From 00ebbffb2829a8f8403139cbf4869192de7cc10b Mon Sep 17 00:00:00 2001 From: iomgaa Date: Fri, 16 Oct 2026 06:47:27 -0400 Subject: [PATCH] =?UTF-8?q?feat(ceus):=20=E4=B8=A4=E5=8D=8A=E5=9D=90?= =?UTF-8?q?=E6=A0=87=E6=8D=A2=E7=AE=97=E6=94=B9=E4=B8=BA=E6=AF=8F=E6=AE=B5?= =?UTF-8?q?=E8=A7=86=E9=A2=91=E9=85=8D=E5=87=86=EF=BC=8C=E5=8A=A0=E7=9B=AE?= =?UTF-8?q?=E8=A7=86=E6=A0=B8=E5=AF=B9=E8=84=9A=E6=9C=AC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 第 1 步分屏后在两半扇形区内用边缘图归一化互相关估 (dx, dy)，再用逐帧多块模板匹配的中位数独立核对； 两者各分量之差超过门槛（15 像素）这一例第 1 步失败。删掉「两个视图左上角相减」的实现与测试。 - 门槛依据：34 例中走完第 1 步的 29 例，差值中位 2、95 分位 6.2、最大 13；参数与门槛在同一批视频上调出。 - scripts/ceus_align_review.py：只做第 1 步，出扇形区棋盘图与边缘叠加图，混入错开 15–30 像素的 阴性对照、随机编号；计分时任一对照未被判「没对齐」即作废。 - ceus_quant 读 review.csv：不通过的病例第 1 步失败；表里没有这段视频或平移对不上即停。 - 配准与紫圈颜色参数进分割缓存键。 Co-Authored-By: Claude Opus 5.5 --- config/ceus.yaml | 39 ++++- pyproject.toml | 7 +- research-wiki/explanation/结构.md | 26 ++- scripts/ceus_align_review.py | 250 +++++++++++++++++++++++++++ scripts/ceus_quant.py | 128 ++++++++++++-- src/cadx/ceus/__init__.py | 26 ++- src/cadx/ceus/_config.py | 48 +++++- src/cadx/ceus/_mapping.py | 51 ++---- src/cadx/ceus/_pipeline.py | 41 +++-- src/cadx/ceus/_polygons.py | 27 ++- src/cadx/ceus/_registration.py | 277 ++++++++++++++++++++++++++++++ tests/fake_ceus_video.py | 57 ++++++ tests/fake_split_scene.py | 92 ++++++++++ tests/test_ceus_align_review.py | 196 +++++++++++++++++++++ tests/test_ceus_contrast_side.py | 4 +- tests/test_ceus_mapping.py | 137 +++++++-------- tests/test_ceus_pipeline.py | 60 ++----- tests/test_ceus_polygons.py | 48 +++--- tests/test_ceus_quant.py | 119 +++++++++---- tests/test_ceus_registration.py | 144 ++++++++++++++++ 20 files changed, 1500 insertions(+), 277 deletions(-) create mode 100644 scripts/ceus_align_review.py create mode 100644 src/cadx/ceus/_registration.py create mode 100644 tests/fake_ceus_video.py create mode 100644 tests/fake_split_scene.py create mode 100644 tests/test_ceus_align_review.py create mode 100644 tests/test_ceus_registration.py diff --git a/config/ceus.yaml b/config/ceus.yaml index 791840e..23c79e9 100644 --- a/config/ceus.yaml +++ b/config/ceus.yaml @@ -1,8 +1,9 @@ # CEUS 定量的配置，给 scripts/ceus_quant.py 用。 # -# 分三段：paths（位置）、segment（分割设置，任务书 02 第二阶段加）、ceus（算法参数）。 +# 分四段：paths（位置）、align_check（配准目视核对的出图，scripts/ceus_align_review.py 用）、 +# segment（分割设置，任务书 02 第二阶段加）、ceus（算法参数）。 # ceus 段的字段定义和校验在 src/cadx/ceus/_config.py，segment 段在 src/cadx/segment/_config.py； -# paths 段和最外层由 scripts/ceus_quant.py 定义。路径都相对仓库根目录。 +# paths、align_check 两段和最外层由 scripts/ceus_quant.py 定义。路径都相对仓库根目录。 # 每个字段必填，可以没有的值写 null。 # # ceus 段的数值逐个照抄旧项目 config/ceus_roi_default.yaml，没有改动（待定问题/02：参数何时冻结）。 @@ -17,6 +18,15 @@ paths: delivery_dir: data/交付包_20260908_40例 # 交付包位置，相对仓库根目录 runs_root: runs/ceus # 每次运行一个子目录，只追加不删 work_root: work # 可删掉重算的中间文件（分割缓存） + align_runs_root: runs/ceus_align # 配准目视核对，每次出图一个子目录，只追加不删 + # 配准目视核对的结果表（scripts/ceus_align_review.py score 写出的 review.csv）。ceus_quant 按视频 + # SHA-256 查它，文件不存在就报错停下。下面是占位路径：目视核对通过后改成实际路径 + align_review: runs/ceus_align/目视核对通过后改成实际运行号/review.csv + +align_check: # 配准目视核对的出图（任务书 02「配准的目视核对」）。本仓库新加 + seed: 20261016 # 挑阴性对照、定错开方向和编号用的随机数种子。本仓库新加，取任务日期，无其他含义 + control_count: 10 # 阴性对照约 10 例。出自任务书 02「配准的目视核对」第 1 步 + control_shift_px: [15, 30] # 对照故意错开的距离（像素，欧氏长度）区间，方向随机。出自同上 segment: # 分割设置，字段定义在 src/cadx/segment/_config.py。出自旧项目 debugs/medsam2_local_track.py medsam2_dir: reference/MedSAM2 # MedSAM2 代码，提交号从它的 .git 读。旧脚本是远程 /home/wangwenhan/MedSAM2 @@ -38,6 +48,31 @@ ceus: split_sample_frames: 25 # 找分屏线只用均匀抽的 25 帧。旧 debug 脚本写死的数，旧项目工程侧定的，无出处 split_search_band: [0.35, 0.65] # 只在宽度的这一段里找分屏线。旧 detect_split_x 的参数默认值，旧项目工程侧定的，无出处 + # 两半配准（第 1 步）：每段视频单独估左半到右半的平移（任务书 02 第一阶段第 1 条）。整段是本仓库 + # 新加的，旧项目用写死的锚点、没有这一步。取值依据 2026-10-16 在 8 段真实视频（覆盖 5 种帧宽）上的 + # 只读实测和同日 34 例的差值分布，见各行 + registration: + sample_frames: 30 # 均匀抽 30 帧求时间平均和标准差。本仓库新加，2026-10-16 实测取的，无外部出处 + top_ignore_ratio: 0.0834 # 顶部信息栏，与 prompt_source.ignore_regions 那条文字带同高。本仓库新加，依据 2026-10-16 实测 + # 时间标准差高于 4（8 位灰度）算扇形区：静止的文字、色标、刻度标准差近 0。本仓库新加，依据 + # 2026-10-16 实测 + fan_std_threshold: 4.0 + fan_open_kernel: 5 # 本仓库新加，依据 2026-10-16 实测 + mark_dilate: 3 # 绿框紫圈像素外扩 3 像素再修补。本仓库新加，依据 2026-10-16 实测 + blur_sigma: 2.0 # 本仓库新加，依据 2026-10-16 实测 + # 搜索范围：|dx| ≤ 0.5 × 左半宽、|dy| ≤ 0.15 × 左半高。8 段实测水平错位 −10 到 +182 像素 + # （最大约 0.38 倍半宽）、竖直 0–5 像素。本仓库新加，依据 2026-10-16 实测 + search_x_ratio: 0.5 + search_y_ratio: 0.15 + template_grid: [4, 4] # 本仓库新加，依据 2026-10-16 实测（3 块时有视频两块误配，改多块取中位数） + template_size_ratio: 0.15 # 本仓库新加，依据 2026-10-16 实测 + template_min_patches: 3 # 本仓库新加，无外部出处 + # 两种方法估计各分量之差的最大值的门槛（像素）。依据：2026-10-16 在交付包检查全部通过的 34 例中 + # 走完第 1 步的 29 例上，两种方法差值最小 0、25% 分位 1、中位 2、75% 分位 2、90% 分位 5、 + # 95% 分位 6.2、最大 13；15 高于观测最大值，远低于唯一明显错配的 66。参数与门槛在同一批视频上 + # 调出，没有独立验证集。负责人 2026-10-16 定 + max_disagreement_px: 15 + prompt_source: # 出自旧项目：findings §3.2 绿框检测 type: color_polygon # 出自旧项目 findings §3.2 hsv_lower: [35, 60, 60] # 出自旧项目 findings §3.2 diff --git a/pyproject.toml b/pyproject.toml index 0166b2c..5daf4c1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -62,12 +62,13 @@ typeCheckingMode = "strict" venvPath = "." venv = ".venv" # medsam2 以 setuptools 的导入钩子方式可编辑安装，pyright 跟不过去，直接指给它源码目录。 -extraPaths = ["tests", "reference/MedSAM2"] +extraPaths = ["tests", "scripts", "reference/MedSAM2"] [tool.pytest.ini_options] testpaths = ["tests"] -# 测试辅助（假交付包生成器）放在 tests/ 下，按普通模块导入。 -pythonpath = ["tests"] +# 测试辅助（假交付包生成器）放在 tests/ 下，按普通模块导入；scripts/ 也在路径上， +# ceus_align_review 引用 ceus_quant 的配置模型（直接运行脚本时 scripts/ 本来就在 sys.path 上）。 +pythonpath = ["tests", "scripts"] addopts = "--strict-markers --import-mode=importlib" markers = [ "live: 真实模型调用，花钱，需要 .env 里的密钥", diff --git a/research-wiki/explanation/结构.md b/research-wiki/explanation/结构.md index 7649843..96c8368 100644 --- a/research-wiki/explanation/结构.md +++ b/research-wiki/explanation/结构.md @@ -193,13 +193,13 @@ library/<构建号>/ ### CEUS 定量的产物（`runs/ceus/` 与 `work/`） -每次运行的产物写进 `runs/ceus/<运行号>/`，只追加；第 3 步的分割结果按下表的键缓存在 `work/`，可以删掉重算。具体布局见 `../任务书/02-CEUS定量.md`。 +每次运行的产物写进 `runs/ceus/<运行号>/`，只追加；全量运行之前的配准目视核对写进 `runs/ceus_align/<运行号>/`（图、`key.csv`、判读、`review.csv`），ceus_quant 按 `review.csv` 拦下核对没通过的病例；第 3 步的分割结果按下表的键缓存在 `work/`，可以删掉重算。具体布局见 `../任务书/02-CEUS定量.md`。 CEUS 定量按下面的阶段进行（旧项目的阶段划分，迁移时逐模块核对）： | 阶段 | 做什么 | 跨运行复用 | | --- | --- | --- | -| 1 解码与分屏 | 读视频；超声机录造影时屏幕左右分两半，一半造影画面、一半普通灰阶画面，找出分界 | 否 | +| 1 解码与分屏 | 读视频；超声机录造影时屏幕左右分两半，一半造影画面、一半普通灰阶画面，找出分界；再配准两半，估出两半局部坐标之间的平移（扇形在两半里的位置不同，每段视频单独估） | 否 | | 2 提示帧 | 在造影画面上找医生画的绿框，挑出绿框位置有变化的帧，换算到另一半画面上 | 否 | | 3 分割 | segment 用 MedSAM2 把提示传播成逐帧病灶掩码；需要显卡，最贵 | **是**，按「视频内容哈希 + 阶段 1–3 的参数 + 模型权重哈希」存 | | 4 多边形 | 掩码转成平滑轮廓，与绿框融合 | 否 | @@ -253,7 +253,7 @@ X 是超声报告工作簿的「推注开始时间s」，含义是视频第 0 | --- | --- | --- | | `build.yaml` | `scripts/build.py` | `delivery_dir`（交付包位置）、`library_root`、`answers_root`（两个库的上级目录）、`alias_seed`（生成化名的随机种子）、`decode_codes`（编码解码开关）、`code_maps`（映射表位置或 null）、`ceus_run`（CEUS 运行目录或 null）；字段定义和校验在 `src/cadx/build/_config.py`，逐项说明见 `../guides/装配两个库.md` | | `code_maps.yaml` | build（经 `build.yaml` 引用） | 编码映射表：交付包里每个用数字编码的列，每个数字代表什么（如性别的 1、2）。还没向临床方拿到，尚未建 | -| `ceus.yaml` | `scripts/ceus_quant.py` | CEUS 定量的全部参数，每个参数用注释写明出处 | +| `ceus.yaml` | `scripts/ceus_quant.py`、`scripts/ceus_align_review.py` | CEUS 定量的全部参数和配准目视核对的出图设置，每个参数用注释写明出处 | | `consult.yaml` | `scripts/consult.py` | 会诊的实验配置与模型用途 | | `evaluate.yaml` | `scripts/evaluate.py` | 评价的指标与统计设置 | @@ -298,7 +298,7 @@ Python 环境用 uv；安装与检查命令见 `../guides/开发环境.md`。 | 模块 | 预期行数 | 实际行数 | 旧项目对应部分的行数 | | --- | --- | --- | --- | -| ceus | 6000–7500（基本照搬旧算法，去掉只为产品层服务的部分） | — | `ceus_roi/` 除 `service.py` 7330，`video_decoder.py` 135 | +| ceus | 6000–7500（基本照搬旧算法，去掉只为产品层服务的部分） | 7085 | `ceus_roi/` 除 `service.py` 7330，`video_decoder.py` 135 | | segment | 150–300 | 233 | `medsam2_local_track.py` 240 | | delivery | 400–700 | 739 | `clinical.py` 266、`modalities.py` 522、`workbook.py` 99、`delivery.py` 134，共 1021，其中入库与数据库部分不迁 | | library | 150–300 | 300 | `input_packages.py` 104 | @@ -310,7 +310,7 @@ Python 环境用 uv；安装与检查命令见 `../guides/开发环境.md`。 | runs | 100–250 | 86（只含快照部分） | — | | consult | 500–1000（会诊口径定后重估） | — | — | | evaluate | 300–600（评价口径定后重估） | — | 已删的 `core/eval` 622 | -| scripts/ 每个 | ≤ 300 | `build.py` 37、`scan_library.py` 88、`ceus_quant.py` 293 | `t1_local_debug.py` 762、`tic_lognormal_replay.py` 627 | +| scripts/ 每个 | ≤ 300 | `build.py` 37、`scan_library.py` 88、`ceus_quant.py` 393（超出，含目视核对表的读取与核对）、`ceus_align_review.py` 250 | `t1_local_debug.py` 762、`tic_lognormal_replay.py` 627 | --- @@ -377,14 +377,16 @@ Python 环境用 uv；安装与检查命令见 `../guides/开发环境.md`。 对外两段接口，中间由 `scripts/ceus_quant.py` 调 segment 做第 3 步： -- `prepare_segmentation(video_path, cfg)`：第 1–2 步，交出 `SegmentationInput`；视频打不开或读出的不合规（fps、尺寸、时间戳、帧数）抛 `VideoDecodeError`、找不到可信分屏线抛 `SplitDetectionError`（都是第 1 步），两半都有或都没有绿框抛 `ContrastSideError`（第 2 步）。 -- `SegmentationInput`：视频的 fps、帧数、宽高、逐帧时间戳，分屏线，`contrast`（`ContrastSide`：哪一半是造影画面、造影视图和分割那一半视图 `target_view`、提示帧），`target_frames`（分割那一半的逐帧 BGR 图像），`prompt_masks`（提示帧帧号 → 绿框挪到 `target_view` 局部坐标的掩膜），造影那一半的第一帧。 +- `prepare_segmentation(video_path, cfg)`：第 1–2 步，交出 `SegmentationInput`；视频打不开或读出的不合规（fps、尺寸、时间戳、帧数）抛 `VideoDecodeError`、找不到可信分屏线抛 `SplitDetectionError`、两半配准失败（找不到扇形区、可用模板匹配太少，或两种方法估出的平移相差超过 `registration.max_disagreement_px`）抛 `RegistrationError`（都是第 1 步），两半都有或都没有绿框抛 `ContrastSideError`（第 2 步）。 +- `register_halves(frames, left_view, right_view, cfg.registration, cfg.prompt_source, cfg.pan_source)`：第 1 步的两半配准，交出 `HalfRegistration`：左半局部坐标 → 右半局部坐标的平移，`correlation`（两半扇形区边缘图的归一化互相关，采用它）、`template`（左半扇形内多块模板逐帧匹配后取中位数，只用来核对）、各次匹配明细、两者各分量之差的最大值；`contrast_to_target(side)` 换成造影那一半 → 另一半的平移。 +- `SegmentationInput`：视频的 fps、帧数、宽高、逐帧时间戳，分屏线，`registration`（`HalfRegistration`），`offset`（造影那一半局部坐标 + `offset` = `target_view` 局部坐标），`contrast`（`ContrastSide`：哪一半是造影画面、造影视图和分割那一半视图 `target_view`、提示帧），`target_frames`（分割那一半的逐帧 BGR 图像），`prompt_masks`（提示帧帧号 → 绿框按 `offset` 挪到 `target_view` 局部坐标的掩膜），造影那一半的第一帧。 - `quantify_case(prepared, target_masks, injection_offset_s, cfg)`：第 4–9 步。`target_masks` 是分割在 `target_view` 局部坐标上的 `(帧数, 高, 宽)` 布尔数组，`injection_offset_s` 是推注时间 X；分割一帧有效轮廓都没有抛 `NoValidContourError`（第 4 步）。拟合没过质量门不抛错。 - `QuantResult`：`trajectory`（`Trajectory`）、`tic`（`TicResult`）、`curves`（`TicCurves`）、推注时间、病灶两条曲线拟合是否都过质量门、`anchors`（7 个槽位的 `AnchorFrameLocation`）、`perfusion`（`roi_all`、`roi_perf` 各自的 `LognormalPerfusionMetrics`，拟合没过为 None）、紫圈检出率。 - `AnchorFrameLocation`：一个关键帧槽位的帧号、视频时间、推注后时刻、状态（`ok` 或未定位的类别）、未定位的原因、锚点时刻。 - 写产物用：`pack_trajectory(trajectory)`（trajectory.npz 字节）、`pack_tic_curves(curves)`（tic_curves.npz 字节）、`render_tic_png(tic, markers=tic_plot_markers(anchors))`（tic.png，标出 D1–D4）、`render_case_overlay(prepared, result)`（overlay.mp4 字节）、`keyframe_images(prepared, anchors)`（定位到的关键帧整幅原图，槽位 → 图像）、`quant_summary(result)`（拟合、灌注指标、紫圈检出率、7 个关键帧定位，可直接写成 JSON）；读回用 `unpack_trajectory`、`unpack_tic_curves`。 -- `CeusConfig` 及各段的配置类：`config/ceus.yaml` 里 `ceus` 段的字段定义（pydantic，字段必填、不许多写、读入后不可改）；`paths` 段由 `scripts/ceus_quant.py` 定义。 -- `CeusCaseError`：只落在这一例上的失败，带失败的步号和原因；`VideoDecodeError`、`SplitDetectionError`、`ContrastSideError`、`NoValidContourError` 是它的子类。脚本据此把这一例记进 `failed.csv`、接着跑下一例。 +- `fan_mask(std, cfg.registration)`：由一半画面的逐像素时间标准差求扇形区掩膜；`color_marks(frames, prompt_cfg, pan_cfg, dilate)`：任一帧里落在绿框或紫圈颜色范围内的像素。两者配准用，目视核对出图也用。 +- `CeusConfig` 及各段的配置类：`config/ceus.yaml` 里 `ceus` 段的字段定义（pydantic，字段必填、不许多写、读入后不可改）；`paths`、`align_check` 两段由 `scripts/ceus_quant.py` 定义。 +- `CeusCaseError`：只落在这一例上的失败，带失败的步号和原因；`VideoDecodeError`、`SplitDetectionError`、`RegistrationError`、`ContrastSideError`、`NoValidContourError` 是它的子类。脚本据此把这一例记进 `failed.csv`、接着跑下一例。 - `VideoFile(path)`：按路径用 cv2 读视频：`probe()` 给 fps、帧数、宽高；`iter_frames()` 每次从第 0 帧逐帧读；`frame_times_s()` 给逐帧显示时间戳。 - 其余导出（几何函数、各步的算法函数与结果类型，如 `split_views`、`locate_contrast_side`、`masks_to_smoothed_polygons`、`extract_dark_region_polygons`、`select_reference_region`、`compute_tic_for_video`、`fit_lognormal_curve`、`locate_anchor_frames`、`compute_lognormal_perfusion_metrics`、`render_overlay`）是两段接口的组成部分，供测试用。 @@ -395,6 +397,12 @@ Python 环境用 uv；安装与检查命令见 `../guides/开发环境.md`。 - `medsam2_commit(medsam2_dir)`：读 MedSAM2 代码目录 `.git` 的 HEAD 提交号，给快照和分割缓存键用。 - `cuda_peak_memory_bytes(device)`：本进程在该设备上分配过的显存峰值（字节），给运行快照用；非 CUDA 设备为 None，本进程没初始化过 CUDA（分割全读了缓存）为 0。放在 segment 里，是因为 `cadx` 里只有它引用 `torch`。 +### CEUS 的两个脚本（`scripts/`） + +- `ceus_quant.py --config config/ceus.yaml`：逐例跑第 1–9 步；第 1 步后用 `check_review` 按视频 SHA-256 查 `paths.align_review`（`read_review` 读入）：没通过的病例第 1 步失败，表里没有这段视频或配准结果与表里不同就整次停下。 +- `ceus_align_review.py render --config config/ceus.yaml`：只做第 1 步，为每例和随机挑出的阴性对照出棋盘拼图和边缘叠加图，按打乱的编号写进 `runs/ceus_align/<运行号>/images/`，对应关系写 `key.csv`。 +- `ceus_align_review.py score --run-dir <运行目录>`：对照 `verdicts_a.csv`、`verdicts_b.csv` 计分；有对照没被两人都判「没对齐」就作废（退出码 1，不写表），否则写 `review.csv`。 + ### build（`src/cadx/build/`） - `run_build(config_path, repo_root, now)`：完整做一次构建：读配置和映射表、读交付包、对每例做 C14（`check_key_image` 返回原因的整例跳过，记进跳过清单）；`ceus_run` 不为 null 时核对运行目录（不合规整次停下），对每例做 C15（CEUS 失败或没有结果的整例跳过）、关键帧也过 C14，按白名单构造 tic.json；取化名、先写进 `.partial` 目录再改名，交出 `BuildSummary`。 diff --git a/scripts/ceus_align_review.py b/scripts/ceus_align_review.py new file mode 100644 index 0000000..4f18592 --- /dev/null +++ b/scripts/ceus_align_review.py @@ -0,0 +1,250 @@ +"""配准的目视核对（任务书 02 第三阶段）：`render --config ` 出图，`score --run-dir` +计分。 + +render 只做第 1 步、不用显卡：第 1 步成功的每例出棋盘拼图和边缘叠加图（只截扇形区），另随机挑 +`align_check` 段定的例数把配准故意错开当阴性对照；图按打乱的编号放进 `images/`，对应关系写 +`key.csv`（不给判读者）。score 对照 `verdicts_a.csv`、`verdicts_b.csv`：任一张对照被任一人判成 +「对齐」或「看不清」就作废（不写 review.csv，退出码 1）；否则写 `review.csv`，两人都判「对齐」才 +通过。终端只打印计数。 +""" + +import argparse +import csv +import sys +from dataclasses import dataclass +from datetime import datetime +from pathlib import Path +from typing import Any + +import cv2 +import numpy as np +import numpy.typing as npt + +from cadx import ceus +from cadx.ceus import CeusConfig, HalfRegistration, ViewRect, crop_frame +from cadx.delivery import read_delivery +from cadx.runs import make_build_id, read_git_state, write_snapshot +from ceus_quant import REPO_ROOT, REVIEW_HEADER, file_sha256, load_config, resolve_paths + +KEY_HEADER = ("id", "video_id", "is_control", "dx", "dy", "video_sha256") +VERDICT_HEADER = ("编号", "判定", "依据") +ALIGNED, MISALIGNED, UNCLEAR = "对齐", "没对齐", "看不清" +_EDGE_FRACTION = 0.1 # 每半取扇形内梯度幅值最高的 10% 像素画成边缘 +type _U8 = npt.NDArray[np.uint8] +type _Gray = npt.NDArray[np.float32] +type _Mask = npt.NDArray[np.bool_] + + +def _indices(count: int, wanted: int) -> set[int]: + """与 `sample_frames_for_split` 相同的均匀抽帧下标。""" + if count <= wanted: + return set(range(count)) + return {int(i) for i in np.linspace(0, count - 1, num=wanted, dtype=np.int32).tolist()} + + +@dataclass(frozen=True) +class Material: + """一例出图要的东西：两个视图；整帧的时间平均灰度、抽样中间那帧的灰度、绿框紫圈像素；两半扇形。""" + + left: ViewRect + right: ViewRect + mean: _Gray + middle: _Gray + marks: _Mask + fan_left: _Mask + fan_right: _Mask + + +def step1(video_path: Path, cfg: CeusConfig) -> tuple[Material, HalfRegistration]: + """第 1 步：逐帧解码只留分屏和配准要的抽样帧，找分屏线，配准两半。""" + video = ceus.VideoFile(video_path) + _, frame_count, _, _ = video.probe() + split_at = _indices(frame_count, cfg.views.split_sample_frames) + register_at = _indices(frame_count, cfg.registration.sample_frames) + split_frames, frames = list[_U8](), list[_U8]() + decoded = 0 + for index, frame in video.iter_frames(): + decoded += 1 + if index in split_at: + split_frames.append(frame) + if index in register_at: + frames.append(frame) + if decoded != frame_count: + raise ceus.VideoDecodeError(f"解码帧数 {decoded} 与容器记录的帧数 {frame_count} 不一致") + _, left, right = ceus.split_views(split_frames, cfg.views) + reg = ceus.register_halves( + frames, left, right, cfg.registration, cfg.prompt_source, cfg.pan_source + ) + gray = np.stack([cv2.cvtColor(f, cv2.COLOR_BGR2GRAY) for f in frames]).astype(np.float32) + std = gray.std(axis=0) + fans = [ceus.fan_mask(crop_frame(std, view), cfg.registration) for view in (left, right)] + marks = ceus.color_marks( + frames, cfg.prompt_source, cfg.pan_source, cfg.registration.mark_dilate + ) + material = Material(left, right, gray.mean(0), gray[len(frames) // 2], marks, *fans) + return material, reg + + +def _shift(image: npt.NDArray[Any], offset: tuple[int, int], size: tuple[int, int]) -> Any: + """右半按平移挪进左半坐标：结果 (x, y) = 右半 (x + dx, y + dy)，出界补 0。""" + matrix = np.array([[1, 0, -offset[0]], [0, 1, -offset[1]]], dtype=np.float32) + return cv2.warpAffine(image, matrix, size, flags=cv2.INTER_NEAREST) + + +def _edges(gray: _Gray, fan: _Mask, marks: _Mask, sigma: float) -> Any: + """扇形内、绿框紫圈之外梯度幅值最高的 `_EDGE_FRACTION` 像素，膨胀 1 像素便于看出重合。""" + smooth = cv2.GaussianBlur(gray, (0, 0), sigma) + magnitude = np.hypot(cv2.Sobel(smooth, cv2.CV_32F, 1, 0), cv2.Sobel(smooth, cv2.CV_32F, 0, 1)) + keep = fan & ~marks # 绿框紫圈只画在一半上，它们的边缘不参与比对 + edge = keep & (magnitude >= np.quantile(magnitude[keep], 1 - _EDGE_FRACTION)) + return cv2.dilate(edge.astype(np.uint8), np.ones((3, 3), np.uint8)) > 0 + + +def render_pair(m: Material, offset: tuple[int, int], sigma: float) -> tuple[_U8, _U8]: + """按左半 → 右半的平移出（棋盘拼图，边缘叠加图），截到两半扇形的外接框，扇形外涂黑。""" + left, right, size = m.left, m.right, (m.left.width, m.left.height) + fan_r = _shift(m.fan_right.astype(np.uint8), offset, size) > 0 + support = m.fan_left | fan_r + ys, xs = np.nonzero(support) + x1, x2, y1, y2 = int(xs.min()), int(xs.max()) + 1, int(ys.min()), int(ys.max()) + 1 + block = max(8, round((x2 - x1) / 12)) + yy, xx = np.mgrid[0 : left.height, 0 : left.width] + first = ((yy - y1) // block + (xx - x1) // block) % 2 == 0 + checker = np.where( + first, crop_frame(m.middle, left), _shift(crop_frame(m.middle, right), offset, size) + ) + checker[~support] = 0 + overlay = np.zeros((left.height, left.width, 3), dtype=np.uint8) + marks_l = crop_frame(m.marks, left) + marks_r = _shift(crop_frame(m.marks, right).astype(np.uint8), offset, size) > 0 + edges_l = _edges(crop_frame(m.mean, left), m.fan_left, marks_l, sigma) + edges_r = _edges(_shift(crop_frame(m.mean, right), offset, size), fan_r, marks_r, sigma) + overlay[..., 2][edges_l] = 255 # 红：左半 + overlay[..., 1][edges_r] = 255 # 绿：右半；重合处红 + 绿 = 黄 + return np.clip(checker[y1:y2, x1:x2], 0, 255).astype(np.uint8), overlay[y1:y2, x1:x2] + + +def render(config_path: Path) -> None: + started_at = datetime.now().astimezone() + cfg = resolve_paths(load_config(config_path)) + git = read_git_state(REPO_ROOT) + run_id = make_build_id(started_at, git) + run_dir = cfg.paths.align_runs_root / run_id + partial = cfg.paths.align_runs_root / f"{run_id}.partial" + for path in (run_dir, partial): + if path.exists(): + raise FileExistsError(f"运行目录已存在，不覆盖：{path}") + (partial / "images").mkdir(parents=True) + rng = np.random.default_rng(cfg.align_check.seed) + + cases: list[tuple[str, str, Material, tuple[int, int]]] = [] + failures: list[tuple[str, int, str]] = [] + for case in read_delivery(cfg.paths.delivery_dir, None).cases: + try: + material, registration = step1(case.video, cfg.ceus) + except ceus.CeusCaseError as error: + failures.append((case.video_id, error.step, error.reason)) + continue + cases.append((case.video_id, file_sha256(case.video), material, registration.correlation)) + # 对照例数多于成功例数时 rng.choice 报错（不放回抽样）。 + controls = rng.choice(len(cases), size=cfg.align_check.control_count, replace=False) + items: list[tuple[int, bool, tuple[int, int]]] = [(i, False, c[3]) for i, c in enumerate(cases)] + for index in sorted(int(i) for i in controls): + length = rng.uniform(*cfg.align_check.control_shift_px) + angle = rng.uniform(0.0, 2 * np.pi) + dx, dy = cases[index][3] + shift = (round(length * np.cos(angle)), round(length * np.sin(angle))) + items.append((index, True, (dx + shift[0], dy + shift[1]))) + numbers = rng.permutation(len(items)) + 1 + rows: list[tuple[str, str, str, int, int, str]] = [] + for number, (index, is_control, offset) in zip(numbers.tolist(), items, strict=True): + video_id, sha, material, _ = cases[index] + name = f"{number:03d}" + checker, edges = render_pair(material, offset, cfg.ceus.registration.blur_sigma) + for kind, image in (("checker", checker), ("edges", edges)): + if not cv2.imwrite(str(partial / "images" / f"{name}_{kind}.png"), image): + raise RuntimeError(f"写 PNG 失败：{name}_{kind}.png") + rows.append((name, video_id, str(is_control).lower(), offset[0], offset[1], sha)) + with (partial / "key.csv").open("w", encoding="utf-8", newline="") as handle: + csv.writer(handle).writerows([KEY_HEADER, *sorted(rows)]) + with (partial / "failed.csv").open("w", encoding="utf-8", newline="") as handle: + csv.writer(handle).writerows([("video_id", "step", "reason"), *failures]) + finished_at, lock = datetime.now().astimezone(), REPO_ROOT / "uv.lock" + configs = {"ceus.yaml": config_path} # align_check 段的随机数种子随 ceus.yaml 进快照 + write_snapshot(partial / "snapshot", run_id, git, started_at, finished_at, lock, configs) + if run_dir.exists(): + raise FileExistsError(f"运行目录已存在，不覆盖：{run_dir}") + partial.rename(run_dir) + print(f"运行目录：{run_dir}") + print(f"第 1 步成功 {len(cases)} 例、失败 {len(failures)} 例；对照 {len(controls)} 组") + + +def _read_csv(path: Path, header: tuple[str, ...]) -> list[list[str]]: + with path.open(encoding="utf-8", newline="") as handle: + rows = list(csv.reader(handle)) + if not rows or tuple(rows[0]) != header or any(len(r) != len(header) for r in rows): + raise ValueError(f"{path} 的表头应为 {header}，且每行 {len(header)} 列") + return rows[1:] + + +def read_verdicts(path: Path, ids: set[str]) -> dict[str, str]: + """一位判读者的判定，编号 → 判定。缺编号、多编号、重复编号、判定取值不合法都报错。""" + rows = _read_csv(path, VERDICT_HEADER) + verdicts = {row[0]: row[1] for row in rows} + if len(verdicts) != len(rows): + raise ValueError(f"{path} 有重复的编号") + if set(verdicts) != ids: + missing, extra = sorted(ids - set(verdicts)), sorted(set(verdicts) - ids) + raise ValueError(f"{path} 的编号与出图不符：缺 {missing}，多 {extra}") + bad = sorted({v for v in verdicts.values() if v not in (ALIGNED, MISALIGNED, UNCLEAR)}) + if bad: + raise ValueError(f"{path} 的判定只能是对齐、没对齐、看不清，收到 {bad}") + return verdicts + + +def score(run_dir: Path) -> int: + """计分，返回退出码：0 写出 review.csv，1 本次核对作废。""" + review_path = run_dir / "review.csv" + if review_path.exists(): + raise FileExistsError(f"review.csv 已存在，不覆盖：{review_path}") + key = _read_csv(run_dir / "key.csv", KEY_HEADER) + ids = {row[0] for row in key} + a = read_verdicts(run_dir / "verdicts_a.csv", ids) + b = read_verdicts(run_dir / "verdicts_b.csv", ids) + controls = [row[0] for row in key if row[2] == "true"] + real = [row for row in key if row[2] == "false"] + caught_a, caught_b = (sum(v[i] == MISALIGNED for i in controls) for v in (a, b)) + print(f"阴性对照 {len(controls)} 张：判读者 A 判「没对齐」{caught_a} 张，B {caught_b} 张") + if caught_a < len(controls) or caught_b < len(controls): + print("本次核对作废：有阴性对照被判成「对齐」或「看不清」，不写 review.csv") + return 1 + if len({row[5] for row in real}) != len(real): + raise ValueError("key.csv 里同一段视频作为真实病例出现了两次") + passed = [a[row[0]] == ALIGNED and b[row[0]] == ALIGNED for row in real] + lines = [ + (r[5], r[3], r[4], str(ok).lower(), a[r[0]], b[r[0]]) + for r, ok in zip(real, passed, strict=True) + ] + with review_path.open("x", encoding="utf-8", newline="") as handle: + csv.writer(handle).writerows([REVIEW_HEADER, *lines]) + unclear = sum(UNCLEAR in (a[row[0]], b[row[0]]) for row in real) + differ = sum(a[row[0]] != b[row[0]] for row in real) + print(f"真实病例 {len(real)} 例：通过 {sum(passed)}，不通过 {len(real) - sum(passed)}") + print(f"其中有人判「看不清」{unclear} 例，两人判定不一致 {differ} 例") + return 0 + + +def main(argv: list[str]) -> None: + parser = argparse.ArgumentParser(description="配准的目视核对：出图、计分") + commands = parser.add_subparsers(dest="command", required=True) + commands.add_parser("render", help="出图").add_argument("--config", type=Path, required=True) + commands.add_parser("score", help="计分").add_argument("--run-dir", type=Path, required=True) + args = parser.parse_args(argv) + if args.command == "render": + render(args.config) + elif score(args.run_dir) != 0: + raise SystemExit(1) + + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/scripts/ceus_quant.py b/scripts/ceus_quant.py index 59da9a5..3aeef98 100644 --- a/scripts/ceus_quant.py +++ b/scripts/ceus_quant.py @@ -1,6 +1,9 @@ """CEUS 定量：对交付包里检查全部通过的病例逐例跑第 1–9 步，写进 `runs/ceus/<运行号>/`。 用法：uv run python scripts/ceus_quant.py --config config/ceus.yaml +第 1 步算出两半配准后，按视频 SHA-256 查配准目视核对的结果表（`paths.align_review`，由 +`scripts/ceus_align_review.py score` 写出）：核对没通过的病例第 1 步失败；表里没有这段视频，或配准 +结果与表里记的不同，说明核对的不是这次跑的东西，整次运行停下。 第 3 步（MedSAM2 分割）的结果按内容键缓存在 `/segment_cache/<键>/masks.npz`，键相同就 直接读缓存、不调显卡。某一例在第 1–6 步失败（`CeusCaseError`）记进 `failed.csv`，接着跑下一例； 其余异常整次运行停下，留下 `.partial` 目录。终端只打印计数，不打印视频号。 @@ -13,7 +16,8 @@ import json import sys import time from collections import Counter -from collections.abc import Callable, Hashable, Sequence +from collections.abc import Callable, Hashable, Mapping, Sequence +from dataclasses import dataclass from datetime import datetime from pathlib import Path from typing import Annotated, Any @@ -22,7 +26,7 @@ import cv2 import numpy as np import numpy.typing as npt import yaml -from pydantic import BaseModel, BeforeValidator, ConfigDict, Field +from pydantic import BaseModel, BeforeValidator, ConfigDict, Field, field_validator from cadx.ceus import ( ANCHOR_SLOTS, @@ -45,14 +49,20 @@ from cadx.segment import SegmentConfig, cuda_peak_memory_bytes, medsam2_commit, REPO_ROOT = Path(__file__).resolve().parents[1] -# 第 1–2 步（`prepare_segmentation`）读的 ceus 段子段：views（分屏）、prompt_source（绿框颜色）、 -# keyframes（提示帧选择）。它们决定交给分割的帧和提示掩码，进分割缓存键。 -_STEP_1_2_CEUS_FIELDS = frozenset({"views", "prompt_source", "keyframes"}) +# 第 1–2 步（`prepare_segmentation`）读的 ceus 段子段：views（分屏）、registration（两半配准）、 +# prompt_source（绿框颜色，配准和选提示帧都用）、pan_source（紫圈颜色，配准时修补掉）、keyframes +# （提示帧选择）。它们决定交给分割的帧和提示掩码，进分割缓存键。 +_STEP_1_2_CEUS_FIELDS = frozenset( + {"views", "registration", "prompt_source", "pan_source", "keyframes"} +) # segment 段里不改变分割结果的开关，不进缓存键；其余字段都进。 _RESULT_NEUTRAL_SEGMENT_FIELDS = frozenset( {"device", "offload_video_to_cpu", "offload_state_to_cpu"} ) _FAILED_HEADER = ("video_id", "step", "reason") +# 配准目视核对结果表 review.csv 的列：(dx, dy) 是左半 → 右半的平移 +# （`HalfRegistration.correlation`），passed 取 true / false，后两列是两位判读者的判定。 +REVIEW_HEADER = ("video_sha256", "dx", "dy", "passed", "verdict_a", "verdict_b") def _reject_empty_path(value: object) -> object: @@ -67,13 +77,37 @@ type _PathField = Annotated[Path, Field(strict=False), BeforeValidator(_reject_e class PathsConfig(BaseModel): - """`paths` 段：交付包、运行目录的上级、中间文件目录。相对路径相对仓库根目录。""" + """`paths` 段：交付包、运行目录的上级、中间文件目录、配准目视核对的运行目录上级和结果表。 + + 相对路径相对仓库根目录。`align_review` 在 ceus_quant 运行时必须存在，出图时用不到。 + """ model_config = ConfigDict(strict=True, extra="forbid", frozen=True) delivery_dir: _PathField runs_root: _PathField work_root: _PathField + align_runs_root: _PathField + align_review: _PathField + + +class AlignCheckConfig(BaseModel): + """`align_check` 段：配准目视核对出图用的随机数种子、阴性对照例数、对照错开距离区间（像素）。""" + + model_config = ConfigDict(strict=True, extra="forbid", frozen=True) + + seed: int + control_count: int = Field(ge=1) + control_shift_px: list[int] + + @field_validator("control_shift_px") + @classmethod + def _shift_range(cls, value: list[int]) -> list[int]: + if len(value) != 2 or not 0 < value[0] <= value[1]: + raise ValueError( + f"align_check.control_shift_px 必须是 [lo, hi] 且 0 < lo ≤ hi，收到 {value}" + ) + return value class QuantConfig(BaseModel): @@ -82,6 +116,7 @@ class QuantConfig(BaseModel): model_config = ConfigDict(strict=True, extra="forbid", frozen=True) paths: PathsConfig + align_check: AlignCheckConfig segment: SegmentConfig ceus: CeusConfig @@ -106,7 +141,7 @@ def load_config(path: Path) -> QuantConfig: return QuantConfig.model_validate(yaml.load(text, Loader=_UniqueKeyLoader)) # noqa: S506 继承 SafeLoader -def _resolved(config: QuantConfig) -> QuantConfig: +def resolve_paths(config: QuantConfig) -> QuantConfig: """相对路径换成仓库根目录下的绝对路径（绝对路径原样保留）。""" paths, segment = config.paths, config.segment return config.model_copy( @@ -116,6 +151,8 @@ def _resolved(config: QuantConfig) -> QuantConfig: "delivery_dir": REPO_ROOT / paths.delivery_dir, "runs_root": REPO_ROOT / paths.runs_root, "work_root": REPO_ROOT / paths.work_root, + "align_runs_root": REPO_ROOT / paths.align_runs_root, + "align_review": REPO_ROOT / paths.align_review, } ), "segment": segment.model_copy( @@ -136,6 +173,49 @@ def file_sha256(path: Path) -> str: return digest.hexdigest() +@dataclass(frozen=True) +class ReviewEntry: + """review.csv 的一行：核对时的配准结果（左半 → 右半）和是否通过。""" + + offset: tuple[int, int] + passed: bool + + +def read_review(path: Path) -> dict[str, ReviewEntry]: + """读配准目视核对结果表，视频 SHA-256 → 这一行。表头、取值不合规或同一视频出现两次就报错。""" + with path.open(encoding="utf-8", newline="") as handle: + rows = list(csv.reader(handle)) + if not rows or tuple(rows[0]) != REVIEW_HEADER: + raise ValueError(f"{path} 的表头应为 {REVIEW_HEADER}") + table: dict[str, ReviewEntry] = {} + for number, row in enumerate(rows[1:], start=2): + if len(row) != len(REVIEW_HEADER) or row[3] not in ("true", "false"): + raise ValueError(f"{path} 第 {number} 行不合规：{row}") + if row[0] in table: + raise ValueError(f"{path} 第 {number} 行：同一段视频出现两次") + table[row[0]] = ReviewEntry(offset=(int(row[1]), int(row[2])), passed=row[3] == "true") + return table + + +def check_review( + prepared: SegmentationInput, video_sha256: str, review: Mapping[str, ReviewEntry] +) -> None: + """第 1 步配准结果对照目视核对表：没通过抛 `CeusCaseError`（第 1 步）；对不上表就报错停下。""" + if video_sha256 not in review: + raise RuntimeError( + f"配准目视核对表里没有这段视频（SHA-256 {video_sha256}），核对的不是这批视频" + ) + entry = review[video_sha256] + offset = prepared.registration.correlation + if offset != entry.offset: + raise RuntimeError( + f"这段视频（SHA-256 {video_sha256}）的配准结果 {offset} 与目视核对表里的 " + f"{entry.offset} 不同，核对的不是这次跑的东西" + ) + if not entry.passed: + raise CeusCaseError(1, "配准目视核对未通过") + + def segment_cache_key( video_sha256: str, config: QuantConfig, checkpoint_sha256: str, commit: str ) -> str: @@ -178,14 +258,21 @@ def _write_png(path: Path, image: npt.NDArray[np.uint8]) -> None: def _run_case( - case: CaseRecord, cfg: QuantConfig, key_of: Callable[[str], str], run_dir: Path + case: CaseRecord, + cfg: QuantConfig, + key_of: Callable[[str], str], + review: Mapping[str, ReviewEntry], + run_dir: Path, ) -> tuple[tuple[str, ...], bool]: """跑一例并写出产物，返回（定位到的关键帧槽位，分割是否命中缓存）。 - `key_of` 把视频的 SHA-256 换成分割缓存键（其余几样跨病例不变）。 + `key_of` 把视频的 SHA-256 换成分割缓存键（其余几样跨病例不变）；`review` 是配准目视核对表。 + 第 1–2 步都在 `prepare_segmentation` 里，所以第 2 步失败的病例不查核对表。 """ prepared = prepare_segmentation(case.video, cfg.ceus) - key = key_of(file_sha256(case.video)) + video_sha256 = file_sha256(case.video) + check_review(prepared, video_sha256, review) + key = key_of(video_sha256) masks, hit = _masks(prepared, key, cfg.paths.work_root / "segment_cache", cfg.segment) result = quantify_case(prepared, masks, float(case.bolus_start_s), cfg.ceus) @@ -193,7 +280,19 @@ def _run_case( (case_dir / "keyframes").mkdir(parents=True, exist_ok=False) (case_dir / "trajectory.npz").write_bytes(pack_trajectory(result.trajectory)) (case_dir / "tic_curves.npz").write_bytes(pack_tic_curves(result.curves)) - summary = {"video_id": case.video_id, "segment_cache_key": key, **quant_summary(result)} + registration = prepared.registration + summary = { + "video_id": case.video_id, + "segment_cache_key": key, + "registration": { + "contrast_side": prepared.contrast.side, + "contrast_to_target_offset": list(prepared.offset), + "left_to_right_correlation": list(registration.correlation), + "left_to_right_template": list(registration.template), + "disagreement_px": registration.disagreement_px, + }, + **quant_summary(result), + } (case_dir / "summary.json").write_text( json.dumps(summary, ensure_ascii=False, indent=2, allow_nan=False) + "\n", encoding="utf-8" ) @@ -214,7 +313,7 @@ def main(argv: Sequence[str]) -> None: clock = time.monotonic() written = load_config(config_path) - cfg = _resolved(written) + cfg = resolve_paths(written) git = read_git_state(REPO_ROOT) run_id = make_build_id(started_at, git) run_dir = cfg.paths.runs_root / run_id @@ -224,6 +323,7 @@ def main(argv: Sequence[str]) -> None: raise FileExistsError(f"运行目录已存在，不覆盖：{path}") partial.mkdir(parents=True, exist_ok=False) + review = read_review(cfg.paths.align_review) delivery = read_delivery(cfg.paths.delivery_dir, None) checkpoint_sha256 = file_sha256(cfg.segment.checkpoint) commit = medsam2_commit(cfg.segment.medsam2_dir) @@ -238,7 +338,7 @@ def main(argv: Sequence[str]) -> None: total = len(delivery.cases) for number, case in enumerate(delivery.cases, start=1): try: - slots, hit = _run_case(case, cfg, key_of, partial) + slots, hit = _run_case(case, cfg, key_of, review, partial) except CeusCaseError as error: failures.append((case.video_id, error.step, error.reason)) print(f"[{number}/{total}] 失败在第 {error.step} 步", flush=True) @@ -259,7 +359,7 @@ def main(argv: Sequence[str]) -> None: started_at=started_at, finished_at=datetime.now().astimezone(), uv_lock=REPO_ROOT / "uv.lock", - configs={"ceus.yaml": config_path}, + configs={"ceus.yaml": config_path, "align_review.csv": cfg.paths.align_review}, ) segment_record = { "checkpoint_sha256": checkpoint_sha256, diff --git a/src/cadx/ceus/__init__.py b/src/cadx/ceus/__init__.py index 2291d5f..3317642 100644 --- a/src/cadx/ceus/__init__.py +++ b/src/cadx/ceus/__init__.py @@ -1,8 +1,8 @@ """ceus：CEUS 定量的全部算法和参数，纯计算，不引用其他 cadx 模块。 -对外两段接口：`prepare_segmentation` 做第 1–2 步（解码、分屏、定造影那一半、选提示帧），交出分割的 -输入；`quantify_case` 拿分割结果做第 4–9 步（轮廓、暗区、背景圆、TIC、拟合、关键帧、灌注指标）。 -第 3 步分割由 segment 做，ceus 不引用它。某一例不满足算法前提时抛 `CeusCaseError` +对外两段接口：`prepare_segmentation` 做第 1–2 步（解码、分屏、配准两半、定造影那一半、选提示帧）， +交出分割的输入；`quantify_case` 拿分割结果做第 4–9 步（轮廓、暗区、背景圆、TIC、拟合、关键帧、 +灌注指标）。第 3 步分割由 segment 做，ceus 不引用它。某一例不满足算法前提时抛 `CeusCaseError` （带失败的步号和原因）。 写产物用：`pack_trajectory`、`pack_tic_curves`、`render_tic_png` + `tic_plot_markers`、 `render_case_overlay`、`keyframe_images`、`quant_summary`。 @@ -30,6 +30,7 @@ from cadx.ceus._config import ( PanSourceConfig, PromptSourceConfig, ReferenceRegionConfig, + RegistrationConfig, RoiShapeConfig, SmoothingConfig, TicAnalysisConfig, @@ -70,7 +71,7 @@ from cadx.ceus._green_prompt import ( extract_color_prompt, prompt_color_mask, ) -from cadx.ceus._mapping import map_source_mask_to_target, target_polygons_to_source, view_offset +from cadx.ceus._mapping import Offset, map_source_mask_to_target, target_polygons_to_source from cadx.ceus._overlay import OverlayConfig, OverlayLayers, draw_polygon_overlay, render_overlay from cadx.ceus._perfusion_metrics import ( LognormalPerfusionMetrics, @@ -98,6 +99,14 @@ from cadx.ceus._polygons import ( ) from cadx.ceus._prompt_frames import PromptFrame, select_prompt_frames from cadx.ceus._reference_region import select_reference_region +from cadx.ceus._registration import ( + HalfRegistration, + RegistrationError, + TemplatePatch, + color_marks, + fan_mask, + register_halves, +) from cadx.ceus._split_screen import ( SplitDetectionError, detect_split_x, @@ -159,6 +168,7 @@ __all__ = [ "ConstraintsConfig", "ContrastSide", "ContrastSideError", + "HalfRegistration", "InteriorDarkRegionConfig", "KeyframePositionConstraint", "KeyframesConfig", @@ -166,6 +176,7 @@ __all__ = [ "LognormalParameters", "LognormalPerfusionMetrics", "NoValidContourError", + "Offset", "OverlayConfig", "OverlayLayers", "PanSourceConfig", @@ -174,10 +185,13 @@ __all__ = [ "PromptSourceConfig", "QuantResult", "ReferenceRegionConfig", + "RegistrationConfig", + "RegistrationError", "RoiShapeConfig", "SegmentationInput", "SmoothingConfig", "SplitDetectionError", + "TemplatePatch", "TicAnalysisConfig", "TicBaselineConfig", "TicCurveSeries", @@ -194,6 +208,7 @@ __all__ = [ "apply_dark_morphology", "apply_prompt_frame_position_constraints", "collect_prompt_observations", + "color_marks", "compute_lognormal_perfusion_metrics", "compute_tic_for_video", "crop_frame", @@ -206,6 +221,7 @@ __all__ = [ "extract_color_prompt", "extract_dark_region_masks", "extract_dark_region_polygons", + "fan_mask", "fill_binary_holes", "fit_lognormal_curve", "from_tic_result", @@ -233,6 +249,7 @@ __all__ = [ "quant_summary", "quantify_case", "rect_from_values", + "register_halves", "render_case_overlay", "render_overlay", "render_tic_png", @@ -248,6 +265,5 @@ __all__ = [ "triangular_smooth", "unpack_tic_curves", "unpack_trajectory", - "view_offset", "view_rect", ] diff --git a/src/cadx/ceus/_config.py b/src/cadx/ceus/_config.py index c112f60..fcb893d 100644 --- a/src/cadx/ceus/_config.py +++ b/src/cadx/ceus/_config.py @@ -7,8 +7,8 @@ null；多写或拼错的键报错；读入后不能修改（`AGENTS.md`「不 旧配置里以下字段删掉了，不留着不用： - `views.source_view` / `views.target_view`：造影画面在哪一半改由两半各自找绿框决定（`_views.py`）； - `views.effective_region`：裁画面黑边只有旧 `service.py` 用，40 例结果走的 debug 路径不裁； -- `mapping`（含锚点）：两半之间的坐标换算改成两个视图矩形左上角之差（`_mapping.py`）， - 不再需要按某一种视频宽度写死的锚点。 +- `mapping`（含锚点）：两半之间的坐标换算改成每段视频单独配准（`registration` 段、 + `_registration.py`），不再需要按某一种视频宽度写死的锚点。 `paths` 段和整份文件最外层的模型由 `scripts/ceus_quant.py` 定义，不在这里。 """ @@ -58,6 +58,49 @@ class ViewsConfig(BaseModel): return value +class RegistrationConfig(BaseModel): + """两半配准（`_registration.py`）：每段视频单独估左半到右半的平移。 + + - `sample_frames`：均匀抽几帧求时间平均和时间标准差； + - `top_ignore_ratio`：顶部信息栏占高度的比例，这几行不算扇形区； + - `fan_std_threshold`：时间标准差（8 位灰度）高于它的像素算扇形区； + - `fan_open_kernel`：扇形掩膜开运算的核边长（像素，正奇数）； + - `mark_dilate`：绿框、紫圈颜色像素修补前外扩的像素数； + - `blur_sigma`：边缘图和模板灰度共用的高斯平滑 σ（像素）； + - `search_x_ratio`、`search_y_ratio`：搜索范围，平移各分量的绝对值上限占左半宽、高的比例； + - `template_grid`：模板块网格 `[行, 列]`；`template_size_ratio`：模板块边长占左半宽度的比例； + - `template_min_patches`：整块落在扇形内、可用的模板块至少几块； + - `max_disagreement_px`：两种方法估计各分量之差的最大值的门槛（像素），超过这一例失败。 + """ + + model_config = _STRICT + + sample_frames: int = Field(ge=2) + top_ignore_ratio: float = Field(ge=0.0, lt=1.0) + fan_std_threshold: float = Field(gt=0.0) + fan_open_kernel: int + mark_dilate: int = Field(ge=0) + blur_sigma: float = Field(gt=0.0) + search_x_ratio: float = Field(gt=0.0, le=1.0) + search_y_ratio: float = Field(gt=0.0, le=1.0) + template_grid: list[int] + template_size_ratio: float = Field(gt=0.0, lt=1.0) + template_min_patches: int = Field(ge=1) + max_disagreement_px: int = Field(ge=0) + + @field_validator("fan_open_kernel") + @classmethod + def _open_kernel_odd(cls, value: int) -> int: + return _require_positive_odd(value, "registration.fan_open_kernel") + + @field_validator("template_grid") + @classmethod + def _grid_pair(cls, value: list[int]) -> list[int]: + if len(value) != 2 or min(value) < 1: + raise ValueError(f"registration.template_grid 必须是两个正整数 [行, 列]，收到 {value}") + return value + + class PromptSourceConfig(BaseModel): """绿框的 HSV 颜色检测。`min_area` 是候选外接框面积 w×h，不是轮廓面积。""" @@ -548,6 +591,7 @@ class CeusConfig(BaseModel): config_version: str = Field(min_length=1) random_seed: int views: ViewsConfig + registration: RegistrationConfig prompt_source: PromptSourceConfig pan_source: PanSourceConfig keyframes: KeyframesConfig diff --git a/src/cadx/ceus/_mapping.py b/src/cadx/ceus/_mapping.py index 4ceaf36..107c5ef 100644 --- a/src/cadx/ceus/_mapping.py +++ b/src/cadx/ceus/_mapping.py @@ -1,63 +1,38 @@ -"""两半画面之间的坐标换算。 +"""两半画面之间的坐标换算：造影那一半局部坐标 + 平移 `(dx, dy)` = 另一半局部坐标。 绿框画在造影那一半；分割在另一半（灰阶画面，病灶边界不随造影剂进出明暗变化）上做，轮廓再挪回 -造影那一半采亮度。两半是同一个探头视野的两路成像，大小相同、只差一个水平平移，所以换算只做平移： +造影那一半采亮度。两半是同一个探头视野的两路成像，大小相同，只差一个平移；平移量由第 1 步 +每段视频单独配准得到（`_registration.py`），这里只做加减和裁边，不缩放、不旋转。 - 平移量 = 目标那一半视图的左上角 − 造影那一半视图的左上角（整帧像素坐标） - -两个视图矩形直接来自第 1 步的分屏结果，不依赖视频宽度。旧代码改用按 1856 像素宽的视频写死的 -锚点，其余宽度的视频会挪偏一两百像素（design 0002 决定 4），这里不再用锚点。 - -两个函数都在各自视图的局部坐标里收发掩膜和多边形：局部 → 整帧（加本视图左上角）→ 平移 → -目标视图局部（减目标视图左上角）。两半宽度若差几个像素，超出目标视图的部分裁掉，不缩放。 +两个函数都在各自视图的局部坐标里收发掩膜和多边形；超出目标视图的部分裁掉。 """ import numpy as np from cadx.ceus._geometry import AnyArray, Mask, Points, ViewRect - -def view_offset(source_view_rect: ViewRect, target_view_rect: ViewRect) -> tuple[int, int]: - """从造影那一半到目标那一半的平移量 `(dx, dy)`：两个视图左上角之差。""" - return ( - target_view_rect.x1 - source_view_rect.x1, - target_view_rect.y1 - source_view_rect.y1, - ) +type Offset = tuple[int, int] def map_source_mask_to_target( - source_mask: AnyArray, - source_view_rect: ViewRect, - target_view_rect: ViewRect, + source_mask: AnyArray, offset: Offset, target_view_rect: ViewRect ) -> Mask: - """把造影那一半局部坐标的掩膜挪到目标那一半的局部坐标，形状为目标视图大小。""" - dx, dy = view_offset(source_view_rect, target_view_rect) - # 源局部 (x, y) → 整帧 (x + sx1, y + sy1) → 平移 → 目标局部 (x + sx1 + dx − tx1, …)。 - shift_x = source_view_rect.x1 + dx - target_view_rect.x1 - shift_y = source_view_rect.y1 + dy - target_view_rect.y1 + """把造影那一半局部坐标的掩膜平移 `offset`，形状为目标视图大小，出界的像素丢掉。""" + dx, dy = offset mapped = np.zeros((target_view_rect.height, target_view_rect.width), dtype=bool) ys, xs = np.nonzero(source_mask) - xs = xs + shift_x - ys = ys + shift_y + xs = xs + dx + ys = ys + dy inside = (xs >= 0) & (xs < target_view_rect.width) & (ys >= 0) & (ys < target_view_rect.height) mapped[ys[inside], xs[inside]] = True return mapped def target_polygons_to_source( - target_polygons: list[AnyArray], - source_view_rect: ViewRect, - target_view_rect: ViewRect, + target_polygons: list[AnyArray], offset: Offset, source_view_rect: ViewRect ) -> list[Points]: - """把目标那一半局部坐标的多边形逐个挪回造影那一半的局部坐标，裁到视图边界内。""" - dx, dy = view_offset(source_view_rect, target_view_rect) - shift = np.array( - [ - target_view_rect.x1 - dx - source_view_rect.x1, - target_view_rect.y1 - dy - source_view_rect.y1, - ], - dtype=np.float32, - ) + """把目标那一半局部坐标的多边形逐个减去 `offset` 挪回造影那一半，裁到视图边界内。""" + shift = np.array([-offset[0], -offset[1]], dtype=np.float32) mapped: list[Points] = [] for polygon in target_polygons: source_poly = polygon.astype(np.float32) + shift[None, :] diff --git a/src/cadx/ceus/_pipeline.py b/src/cadx/ceus/_pipeline.py index 1ed0aff..5782255 100644 --- a/src/cadx/ceus/_pipeline.py +++ b/src/cadx/ceus/_pipeline.py @@ -1,7 +1,8 @@ """ceus 对外的两段接口，中间由 `scripts/ceus_quant.py` 调 segment 做第 3 步分割。 -- `prepare_segmentation`（第 1–2 步）：读视频、找分屏线、判断造影画面在哪一半、选提示帧，交出 - 要交给分割的那一半（灰阶那一半，`target_view`）的逐帧图像和提示帧掩膜（`target_view` 局部坐标）。 +- `prepare_segmentation`（第 1–2 步）：读视频、找分屏线、配准两半、判断造影画面在哪一半、选提示帧， + 交出要交给分割的那一半（灰阶那一半，`target_view`）的逐帧图像和提示帧掩膜（`target_view` + 局部坐标）。 - `quantify_case`（第 4–9 步）：拿分割得到的逐帧掩膜（`target_view` 局部坐标）整理成平滑轮廓并与 绿框融合、挪回造影那一半，找暗区、选背景圆、逐帧采 TIC、拟合、定位 7 个关键帧、算灌注指标。 @@ -33,7 +34,7 @@ from cadx.ceus._config import CeusConfig from cadx.ceus._contrast_side import ContrastSide, locate_contrast_side from cadx.ceus._dark_region import extract_dark_region_polygons from cadx.ceus._geometry import AnyArray, Frame, Mask, crop_frame -from cadx.ceus._mapping import map_source_mask_to_target, target_polygons_to_source +from cadx.ceus._mapping import Offset, map_source_mask_to_target, target_polygons_to_source from cadx.ceus._overlay import OverlayConfig, OverlayLayers, render_overlay from cadx.ceus._perfusion_metrics import ( LognormalPerfusionMetrics, @@ -45,6 +46,7 @@ from cadx.ceus._polygons import ( masks_to_smoothed_polygons, ) from cadx.ceus._reference_region import select_reference_region +from cadx.ceus._registration import HalfRegistration, register_halves from cadx.ceus._split_screen import split_views from cadx.ceus._tic import CURVE_NAMES, TicResult, compute_tic_for_video from cadx.ceus._tic_curves import TicCurves, from_tic_result @@ -80,8 +82,10 @@ class SegmentationInput: """第一段（第 1–2 步）的结果。 - `video_path`、`fps`、`frame_count`、`frame_size`（宽, 高）、`frame_times_s`：视频本身； - - `split_x`：分屏线的列；`contrast`：哪一半是造影画面、两个视图矩形（整帧坐标）、提示帧 - （造影那一半局部坐标的绿框掩膜）； + - `split_x`：分屏线的列；`registration`：两半配准（左半 → 右半，两种方法各自的估计）； + `offset`：造影那一半局部坐标 → `target_view` 局部坐标的平移 `(dx, dy)`，取互相关的估计； + - `contrast`：哪一半是造影画面、两个视图矩形（整帧坐标）、提示帧（造影那一半局部坐标的 + 绿框掩膜）； - `target_frames`：要交给分割的那一半（`contrast.target_view`，灰阶画面）的逐帧 BGR 图像， 与帧号一一对应； - `prompt_masks`：提示帧帧号 → 这一帧绿框挪到 `target_view` 局部坐标后的布尔掩膜； @@ -94,6 +98,8 @@ class SegmentationInput: frame_size: tuple[int, int] frame_times_s: AnyArray split_x: int + registration: HalfRegistration + offset: Offset contrast: ContrastSide target_frames: list[Frame] prompt_masks: dict[int, Mask] @@ -124,11 +130,11 @@ class QuantResult: def prepare_segmentation(video_path: Path, cfg: CeusConfig) -> SegmentationInput: - """第 1–2 步：解码、分屏、找绿框定造影那一半、选提示帧，交出分割的输入。 + """第 1–2 步：解码、分屏、配准两半、找绿框定造影那一半、选提示帧，交出分割的输入。 - 视频读不了或读出的不合规抛 `VideoDecodeError`、找不到可信分屏线抛 `SplitDetectionError`（都是 - 第 1 步），两半都有或都没有绿框抛 `ContrastSideError` - （第 2 步），两者都是 `CeusCaseError`。 + 视频读不了或读出的不合规抛 `VideoDecodeError`、找不到可信分屏线抛 `SplitDetectionError`、两半 + 配准失败（含两种方法不一致）抛 `RegistrationError`（都是第 1 步），两半都有或都没有绿框抛 + `ContrastSideError`（第 2 步），都是 `CeusCaseError`。 """ video = VideoFile(video_path) fps, frame_count, width, height = video.probe() @@ -140,15 +146,17 @@ def prepare_segmentation(video_path: Path, cfg: CeusConfig) -> SegmentationInput ) split_x, left_view, right_view = split_views(full_frames, cfg.views) + registration = register_halves( + full_frames, left_view, right_view, cfg.registration, cfg.prompt_source, cfg.pan_source + ) contrast = locate_contrast_side( full_frames, left_view, right_view, cfg.prompt_source, cfg.keyframes ) + offset = registration.contrast_to_target(contrast.side) target_frames = [crop_frame(frame, contrast.target_view).copy() for frame in full_frames] first_contrast_frame = crop_frame(full_frames[0], contrast.contrast_view).copy() prompt_masks = { - prompt.frame_index: map_source_mask_to_target( - prompt.mask, contrast.contrast_view, contrast.target_view - ) + prompt.frame_index: map_source_mask_to_target(prompt.mask, offset, contrast.target_view) for prompt in contrast.prompt_frames } return SegmentationInput( @@ -158,6 +166,8 @@ def prepare_segmentation(video_path: Path, cfg: CeusConfig) -> SegmentationInput frame_size=(width, height), frame_times_s=frame_times_s, split_x=split_x, + registration=registration, + offset=offset, contrast=contrast, target_frames=target_frames, prompt_masks=prompt_masks, @@ -212,20 +222,20 @@ def quantify_case( target_polygons = apply_prompt_frame_position_constraints( target_polygons, prompt_frame_dicts, - contrast.contrast_view, + prepared.offset, contrast.target_view, cfg.constraints.keyframe_position, ) target_polygons = fuse_prompt_prior_with_medsam( target_polygons, prompt_frame_dicts, - contrast.contrast_view, + prepared.offset, contrast.target_view, cfg.roi_shape, cfg.smoothing, ) source_polygons = target_polygons_to_source( - target_polygons, contrast.contrast_view, contrast.target_view + target_polygons, prepared.offset, contrast.contrast_view ) # 第 5 步：暗区（在造影那一半上找）。 @@ -296,6 +306,7 @@ def quantify_case( "split_x": prepared.split_x, "source": "auto_detected", "contrast_side": contrast.side, + "contrast_to_target_offset": list(prepared.offset), }, ) curves = from_tic_result( diff --git a/src/cadx/ceus/_polygons.py b/src/cadx/ceus/_polygons.py index c45906f..102ec40 100644 --- a/src/cadx/ceus/_polygons.py +++ b/src/cadx/ceus/_polygons.py @@ -13,8 +13,9 @@ 其余帧用 `propagation_shape_strength`），再把面积收紧到先验面积的 `[min, max]` 倍 （绿框在最终轮廓里占多少，见 `待定问题/01`）。 -照搬旧项目 `polygons.py`，去掉产品层；坐标换算改用两个视图矩形左上角之差（`_mapping.py`）， -不再需要旧的有效区矩形和映射锚点。 +照搬旧项目 `polygons.py`，去掉产品层；绿框挪到 target_view 用第 1 步配准得到的平移 `offset` +（造影那一半局部坐标 + offset = target_view 局部坐标，`_mapping.py`），不再需要旧的有效区矩形和 +映射锚点。 """ from typing import Any, Protocol @@ -32,7 +33,7 @@ from cadx.ceus._geometry import ( resample_closed_contour, triangular_smooth, ) -from cadx.ceus._mapping import map_source_mask_to_target +from cadx.ceus._mapping import Offset, map_source_mask_to_target class _SmoothingConfigLike(Protocol): @@ -205,15 +206,13 @@ def masks_to_smoothed_polygons( def _prompt_frame_target_centers( prompt_frames: list[dict[str, Any]], - source_view_rect: ViewRect, + offset: Offset, target_view_rect: ViewRect, ) -> dict[int, AnyArray]: """各提示帧的绿框掩膜挪到 target_view 后的质心，`帧号 → (cx, cy)`。""" centers: dict[int, AnyArray] = {} for prompt_frame in prompt_frames: - target_mask = map_source_mask_to_target( - prompt_frame["mask"], source_view_rect, target_view_rect - ) + target_mask = map_source_mask_to_target(prompt_frame["mask"], offset, target_view_rect) centers[int(prompt_frame["frame_index"])] = mask_center(target_mask) return centers @@ -221,7 +220,7 @@ def _prompt_frame_target_centers( def apply_prompt_frame_position_constraints( polygons: list[AnyArray], prompt_frames: list[dict[str, Any]], - source_view_rect: ViewRect, + offset: Offset, target_view_rect: ViewRect, constraint_cfg: _KeyframePositionConstraintLike, ) -> list[AnyArray]: @@ -241,7 +240,7 @@ def apply_prompt_frame_position_constraints( ) strength = float(np.clip(constraint_cfg.strength, 0.0, 1.0)) - anchor_centers = _prompt_frame_target_centers(prompt_frames, source_view_rect, target_view_rect) + anchor_centers = _prompt_frame_target_centers(prompt_frames, offset, target_view_rect) center_track = _interpolate_centers(anchor_centers, len(polygons)) corrected: list[AnyArray] = [] @@ -254,7 +253,7 @@ def apply_prompt_frame_position_constraints( def _prompt_prior_polygons( prompt_frames: list[dict[str, Any]], frame_count: int, - source_view_rect: ViewRect, + offset: Offset, target_view_rect: ViewRect, smoothing_cfg: _SmoothingConfigLike, ) -> list[AnyArray]: @@ -270,9 +269,7 @@ def _prompt_prior_polygons( raw_key_polygons: dict[int, AnyArray] = {} previous: AnyArray | None = None for prompt_frame in prompt_frames: - target_mask = map_source_mask_to_target( - prompt_frame["mask"], source_view_rect, target_view_rect - ) + target_mask = map_source_mask_to_target(prompt_frame["mask"], offset, target_view_rect) resampled = resample_closed_contour(largest_contour(target_mask), num_points) if resampled is None: continue @@ -316,7 +313,7 @@ def _tighten_area_to_prior_band( def fuse_prompt_prior_with_medsam( medsam_polygons: list[AnyArray], prompt_frames: list[dict[str, Any]], - source_view_rect: ViewRect, + offset: Offset, target_view_rect: ViewRect, shape_cfg: _RoiShapeConfigLike, smoothing_cfg: _SmoothingConfigLike, @@ -332,7 +329,7 @@ def fuse_prompt_prior_with_medsam( prior_polygons = _prompt_prior_polygons( prompt_frames, len(medsam_polygons), - source_view_rect, + offset, target_view_rect, smoothing_cfg, ) diff --git a/src/cadx/ceus/_registration.py b/src/cadx/ceus/_registration.py new file mode 100644 index 0000000..3f4a651 --- /dev/null +++ b/src/cadx/ceus/_registration.py @@ -0,0 +1,277 @@ +"""第 1 步的两半配准：估出左半局部坐标到右半局部坐标的平移 `(dx, dy)`。 + +两半是同一个探头视野的两路成像（造影和灰阶），大小相同，但扇形在各自那一半里的位置不同： +造影那一半的扇形左边常有一段黑边，灰阶那一半的扇形贴着分屏线，而且同一帧宽的两段视频错位也 +可能不同。所以每段视频单独配准，不写死、也不用两个视图左上角相减（任务书 02 第一阶段第 1 条）。 + +做法：均匀抽若干帧，两半各自求逐像素的时间平均灰度和时间标准差。 + +1. 扇形区：时间标准差高于阈值的像素（画面在动的地方；文字、色标、刻度这些静止的界面元素 + 标准差近 0），去掉顶部信息栏那几行，开运算断开细连接，取最大连通域的凸包。 +2. 绿框、紫圈的颜色（`prompt_source`、`pan_source` 的 HSV 范围）在任一抽样帧出现过的像素 + 膨胀后用周围像素修补掉：它们只画在一半上，会被当成边缘或纹理去匹配。 +3. 方法一（互相关，结果取它）：两半的边缘图（高斯平滑后的梯度幅值，只留扇形区和它的边）做 + 归一化互相关，在搜索范围内穷举平移取最大。扇形的轮廓和组织边界都参与。 +4. 方法二（模板匹配，只用来核对）：在左半扇形内按网格定若干块，每一抽样帧上各切一次平滑后的 + 灰度，在右半搜索范围内、块中心落在右半扇形内的位置里找最相似的（归一化相关系数），全部 + 匹配的平移各分量取中位数。它用的是逐帧灰度纹理、局部块和中位数，不用边缘和扇形轮廓，和 + 方法一相互独立。 +5. 两种估计各分量之差的最大值（切比雪夫距离：水平、竖直各自都不许错开超过门槛）超过门槛， + 这一例第 1 步失败（`RegistrationError`）。 + +只做平移，不做缩放和旋转。 +""" + +from collections.abc import Sequence +from dataclasses import dataclass +from typing import Literal + +import cv2 +import numpy as np +import numpy.typing as npt + +from cadx.ceus._config import PanSourceConfig, PromptSourceConfig, RegistrationConfig +from cadx.ceus._errors import CeusCaseError +from cadx.ceus._geometry import Frame, Mask, ViewRect, crop_frame +from cadx.ceus._mapping import Offset +from cadx.ceus._split_screen import sample_frames_for_split + +REGISTRATION_STEP = 1 + +type _Gray = npt.NDArray[np.float32] + + +class RegistrationError(CeusCaseError): + """两半配准失败：找不到扇形区、可用的模板块不够，或两种方法估出的平移不一致。""" + + def __init__(self, reason: str) -> None: + super().__init__(REGISTRATION_STEP, reason) + + +@dataclass(frozen=True) +class TemplatePatch: + """方法二的一次匹配：模板块在左半局部坐标的左上角、第几张抽样帧、在右半找到的平移、相关系数。""" + + x: int + y: int + frame: int + offset: Offset + score: float + + +@dataclass(frozen=True) +class HalfRegistration: + """左半局部坐标 → 右半局部坐标的平移：`右半局部 = 左半局部 + (dx, dy)`。 + + `correlation` 是方法一（边缘图互相关）的估计，采用它；`correlation_score` 是它的峰值； + `template` 是方法二（多块模板逐帧匹配取中位数）的估计，`patches` 是每次匹配的明细； + `disagreement_px` 是两者各分量之差的最大值。 + """ + + correlation: Offset + correlation_score: float + template: Offset + patches: tuple[TemplatePatch, ...] + disagreement_px: int + + def contrast_to_target(self, contrast_side: Literal["left", "right"]) -> Offset: + """造影那一半局部坐标 → 另一半局部坐标的平移。""" + dx, dy = self.correlation + if contrast_side == "left": + return (dx, dy) + return (-dx, -dy) + + +@dataclass(frozen=True) +class _Half: + """一半画面的配准素材：扇形区掩膜，修补绿框紫圈并平滑后的时间平均灰度和逐帧灰度。""" + + fan: Mask + mean: _Gray + frames: list[_Gray] + + +def color_marks( + frames: Sequence[Frame], prompt_cfg: PromptSourceConfig, pan_cfg: PanSourceConfig, dilate: int +) -> Mask: + """任一帧里落在绿框或紫圈 HSV 范围内的像素，膨胀 `dilate` 像素。""" + marks = np.zeros(frames[0].shape[:2], dtype=np.uint8) + for frame in frames: + hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV) + for cfg in (prompt_cfg, pan_cfg): + lower = np.asarray(cfg.hsv_lower, dtype=np.uint8) + upper = np.asarray(cfg.hsv_upper, dtype=np.uint8) + marks |= np.asarray(cv2.inRange(hsv, lower, upper), dtype=np.uint8) + if dilate > 0: + kernel = np.ones((2 * dilate + 1, 2 * dilate + 1), dtype=np.uint8) + marks = np.asarray(cv2.dilate(marks, kernel), dtype=np.uint8) + return marks > 0 + + +def fan_mask(std: _Gray, cfg: RegistrationConfig) -> Mask: + """扇形区：时间标准差高于阈值、去掉顶部信息栏、开运算后最大连通域的凸包。""" + moving = (std > cfg.fan_std_threshold).astype(np.uint8) + moving[: round(cfg.top_ignore_ratio * std.shape[0])] = 0 + size = cfg.fan_open_kernel + moving = cv2.morphologyEx(moving, cv2.MORPH_OPEN, np.ones((size, size), dtype=np.uint8)) + count, labels, stats, _ = cv2.connectedComponentsWithStats(moving, connectivity=8) + if count < 2: + raise RegistrationError("这一半画面里找不到扇形区（没有随时间变化的区域）") + largest = 1 + int(np.argmax(stats[1:, cv2.CC_STAT_AREA])) + component = np.asarray(labels, dtype=np.int32) == largest + points = np.argwhere(component)[:, ::-1].astype(np.int32) # (x, y) + hull = np.asarray(cv2.convexHull(points), dtype=np.int32) + filled = np.zeros(std.shape, dtype=np.uint8) + cv2.fillConvexPoly(filled, hull, 1) + return filled > 0 + + +def _repair_smooth(gray: _Gray, marks: Mask, sigma: float) -> _Gray: + """绿框紫圈像素用周围像素修补（Telea 法，半径 3），再高斯平滑。""" + repaired = cv2.inpaint( + np.clip(gray, 0, 255).astype(np.uint8), marks.astype(np.uint8), 3, cv2.INPAINT_TELEA + ) + return np.asarray(cv2.GaussianBlur(repaired.astype(np.float32), (0, 0), sigma), np.float32) + + +def _prepare_half(gray: _Gray, view: ViewRect, marks: Mask, cfg: RegistrationConfig) -> _Half: + """从整帧灰度栈 `(N, H, W)` 切出一半，求扇形区、平均灰度和逐帧灰度。""" + stack = gray[:, view.y1 : view.y2, view.x1 : view.x2] + local_marks = crop_frame(marks, view) + return _Half( + fan=fan_mask(stack.std(axis=0), cfg), + mean=_repair_smooth(stack.mean(axis=0), local_marks, cfg.blur_sigma), + frames=[_repair_smooth(frame, local_marks, cfg.blur_sigma) for frame in stack], + ) + + +def _edge_map(half: _Half, sigma: float) -> _Gray: + """时间平均灰度的梯度幅值，只留扇形区和它的边（扇形掩膜外扩 3σ，轮廓上的梯度也算进来）。""" + gx = cv2.Sobel(half.mean, cv2.CV_32F, 1, 0) + gy = cv2.Sobel(half.mean, cv2.CV_32F, 0, 1) + margin = int(np.ceil(3 * sigma)) + kernel = np.ones((2 * margin + 1, 2 * margin + 1), dtype=np.uint8) + support = cv2.dilate(half.fan.astype(np.uint8), kernel) > 0 + edges = np.hypot(gx, gy).astype(np.float32) + edges[~support] = 0 + return edges + + +def _bbox(mask: Mask) -> tuple[int, int, int, int]: + """掩膜的外接框 `(x1, y1, x2, y2)`，右下不含。""" + ys, xs = np.nonzero(mask) + return int(xs.min()), int(ys.min()), int(xs.max()) + 1, int(ys.max()) + 1 + + +def _correlate_edges( + left: _Half, right: _Half, search: Offset, sigma: float +) -> tuple[Offset, float]: + """方法一：左半边缘图（裁到扇形外接框）在补零后的右半边缘图上做归一化互相关。""" + sx, sy = search + left_edges = _edge_map(left, sigma) + right_edges = _edge_map(right, sigma) + x1, y1, x2, y2 = _bbox(left_edges > 0) + template = left_edges[y1:y2, x1:x2] + height, width = left_edges.shape + # 右半四周补零，使平移 dx ∈ [−sx, sx]、dy ∈ [−sy, sy] 都落在结果图里；右、下多补一个左半大小， + # 右半比左半窄几个像素时也够。 + padded = cv2.copyMakeBorder( + right_edges, sy, sy + height, sx, sx + width, cv2.BORDER_CONSTANT, value=0.0 + ) + scores = cv2.matchTemplate(padded, template, cv2.TM_CCORR_NORMED) + # 结果图里位置 (u, v) 对应平移 dx = u − sx − x1、dy = v − sy − y1。 + window = scores[y1 : y1 + 2 * sy + 1, x1 : x1 + 2 * sx + 1] + _, best, _, (u, v) = cv2.minMaxLoc(window) + return (u - sx, v - sy), float(best) + + +def _match_templates( + left: _Half, right: _Half, search: Offset, cfg: RegistrationConfig +) -> tuple[TemplatePatch, ...]: + """方法二：左半扇形外接框上按网格切块，整块都在扇形内的才用；每块在每一抽样帧上各找一次。 + + 逐帧匹配而不是在时间平均图上匹配：探头在视频里移动时，平均图上的组织纹理糊掉了。 + """ + sx, sy = search + rows, cols = cfg.template_grid + side = round(cfg.template_size_ratio * left.fan.shape[1]) + x1, y1, x2, y2 = _bbox(left.fan) + right_h, right_w = right.fan.shape + patches: list[TemplatePatch] = [] + for row in range(rows): + for col in range(cols): + px = round(x1 + (col + 0.5) * (x2 - x1) / cols - side / 2) + py = round(y1 + (row + 0.5) * (y2 - y1) / rows - side / 2) + inside_view = px >= 0 and py >= 0 and px + side <= left.fan.shape[1] + inside_view = inside_view and py + side <= left.fan.shape[0] + if not inside_view or not left.fan[py : py + side, px : px + side].all(): + continue + wx1, wy1 = max(0, px - sx), max(0, py - sy) + wx2, wy2 = min(right_w, px + side + sx), min(right_h, py + side + sy) + if wx2 - wx1 < side or wy2 - wy1 < side: + continue + # 只认块中心落在右半扇形内的位置：扇形外是纯黑，零方差窗口的相关系数 OpenCV 记为 1。 + half = side // 2 + centers = right.fan[ + wy1 + half : wy2 - side + half + 1, wx1 + half : wx2 - side + half + 1 + ] + if not centers.any(): + continue + for index, (left_frame, right_frame) in enumerate( + zip(left.frames, right.frames, strict=True) + ): + template = left_frame[py : py + side, px : px + side] + scores = cv2.matchTemplate( + right_frame[wy1:wy2, wx1:wx2], template, cv2.TM_CCOEFF_NORMED + ) + scores[~centers] = -1.0 + _, best, _, (u, v) = cv2.minMaxLoc(scores) + offset = (wx1 + u - px, wy1 + v - py) + patches.append(TemplatePatch(px, py, index, offset, float(best))) + return tuple(patches) + + +def register_halves( + frames: Sequence[Frame], + left_view: ViewRect, + right_view: ViewRect, + cfg: RegistrationConfig, + prompt_cfg: PromptSourceConfig, + pan_cfg: PanSourceConfig, +) -> HalfRegistration: + """两种方法各估一次左半 → 右半的平移，相差超过门槛抛 `RegistrationError`。""" + sampled = sample_frames_for_split(frames, cfg.sample_frames) + gray = np.stack([cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY) for frame in sampled]) + gray = gray.astype(np.float32) + marks = color_marks(sampled, prompt_cfg, pan_cfg, cfg.mark_dilate) + left = _prepare_half(gray, left_view, marks, cfg) + right = _prepare_half(gray, right_view, marks, cfg) + search = ( + round(cfg.search_x_ratio * left_view.width), + round(cfg.search_y_ratio * left_view.height), + ) + + correlation, correlation_score = _correlate_edges(left, right, search, cfg.blur_sigma) + patches = _match_templates(left, right, search, cfg) + blocks = len({(patch.x, patch.y) for patch in patches}) + if blocks < cfg.template_min_patches: + raise RegistrationError( + f"左半扇形内可用的模板块只有 {blocks} 块，少于 {cfg.template_min_patches} 块" + ) + template = ( + round(float(np.median([patch.offset[0] for patch in patches]))), + round(float(np.median([patch.offset[1] for patch in patches]))), + ) + disagreement = max(abs(correlation[0] - template[0]), abs(correlation[1] - template[1])) + if disagreement > cfg.max_disagreement_px: + raise RegistrationError( + f"两半配准不一致：互相关估计 (dx, dy) = {correlation}，模板匹配估计 {template}，" + f"相差 {disagreement} 像素，超过门槛 {cfg.max_disagreement_px}" + ) + return HalfRegistration( + correlation=correlation, + correlation_score=correlation_score, + template=template, + patches=patches, + disagreement_px=disagreement, + ) diff --git a/tests/fake_ceus_video.py b/tests/fake_ceus_video.py new file mode 100644 index 0000000..7ab77eb --- /dev/null +++ b/tests/fake_ceus_video.py @@ -0,0 +1,57 @@ +"""测试用合成 CEUS 视频：400×200、10 fps、40 秒的分屏造影视频（test_ceus_pipeline、test_ceus_quant、 +test_ceus_align_review 共用）。 + +两半各有一块位置相同、带纹理、随时间轻微起伏的扇形区（`fake_split_scene`，配准平移 (0, 0)）。 +造影画面在右半：病灶是整帧 `LESION` 方块，亮度按 Lognormal 团注曲线 `BOLUS` 变化，内部有一块一直 +发黑的囊 `CYST`；`green_box` 为真时画医生的绿框（病灶的空心描边）。 +""" + +from pathlib import Path + +import cv2 +import numpy as np + +from cadx.ceus import LognormalParameters, evaluate_lognormal +from fake_split_scene import make_scene, scene_gray, to_bgr + +WIDTH, HEIGHT, SPLIT = 400, 200, 200 +FPS, FRAMES = 10.0, 400 +LESION = (250, 60, 369, 159) # 整帧坐标 x1, y1, x2, y2（含端点） +CYST = (290, 95, 320, 125) +BOLUS = LognormalParameters(area=1300.0, mu=2.4, sigma=0.45, arrival_offset_s=4.0, baseline=40.0) + + +def write_video(path: Path, *, green_box: bool) -> None: + rng = np.random.default_rng(3) + intensity = evaluate_lognormal(np.arange(FRAMES) / FPS, BOLUS) + # 两半是同一幅带纹理的扇形，平移 (0, 0)：配准应估出 (0, 0)，绿框原位挪到左半。 + scene = make_scene( + height=HEIGHT, + half_width=SPLIT, + fan_polygon=[(30, 22), (170, 22), (197, 192), (3, 192)], + offset=(0, 0), + interior_shift=None, + texture_std=20.0, + # 造影那一半调暗（扇形内平均约 32，与改动前的 30 相当）：mp4 压缩把绿框描边的颜色渗到 + # 相邻像素，背景亮时渗出的像素更容易落进绿框的 HSV 范围，绿框会被认大。 + right_gain=0.15, + right_bias=0.0, + seed=7, + ) + x1, y1, x2, y2 = LESION + # 病灶里保留组织纹理（减去它在病灶方块内的均值），病灶平均亮度仍等于团注曲线。 + detail = scene.right[y1 : y2 + 1, x1 - SPLIT : x2 - SPLIT + 1] + detail = detail - detail.mean() + writer = cv2.VideoWriter(str(path), cv2.VideoWriter.fourcc(*"mp4v"), FPS, (WIDTH, HEIGHT)) + for index in range(FRAMES): + gray = scene_gray(scene, index, rng) + gray[:, SPLIT] = 0 # 分屏缝 + lesion = intensity[index] + detail + rng.normal(0.0, 6.0, detail.shape) + gray[y1 : y2 + 1, x1 : x2 + 1] = lesion + cx1, cy1, cx2, cy2 = CYST + gray[cy1 : cy2 + 1, cx1 : cx2 + 1] = 10 + frame = np.ascontiguousarray(to_bgr(gray)) + if green_box: + cv2.rectangle(frame, (x1, y1), (x2, y2), (0, 255, 0), thickness=3) + writer.write(frame) + writer.release() diff --git a/tests/fake_split_scene.py b/tests/fake_split_scene.py new file mode 100644 index 0000000..b0f7f01 --- /dev/null +++ b/tests/fake_split_scene.py @@ -0,0 +1,92 @@ +"""测试用合成分屏画面：左右两半各有一块带纹理的扇形区，右半内容是左半平移 `(dx, dy)` 后的样子。 + +模仿真实录像里配准要用到的几样东西：扇形区里的组织纹理随时间轻微起伏（扇形外的界面是静止的）， +两半的成像亮度不同（右半扇形内做一次线性变换 `gain × 灰度 + bias`），各帧噪声独立。 + +- `right_local(x, y) = left_local(x − dx, y − dy)`：右半局部坐标 = 左半局部坐标 + (dx, dy)； +- `interior_shift` 不为 None 时，右半扇形的轮廓仍按 `offset` 平移，里面的纹理却按 `interior_shift` + 平移，用来造「扇形轮廓和组织纹理给出不同平移」的画面。 +""" + +from dataclasses import dataclass + +import cv2 +import numpy as np +import numpy.typing as npt + +type Gray = npt.NDArray[np.float32] + +OUTSIDE = 20.0 # 扇形外静止界面的灰度 +BASE = 100.0 # 扇形内的平均灰度 + + +@dataclass(frozen=True) +class SplitScene: + """一幅分屏画面的静止部分：两半的灰度内容和扇形掩膜（各自局部坐标）。""" + + left: Gray + right: Gray + left_fan: npt.NDArray[np.bool_] + right_fan: npt.NDArray[np.bool_] + + +def _fan(shape: tuple[int, int], polygon: list[tuple[int, int]], shift: tuple[int, int]) -> Gray: + points = np.asarray(polygon, dtype=np.int32) + np.asarray(shift, dtype=np.int32)[None, :] + mask = np.zeros(shape, dtype=np.uint8) + cv2.fillConvexPoly(mask, points, 1) + return mask.astype(np.float32) + + +def make_scene( + *, + height: int, + half_width: int, + fan_polygon: list[tuple[int, int]], + offset: tuple[int, int], + interior_shift: tuple[int, int] | None, + texture_std: float, + right_gain: float, + right_bias: float, + seed: int, +) -> SplitScene: + """左半扇形取 `fan_polygon`（左半局部坐标），右半按 `offset` 平移；纹理是平滑过的随机场。""" + rng = np.random.default_rng(seed) + margin = 64 + canvas = cv2.GaussianBlur( + rng.normal(0.0, 1.0, (height + 2 * margin, half_width + 2 * margin)).astype(np.float32), + (0, 0), + 4.0, + ) + canvas *= texture_std / float(canvas.std()) + shape = (height, half_width) + left_texture = canvas[margin : margin + height, margin : margin + half_width] + tdx, tdy = offset if interior_shift is None else interior_shift + right_texture = canvas[ + margin - tdy : margin - tdy + height, margin - tdx : margin - tdx + half_width + ] + left_fan = _fan(shape, fan_polygon, (0, 0)) + right_fan = _fan(shape, fan_polygon, offset) + left = OUTSIDE + left_fan * (BASE - OUTSIDE + left_texture) + # 右半是另一种成像：同一组织，亮度做一次线性变换。 + right = OUTSIDE + right_fan * (right_gain * (BASE - OUTSIDE + right_texture) + right_bias) + return SplitScene( + left=left.astype(np.float32), + right=right.astype(np.float32), + left_fan=left_fan > 0, + right_fan=right_fan > 0, + ) + + +def scene_gray(scene: SplitScene, index: int, rng: np.random.Generator) -> Gray: + """第 `index` 帧的整帧灰度（左右拼接）：扇形内加缓慢起伏和独立噪声，扇形外不动。""" + wave = 10.0 * np.sin(2 * np.pi * index / 7.0) + halves: list[Gray] = [] + for content, fan in ((scene.left, scene.left_fan), (scene.right, scene.right_fan)): + noise = rng.normal(0.0, 3.0, content.shape).astype(np.float32) + halves.append(content + fan * (wave + noise)) + return np.hstack(halves) + + +def to_bgr(gray: Gray) -> npt.NDArray[np.uint8]: + """灰度裁到 0–255 转成三通道 uint8。""" + return np.repeat(np.clip(gray, 0, 255).astype(np.uint8)[:, :, None], 3, axis=2) diff --git a/tests/test_ceus_align_review.py b/tests/test_ceus_align_review.py new file mode 100644 index 0000000..b80c43f --- /dev/null +++ b/tests/test_ceus_align_review.py @@ -0,0 +1,196 @@ +"""`scripts/ceus_align_review.py`：配准目视核对的出图和计分（任务书 02「配准的目视核对」）。 + +出图：假交付包两例合成分屏视频（`fake_ceus_video`，两半扇形对齐，配准 (0, 0)），对照取 1 例。 +计分：key.csv 和两份判读手写，预期按任务书的规则手算。表头按任务书和协调者给的列名独立写。 +""" + +import csv +import hashlib +import importlib.util +import math +import re +from pathlib import Path +from types import ModuleType +from typing import Any + +import cv2 +import pytest +import yaml + +from conftest import CEUS_YAML +from fake_ceus_video import write_video +from fake_delivery import make_delivery + +_SCRIPT = Path(__file__).resolve().parents[1] / "scripts" / "ceus_align_review.py" +_KEY_HEADER = ["id", "video_id", "is_control", "dx", "dy", "video_sha256"] +_VERDICT_HEADER = ["编号", "判定", "依据"] +_REVIEW_HEADER = ["video_sha256", "dx", "dy", "passed", "verdict_a", "verdict_b"] + + +def _load_script() -> ModuleType: + spec = importlib.util.spec_from_file_location("ceus_align_review", _SCRIPT) + assert spec is not None and spec.loader is not None + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def _read(path: Path) -> list[list[str]]: + with path.open(encoding="utf-8", newline="") as handle: + return list(csv.reader(handle)) + + +def _write(path: Path, rows: list[list[str]]) -> None: + with path.open("w", encoding="utf-8", newline="") as handle: + csv.writer(handle).writerows(rows) + + +def test_render_writes_shuffled_images_controls_and_key( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch, capsys: pytest.CaptureFixture[str] +) -> None: + delivery_dir = make_delivery(2).write(tmp_path / "delivery") + videos = { + "U001": delivery_dir / "CEUS_agent" / "001" / "U001_roi2.mp4", + "U002": delivery_dir / "CEUS_agent" / "002" / "U002_roi2.mp4", + } + write_video(videos["U001"], green_box=True) + write_video(videos["U002"], green_box=False) # 第 2 步才会失败，第 1 步照常出图 + document: dict[str, Any] = yaml.safe_load(CEUS_YAML.read_text(encoding="utf-8")) + document["paths"].update( + delivery_dir=str(delivery_dir), + align_runs_root=str(tmp_path / "align"), + align_review=str(tmp_path / "unused_review.csv"), # 出图用不到，文件不存在也行 + ) + document["align_check"]["control_count"] = 1 + config = tmp_path / "ceus.yaml" + config.write_text(yaml.safe_dump(document, allow_unicode=True), encoding="utf-8") + module = _load_script() + + def fake_run_id(now: object, git: object) -> str: + return "run-a" + + monkeypatch.setattr(module, "make_build_id", fake_run_id) + module.main(["render", "--config", str(config)]) + + run_dir = tmp_path / "align" / "run-a" + assert not (tmp_path / "align" / "run-a.partial").exists() + key = _read(run_dir / "key.csv") + assert key[0] == _KEY_HEADER + rows = key[1:] + # 2 例真实 + 1 例对照，编号是 001–003 的一个排列，每个编号两张图。 + assert sorted(row[0] for row in rows) == ["001", "002", "003"] + names = sorted(p.name for p in (run_dir / "images").iterdir()) + assert names == sorted(f"{row[0]}_{kind}.png" for row in rows for kind in ("checker", "edges")) + assert all(re.fullmatch(r"\d{3}_(checker|edges)\.png", name) for name in names) + real = [row for row in rows if row[2] == "false"] + controls = [row for row in rows if row[2] == "true"] + assert sorted(row[1] for row in real) == ["U001", "U002"] + assert [(row[3], row[4]) for row in real] == [("0", "0"), ("0", "0")] + (control,) = controls + # 对照的平移 = 真实配准 (0, 0) + 长度 15–30 像素的随机错开（取整后容差 1 像素）。 + assert 14 <= math.hypot(int(control[3]), int(control[4])) <= 31 + for row in rows: + assert row[5] == hashlib.sha256(videos[row[1]].read_bytes()).hexdigest() + checker = cv2.imread(str(run_dir / "images" / f"{control[0]}_checker.png")) + edges = cv2.imread(str(run_dir / "images" / f"{control[0]}_edges.png")) + assert checker is not None and edges is not None and checker.shape == edges.shape + assert _read(run_dir / "failed.csv") == [["video_id", "step", "reason"]] + assert (run_dir / "snapshot" / "ceus.yaml").read_bytes() == config.read_bytes() + out = capsys.readouterr().out + assert "U001" not in out and "U002" not in out + + +# 手写的 key.csv：两例真实（001、003）、一例对照（002）。 +_KEY = [ + _KEY_HEADER, + ["001", "V1", "false", "-24", "0", "a" * 64], + ["002", "V1", "true", "-4", "17", "a" * 64], + ["003", "V2", "false", "183", "-2", "b" * 64], +] + + +def _score(module: ModuleType, run_dir: Path, a: list[list[str]], b: list[list[str]]) -> None: + _write(run_dir / "key.csv", _KEY) + _write(run_dir / "verdicts_a.csv", [_VERDICT_HEADER, *a]) + _write(run_dir / "verdicts_b.csv", [_VERDICT_HEADER, *b]) + module.main(["score", "--run-dir", str(run_dir)]) + + +def test_score_writes_review_when_controls_caught( + tmp_path: Path, capsys: pytest.CaptureFixture[str] +) -> None: + """对照两人都判「没对齐」；001 两人都判对齐 → 通过，003 有一人看不清 → 不通过。""" + module = _load_script() + a = [["001", "对齐", "血管连续"], ["002", "没对齐", "边缘错开"], ["003", "对齐", "扇形重合"]] + b = [["001", "对齐", "血管连续"], ["002", "没对齐", "边缘错开"], ["003", "看不清", "太暗"]] + _score(module, tmp_path, a, b) + + assert _read(tmp_path / "review.csv") == [ + _REVIEW_HEADER, + ["a" * 64, "-24", "0", "true", "对齐", "对齐"], + ["b" * 64, "183", "-2", "false", "对齐", "看不清"], + ] + out = capsys.readouterr().out + assert "通过 1，不通过 1" in out + assert "「看不清」1 例" in out and "不一致 1 例" in out + # ceus_quant 读得回这份表。 + review = _load_quant().read_review(tmp_path / "review.csv") + assert review["a" * 64].offset == (-24, 0) and review["a" * 64].passed is True + assert review["b" * 64].passed is False + + +def _load_quant() -> ModuleType: + spec = importlib.util.spec_from_file_location("ceus_quant", _SCRIPT.parent / "ceus_quant.py") + assert spec is not None and spec.loader is not None + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +@pytest.mark.parametrize("control_verdict", ["对齐", "看不清"]) +def test_score_is_void_when_a_control_is_missed( + tmp_path: Path, capsys: pytest.CaptureFixture[str], control_verdict: str +) -> None: + module = _load_script() + a = [["001", "对齐", "x"], ["002", "没对齐", "x"], ["003", "对齐", "x"]] + b = [["001", "对齐", "x"], ["002", control_verdict, "x"], ["003", "对齐", "x"]] + with pytest.raises(SystemExit) as raised: + _score(module, tmp_path, a, b) + assert raised.value.code == 1 + assert not (tmp_path / "review.csv").exists() + assert "作废" in capsys.readouterr().out + + +@pytest.mark.parametrize( + ("rows", "message"), + [ + ([["001", "对齐", "x"], ["002", "没对齐", "x"]], "缺"), # 缺 003 + ( + [ + ["001", "对齐", "x"], + ["002", "没对齐", "x"], + ["003", "对齐", "x"], + ["004", "对齐", "x"], + ], + "多", + ), + ([["001", "对齐", "x"], ["002", "没对齐", "x"], ["003", "差不多", "x"]], "判定只能是"), + ( + [ + ["001", "对齐", "x"], + ["001", "对齐", "x"], + ["002", "没对齐", "x"], + ["003", "对齐", "x"], + ], + "重复", + ), + ], +) +def test_score_rejects_malformed_verdicts( + tmp_path: Path, rows: list[list[str]], message: str +) -> None: + module = _load_script() + good = [["001", "对齐", "x"], ["002", "没对齐", "x"], ["003", "对齐", "x"]] + with pytest.raises(ValueError, match=message): + _score(module, tmp_path, rows, good) + assert not (tmp_path / "review.csv").exists() diff --git a/tests/test_ceus_contrast_side.py b/tests/test_ceus_contrast_side.py index d88f143..f568b65 100644 --- a/tests/test_ceus_contrast_side.py +++ b/tests/test_ceus_contrast_side.py @@ -65,8 +65,8 @@ def test_green_box_in_right_half(ceus_config: CeusConfig) -> None: # 整帧 x=250 在右半局部是 250 − 200 = 50。 prompt = side.prompt_frames[0] assert abs(prompt.bbox_xyxy[0] - 50) <= 3 - # 挪到左半（灰阶那一半）：整帧 x = 250 − 200 = 50，局部也是 50。 - mapped = map_source_mask_to_target(prompt.mask, side.contrast_view, side.target_view) + # 挪到左半（灰阶那一半），两半扇形对齐、配准平移为 (0, 0)：局部 x 仍是 50。 + mapped = map_source_mask_to_target(prompt.mask, (0, 0), side.target_view) xs = np.where(mapped)[1] assert abs(int(xs.min()) - 50) <= 3 assert abs(int(xs.max()) - 169) <= 3 diff --git a/tests/test_ceus_mapping.py b/tests/test_ceus_mapping.py index 375d1c6..2972b23 100644 --- a/tests/test_ceus_mapping.py +++ b/tests/test_ceus_mapping.py @@ -1,7 +1,8 @@ -"""两半画面之间的坐标换算：平移量 = 目标视图左上角 − 造影视图左上角（任务书 02 第一阶段第 1 条）。 +"""两半画面之间的坐标换算：造影那一半局部坐标 + 配准平移 (dx, dy) = 另一半局部坐标。 -后两条照搬旧 test_ceus_mapping（去掉锚点配置参数，断言不变）；其余是新写的手算测试。 -旧文件里依赖锚点、有效区或缩放的用例不再成立，已删（见任务书 02 执行报告）。 +前两条照搬旧 test_ceus_mapping（去掉锚点配置参数，平移取 (0, 0)，断言不变）；其余是新写的手算 +测试，平移量直接给定（怎么估平移见 test_ceus_registration）。旧文件里依赖锚点、有效区或缩放的 +用例不再成立，已删（见任务书 02 执行报告）。 """ import cv2 @@ -9,15 +10,11 @@ import numpy as np import pytest from cadx.ceus import ( - CeusConfig, ViewRect, map_source_mask_to_target, polygon_to_mask, - split_views, target_polygons_to_source, - view_offset, ) -from conftest import load_ceus_section def _largest_external_contour(mask: np.ndarray) -> np.ndarray: @@ -49,14 +46,13 @@ def _bbox(points: np.ndarray) -> tuple[float, float, float, float]: def test_map_source_mask_to_target_translates_into_target_effective() -> None: """source 视图局部 mask 经映射落入 target 有效区局部坐标，形状为 target_eff 尺寸。""" - # source view = 整个左半区 [0,100)×[0,100)；有效区同 view（无黑边）。 - source_view = ViewRect(0, 0, 100, 100) + # source view = 整个左半区 [0,100)×[0,100)；有效区同 view（无黑边），平移 (0, 0)。 target_eff = ViewRect(100, 0, 200, 100) source_mask = np.zeros((100, 100), dtype=bool) source_mask[30:50, 20:40] = True # source 局部一个方块 - target_mask = map_source_mask_to_target(source_mask, source_view, target_eff) + target_mask = map_source_mask_to_target(source_mask, (0, 0), target_eff) assert target_mask.dtype == bool assert target_mask.shape == (target_eff.height, target_eff.width) @@ -83,97 +79,94 @@ def test_real_roundtrip_same_size_no_offset_restores_polygon() -> None: # 正向：source 局部多边形 → 栅格化 mask → 映射到 target 有效区局部 mask。 source_mask = polygon_to_mask(source_polygon, source_view.height, source_view.width) - target_mask = map_source_mask_to_target(source_mask, source_view, target_eff) + target_mask = map_source_mask_to_target(source_mask, (0, 0), target_eff) assert target_mask.sum() > 0 # 从真实 target mask 提最大外轮廓得 target 多边形，再逆向映回 source 视图局部。 target_polygon = _largest_external_contour(target_mask) - [back_to_source] = target_polygons_to_source([target_polygon], source_view, target_eff) + [back_to_source] = target_polygons_to_source([target_polygon], (0, 0), source_view) # 形状/包围盒近似还原（栅格化+轮廓量化误差 ≤2px）。 np.testing.assert_allclose(_bbox(back_to_source), _bbox(source_polygon), atol=2.0) -def _seam_frames(width: int, height: int, split_x: int) -> list[np.ndarray]: - """合成几帧左右分屏画面：两半是不同的常量灰度，`split_x` 处一列纯黑缝。""" - frame = np.zeros((height, width, 3), dtype=np.uint8) - frame[:, :split_x] = 180 - frame[:, split_x + 1 :] = 120 - return [frame.copy() for _ in range(5)] - - -# 三种宽度（40 例里有 5 种），分屏线在正中：右半左上角 x = W/2，平移量手算为 (W/2, 0)。 @pytest.mark.parametrize( - ("width", "expected_dx"), - [(1856, 928), (1280, 640), (1024, 512)], -) -def test_view_offset_is_difference_of_view_top_left_corners(width: int, expected_dx: int) -> None: - config = CeusConfig.model_validate(load_ceus_section()) - _, left, right = split_views(_seam_frames(width, 200, expected_dx), config.views) - assert view_offset(left, right) == (expected_dx, 0) - assert view_offset(right, left) == (-expected_dx, 0) - - -@pytest.mark.parametrize( - ("width", "split_x", "box_x1", "expected_target_x1"), + ("offset", "expected_x", "expected_y"), [ - # 造影在左半：整帧 x=300 的框挪到右半，整帧 x = 300 + 928 = 1228。 - (1856, 928, 300, 1228), - # 造影在左半：整帧 x=300 → 300 + 640 = 940。 - (1280, 640, 300, 940), - # 造影在左半：整帧 x=100 → 100 + 512 = 612。 - (1024, 512, 100, 612), + # 局部方块 x 300..339、y 50..79。平移 (−37, 4)：x 263..302、y 54..83。 + ((-37, 4), (263, 302), (54, 83)), + # 平移 (25, −6)：x 325..364、y 44..73。 + ((25, -6), (325, 364), (44, 73)), + # 平移 (0, 0)：原位不动。 + ((0, 0), (300, 339), (50, 79)), ], ) -def test_mask_moves_by_view_offset_in_frame_coordinates( - width: int, split_x: int, box_x1: int, expected_target_x1: int +def test_mask_moves_by_offset( + offset: tuple[int, int], expected_x: tuple[int, int], expected_y: tuple[int, int] ) -> None: - left = ViewRect(0, 0, split_x, 200) - right = ViewRect(split_x, 0, width, 200) - source_mask = np.zeros((left.height, left.width), dtype=bool) - source_mask[50:80, box_x1 - left.x1 : box_x1 - left.x1 + 40] = True + target = ViewRect(640, 0, 1280, 200) + source_mask = np.zeros((200, 640), dtype=bool) + source_mask[50:80, 300:340] = True - target_mask = map_source_mask_to_target(source_mask, left, right) + target_mask = map_source_mask_to_target(source_mask, offset, target) + + assert target_mask.shape == (200, 640) + ys, xs = np.where(target_mask) + assert (int(xs.min()), int(xs.max())) == expected_x + assert (int(ys.min()), int(ys.max())) == expected_y + + +def test_mask_crops_what_leaves_target_view() -> None: + """平移 (−310, 0)：局部 x 300..339 落到 −10..29，出界的 10 列丢掉，不缩放。""" + target = ViewRect(0, 0, 640, 200) + source_mask = np.zeros((200, 640), dtype=bool) + source_mask[10:20, 300:340] = True + + target_mask = map_source_mask_to_target(source_mask, (-310, 0), target) ys, xs = np.where(target_mask) - assert int(xs.min()) + right.x1 == expected_target_x1 - assert int(xs.max()) + right.x1 == expected_target_x1 + 39 - assert (int(ys.min()), int(ys.max())) == (50, 79) - - -def test_mask_moves_from_right_half_to_left_half() -> None: - """造影在右半：1856 宽、分屏线 928，整帧 x=1000 的框挪到左半，整帧 x = 1000 − 928 = 72。""" - left = ViewRect(0, 0, 928, 200) - right = ViewRect(928, 0, 1856, 200) - source_mask = np.zeros((right.height, right.width), dtype=bool) - source_mask[10:20, 1000 - 928 : 1040 - 928] = True - - target_mask = map_source_mask_to_target(source_mask, right, left) - - ys, xs = np.where(target_mask) - assert (int(xs.min()), int(xs.max())) == (72, 111) + assert (int(xs.min()), int(xs.max())) == (0, 29) assert (int(ys.min()), int(ys.max())) == (10, 19) -def test_polygons_move_back_by_negative_view_offset() -> None: - """右半局部 (40, 30) 的点挪回左半：整帧 (968, 30) − 928 = (40, 30)；超出左半的点裁到边界。""" - left = ViewRect(0, 0, 928, 200) - right = ViewRect(928, 0, 1856, 200) - polygon = np.array([[40.0, 30.0], [60.0, 30.0], [927.0, 250.0]], dtype=np.float32) +def test_polygons_move_back_by_negative_offset() -> None: + """平移 (25, 4)：目标局部 (40, 30) 挪回 (15, 26)；超出造影视图的点裁到边界。""" + source = ViewRect(0, 0, 928, 200) + polygon = np.array([[40.0, 30.0], [60.0, 30.0], [10.0, 250.0]], dtype=np.float32) - [back] = target_polygons_to_source([polygon], left, right) + [back] = target_polygons_to_source([polygon], (25, 4), source) - np.testing.assert_allclose(back, [[40.0, 30.0], [60.0, 30.0], [927.0, 199.0]]) + np.testing.assert_allclose(back, [[15.0, 26.0], [35.0, 26.0], [0.0, 199.0]]) def test_mask_crops_when_target_half_is_narrower() -> None: - """宽度为奇数时右半比左半窄 1 像素：不缩放，超出右半的那一列裁掉。""" - left = ViewRect(0, 0, 101, 50) + """宽度为奇数时右半比左半窄 1 像素：平移 (0, 0)，不缩放，超出右半的那一列裁掉。""" right = ViewRect(101, 0, 201, 50) source_mask = np.zeros((50, 101), dtype=bool) source_mask[:, 99:101] = True - target_mask = map_source_mask_to_target(source_mask, left, right) + target_mask = map_source_mask_to_target(source_mask, (0, 0), right) assert target_mask.shape == (50, 100) assert np.where(target_mask)[1].tolist() == [99] * 50 + + +def test_real_roundtrip_with_offset_restores_polygon() -> None: + """平移 (−37, 4) 下真往返：多边形 → 栅格化 → 正向 → 提轮廓 → 逆向，近似还原（误差 ≤2px）。""" + source_view = ViewRect(0, 0, 200, 150) + target_view = ViewRect(200, 0, 400, 150) + source_polygon = np.array( + [[60.0, 30.0], [120.0, 25.0], [130.0, 90.0], [70.0, 85.0]], dtype=np.float32 + ) + source_mask = polygon_to_mask(source_polygon, source_view.height, source_view.width) + target_mask = map_source_mask_to_target(source_mask, (-37, 4), target_view) + target_polygon = _largest_external_contour(target_mask) + np.testing.assert_allclose( + _bbox(target_polygon), + np.array(_bbox(source_polygon)) + np.array([-37, -37, 4, 4]), + atol=2.0, + ) + + [back] = target_polygons_to_source([target_polygon], (-37, 4), source_view) + + np.testing.assert_allclose(_bbox(back), _bbox(source_polygon), atol=2.0) diff --git a/tests/test_ceus_pipeline.py b/tests/test_ceus_pipeline.py index 5c134bd..c21353e 100644 --- a/tests/test_ceus_pipeline.py +++ b/tests/test_ceus_pipeline.py @@ -3,7 +3,8 @@ 视频 400×200、10 fps、40 秒（400 帧），分屏线在 x=200（这一列全黑）。造影画面在右半：病灶是整帧 [250, 369] × [60, 159] 的方块，医生的绿框是它的空心描边；病灶内部亮度按 Lognormal 团注曲线从 基线 40 升到约 170 再回落，内部另有一块一直发黑的囊（暗区）；每帧加一点随机噪声，免得相邻帧 -被当成重复帧。左半是不变的灰阶画面。 +被当成重复帧。两半各有一块位置相同、带纹理、随时间轻微起伏的扇形区（`fake_split_scene`， +配准要靠它），左半是灰阶画面，右半（造影）调暗。 手算的预期： - 绿框在右半 → 造影画面是右半，分割那一半是左半 [0, 200)； @@ -14,19 +15,15 @@ import json import math -from pathlib import Path -import cv2 import numpy as np import pytest from cadx.ceus import ( ANCHOR_SLOTS, CeusConfig, - LognormalParameters, NoValidContourError, SegmentationInput, - evaluate_lognormal, keyframe_images, pack_tic_curves, pack_trajectory, @@ -39,36 +36,11 @@ from cadx.ceus import ( unpack_tic_curves, unpack_trajectory, ) +from fake_ceus_video import BOLUS, FPS, FRAMES, HEIGHT, LESION, SPLIT, WIDTH, write_video -_WIDTH, _HEIGHT, _SPLIT = 400, 200, 200 -_FPS, _FRAMES = 10.0, 400 -_LESION = (250, 60, 369, 159) # 整帧坐标 x1, y1, x2, y2（含端点） -_CYST = (290, 95, 320, 125) -_BOLUS = LognormalParameters(area=1300.0, mu=2.4, sigma=0.45, arrival_offset_s=4.0, baseline=40.0) _INJECTION_OFFSET_S = 5.0 -def _write_video(path: Path) -> None: - rng = np.random.default_rng(3) - times = np.arange(_FRAMES) / _FPS - intensity = evaluate_lognormal(times, _BOLUS) - writer = cv2.VideoWriter(str(path), cv2.VideoWriter.fourcc(*"mp4v"), _FPS, (_WIDTH, _HEIGHT)) - x1, y1, x2, y2 = _LESION - for index in range(_FRAMES): - frame = np.zeros((_HEIGHT, _WIDTH, 3), dtype=np.uint8) - frame[:, :_SPLIT] = 120 - frame[:, _SPLIT + 1 :] = 30 - lesion = np.clip( - intensity[index] + rng.normal(0.0, 6.0, (y2 - y1 + 1, x2 - x1 + 1)), 0, 255 - ) - frame[y1 : y2 + 1, x1 : x2 + 1] = lesion.astype(np.uint8)[:, :, None] - cx1, cy1, cx2, cy2 = _CYST - frame[cy1 : cy2 + 1, cx1 : cx2 + 1] = 10 - cv2.rectangle(frame, (x1, y1), (x2, y2), (0, 255, 0), thickness=3) - writer.write(frame) - writer.release() - - @pytest.fixture(scope="module") def config() -> CeusConfig: from conftest import load_ceus_section @@ -79,7 +51,7 @@ def config() -> CeusConfig: @pytest.fixture(scope="module") def prepared(tmp_path_factory: pytest.TempPathFactory, config: CeusConfig) -> SegmentationInput: path = tmp_path_factory.mktemp("ceus") / "synthetic.mp4" - _write_video(path) + write_video(path, green_box=True) return prepare_segmentation(path, config) @@ -87,21 +59,25 @@ def _lesion_masks(prepared: SegmentationInput) -> np.ndarray: """合成「分割结果」：每帧都是病灶方块挪到左半局部坐标后的位置（略往里收 3 像素）。""" view = prepared.contrast.target_view masks = np.zeros((prepared.frame_count, view.height, view.width), dtype=bool) - x1, y1, x2, y2 = _LESION - masks[:, y1 + 3 : y2 - 2, x1 - _SPLIT + 3 : x2 - _SPLIT - 2] = True + x1, y1, x2, y2 = LESION + masks[:, y1 + 3 : y2 - 2, x1 - SPLIT + 3 : x2 - SPLIT - 2] = True return masks def test_prepare_segmentation_finds_contrast_half_and_maps_prompts( prepared: SegmentationInput, ) -> None: - assert prepared.frame_count == _FRAMES - assert prepared.frame_size == (_WIDTH, _HEIGHT) + assert prepared.frame_count == FRAMES + assert prepared.frame_size == (WIDTH, HEIGHT) assert prepared.contrast.side == "right" + # 两半扇形对齐：两种配准都估出 (0, 0)；造影在右半，平移取反后仍是 (0, 0)。 + assert prepared.registration.correlation == (0, 0) + assert prepared.registration.disagreement_px <= 1 + assert prepared.offset == (0, 0) target = prepared.contrast.target_view assert (target.x1, target.x2) == (0, prepared.split_x) - assert abs(prepared.split_x - _SPLIT) <= 1 - assert len(prepared.target_frames) == _FRAMES + assert abs(prepared.split_x - SPLIT) <= 1 + assert len(prepared.target_frames) == FRAMES assert all(frame.shape == (target.height, target.width, 3) for frame in prepared.target_frames) assert sorted(prepared.prompt_masks) == [p.frame_index for p in prepared.contrast.prompt_frames] assert 0 in prepared.prompt_masks @@ -114,7 +90,7 @@ def test_prepare_segmentation_finds_contrast_half_and_maps_prompts( def test_quantify_case_end_to_end(prepared: SegmentationInput, config: CeusConfig) -> None: result = quantify_case(prepared, _lesion_masks(prepared), _INJECTION_OFFSET_S, config) - assert result.trajectory.source_polygons.shape == (_FRAMES, 64, 2) + assert result.trajectory.source_polygons.shape == (FRAMES, 64, 2) assert result.lesion_fit_ready is True assert list(result.anchors) == list(ANCHOR_SLOTS) for slot in ("D1", "D2", "D3"): @@ -124,13 +100,13 @@ def test_quantify_case_end_to_end(prepared: SegmentationInput, config: CeusConfi assert result.anchors["F2"].status == "not_covered" assert result.anchors["F2"].reason is not None assert result.perfusion["roi_perf"] is not None - peak_time = _BOLUS.arrival_offset_s + math.exp(_BOLUS.mu - _BOLUS.sigma**2) + peak_time = BOLUS.arrival_offset_s + math.exp(BOLUS.mu - BOLUS.sigma**2) assert result.perfusion["roi_perf"].ttp_s == pytest.approx(peak_time, abs=0.5) # 画面上没有紫圈。 assert result.pan_detected_fraction == 0.0 assert result.tic.reference_valid is False # 暗区（囊）被找出来，去暗区曲线比病灶整体亮。 - peak_row = result.tic.rows[round(peak_time * _FPS)] + peak_row = result.tic.rows[round(peak_time * FPS)] assert peak_row["dark_area_px"] > 0 assert peak_row["roi_perf_raw"] > peak_row["roi_all_raw"] @@ -150,7 +126,7 @@ def test_quantify_case_end_to_end(prepared: SegmentationInput, config: CeusConfi assert len(render_case_overlay(prepared, result)) > 0 images = keyframe_images(prepared, result.anchors) assert "F1" in images and "F2" not in images - assert all(image.shape == (_HEIGHT, _WIDTH, 3) for image in images.values()) + assert all(image.shape == (HEIGHT, WIDTH, 3) for image in images.values()) def test_empty_segmentation_fails_case_at_step_4( diff --git a/tests/test_ceus_polygons.py b/tests/test_ceus_polygons.py index 0d88f59..d08225b 100644 --- a/tests/test_ceus_polygons.py +++ b/tests/test_ceus_polygons.py @@ -1,5 +1,5 @@ -"""polygons 的单元测试（照搬旧项目 test_ceus_polygons.py；坐标换算改为两个视图矩形， -去掉旧的有效区矩形和映射锚点参数）。 +"""polygons 的单元测试（照搬旧项目 test_ceus_polygons.py；坐标换算改为配准得到的平移 `offset`， +去掉旧的有效区矩形和映射锚点参数。旧用例里两个视图左上角相同，平移取 (0, 0)，断言不变）。 旧标题：B6 polygons.py 的单元测试：mask→平滑多边形 + 位置约束 + 绿框先验融合。 @@ -104,7 +104,7 @@ def test_masks_to_smoothed_polygons_all_empty_raises() -> None: def test_position_constraint_pulls_center_to_keyframe_track_without_area_change() -> None: """strength=1.0 把每帧质心拉到 keyframe 中心插值轨迹，且整体平移不改面积。""" target_eff = ViewRect(0, 0, 200, 200) - source_view = ViewRect(0, 0, 200, 200) + offset = (0, 0) # 两个 keyframe：source mask 中心分别在 (40, 100) 与 (160, 100) keyframes = [ @@ -118,7 +118,7 @@ def test_position_constraint_pulls_center_to_keyframe_track_without_area_change( corrected = apply_prompt_frame_position_constraints( polygons, keyframes, - source_view, + offset, target_eff, _CONFIG.constraints.keyframe_position, # enabled、center_interpolation、strength=1.0 ) @@ -138,14 +138,14 @@ def test_position_constraint_pulls_center_to_keyframe_track_without_area_change( def test_position_constraint_disabled_returns_input() -> None: """enabled=False 时原样返回输入（门关）。""" target_eff = ViewRect(0, 0, 100, 100) - source_view = ViewRect(0, 0, 100, 100) + offset = (0, 0) polygons = [_circle_polygon((50.0, 50.0), 10.0, 64)] keyframes = [{"frame_index": 0, "mask": _disk_mask(100, 100, (50, 50), 10)}] out = apply_prompt_frame_position_constraints( polygons, keyframes, - source_view, + offset, target_eff, _CONFIG.constraints.keyframe_position.model_copy(update={"enabled": False}), ) @@ -156,14 +156,14 @@ def test_position_constraint_disabled_returns_input() -> None: def test_fusion_gate_closed_returns_pure_medsam() -> None: """roi_shape.source != 'prompt_prior' → 直接返回纯 medsam（门关，findings §3.7）。""" target_eff = ViewRect(0, 0, 100, 100) - source_view = ViewRect(0, 0, 100, 100) + offset = (0, 0) medsam = [_circle_polygon((50.0, 50.0), 10.0, 64)] keyframes = [{"frame_index": 0, "mask": _disk_mask(100, 100, (50, 50), 10)}] out = fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, _CONFIG.roi_shape.model_copy(update={"source": "medsam2"}), _CONFIG.smoothing, @@ -175,7 +175,7 @@ def test_fusion_gate_closed_returns_pure_medsam() -> None: def test_fusion_keyframe_strength_one_yields_prior_shape() -> None: """keyframe 帧 strength=1.0：融合结果几何上贴近绿框先验（逐顶点混合 s=1）。""" target_eff = ViewRect(0, 0, 200, 200) - source_view = ViewRect(0, 0, 200, 200) + offset = (0, 0) # 绿框先验中心 (100, 100) 半径 30；medsam 偏移到 (130, 130) 半径 30 keyframes = [{"frame_index": 0, "mask": _disk_mask(200, 200, (100, 100), 30)}] @@ -184,7 +184,7 @@ def test_fusion_keyframe_strength_one_yields_prior_shape() -> None: fused = fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, _CONFIG.roi_shape.model_copy(update={"keyframe_shape_strength": 1.0}), _CONFIG.smoothing, @@ -200,7 +200,7 @@ def test_fusion_keyframe_strength_one_yields_prior_shape() -> None: def test_fusion_area_clamped_to_prior_band() -> None: """融合面积超过 max_area_ratio_to_prompt 时按 sqrt 绕先验质心缩放回界（findings §3.7）。""" target_eff = ViewRect(0, 0, 400, 400) - source_view = ViewRect(0, 0, 400, 400) + offset = (0, 0) # 关键帧只在 frame 0；frame 1 为传播帧（非 keyframe），用 propagation_shape_strength=0 # 让 medsam 形状主导，从而触发面积收紧。先验半径 30，medsam 半径 80（面积远超 1.25×）。 @@ -222,7 +222,7 @@ def test_fusion_area_clamped_to_prior_band() -> None: fused = fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, shape_cfg, _CONFIG.smoothing, @@ -246,7 +246,7 @@ def test_fusion_propagation_frame_linear_blend_without_tightening() -> None: 使逐顶点线性混合的几何结果可被精确预测（绕开起点对齐 / 收紧的扰动）。 """ target_eff = ViewRect(0, 0, 400, 400) - source_view = ViewRect(0, 0, 400, 400) + offset = (0, 0) # 单关键帧（frame 0），frame 1 为传播帧；先验半径 50，medsam 半径 60（面积比≈1.44， # 但 max_ratio 放宽到 2.0 → 不触发收紧），中心同在 (200, 200) @@ -267,7 +267,7 @@ def test_fusion_propagation_frame_linear_blend_without_tightening() -> None: fused = fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, shape_cfg, _CONFIG.smoothing, @@ -297,7 +297,7 @@ def test_fusion_propagation_frame_linear_blend_without_tightening() -> None: def test_fusion_lower_bound_tightening_expands_area() -> None: """下界收紧：medsam 面积 < min_ratio × prior 面积 → 被放大到 ≈ min_ratio × prior。""" target_eff = ViewRect(0, 0, 400, 400) - source_view = ViewRect(0, 0, 400, 400) + offset = (0, 0) # 先验半径 60；medsam 传播帧半径 20（面积比≈0.11 < min_ratio 0.70 → 触发下界放大） keyframes = [{"frame_index": 0, "mask": _disk_mask(400, 400, (200, 200), 60)}] @@ -316,7 +316,7 @@ def test_fusion_lower_bound_tightening_expands_area() -> None: fused = fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, shape_cfg, _CONFIG.smoothing, @@ -336,7 +336,7 @@ def test_fusion_start_point_misalignment_does_not_distort_shape() -> None: 若不先对齐起点，逐顶点混合会把不同角度的点平均出来导致半径塌缩；align 后应保形。 """ target_eff = ViewRect(0, 0, 400, 400) - source_view = ViewRect(0, 0, 400, 400) + offset = (0, 0) radius = 50.0 keyframes = [{"frame_index": 0, "mask": _disk_mask(400, 400, (200, 200), int(radius))}] @@ -357,7 +357,7 @@ def test_fusion_start_point_misalignment_does_not_distort_shape() -> None: fused = fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, shape_cfg, _CONFIG.smoothing, @@ -374,7 +374,7 @@ def test_fusion_start_point_misalignment_does_not_distort_shape() -> None: def test_fusion_propagation_strength_clipped_to_unit_interval() -> None: """传播权重越界（1.5 / −0.5）经函数内部 np.clip → 行为等价 1.0 / 0.0。""" target_eff = ViewRect(0, 0, 400, 400) - source_view = ViewRect(0, 0, 400, 400) + offset = (0, 0) keyframes = [{"frame_index": 0, "mask": _disk_mask(400, 400, (200, 200), 40)}] medsam = [ @@ -392,7 +392,7 @@ def test_fusion_propagation_strength_clipped_to_unit_interval() -> None: return fuse_prompt_prior_with_medsam( medsam, keyframes, - source_view, + offset, target_eff, shape_cfg, _CONFIG.smoothing, @@ -465,7 +465,7 @@ def test_temporal_smoothing_no_vertex_crossing_distortion() -> None: def test_position_constraint_invalid_mode_raises() -> None: """mode != 'center_interpolation' → RuntimeError（当前仅支持该模式，纵深防御）。""" target_eff = ViewRect(0, 0, 200, 200) - source_view = ViewRect(0, 0, 200, 200) + offset = (0, 0) polygons = [_circle_polygon((100.0, 100.0), 20.0, 64)] keyframes = [{"frame_index": 0, "mask": _disk_mask(200, 200, (100, 100), 20)}] @@ -473,7 +473,7 @@ def test_position_constraint_invalid_mode_raises() -> None: apply_prompt_frame_position_constraints( polygons, keyframes, - source_view, + offset, target_eff, _RawConstraintCfg(mode="affine"), ) @@ -482,7 +482,7 @@ def test_position_constraint_invalid_mode_raises() -> None: def test_position_constraint_strength_clipped_to_unit_interval() -> None: """strength 越界（1.5 / −0.5）经函数内部 np.clip → 行为等价 1.0 / 0.0。""" target_eff = ViewRect(0, 0, 200, 200) - source_view = ViewRect(0, 0, 200, 200) + offset = (0, 0) keyframes = [ {"frame_index": 0, "mask": _disk_mask(200, 200, (40, 100), 15)}, {"frame_index": 2, "mask": _disk_mask(200, 200, (160, 100), 15)}, @@ -493,7 +493,7 @@ def test_position_constraint_strength_clipped_to_unit_interval() -> None: corrected = apply_prompt_frame_position_constraints( [poly.copy() for poly in base_polys], keyframes, - source_view, + offset, target_eff, _RawConstraintCfg(strength=strength), ) diff --git a/tests/test_ceus_quant.py b/tests/test_ceus_quant.py index b415976..181bf19 100644 --- a/tests/test_ceus_quant.py +++ b/tests/test_ceus_quant.py @@ -10,6 +10,7 @@ """ import csv +import hashlib import importlib.util import json from collections.abc import Callable, Mapping, Sequence @@ -26,17 +27,15 @@ import yaml from cadx.answers import read_aliases from cadx.build import run_build -from cadx.ceus import ANCHOR_SLOTS, LognormalParameters, evaluate_lognormal +from cadx.ceus import ANCHOR_SLOTS from cadx.library import list_aliases, read_ultrasound from conftest import CEUS_YAML +from fake_ceus_video import HEIGHT, LESION, SPLIT, WIDTH, write_video from fake_delivery import REPO_ROOT, make_delivery, write_config _SCRIPT = Path(__file__).resolve().parents[1] / "scripts" / "ceus_quant.py" -_WIDTH, _HEIGHT, _SPLIT = 400, 200, 200 -_FPS, _FRAMES = 10.0, 400 -_LESION = (250, 60, 369, 159) # 整帧坐标 x1, y1, x2, y2（含端点） -_CYST = (290, 95, 320, 125) -_BOLUS = LognormalParameters(area=1300.0, mu=2.4, sigma=0.45, arrival_offset_s=4.0, baseline=40.0) +# review.csv 的列，按任务书 02「配准的目视核对」第 3 步独立写，不取脚本里的常量。 +_REVIEW_HEADER = ("video_sha256", "dx", "dy", "passed", "verdict_a", "verdict_b") def _load_script() -> ModuleType: @@ -47,33 +46,14 @@ def _load_script() -> ModuleType: return module -def _write_video(path: Path, green_box: bool) -> None: - rng = np.random.default_rng(3) - intensity = evaluate_lognormal(np.arange(_FRAMES) / _FPS, _BOLUS) - writer = cv2.VideoWriter(str(path), cv2.VideoWriter.fourcc(*"mp4v"), _FPS, (_WIDTH, _HEIGHT)) - x1, y1, x2, y2 = _LESION - for index in range(_FRAMES): - frame = np.zeros((_HEIGHT, _WIDTH, 3), dtype=np.uint8) - frame[:, :_SPLIT] = 120 - frame[:, _SPLIT + 1 :] = 30 - lesion = np.clip( - intensity[index] + rng.normal(0.0, 6.0, (y2 - y1 + 1, x2 - x1 + 1)), 0, 255 - ) - frame[y1 : y2 + 1, x1 : x2 + 1] = lesion.astype(np.uint8)[:, :, None] - cx1, cy1, cx2, cy2 = _CYST - frame[cy1 : cy2 + 1, cx1 : cx2 + 1] = 10 - if green_box: - cv2.rectangle(frame, (x1, y1), (x2, y2), (0, 255, 0), thickness=3) - writer.write(frame) - writer.release() - - def _config_dict(tmp_path: Path, delivery_dir: Path) -> dict[str, Any]: document: dict[str, Any] = yaml.safe_load(CEUS_YAML.read_text(encoding="utf-8")) document["paths"] = { "delivery_dir": str(delivery_dir), "runs_root": str(tmp_path / "runs"), "work_root": str(tmp_path / "work"), + "align_runs_root": str(tmp_path / "runs_align"), + "align_review": str(tmp_path / "review.csv"), } checkpoint = tmp_path / "fake_checkpoint.pt" checkpoint.write_bytes(b"fake weights") @@ -89,6 +69,17 @@ def _write_config(tmp_path: Path, document: dict[str, Any], name: str) -> Path: return path +def _write_review(document: dict[str, Any], entries: dict[str, tuple[int, int, bool]]) -> None: + """按视频号写配准目视核对表：(dx, dy, 是否通过)，视频 SHA-256 现算。""" + videos = Path(document["paths"]["delivery_dir"]) / "CEUS_agent" + lines = [",".join(_REVIEW_HEADER)] + for video_id, (dx, dy, passed) in entries.items(): + video = videos / f"{int(video_id[1:]):03d}" / f"{video_id}_roi2.mp4" + sha = hashlib.sha256(video.read_bytes()).hexdigest() + lines.append(f"{sha},{dx},{dy},{str(passed).lower()},对齐,对齐") + Path(document["paths"]["align_review"]).write_text("\n".join(lines) + "\n", encoding="utf-8") + + class _CountingSegment: """假分割：每帧都是病灶方块挪到左半局部坐标后的位置（略往里收 3 像素）。""" @@ -104,8 +95,8 @@ class _CountingSegment: self.calls += 1 height, width = frames[0].shape[:2] masks = np.zeros((len(frames), height, width), dtype=bool) - x1, y1, x2, y2 = _LESION - masks[:, y1 + 3 : y2 - 2, x1 - _SPLIT + 3 : x2 - _SPLIT - 2] = True + x1, y1, x2, y2 = LESION + masks[:, y1 + 3 : y2 - 2, x1 - SPLIT + 3 : x2 - SPLIT - 2] = True return masks @@ -115,8 +106,8 @@ def setup( ) -> tuple[ModuleType, _CountingSegment, dict[str, Any], Callable[[str], None]]: """写假交付包（两例真视频），载入脚本并换掉分割、MedSAM2 提交号和运行号。""" delivery_dir = make_delivery(2).write(tmp_path / "delivery") - _write_video(delivery_dir / "CEUS_agent" / "001" / "U001_roi2.mp4", green_box=True) - _write_video(delivery_dir / "CEUS_agent" / "002" / "U002_roi2.mp4", green_box=False) + write_video(delivery_dir / "CEUS_agent" / "001" / "U001_roi2.mp4", green_box=True) + write_video(delivery_dir / "CEUS_agent" / "002" / "U002_roi2.mp4", green_box=False) module = _load_script() segment = _CountingSegment() monkeypatch.setattr(module, "segment_video", segment) @@ -130,7 +121,10 @@ def setup( monkeypatch.setattr(module, "medsam2_commit", fake_commit) monkeypatch.setattr(module, "make_build_id", fake_run_id) - return module, segment, _config_dict(tmp_path, delivery_dir), run_ids.append + document = _config_dict(tmp_path, delivery_dir) + # 两段合成视频两半扇形对齐，配准估出 (0, 0)（test_ceus_pipeline 断言过），目视核对都通过。 + _write_review(document, {"U001": (0, 0, True), "U002": (0, 0, True)}) + return module, segment, document, run_ids.append def test_run_writes_outputs_and_failed_case( @@ -157,6 +151,8 @@ def test_run_writes_outputs_and_failed_case( assert (case_dir / name).stat().st_size > 0, name summary = json.loads((case_dir / "summary.json").read_text(encoding="utf-8")) assert summary["video_id"] == "U001" + assert summary["registration"]["contrast_side"] == "right" + assert summary["registration"]["contrast_to_target_offset"] == [0, 0] assert summary["injection_offset_s"] == 5.0 assert list(summary["anchor_frames"]) == list(ANCHOR_SLOTS) assert summary["perfusion_metrics"]["roi_perf"] is not None @@ -167,10 +163,11 @@ def test_run_writes_outputs_and_failed_case( keyframes = {p.stem for p in (case_dir / "keyframes").iterdir()} assert keyframes == located image = cv2.imread(str(case_dir / "keyframes" / "F1.png")) - assert image is not None and image.shape == (_HEIGHT, _WIDTH, 3) # 整幅原始画面 + assert image is not None and image.shape == (HEIGHT, WIDTH, 3) # 整幅原始画面 snapshot = run_dir / "snapshot" assert (snapshot / "ceus.yaml").read_bytes() == (tmp_path / "ceus.yaml").read_bytes() + assert (snapshot / "align_review.csv").read_bytes() == (tmp_path / "review.csv").read_bytes() segment_record = json.loads((snapshot / "segment.json").read_text(encoding="utf-8")) assert segment_record["medsam2_commit"] == "c" * 40 assert segment_record["cuda_peak_memory_bytes"] is None # device 是 cpu @@ -220,7 +217,8 @@ def test_unreadable_video_fails_at_step_1_and_next_case_runs( # 第 1 例换成损坏的视频（cv2 打不开），第 2 例换成正常视频。 videos = Path(document["paths"]["delivery_dir"]) / "CEUS_agent" (videos / "001" / "U001_roi2.mp4").write_bytes(b"not a real video") - _write_video(videos / "002" / "U002_roi2.mp4", green_box=True) + write_video(videos / "002" / "U002_roi2.mp4", green_box=True) + _write_review(document, {"U002": (0, 0, True)}) # 第 1 步失败的视频不进核对表 next_run_id("run-a") module.main(["--config", str(_write_config(tmp_path, document, "ceus.yaml"))]) @@ -234,6 +232,56 @@ def test_unreadable_video_fails_at_step_1_and_next_case_runs( assert segment.calls == 1 +def test_review_rejected_case_fails_at_step_1( + tmp_path: Path, + setup: tuple[ModuleType, _CountingSegment, dict[str, Any], Callable[[str], None]], +) -> None: + """目视核对没通过的病例第 1 步失败、不调分割；另一例照常（这里是第 2 步失败）。""" + module, segment, document, next_run_id = setup + _write_review(document, {"U001": (0, 0, False), "U002": (0, 0, True)}) + next_run_id("run-a") + module.main(["--config", str(_write_config(tmp_path, document, "ceus.yaml"))]) + + with (tmp_path / "runs" / "run-a" / "failed.csv").open(encoding="utf-8", newline="") as handle: + rows = list(csv.reader(handle)) + assert rows[1] == ["U001", "1", "配准目视核对未通过"] + assert rows[2][:2] == ["U002", "2"] + assert segment.calls == 0 + + +@pytest.mark.parametrize( + ("entries", "message"), + [ + ({"U002": (0, 0, True)}, "没有这段视频"), # U001 不在表里 + ({"U001": (3, 0, True), "U002": (0, 0, True)}, "不同"), # 表里的平移与这次算出的不同 + ], +) +def test_review_table_mismatch_stops_run( + tmp_path: Path, + setup: tuple[ModuleType, _CountingSegment, dict[str, Any], Callable[[str], None]], + entries: dict[str, tuple[int, int, bool]], + message: str, +) -> None: + module, segment, document, next_run_id = setup + _write_review(document, entries) + next_run_id("run-a") + with pytest.raises(RuntimeError, match=message): + module.main(["--config", str(_write_config(tmp_path, document, "ceus.yaml"))]) + assert not (tmp_path / "runs" / "run-a").exists() + assert segment.calls == 0 + + +def test_missing_review_file_stops_run( + tmp_path: Path, + setup: tuple[ModuleType, _CountingSegment, dict[str, Any], Callable[[str], None]], +) -> None: + module, _, document, next_run_id = setup + Path(document["paths"]["align_review"]).unlink() + next_run_id("run-a") + with pytest.raises(FileNotFoundError): + module.main(["--config", str(_write_config(tmp_path, document, "ceus.yaml"))]) + + def test_segment_cache_reused_and_keyed_on_step_1_3_params( tmp_path: Path, setup: tuple[ModuleType, _CountingSegment, dict[str, Any], Callable[[str], None]], @@ -292,6 +340,9 @@ def _key(module: ModuleType, document: dict[str, Any]) -> str: ("ceus", ("views", "split_search_band"), [0.3, 0.7], True), ("ceus", ("prompt_source", "min_area"), 9000, True), ("ceus", ("keyframes", "iou_threshold"), 0.9, True), + ("ceus", ("registration", "fan_std_threshold"), 5.0, True), # 第 1 步配准 + ("ceus", ("registration", "max_disagreement_px"), 9, True), + ("ceus", ("pan_source", "hsv_lower"), [120, 60, 60], True), # 配准时修补紫圈颜色 ("segment", ("sam2_config",), "configs/other.yaml", True), ("segment", ("checkpoint",), "models/other.pt", True), ("segment", ("device",), "cuda:1", False), diff --git a/tests/test_ceus_registration.py b/tests/test_ceus_registration.py new file mode 100644 index 0000000..2e5f02b --- /dev/null +++ b/tests/test_ceus_registration.py @@ -0,0 +1,144 @@ +"""两半配准：每段视频单独估左半 → 右半的平移（任务书 02 第一阶段第 1 条）。 + +合成分屏画面（`fake_split_scene`）：每半 240×200，左半扇形是一个梯形，右半的扇形和纹理是左半平移 +(dx, dy) 后的样子，再换一种亮度、加独立噪声。预期值就是造图时设定的平移，不从实现里推。 +""" + +import numpy as np +import pytest + +from cadx.ceus import ( + CeusCaseError, + CeusConfig, + RegistrationError, + ViewRect, + register_halves, +) +from fake_split_scene import SplitScene, make_scene, scene_gray, to_bgr + +_HEIGHT, _HALF = 200, 240 +_LEFT, _RIGHT = ViewRect(0, 0, _HALF, _HEIGHT), ViewRect(_HALF, 0, 2 * _HALF, _HEIGHT) +# 左半局部坐标的梯形扇形：上窄下宽，顶边在顶部信息栏（200 × 0.0834 ≈ 17 行）之下。 +_FAN = [(80, 30), (170, 30), (215, 185), (35, 185)] + + +def _frames(scene: SplitScene, *, green_box: bool) -> list[np.ndarray]: + rng = np.random.default_rng(11) + frames: list[np.ndarray] = [] + for index in range(12): + frame = to_bgr(scene_gray(scene, index, rng)) + if green_box: + # 绿框只画在左半：配准前要修补掉，否则它的边缘只在一半上有。 + frame[90:140, 100] = frame[90:140, 160] = (0, 255, 0) + frame[90, 100:161] = frame[140, 100:161] = (0, 255, 0) + frames.append(frame) + return frames + + +@pytest.mark.parametrize( + ("offset", "green_box"), + [ + ((-37, 0), False), # 负的水平偏移（真实视频里最常见：造影那一半扇形左边有黑边） + ((25, 4), False), # 正的水平偏移、非零竖直偏移 + ((-12, -6), True), # 两个分量都为负，左半画着绿框 + ((0, 0), False), + ], +) +def test_known_offset_is_recovered_by_both_methods( + ceus_config: CeusConfig, offset: tuple[int, int], green_box: bool +) -> None: + scene = make_scene( + height=_HEIGHT, + half_width=_HALF, + fan_polygon=_FAN, + offset=offset, + interior_shift=None, + texture_std=25.0, + right_gain=0.8, + right_bias=10.0, + seed=5, + ) + registration = register_halves( + _frames(scene, green_box=green_box), + _LEFT, + _RIGHT, + ceus_config.registration, + ceus_config.prompt_source, + ceus_config.pan_source, + ) + assert registration.correlation == offset + assert registration.template == offset + assert registration.disagreement_px == 0 + # 造影在左半时平移照用；造影在右半时（把右半的东西挪到左半）取反。 + assert registration.contrast_to_target("left") == offset + assert registration.contrast_to_target("right") == (-offset[0], -offset[1]) + + +# 冲突画面用更宽的扇形：模板块落在扇形中部，纹理平移后真正的匹配位置仍在右半扇形内。 +_WIDE_FAN = [(40, 25), (200, 25), (235, 195), (5, 195)] + + +def _conflicting_scene() -> SplitScene: + """扇形轮廓按 (15, 0) 平移，里面的纹理按 (−15, 0) 平移：两种方法各信一样。 + + 纹理起伏（标准差 5）远小于扇形内外的亮度差（80），边缘图上占主导的是扇形轮廓，互相关跟着 + 轮廓走；模板块整块落在扇形里，只看得到纹理，跟着纹理走。 + """ + return make_scene( + height=_HEIGHT, + half_width=_HALF, + fan_polygon=_WIDE_FAN, + offset=(15, 0), + interior_shift=(-15, 0), + texture_std=5.0, + right_gain=0.8, + right_bias=10.0, + seed=5, + ) + + +def test_methods_disagree_on_conflicting_scene(ceus_config: CeusConfig) -> None: + """门槛放到极大时两种方法各自的估计：互相关跟轮廓 (15, 0)，模板匹配跟纹理 (−15, 0)。""" + loose = ceus_config.registration.model_copy(update={"max_disagreement_px": 10_000}) + registration = register_halves( + _frames(_conflicting_scene(), green_box=False), + _LEFT, + _RIGHT, + loose, + ceus_config.prompt_source, + ceus_config.pan_source, + ) + assert registration.correlation == (15, 0) + assert registration.template == (-15, 0) + assert registration.disagreement_px == 30 + + +def test_disagreement_beyond_threshold_fails_case_at_step_1(ceus_config: CeusConfig) -> None: + with pytest.raises(RegistrationError, match="两半配准不一致") as raised: + register_halves( + _frames(_conflicting_scene(), green_box=False), + _LEFT, + _RIGHT, + ceus_config.registration, + ceus_config.prompt_source, + ceus_config.pan_source, + ) + assert isinstance(raised.value, CeusCaseError) + assert raised.value.step == 1 + # 原因里附两个估计值。 + assert "(15, 0)" in raised.value.reason and "(-15, 0)" in raised.value.reason + + +def test_static_half_has_no_fan_and_fails(ceus_config: CeusConfig) -> None: + """一半画面完全不动（没有扇形区）：第 1 步失败，不猜一个平移。""" + frames = [np.full((_HEIGHT, 2 * _HALF, 3), 60, dtype=np.uint8) for _ in range(12)] + with pytest.raises(RegistrationError, match="找不到扇形区") as raised: + register_halves( + frames, + _LEFT, + _RIGHT, + ceus_config.registration, + ceus_config.prompt_source, + ceus_config.pan_source, + ) + assert raised.value.step == 1