feat(pipelines): integrate niftyreg and adopt Mattes MI as the benchmark headline score

Add `register_niftyreg` as the sixth rigid CT→T1c pipeline using `reg_aladin` (symmetric block matching + trimmed least squares, NCC similarity). Runs natively via `conda install -c conda-forge niftyreg`; no Docker required. Change the benchmark's primary reported metric from geometric skull-ring QA to the achieved Mattes MI (`mi_ct_mr`, nats) between the warped CT and the original MR over the MR head. Geometric QA fields are still recorded per case in the JSON output. Raise the SimpleITK benchmark-mode optimizer budget from 150 to 2000 iterations per pyramid level so all toolkits are compared on an equal iteration footing. Update `REGISTRATION_BENCHMARK.md` (five→six pipelines, niftyreg row, elastix bin count 10→32), `AGENTS.md` (toolkit list), and `bench_registration.py` docstrings/usage accordingly.
feat(benchmark): switch to per-toolkit native init and full-volume registration
2026-09-29 09:47:28 +08:00 · 2026-09-29 07:43:27 +08:00
4 changed files with 564 additions and 452 deletions
--- 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
--- 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
@ -22,34 +22,44 @@ sitk.Resample(tv,  t1c, T, sitk.sitkNearestNeighbor, 0) # GTV onto the MR grid
 ```
 places the CT / GTV concentric with the MR head. In `Resample`'s
-destination→source convention this `T` maps **T1c → CT**. All pipelines use
+destination→source convention this `T` maps **T1c → CT**. All pipelines
-**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c,
+register on the **full CT / T1c volumes** (no head box) and use each toolkit's
-moving = CT** (swapping its roles breaks the TotalSegmentator-based metric,
+**native automatic transform initialization** (no manual centroid offset):
-which aborts with "Error in metric" or converges to a wrong rotation). The
+elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
-intensity-metric pipelines register on the same CT-space head box (MI computed
+antspyx its default center-of-mass alignment, simpleitk starts from identity,
-on head tissue, not neck/air); `synthmorph` uses the full volumes (it
+niftyreg (reg_aladin) initializes from the image centres, synthmorph needs
-normalizes each one to an isotropic network space internally). The local
+none. All use **fixed = CT, moving = T1c** except `impact`, which keeps
-optimizers that cannot recover the gross inter-session offset from identity
+**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than
-(`antspyx`, `impact`) start from a centroid-aligned translation
+the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
-(`centroid_offset`); the others start from identity, so the comparison is on
+samples map outside the T1c buffer and the IMPACT metric aborts ("too many
-each toolkit's native optimizer / metric / global-search capability.
+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 use **fixed = CT, moving = T1c** except `impact`,
+Listed alphabetically. All register on the **full volumes** (no head box) and
-which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator
+use each toolkit's **native automatic transform initialization** (no manual
-metric). A pipeline that is a local optimizer and cannot recover the gross
+centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
-inter-session offset from identity (`antspyx`, `impact`) starts from a
+(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
-centroid-aligned translation; the others start from identity / their own global
+starts from identity, niftyreg initializes from the image centres, and
-search.
+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 | centroid init | metric | transform | status in this env |
+| # | pipeline | fixed / moving | init | metric | transform | status in this env |
-|---|----------|----------------|---------------|--------|-----------|--------------------|
+|---|----------|----------------|------|--------|-----------|--------------------|
-| 1 | `antspyx`    | CT / T1c | **yes** (required) | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
+| 1 | `antspyx`    | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
-| 2 | `elastix`    | CT / T1c | no (identity) | 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** | **yes** | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) |
+| 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 | no (identity) | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
+| 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 | `synthmorph` | CT / T1c | no (global deep) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
+| 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
@ -73,135 +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
+  (SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are
-  off as a standalone alignment score: higher = the warped-CT and MR intensities
+  better in-sync = better aligned; it collapses toward 0 when the CT is
-  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
-  rotated/translated away from the MR.
+  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) |
+Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
-|---|---|---|---|---|---|---|---|---|
+the MR head; higher = better).
 | **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 18.9 / 20.1 | 31.7 / 31.3 | 10.9 | 0.94 | **0.436 / 0.393** |
 | **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 16.6 / 17.6 | 16.9 / 18.8 | 2.4 | 0.93 | **0.641 / 0.583** |
 | **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 13.1 / 16.8 | 13.2 / 18.6 | **576.0** | 0.93 | **0.469 / 0.482** |
 | **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.7 / 16.0 | 13.5 / 15.7 | 7.6 | 0.93 | **0.645 / 0.632** |
 | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 22.9 / 21.7 | 26.0 | 0.93 | **0.241 / 0.308** |
-*(antspyx and impact are stochastic (random metric/patch subsampling), so their
+| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
 accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.6
 min/case on GPU, 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 |
 |---|---|---|---|---|---|
-| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
+| **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** |
-| m6_22IRRKU6_20240115 | 18.9 / 43.8 | 27.0 / 29.3 | 8.3 / 13.2 | 26.2 / 29.6 | 10.4 / 23.7 |
+| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** |
-| m6_23XIMO4M_20230222 | 22.9 / 24.5 | 22.6 / 23.9 | 23.1 / 25.7 | 22.6 / 24.0 | 20.3 / 23.3 |
+| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** |
-| m6_24V62QQM_20230223 | 18.7 / 47.5 | 11.3 / 16.9 | 9.9 / 7.9 **OK** | 10.2 / 16.2 | 10.4 / 22.9 |
+| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** |
-| m6_24V62QQM_20231017 | 15.3 / 37.9 | 16.6 / 11.1 **OK** | 6.6 / 6.9 **OK** | 11.9 / 7.7 **OK** | 11.5 / 8.4 **OK** |
+| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** |
-| m6_24V62QQM_20240820 | 38.7 / 54.0 | 10.2 / 8.2 **OK** | 37.1 / 45.6 | 10.8 / 9.0 **OK** | 43.4 / 15.8 |
+| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** |
-| m6_26TYQI3C_20221027 | 11.2 / 3.7 **OK** | 11.7 / 3.7 **OK** | 23.1 / 13.0 | 11.1 / 2.6 **OK** | 11.8 / 2.9 **OK** |
+
-| m6_2BBN5XA4_20250411 | 20.6 / 24.2 | 25.3 / 30.8 | 10.1 / 6.5 **OK** | 19.1 / 8.8 **OK** | 52.4 / 53.0 |
+\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default
-| m6_2C4TKI6H_20240506 | 15.3 / 14.2 | 14.2 / 15.3 | 13.1 / 17.4 | 12.7 / 13.5 | 16.7 / 28.2 |
+(MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
-| m6_2CDNN7BH_20240102 | 19.4 / 31.7 | 19.5 / 30.3 | 19.6 / 30.9 | 19.3 / 30.1 | 18.9 / 16.8 |
+(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 6–20 mm,
+- **SimpleITK and ANTsPy tie for best** (median MI 0.651 / 0.652): on the
-  within the `≤ 20 mm` band, except the grossly-drifted cases below).
+  well-positioned cases they cluster at 0.60–0.70 and both survive the harder
-  `tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for
+  repositioned cases (0.58–0.65 on `m6_2BBN5XA4` / `m6_22IRRKU6`). SimpleITK's
-  every pipeline, so the target is in the right place.
+  per-case MI is **unchanged** from the previous 150-iteration budget — RSGD
- **The binding failure is rotation** — `std > 12 mm` on the repositioned cases
+  converges well before the 2000-iteration cap, so the cap is now headroom
-  (the ~30° inter-session tilt produces local minima that a single start
+  rather than a constraint.
-  (identity or centroid) can lock onto). This is exactly why production
+- **The grossly repositioned cases define the ranking.** On
-  `02_build_m6_dataset.py` uses a **multi-start tilt grid**.
+  `m6_2BBN5XA4` and `m6_24V62QQM_20240820` (~30° inter-session tilt) the tools
- **SimpleITK is the most reliable intensity-metric tool in this draw**
+  split: ANTs/SimpleITK hold 0.58–0.65, elastix drops to 0.29–0.64, niftyreg
-  (4/10 accepted, best median mid 12.7), and **Elastix is the fastest**
+  0.07–0.12, synthmorph 0.06–0.10 — single-start optimizers locking onto local
-  (2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best
+  minima. This is exactly why production `02_build_m6_dataset.py` uses a
-  median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c
+  **multi-start tilt grid**.
-  convention the two settle into different local minima on the hard tilted
+- **Elastix is second at the new budget** (2000 iters per resolution, 32 bins,
-  cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017`
+  5000 samples): median 0.516, below the 0.641 recorded at the old 500/500/250
-  (std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles
+  + 10-bin setting. The long ASGD run nails the clean cases (0.64–0.70) but its
-  better, and vice-versa on the drift cases. The fixed/moving direction
+  single CenterOfGravity start still locks onto local minima on the
-  materially affects where the ASGD/RSGD optimizers converge on the tilted cases.
+  repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
- **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this
+  0.416).
-  draw; it is stochastic, so the count varies): it frequently locks onto the
+- **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI
-  tilted local minimum (median std 31.7, worst of all). It accepts only
+  on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618),
-  `m6_26TYQI3C` (mid 11.2 / std 3.7) and misses `m6_24V62QQM_20230223`
+  where anatomical-segmentation similarity out-resolves intensity MI — but its
-  (std 47.5) and `m6_24V62QQM_20231017` (std 37.9) that the other tools accept —
+  median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is
-  textbook tilted local minima.
+  also by far the slowest (~14 min/case at the shown budget; the new 2000/5000
- **IMPACT (semantic metric) is the standout on the hardest repositioned case**
+  budget would be on the order of hours/case).
-  and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every
+- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
-  intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of
+  is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class;
-  any pipeline and within ~1 mm of the acceptance threshold — because matching
+  `m6_26TYQI3C` 0.608) but falls into the same tilt local minima
-  *anatomical semantics* (TotalSegmentator features) is far more discriminating
+  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower
-  for gross pose than raw cross-modality intensity. But the semantic landscape
+  than elastix at full volumes.
-  is rugged: single-init ASGD lands in bad local minima on several cases, and
+- **SynthMorph is weakest of the six on these cross-modality CT pairs**
-  two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0).
+  (median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30°
-  Net 3/10 accepted — between ANTs and SimpleITK — at ~9.6 min/case, two orders
+  tilt defeats it on the repositioned cases (0.06–0.10) even though it matches
-  of magnitude slower than the intensity tools.
+  the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.4 s
+  centroid init** (global deep method) and is deterministic.
-   (warm workers), SimpleITK ~7.6 s, ANTs ~10.9 s (CPU-bound),
+- **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the
-   **SynthMorph ~26 s** (GPU: a single pretrained-model forward pass per case),
+  CT head everywhere, so coverage is not the discriminator — MI is.
-   **IMPACT ~576 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator
+- **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s,
-   extractor on the sampled patches with backprop).
+  SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s,
-  - **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**:
+  **NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous
-    medians 0.645 (SimpleITK) / 0.567 (Elastix) / 0.469 (IMPACT) / 0.436 (ANTs) /
+  budget; GPU feature extraction with backprop per iteration).
-    0.241 (SynthMorph),
+- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
   and it diverges sharply on the tilted cases — ANTs `m6_24V62QQM_20230223`
   drops to ~0.06 nats (rotation ~47 mm off) while SimpleITK holds 0.703. A high
   MI means the CT skull/tissue intensities fall on the matching MR structures
   after the inverse transform. ANTs' and SynthMorph's low median MI is
   consistent with their tilted-case convergence failures.
 - **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.241)
    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
  per-case numbers carry run-to-run variance (this table is one draw);
-  SimpleITK's MI (and registration) is deterministic at fixed thread count
+  SimpleITK's MI (and registration) is deterministic at fixed thread count;
-  (RSGD wall time is thread-sensitive, convergence is not); the in-process
+  the in-process Elastix ASGD is deterministic given its fixed seed;
-  Elastix ASGD is deterministic given its fixed seed. **IMPACT is the most
+  SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**:
-  stochastic**: its Jacobian mode draws random patches every iteration, so the
+  its Jacobian mode draws random patches every iteration, so the final pose
-  final pose — and even *which* cases get accepted — changes run to run (two
+  changes run to run. Treat its per-case numbers as one sample from a
-  independent 10-case draws both accepted 3/10 but on different cases). Treat
+  distribution.
  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
@ -217,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 <int|-1>`,
 `--synthmorph-extent <int>`, `--impact-model <rel|abs>`,
-`--impact-gpu <int|-1>`, `--impact-iters N`, `--impact-samples N`,
+`--impact-gpu <int|-1>`, `--impact-iters N` (default 2000),
 `--impact-samples N` (default 5000),
 `--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
 ## Visual QA images (checkerboard + color overlay)
@ -242,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
+The title carries the pipeline, case key and the achieved Mattes MI, so a
-(mid/std), so a glance tells you both the achieved-metric score and whether the
+glance tells you whether the visuals agree with the score. Generate for the
-visuals agree. Generate for the available pipelines with:
+available pipelines with:
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
@ -296,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)
@ -303,8 +338,10 @@ image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121)
 *semantic* similarity metric compiled in as a native Elastix metric, plus a
 TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
 Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
-**fixed = T1c, moving = CT** (the reverse of the others): swapping its roles
+**fixed = T1c, moving = CT** (the reverse of the others): the CT FOV
-breaks the TotalSegmentator metric ("Error in metric" / wrong rotation).
+(~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").
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
@ -330,49 +367,50 @@ repo for the full list and FOV/layer tables.
 ### How the adapter drives the image
-The image's `CMD.sh` only runs `elastix -f -m -p -out` (no `-t0`), so the
+The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
-adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`)
+binary directly (`/lib/elastix-install/bin/elastix`):
 to add the shared centered start:
-1. **Head box** — resamples fixed=T1c and moving=CT into the same CT-space head
+1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box).
-   box (physical frame) the other pipelines use, so MI/IMPACT is computed on
+   T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too
-   head tissue, not neck/air.
+   many CT samples map outside the T1c buffer and the metric aborts.
-2. **Shared init via `-t0`** — the `centroid_offset` translation is written to a
+2. **Auto init (no `-t0`)** — the parameter map sets
-   `TranslationTransform` parameter file and passed with `-t0`; with
+   `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
-   `HowToCombineTransforms "Compose"` the final transform is
+   elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
-   `T = T_euler(T_init(p))`. (Elastix 5.x removed in-file multi-`Transform`
+   optimizer scales.
   chains and `UseScalingFactors`, so the initial transform must be a separate
   `-t0` file.)
 3. **Parameter map** — rigid `EulerTransform` (6-DOF), `Impact` metric in
   **Jacobian** mode (differentiable, per the recommended config), 4 semantic
-   scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the box's base spacing) with
+   scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the T1c base spacing) with
   **no image resampling** (all-1 pyramid schedules, as the IMPACT docs require),
   ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid
   (`--impact-gpu -1` forces CPU).
 4. **Read-back** — the result `TransformParameters.0.itk.txt` is an
-   `Euler3DTransform` in the box/world physical frame; the adapter composes it
+    `Euler3DTransform` in the world physical frame (the
-   with the `-t0` translation and runs the usual `fix_orientation` gate.
+    `AutomaticTransformInitialization` is folded in by elastix); the adapter
    runs it through the usual `fix_orientation` gate.
 ## Notes / caveats
- The benchmark compares each toolkit in its **native configuration**: the
+- The benchmark compares each toolkit in its **native configuration**: all
-  local optimizers that cannot recover the gross offset from identity
+  register on the full volumes (no head box) and use each toolkit's native
-  (`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`,
+  automatic transform initialization (no manual `centroid_offset`) — elastix /
-  `elastix` and `synthmorph` start from identity / their own global search.
+  impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's
-  Production uses `--simpleitk-mode production`'s multi-start tilt grid, which
+  center-of-mass alignment, simpleitk identity, synthmorph none. Production
-  raises SimpleITK's acceptance at ~20× the cost; that is out of scope for the
+  uses `--simpleitk-mode production`'s multi-start tilt grid, which raises
  SimpleITK's acceptance at ~20× the cost; that is out of scope for the
  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
+  IMPACT. Its low MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
-  weakness, not variance.
+  genuine cross-modality / tilt weakness, not variance.
- `fix_orientation` (elastix/ants/synthmorph/impact) chooses between a transform
+- `fix_orientation` (elastix/ants/niftyreg/synthmorph/impact) chooses between a
-  and its inverse by skull-ring concentricity — a correctness gate that makes
+  transform and its inverse by skull-ring concentricity — a correctness gate
-  the result independent of each toolkit's transform-direction convention. With
+  that makes the result independent of each toolkit's transform-direction
-  the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the
+  convention. With the fixed=CT/moving=T1c swap, `elastix`'s combination
-  inverse of the contract, so it reports `auto_inverted: True` (the gate
+  transform is now the inverse of the contract, so it reports
-  corrects it); antspyx/synthmorph/impact return the contract direction
+  `auto_inverted: True` (the gate corrects it); antspyx/synthmorph/impact
-  directly (`auto_inverted: False`).
+  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
@ -385,8 +423,27 @@ to add the shared centered start:
      path).
    - **GPU sharing.** The container runs as the host uid with
      `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
-    - **~26 s/case** (GPU: one pretrained forward pass), between ANTs (~11 s)
+    - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK
-      and IMPACT (~576 s).
+      (~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
@ -394,11 +451,11 @@ to add the shared centered start:
    the component set), the metric key is `Metric` (not `ImageMetric`), the
    random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x
    names and is not compatible with elastix ≤ 4.x parameter files.
-  - No initial transform is passed (identity start): the ASGD + 3-resolution
+  - `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the
-    Mattes-MI pyramid recovers the gross inter-session offset on its own, so
+    start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the
-    the shared `centroid_offset` start is not required. An **output directory**
+    optimizer scales. An **output directory** is still set
-    is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter
+    (`SetOutputDirectory`) for the C++ filter to run; the adapter uses a
-    uses a scratch dir under `data/m6reg/elastix/itkelastix/`.
+    scratch dir under `data/m6reg/elastix/itkelastix/`.
  - `WriteResultImage` must be the string `"false"` (Python `False`); the
    default (absent) makes the filter raise "Could not read result image"
    after a *successful* registration.
@ -408,32 +465,40 @@ to add the shared centered start:
    (`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
+  - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
-    500/500/250 over 3 resolutions).
+    (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
+      extractor on the sampled patches with backpropagation. The benchmark
-     500-iter / 2000-sample budget is ~9.6 min/case (vs 2–11 s for the other
+      budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
-     tools). `--impact-iters` / `--impact-samples` trade accuracy for speed.
+      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: swapping its roles
+     follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251
-     breaks the TotalSegmentator metric (elastix aborts with "Error in metric"
+     mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the
-     or converges to a wrong rotation). The returned transform still satisfies
+     fixed image too many CT samples map outside the T1c buffer and the metric
-     the contract (T1c→CT).
+     aborts ("too many samples map outside moving image buffer"). The returned
     transform still satisfies the contract (T1c→CT).
   - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns
     the intensity centroids; no manual `-t0` translation is used.
   - **Rugged single-init convergence.** The semantic landscape has more local
-     minima than Mattes-MI: the pipeline nails the hardest repositioned case
+     minima than Mattes-MI: on full volumes it no longer nails the hardest
-     (`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some
+     repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box
-     others. A multi-start or higher-iteration budget would likely raise its
+     used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration
-     acceptance; that is out of scope for the single-init head-to-head.
+     budget would likely raise its acceptance; that is out of scope for the
     single-init head-to-head.
   - **Stochastic.** Jacobian mode resamples patches each iteration, so results
     vary run to run (see Stochasticity).
   - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image
     bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors`
-    are gone — the shared centered start is a separate `-t0` translation file,
+     are gone. The result is read back as an `Euler3DTransform` (`.itk.txt`).
-    and the result is read back as an `Euler3DTransform` (`.itk.txt`) and
+     The image's own `CMD.sh` is a generic launcher, so the adapter invokes
-    composed with `-t0`. The image's own `CMD.sh` doesn't pass `-t0`, so the
+     `/lib/elastix-install/bin/elastix` directly.
    adapter invokes `/lib/elastix-install/bin/elastix` directly.
  - **Models are not in the image.** The TorchScript model(s) must be present
    under the mounted `/Data/Models` (the repo caches them in
    `data/impact_models/`, bind-mounted read-only).
--- 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
+the T1c native grid.  The headline alignment score is the achieved Mattes MI
-geometric QA used in production (scripts/02_build_m6_dataset.py):
+(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
+  - mi_ct_mr           : Mattes MI (nats), warped-CT ~ MR over the MR head
-                         head silhouette (mid/low level + angular std).
+                         (the reported alignment score).
                         Aligned = 0 <= mid <= 20 mm and std <= 12 mm.
  - 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
+Because some toolkits are not pip-installable into every
-`longitudinal-313` env (antspyx is kept in an isolated env, see
+`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
+ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively
-``impact`` run in Docker images (``freesurfer/synthmorph`` and
+(``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
-``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see
+``synthmorph`` and ``impact`` run in Docker images (``freesurfer/synthmorph``
-REGISTRATION_BENCHMARK.md).  The `--merge` mode then combines the
+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
@ -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')}"
        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}")
        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']} "
+                print(f"   {c['key']}: mi={r.get('mi_ct_mr')} "
                      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')} "
                      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"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} "
              f"{_fmt(gm, '{:.1f}'):>7s} {_fmt(gs, '{:.1f}'):>7s} {_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))
+    fig, axes = plt.subplots(1, 2, figsize=(10, 4))
-    colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD", "#E08E45"]
+    colors = ["#4878CF", "#6ACC65"]
-    for ax, (title, vals, fmt, ylab, inv) in zip(axes, [
+    for ax, (title, vals, fmt, ylab) in zip(axes, [
-        ("Skull-ring gap mid (mm, median)", gap_mid, "{:.1f}", "mm", True),
+        ("Registration time (s, median)", regt, "{:.0f}", "s"),
-        ("Skull-ring gap std (mm, median)", gap_std, "{:.1f}", "mm", True),
+        ("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]):
        ("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)]):
        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 "
+        "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI "
-                     "concentricity (aligned: 0<=mid<=20mm, std<=12mm) + GTV-in-CT fraction "
+                      "(nats, warped-CT ~ MR over the MR head, higher = better); reg time "
-                     "(>=0.5); reg time excludes I/O+scoring; per-pipeline best-of-T vs T^-1 "
+                      "excludes I/O+scoring; per-case JSON retains the production skull-ring "
-                     "orientation gate for elastix/antspyx"),
+                      "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['n_error']} | "
            f"{(r.get('skull_ring_gap_mid_mm') or {}).get('median')} | "
            f"{(r.get('skull_ring_gap_std_mm') or {}).get('median')} | "
            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,
+    ap.add_argument("--impact-iters", type=int, default=2000,
-                    help="IMPACT per-resolution max iterations (default 500, the authors' "
+                    help="IMPACT MaximumNumberOfIterations (default 2000, the benchmark "
-                         "recommended Jacobian budget; lower to trade accuracy for speed)")
+                         "budget; ~10x the authors' 500-iter Jacobian budget, much slower)")
-    ap.add_argument("--impact-samples", type=int, default=2000,
+    ap.add_argument("--impact-samples", type=int, default=5000,
-                    help="IMPACT random spatial samples per iteration (default 2000, "
+                    help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
-                         "the recommended budget)")
+                         "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")
--- 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``)
 --------------------------------------------------------
@ -23,34 +23,41 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
 Pipelines (listed alphabetically)
 ---------------------------------
-All pipelines use the shared convention **fixed = CT, moving = T1c**, except
+All pipelines register on the **full volumes** (no head box) and use each
-IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its
+toolkit's **native automatic transform initialization** (no manual
-TotalSegmentator metric -- see ``register_impact``).  A pipeline that is a
+``centroid_offset``): elastix / impact use elastix's
-local optimizer and cannot recover the gross inter-session offset from
+``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its
-identity (antspyx, impact) starts from a centroid-aligned translation; the
+default center-of-mass alignment, simpleitk starts from identity, niftyreg
-others (simpleitk, elastix, synthmorph) do not require it and start from
+(reg_aladin) initializes from the image centres, and synthmorph needs none
-identity / their own global search.
+(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).
-  Centroid pre-shift of the moving (T1c) image is *required* (ANTs
+  Default ``initial_transform=None`` computes a translation aligning the image
-  Rigid-from-identity collapses on the gross offset).
+  centers of mass.
 - ``register_elastix``    : ITK-Elastix Python API (pip ``itk`` +
  ``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
-  Identity start (the ASGD + multi-resolution pyramid recovers the offset).
+  ``AutomaticTransformInitialization`` (CenterOfGravity) start.
 - ``register_impact``     : IMPACT (semantic-similarity metric) inside the
  ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript
  TotalSegmentator feature extractor as a native Elastix metric).  Rigid 6-DOF,
-  GPU (Jacobian mode).  Kept as fixed=T1c/moving=CT (see above); centroid init
+  GPU (Jacobian mode).  Kept as fixed=T1c/moving=CT (see above); elastix
-  kept (local ASGD on the semantic metric benefits from the gross alignment).
+  ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
  Needs Docker + the image + a model asset.
- ``register_simpleitk``  : SimpleITK ImageRegistrationMethod, Mattes MI on a
+- ``register_niftyreg``   : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
-  CT-space head box.  Reuses the *proven* production ``register_ct_to_t1c`` in
+  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
  multi-resolution pyramid run.
 - ``register_synthmorph`` : FreeSurfer SynthMorph (TensorFlow, pretrained
  anatomy-aware contrast-invariant model) inside the ``freesurfer/synthmorph``
  Docker image.  Rigid 6-DOF matrix model, GPU when available.  Global deep
-  method -- no centroid init.  Needs Docker + the pulled image.
+  method -- no init.  Needs Docker + the pulled image.
 A toolkit that is not installed (or, for the Docker-based pipelines, missing
 its image)
@ -70,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
@ -157,38 +165,18 @@ def centroid_offset(ct, t1c, ct_box=None):
    return (float(t1c_c[0] - c[0]), float(t1c_c[1] - c[1]), float(t1c_c[2] - c[2]))
 def _ct_head_box_mm(ct):
    """CT-space head box in physical mm: (box_lo, box_hi) tuples, or None."""
    m6 = _load_m6()
    ca = sitk.GetArrayFromImage(ct)
    tissue = ca > -150
    csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
    vz = m6._ct_vertex_slice(ct)
    if vz is None:
        return None, None
    ztop = co[2] + vz * csp[2]
    band = tissue[max(0, vz - int(round(150 / csp[2]))):vz + 1]
    cols = np.nonzero(band.any(0))
    box_lo = (co[0] + cols[1].min() * csp[0] - 15,
              co[1] + cols[0].min() * csp[1] - 15, ztop - 215)
    box_hi = (co[0] + cols[1].max() * csp[0] + 15,
              co[1] + cols[0].max() * csp[1] + 15, ztop + 10)
    return box_lo, box_hi
 # --------------------------------------------------------------------------- #
 # 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 on a CT-space head box.
+    """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.
+    identity start, with ``iters_full`` (=2000) optimizer iterations at every
-        This is the *standard* SimpleITK algorithm and matches the budget the
+    pyramid level.  This is the *standard* SimpleITK algorithm and matches the
-        other toolkits get in the benchmark (one rigid MI registration), so the
+    2000-iteration budget the other toolkits get in the benchmark, so the
-        pipelines are compared on equal footing.  The pyramid recovers the gross
+    pipelines are compared on equal footing.  No centroid init is used.
        inter-session offset, so 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).
@ -215,20 +203,16 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
    except AttributeError:
        pass
    ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
-    lo, hi = _ct_head_box_mm(ct)
+    # full volumes, fixed = CT, moving = T1c (no head box)
-    if lo is None:
+    fixed, moving = ct, t1c
-        raise PipelineUnavailable("simpleitk: no CT head box found (corrupt CT?)")
+    shape = np.array(ct.GetSize())
-    tsp = np.array(t1c.GetSpacing())
+    # Identity start (no centroid init); the multi-resolution Mattes-MI pyramid
-    shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
+    # recovers the gross inter-session offset.  The Euler center is the fixed
-    ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
+    # (CT) geometric center.
-    ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
+    init = sitk.Euler3DTransform()
-    ref.SetOrigin(tuple(float(x) for x in lo))
+    init.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint(
-    fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
+        ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2,
-    moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
+         (fixed.GetSize()[2] - 1) // 2))))
    # Identity start: the multi-resolution Mattes-MI pyramid recovers the gross
    # inter-session offset on its own, so no centroid init is required here.
    ct_c = np.array([(lo[i] + hi[i]) / 2.0 for i in range(3)])
    init = sitk.Euler3DTransform(); init.SetCenter(tuple(float(x) for x in ct_c))
    init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
    import time
    t0 = time.time()
@ -250,9 +234,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
    meta = {
        "pipeline": "simpleitk", "mode": "benchmark",
        "toolkit": f"SimpleITK {sitk.Version_VersionString()}",
-        "metric": "MattesMutualInformation (head box, CT space)",
+        "metric": "MattesMutualInformation (full volumes)",
        "transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
-                      f"{iters_full} iters/level)"),
+                      f"{iters_full} iters/level, identity start)"),
        "convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
        "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
        "skull_ring_gap_mm": gap,
@ -311,34 +295,29 @@ def _itk_transform_to_sitk_rigid(itk_tf):
 def register_elastix(ct, t1c, workdir=None, threads=8):
-    """ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), Mattes MI.
+    """ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), 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
-    (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common
+    (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the full
-    CT-space head box, from the identity start (the ASGD + multi-resolution
+    CT / T1c volumes (no head box).  No manual centroid init: elastix's
-    pyramid recovers the gross inter-session offset, so no centroid init is
+    ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity
-    required).  Elastix 5.x component names: sampler ``Random`` (not
+    centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales.
-    ``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``).  An output
+    Budget: ``MaximumNumberOfIterations`` 2000 per resolution; the sampler /
-    directory is required for the C++ filter to run, and ``WriteResultImage``
+    metric run at their defaults (``NumberOfSpatialSamples`` 5000,
-    must be the string ``"false"`` (a Python bool False) or the filter raises
+    ``NumberOfHistogramBins`` 32).
-    "Could not read result image" after a successful registration.
+    Elastix 5.x
    component names: sampler ``Random`` (not ``RandomSampler``), metric key
    ``Metric`` (not ``ImageMetric``).  An output directory is required for the
    C++ filter to run, and ``WriteResultImage`` must be the string ``"false"``
    (a Python bool False) or the filter raises "Could not read result image"
    after a successful registration.
    """
    itk = _elastix_api()
-    lo, hi = _ct_head_box_mm(ct)
+    # full volumes, fixed = CT, moving = T1c (no head box)
-    if lo is None:
+    fixed = ct
-        raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)")
+    moving = t1c
    # Resample both into a common CT-space head box (mirrors the other pipelines)
    # so MI is computed on head tissue, not neck/air.
    tsp = np.array(t1c.GetSpacing())
    shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
    ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
    ref.SetSpacing(t1c.GetSpacing())
    ref.SetDirection(t1c.GetDirection())
    ref.SetOrigin(tuple(float(x) for x in lo))
    fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)     # CT  = reference
    moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)   # T1c = moving
    params = {
        "ImageDimension": 3,
        "FixedInternalImagePixelType": "float",
@ -348,10 +327,13 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
        "FixedCenteredTransform": True,
        "TransformParameterReadoutFormat": 1,
        "HowToCombineTransforms": "Compose",
        "AutomaticTransformInitialization": True,
        "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,
@ -388,8 +370,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
        version = "ITK-Elastix (python API)"
    meta = {
        "pipeline": "elastix", "toolkit": version,
-        "metric": "MattesMutualInformation (head box, CT space)",
+        "metric": "MattesMutualInformation (full volumes)",
-        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)",
+        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)",
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
        "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
    }
@ -402,18 +384,16 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
 def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
                     aff_shrink_factors=None, aff_smoothing_sigmas=None,
-                     aff_metric="mattes", use_head_box=True, center_init=True):
+                     aff_metric="mattes"):
-    """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c).
+    """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c), full volumes.
    antspyx 0.6.3: ``ants.registration(fixed, moving, type_of_transform='Rigid',
    aff_metric='mattes', aff_iterations=..., aff_shrink_factors=...,
    aff_smoothing_sigmas=...)``.  ``aff_metric='mattes'`` (the default) gives a
    multi-resolution Mattes-MI rigid registration -- the cross-modality metric
-    this CT/MR pair needs.  To keep the benchmark tractable and consistent with
+    this CT/MR pair needs.  Registers on the full CT / T1c volumes (no head box).
-    the other pipelines (and to avoid neck/air in the metric), both images are
+    No manual centroid init: antspyx's default ``initial_transform=None``
-    first resampled into the same CT-space head box.  ANTs Rigid-from-identity
+    computes a translation that aligns the image centers of mass.  Returns the
    does NOT recover the gross inter-session offset, so a centroid pre-shift of
    the moving (T1c) image is required (see ``center_init``).  Returns the
    transform that warps CT into the T1c grid (orientation verified by
    :func:`fix_orientation`).
    """
@ -428,36 +408,12 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
    aff_iterations = tuple(aff_iterations) if aff_iterations else (500, 500, 200)
    aff_shrink_factors = tuple(aff_shrink_factors) if aff_shrink_factors else (4, 2, 1)
    aff_smoothing_sigmas = tuple(aff_smoothing_sigmas) if aff_smoothing_sigmas else (3, 2, 0)
-    if use_head_box:
+    # full volumes, fixed = CT, moving = T1c (no head box)
        lo, hi = _ct_head_box_mm(ct)
        if lo is None:
            raise PipelineUnavailable("antspyx: no CT head box found (corrupt CT?)")
        tsp = np.array(t1c.GetSpacing())
        shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
        ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
        ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
        ref.SetOrigin(tuple(float(x) for x in lo))
        fx = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)      # CT  box (fixed)
        mv = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)     # T1c box (moving)
    else:
    fx, mv = ct, t1c
    import time
    t0 = time.time()
-    # Centered start (REQUIRED for ANTs): antspyx cannot use a
+    # No manual centroid init: antspyx's default initial_transform=None computes
-    # TranslationTransform as initial_transform (the ANTs CLI rejects it), and
+    # a translation aligning the image centers of mass.
    # Rigid-from-identity does not recover the large (~100-130 mm) inter-session
    # offset -- without this init ANTs collapses on 8/9 cases.  So we pre-shift the
    # T1c (moving) box so its head centroid coincides with the CT (fixed) head
    # centroid, run Rigid from identity, then compose the shift + inverse back so
    # the result is expressed in the original T1c/CT coordinates.
    off = None
    if center_init:
        off = centroid_offset(ct, t1c, ct_box=fx)
        if off is not None:
            # shift the T1c content by -off (head c_t1c -> c_ct): resample with +off
            tr_off = sitk.TranslationTransform(3)
            tr_off.SetOffset((off[0], off[1], off[2]))
            mv = sitk.Resample(mv, mv, tr_off, sitk.sitkLinear, 0.0)
    af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32))
    am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32))
    r = ants.registration(fixed=af, moving=am, type_of_transform=transform_type,
@ -466,20 +422,15 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
                          aff_smoothing_sigmas=aff_smoothing_sigmas)
    dt = time.time() - t0
    # r['fwdtransforms'][0] is the ANTs forward transform: fixed(CT) -> moving
-    # (pre-shifted T1c).  We need the contract T1c(orig) -> CT: map T1c(orig) to
+    # (T1c).  The contract is T1c -> CT, i.e. its inverse; fix_orientation
-    # the pre-shifted T1c (translation by -off), then invert the ANTs transform.
+    # verifies the direction against the skull ring.
    T_ants = sitk.ReadTransform(r["fwdtransforms"][0])
    if off is not None:
        tr_shift = sitk.TranslationTransform(3)
        tr_shift.SetOffset((-off[0], -off[1], -off[2]))          # T1c orig -> pre
        cand = sitk.CompositeTransform([tr_shift, T_ants.GetInverse()])
    else:
    cand = T_ants.GetInverse()
    T, inverted = fix_orientation(ct, t1c, cand)
    meta = {
        "pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
-        "metric": f"ANTs {aff_metric} (multi-resolution)",
+        "metric": f"ANTs {aff_metric} (multi-resolution, full volumes)",
-        "transform": f"ANTs {transform_type} (head box)", "convention": "T1c->CT (verified)",
+        "transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)",
        "auto_inverted": inverted, "aff_iterations": list(aff_iterations),
        "aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
    }
@ -740,12 +691,18 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
        "(ImpactFeaturesMapUpdateInterval -1)",
        '(ImpactWriteFeatureMaps "false")',
        "",
-        "// rigid 6-DOF to be optimized (the shared centered start comes via -t0)",
+        "// rigid 6-DOF to be optimized (auto center init, no -t0)",
        '(Transform "EulerTransform")',
        '(TransformBase "ElastixTransformBase")',
        '(TransformParameterReadoutFormat "AllParametersOneLine")',
        '(FixedCenteredTransform "true")',
        '(HowToCombineTransforms "Compose")',
        '(AutomaticTransformInitialization "true")',
        # align intensity centroids (heads), not buffer corners -- the CT FOV is
        # much larger than the T1c crop, so GeometricalCenter leaves the heads
        # ~100 mm apart and the metric aborts ("too many samples outside buffer")
        '(AutomaticTransformInitializationMethod "CenterOfGravity")',
        '(AutomaticScalesEstimation "true")',
        "",
        "(ImageDimension 3)",
        '(FixedInternalImagePixelType "float")',
@ -788,19 +745,12 @@ 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
    ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
    NOTE on the fixed/moving convention: every other pipeline here uses
    fixed=CT, moving=T1c.  IMPACT is the exception -- swapping its fixed/moving
    to CT/T1c breaks the TotalSegmentator-based semantic metric (elastix aborts
    with "Error in metric" or converges to a wrong rotation), so IMPACT keeps
    the working MRI-as-fixed (T1c) / CT-as-moving convention.  The returned
    transform still satisfies the unified contract (T1c->CT).
    IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
    one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run
    on randomly sampled patches of the fixed and moving images and the
@ -809,21 +759,26 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
    the Mattes-MI pipelines above -- it is the cross-modality tool the image was
    built for.
-    To stay on the benchmark's equal footing the pipeline:
+    The pipeline registers on the full CT / T1c volumes (no head box) with no
-      - registers on the same CT-space head box (both images resampled into it),
+    manual centroid init: elastix's ``AutomaticTransformInitialization``
-      - starts from the *shared* centroid-aligned translation (passed to Elastix
+    (CenterOfGravity) aligns the intensity centroids.  It uses a rigid 6-DOF
-        as a fixed initial transform via ``-t0``; with ``HowToCombineTransforms
+    ``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level
-        Compose`` the final transform = ``T_main(T_init(p))``),
+    scale pyramid (no image resampling).
-      - uses a rigid 6-DOF ``EulerTransform`` (the 5.x name for the classic
+
-        3D rigid), ASGD, a 4-level scale pyramid (no image resampling).
+    NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other
-    Budget: the default ``iters=500`` / ``samples=2000`` is the authors'
+    pipelines).  The CT FOV (~300x300x251 mm) is much larger than the T1c crop
-    recommended Jacobian budget.  It is *much* slower than the intensity-metric
+    (~220x220x150 mm); with CT as the fixed image, too many CT samples map
-    pipelines (minutes per case on GPU vs seconds) because every iteration runs
+    outside the T1c buffer and the IMPACT metric aborts ("too many samples map
-    the TotalSegmentator feature extractor on the sampled patches with
+    outside moving image buffer").  The returned transform still satisfies the
-    backpropagation.  Lower ``iters``/``samples`` to trade accuracy for speed.
+    unified contract (T1c->CT).
    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 box/world physical frame) and composed with
+    (an ``Euler3DTransform`` in the world physical frame) and checked by
-    the ``-t0`` translation, then checked by :func:`fix_orientation`.
+    :func:`fix_orientation`.
    Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, the
    ``vboussot/elastix_impact`` image, and the model asset under the cache.  A
@ -837,25 +792,13 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
        raise PipelineUnavailable(f"impact: {note}")
    model_host, model_rel = _impact_model_host(model)
-    lo, hi = _ct_head_box_mm(ct)
+    # full volumes, fixed = T1c, moving = CT (no head box, no manual init --
-    if lo is None:
+    # elastix's AutomaticTransformInitialization aligns the image centers).
-        raise PipelineUnavailable("impact: no CT head box found (corrupt CT?)")
+    # T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger
-
+    # than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples
-    # resample both into the shared CT-space head box (mirrors the other pipelines)
+    # map outside the T1c buffer and the IMPACT metric aborts.
-    tsp = np.array(t1c.GetSpacing())
+    fixed = t1c
-    shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
+    moving = ct
    ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
    ref.SetSpacing(t1c.GetSpacing())
    ref.SetDirection(t1c.GetDirection())
    ref.SetOrigin(tuple(float(x) for x in lo))
    fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)    # T1c = reference
    moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)     # CT  = moving
    # shared centered start (REQUIRED: the local ASGD optimizer does not recover
    # the gross offset from identity).  Elastix maps fixed(CT)->moving(T1c), the
    # same direction as centroid_offset, so no negation.
    off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0)
    init = (-float(off[0]), -float(off[1]), -float(off[2]))
    # GPU: auto-spread across the host GPUs by pid, or forced value
    if gpu is None:
@ -901,13 +844,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
            f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
                                       gpu, iters, samples, voxel_sizes,
                                       rescale=schedule, smoothing=smoothe))
        with open(os.path.join(datadir, "InitialTransform.txt"), "w") as f:
            f.write(
                '(Transform "TranslationTransform")\n'
                '(TransformBase "ElastixTransformBase")\n'
                '(TransformParameterReadoutFormat "AllParametersOneLine")\n'
                "(NumberOfParameters 3)\n"
                f"(TransformParameters {init[0]:g} {init[1]:g} {init[2]:g})\n")
        cmd = ["docker", "run", "--rm"]
        if gpu >= 0:
@ -919,7 +855,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
            IMPACT_IMAGE, IMPACT_ELASTIX_BIN,
            "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha",
            "-p", "/Data/ParameterMap.txt",
            "-t0", "/Data/InitialTransform.txt",
            "-out", "/Out/", "-threads", str(threads),
        ]
        t0 = time.time()
@ -934,14 +869,12 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
            tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
            raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
-        # the optimized Euler (box/world frame) composed with the fixed -t0
+        # the optimized Euler (world frame) is the full combination transform
-        # translation: T = T_euler(T_init(p)).  Both map fixed(CT)->moving(T1c),
+        # (AutomaticTransformInitialization is folded in by elastix); it maps
-        # i.e. the inverse of the contract (T1c->CT); fix_orientation resolves it.
+        # fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT);
        # fix_orientation resolves the direction.
        T_main = sitk.ReadTransform(itk_tf)
-        T_init = sitk.TranslationTransform(3)
+        T, inverted = fix_orientation(ct, t1c, T_main)
        T_init.SetOffset(init)
        cand = sitk.CompositeTransform([T_init, T_main])
        T, inverted = fix_orientation(ct, t1c, cand)
    finally:
        if not keep_scratch:
            shutil.rmtree(tmp, ignore_errors=True)
@ -950,7 +883,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
        "pipeline": "impact",
        "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
        "metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)",
-        "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, centroid init; "
+        "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; "
                      f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
        "fixed": "T1c", "moving": "CT",
@ -963,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,
 }
@ -994,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: