4 changed files with 452 additions and 564 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. 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`.
+Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK.
 ## Project
--- a/REGISTRATION_BENCHMARK.md
+++ b/REGISTRATION_BENCHMARK.md
@ -1,6 +1,6 @@
 # Registration Pipeline Benchmark — CT → T1c (M6)
-Benchmarks **six rigid, cross-modality registration pipelines** that align a
+Benchmarks **five 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,44 +22,34 @@ 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
+destination→source convention this `T` maps **T1c → CT**. All pipelines use
-register on the **full CT / T1c volumes** (no head box) and use each toolkit's
+**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c,
-**native automatic transform initialization** (no manual centroid offset):
+moving = CT** (swapping its roles breaks the TotalSegmentator-based metric,
-elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
+which aborts with "Error in metric" or converges to a wrong rotation). The
-antspyx its default center-of-mass alignment, simpleitk starts from identity,
+intensity-metric pipelines register on the same CT-space head box (MI computed
-niftyreg (reg_aladin) initializes from the image centres, synthmorph needs
+on head tissue, not neck/air); `synthmorph` uses the full volumes (it
-none. All use **fixed = CT, moving = T1c** except `impact`, which keeps
+normalizes each one to an isotropic network space internally). The local
-**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than
+optimizers that cannot recover the gross inter-session offset from identity
-the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
+(`antspyx`, `impact`) start from a centroid-aligned translation
-samples map outside the T1c buffer and the IMPACT metric aborts ("too many
+(`centroid_offset`); the others start from identity, so the comparison is on
-samples map outside moving image buffer"). The comparison is on each
+each toolkit's native optimizer / metric / global-search capability.
 toolkit's native optimizer / metric / global-search capability.
-## The six pipelines
+## The five pipelines
-Listed alphabetically. All register on the **full volumes** (no head box) and
+Listed alphabetically. All use **fixed = CT, moving = T1c** except `impact`,
-use each toolkit's **native automatic transform initialization** (no manual
+which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator
-centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
+metric). A pipeline that is a local optimizer and cannot recover the gross
-(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
+inter-session offset from identity (`antspyx`, `impact`) starts from a
-starts from identity, niftyreg initializes from the image centres, and
+centroid-aligned translation; the others start from identity / their own global
-synthmorph needs none (global deep method). All use **fixed = CT, moving =
+search.
 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 |
+| # | pipeline | fixed / moving | centroid 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) |
+| 1 | `antspyx`    | CT / T1c | **yes** (required) | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
-| 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) |
+| 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) |
-| 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) |
+| 3 | `impact`     | **T1c / CT** | **yes** | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ 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) |
+| 4 | `simpleitk`  | CT / T1c | no (identity) | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
-| 5 | `simpleitk`  | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
+| 5 | `synthmorph` | CT / T1c | no (global deep) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
 | 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
@ -83,123 +73,135 @@ m6_2BBN5XA4_20250411  m6_2C4TKI6H_20240506  m6_2CDNN7BH_20240102
 ## Metrics (identical for every pipeline)
-The headline alignment score is **`mi_ct_mr`** — the registration's *achieved*
+Production QA, crop-independent:
 Mattes MI, read off as a standalone score:
 - **`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). 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
+  CT) at mid and low levels, plus the angular std. **Aligned ⇔
-  translational concentricity; `std` measures rotational concentricity.
+  `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; aligned ⇔ `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`).
+  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.
 - **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring).
 ## Results (10 cases)
-Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
+| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
-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** |
-| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
+*(antspyx and impact are stochastic (random metric/patch subsampling), so their
 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 |
 |---|---|---|---|---|---|
-| **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** |
+| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
-| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** |
+| m6_22IRRKU6_20240115 | 18.9 / 43.8 | 27.0 / 29.3 | 8.3 / 13.2 | 26.2 / 29.6 | 10.4 / 23.7 |
-| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** |
+| m6_23XIMO4M_20230222 | 22.9 / 24.5 | 22.6 / 23.9 | 23.1 / 25.7 | 22.6 / 24.0 | 20.3 / 23.3 |
-| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** |
+| m6_24V62QQM_20230223 | 18.7 / 47.5 | 11.3 / 16.9 | 9.9 / 7.9 **OK** | 10.2 / 16.2 | 10.4 / 22.9 |
-| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** |
+| 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** |
-| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** |
+| 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 |
-
+| 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** |
-\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default
+| 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 |
-(MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
+| m6_2C4TKI6H_20240506 | 15.3 / 14.2 | 14.2 / 15.3 | 13.1 / 17.4 | 12.7 / 13.5 | 16.7 / 28.2 |
-(hours/case) and was not re-run.
+| m6_2CDNN7BH_20240102 | 19.4 / 31.7 | 19.5 / 30.3 | 19.6 / 30.9 | 19.3 / 30.1 | 18.9 / 16.8 |
 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
- **SimpleITK and ANTsPy tie for best** (median MI 0.651 / 0.652): on the
+- **Translation is recovered well by all five** (mid gaps mostly 6–20 mm,
-  well-positioned cases they cluster at 0.60–0.70 and both survive the harder
+  within the `≤ 20 mm` band, except the grossly-drifted cases below).
-  repositioned cases (0.58–0.65 on `m6_2BBN5XA4` / `m6_22IRRKU6`). SimpleITK's
+  `tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for
-  per-case MI is **unchanged** from the previous 150-iteration budget — RSGD
+  every pipeline, so the target is in the right place.
-  converges well before the 2000-iteration cap, so the cap is now headroom
+- **The binding failure is rotation** — `std > 12 mm` on the repositioned cases
-  rather than a constraint.
+  (the ~30° inter-session tilt produces local minima that a single start
- **The grossly repositioned cases define the ranking.** On
+  (identity or centroid) can lock onto). This is exactly why production
-  `m6_2BBN5XA4` and `m6_24V62QQM_20240820` (~30° inter-session tilt) the tools
+  `02_build_m6_dataset.py` uses a **multi-start tilt grid**.
-  split: ANTs/SimpleITK hold 0.58–0.65, elastix drops to 0.29–0.64, niftyreg
+- **SimpleITK is the most reliable intensity-metric tool in this draw**
-  0.07–0.12, synthmorph 0.06–0.10 — single-start optimizers locking onto local
+  (4/10 accepted, best median mid 12.7), and **Elastix is the fastest**
-  minima. This is exactly why production `02_build_m6_dataset.py` uses a
+  (2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best
-  **multi-start tilt grid**.
+  median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c
- **Elastix is second at the new budget** (2000 iters per resolution, 32 bins,
+  convention the two settle into different local minima on the hard tilted
-  5000 samples): median 0.516, below the 0.641 recorded at the old 500/500/250
+  cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017`
-  + 10-bin setting. The long ASGD run nails the clean cases (0.64–0.70) but its
+  (std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles
-  single CenterOfGravity start still locks onto local minima on the
+  better, and vice-versa on the drift cases. The fixed/moving direction
-  repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
+  materially affects where the ASGD/RSGD optimizers converge on the tilted cases.
-  0.416).
+- **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this
- **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI
+  draw; it is stochastic, so the count varies): it frequently locks onto the
-  on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618),
+  tilted local minimum (median std 31.7, worst of all). It accepts only
-  where anatomical-segmentation similarity out-resolves intensity MI — but its
+  `m6_26TYQI3C` (mid 11.2 / std 3.7) and misses `m6_24V62QQM_20230223`
-  median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is
+  (std 47.5) and `m6_24V62QQM_20231017` (std 37.9) that the other tools accept —
-  also by far the slowest (~14 min/case at the shown budget; the new 2000/5000
+  textbook tilted local minima.
-  budget would be on the order of hours/case).
+- **IMPACT (semantic metric) is the standout on the hardest repositioned case**
- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
+  and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every
-  is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class;
+  intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of
-  `m6_26TYQI3C` 0.608) but falls into the same tilt local minima
+  any pipeline and within ~1 mm of the acceptance threshold — because matching
-  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower
+  *anatomical semantics* (TotalSegmentator features) is far more discriminating
-  than elastix at full volumes.
+  for gross pose than raw cross-modality intensity. But the semantic landscape
- **SynthMorph is weakest of the six on these cross-modality CT pairs**
+  is rugged: single-init ASGD lands in bad local minima on several cases, and
-  (median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30°
+  two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0).
-  tilt defeats it on the repositioned cases (0.06–0.10) even though it matches
+  Net 3/10 accepted — between ANTs and SimpleITK — at ~9.6 min/case, two orders
-  the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no
+  of magnitude slower than the intensity tools.
-  centroid init** (global deep method) and is deterministic.
+- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.4 s
- **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the
+   (warm workers), SimpleITK ~7.6 s, ANTs ~10.9 s (CPU-bound),
-  CT head everywhere, so coverage is not the discriminator — MI is.
+   **SynthMorph ~26 s** (GPU: a single pretrained-model forward pass per case),
- **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s,
+   **IMPACT ~576 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator
-  SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s,
+   extractor on the sampled patches with backprop).
-  **NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous
+  - **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**:
-  budget; GPU feature extraction with backprop per iteration).
+    medians 0.645 (SimpleITK) / 0.567 (Elastix) / 0.469 (IMPACT) / 0.436 (ANTs) /
- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
+    0.241 (SynthMorph),
   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
-  the in-process Elastix ASGD is deterministic given its fixed seed;
+  (RSGD wall time is thread-sensitive, convergence is not); the in-process
-  SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**:
+  Elastix ASGD is deterministic given its fixed seed. **IMPACT is the most
-  its Jacobian mode draws random patches every iteration, so the final pose
+  stochastic**: its Jacobian mode draws random patches every iteration, so the
-  changes run to run. Treat its per-case numbers as one sample from a
+  final pose — and even *which* cases get accepted — changes run to run (two
-  distribution.
+  independent 10-case draws both accepted 3/10 but on different cases). Treat
  its per-case numbers as one sample from a distribution.
 ## Reproducing
 ```bash
-# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (pip/conda-installed here)
+# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here)
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
-python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg \
+python scripts/bench_registration.py --pipelines simpleitk elastix antspyx \
    --n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
 # Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
@ -215,12 +217,10 @@ 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` (also the OpenMP thread count for `reg_aladin`),
+`--threads T`, `--simpleitk-mode {benchmark,production}`,
 `--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` (default 2000),
+`--impact-gpu <int|-1>`, `--impact-iters N`, `--impact-samples N`,
 `--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 and the achieved Mattes MI, so a
+The title carries the pipeline, case key, Mattes MI, and the skull-ring gap
-glance tells you whether the visuals agree with the score. Generate for the
+(mid/std), so a glance tells you both the achieved-metric score and whether the
-available pipelines with:
+visuals agree. Generate for the available pipelines with:
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
@ -296,41 +296,6 @@ 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)
@ -338,10 +303,8 @@ 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): the CT FOV
+**fixed = T1c, moving = CT** (the reverse of the others): swapping its roles
-(~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with
+breaks the TotalSegmentator metric ("Error in metric" / wrong rotation).
 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):
@ -367,50 +330,49 @@ repo for the full list and FOV/layer tables.
 ### How the adapter drives the image
-The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
+The image's `CMD.sh` only runs `elastix -f -m -p -out` (no `-t0`), so the
-binary directly (`/lib/elastix-install/bin/elastix`):
+adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`)
 to add the shared centered start:
-1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box).
+1. **Head box** — resamples fixed=T1c and moving=CT into the same CT-space head
-   T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too
+   box (physical frame) the other pipelines use, so MI/IMPACT is computed on
-   many CT samples map outside the T1c buffer and the metric aborts.
+   head tissue, not neck/air.
-2. **Auto init (no `-t0`)** — the parameter map sets
+2. **Shared init via `-t0`** — the `centroid_offset` translation is written to a
-   `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
+   `TranslationTransform` parameter file and passed with `-t0`; with
-   elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
+   `HowToCombineTransforms "Compose"` the final transform is
-   optimizer scales.
+   `T = T_euler(T_init(p))`. (Elastix 5.x removed in-file multi-`Transform`
   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 T1c base spacing) with
+   scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the box's 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 world physical frame (the
+   `Euler3DTransform` in the box/world physical frame; the adapter composes it
-    `AutomaticTransformInitialization` is folded in by elastix); the adapter
+   with the `-t0` translation and runs the usual `fix_orientation` gate.
    runs it through the usual `fix_orientation` gate.
 ## Notes / caveats
- The benchmark compares each toolkit in its **native configuration**: all
+- The benchmark compares each toolkit in its **native configuration**: the
-  register on the full volumes (no head box) and use each toolkit's native
+  local optimizers that cannot recover the gross offset from identity
-  automatic transform initialization (no manual `centroid_offset`) — elastix /
+  (`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`,
-  impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's
+  `elastix` and `synthmorph` start from identity / their own global search.
-  center-of-mass alignment, simpleitk identity, synthmorph none. Production
+  Production uses `--simpleitk-mode production`'s multi-start tilt grid, which
-  uses `--simpleitk-mode production`'s multi-start tilt grid, which raises
+  raises SimpleITK's acceptance at ~20× the cost; that is out of scope for the
  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 MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
+  IMPACT. Its low acceptance (2/10) is a genuine cross-modality / tilt
-  genuine cross-modality / tilt weakness, not variance.
+  weakness, not variance.
- `fix_orientation` (elastix/ants/niftyreg/synthmorph/impact) chooses between a
+- `fix_orientation` (elastix/ants/synthmorph/impact) chooses between a transform
-  transform and its inverse by skull-ring concentricity — a correctness gate
+  and its inverse by skull-ring concentricity — a correctness gate that makes
-  that makes the result independent of each toolkit's transform-direction
+  the result independent of each toolkit's transform-direction convention. With
-  convention. With the fixed=CT/moving=T1c swap, `elastix`'s combination
+  the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the
-  transform is now the inverse of the contract, so it reports
+  inverse of the contract, so it reports `auto_inverted: True` (the gate
-  `auto_inverted: True` (the gate corrects it); antspyx/synthmorph/impact
+  corrects it); antspyx/synthmorph/impact return the contract direction
-  return the contract direction directly (`auto_inverted: False`). `niftyreg`'s
+  directly (`auto_inverted: False`).
  `-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
@ -423,27 +385,8 @@ binary directly (`/lib/elastix-install/bin/elastix`):
      path).
    - **GPU sharing.** The container runs as the host uid with
      `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
-    - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK
+    - **~26 s/case** (GPU: one pretrained forward pass), between ANTs (~11 s)
-      (~10 s) / ANTs (~37 s) and IMPACT (~834 s).
+      and IMPACT (~576 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
@ -451,11 +394,11 @@ binary directly (`/lib/elastix-install/bin/elastix`):
    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.
-  - `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the
+  - No initial transform is passed (identity start): the ASGD + 3-resolution
-    start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the
+    Mattes-MI pyramid recovers the gross inter-session offset on its own, so
-    optimizer scales. An **output directory** is still set
+    the shared `centroid_offset` start is not required. An **output directory**
-    (`SetOutputDirectory`) for the C++ filter to run; the adapter uses a
+    is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter
-    scratch dir under `data/m6reg/elastix/itkelastix/`.
+    uses a 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.
@ -465,40 +408,32 @@ 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)`.
-  - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
+  - Registration wall time is ~2.4 s/case median at full budget (ASGD
-    (ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their
+    500/500/250 over 3 resolutions).
    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. The benchmark
+     extractor on the sampled patches with backpropagation, so the recommended
-      budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
+     500-iter / 2000-sample budget is ~9.6 min/case (vs 2–11 s for the other
-      the elastix default) is ~10× the authors' 500-iter / 2000-sample
+     tools). `--impact-iters` / `--impact-samples` trade accuracy for speed.
      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
+     follow the shared fixed=CT/moving=T1c convention: swapping its roles
-     mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the
+     breaks the TotalSegmentator metric (elastix aborts with "Error in metric"
-     fixed image too many CT samples map outside the T1c buffer and the metric
+     or converges to a wrong rotation). The returned transform still satisfies
-     aborts ("too many samples map outside moving image buffer"). The returned
+     the contract (T1c→CT).
     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: on full volumes it no longer nails the hardest
+     minima than Mattes-MI: the pipeline nails the hardest repositioned case
-     repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box
+     (`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some
-     used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration
+     others. A multi-start or higher-iteration budget would likely raise its
-     budget would likely raise its acceptance; that is out of scope for the
+     acceptance; that is out of scope for the single-init head-to-head.
-     single-init head-to-head.
+  - **Stochastic.** Jacobian mode resamples patches each iteration, so results
-   - **Stochastic.** Jacobian mode resamples patches each iteration, so results
+    vary run to run (see Stochasticity).
-     vary run to run (see Stochasticity).
+  - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image
-   - **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`
-     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,41 +1,37 @@
-"""Benchmark the 6 rigid CT->T1c registration pipelines on N M6 test cases.
+"""Benchmark the 5 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 | niftyreg
+  simpleitk | elastix | antspyx | synthmorph | impact
 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.  The headline alignment score is the achieved Mattes MI
+the T1c native grid.  Every pipeline is scored with the SAME crop-independent
-(mi_ct_mr, nats) between the warped-CT and the original MR over the MR head --
+geometric QA used in production (scripts/02_build_m6_dataset.py):
 higher = better in-sync = better aligned.  The production geometric QA
 (scripts/02_build_m6_dataset.py) is still recorded per case in the JSON:
-  - mi_ct_mr           : Mattes MI (nats), warped-CT ~ MR over the MR head
+  - skull_ring_gap_mm  : concentricity of the warped-CT skull ring vs the MR
-                         (the reported alignment score).
+                         head silhouette (mid/low level + angular std).
                         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 every
+Because some toolkits are not pip-installable into the main
-`longitudinal-313` env (antspyx needs a cp310-cp313 wheel, see
+`longitudinal-313` env (antspyx is kept in an isolated env, 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``); niftyreg runs natively
+ITK-Elastix python API (pip ``itk`` + ``itk-elastix``).  ``synthmorph`` and
-(``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
+``impact`` run in Docker images (``freesurfer/synthmorph`` and
-``synthmorph`` and ``impact`` run in Docker images (``freesurfer/synthmorph``
+``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see
-and ``vboussot/elastix_impact`` respectively; needs Docker + the pulled image,
+REGISTRATION_BENCHMARK.md).  The `--merge` mode then combines the
 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 + niftyreg)
+  # main env (SimpleITK + elastix + antspyx)
  source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
-  python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg
+  python scripts/bench_registration.py --pipelines simpleitk elastix antspyx
  # synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
  python scripts/bench_registration.py --pipelines synthmorph impact
@ -86,7 +82,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():
@ -236,8 +232,6 @@ 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)
@ -261,7 +255,10 @@ def _run_case(pipeline, case, args):
    out.update(meta)
    out.update(evaluate(ct, t1c, tv, T))
    if getattr(args, "qa_images", False):
-        title = f"{pipeline} | {key} | MI={out.get('mi_ct_mr')}"
+        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}")
        qp = os.path.join(args.qa_out, f"{pipeline}__{key}.png")
        try:
            wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
@ -290,8 +287,12 @@ 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']}: mi={r.get('mi_ct_mr')} "
+                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')} "
                      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)})
@ -347,24 +348,27 @@ 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" + "=" * 92)
+    print("\n" + "=" * 100)
    hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} "
-           f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
+            f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
    print(hdr)
-    print("-" * 92)
+    print("-" * 100)
    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}   [{r['reason']}]")
+                  f"{'-':>7s} {'-':>7s} {'-':>7s} {'-':>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} {_fmt(rt, '{:.1f}'):>7s} "
+              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"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
-    print("=" * 92)
+    print("=" * 100)
 def write_plot(rows, out_png):
@ -380,13 +384,19 @@ 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, 2, figsize=(10, 4))
+    fig, axes = plt.subplots(1, 5, figsize=(20, 4))
-    colors = ["#4878CF", "#6ACC65"]
+    colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD", "#E08E45"]
-    for ax, (title, vals, fmt, ylab) in zip(axes, [
+    for ax, (title, vals, fmt, ylab, inv) in zip(axes, [
-        ("Registration time (s, median)", regt, "{:.0f}", "s"),
+        ("Skull-ring gap mid (mm, median)", gap_mid, "{:.1f}", "mm", True),
-        ("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]):
+        ("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)]):
        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()),
@ -416,11 +426,10 @@ def merge(args):
    print_table(rows)
    summary = {
        "timestamp": datetime.now().isoformat(timespec="seconds"),
-        "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI "
+        "protocol": ("rigid CT->T1c; cross-modality; scored by production skull-ring "
-                      "(nats, warped-CT ~ MR over the MR head, higher = better); reg time "
+                     "concentricity (aligned: 0<=mid<=20mm, std<=12mm) + GTV-in-CT fraction "
-                      "excludes I/O+scoring; per-case JSON retains the production skull-ring "
+                     "(>=0.5); reg time excludes I/O+scoring; per-pipeline best-of-T vs T^-1 "
-                      "gap / acceptance QA; best-of-T vs T^-1 orientation gate for the "
+                     "orientation gate for elastix/antspyx"),
                      "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,
@ -431,16 +440,18 @@ 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 | s/reg (median) | tv-in-CT | MI (nats, median) |",
+    md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT | MI (nats) |",
-          "|---|---|---|---|---|---|---|---|---|"]
+          "|---|---|---|---|---|---|---|---|---|---|---|---|"]
    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']} | "
+            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.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')} |")
@ -481,12 +492,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=2000,
+    ap.add_argument("--impact-iters", type=int, default=500,
-                    help="IMPACT MaximumNumberOfIterations (default 2000, the benchmark "
+                    help="IMPACT per-resolution max iterations (default 500, the authors' "
-                         "budget; ~10x the authors' 500-iter Jacobian budget, much slower)")
+                         "recommended Jacobian budget; lower to trade accuracy for speed)")
-    ap.add_argument("--impact-samples", type=int, default=5000,
+    ap.add_argument("--impact-samples", type=int, default=2000,
-                    help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
+                    help="IMPACT random spatial samples per iteration (default 2000, "
-                         "elastix default)")
+                         "the recommended budget)")
    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 @@
-"""Six rigid CT -> T1c registration pipelines behind one unified contract.
+"""Five rigid CT -> T1c registration pipelines behind one unified contract.
 Problem (shared with ``scripts/02_build_m6_dataset.py``)
 --------------------------------------------------------
@ -23,41 +23,34 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
 Pipelines (listed alphabetically)
 ---------------------------------
-All pipelines register on the **full volumes** (no head box) and use each
+All pipelines use the shared convention **fixed = CT, moving = T1c**, except
-toolkit's **native automatic transform initialization** (no manual
+IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its
-``centroid_offset``): elastix / impact use elastix's
+TotalSegmentator metric -- see ``register_impact``).  A pipeline that is a
-``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its
+local optimizer and cannot recover the gross inter-session offset from
-default center-of-mass alignment, simpleitk starts from identity, niftyreg
+identity (antspyx, impact) starts from a centroid-aligned translation; the
-(reg_aladin) initializes from the image centres, and synthmorph needs none
+others (simpleitk, elastix, synthmorph) do not require it and start from
-(global deep method).  All use **fixed = CT, moving = T1c** except IMPACT,
+identity / their own global search.
 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
+  Centroid pre-shift of the moving (T1c) image is *required* (ANTs
-  centers of mass.
+  Rigid-from-identity collapses on the gross offset).
 - ``register_elastix``    : ITK-Elastix Python API (pip ``itk`` +
  ``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
-  ``AutomaticTransformInitialization`` (CenterOfGravity) start.
+  Identity start (the ASGD + multi-resolution pyramid recovers the offset).
 - ``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); elastix
+  GPU (Jacobian mode).  Kept as fixed=T1c/moving=CT (see above); centroid init
-  ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
+  kept (local ASGD on the semantic metric benefits from the gross alignment).
  Needs Docker + the image + a model asset.
- ``register_niftyreg``   : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
+- ``register_simpleitk``  : SimpleITK ImageRegistrationMethod, Mattes MI on a
-  symmetric block-matching + trimmed least squares (Modat et al., J. Med.
+  CT-space head box.  Reuses the *proven* production ``register_ct_to_t1c`` in
  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 init.  Needs Docker + the pulled image.
+  method -- no centroid init.  Needs Docker + the pulled image.
 A toolkit that is not installed (or, for the Docker-based pipelines, missing
 its image)
@ -77,7 +70,6 @@ applied to elastix / antspyx / synthmorph / impact, not to the trusted
 production SimpleITK path.
 """
 import os
 import re
 import shutil
 import subprocess
 import importlib.util
@ -165,18 +157,38 @@ 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=2000):
+def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
-    """SimpleITK rigid CT->T1c, Mattes MI, full volumes (no head box).
+    """SimpleITK rigid CT->T1c, Mattes MI on a CT-space head box.
    mode="benchmark" (default): a single multi-resolution pyramid run from the
-    identity start, with ``iters_full`` (=2000) optimizer iterations at every
+    identity start, capped at ``iters_full`` iterations at full resolution.
-    pyramid level.  This is the *standard* SimpleITK algorithm and matches the
+        This is the *standard* SimpleITK algorithm and matches the budget the
-    2000-iteration budget the other toolkits get in the benchmark, so the
+        other toolkits get in the benchmark (one rigid MI registration), so the
-    pipelines are compared on equal footing.  No centroid init is used.
+        pipelines are compared on equal footing.  The pyramid recovers the gross
        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).
@ -203,16 +215,20 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
    except AttributeError:
        pass
    ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
-    # full volumes, fixed = CT, moving = T1c (no head box)
+    lo, hi = _ct_head_box_mm(ct)
-    fixed, moving = ct, t1c
+    if lo is None:
-    shape = np.array(ct.GetSize())
+        raise PipelineUnavailable("simpleitk: no CT head box found (corrupt CT?)")
-    # Identity start (no centroid init); the multi-resolution Mattes-MI pyramid
+    tsp = np.array(t1c.GetSpacing())
-    # recovers the gross inter-session offset.  The Euler center is the fixed
+    shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
-    # (CT) geometric center.
+    ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
-    init = sitk.Euler3DTransform()
+    ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
-    init.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint(
+    ref.SetOrigin(tuple(float(x) for x in lo))
-        ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2,
+    fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
-         (fixed.GetSize()[2] - 1) // 2))))
+    moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
    # 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()
@ -234,9 +250,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
    meta = {
        "pipeline": "simpleitk", "mode": "benchmark",
        "toolkit": f"SimpleITK {sitk.Version_VersionString()}",
-        "metric": "MattesMutualInformation (full volumes)",
+        "metric": "MattesMutualInformation (head box, CT space)",
        "transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
-                      f"{iters_full} iters/level, identity start)"),
+                      f"{iters_full} iters/level)"),
        "convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
        "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
        "skull_ring_gap_mm": gap,
@ -295,29 +311,34 @@ 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, full images.
+    """ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), Mattes MI.
    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 full
+    (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common
-    CT / T1c volumes (no head box).  No manual centroid init: elastix's
+    CT-space head box, from the identity start (the ASGD + multi-resolution
-    ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity
+    pyramid recovers the gross inter-session offset, so no centroid init is
-    centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales.
+    required).  Elastix 5.x component names: sampler ``Random`` (not
-    Budget: ``MaximumNumberOfIterations`` 2000 per resolution; the sampler /
+    ``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``).  An output
-    metric run at their defaults (``NumberOfSpatialSamples`` 5000,
+    directory is required for the C++ filter to run, and ``WriteResultImage``
-    ``NumberOfHistogramBins`` 32).
+    must be the string ``"false"`` (a Python bool False) or the filter raises
-    Elastix 5.x
+    "Could not read result image" after a successful registration.
    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()
-    # full volumes, fixed = CT, moving = T1c (no head box)
+    lo, hi = _ct_head_box_mm(ct)
-    fixed = ct
+    if lo is None:
-    moving = t1c
+        raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)")
    # 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",
@ -327,13 +348,10 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
        "FixedCenteredTransform": True,
        "TransformParameterReadoutFormat": 1,
        "HowToCombineTransforms": "Compose",
        "AutomaticTransformInitialization": True,
        "AutomaticTransformInitializationMethod": "CenterOfGravity",
        "AutomaticScalesEstimation": True,
        "Optimizer": "AdaptiveStochasticGradientDescent",
-        "MaximumNumberOfIterations": [2000, 2000, 2000],
+        "MaximumNumberOfIterations": [500, 500, 250],
        "Metric": "AdvancedMattesMutualInformation",
-        "NumberOfHistogramBins": 32,
+        "NumberOfHistogramBins": 10,
        "ImageSampler": "Random",
        "ImageSamplerPatchSize": [25, 25, 25],
        "WriteResultImage": False,
@ -370,8 +388,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
        version = "ITK-Elastix (python API)"
    meta = {
        "pipeline": "elastix", "toolkit": version,
-        "metric": "MattesMutualInformation (full volumes)",
+        "metric": "MattesMutualInformation (head box, CT space)",
-        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)",
+        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)",
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
        "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
    }
@ -384,16 +402,18 @@ 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"):
+                     aff_metric="mattes", use_head_box=True, center_init=True):
-    """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c), full volumes.
+    """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c).
    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.  Registers on the full CT / T1c volumes (no head box).
+    this CT/MR pair needs.  To keep the benchmark tractable and consistent with
-    No manual centroid init: antspyx's default ``initial_transform=None``
+    the other pipelines (and to avoid neck/air in the metric), both images are
-    computes a translation that aligns the image centers of mass.  Returns the
+    first resampled into the same CT-space head box.  ANTs Rigid-from-identity
    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`).
    """
@ -408,12 +428,36 @@ 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)
-    # full volumes, fixed = CT, moving = T1c (no head box)
+    if use_head_box:
-    fx, mv = ct, t1c
+        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()
-    # No manual centroid init: antspyx's default initial_transform=None computes
+    # Centered start (REQUIRED for ANTs): antspyx cannot use a
-    # a translation aligning the image centers of mass.
+    # TranslationTransform as initial_transform (the ANTs CLI rejects it), and
    # 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,
@ -422,15 +466,20 @@ 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
-    # (T1c).  The contract is T1c -> CT, i.e. its inverse; fix_orientation
+    # (pre-shifted T1c).  We need the contract T1c(orig) -> CT: map T1c(orig) to
-    # verifies the direction against the skull ring.
+    # the pre-shifted T1c (translation by -off), then invert the ANTs transform.
    T_ants = sitk.ReadTransform(r["fwdtransforms"][0])
-    cand = T_ants.GetInverse()
+    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, full volumes)",
+        "metric": f"ANTs {aff_metric} (multi-resolution)",
-        "transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)",
+        "transform": f"ANTs {transform_type} (head box)", "convention": "T1c->CT (verified)",
        "auto_inverted": inverted, "aff_iterations": list(aff_iterations),
        "aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
    }
@ -691,18 +740,12 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
        "(ImpactFeaturesMapUpdateInterval -1)",
        '(ImpactWriteFeatureMaps "false")',
        "",
-        "// rigid 6-DOF to be optimized (auto center init, no -t0)",
+        "// rigid 6-DOF to be optimized (the shared centered start comes via -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")',
@ -745,12 +788,19 @@ 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=2000, samples=5000,
+def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
                    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
@ -759,26 +809,21 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
    the Mattes-MI pipelines above -- it is the cross-modality tool the image was
    built for.
-    The pipeline registers on the full CT / T1c volumes (no head box) with no
+    To stay on the benchmark's equal footing the pipeline:
-    manual centroid init: elastix's ``AutomaticTransformInitialization``
+      - registers on the same CT-space head box (both images resampled into it),
-    (CenterOfGravity) aligns the intensity centroids.  It uses a rigid 6-DOF
+      - starts from the *shared* centroid-aligned translation (passed to Elastix
-    ``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level
+        as a fixed initial transform via ``-t0``; with ``HowToCombineTransforms
-    scale pyramid (no image resampling).
+        Compose`` the final transform = ``T_main(T_init(p))``),
-
+      - uses a rigid 6-DOF ``EulerTransform`` (the 5.x name for the classic
-    NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other
+        3D rigid), ASGD, a 4-level scale pyramid (no image resampling).
-    pipelines).  The CT FOV (~300x300x251 mm) is much larger than the T1c crop
+    Budget: the default ``iters=500`` / ``samples=2000`` is the authors'
-    (~220x220x150 mm); with CT as the fixed image, too many CT samples map
+    recommended Jacobian budget.  It is *much* slower than the intensity-metric
-    outside the T1c buffer and the IMPACT metric aborts ("too many samples map
+    pipelines (minutes per case on GPU vs seconds) because every iteration runs
-    outside moving image buffer").  The returned transform still satisfies the
+    the TotalSegmentator feature extractor on the sampled patches with
-    unified contract (T1c->CT).
+    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
+    (an ``Euler3DTransform`` in the box/world physical frame) and composed with
-    :func:`fix_orientation`.
+    the ``-t0`` translation, then checked by :func:`fix_orientation`.
    Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, the
    ``vboussot/elastix_impact`` image, and the model asset under the cache.  A
@ -792,13 +837,25 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
        raise PipelineUnavailable(f"impact: {note}")
    model_host, model_rel = _impact_model_host(model)
-    # full volumes, fixed = T1c, moving = CT (no head box, no manual init --
+    lo, hi = _ct_head_box_mm(ct)
-    # elastix's AutomaticTransformInitialization aligns the image centers).
+    if lo is None:
-    # T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger
+        raise PipelineUnavailable("impact: no CT head box found (corrupt CT?)")
-    # than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples
+
-    # map outside the T1c buffer and the IMPACT metric aborts.
+    # resample both into the shared CT-space head box (mirrors the other pipelines)
-    fixed = t1c
+    tsp = np.array(t1c.GetSpacing())
-    moving = ct
+    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(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:
@ -842,8 +899,15 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
                        os.path.join(datadir, "Moving_image.mha"), True)
        with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
            f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
-                                       gpu, iters, samples, voxel_sizes,
+                                      gpu, iters, samples, voxel_sizes,
-                                       rescale=schedule, smoothing=smoothe))
+                                      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:
@ -855,6 +919,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
            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()
@ -869,12 +934,14 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
            tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
            raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
-        # the optimized Euler (world frame) is the full combination transform
+        # the optimized Euler (box/world frame) composed with the fixed -t0
-        # (AutomaticTransformInitialization is folded in by elastix); it maps
+        # translation: T = T_euler(T_init(p)).  Both map fixed(CT)->moving(T1c),
-        # fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT);
+        # i.e. the inverse of the contract (T1c->CT); fix_orientation resolves it.
        # fix_orientation resolves the direction.
        T_main = sitk.ReadTransform(itk_tf)
-        T, inverted = fix_orientation(ct, t1c, T_main)
+        T_init = sitk.TranslationTransform(3)
        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)
@ -883,7 +950,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
        "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, auto center init; "
+        "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, centroid init; "
                      f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
        "fixed": "T1c", "moving": "CT",
@ -896,128 +963,10 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
    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,
 }
@ -1045,13 +994,6 @@ 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: