diff --git a/AGENTS.md b/AGENTS.md index 5f5b11a..d90fb4c 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -10,7 +10,7 @@ Activate with: source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 ``` -Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. +Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), plus Docker images `freesurfer/synthmorph` and `vboussot/elastix_impact`. ## Project diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index 4c93ec6..7f16871 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/REGISTRATION_BENCHMARK.md @@ -1,6 +1,6 @@ # Registration Pipeline Benchmark — CT → T1c (M6) -Benchmarks **five rigid, cross-modality registration pipelines** that align a +Benchmarks **six rigid, cross-modality registration pipelines** that align a head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be warped onto the T1c (MR) native grid. This is the registration step of `scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT @@ -27,30 +27,39 @@ register on the **full CT / T1c volumes** (no head box) and use each toolkit's **native automatic transform initialization** (no manual centroid offset): elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity), antspyx its default center-of-mass alignment, simpleitk starts from identity, -synthmorph needs none. All use **fixed = CT, moving = T1c** except `impact`, -which keeps **fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much -larger than the T1c crop (~220×220×150 mm), so with CT as the fixed image too -many CT samples map outside the T1c buffer and the IMPACT metric aborts -("too many samples map outside moving image buffer"). The comparison is on -each toolkit's native optimizer / metric / global-search capability. +niftyreg (reg_aladin) initializes from the image centres, synthmorph needs +none. All use **fixed = CT, moving = T1c** except `impact`, which keeps +**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than +the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT +samples map outside the T1c buffer and the IMPACT metric aborts ("too many +samples map outside moving image buffer"). The comparison is on each +toolkit's native optimizer / metric / global-search capability. -## The five pipelines +## The six pipelines Listed alphabetically. All register on the **full volumes** (no head box) and use each toolkit's **native automatic transform initialization** (no manual centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk -starts from identity, and synthmorph needs none (global deep method). All use -**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c, moving -= CT** (swapping breaks its TotalSegmentator metric — see below). +starts from identity, niftyreg initializes from the image centres, and +synthmorph needs none (global deep method). All use **fixed = CT, moving = +T1c** except `impact`, which keeps **fixed = T1c, moving = CT** (swapping +breaks its TotalSegmentator metric — see below). | # | pipeline | fixed / moving | init | metric | transform | status in this env | |---|----------|----------------|------|--------|-----------|--------------------| | 1 | `antspyx` | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | -| 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (10 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) | +| 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (32 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) | | 3 | `impact` | **T1c / CT** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) | -| 4 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | -| 5 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) | +| 4 | `niftyreg` | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (main env, conda-forge `reg_aladin` 1.5.69) | +| 5 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | +| 6 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) | + +Benchmark budget (best-results settings): `elastix`, `impact` and `simpleitk` +run with **MaximumNumberOfIterations = 2000** (per resolution / per pyramid +level), with the sampler and metric left at their defaults +(**NumberOfSpatialSamples 5000**, **NumberOfHistogramBins 32**). `antspyx`, +`niftyreg` and `synthmorph` keep their native defaults. `simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid, centered start — the default, comparable to the other tools) and @@ -74,128 +83,123 @@ m6_2BBN5XA4_20250411 m6_2C4TKI6H_20240506 m6_2CDNN7BH_20240102 ## Metrics (identical for every pipeline) -Production QA, crop-independent: +The headline alignment score is **`mi_ct_mr`** — the registration's *achieved* +Mattes MI, read off as a standalone score: -- **`skull_ring_gap_mm`** — march rays from the MR head center; the gap between - the MR head silhouette boundary and the warped-CT skull ring (p97 of positive - CT) at mid and low levels, plus the angular std. **Aligned ⇔ - `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`.** `mid` measures translational - concentricity; `std` measures rotational concentricity. -- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT - head tissue (the slab must straddle the lesion). Acceptance needs `≥ 0.5`. -- **`tv_nvox`** — warped GTV volume. -- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the - production acceptance rule). - **`mi_ct_mr`** — mean Mattes MI (nats) between the **CT inverse-transformed into MR coordinates** (the warped CT on the T1c grid) and the **original MR**, restricted to the MR head and each image percentile-binned independently - (SimpleITK Mattes-MI style). This is the registration's *achieved* metric read - off as a standalone alignment score: higher = the warped-CT and MR intensities - are better in-sync = better aligned; it collapses toward 0 when the CT is - rotated/translated away from the MR. + (SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are + better in-sync = better aligned; it collapses toward 0 when the CT is + rotated/translated away from the MR. Tables and plots report this (plus + registration time and GTV-in-CT coverage). +- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT + head tissue (the slab must straddle the lesion). - **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring). +The production geometric QA is still computed per case and kept in the JSON +(it drives dataset acceptance in `02_build_m6_dataset.py`, but its summary +numbers were found to correlate poorly with perceived alignment, so they are no +longer reported in the tables/plots): + +- **`skull_ring_gap_mm`** — march rays from the MR head center; the gap between + the MR head silhouette boundary and the warped-CT skull ring (p97 of positive + CT) at mid and low levels, plus the angular std. `mid` measures + translational concentricity; `std` measures rotational concentricity. +- **`tv_nvox`** — warped GTV volume. +- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the + production acceptance rule; aligned ⇔ `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`). + ## Results (10 cases) -| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) | -|---|---|---|---|---|---|---|---|---| -| **antspyx** | ANTsPy 0.6.3 | 9 | **4/10 (0.40)** | 12.8 / 15.6 | 15.6 / 16.5 | 36.7 | 0.93 | **0.652 / 0.644** | -| **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 13.6 / 15.7 | 14.8 / 17.0 | 2.6 | 0.93 | **0.641 / 0.640** | -| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **4/10 (0.40)** | 11.6 / 11.2 | 16.2 / 17.3 | **833.5** | 0.93 | **0.511 / 0.493** | -| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 14.4 / 16.6 | 17.9 / 17.5 | 9.7 | 0.93 | **0.651 / 0.641** | -| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 21.8 / 21.5 | 25.0 | 0.93 | **0.238 / 0.308** | +Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over +the MR head; higher = better). -*(antspyx and impact are stochastic (random metric/patch subsampling), so their -accepted counts are one draw each — see Interpretation. impact's s/reg is ~14 -min/case on GPU (full volumes), orders of magnitude above the intensity-metric -tools. elastix, simpleitk and synthmorph are deterministic at fixed threads.)* - -Per-case (mid / std, `OK` = accepted): - -| case | antspyx | elastix | impact | simpleitk | synthmorph | +| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) | |---|---|---|---|---|---| -| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | -| m6_22IRRKU6_20240115 | 24.6 / 28.6 | 24.0 / 28.3 | 22.7 / 27.0 | 26.3 / 29.6 | 10.4 / 23.7 | -| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 22.3 / 23.5 | 9.1 / 16.2 | 21.9 / 23.2 | 20.3 / 23.3 | -| m6_24V62QQM_20230223 | 9.5 / 15.6 | 11.6 / 14.8 | 8.6 / 7.2 **OK** | 14.4 / 23.2 | 11.1 / 21.8 | -| m6_24V62QQM_20231017 | 11.7 / 7.3 **OK** | 12.2 / 7.4 **OK** | 4.6 / 6.6 **OK** | 11.7 / 7.5 **OK** | 11.5 / 8.4 **OK** | -| m6_24V62QQM_20240820 | 10.8 / 9.1 **OK** | 11.4 / 9.0 **OK** | 11.6 / 4.9 **OK** | 11.5 / 9.1 **OK** | 43.4 / 15.8 | -| m6_26TYQI3C_20221027 | 12.3 / 4.4 **OK** | 10.6 / 3.1 **OK** | 11.9 / 4.5 **OK** | 11.3 / 4.6 **OK** | 11.8 / 2.9 **OK** | -| m6_2BBN5XA4_20250411 | 19.1 / 7.7 **OK** | 18.2 / 14.1 | −1.7 / 29.3 | 18.8 / 7.5 **OK** | 52.4 / 53.0 | -| m6_2C4TKI6H_20240506 | 12.8 / 17.2 | 13.6 / 20.0 | 15.4 / 23.2 | 14.0 / 17.9 | 16.7 / 28.2 | -| m6_2CDNN7BH_20240102 | 17.5 / 34.6 | 17.6 / 32.8 | 18.9 / 37.1 | 19.3 / 35.0 | 18.9 / 16.8 | +| **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** | +| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** | +| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** | +| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** | +| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** | +| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** | + +\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default +(MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower +(hours/case) and was not re-run. + +Per-case MI (nats; **bold** = best for the case): + +| case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph | +|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | 0.649 | 0.472 | **0.624** | 0.296 | 0.642 | 0.235 | +| m6_23XIMO4M_20230222 | 0.655 | 0.416 | 0.474 | **0.634** | 0.652 | 0.510 | +| m6_24V62QQM_20230223 | 0.703 | **0.701** | 0.554 | 0.260 | 0.663 | 0.238 | +| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.419 | 0.493 | **0.602** | 0.516 | +| m6_24V62QQM_20240820 | 0.652 | 0.643 | 0.458 | 0.118 | **0.652** | 0.096 | +| m6_26TYQI3C_20221027 | 0.657 | 0.653 | 0.511 | 0.608 | **0.695** | 0.604 | +| m6_2BBN5XA4_20250411 | 0.576 | 0.287 | 0.155 | 0.069 | **0.577** | 0.060 | +| m6_2C4TKI6H_20240506 | 0.645 | 0.401 | **0.626** | 0.231 | 0.634 | 0.209 | +| m6_2CDNN7BH_20240102 | 0.655 | 0.516 | **0.618** | 0.344 | 0.651 | 0.301 | ### Interpretation -- **Translation is recovered well by all five** (mid gaps mostly 4–16 mm, - within the `≤ 20 mm` band, except the grossly-drifted cases). - `tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for - every pipeline, so the target is in the right place. -- **The binding failure is rotation, and removing the head box makes it - harder.** On the repositioned cases `std > 12 mm` — the ~30° inter-session - tilt produces local minima a single auto-start locks onto, and without the - head-box crop the full-volume metric sees more neck/air confounders. Notably, - the hardest tilt case `m6_22IRRKU6` is now **rejected by all five** (std - 23–30 mm), which the head-box config got through with IMPACT. This is exactly - why production `02_build_m6_dataset.py` uses a **multi-start tilt grid**. -- **The easy cases are now accepted by nearly all pipelines.** - `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820` and `m6_26TYQI3C` are accepted - by ANTs/Elastix/IMPACT/SimpleITK (SynthMorph accepts the two cleanest). The - native automatic inits (center-of-gravity / center-of-mass) place the heads - close enough that the optimizers converge on the well-positioned cases. -- **SimpleITK, ANTs and IMPACT tie at 4/10; Elastix 3/10; SynthMorph 2/10.** - ANTs is no longer the weakest — on the full volumes its native center-of-mass - init works well (median std 15.6, up from the head-box+manual-pre-shift - 31.7). It is stochastic (random MI subsampling), so the count varies by draw. -- **IMPACT (semantic metric) is 4/10 but the slowest (~14 min/case on full - volumes) and lost its standout on the hardest tilt** (`m6_22IRRKU6` now - mid 22.7 / std 27.0 vs mid 8.3 / std 13.2 with the head box). Matching - *anatomical semantics* still gives it the best median mid (11.6) and a - competitive MI (0.511), but the full-volume registration is a larger problem - for the local ASGD. It also drifts grossly on `m6_2BBN5XA4` (mid −1.7, the - transform flips to the far side). -- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.6 s - (warm workers), SimpleITK ~9.7 s, ANTs ~36.7 s (CPU-bound, full volumes), - **SynthMorph ~25 s** (GPU: a single pretrained-model forward pass per case), - **IMPACT ~834 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator - extractor on the sampled patches with backprop). - - **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**: - medians 0.651 (SimpleITK) / 0.641 (Elastix) / 0.511 (IMPACT) / 0.652 (ANTs) / - 0.238 (SynthMorph). - The intensity-metric tools cluster at ~0.64 (a genuine alignment signal), - while SynthMorph's 0.238 reflects its tilt failures. A high MI means the CT - skull/tissue intensities fall on the matching MR structures after the inverse - transform. -- **SynthMorph (deep, pretrained) is near-weakest of the five on these - cross-modality CT pairs.** Its rigid model is trained on brain *MRI* - (T1/T2, anatomy-aware but intensity-domain MRI), so it recovers gross - translation well (mid 10–20 mm on most cases) but is repeatedly defeated by - the ~30° inter-session tilt: 7 of 9 cases have `std > 15 mm` (rotation error) - and two drift grossly (`m6_2BBN5XA4` mid 52 / std 53, `m6_24V62QQM_20240820` - mid 43). It accepts the clean `m6_26TYQI3C` (mid 11.8 / std 2.9, MI 0.60) and - `m6_24V62QQM_20231017` (mid 11.5 / std 8.4) — where it is competitive with the - intensity tools. Notably it needs **no centroid init** (it is a global deep - method that recovers the gross offset on its own). Its low median MI (0.238) - reflects the tilt failures. On same-modality MRI it would likely rank much - higher; this CT benchmark is an out-of-distribution stress test for it. It is - deterministic (one fixed pretrained model, no random subsampling), so unlike - ANTs/IMPACT the per-case numbers are stable run to run. - - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its +- **SimpleITK and ANTsPy tie for best** (median MI 0.651 / 0.652): on the + well-positioned cases they cluster at 0.60–0.70 and both survive the harder + repositioned cases (0.58–0.65 on `m6_2BBN5XA4` / `m6_22IRRKU6`). SimpleITK's + per-case MI is **unchanged** from the previous 150-iteration budget — RSGD + converges well before the 2000-iteration cap, so the cap is now headroom + rather than a constraint. +- **The grossly repositioned cases define the ranking.** On + `m6_2BBN5XA4` and `m6_24V62QQM_20240820` (~30° inter-session tilt) the tools + split: ANTs/SimpleITK hold 0.58–0.65, elastix drops to 0.29–0.64, niftyreg + 0.07–0.12, synthmorph 0.06–0.10 — single-start optimizers locking onto local + minima. This is exactly why production `02_build_m6_dataset.py` uses a + **multi-start tilt grid**. +- **Elastix is second at the new budget** (2000 iters per resolution, 32 bins, + 5000 samples): median 0.516, below the 0.641 recorded at the old 500/500/250 + + 10-bin setting. The long ASGD run nails the clean cases (0.64–0.70) but its + single CenterOfGravity start still locks onto local minima on the + repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M` + 0.416). +- **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI + on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618), + where anatomical-segmentation similarity out-resolves intensity MI — but its + median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is + also by far the slowest (~14 min/case at the shown budget; the new 2000/5000 + budget would be on the order of hours/case). +- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching + is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class; + `m6_26TYQI3C` 0.608) but falls into the same tilt local minima + (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower + than elastix at full volumes. +- **SynthMorph is weakest of the six on these cross-modality CT pairs** + (median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30° + tilt defeats it on the repositioned cases (0.06–0.10) even though it matches + the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no + centroid init** (global deep method) and is deterministic. +- **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the + CT head everywhere, so coverage is not the discriminator — MI is. +- **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s, + SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s, + **NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous + budget; GPU feature extraction with backprop per iteration). +- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its per-case numbers carry run-to-run variance (this table is one draw); - SimpleITK's MI (and registration) is deterministic at fixed thread count - (RSGD wall time is thread-sensitive, convergence is not); the in-process - Elastix ASGD is deterministic given its fixed seed. **IMPACT is the most - stochastic**: its Jacobian mode draws random patches every iteration, so the - final pose — and even *which* cases get accepted — changes run to run (two - independent 10-case draws both accepted 3/10 but on different cases). Treat - its per-case numbers as one sample from a distribution. + SimpleITK's MI (and registration) is deterministic at fixed thread count; + the in-process Elastix ASGD is deterministic given its fixed seed; + SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**: + its Jacobian mode draws random patches every iteration, so the final pose + changes run to run. Treat its per-case numbers as one sample from a + distribution. ## Reproducing ```bash -# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here) +# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (pip/conda-installed here) source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 -python scripts/bench_registration.py --pipelines simpleitk elastix antspyx \ +python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg \ --n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark # Docker pipelines (need the pulled image; see the SynthMorph / IMPACT @@ -211,10 +215,12 @@ python scripts/bench_registration.py --merge Other flags: `--sample random --seed 0` (reproducible random 10), `--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`, -`--threads T`, `--simpleitk-mode {benchmark,production}`, +`--threads T` (also the OpenMP thread count for `reg_aladin`), +`--simpleitk-mode {benchmark,production}`, `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu `, `--synthmorph-extent `, `--impact-model `, -`--impact-gpu `, `--impact-iters N`, `--impact-samples N`, +`--impact-gpu `, `--impact-iters N` (default 2000), +`--impact-samples N` (default 5000), `--out-dir `, `--qa-images` / `--qa-out ` (visual QA, below). ## Visual QA images (checkerboard + color overlay) @@ -236,9 +242,9 @@ parenchyma is visible and comparable to the MR (the skull, ~800+ HU, clips to white); the MR uses a 2–98 % head window. The coronal and sagittal views are oriented superior-up to match the axial view. -The title carries the pipeline, case key, Mattes MI, and the skull-ring gap -(mid/std), so a glance tells you both the achieved-metric score and whether the -visuals agree. Generate for the available pipelines with: +The title carries the pipeline, case key and the achieved Mattes MI, so a +glance tells you whether the visuals agree with the score. Generate for the +available pipelines with: ```bash source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 @@ -290,6 +296,41 @@ No extra model assets: the rigid weights ship inside the image removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so concurrent benchmark workers can share a GPU. +## NiftyReg pipeline — one-time setup + +The `niftyreg` pipeline runs `reg_aladin` from the [NiftyReg](https://github.com/KCL-BMEIS/niftyreg) +suite (KCL-BMEIS, the actively-maintained NiftyReg — the release zips for +v2.x need a newer glibc, so this env uses the conda-forge build): +**symmetric block matching + trimmed least squares** (Modat et al., +J. Med. Img. 1(2) 024003, 2014). Instead of a histogram MI objective it +matches 4-voxel image blocks by **normalised cross-correlation** at each level +of a 3-level Gaussian pyramid and solves a robust (trimmed) least-squares for +the rigid transform — a genuinely different algorithmic family, and +cross-modality tolerant because NCC is invariant to intensity scaling. + +Prerequisites (check with `RP.availability()`; a missing piece raises +`PipelineUnavailable` with the exact command): + +```bash +# one-time, in the benchmark env (provides reg_aladin + reg_f3d, NiftyReg 1.5.69) +conda install -c conda-forge niftyreg +``` + +### How the adapter drives the tool + +1. **Rigid 6-DOF, full volumes** — `reg_aladin -ref CT -flo T1c -rigOnly` + (symmetric variant on by default, which matters: the non-symmetric mode + drifts grossly on these repositioned pairs). `-ln 3` pyramid, + `-maxit 15` TLS iterations per level (the tool's default of 5 is too + coarse), `-interp 3` cubic internal warping, default initialisation + (image centres). +2. **Read-back** — the output affine (`-aff`) is a 4×4 row-major matrix, + documented as `Affine * Reference = Floating` (CT→T1c); the adapter parses + it into a SimpleITK `AffineTransform` and runs it through the usual + `fix_orientation` gate (skull-ring concentricity). On the M6 cases the gate + keeps the file's affine as-is (`auto_inverted: False`), i.e. it is already + the contract transform (T1c→CT). + ## IMPACT pipeline (Docker) — one-time setup The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact) @@ -359,15 +400,17 @@ binary directly (`/lib/elastix-install/bin/elastix`): head-to-head but is the config the dataset actually ships with. - `synthmorph` is deterministic (one fixed pretrained model, no random subsampling), so its per-case numbers are stable run to run — unlike ANTs / - IMPACT. Its low acceptance (2/10) is a genuine cross-modality / tilt - weakness, not variance. -- `fix_orientation` (elastix/ants/synthmorph/impact) chooses between a transform - and its inverse by skull-ring concentricity — a correctness gate that makes - the result independent of each toolkit's transform-direction convention. With - the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the - inverse of the contract, so it reports `auto_inverted: True` (the gate - corrects it); antspyx/synthmorph/impact return the contract direction - directly (`auto_inverted: False`). + IMPACT. Its low MI (median 0.238; 0.06–0.10 on the repositioned cases) is a + genuine cross-modality / tilt weakness, not variance. +- `fix_orientation` (elastix/ants/niftyreg/synthmorph/impact) chooses between a + transform and its inverse by skull-ring concentricity — a correctness gate + that makes the result independent of each toolkit's transform-direction + convention. With the fixed=CT/moving=T1c swap, `elastix`'s combination + transform is now the inverse of the contract, so it reports + `auto_inverted: True` (the gate corrects it); antspyx/synthmorph/impact + return the contract direction directly (`auto_inverted: False`). `niftyreg`'s + `-aff` file is documented as CT→T1c but the gate keeps it as-is on these + cases (`auto_inverted: False`) — the empirical gate, not the docs, decides. - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU): - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly @@ -382,6 +425,25 @@ binary directly (`/lib/elastix-install/bin/elastix`): `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU. - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK (~10 s) / ANTs (~37 s) and IMPACT (~834 s). +- **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`, + symmetric block matching + TLS, CPU/OpenMP): + - **Different algorithmic family.** No histogram MI objective: 4-voxel + blocks are matched by normalised cross-correlation and the rigid + transform is solved by a trimmed least squares, iterated per pyramid + level (`-ln 3 -maxit 15`). NCC's scale invariance is what makes it + cross-modality tolerant. + - **Symmetric mode matters.** The default symmetric variant (optimising + forward and backward transforms concurrently) is the robust one on these + repositioned pairs; `-noSym` drifts grossly (gap std 45 mm on a clean + case every other tool passes). + - **`-maxit 5` (the tool default) is too coarse** — on the clean case + `m6_26TYQI3C` it reaches MI 0.41 vs 0.61 at 15 TLS iterations/level; + 15 is the benchmark setting. + - **CPU-only, thread-bound** (`-omp`); ~4 min/case median (120–440 s range) + at full volumes with 24 threads. + - **Version note.** The upstream KCL-BMEIS release zips (v2.x) need + glibc ≥ 2.33; this host (Ubuntu 20.04, glibc 2.31) uses the conda-forge + build (1.5.69) instead. - **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python API — no CLI binary involved): - Elastix 5.x renamed components vs the classic parameter files: the rigid @@ -403,15 +465,19 @@ binary directly (`/lib/elastix-install/bin/elastix`): (`GetCombinationTransform().TransformPoint`) into an affine — independent of elastix's Euler-angle parameterization — and wrapped in a SimpleITK `AffineTransform(3)`. - - Registration wall time is ~2.4 s/case median at full budget (ASGD - 500/500/250 over 3 resolutions). + - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution + (ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their + defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32). - **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0, GPU, Jacobian mode): - - **Slow, by design.** Every ASGD iteration runs the TotalSegmentator - extractor on the sampled patches with backpropagation, so the recommended - 500-iter / 2000-sample budget is ~14 min/case on full volumes (vs 2–37 s - for the other tools). `--impact-iters` / `--impact-samples` trade accuracy - for speed. +- **Slow, by design.** Every ASGD iteration runs the TotalSegmentator + extractor on the sampled patches with backpropagation. The benchmark + budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 = + the elastix default) is ~10× the authors' 500-iter / 2000-sample + Jacobian budget, i.e. on the order of *hours* per case on full volumes — + the results shown for `impact` are from the previous 500/2000 budget and + were not re-run at the new budget. `--impact-iters` / `--impact-samples` + trade accuracy for speed. - **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index 3b2313b..aec7571 100644 --- a/scripts/bench_registration.py +++ b/scripts/bench_registration.py @@ -1,37 +1,41 @@ -"""Benchmark the 5 rigid CT->T1c registration pipelines on N M6 test cases. +"""Benchmark the 6 rigid CT->T1c registration pipelines on N M6 test cases. Pipelines (see src/registration_pipelines.py for the shared transform contract): - simpleitk | elastix | antspyx | synthmorph | impact + simpleitk | elastix | antspyx | synthmorph | impact | niftyreg Each pipeline registers a head-and-neck CT to a tightly-cropped T1c MRI (cross-modality, inter-session tilt up to ~30 deg) and the GTV is warped onto -the T1c native grid. Every pipeline is scored with the SAME crop-independent -geometric QA used in production (scripts/02_build_m6_dataset.py): +the T1c native grid. The headline alignment score is the achieved Mattes MI +(mi_ct_mr, nats) between the warped-CT and the original MR over the MR head -- +higher = better in-sync = better aligned. The production geometric QA +(scripts/02_build_m6_dataset.py) is still recorded per case in the JSON: - - skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR - head silhouette (mid/low level + angular std). - Aligned = 0 <= mid <= 20 mm and std <= 12 mm. + - mi_ct_mr : Mattes MI (nats), warped-CT ~ MR over the MR head + (the reported alignment score). - tv_in_ct_frac : fraction of the (dilated) warped GTV that falls inside the warped-CT head tissue (the slab must straddle it). - tv_nvox : warped GTV volume (>0.5 voxels). + - skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR + head silhouette (mid/low level + angular std). - accepted : production rule (nvox>0 and tv_in_ct and concentric). - reg_wall_s : wall-clock registration time (excludes I/O + scoring). -Because some toolkits are not pip-installable into the main -`longitudinal-313` env (antspyx is kept in an isolated env, see +Because some toolkits are not pip-installable into every +`longitudinal-313` env (antspyx needs a cp310-cp313 wheel, see REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under whichever conda env has a given toolkit installed, and it records only the pipelines that are importable there. elastix runs in the main env through the -ITK-Elastix python API (pip ``itk`` + ``itk-elastix``). ``synthmorph`` and -``impact`` run in Docker images (``freesurfer/synthmorph`` and -``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see -REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the +ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively +(``conda install -c conda-forge niftyreg`` provides ``reg_aladin``); +``synthmorph`` and ``impact`` run in Docker images (``freesurfer/synthmorph`` +and ``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, +see REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the per-pipeline JSON files (accumulated across envs) into one summary + plot. Usage (repo root): - # main env (SimpleITK + elastix + antspyx) + # main env (SimpleITK + elastix + antspyx + niftyreg) source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 - python scripts/bench_registration.py --pipelines simpleitk elastix antspyx + python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg # synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md python scripts/bench_registration.py --pipelines synthmorph impact @@ -82,7 +86,7 @@ def _mi_masked(fixed, moving, mask, bins=32): px, py = p.sum(axis=1, keepdims=True), p.sum(axis=0, keepdims=True) with np.errstate(divide="ignore", invalid="ignore"): ll = np.log(p / (px * py)) - return float(np.where((p > 0) & (px > 0) & (py > 0), p * ll, 0.0).sum()) + return float(np.where((p > 0) & (px > 0) & (py > 0), p * ll, 0.0).sum()) def _m6(): @@ -232,6 +236,8 @@ def _call(pipeline, ct, t1c, tv, args): if pipeline == "impact": return RP.register_impact(ct, t1c, model=args.impact_model, gpu=args.impact_gpu, iters=args.impact_iters, samples=args.impact_samples) + if pipeline == "niftyreg": + return RP.register_niftyreg(ct, t1c, threads=args.threads) raise KeyError(pipeline) @@ -255,10 +261,7 @@ def _run_case(pipeline, case, args): out.update(meta) out.update(evaluate(ct, t1c, tv, T)) if getattr(args, "qa_images", False): - g = out.get("skull_ring_gap_mm") - title = (f"{pipeline} | {key} | MI={out.get('mi_ct_mr')} | " - f"gap mid {g['mid_mm']}mm std {g['std_mm']}mm" if g - else f"{pipeline} | {key}") + title = f"{pipeline} | {key} | MI={out.get('mi_ct_mr')}" qp = os.path.join(args.qa_out, f"{pipeline}__{key}.png") try: wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear)) @@ -287,12 +290,8 @@ def run_pipeline(pipeline, cases, args): try: results.append(fu.result()) r = results[-1] - g = r.get("skull_ring_gap_mm") reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT" - print(f" {c['key']}: accepted={r['accepted']} " - f"gap={g['mid_mm'] if g else 'None'}mm " - f"std={g['std_mm'] if g else 'None'}mm " - f"mi={r.get('mi_ct_mr')} " + print(f" {c['key']}: mi={r.get('mi_ct_mr')} " f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True) except Exception as e: # noqa errors.append({"key": c["key"], "error": repr(e)}) @@ -348,27 +347,24 @@ def _fmt(v, spec="{}", default="-"): def print_table(rows): order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)] - print("\n" + "=" * 100) + print("\n" + "=" * 92) hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} " - f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}") + f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}") print(hdr) - print("-" * 100) + print("-" * 92) for p in order: r = next(x for x in rows if x["pipeline"] == p) if not r.get("available"): print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} " - f"{'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]") + f"{'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]") continue - gm = (r.get("skull_ring_gap_mid_mm") or {}).get("median") - gs = (r.get("skull_ring_gap_std_mm") or {}).get("median") rt = (r.get("reg_wall_s") or {}).get("median") tv = (r.get("tv_in_ct_frac") or {}).get("mean") mi = (r.get("mi_ct_mr") or {}).get("median") print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} " - f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {r['accepted_rate']:>7.2f} " - f"{_fmt(gm, '{:.1f}'):>7s} {_fmt(gs, '{:.1f}'):>7s} {_fmt(rt, '{:.1f}'):>7s} " + f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} " f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}") - print("=" * 100) + print("=" * 92) def write_plot(rows, out_png): @@ -384,19 +380,13 @@ def write_plot(rows, out_png): print("no available pipelines to plot") return labels = list(ok_rows) - gap_mid = [ok_rows[p].get("skull_ring_gap_mid_mm", {}).get("median") or 0 for p in labels] - gap_std = [ok_rows[p].get("skull_ring_gap_std_mm", {}).get("median") or 0 for p in labels] regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels] - acc = [ok_rows[p].get("accepted_rate") or 0 for p in labels] mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels] - fig, axes = plt.subplots(1, 5, figsize=(20, 4)) - colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD", "#E08E45"] - for ax, (title, vals, fmt, ylab, inv) in zip(axes, [ - ("Skull-ring gap mid (mm, median)", gap_mid, "{:.1f}", "mm", True), - ("Skull-ring gap std (mm, median)", gap_std, "{:.1f}", "mm", True), - ("Registration time (s, median)", regt, "{:.0f}", "s", True), - ("Acceptance rate", acc, "{:.2f}", "frac", False), - ("Mean Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats", False)]): + fig, axes = plt.subplots(1, 2, figsize=(10, 4)) + colors = ["#4878CF", "#6ACC65"] + for ax, (title, vals, fmt, ylab) in zip(axes, [ + ("Registration time (s, median)", regt, "{:.0f}", "s"), + ("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]): bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)]) for b, v in zip(bars, vals): ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()), @@ -426,10 +416,11 @@ def merge(args): print_table(rows) summary = { "timestamp": datetime.now().isoformat(timespec="seconds"), - "protocol": ("rigid CT->T1c; cross-modality; scored by production skull-ring " - "concentricity (aligned: 0<=mid<=20mm, std<=12mm) + GTV-in-CT fraction " - "(>=0.5); reg time excludes I/O+scoring; per-pipeline best-of-T vs T^-1 " - "orientation gate for elastix/antspyx"), + "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI " + "(nats, warped-CT ~ MR over the MR head, higher = better); reg time " + "excludes I/O+scoring; per-case JSON retains the production skull-ring " + "gap / acceptance QA; best-of-T vs T^-1 orientation gate for the " + "elastix-style toolkits"), "n_cases": next((r["n_cases"] for r in rows if r.get("available")), None), "case_keys": sorted({c["key"] for r in rows if r.get("available") for c in r.get("cases", [])}), "results": rows, @@ -440,18 +431,16 @@ def merge(args): print("summary:", sf) write_plot(rows, os.path.join(out_dir, "bench_registration.png")) # markdown table for the report - md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT | MI (nats) |", - "|---|---|---|---|---|---|---|---|---|---|---|---|"] + md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | tv-in-CT | MI (nats, median) |", + "|---|---|---|---|---|---|---|---|---|"] for p in [x["pipeline"] for x in rows]: r = next(x for x in rows if x["pipeline"] == p) if not r.get("available"): - md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | - | - | - | ({r['reason']}) |") + md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | ({r['reason']}) |") continue md.append( f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | " - f"{r['n_error']} | {r['accepted_rate']:.2f} | " - f"{(r.get('skull_ring_gap_mid_mm') or {}).get('median')} | " - f"{(r.get('skull_ring_gap_std_mm') or {}).get('median')} | " + f"{r['n_error']} | " f"{(r.get('reg_wall_s') or {}).get('median')} | " f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | " f"{(r.get('mi_ct_mr') or {}).get('median')} |") @@ -492,12 +481,12 @@ def parse_args(): ap.add_argument("--impact-gpu", type=int, default=None, help="GPU device index for IMPACT (default: auto-spread over " "the host GPUs by pid; -1 forces CPU)") - ap.add_argument("--impact-iters", type=int, default=500, - help="IMPACT per-resolution max iterations (default 500, the authors' " - "recommended Jacobian budget; lower to trade accuracy for speed)") - ap.add_argument("--impact-samples", type=int, default=2000, - help="IMPACT random spatial samples per iteration (default 2000, " - "the recommended budget)") + ap.add_argument("--impact-iters", type=int, default=2000, + help="IMPACT MaximumNumberOfIterations (default 2000, the benchmark " + "budget; ~10x the authors' 500-iter Jacobian budget, much slower)") + ap.add_argument("--impact-samples", type=int, default=5000, + help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the " + "elastix default)") ap.add_argument("--qa-images", action="store_true", help="also write a checkreg-style PNG (checkerboard blend + color " "overlay of the inverse-warped CT vs the original MR) per case") diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index 243f722..26f23fa 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -1,4 +1,4 @@ -"""Five rigid CT -> T1c registration pipelines behind one unified contract. +"""Six rigid CT -> T1c registration pipelines behind one unified contract. Problem (shared with ``scripts/02_build_m6_dataset.py``) -------------------------------------------------------- @@ -27,10 +27,11 @@ All pipelines register on the **full volumes** (no head box) and use each toolkit's **native automatic transform initialization** (no manual ``centroid_offset``): elastix / impact use elastix's ``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its -default center-of-mass alignment, simpleitk starts from identity, and -synthmorph needs none (global deep method). All use **fixed = CT, moving = -T1c** except IMPACT, which keeps **fixed = T1c, moving = CT** (swapping breaks -its TotalSegmentator metric -- see ``register_impact``). +default center-of-mass alignment, simpleitk starts from identity, niftyreg +(reg_aladin) initializes from the image centres, and synthmorph needs none +(global deep method). All use **fixed = CT, moving = T1c** except IMPACT, +which keeps **fixed = T1c, moving = CT** (swapping breaks its +TotalSegmentator metric -- see ``register_impact``). - ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). Default ``initial_transform=None`` computes a translation aligning the image @@ -44,6 +45,11 @@ its TotalSegmentator metric -- see ``register_impact``). GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start. Needs Docker + the image + a model asset. +- ``register_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` -- + symmetric block-matching + trimmed least squares (Modat et al., J. Med. + Img. 2014), NCC-based block similarity (cross-modality tolerant). Rigid + 6-DOF (``-rigOnly``), 3-level pyramid, 15 TLS iterations/level. Needs + ``conda install -c conda-forge niftyreg`` (provides ``reg_aladin``). - ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on the full volumes. Reuses the *proven* production ``register_ct_to_t1c`` in ``mode="production"``; the benchmark mode is a single identity-start @@ -71,6 +77,7 @@ applied to elastix / antspyx / synthmorph / impact, not to the trusted production SimpleITK path. """ import os +import re import shutil import subprocess import importlib.util @@ -162,14 +169,14 @@ def centroid_offset(ct, t1c, ct_box=None): # 2. SimpleITK (production baseline -- reused verbatim) # --------------------------------------------------------------------------- # -def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150): +def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000): """SimpleITK rigid CT->T1c, Mattes MI, full volumes (no head box). mode="benchmark" (default): a single multi-resolution pyramid run from the - identity start, capped at ``iters_full`` iterations at full resolution. - This is the *standard* SimpleITK algorithm and matches the budget the - other toolkits get in the benchmark (one rigid MI registration), so the - pipelines are compared on equal footing. No centroid init is used. + identity start, with ``iters_full`` (=2000) optimizer iterations at every + pyramid level. This is the *standard* SimpleITK algorithm and matches the + 2000-iteration budget the other toolkits get in the benchmark, so the + pipelines are compared on equal footing. No centroid init is used. mode="production": the exact multi-start tilt-grid ``register_ct_to_t1c`` from scripts/02_build_m6_dataset.py (robust to the ~30 deg inter-session tilt but much slower -- not the default for the benchmark). @@ -297,6 +304,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8): CT / T1c volumes (no head box). No manual centroid init: elastix's ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales. + Budget: ``MaximumNumberOfIterations`` 2000 per resolution; the sampler / + metric run at their defaults (``NumberOfSpatialSamples`` 5000, + ``NumberOfHistogramBins`` 32). Elastix 5.x component names: sampler ``Random`` (not ``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``). An output directory is required for the @@ -321,9 +331,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8): "AutomaticTransformInitializationMethod": "CenterOfGravity", "AutomaticScalesEstimation": True, "Optimizer": "AdaptiveStochasticGradientDescent", - "MaximumNumberOfIterations": [500, 500, 250], + "MaximumNumberOfIterations": [2000, 2000, 2000], "Metric": "AdvancedMattesMutualInformation", - "NumberOfHistogramBins": 10, + "NumberOfHistogramBins": 32, "ImageSampler": "Random", "ImageSamplerPatchSize": [25, 25, 25], "WriteResultImage": False, @@ -735,7 +745,7 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes, return "\n".join(L) + "\n" -def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, +def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000, threads=12, keep_scratch=False, scales=(6.0, 3.0, 1.5, 1.0), downsample=False): """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image @@ -761,11 +771,11 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, outside the T1c buffer and the IMPACT metric aborts ("too many samples map outside moving image buffer"). The returned transform still satisfies the unified contract (T1c->CT). - Budget: the default ``iters=500`` / ``samples=2000`` is the authors' - recommended Jacobian budget. It is *much* slower than the intensity-metric - pipelines (minutes per case on GPU vs seconds) because every iteration runs - the TotalSegmentator feature extractor on the sampled patches with - backpropagation. Lower ``iters``/``samples`` to trade accuracy for speed. + Budget: the default ``iters=2000`` (MaximumNumberOfIterations) / + ``samples=5000`` (NumberOfSpatialSamples, the elastix default) runs the + full benchmark budget for best results -- roughly 10x the authors' + 500/2000 Jacobian budget, i.e. hours per case on full volumes. Lower + ``iters``/``samples`` to trade accuracy for speed. The result is read back from ``TransformParameters.0.itk.txt`` (an ``Euler3DTransform`` in the world physical frame) and checked by :func:`fix_orientation`. @@ -886,10 +896,128 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, return T, meta +# --------------------------------------------------------------------------- # +# 6. NiftyReg (KCL-BMEIS suite, reg_aladin) +# --------------------------------------------------------------------------- # + + +def _niftyreg_bin(): + """Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable.""" + b = shutil.which("reg_aladin") + if b: + return b + raise PipelineUnavailable( + "niftyreg: reg_aladin not found on PATH. Install NiftyReg with " + "`conda install -c conda-forge niftyreg` (NiftyReg 1.5.69, the " + "KCL-BMEIS suite: reg_aladin rigid/affine + reg_f3d non-linear) and " + "re-run in the activated env.") + + +def _niftyreg_version(binpath): + try: + r = subprocess.run([binpath, "--version"], capture_output=True, + text=True, timeout=60) + m = re.search(r"(\d+\.\d+\.\d+)", r.stdout + r.stderr) + if m: + return f"NiftyReg {m.group(1)} (reg_aladin)" + except Exception: # noqa + pass + return "NiftyReg (reg_aladin)" + + +def _parse_niftyreg_affine(path): + """Parse reg_aladin's output affine (4x4, row-major) into a SimpleITK + AffineTransform. The file is *documented* as + ``Affine * Reference = Floating``; the direction that actually satisfies + the contract is verified empirically by :func:`fix_orientation` (on the + M6 cases the gate keeps the file's affine as-is).""" + M = np.array([float(x) for x in open(path).read().split()]).reshape(4, 4) + T = sitk.AffineTransform(3) + T.SetMatrix(M[:3, :3].ravel()) + T.SetTranslation(M[:3, 3]) + return T + + +def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False): + """NiftyReg rigid CT->T1c (reference=CT, floating=T1c) via ``reg_aladin``. + + reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global + block-matching registration: at each level of a ``levels``-level Gaussian + pyramid it matches 4-voxel blocks by normalised cross-correlation + (NCC -- a cross-modality-tolerant block similarity, unlike the histogram + MI objectives of the other pipelines), then solves a trimmed least squares + (TLS) problem for the rigid transform, iterating ``maxit`` times per + level. The symmetric variant (default) optimises forward and backward + transforms concurrently; ``-maxit 5`` (the tool default) is too coarse for + these repositioned pairs -- the benchmark uses 15. + + Full volumes (no head box); default initialisation (image centres); + ``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3`` (cubic) internal + warping. The output affine is documented as ``Affine * Reference = + Floating``; :func:`fix_orientation` resolves the direction that actually + satisfies the contract against the skull ring (on the M6 cases it keeps + the file's affine as-is, ``auto_inverted: False``) and returns the T1c -> + CT transform. + + Preconditions (else :class:`PipelineUnavailable`): ``reg_aladin`` on PATH + (``conda install -c conda-forge niftyreg``). + """ + import time + + binpath = _niftyreg_bin() + # full volumes, reference = CT, floating = T1c (no head box) + tmp = tempfile.mkdtemp(prefix="niftyreg_reg_") + dt = 0.0 + try: + sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32), + os.path.join(tmp, "ref.nii.gz"), True) + sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32), + os.path.join(tmp, "flo.nii.gz"), True) + cmd = [binpath, "-ref", "ref.nii.gz", "-flo", "flo.nii.gz", + "-rigOnly", "-aff", "aff.txt", + "-ln", str(levels), "-maxit", str(maxit), + "-interp", "3", "-omp", str(threads), "-voff"] + t0 = time.time() + proc = subprocess.run(cmd, capture_output=True, text=True, + cwd=tmp, timeout=7200) + dt = time.time() - t0 + if proc.returncode != 0: + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"niftyreg: reg_aladin failed (rc={proc.returncode}):\n{tail}") + aff = os.path.join(tmp, "aff.txt") + if not os.path.exists(aff): + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"niftyreg: no affine output produced:\n{tail}") + T = _parse_niftyreg_affine(aff) + # aff.txt is documented as CT -> T1c (Affine*Reference=Floating); + # the contract is T1c -> CT. fix_orientation resolves the direction + # empirically (on these cases it keeps the file's affine as-is). + T, inverted = fix_orientation(ct, t1c, T) + finally: + if not keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + + meta = { + "pipeline": "niftyreg", + "toolkit": _niftyreg_version(binpath), + "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)"), + "convention": "T1c->CT (verified)", "auto_inverted": inverted, + "fixed": "CT", "moving": "T1c", + "seconds_reg": round(dt, 1), + "levels": levels, "maxit": maxit, + } + if keep_scratch: + meta["scratch"] = tmp + return T, meta + + REGISTERS = { "antspyx": register_antspyx, "elastix": register_elastix, "impact": register_impact, + "niftyreg": register_niftyreg, "simpleitk": register_simpleitk, "synthmorph": register_synthmorph, } @@ -917,6 +1045,13 @@ def availability(): out[name] = (True, f"ants importable v{getattr(ants, '__version__', '?')}") except Exception as e: # noqa out[name] = (False, f"ants(antspyx) not importable ({type(e).__name__})") + elif name == "niftyreg": + b = shutil.which("reg_aladin") + if b: + out[name] = (True, _niftyreg_version(b)) + else: + out[name] = (False, "reg_aladin not on PATH " + "(conda install -c conda-forge niftyreg)") elif name == "synthmorph": ok, note = _synthmorph_docker_ok() if not ok: