diff --git a/REGISTRATION_BENCHMARK_AFFINE.md b/REGISTRATION_BENCHMARK_AFFINE.md new file mode 100644 index 0000000..f94efb6 --- /dev/null +++ b/REGISTRATION_BENCHMARK_AFFINE.md @@ -0,0 +1,231 @@ +# Registration Pipeline Benchmark — CT → T1c (M6) — Affine + +Companion to [REGISTRATION_BENCHMARK_RIGID.md](REGISTRATION_BENCHMARK_RIGID.md) +and [REGISTRATION_BENCHMARK_DEFORMABLE.md](REGISTRATION_BENCHMARK_DEFORMABLE.md): +the **affine** (12-DOF) variant of the same benchmark. The rigid benchmark +asked "which toolkit aligns the head best with 6 DOF?"; this one asks "how +much does the full 12-DOF affine buy on top of rigid — and do the toolkits' +affine workflows still rank the same way?" — on the same 10 M6 CT→T1c test +cases, the same unified contract, and the same metrics. + +Run date: 2026-10-03, 13:50–14:40 local (single run, 6 workers, 24 threads, +2 workers per GPU; ~50 min wall, dominated by `impact`'s ~6 min phase). +10 M6 test cases, 3× NVIDIA A100 80GB PCIe, 96-core host, conda env +`longitudinal-313`. + +## Unified contract + +Identical to the rigid benchmark: the returned transform `T` must satisfy +`sitk.Resample(ct, t1c, T, ...)` = CT on the T1c grid, i.e. **T maps T1c → +CT coordinates**. fixed = T1c, moving = CT, full volumes (no head box, no +mask, no crop, no pad), each toolkit's native automatic transform +initialization, CT clipped to [0, 80] HU only for the pipelines where the +rigid A/B showed clipping helps (niftyreg, synthmorph). + +One change vs the rigid benchmark: every pipeline estimates a **full affine +(12-DOF)** transform instead of a rigid (6-DOF) one — each toolkit's native +affine workflow (or, for the antspyx rows, a stock linear-only preset). + +## The nine pipelines (affine configurations) + +| pipeline | affine configuration | metric | +|---|---|---| +| **antspyx** | antspyx stock `Affine` preset (the toolkit's native affine): Affine[0.25], 2100×1200×1200×10 iters, shrink 6×4×2×1, smoothing 3×2×1×0, center-of-mass auto init | Mattes MI (32 bins, 20 % regular sampling) | +| **ants_trsaa** | antspyx stock `TRSAA` preset: Translation[1] → Rigid[1] → Similarity[1] (each 2000×2000×0 iters, 4×2×1, 2×1×0) → Affine[1] × 2 (each 40×20×0, 4×2×1, 2×1×0), COM init — the full TRSA cascade | Mattes MI (32 bins), all stages | +| **ants_reproa** | `antsRegistrationSyNRepro[a]`: Rigid[0.1] + Affine[0.1] (each 1000×500×250×100, 8×4×2×1, 3×2×1×0vox), COM init — the linear-only half of the Repro workflow | GC (1/1/Regular/0.25), both stages | +| **ants_quickreproa** | `antsRegistrationSyNQuickRepro[a]`: Rigid[0.1] + Affine[0.1] (each 1000×500×250×0 — no full-res level —, 8×4×2×1, 3×2×1×0vox), COM init | GC (1/1/Regular/0.25), both stages | +| **elastix** | AffineTransform (12-DOF, the rigid benchmark's map with Euler→Affine), ASGD 2000 iters × 3 resolutions, CenterOfGravity auto init | Mattes MI (32 bins), 5000 samples/level | +| **impact** | the ImpactReg `MR_CT_TS` reference map with **AffineTransform** in place of EulerTransform (same MIND + TotalSegmentator static features, 4 scales 400/300/200/200 iters, 2000 samples, same CoG-init / 0.1 sample-gate adaptations) | IMPACT semantic (L1 + Dice) + Mattes MI | +| **niftyreg** | reg_aladin **full affine** (no `-rigOnly`), 3-level pyramid, 15 TLS iters/level | NCC block matching (init) + TLS | +| **simpleitk** | two-stage (the toolkit's canonical affine workflow): Euler3D stage (2000 iters/level, pyramid [4,2,1], identity start) → AffineTransform stage (same budget) initialized to reproduce the rigid stage exactly — a single-stage affine from identity collapses its scale on these pairs | Mattes MI (32 bins), both stages | +| **synthmorph** | `affine` matrix model (the toolkit's native 12-DOF deep model) | anatomy-aware contrast-invariant network (pretrained) | + +The `ants_trsaa` / `ants_reproa` / `ants_quickreproa` rows are the linear-only +(`[a]`) antspyx presets — the same stage configurations the +`antsRegistrationSyNRepro*.sh` / `antsRegistrationSyNQuickRepro*.sh` +workflows use for their Rigid + Affine head. + +## Test cases + +Same 10 as the rigid / deformable benchmarks (first 10 M6 labelled cases of +`m6.discover()`, all with a GTV label): +`m6_22GD5VL2_20230724`, `m6_22IRRKU6_20240115`, `m6_23XIMO4M_20230222`, +`m6_24V62QQM_20230223`, `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820`, +`m6_26TYQI3C_20221027`, `m6_2BBN5XA4_20250411`, `m6_2C4TKI6H_20240506`, +`m6_2CDNN7BH_20240102`. + +(`m6_22GD5VL2_20230724` has an incomplete CT head — a known data issue, +not a pipeline failure; it is scored as ct_incomplete in every pipeline.) + +## Metrics (identical for every pipeline) + +- **Mattes MI** (nats) between the warped CT and the T1c over the MR head + (largest connected component of the T1c mask) — achieved alignment, higher + = better. +- **IMPACT loss** (vboussot/ImpactLoss; pretrained TotalSegmentator + `TS/M730` deep features — weighted L1 over the two deepest layers, full + volumes; the same scorer as the other benchmarks) between the warped CT + and the T1c — semantic-anatomy similarity, lower = better. +- **Registration wall time** (seconds, excludes I/O + scoring). + +Per-case JSONs additionally retain the production skull-ring QA (gap + +concentricity, and the acceptance verdict) for reference. + +### Result + +Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ T1c +over the MR head; higher = better) **and** the **ImpactLoss** semantic score +(warped-CT ~ T1c over the full volumes; lower = better), on the same +unclipped scoring volumes as the other benchmarks, so all rows are +comparable on the same scales and against the rigid / deformable numbers. +`niftyreg` and `synthmorph` clip the CT (the moving image) to [0, 80] HU +*before registration* (per the rigid benchmark's clipping A/B); the others +register the CT as-is. All nine transforms are affine (12-DOF). + +All nine pipelines registered 9/9 cases (the 10th is the known +incomplete-CT data issue, unscored in every pipeline). **Bold** = best for +the column (MI higher is better, ImpactLoss lower is better). + +| pipeline | toolkit | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) | +|---|---|---|---|---|---| +| **antspyx** | ANTsPy 0.6.3 (`Affine` preset) | 9 | 12.6 | 0.156 / 0.155 | **0.660 / 0.660** | +| **ants_trsaa** | ANTsPy 0.6.3 (`TRSAA` preset) | 9 | 30.4 | 0.156 / 0.155 | 0.658 / 0.659 | +| **ants_reproa** | ANTsPy 0.6.3 (`Repro[a]`) | 9 | 63.1 | 0.160 / 0.159 | 0.525 / 0.514 | +| **ants_quickreproa** | ANTsPy 0.6.3 (`QuickRepro[a]`) | 9 | 28.4 | 0.161 / 0.159 | 0.523 / 0.514 | +| **elastix** | ITK-Elastix 0.25.4 | 9 | 12.3 | 0.156 / 0.156 | 0.651 / 0.645 | +| **impact** | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT, affine MIND+TS reference map) | 9 | 215.5 | 0.159 / 0.158 | 0.623 / 0.567 | +| **niftyreg** | NiftyReg 2.1.1 (reg_aladin full affine, Docker + CUDA) | 9 | 19.8 | 0.167 / 0.170 | 0.278 / 0.330 | +| **simpleitk** | SimpleITK 2.5.6 | 9 | 23.3 | 0.156 / 0.156 | 0.654 / 0.641 | +| **synthmorph** | Docker `freesurfer/synthmorph` (TensorFlow, pretrained, `affine` model) | 9 | 27.1 | 0.167 / 0.171 | 0.282 / 0.320 | + +Per-case MI (nats, unclipped score; **bold** = best for the case): + +| case | antspyx | ants_trsaa | ants_reproa | ants_quickreproa | elastix | impact | niftyreg | simpleitk | synthmorph | +|---|---|---|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | **0.655** | 0.647 | 0.547 | 0.546 | 0.638 | 0.391 | 0.278 | 0.654 | 0.282 | +| m6_23XIMO4M_20230222 | **0.661** | 0.660 | 0.519 | 0.521 | 0.656 | 0.657 | 0.611 | 0.660 | 0.520 | +| m6_24V62QQM_20230223 | **0.710** | 0.708 | 0.568 | 0.569 | 0.684 | 0.708 | 0.254 | **0.710** | 0.252 | +| m6_24V62QQM_20231017 | 0.628 | **0.635** | 0.425 | 0.423 | 0.622 | 0.623 | 0.473 | 0.619 | 0.479 | +| m6_24V62QQM_20240820 | 0.670 | **0.673** | 0.459 | 0.464 | 0.670 | 0.659 | 0.121 | 0.545 | 0.122 | +| m6_26TYQI3C_20221027 | **0.719** | 0.717 | 0.525 | 0.523 | 0.675 | 0.676 | 0.593 | 0.699 | 0.570 | +| m6_2BBN5XA4_20250411 | **0.590** | 0.588 | 0.505 | 0.507 | 0.568 | 0.360 | 0.069 | **0.590** | 0.071 | +| m6_2C4TKI6H_20240506 | **0.646** | 0.645 | 0.546 | 0.542 | 0.636 | 0.451 | 0.241 | 0.635 | 0.256 | +| m6_2CDNN7BH_20240102 | **0.660** | 0.658 | 0.535 | 0.534 | 0.651 | 0.575 | 0.327 | 0.656 | 0.328 | + +Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) +for the case): + +| case | antspyx | ants_trsaa | ants_reproa | ants_quickreproa | elastix | impact | niftyreg | simpleitk | synthmorph | +|---|---|---|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | 0.1514 | 0.1516 | 0.1527 | 0.1528 | **0.1509** | 0.1607 | 0.1667 | 0.1513 | 0.1666 | +| m6_23XIMO4M_20230222 | 0.1565 | 0.1564 | 0.1598 | 0.1599 | 0.1565 | 0.1564 | 0.1585 | **0.1563** | 0.1616 | +| m6_24V62QQM_20230223 | **0.1466** | 0.1468 | 0.1527 | 0.1529 | 0.1470 | 0.1468 | 0.1670 | **0.1466** | 0.1672 | +| m6_24V62QQM_20231017 | 0.1519 | **0.1517** | 0.1599 | 0.1599 | 0.1527 | 0.1533 | 0.1587 | 0.1520 | 0.1581 | +| m6_24V62QQM_20240820 | 0.1560 | **0.1559** | 0.1618 | 0.1617 | 0.1564 | 0.1570 | 0.1840 | 0.1599 | 0.1843 | +| m6_26TYQI3C_20221027 | **0.1563** | 0.1570 | 0.1604 | 0.1615 | 0.1615 | 0.1657 | 0.1621 | 0.1585 | 0.1624 | +| m6_2BBN5XA4_20250411 | 0.1617 | 0.1617 | 0.1622 | 0.1622 | **0.1593** | 0.1662 | 0.1845 | 0.1617 | 0.1849 | +| m6_2C4TKI6H_20240506 | **0.1559** | 0.1561 | 0.1612 | 0.1613 | 0.1561 | 0.1594 | 0.1759 | **0.1559** | 0.1766 | +| m6_2CDNN7BH_20240102 | 0.1591 | 0.1590 | 0.1609 | 0.1617 | **0.1587** | **0.1587** | 0.1740 | **0.1587** | 0.1741 | + +### Interpretation + +- **The four intensity tools form a statistical tie, and affine barely + moves them.** `antspyx` (0.660), `ants_trsaa` (0.658), `simpleitk` + (0.654) and `elastix` (0.651) all sit within 0.009 nats of each other on + the median (ImpactLoss medians 0.156), and per-case MI bests split + antspyx 5/9, trsaa 2/9, simpleitk 2/9 (ties). Against the rigid + benchmark's per-pipeline medians, affine buys only + `antspyx` 0.652 → 0.660 (+0.008), `simpleitk` 0.647 → 0.654 (+0.007), + `elastix` 0.646 → 0.651 (+0.005) — the M6 head is essentially rigid + between sessions, so the extra 6 DOF mostly absorb inter-session + anisotropic scaling (≤ a few %), not misalignment. +- **`ants_trsaa` is the full TRSA cascade for free quality — and it is not + free.** Translation → Rigid → Similarity → Affine → Affine (5 linear + stages) reaches 0.658, i.e. the same as the single-stage `Affine` preset + (0.660) at 2.4× the wall time (30.4 s vs 12.6 s): the cascaded 6-DOF and + 9-DOF stages just re-solve what the one 12-DOF stage finds. +- **`ants_quickreproa` ≈ `ants_reproa` quality at 45 % of the time** + (0.523 vs 0.525 MI; 28.4 s vs 63.1 s): dropping the full-resolution rigid + level costs nothing. But both are clearly the weakest of the *top-tier + eligible* antspyx presets — the GC-based Rigid+Affine head (the linear + half of the Repro T1w-T1w workflow) locks a cross-modality-weak pose, the + same story as the Repro family in the deformable benchmark (0.530–0.549 + MI there). +- **`impact`'s affine run degrades vs its own rigid.** Median 0.623 + (mean 0.567) vs 0.652 for rigid: the semantic 12-DOF landscape has more + local minima, and the single CenterOfGravity start falls into two of them + (`m6_22IRRKU6` 0.391, `m6_2BBN5XA4` 0.360 — both 0.645 / 0.572 in the + rigid run). It also roughly doubles in wall time (215.5 s vs 111.2 s) + for the extra 6 DOF. +- **`niftyreg` full affine does not fix its gross-pose failures** (0.278 vs + 0.292 rigid; the two repositioned cases still at 0.069 / 0.121): reg_ + aladin's block-matching + TLS is a local optimizer either way, and the + output affine is still flipped by the orientation gate on 8/9 cases. +- **`synthmorph`'s `affine` model confirms the rigid-benchmark + investigation** (0.282 vs 0.277; within ±0.03 per case): the linear part + of its output is a near-pure rotation — the extra 6 DOF are essentially + unused, and the pretrained matrix head's gross-pose estimation weakness + stands for 12 DOF as for 6. +- **The deformable stage is still worth it.** The best affine median here + is 0.660 (`antspyx`); the best deformable in the companion benchmark is + 0.753 (`ants_rigid_synonly`; 0.723 for the stock `SyN` preset) — the + deformable stage adds +0.06 to +0.09 MI on top of the best affine, which + is fine-tissue alignment, not gross pose. +- **Production skull-ring QA (reference only, kept in the per-case JSONs):** + acceptance 4/9 `antspyx` / `ants_reproa` / `ants_quickreproa` / + `elastix` / `simpleitk`, 3/9 `ants_trsaa`, 2/9 `impact`, 1/9 `niftyreg` / + `synthmorph`. Concentricity: `ants_reproa` tightest mid-gap (5.6 mm), + `niftyreg` / `impact` / `synthmorph` worst on the rotational term + (23.5–25.3 mm std) — the same weak correlation with perceived alignment + noted in the other benchmarks, so it is reported for reference only. + +## Transform conventions (affine-specific) + +- **antspyx linear-only presets write one merged affine per consecutive + linear run.** The `TRSAA` cascade (5 linear stages) comes back as a single + *binary* `0GenericAffine.mat` that SimpleITK cannot parse; the pipeline + recovers it by point-probing the transform (exact for affine maps, + verified to NCC 1.0 against ANTs' own `warpedmovout` on the pilot cases). + `Repro[a]` / `QuickRepro[a]` (Rigid + Affine) merge the same way. +- **`niftyreg`'s `-aff` file** is documented as `Affine * Reference = + Floating`; the `fix_orientation` gate resolves the direction empirically + and inverts it on 8/9 M6 cases (as in the rigid benchmark). +- **`elastix` / `impact`** return the full combination transform (auto-init + folded in) in the contract direction (`auto_inverted: False`); the impact + affine read-back is `TransformParameters.0.itk.txt` (parsed by + SimpleITK, with a manual 12-parameter fallback). +- **`simpleitk`** returns the two-stage composite (rigid-init affine o + optimized affine) condensed to a single AffineTransform by point mapping. +- **`synthmorph`**'s `-T` LTA (LINEAR_RAS_TO_RAS 4×4) is the contract + transform directly, direction-checked by the usual gate. + +## Reproducing + +```bash +source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 +cd /mnt/b4/xfr/git26/longitudinal + +# all nine affine pipelines (results go to results/bench_registration_affine/; +# 2026-10-03 run, 6 workers, ~50 min wall dominated by impact's ~6 min phase) +python scripts/bench_registration.py --affine --pipelines "antspyx ants_trsaa ants_reproa ants_quickreproa elastix impact niftyreg simpleitk synthmorph" --workers 6 --threads 24 --qa-images + +# single pipeline / case subset +python scripts/bench_registration.py --affine --pipelines antspyx simpleitk --n-cases 3 + +# combine the per-pipeline JSONs into summary.json / summary.md / plot +python scripts/bench_registration.py --affine --merge +``` + +## Visual QA images + +`--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to +`results/bench_registration_affine/qa/__.png` (9 per +pipeline; the incomplete-CT case is skipped) — same format as the other +benchmarks: a 3×2 grid of the mid-axial / mid-coronal / mid-sagittal slices +through the head center, **left column** checkerboard blend of the warped CT +vs the T1c (clean registration = seamless mosaic), **right column** additive +color overlay (CT red, T1c green; agreement = yellow, misregistration = +red/green fringing). CT is brain-windowed 0–80 HU, T1c 2–98 % head window; +the title carries the pipeline, case key and the achieved Mattes MI. \ No newline at end of file diff --git a/REGISTRATION_BENCHMARK_DEFORMABLE.md b/REGISTRATION_BENCHMARK_DEFORMABLE.md index 843973d..51b29cc 100644 --- a/REGISTRATION_BENCHMARK_DEFORMABLE.md +++ b/REGISTRATION_BENCHMARK_DEFORMABLE.md @@ -245,10 +245,14 @@ for the case): ## ANTs variants: direct SyNRA / SyN and the Repro / QuickRepro workflows -The current `antspyx` row above is antspyx's stock **`SyN` preset** (Affine -2100×1200×1200×0 + SyN[0.2,3,0] 40×20×0, Mattes MI 32, center-of-mass init). -As a follow-up (2026-10-03, 6 workers as the main run) seven further -ANTs configurations were benchmarked on the same 10 cases and the same +The `ants_syn` row below is antspyx's stock **`SyN` preset** (Affine +2100×1200×1200×0 + SyN[0.2,3,0] 40×20×0, Mattes MI 32, center-of-mass +init) — the same configuration as the `antspyx` row of the main table +above, re-run in this follow-up batch and reproduced case-for-case within +0.002 nats (results in `ants_syn.json`), listed under its own name to +avoid confusing it with the `ants_synra` / `ants_rigid_synonly` variants. +As a follow-up (2026-10-03, 6 workers as the main run) six further ANTs +configurations were benchmarked on the same 10 cases and the same contract — the question: is there a better stock ANTs deformable workflow for CT→T1c than the `SyN` preset? @@ -256,7 +260,6 @@ for CT→T1c than the `SyN` preset? |---|---|---| | **ants_rigid_synonly** | legacy two-call (`Rigid` + `SyNOnly`) | Rigid[0.25] (Mattes 32, 500×500×200, 4×2×1, 3×2×0) → SyNOnly SyN[0.2,3,0] (Mattes 32, 100×100×0, 4×2×1, 2×1×0) initialized from the Rigid stage; `Composite.h5` read-back | | **ants_synra** | `SyNRA` | Rigid[0.25] + Affine[0.25] (each Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1) | -| **ants_syn** | `SyN` | Affine[0.25] (Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1); COM init — the stock preset, identical configuration to the existing `antspyx` row | | **ants_synreprosr** | `antsRegistrationSyNRepro[sr]` | Rigid[0.1] (GC 1/1/Regular/0.25, 1000×500×250×100, 8×4×2×1, 3×2×1×0vox) + SyN[0.2,3,0] (CC 1/2, 50×20, 2×1, 1×0vox) | | **ants_synreprobr** | `antsRegistrationSyNRepro[br]` | the `Repro` Rigid[0.1] stage + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×20, 2×1, 1×0vox) | | **ants_synquickreprosr** | `antsRegistrationSyNQuickRepro[sr]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + SyN[0.2,3,0] (CC 1/2, 50×0, 2×1) | @@ -280,16 +283,14 @@ linear stages (`SyNRA`), it writes the merged linear stage as a *binary* by point-probing the transform (exact for affine maps, verified to 2×10⁻⁴ mm). -All six variants registered 9/9 (the 10th is the known incomplete-CT data -issue). **Bold** = best for the column; `antspyx` is the existing row, -re-listed for reference. +All seven rows registered 9/9 (the 10th is the known incomplete-CT data +issue). **Bold** = best for the column. | pipeline | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) | |---|---|---|---|---| | **ants_rigid_synonly** | 9 | 208.4 | 0.156 / 0.155 | **0.753 / 0.763** | -| **antspyx** (existing) | 9 | 31.7 | 0.158 / 0.156 | 0.723 / 0.729 | -| **ants_synra** | 9 | 43.6 | 0.158 / 0.156 | 0.722 / 0.729 | | **ants_syn** | 9 | 33.8 | 0.158 / 0.156 | 0.722 / 0.729 | +| **ants_synra** | 9 | 43.6 | 0.158 / 0.156 | 0.722 / 0.729 | | **ants_synreprosr** | 9 | 294.7 | **0.151 / 0.153** | 0.530 / 0.529 | | **ants_synreprobr** | 9 | 370.3 | 0.156 / 0.157 | 0.541 / 0.537 | | **ants_synquickreprosr** | 9 | 83.7 | 0.154 / 0.156 | 0.542 / 0.534 | @@ -297,61 +298,56 @@ re-listed for reference. Per-case MI (nats, unclipped score; **bold** = best for the case): -| case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | -|---|---|---|---|---|---|---|---|---| -| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | -| m6_22IRRKU6_20240115 | 0.679 | **0.741** | 0.681 | 0.681 | 0.623 | 0.646 | 0.630 | 0.627 | -| m6_23XIMO4M_20230222 | 0.723 | **0.760** | 0.722 | 0.722 | 0.530 | 0.544 | 0.529 | 0.530 | -| m6_24V62QQM_20230223 | 0.789 | **0.824** | 0.789 | 0.787 | 0.520 | 0.528 | 0.542 | 0.549 | -| m6_24V62QQM_20231017 | **0.720** | 0.716 | **0.720** | **0.720** | 0.473 | 0.475 | 0.469 | 0.462 | -| m6_24V62QQM_20240820 | 0.738 | **0.753** | 0.738 | 0.738 | 0.524 | 0.525 | 0.517 | 0.524 | -| m6_26TYQI3C_20221027 | 0.799 | **0.857** | 0.798 | 0.799 | 0.557 | 0.541 | 0.599 | 0.632 | -| m6_2BBN5XA4_20250411 | 0.665 | **0.696** | 0.665 | 0.666 | 0.230 | 0.303 | 0.235 | 0.297 | -| m6_2C4TKI6H_20240506 | 0.719 | **0.751** | 0.719 | 0.719 | 0.672 | 0.649 | 0.654 | 0.621 | -| m6_2CDNN7BH_20240102 | 0.726 | **0.765** | 0.726 | 0.725 | 0.629 | 0.622 | 0.631 | 0.599 | +| case | ants_syn | ants_rigid_synonly | ants_synra | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | +|---|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | 0.681 | **0.741** | 0.681 | 0.623 | 0.646 | 0.630 | 0.627 | +| m6_23XIMO4M_20230222 | 0.722 | **0.760** | 0.722 | 0.530 | 0.544 | 0.529 | 0.530 | +| m6_24V62QQM_20230223 | 0.787 | **0.824** | 0.789 | 0.520 | 0.528 | 0.542 | 0.549 | +| m6_24V62QQM_20231017 | **0.720** | 0.716 | **0.720** | 0.473 | 0.475 | 0.469 | 0.462 | +| m6_24V62QQM_20240820 | 0.738 | **0.753** | 0.738 | 0.524 | 0.525 | 0.517 | 0.524 | +| m6_26TYQI3C_20221027 | 0.799 | **0.857** | 0.798 | 0.557 | 0.541 | 0.599 | 0.632 | +| m6_2BBN5XA4_20250411 | 0.666 | **0.696** | 0.665 | 0.230 | 0.303 | 0.235 | 0.297 | +| m6_2C4TKI6H_20240506 | 0.719 | **0.751** | 0.719 | 0.672 | 0.649 | 0.654 | 0.621 | +| m6_2CDNN7BH_20240102 | 0.725 | **0.765** | 0.726 | 0.629 | 0.622 | 0.631 | 0.599 | Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) for the case): -| case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | -|---|---|---|---|---|---|---|---|---| -| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | -| m6_22IRRKU6_20240115 | 0.1541 | 0.1530 | 0.1542 | 0.1540 | **0.1442** | 0.1482 | 0.1474 | 0.1493 | -| m6_23XIMO4M_20230222 | 0.1581 | 0.1565 | 0.1581 | 0.1581 | **0.1505** | 0.1549 | 0.1535 | 0.1585 | -| m6_24V62QQM_20230223 | 0.1473 | **0.1464** | 0.1472 | 0.1472 | 0.1515 | 0.1553 | 0.1535 | 0.1560 | -| m6_24V62QQM_20231017 | 0.1518 | **0.1512** | 0.1516 | 0.1519 | 0.1559 | 0.1575 | 0.1577 | 0.1581 | -| m6_24V62QQM_20240820 | **0.1545** | **0.1545** | 0.1547 | 0.1546 | 0.1553 | 0.1578 | 0.1586 | 0.1586 | -| m6_26TYQI3C_20221027 | 0.1580 | 0.1562 | 0.1578 | 0.1579 | **0.1561** | 0.1574 | 0.1602 | 0.1610 | -| m6_2BBN5XA4_20250411 | 0.1624 | **0.1611** | 0.1624 | 0.1624 | 0.1680 | 0.1695 | 0.1690 | 0.1707 | -| m6_2C4TKI6H_20240506 | 0.1594 | 0.1582 | 0.1594 | 0.1593 | **0.1490** | 0.1526 | 0.1517 | 0.1545 | -| m6_2CDNN7BH_20240102 | 0.1590 | 0.1582 | 0.1589 | 0.1589 | **0.1493** | 0.1560 | 0.1538 | 0.1580 | +| case | ants_syn | ants_rigid_synonly | ants_synra | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | +|---|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | 0.1540 | 0.1530 | 0.1542 | **0.1442** | 0.1482 | 0.1474 | 0.1493 | +| m6_23XIMO4M_20230222 | 0.1581 | 0.1565 | 0.1581 | **0.1505** | 0.1549 | 0.1535 | 0.1585 | +| m6_24V62QQM_20230223 | 0.1472 | **0.1464** | 0.1472 | 0.1515 | 0.1553 | 0.1535 | 0.1560 | +| m6_24V62QQM_20231017 | 0.1519 | **0.1512** | 0.1516 | 0.1559 | 0.1575 | 0.1577 | 0.1581 | +| m6_24V62QQM_20240820 | 0.1546 | **0.1545** | 0.1547 | 0.1553 | 0.1578 | 0.1586 | 0.1586 | +| m6_26TYQI3C_20221027 | 0.1579 | 0.1562 | 0.1578 | **0.1561** | 0.1574 | 0.1602 | 0.1610 | +| m6_2BBN5XA4_20250411 | 0.1624 | **0.1611** | 0.1624 | 0.1680 | 0.1695 | 0.1690 | 0.1707 | +| m6_2C4TKI6H_20240506 | 0.1593 | 0.1582 | 0.1594 | **0.1490** | 0.1526 | 0.1517 | 0.1545 | +| m6_2CDNN7BH_20240102 | 0.1589 | 0.1582 | 0.1589 | **0.1493** | 0.1560 | 0.1538 | 0.1580 | ### Interpretation (ANTs variants) -- **`ants_rigid_synonly` is the best ANTs configuration of the eight** — - the only one that beats the current `antspyx` row: median MI **0.753** +- **`ants_rigid_synonly` is the best ANTs configuration of the seven** — + the only one that beats the current `ants_syn` row: median MI **0.753** (per-case best on 8 of 9; its 0.857 on `m6_26TYQI3C` is the highest MI of the whole benchmark, ahead of the previous 0.801) at ImpactLoss 0.156 (second-best median after `synthmorph`'s 0.156). The gain over the `SyN` preset comes from the *deeper* SyN stage: 100 iterations at every level **including full resolution** (100×100×0 vs 40×20×0) - refining the Rigid stage's pose. The cost is 208.4 s/case (~6.6× the - preset) — the full-resolution SyN iterations dominate. -- **`ants_synra` is indistinguishable from the current `antspyx` row** on - quality (median MI 0.722 vs 0.723, ImpactLoss 0.158 vs 0.158; per-case - deltas within 0.002 nats, ties on 4 of 9) at **+38 % wall time** - (43.6 s vs 31.7 s): the explicit Rigid stage ahead of the Affine stage + refining the Rigid stage's pose. The cost is 208.4 s/case (~6.2× the + `ants_syn` time) — the full-resolution SyN iterations dominate. +- **`ants_synra` is indistinguishable from the `ants_syn` row** on + quality (median MI 0.722 vs 0.722, ImpactLoss 0.158 vs 0.158; per-case + deltas within 0.002 nats, identical on 6 of 9) at **+29 % wall time** + (43.6 s vs 33.8 s): the explicit Rigid stage ahead of the Affine stage buys nothing — on these pairs the 12-DOF stage from the center-of-mass init already finds the gross pose, so the extra 6-DOF stage just re-solves it. -- **`ants_syn` reproduces the existing `antspyx` row** (median MI 0.722 - vs 0.723, ImpactLoss 0.158, s/reg 33.8 vs 31.7; per-case deltas within - 0.002 nats) — as expected, since it is the *same stock `SyN` preset* - (Affine + SyN[0.2,3,0], COM init) re-run in this follow-up batch: it - serves as a consistency check that the shared read-back / scoring path - is configuration-stable. - **The Repro / QuickRepro family is semantically competitive but weaker on - intensity:** median MI 0.530–0.549 (vs 0.723 for the `SyN` preset), + intensity:** median MI 0.530–0.549 (vs 0.722 for the `SyN` preset, i.e. + `ants_syn`), driven by the GC-based Rigid stage and the short 2-level SyN/BSplineSyN refinement (50×20 or 50×0 iters, 2×1 pyramid, CC 2-bin metric) — a workflow tuned for intra-modality T1w-T1w. `ants_synreprosr` still takes @@ -363,10 +359,10 @@ the case): slowest ANTs configuration (370 s/case) at the same (low) MI as the 84–96 s `QuickRepro` variants (which drop the 4th rigid level and the 2nd SyN level, changing the MI by < 0.02). Only `ants_rigid_synonly` - challenges the `SyN` preset's quality (better MI at 6.6× the time); - every other variant is worse for more or for the same. + challenges `ants_syn`'s quality (better MI at ≈6× the time); every + other variant is worse for more or for the same. - **Production skull-ring QA (reference only):** acceptance 4/9 - `ants_syn` / `ants_synra` / `ants_rigid_synonly` (= baseline), 2/9 + `ants_syn` / `ants_synra` / `ants_rigid_synonly`, 2/9 `ants_synreprobr`, 1/9 the rest; rotational dispersion worst in the `Repro` / `QuickRepro` family (24–29 mm std) vs 14–15 mm for `ants_syn` / `ants_synra` / `ants_rigid_synonly` — consistent with diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index d7e91c7..04365b3 100644 --- a/scripts/bench_registration.py +++ b/scripts/bench_registration.py @@ -293,6 +293,28 @@ def make_qa_png(warped_ct, mr, out_png, title, block=24, ct_brain_window=(0.0, 8 def _call(pipeline, ct, t1c, tv, args): + if args.affine: + if pipeline in ("antspyx", "ants_trsaa", "ants_reproa", "ants_quickreproa"): + return RP.register_ants_affine_variant(ct, t1c, preset_name=pipeline) + if pipeline == "elastix": + return RP.register_elastix_affine(ct, t1c, workdir=args.elastix_dir, + threads=args.threads) + if pipeline == "impact": + return RP.register_impact_affine(ct, t1c, gpu=args.impact_gpu, + iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS, + samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES, + threads=args.threads) + if pipeline == "niftyreg": + return RP.register_niftyreg_affine(ct, t1c, threads=args.threads, + backend=args.niftyreg_backend, + gpu=args.niftyreg_gpu) + if pipeline == "simpleitk": + return RP.register_simpleitk_affine(ct, t1c, threads=args.threads) + if pipeline == "synthmorph": + return RP.register_synthmorph_affine(ct, t1c, gpu=args.synthmorph_gpu, + threads=args.threads, + extent=args.synthmorph_extent) + raise KeyError(pipeline) if args.deformable: if pipeline == "simpleitk": return RP.register_simpleitk_deformable(ct, t1c, threads=args.threads) @@ -526,9 +548,24 @@ def merge(args): raise SystemExit(f"no pipeline results found in {out_dir}") print(f"[merge] {len(rows)} pipeline result files in {out_dir}") print_table(rows) - deformable = (os.path.basename(os.path.normpath(args.out_dir)) - == "bench_registration_deformable") - if deformable: + out_base = os.path.basename(os.path.normpath(args.out_dir)) + affine = out_base == "bench_registration_affine" + deformable = out_base == "bench_registration_deformable" + if affine: + protocol = ("affine CT->T1c; cross-modality; each toolkit estimates a " + "full affine (12-DOF) transform with its native workflow " + "(antspyx: stock Affine preset + TRSAA / Repro[a] / " + "QuickRepro[a] linear-only presets; elastix AffineTransform; " + "impact MIND+TS reference map with AffineTransform; " + "niftyreg reg_aladin full affine; simpleitk AffineTransform; " + "synthmorph affine matrix model); alignment scored by " + "achieved Mattes MI (nats, warped-CT ~ MR over the MR head, " + "higher = better) and the IMPACT semantic loss " + "(vboussot/ImpactLoss, TS/M730 features, warped-CT ~ MR over " + "the full volumes, lower = better); reg time excludes " + "I/O+scoring; per-case JSON retains the production " + "skull-ring gap / acceptance QA") + elif deformable: protocol = ("deformable CT->T1c; cross-modality; each pipeline runs its " "toolkit's native deformable workflow (linear init + deformable " "stage, or a single deformable stage where that is native); the " @@ -559,7 +596,8 @@ def merge(args): json.dump(summary, f, indent=1) print("summary:", sf) write_plot(rows, os.path.join(out_dir, "bench_registration.png"), - label="deformable" if deformable else "rigid") + label="deformable" if deformable + else ("affine" if affine else "rigid")) # markdown table for the report md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |", "|---|---|---|---|---|---|---|---|---|"] @@ -582,6 +620,12 @@ def merge(args): def parse_args(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--affine", action="store_true", + help="run the affine (12-DOF) variants of the pipelines " + "(each toolkit's affine workflow + the antspyx " + "Affine / TRSAA / Repro[a] / QuickRepro[a] presets; " + "results go to results/bench_registration_affine " + "by default)") ap.add_argument("--deformable", action="store_true", help="run the deformable variants of the pipelines " "(each toolkit's native deformable workflow; results go " @@ -650,6 +694,9 @@ def main(): if args.deformable: if args.out_dir == os.path.join(ROOT, "results", "bench_registration"): args.out_dir = os.path.join(ROOT, "results", "bench_registration_deformable") + if args.affine: + if args.out_dir == os.path.join(ROOT, "results", "bench_registration"): + args.out_dir = os.path.join(ROOT, "results", "bench_registration_affine") if args.merge: merge(args) return diff --git a/scripts/qa_label_interpolation.py b/scripts/qa_label_interpolation.py new file mode 100644 index 0000000..e856ea1 --- /dev/null +++ b/scripts/qa_label_interpolation.py @@ -0,0 +1,436 @@ +"""Visual QA of TV label warping: registration method x label interpolation method. + +For one M6 case (default: 2BBN5XA4 / 20250411): + - fixed = T1c (MR, picked by scripts/02_build_m6_dataset.pick_t1c), + moving = CT (RT/ct_image.nii.gz); the TV (RT/TV-*.nii.gz) lives in the CT grid. + - registration via ANTsPy presets (src/registration_pipelines.py), each returning + a transform T that maps T1c -> CT, so sitk.Resample(tv, t1c, T, ...) places + the TV in MR (T1c) coordinates: + ants_rigid : Rigid[0.25] (Mattes, 500x500x200, 4x2x1) + ants_affine : Affine[0.25] (Mattes, 2100x1200x1200x10, 6x4x2x1) + ants_syn : SyN (Affine[0.25] + SyN[0.2,3,0] stages) + - the binary TV is then warped to the T1c grid by nearest-neighbor (nn) and + by SignedMaurerDistanceMap (mm, inside<0) -> {linear, B-spline order-3} + interp -> threshold (smd <= thr, default 0.5) + - one PNG per registration: 3 rows (axial / coronal / sagittal) x 4 columns + (leftmost = bare MR, then nn + each SMD method at its threshold); all + columns show the SAME slice + (the nn warp's largest-TV cross-section per view), ZOOMED to the TV + bounding box + margin (the NN box, shared by all columns); only the + 1-voxel TV contour is drawn (yellow, no fill), plus the volume + + % of the original CT-grid TV volume in the column header. + +Outputs (default results/label_interp_qa/_/): + _T.h5 registration transform (cached; re-runs skip it) + ___tv.nii.gz warped TV on the T1c grid (warp_key includes thr) + __label_qa_.png one 3x5 figure per registration + summary.json per-combination metrics (nvox, mm3, dice, shifts) + + the SMD thresholding volume-preservation check + +Usage (repo root): + source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 + python scripts/qa_label_interpolation.py + python scripts/qa_label_interpolation.py --subject 2BBN5XA4 --date 20250411 \ + --regs ants_rigid ants_affine ants_syn +""" +import argparse +import json +import os +import sys +import time + +ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, ROOT) + +import numpy as np +import SimpleITK as sitk +from scipy import ndimage + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.colors import ListedColormap + +from src.common import d, read_nii_img +from src import registration_pipelines as RP + +M6 = RP._load_m6() + +REGS = ["ants_rigid", "ants_affine", "ants_syn"] +SMD_THRESHOLDS = [0.5] # SMD mask = (resampled smd <= thr); 0.0 / 1.0 were +# tested earlier (36% under / 40% over volume) and dropped from the figure +METHOD_LABELS = {"nn": "NN", "smd_linear": "SMD + linear", "smd_bspline": "SMD + B-spline(3)"} +VIEWS = ["axial", "coronal", "sagittal"] +SMD_DEFAULT = 1e5 # out-of-grid fill for the SMD volume: far OUTSIDE (inside < 0) +# figure / warp columns: NN, then each SMD method at each threshold +COLS = ([("nn", None)] + + [(m, t) for m in ("smd_linear", "smd_bspline") for t in SMD_THRESHOLDS]) +ZOOM_MARGIN = 12 # voxels of context around the TV bbox in the displayed slices + + +def warp_key(method, thr): + return method if thr is None else f"{method}_thr{thr:g}" + + +def col_label(method, thr): + if thr is None: + return METHOD_LABELS[method] + return f"{METHOD_LABELS[method]} (smd<={thr:g})" + + +def reg_call(reg, ct, t1c): + """ANTsPy registration, fixed = T1c / moving = CT (the shared RP rule).""" + if reg == "ants_rigid": + return RP.register_antspyx(ct, t1c) + if reg == "ants_affine": + return RP.register_ants_affine_variant(ct, t1c, "antspyx") + if reg == "ants_syn": + return RP.register_ants_variant_deformable(ct, t1c, "ants_syn") + raise KeyError(f"unknown registration {reg!r} (one of {REGS})") + + +def find_case(subject, date): + labeled, _unlabeled = M6.discover() + for r in labeled: + if r.get("subject") == subject and r.get("date") == date: + return r + raise SystemExit(f"case {subject}/{date} not found among m6 labeled rows") + + +def get_transform(reg, ct, t1c, case_dir, use_cache): + """Register (or load the cached transform), return (T, meta).""" + tf = os.path.join(case_dir, reg + "_T.h5") + if use_cache and os.path.exists(tf): + print(f"[{reg}] loading cached transform {tf}", flush=True) + return sitk.ReadTransform(tf), {"cached": True} + t0 = time.time() + T, meta = reg_call(reg, ct, t1c) + dt = time.time() - t0 + sitk.WriteTransform(T, tf) + print(f"[{reg}] registered in {dt:.0f}s ({meta.get('toolkit')}) " + f"inverted={meta.get('auto_inverted')} gap={meta.get('skull_ring_gap_mm')}", + flush=True) + return T, {**meta, "seconds_wall": round(dt, 1), "cached": False} + + +def make_tv_bin(tv): + """Binary 0/1 TV image on the TV's own (CT) grid.""" + a = (sitk.GetArrayFromImage(tv) > 0.5).astype(np.uint8) + img = sitk.GetImageFromArray(a) + img.CopyInformation(tv) + return img + + +def warp_tv(tv_bin, t1c, T, method, thr=0.0): + """Warp the binary TV onto the t1c grid with T (T1c -> CT). Returns a + (z, y, x) uint8 0/1 mask (the interpolated label in MR coordinates). + For the SMD methods the threshold is mask = (smd <= thr).""" + rs = sitk.ResampleImageFilter() + rs.SetReferenceImage(t1c) + rs.SetTransform(T) + if method == "nn": + rs.SetInterpolator(sitk.sitkNearestNeighbor) + rs.SetDefaultPixelValue(0) + rs.SetOutputPixelType(sitk.sitkUInt8) + return (sitk.GetArrayFromImage(rs.Execute(tv_bin)) > 0).astype(np.uint8) + # signed distance to the TV boundary in mm: negative inside, positive + # outside. Note the SimpleITK 2.5 signature is + # (image1, insideIsPositive=False, squaredDistance=True, useImageSpacing=False, + # backgroundValue=0.0) -- keywords are mandatory: the defaults would give + # a SQUARED distance map, which interpolates wrongly. On this build the + # boundary (outer shell) voxels of the foreground get exactly 0, so the + # mask is recovered with <= 0, not < 0. + smd = sitk.SignedMaurerDistanceMap( + tv_bin, insideIsPositive=False, squaredDistance=False, useImageSpacing=True) + smd = sitk.Cast(smd, sitk.sitkFloat32) + if method == "smd_linear": + rs.SetInterpolator(sitk.sitkLinear) + else: + rs.SetInterpolator(sitk.sitkBSpline3) # B-spline order 3 + rs.SetDefaultPixelValue(SMD_DEFAULT) + rs.SetOutputPixelType(sitk.sitkFloat32) + return (sitk.GetArrayFromImage(rs.Execute(smd)) <= thr).astype(np.uint8) + + +def best_slices(mask, shape): + """(idx, area) of the largest TV cross-section per orthogonal view.""" + if not mask.any(): + return {"axial": (shape[0] // 2, 0), + "coronal": (shape[1] // 2, 0), + "sagittal": (shape[2] // 2, 0)} + az = mask.sum(axis=(1, 2)) # per z (axial) + ay = mask.sum(axis=(0, 2)) # per y (coronal) + ax = mask.sum(axis=(0, 1)) # per x (sagittal) + return {"axial": (int(np.argmax(az)), int(az.max())), + "coronal": (int(np.argmax(ay)), int(ay.max())), + "sagittal": (int(np.argmax(ax)), int(ax.max()))} + + +def tv_crop_box(mask, shape, margin): + """(lo, hi) 3-D index box around the mask (per axis of (z, y, x)) with a + margin; centered box when the mask is empty.""" + shape = np.array(shape) + if mask.any(): + idx = np.argwhere(mask) + lo, hi = idx.min(0), idx.max(0) + else: + lo = hi = shape // 2 + lo = np.clip(lo - margin, 0, shape - 1) + hi = np.clip(hi + margin, 0, shape - 1) + return lo, hi + + +def view_crop(arr, kind, idx, lo_b, hi_b): + """Zoomed slice display around the TV box (bench_registration.make_qa_png + convention: coronal / sagittal flipped so superior is up).""" + if kind == "axial": + s = arr[idx, lo_b[1]:hi_b[1] + 1, lo_b[2]:hi_b[2] + 1] + elif kind == "coronal": + s = arr[lo_b[0]:hi_b[0] + 1, idx, lo_b[2]:hi_b[2] + 1][::-1, :] + else: + s = arr[lo_b[0]:hi_b[0] + 1, lo_b[1]:hi_b[1] + 1, idx][::-1, :] + return np.ascontiguousarray(s) + + +def dice3(a, b): + a, b = a > 0, b > 0 + if not a.any() and not b.any(): + return 1.0 + return float(2 * (a & b).sum() / max(a.sum() + b.sum(), 1)) + + +def centroid_mm(mask, spacing): + c = np.array(ndimage.center_of_mass(mask)) * np.array(spacing) + return c + + +def check_volume(tv_bin, t1c, source_vol_mm3): + """Does thresholding the SMD preserve the TV volume? Resampled with the + IDENTITY transform (CT grid -> MR grid, no registration), which isolates + the interpolation + thresholding step from registration error. One row + per figure column (NN + each SMD method x threshold).""" + T = sitk.Transform(3, sitk.sitkIdentity) + sp = [float(s) for s in t1c.GetSpacing()] + vv = sp[0] * sp[1] * sp[2] + out = {} + for method, thr in COLS: + mask = warp_tv(tv_bin, t1c, T, method, thr or 0.0) + out[warp_key(method, thr)] = round(float(mask.sum() * vv), 1) + print(f"volume preservation (identity resample CT grid -> MR grid): " + f"source TV = {source_vol_mm3:.0f} mm3", flush=True) + for tag, v in out.items(): + print(f" {tag:20s}: {v:8.1f} mm3 ({100.0 * v / source_vol_mm3:5.1f}%)", + flush=True) + return out + + +def make_png(mr, combos, out_png, title, spacing, window, source_vol): + """One registration: 3 rows (views) x (len(combos)+1) columns; the + leftmost column is the bare MR (no overlay). All columns show the SAME + slices (the NN warp's largest-area slice per view); volumes go in the + column headers (row-independent): the original CT TV volume under the + bare-MR column, "% of it" under each warp column.""" + lo, hi = window + ncol, nrow = len(combos) + 1, len(VIEWS) + cw, ch, gap, m = 3.6, 3.4, 0.15, 0.62 + top = 1.5 + fig_w = ncol * cw + (ncol - 1) * gap + 2 * m + fig_h = nrow * ch + (nrow - 1) * gap + m + top + fig = plt.figure(figsize=(fig_w, fig_h)) + fig.patch.set_facecolor("k") + aspects = {"axial": spacing[1] / spacing[0], + "coronal": spacing[2] / spacing[0], + "sagittal": spacing[2] / spacing[1]} + cmap_edge = ListedColormap(["#ffee00"]) + + def w(a): + return np.clip((a - lo) / max(hi - lo, 1e-9), 0.0, 1.0) + + for r, kind in enumerate(VIEWS): + for c in range(ncol): + x0 = (m + c * (cw + gap)) / fig_w + y0 = (m + (nrow - 1 - r) * (ch + gap)) / fig_h + axi = fig.add_axes([x0, y0, cw / fig_w, ch / fig_h]) + # every column (incl. the bare-MR one) shows the same slice/box + ref = combos[0]["views"][kind] + axi.imshow(w(ref["mr"]), cmap="gray", aspect=aspects[kind]) + if c > 0: + v = combos[c - 1]["views"][kind] + m2 = (v["mask"] > 0) + if m2.any(): + cont = m2 & ~ndimage.binary_erosion(m2) + axi.imshow(np.ma.masked_where(~cont, cont), cmap=cmap_edge, + aspect=aspects[kind]) + else: + axi.text(0.5, 0.5, "empty warp", transform=axi.transAxes, + ha="center", va="center", color="r", fontsize=8) + axi.text(0.02, 0.02, f"{kind}={ref['idx']}", transform=axi.transAxes, + color="w", fontsize=7, + bbox=dict(facecolor="k", alpha=0.55, pad=1.0)) + axi.set_xticks([]); axi.set_yticks([]) + for s_ in axi.spines.values(): + s_.set_edgecolor("0.4") + # column headers: bare MR + original CT TV volume, then method + % volume + cx0 = (m + cw / 2) / fig_w + fig.text(cx0, 1 - 0.74 / fig_h, "MR (no overlay)", color="#ffffff", + fontsize=9.5, ha="center") + fig.text(cx0, 1 - 1.02 / fig_h, f"CT TV: {source_vol:.0f} mm3", + color="#b6f0b6", fontsize=9, ha="center") + for c, comb in enumerate(combos): + cx = (m + (c + 1) * (cw + gap) + cw / 2) / fig_w + fig.text(cx, 1 - 0.74 / fig_h, comb["ilabel"], color="#9fd0ff", + fontsize=9.5, ha="center") + fig.text(cx, 1 - 1.02 / fig_h, + f"{comb['vol_mm3']:.0f} mm3 ({comb['pct_of_source']:.0f}%)", + color="#b6f0b6", fontsize=9, ha="center") + # view labels (left) + for r, kind in enumerate(VIEWS): + y0 = (m + (nrow - 1 - r) * (ch + gap)) / fig_h + fig.text(0.012, y0 + ch / 2 / fig_h, kind.capitalize(), color="w", fontsize=10, + ha="center", va="center", rotation=90) + fig.suptitle(title, color="w", fontsize=11, y=1 - 0.16 / fig_h) + fig.savefig(out_png, facecolor="k", dpi=100) + plt.close(fig) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--subject", default="2BBN5XA4") + ap.add_argument("--date", default="20250411") + ap.add_argument("--regs", nargs="*", default=REGS, + help="subset of: " + " ".join(REGS)) + ap.add_argument("--tv", default=None, + help="explicit TV file (default: first TV-*.nii.gz in RT)") + ap.add_argument("--out", default=None, + help="output dir (default results/label_interp_qa/_)") + ap.add_argument("--no-cache", action="store_true", + help="re-register and re-warp even if cached outputs exist") + args = ap.parse_args() + + row = find_case(args.subject, args.date) + case_dir = args.out or os.path.join(ROOT, "results", "label_interp_qa", + f"{args.subject}_{args.date}") + os.makedirs(case_dir, exist_ok=True) + ct_path = row["ct"] + t1c_path = row["img"] + tv_path = args.tv or row["label"] + + print(f"case {row['key']}: ct={os.path.basename(ct_path)} " + f"t1c={os.path.basename(t1c_path)} tv={os.path.basename(tv_path)}", flush=True) + ct = read_nii_img(ct_path) + t1c = read_nii_img(t1c_path) + tv = read_nii_img(tv_path) + tv_bin = make_tv_bin(tv) + mr = sitk.GetArrayFromImage(t1c).astype("float32") + spacing = [float(s) for s in t1c.GetSpacing()] # (x, y, z) + print(f"t1c {t1c.GetSize()} @ {spacing}mm, tv nvox={int((sitk.GetArrayFromImage(tv_bin) > 0).sum())} " + f"({int((sitk.GetArrayFromImage(tv_bin) > 0).sum()) * ct.GetSpacing()[0] * ct.GetSpacing()[1] * ct.GetSpacing()[2]:.0f} mm3 on the CT grid)", + flush=True) + + # MR display window from the head + pos = mr > 0 + thr = np.percentile(mr[pos], 2) if pos.any() else 0.0 + head = mr > max(thr, 1e-6) + lo, hi = (np.percentile(mr[head], (2, 98)) if head.any() else (0.0, 1.0)) + + # volume-preservation check for the SMD thresholding (no registration) + sp_ct = [float(s) for s in tv_bin.GetSpacing()] + source_vol = float((sitk.GetArrayFromImage(tv_bin) > 0).sum()) * sp_ct[0] * sp_ct[1] * sp_ct[2] + vol_check = check_volume(tv_bin, t1c, source_vol) + + combos = [] + reg_meta = {} + for reg in args.regs: + T, meta = get_transform(reg, ct, t1c, case_dir, use_cache=not args.no_cache) + reg_meta[reg] = meta + reg_combos = [] + # warp every column (NN + each SMD method x threshold) + for method, thr in COLS: + key = warp_key(method, thr) + p = os.path.join(case_dir, f"{reg}__{key}_tv.nii.gz") + if not args.no_cache and os.path.exists(p): + print(f"[{reg}/{key}] loading cached warp {p}", flush=True) + mask = (sitk.GetArrayFromImage(read_nii_img(p)) > 0).astype(np.uint8) + dt = 0.0 + else: + t0 = time.time() + mask = warp_tv(tv_bin, t1c, T, method, thr or 0.0) + dt = time.time() - t0 + mimg = sitk.GetImageFromArray(mask) + mimg.CopyInformation(t1c) + sitk.WriteImage(mimg, p, True) + nvox = int(mask.sum()) + reg_combos.append({"reg": reg, "method": method, "thr": thr, + "key": key, "ilabel": col_label(method, thr), + "mask": mask, "nvox": nvox, + "vol_mm3": round(nvox * spacing[0] * spacing[1] * spacing[2], 1), + "seconds_warp": round(dt, 1)}) + print(f"[{reg}/{key}] nvox={nvox} ({reg_combos[-1]['vol_mm3']} mm3) ({dt:.1f}s)", + flush=True) + # all columns show the SAME slices (the NN warp's largest-area slice + # per view) and the SAME zoom box (NN mask bbox + margin) + nn = next(c for c in reg_combos if c["method"] == "nn") + sl = best_slices(nn["mask"], mr.shape) + lo_b, hi_b = tv_crop_box(nn["mask"], mr.shape, ZOOM_MARGIN) + for c in reg_combos: + c["slices"] = [sl[k][0] for k in VIEWS] + c["views"] = {kind: {"mr": view_crop(mr, kind, sl[kind][0], lo_b, hi_b), + "mask": view_crop(c["mask"], kind, sl[kind][0], lo_b, hi_b), + "idx": sl[kind][0], "area": sl[kind][1]} for kind in VIEWS} + c["pct_of_source"] = 100.0 * c["vol_mm3"] / source_vol + c["pct_vs_nn"] = (100.0 * c["vol_mm3"] / nn["vol_mm3"] + if nn["vol_mm3"] else 0.0) + if c["method"] == "nn": + c["dice_vs_nn"] = 1.0 + c["centroid_shift_mm_vs_nn"] = 0.0 + else: + c["dice_vs_nn"] = round(dice3(nn["mask"], c["mask"]), 4) + cs = centroid_mm(c["mask"], spacing) - centroid_mm(nn["mask"], spacing) + c["centroid_shift_mm_vs_nn"] = [round(float(x), 2) for x in cs] + combos.extend(reg_combos) + # warped-volume table per registration (% of the original CT-grid TV volume) + parts = " | ".join(f"{c['key']}={c['vol_mm3']:.0f} mm3 ({c['pct_of_source']:.0f}%)" + for c in reg_combos) + print(f"warped volumes {reg} (CT TV={source_vol:.0f} mm3, " + f"slices ax/cor/sag={nn['slices']}): {parts}", flush=True) + + pngs = {} + for reg in args.regs: + title = (f"{args.subject} {args.date} -- {reg}: TV warping QA " + f"(same slices = NN largest-area per view)\n" + f"fixed=T1c ({os.path.basename(t1c_path)[:40]}...) moving=CT | " + f"yellow = TV contour | MR window {lo:.0f}..{hi:.0f}") + png = os.path.join(case_dir, f"{args.subject}_{args.date}_label_qa_{reg}.png") + make_png(mr, [c for c in combos if c["reg"] == reg], png, title, spacing, + (float(lo), float(hi)), source_vol) + pngs[reg] = png + print(f"png: {png}", flush=True) + + summary = { + "case": row["key"], "subject": args.subject, "date": args.date, + "inputs": {"ct": ct_path, "t1c": t1c_path, "tv": tv_path}, + "t1c_geometry": {"size": list(t1c.GetSize()), "spacing": spacing, + "origin": [float(o) for o in t1c.GetOrigin()]}, + "regs": {k: {kk: vv for kk, vv in v.items() if kk != "stages"} + for k, v in reg_meta.items()}, + "volume_check": {"source_tv_mm3": round(source_vol, 1), + "identity_resample_mm3": vol_check}, + "pngs": pngs, + "combos": {f"{c['reg']}/{c['key']}": { + "nvox": c["nvox"], "vol_mm3": c["vol_mm3"], + "pct_of_source_tv": round(c["pct_of_source"], 1), + "slices_ax_cor_sag": c["slices"], + "dice_vs_nn": c["dice_vs_nn"], + "centroid_shift_mm_vs_nn": c["centroid_shift_mm_vs_nn"], + "seconds_warp": c["seconds_warp"]} for c in combos}, + } + sf = os.path.join(case_dir, "summary.json") + with open(sf, "w") as f: + json.dump(summary, f, indent=1) + print(f"summary: {sf}", flush=True) + print("done", flush=True) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index 3f394e9..c4ad52e 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -359,12 +359,15 @@ def _itk_transform_to_sitk_rigid(itk_tf): return t -def register_elastix(ct, t1c, workdir=None, threads=8): - """ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI, full images. +def register_elastix(ct, t1c, workdir=None, threads=8, + transform="EulerTransform"): + """ITK-Elastix rigid/affine CT->T1c (fixed=T1c, moving=CT), Mattes MI, + full images. Runs in-process through the ITK-Elastix python API (pip ``itk`` + - ``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the - 5.x replacement for ``Rigid3DTransform``) + Mattes MI + ``itk-elastix``, elastix 5.x). ``transform`` = ``EulerTransform`` (rigid + 6-DOF, the 5.x replacement for ``Rigid3DTransform``) or + ``AffineTransform`` (12-DOF, the affine benchmark) + Mattes MI (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the full CT / T1c volumes (no head box). No manual centroid init: elastix's ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity @@ -386,7 +389,7 @@ def register_elastix(ct, t1c, workdir=None, threads=8): "ImageDimension": 3, "FixedInternalImagePixelType": "float", "MovingInternalImagePixelType": "float", - "Transform": "EulerTransform", + "Transform": transform, "TransformBase": "ElastixTransformBase", "FixedCenteredTransform": True, "TransformParameterReadoutFormat": 1, @@ -402,7 +405,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8): "ImageSamplerPatchSize": [25, 25, 25], "WriteResultImage": False, } - outdir = os.path.join(workdir or d("data/m6reg/elastix"), "itkelastix") + outdir = os.path.join(workdir or d("data/m6reg/elastix"), + "itkelastix" if transform == "EulerTransform" + else "itkelastix_aff") os.makedirs(outdir, exist_ok=True) import time, io, contextlib t0 = time.time() @@ -435,7 +440,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8): meta = { "pipeline": "elastix", "toolkit": version, "metric": "MattesMutualInformation (full volumes)", - "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)", + "transform": (f"{'Rigid' if transform == 'EulerTransform' else 'Affine'} " + f"({transform}, ASGD, 3 resolutions, auto center init)"), "convention": "T1c->CT (verified)", "auto_inverted": inverted, "fixed": _modality_name(fixed), "moving": _modality_name(moving), "seconds_reg": round(dt, 1), @@ -443,6 +449,13 @@ def register_elastix(ct, t1c, workdir=None, threads=8): return T, meta +def register_elastix_affine(ct, t1c, workdir=None, threads=8): + """Affine (12-DOF) variant of :func:`register_elastix` — same budget, + ``AffineTransform`` instead of ``EulerTransform``.""" + return register_elastix(ct, t1c, workdir=workdir, threads=threads, + transform="AffineTransform") + + # --------------------------------------------------------------------------- # # 3. ANTsPy # --------------------------------------------------------------------------- # @@ -738,14 +751,40 @@ def _impact_models_host(): return list(dict.fromkeys(rels)) -def _impact_param_map(gpu, iters, samples, voxel_sizes): - """Build the Elastix-5.2 parameter map for a RIGID (Euler) IMPACT - registration: the MIND + TotalSegmentator static feature combination of - the ImpactReg ``MR_CT_TS`` reference map (one model pair per level, - coarse -> fine, see :data:`IMPACT_LEVELS`) + a secondary Mattes-MI metric - (:data:`IMPACT_MI_WEIGHT`). The image pyramid is NOT resampled (all-1 - schedules); the scale comes from the per-level ImpactVoxelSize (as the - IMPACT docs recommend). ``iters`` is a per-level tuple or a single int. +def _itk_affine_txt_to_sitk(path): + """Parse an ITK ``AffineTransform`` .txt file into a SimpleITK + AffineTransform (fallback for when SimpleITK's own reader does not accept + the file). ``Parameters`` = 9 matrix values (row-major) + 3 translation, + plus a ``Center`` line; the mapping is p -> M (p - c) + c + t.""" + text = open(path).read() + if "AffineTransform" not in text: + raise RuntimeError(f"{path} is not an AffineTransform file") + body = text.split("Parameters:", 1)[1].split("FixedParameters:", 1)[0] + params = [float(x) for x in body.split()] + if len(params) != 12: + raise RuntimeError(f"expected 12 affine parameters in {path}, " + f"got {len(params)}") + m = re.search(r"^Center:\s*([\d\.\-e+ ]+)", text, re.M) + center = ([float(x) for x in m.group(1).split()[:3]] if m + else [0.0, 0.0, 0.0]) + t = sitk.AffineTransform(3) + t.SetCenter(tuple(center)) + t.SetMatrix(params[0:9]) + t.SetTranslation(params[9:12]) + return t + + +def _impact_param_map(gpu, iters, samples, voxel_sizes, + transform="EulerTransform"): + """Build the Elastix-5.2 parameter map for a RIGID (Euler) or AFFINE + IMPACT registration: the MIND + TotalSegmentator static feature + combination of the ImpactReg ``MR_CT_TS`` reference map (one model pair + per level, coarse -> fine, see :data:`IMPACT_LEVELS`) + a secondary + Mattes-MI metric (:data:`IMPACT_MI_WEIGHT`). The image pyramid is NOT + resampled (all-1 schedules); the scale comes from the per-level + ImpactVoxelSize (as the IMPACT docs recommend). ``iters`` is a per-level + tuple or a single int; ``transform`` = ``EulerTransform`` (6-DOF) or + ``AffineTransform`` (12-DOF, the affine benchmark). """ n = len(voxel_sizes) if isinstance(iters, (int,)): @@ -790,8 +829,9 @@ def _impact_param_map(gpu, iters, samples, voxel_sizes): '(CheckNumberOfSamples "true")', "(RequiredRatioOfValidSamples 0.1)", "", - "// rigid 6-DOF to be optimized (auto center init, no -t0)", - '(Transform "EulerTransform")', + f"// {'rigid 6-DOF' if transform == 'EulerTransform' else 'affine 12-DOF'} " + "to be optimized (auto center init, no -t0)", + f'(Transform "{transform}")', '(TransformBase "ElastixTransformBase")', '(TransformParameterReadoutFormat "AllParametersOneLine")', '(FixedCenteredTransform "true")', @@ -849,8 +889,8 @@ def _impact_param_map(gpu, iters, samples, voxel_sizes): def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, samples=IMPACT_DEFAULT_SAMPLES, threads=12, - keep_scratch=False): - """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image + keep_scratch=False, transform="EulerTransform"): + """IMPACT rigid/affine CT->T1c (fixed=T1c, moving=CT) in the Docker image ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric). IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic* @@ -930,7 +970,7 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, os.path.join(datadir, "Fixed_image.mha"), True) sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), os.path.join(datadir, "Moving_image.mha"), True) - param_text = _impact_param_map(gpu, iters, samples, voxel_sizes) + param_text = _impact_param_map(gpu, iters, samples, voxel_sizes, transform) with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: f.write(param_text) @@ -968,11 +1008,14 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}") - # the optimized Euler (world frame) is the full combination transform + # the optimized transform (world frame) is the full combination # (AutomaticTransformInitialization is folded in by elastix); with - # fixed=T1c / moving=CT the direction it holds is resolved empirically - # by fix_orientation against the skull ring. - T_main = sitk.ReadTransform(itk_tf) + # fixed=T1c / moving=CT the direction it holds is resolved + # empirically by fix_orientation against the skull ring. + try: + T_main = sitk.ReadTransform(itk_tf) + except Exception: # noqa: SimpleITK reader lacks this ITK format + T_main = _itk_affine_txt_to_sitk(itk_tf) T, inverted = fix_orientation(ct, t1c, T_main) finally: if not keep_scratch: @@ -984,7 +1027,8 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", "metric": (f"IMPACT semantic ({models}, L1+Dice, Static) + " f"Mattes MI (w {IMPACT_MI_WEIGHT})"), - "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; " + "transform": (f"{'Rigid' if transform == 'EulerTransform' else 'Affine'} " + f"({transform}, ASGD, {n} scales, auto center init; " f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), "convention": "T1c->CT (verified)", "auto_inverted": inverted, "fixed": _modality_name(fixed), "moving": _modality_name(moving), @@ -998,6 +1042,17 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, return T, meta +def register_impact_affine(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, + samples=IMPACT_DEFAULT_SAMPLES, threads=12, + keep_scratch=False): + """Affine (12-DOF) variant of :func:`register_impact` — the MIND+TS + reference map with ``AffineTransform`` instead of ``EulerTransform``, + same budgets and M6 adaptations.""" + return register_impact(ct, t1c, gpu=gpu, iters=iters, samples=samples, + threads=threads, keep_scratch=keep_scratch, + transform="AffineTransform") + + # --------------------------------------------------------------------------- # # 6. NiftyReg (KCL-BMEIS suite, reg_aladin) # --------------------------------------------------------------------------- # @@ -1064,7 +1119,8 @@ def _parse_niftyreg_affine(path): def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False, - backend="auto", gpu=None, clip_ct=(0.0, 80.0)): + backend="auto", gpu=None, clip_ct=(0.0, 80.0), + rigid_only=True): """NiftyReg rigid CT->T1c (reference=T1c, floating=CT) via ``reg_aladin``. reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global @@ -1147,10 +1203,12 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals os.path.join(tmp, "ref.nii.gz"), True) sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), os.path.join(tmp, "flo.nii.gz"), True) - flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz", - "-rigOnly", "-aff", "aff.txt", - "-ln", str(levels), "-maxit", str(maxit), - "-interp", interp, "-omp", str(threads)] + flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz"] + if rigid_only: + flags.append("-rigOnly") + flags += ["-aff", "aff.txt", + "-ln", str(levels), "-maxit", str(maxit), + "-interp", interp, "-omp", str(threads)] if platf is not None: flags += ["-platf", platf] if use_gpu: @@ -1194,8 +1252,9 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals "toolkit": toolkit, "backend": "docker" if use_docker else "conda", "metric": "NCC block matching + trimmed least squares (symmetric, full volumes)", - "transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, " - f"{maxit} TLS iters/level, cubic interp)"), + "transform": (f"{'Rigid' if rigid_only else 'Affine'} (reg_aladin" + f"{' -rigOnly' if rigid_only else ''}, {levels}-level " + f"pyramid, {maxit} TLS iters/level)"), "convention": "T1c->CT (verified)", "auto_inverted": inverted, "fixed": fixed_name, "moving": moving_name, "seconds_reg": round(dt, 1), @@ -1207,6 +1266,17 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals return T, meta +def register_niftyreg_affine(ct, t1c, threads=24, levels=3, maxit=15, + keep_scratch=False, backend="auto", gpu=None, + clip_ct=(0.0, 80.0)): + """Affine (12-DOF) variant of :func:`register_niftyreg` — ``reg_aladin`` + full affine (no ``-rigOnly``), same pyramid / TLS budget.""" + return register_niftyreg(ct, t1c, threads=threads, levels=levels, + maxit=maxit, keep_scratch=keep_scratch, + backend=backend, gpu=gpu, clip_ct=clip_ct, + rigid_only=False) + + REGISTERS = { "antspyx": register_antspyx, "elastix": register_elastix, @@ -1218,6 +1288,203 @@ REGISTERS = { ALL_PIPELINES = list(REGISTERS) +# --------------------------------------------------------------------------- # +# Affine variants (see REGISTRATION_BENCHMARK_AFFINE.md) +# --------------------------------------------------------------------------- # +# +# Same unified contract as the rigid pipelines, but every toolkit estimates a +# full affine (12-DOF) transform. The antspyx rows run stock linear-only +# presets; the other toolkits keep their rigid-benchmark budgets with an +# affine transform in place of the rigid one. All return (T, meta) with T an +# affine (or composite of affines) mapping T1c -> CT, direction resolved by +# fix_orientation like the rigid pipelines. + +ANT_AFFINE_PRESETS = { + # pipeline name -> (antspyx type_of_transform, description) + "antspyx": ( + "Affine", + "direct Affine preset: Affine[0.25] (Mattes 32, 2100x1200x1200x10, " + "6x4x2x1, 3x2x1x0); COM init"), + "ants_trsaa": ( + "TRSAA", + "TRSA cascade: Translation[1] -> Rigid[1] -> Similarity[1] (each " + "Mattes 32, 2000x2000x0, 4x2x1, 2x1x0) -> Affine[1] x2 (Mattes 32, " + "40x20x0, 4x2x1, 2x1x0); COM init"), + "ants_reproa": ( + "antsRegistrationSyNRepro[a]", + "Repro (a): Rigid[0.1] + Affine[0.1] (each GC 1/1/Regular/0.25, " + "1000x500x250x100, 8x4x2x1, 3x2x1x0vox); COM init"), + "ants_quickreproa": ( + "antsRegistrationSyNQuickRepro[a]", + "QuickRepro (a): Rigid[0.1] + Affine[0.1] (each GC 1/1/Regular/0.25, " + "1000x500x250x0, 8x4x2x1, 3x2x1x0vox); COM init"), +} + + +def register_ants_affine_variant(ct, t1c, preset_name): + """Run one stock antspyx linear-only preset (:data:`ANT_AFFINE_PRESETS`) + as an affine pipeline under the unified contract (fixed = T1c, moving = + CT, full volumes, no clip). All stages are linear, so the stage files + ANTs writes (in numeric order = the order it applies them) are composed + into the contract affine transform; binary merged-linear .mat files + (Ants combines consecutive linear stages) are recovered by point-probing + (:func:`_ants_stage_to_sitk`, exact for affine maps). OMP threads are + left at the system default, as for the other antspyx rows.""" + if preset_name not in ANT_AFFINE_PRESETS: + raise KeyError(f"unknown antspyx affine preset {preset_name!r} " + f"(one of {list(ANT_AFFINE_PRESETS)})") + tx_type, desc = ANT_AFFINE_PRESETS[preset_name] + try: + import ants + except Exception as e: # noqa + raise PipelineUnavailable(f"antspyx: {e!r}") + import glob + import time + fx, mv = fixed_moving(ct, t1c) + af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)) + am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) + tmp = tempfile.mkdtemp(prefix=f"ants_aff_{preset_name}_") + dt = 0.0 + stage_names = [] + try: + t0 = time.time() + r = ants.registration(fixed=af, moving=am, type_of_transform=tx_type, + outprefix=tmp + "/") + dt = time.time() - t0 + stage_files = [f for f in glob.glob(os.path.join(tmp, "*")) + if f.endswith(".mat")] + if not stage_files: + raise RuntimeError(f"{preset_name}: no .mat stage files in {tmp}") + stage_files.sort(key=lambda p: int(re.match( + r"(\d+)", os.path.basename(p)).group(1))) + stage_names = [os.path.basename(p) for p in stage_files] + stages = [_ants_stage_to_sitk(p, fx) for p in stage_files] + T = stages[0] if len(stages) == 1 else sitk.CompositeTransform(stages) + T, inverted = fix_orientation(ct, t1c, T) + finally: + shutil.rmtree(tmp, ignore_errors=True) + meta = { + "pipeline": preset_name, + "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", + "metric": "ANTs stock preset metrics (Mattes MI / GC per stage)", + "transform": f"Affine (antspyx preset '{tx_type}'): {desc}", + "convention": "T1c->CT (verified)", "auto_inverted": bool(inverted), + "fixed": _modality_name(fx), "moving": _modality_name(mv), + "seconds_reg": round(dt, 1), + "stage_files": stage_names, + } + return T, meta + + +def register_simpleitk_affine(ct, t1c, threads=12, iters_full=2000): + """SimpleITK affine CT->T1c, Mattes MI, full volumes: the rigid + benchmark's Euler3D stage (``iters_full`` iters/level on the [4,2,1] + pyramid, identity start) followed by an ``AffineTransform`` stage with + the same budget, initialised to reproduce the rigid stage exactly so the + affine refines the pose instead of re-estimating it (a single-stage affine + from the identity collapses its scale on these pairs). The returned + transform is the composite condensed to a single affine by point mapping + (T1c -> CT).""" + m6 = _load_m6() + _sitk_set_threads(threads) + ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32) + fixed, moving = fixed_moving(ct, t1c) + shape = np.array(ct.GetSize()) + sh = [s for s in [4, 2, 1] if min(shape) / s >= 6] or [1] + sm = [2, 1, 0][:len(sh)] + center = tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint( + ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2, + (fixed.GetSize()[2] - 1) // 2))) + + def _method(): + rm = sitk.ImageRegistrationMethod() + rm.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32) + rm.SetOptimizerAsRegularStepGradientDescent( + learningRate=10.0, minStep=0.2, numberOfIterations=iters_full, + relaxationFactor=0.5, gradientMagnitudeTolerance=1e-6) + rm.SetOptimizerScalesFromPhysicalShift() + rm.SetShrinkFactorsPerLevel(sh); rm.SetSmoothingSigmasPerLevel(sm) + rm.SetInterpolator(sitk.sitkLinear) + return rm + + def _condense_affine(tr): + """Condense an affine (or composite of affines) to a single + AffineTransform by point mapping (exact for affine maps).""" + o = np.array(tr.TransformPoint((0.0, 0.0, 0.0)), dtype=float) + M = np.empty((3, 3), dtype=float) + for i in range(3): + p = [0.0, 0.0, 0.0] + p[i] = 1.0 + M[:, i] = np.array(tr.TransformPoint(tuple(p)), dtype=float) - o + a = sitk.AffineTransform(3) + a.SetMatrix(M.ravel()) + a.SetTranslation(o) + return a + + import time + t0 = time.time() + init1 = sitk.Euler3DTransform() + init1.SetCenter(center) + init1.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) + r1 = _method() + r1.SetInitialTransform(init1, inPlace=False) + tr1 = _condense_affine(r1.Execute(fixed, moving)) + + init2 = sitk.AffineTransform(3) + init2.SetCenter(tr1.GetCenter()) + init2.SetMatrix(tr1.GetMatrix()) + init2.SetTranslation(tr1.GetTranslation()) + r2 = _method() + r2.SetInitialTransform(init2, inPlace=False) + tr = _condense_affine(r2.Execute(fixed, moving)) + dt = time.time() - t0 + ha = sitk.GetArrayFromImage(t1c) + gap = m6._skull_ring_gap(_warp_ct(ct, t1c, tr), ha) + meta = { + "pipeline": "simpleitk", "mode": "benchmark", + "toolkit": f"SimpleITK {sitk.Version_VersionString()}", + "metric": "MattesMutualInformation (both stages, full volumes)", + "transform": (f"Affine (Euler3D rigid init + AffineTransform, " + f"{iters_full} iters/level on pyramid {sh} each, " + f"identity start on the rigid, affine initialised to " + f"reproduce the rigid stage)"), + "convention": "T1c->CT (benchmark budget)", "auto_inverted": False, + "fixed": _modality_name(fixed), "moving": _modality_name(moving), + "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm, + "skull_ring_gap_mm": gap, + } + return tr, meta + + +def _ants_affine_register(name): + def _reg(ct, t1c): + return register_ants_affine_variant(ct, t1c, preset_name=name) + return _reg + + +def register_synthmorph_affine(ct, t1c, gpu=None, threads=None, extent=256, + keep_scratch=False, clip_ct=(0.0, 80.0)): + """Affine (12-DOF) variant of :func:`register_synthmorph` — the + ``-m affine`` matrix model.""" + return register_synthmorph(ct, t1c, model="affine", gpu=gpu, threads=threads, + extent=extent, keep_scratch=keep_scratch, + clip_ct=clip_ct) + + +REGISTERS_AFFINE = { + "antspyx": _ants_affine_register("antspyx"), + "ants_trsaa": _ants_affine_register("ants_trsaa"), + "ants_reproa": _ants_affine_register("ants_reproa"), + "ants_quickreproa": _ants_affine_register("ants_quickreproa"), + "elastix": register_elastix_affine, + "impact": register_impact_affine, + "niftyreg": register_niftyreg_affine, + "simpleitk": register_simpleitk_affine, + "synthmorph": register_synthmorph_affine, +} +ALL_PIPELINES_AFFINE = list(REGISTERS_AFFINE) + + # --------------------------------------------------------------------------- # # Deformable variants (see REGISTRATION_BENCHMARK_DEFORMABLE.md) # --------------------------------------------------------------------------- # @@ -2311,8 +2578,10 @@ REGISTERS_DEFORMABLE.update({n: _ants_variant_register(n) for n in ANT_VARIANTS} REGISTERS_DEFORMABLE["ants_rigid_synonly"] = register_ants_rigid_synonly_deformable ALL_PIPELINES_DEFORMABLE = list(REGISTERS_DEFORMABLE) # every pipeline name the harness may dispatch / merge on (rigid + variants) -ALL_PIPELINE_NAMES = (ALL_PIPELINES + list(ANT_VARIANTS) - + ["ants_rigid_synonly"]) +ALL_PIPELINE_NAMES = ( + ALL_PIPELINES + list(ANT_VARIANTS) + ["ants_rigid_synonly"] + + [n for n in ALL_PIPELINES_AFFINE if n not in ALL_PIPELINES] +) def availability(): @@ -2374,6 +2643,7 @@ def availability(): except Exception as e: # noqa out[name] = (False, repr(e)) # the ANTs variants all run through antspyx (same availability as antspyx) - for name in list(ANT_VARIANTS) + ["ants_rigid_synonly"]: + for name in (list(ANT_VARIANTS) + ["ants_rigid_synonly"] + + [n for n in ALL_PIPELINES_AFFINE if n not in ALL_PIPELINES]): out[name] = out.get("antspyx", (False, "antspyx not checked")) return out \ No newline at end of file