6 changed files with 267 additions and 698 deletions
--- a/.gitignore
+++ b/.gitignore
@ -9,5 +9,3 @@ runs_monai/
 __pycache__/
 elastix/
 niftyreg/NiftyReg-Ubuntu-v2.1.1.zip
 niftyreg/NiftyReg-Ubuntu-v2.1.1/
--- a/REGISTRATION_BENCHMARK.md
+++ b/REGISTRATION_BENCHMARK.md
@ -22,50 +22,36 @@ 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**. Shared rules for
+destination→source convention this `T` maps **T1c → CT**. All pipelines
-every pipeline:
+register on the **full CT / T1c volumes** (no head box) and use each toolkit's
-
+**native automatic transform initialization** (no manual centroid offset):
- **No mask, no crop / pad**: registration runs on the full input volumes
+elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
-  exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
+antspyx its default center-of-mass alignment, simpleitk starts from identity,
-  padding).
+niftyreg (reg_aladin) initializes from the image centres, synthmorph needs
- **CT is always the fixed image** (moving = T1c); when both inputs are MR
+none. All use **fixed = CT, moving = T1c** except `impact`, which keeps
-  the first argument is kept as fixed.
+**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than
-
+the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
-All pipelines also use each toolkit's **native automatic transform
+samples map outside the T1c buffer and the IMPACT metric aborts ("too many
-initialization** (no manual centroid offset): elastix / impact use
+samples map outside moving image buffer"). The comparison is on each
-`AutomaticTransformInitialization` (CenterOfGravity), antspyx its default
+toolkit's native optimizer / metric / global-search capability.
 center-of-mass alignment, simpleitk starts from identity, niftyreg
 (reg_aladin) initializes from the image centres, synthmorph needs none. One
 consequence of the no-mask rule: the IMPACT metric only evaluates samples that
 map inside the moving buffer, so with the CT (head-neck, ~300×300×251 mm) as
 the fixed image and the tight T1c head crop (~220×220×150 mm) as moving, most
 of the CT buffer maps outside the T1c buffer (~14% of samples in-buffer at the
 auto-init). Elastix's sample-validity gate (build default
 `RequiredRatioOfValidSamples 0.25`) used to abort the M6 pairs on that ("too
 many samples map outside moving image buffer"); the parameter map now sets
 `RequiredRatioOfValidSamples 0.1` so the FOV-mismatched runs proceed — see the
 IMPACT section. The comparison is on each toolkit's native optimizer / metric
 / global-search capability.
 ## The six pipelines
-Listed alphabetically. All register on the **full volumes** (no head box, no
+Listed alphabetically. All register on the **full volumes** (no head box) and
-mask, no crop / pad) and use each toolkit's **native automatic transform
+use each toolkit's **native automatic transform initialization** (no manual
-initialization** (no manual centroid offset): elastix / impact use elastix's
+centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
-`AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default
+(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
-center-of-mass alignment, simpleitk starts from identity, niftyreg
+starts from identity, niftyreg initializes from the image centres, and
-initializes from the image centres, and synthmorph needs none (global deep
+synthmorph needs none (global deep method). All use **fixed = CT, moving =
-method). **All use fixed = CT, moving = T1c** (shared rule above; when both
+T1c** except `impact`, which keeps **fixed = T1c, moving = CT** (swapping
-inputs are MR the first argument is fixed). `niftyreg` and `synthmorph`
+breaks its TotalSegmentator metric — see below).
 clip the CT to [0, 80] HU before registration (per the CT-clipping A/B
 below; `clip_ct=None` disables it).
 | # | pipeline | fixed / moving | init | metric | transform | status in this env |
 |---|----------|----------------|------|--------|-----------|--------------------|
 | 1 | `antspyx`    | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
 | 2 | `elastix`    | CT / T1c | auto CenterOfGravity | Mattes MI (32 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
-| 3 | `impact`     | **CT / T1c** | auto CenterOfGravity | **IMPACT semantic** (MRSegmentator `MRSeg/MRSeg.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (Docker + GPU; sample gate 0.1 — see IMPACT section) |
+| 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 | `niftyreg`   | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (Docker `niftyreg:2.1.1-ubuntu`, GPU `-platf 1`) |
+| 4 | `niftyreg`   | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (main env, conda-forge `reg_aladin` 1.5.69) |
 | 5 | `simpleitk`  | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
 | 6 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
@ -76,10 +62,10 @@ level), with the sampler and metric left at their defaults
 `niftyreg` and `synthmorph` keep their native defaults.
 `simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid,
-identity start — the default, comparable to the other tools) and
+centered start — the default, comparable to the other tools) and
 `--simpleitk-mode production` (the full multi-start tilt grid from
-`02_build_m6_dataset.py`, now also on the full volumes per the no-crop rule,
+`02_build_m6_dataset.py`, ~20× slower, the most robust but not used for the
-~20× slower, the most robust but not used for the head-to-head).
+head-to-head).
 ## Test cases
@ -122,45 +108,40 @@ longer reported in the tables/plots):
  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 ⇔ `|mid| ≤ 20 mm` and `std ≤ 12 mm` —
+  production acceptance rule; aligned ⇔ `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`).
  a small negative `mid` is fine: on full-volume registration the CT skull
  ring can land slightly inside the MR silhouette, which is equally
  concentric).
-## Results (10 cases, no clip)
+## Results (10 cases)
 Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
-the MR head; higher = better) with the **unclipped** CT, so all rows are
+the MR head; higher = better).
 comparable on the same score scale. Full volumes, no mask / crop / pad,
 fixed = CT. `niftyreg` and `synthmorph` clip the CT to [0, 80] HU *before
 registration* (per the CT-clipping A/B below; `clip_ct=None` disables it);
 `antspyx`, `elastix` and `simpleitk` register the CT as-is. `niftyreg` is the
 v2.1.1 Docker build on the GPU (`-platf 1`); `impact` cannot run on these
 pairs under the no-mask rule (see the IMPACT section).
 | pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
 |---|---|---|---|---|---|
-| **antspyx** | ANTsPy 0.6.3 | 9 | 36.6 | 0.93 | **0.653 / 0.644** |
+| **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** |
-| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.517 / 0.520** |
+| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** |
-| **impact** | Docker Elastix 5.2 + IMPACT (MRSegmentator) | 8 | 2534 | 0.93 | **0.591 / 0.575** |
+| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** |
-| **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 20.8 | 0.93 | **0.292 / 0.338** |
+| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** |
-| **simpleitk** | SimpleITK 2.5.6 | 9 | 10.6 | 0.93 | **0.651 / 0.641** |
+| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** |
-| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 24.5 | 0.93 | **0.277 / 0.328** |
+| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** |
-Per-case MI (nats, unclipped score; **bold** = best for the case):
+\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default
 (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
 (hours/case) and was not re-run.
 Per-case MI (nats; **bold** = best for the case):
 | case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
 |---|---|---|---|---|---|---|
 | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
-| m6_22IRRKU6_20240115 | **0.649** | 0.472 | 0.635 | 0.292 | 0.642 | 0.277 |
+| m6_22IRRKU6_20240115 | 0.649 | 0.472 | **0.624** | 0.296 | 0.642 | 0.235 |
-| m6_23XIMO4M_20230222 | **0.655** | 0.416 | 0.478 | 0.632 | 0.652 | 0.616 |
+| m6_23XIMO4M_20230222 | 0.655 | 0.416 | 0.474 | **0.634** | 0.652 | 0.510 |
-| m6_24V62QQM_20230223 | **0.704** | 0.701 | 0.697 | 0.260 | 0.663 | 0.256 |
+| m6_24V62QQM_20230223 | 0.703 | **0.701** | 0.554 | 0.260 | 0.663 | 0.238 |
-| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.556 | 0.493 | **0.602** | 0.484 |
+| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.419 | 0.493 | **0.602** | 0.516 |
-| m6_24V62QQM_20240820 | **0.653** | 0.644 | 0.531 | 0.118 | 0.652 | 0.111 |
+| m6_24V62QQM_20240820 | 0.652 | 0.643 | 0.458 | 0.118 | **0.652** | 0.096 |
-| m6_26TYQI3C_20221027 | 0.659 | 0.653 | gate-abort | 0.609 | **0.695** | 0.590 |
+| m6_26TYQI3C_20221027 | 0.657 | 0.653 | 0.511 | 0.608 | **0.695** | 0.604 |
-| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | 0.450 | 0.069 | **0.577** | 0.061 |
+| m6_2BBN5XA4_20250411 | 0.576 | 0.287 | 0.155 | 0.069 | **0.577** | 0.060 |
-| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | 0.629 | 0.233 | 0.634 | 0.231 |
+| m6_2C4TKI6H_20240506 | 0.645 | 0.401 | **0.626** | 0.231 | 0.634 | 0.209 |
-| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | 0.625 | 0.340 | 0.651 | 0.323 |
+| m6_2CDNN7BH_20240102 | 0.655 | 0.516 | **0.618** | 0.344 | 0.651 | 0.301 |
 ### Interpretation
@ -182,138 +163,52 @@ Per-case MI (nats, unclipped score; **bold** = best for the case):
  single CenterOfGravity start still locks onto local minima on the
  repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
  0.416).
- **IMPACT (semantic metric) now runs on these pairs** — the no-mask abort was
+- **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI
-  elastix's sample-validity gate (build default `RequiredRatioOfValidSamples`
+  on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618),
-  0.25) rejecting the ~14% in-buffer fraction of the CT head-neck FOV vs the
+  where anatomical-segmentation similarity out-resolves intensity MI — but its
-  tight T1c crop; the parameter map sets it to 0.1 (see the IMPACT section).
+  median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is
-  With MRSegmentator at the 2000/5000 budget it is top-three in MI on 7 of the
+  also by far the slowest (~14 min/case at the shown budget; the new 2000/5000
-  8 registered cases (median 0.591): 0.635 on the hardest `m6_22IRRKU6`
+  budget would be on the order of hours/case).
  (masked-era 0.624), 0.697 on `m6_24V62QQM_20230223` (near-best), 0.629 on
  `m6_2C4TKI6H` — and it shares the production-acceptance of the grossly
  repositioned `m6_2BBN5XA4` (std 10.7 mm) with antspyx/simpleitk, the only
  non-intensity tool to do so. One run diverges mid-registration
  (`m6_26TYQI3C`: the in-buffer fraction falls below 10% at iteration 256 and
  the gate aborts — a stochastic-draw outcome). It remains by far the slowest
  (~42 min/case at 2000/5000).
 - **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
-  is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class;
+  is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class;
-  `m6_26TYQI3C` 0.609) but falls into the same tilt local minima
+  `m6_26TYQI3C` 0.608) but falls into the same tilt local minima
-  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). The v2.1.1 Docker build
+  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower
-  on the GPU (`-platf 1`) is ~21 s/case, ~4× slower than elastix (vs ~4 min
+  than elastix at full volumes.
  CPU for the conda 1.5.69 build). It clips the CT to [0,80] HU by default
  (per the A/B below) — neutral on the unclipped score scale
  (0.296→0.292 median).
 - **SynthMorph is weakest of the six on these cross-modality CT pairs**
-  (median 0.277 with its default CT clip; 0.241 unclipped): its `rigid` model
+  (median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30°
-  is trained on brain *MRI*, so the ~30° tilt defeats it on the repositioned
+  tilt defeats it on the repositioned cases (0.06–0.10) even though it matches
-  cases (0.06–0.11) even though it matches the intensity tools on the
+  the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no
-  cleanest (`m6_26TYQI3C` 0.590). The [0,80] brain-window clip of the CT is
+  centroid init** (global deep method) and is deterministic.
  its only genuine A/B win (+0.036 median, 7/9 cases), so it is the default.
  It needs **no centroid init** (global deep method).
 - **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the
  CT head everywhere, so coverage is not the discriminator — MI is.
 - **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s,
-  SimpleITK 10.6 s, **NiftyReg 20.8 s** (v2.1.1, GPU block matching),
+  SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s,
-  SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT is by far the
+  **NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous
-  slowest — 2534 s/case median at the 2000/5000 benchmark budget, ~276 s/case
+  budget; GPU feature extraction with backprop per iteration).
  at the authors' 500/2000; GPU feature extraction with backprop per
  iteration. No-mask runs are faster per iteration than the masked-era
  ~833 s/case @500/2000 because only the ~14% in-buffer samples are scored.)
 - **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;
  the in-process Elastix ASGD is deterministic given its fixed seed;
-  SynthMorph is deterministic; NiftyReg is deterministic on the CPU
+  SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**:
  (conda 1.5.69) but shows small run-to-run variation on the GPU (v2.1.1,
  ~±0.01–0.05 MI on the weak cases). **IMPACT is the most stochastic**:
  its Jacobian mode draws random patches every iteration, so the final pose
  changes run to run. Treat its per-case numbers as one sample from a
  distribution.
 ## CT clipping A/B (`--clip-ct 0 80`)
 Question: does clipping the CT to a brain window ([0, 80] HU — skull and air
 removed) **before registration** improve the result? The harness clips the CT
 before registration AND before scoring (percentile-binned MI is scale-free,
 so the A/B compares like for like). Same 10 cases, same budgets.
 Median achieved MI (nats) and wall-clock (9 complete-CT cases):
 | pipeline | no clip | clip [0,80] | Δ (median MI) | clip wins (cases by MI) |
 |---|---|---|---|---|
 | antspyx | 0.653 / 36.6 s | 0.466 / 34.7 s | **−0.186** | 0/9 |
 | elastix | 0.517 / 6.0 s | 0.425 / 5.9 s | **−0.092** | 1/9 |
 | niftyreg | 0.296 / 23.8 s | 0.320 / 20.2 s | +0.024 | 5/9 |
 | simpleitk | 0.651 / 10.6 s | 0.436 / 11.1 s | **−0.215** | 0/9 |
 | synthmorph | 0.241 / 26.0 s | 0.318 / 25.6 s | +0.077 | 5/9 |
 | impact | abort | abort | - | - (historical: at the time both aborted on the old 0.25 sample-validity gate; the gate is now 0.1, A/B to be re-done) |
 **Verdict: no — clipping the CT to [0, 80] HU does not improve the result; do
 not clip.**
 - The headline alignment score drops for **all three leading intensity
  pipelines**: antspyx 0.653→0.466, simpleitk 0.651→0.436, elastix
  0.517→0.425 (median nats; simpleitk and antspyx lose on **all 9** cases,
  elastix on 8/9).
 - The two weaker pipelines show small median-MI gains (niftyreg +0.024,
  synthmorph +0.077) but only mixed case-by-case (5/9), and their geometric
  QA does not support a real improvement.
 - **Geometric concentricity (skull-ring gap) is a wash**: mean mid/std move by
  ~1–5 mm in either direction (niftyreg mid 23.6→18.1 the only clear gain;
  antspyx std 16.4→18.0 the clear loss). Production acceptance even falls
  (no-clip 4/4/2/4/2 → clip 1/0/1/0/2).
 - **Mechanism.** Clipping is a saturating (non-invertible) intensity map: the
  skull atom (800+ HU → 80) removes the strongest intensity anchor
  (skull ↔ MR scalp) and compresses the CT rank distribution, so the
  percentile-binned MI between the warped CT and the MR is inherently lower
  *for the same pose*. The pose itself does not systematically degrade — the
  head is rigid, so skull-ring concentricity is its fingerprint, and that is
  unchanged. The MI drop is therefore mostly a score-scale effect, but it is
  exactly the score the benchmark reports, and there is no geometric benefit
  to pay it for.
 - Scoring note: `_skull_ring_gap` tests the ring with `>= p97` (not `>`) so
  the orientation gate and the gap QA also work on a clipped CT (where p97
  sits exactly on the clipped-skull atom).
 - **Decision / adopted default.** The pipelines were updated: `niftyreg` and
  `synthmorph` clip the CT to [0, 80] HU before registration
  (`clip_ct=(0.0, 80.0)` default, `clip_ct=None` disables it; MR inputs are
  never clipped); `antspyx` / `elastix` / `simpleitk` register the CT as-is.
  Re-checked on the **unclipped score scale** (the protocol used by the table
  above): synthmorph keeps a genuine gain (0.241→0.277 median, 7/9 cases),
  niftyreg is neutral (0.296→0.292; its apparent +0.024 above is mostly the
  clip-scoring effect). All rows of the main table are scored on the
  unclipped warped CT, so pipelines remain comparable.
 ## Reproducing
 ```bash
-# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1)
+# 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 niftyreg \
    --n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
 # force the conda 1.5.69 CPU build: --niftyreg-backend conda
 # Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
 # sections below for the one-time setup)
 python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \
    --workers 4 --threads 24
 # NOTE: impact now runs on the M6 pairs (the parameter map relaxes the
 # sample-validity gate to RequiredRatioOfValidSamples 0.1) but is by far the
 # slowest: ~5 min/case at 500/2000, ~42 min/case (2534 s median) at 2000/5000
 # (--impact-iters / --impact-samples) — see the IMPACT section
 python scripts/bench_registration.py --pipelines impact --n-cases 10 \
    --workers 4 --threads 24
 # A/B: does clipping the CT to a brain window help? (--clip-ct clips the CT
 # for ALL pipelines, before registration and scoring; niftyreg/synthmorph
 # already clip by default, so this only affects antspyx/elastix/simpleitk)
 python scripts/bench_registration.py --n-cases 10 --workers 4 --threads 24 \
    --out-dir results/bench_reg_noclip
 python scripts/bench_registration.py --n-cases 10 --workers 4 --threads 24 \
    --clip-ct 0 80 --out-dir results/bench_reg_clip0_80
 python scripts/bench_registration.py --merge --out-dir results/bench_reg_noclip
 python scripts/bench_registration.py --merge --out-dir results/bench_reg_clip0_80
 # combine the per-pipeline JSONs accumulated across envs into one table + plot
 python scripts/bench_registration.py --merge
 ```
@ -354,9 +249,11 @@ available pipelines with:
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 python scripts/bench_registration.py \
-    --pipelines simpleitk elastix antspyx niftyreg --n-cases 10 --qa-images
+    --pipelines simpleitk elastix antspyx --n-cases 10 --qa-images
 python scripts/bench_registration.py \
    --pipelines synthmorph --n-cases 10 --qa-images
 python scripts/bench_registration.py \
    --pipelines impact --n-cases 10 --qa-images
 ```
 (9 images each; the incomplete-CT case is skipped.)
@ -402,33 +299,20 @@ No extra model assets: the rigid weights ship inside the image
 ## 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): **symmetric block
+suite (KCL-BMEIS, the actively-maintained NiftyReg — the release zips for
-matching + trimmed least squares** (Modat et al., J. Med. Img. 1(2) 024003,
+v2.x need a newer glibc, so this env uses the conda-forge build):
-2014). Instead of a histogram MI objective it matches 4-voxel image blocks by
+**symmetric block matching + trimmed least squares** (Modat et al.,
-**normalised cross-correlation** at each level of a 3-level Gaussian pyramid
+J. Med. Img. 1(2) 024003, 2014). Instead of a histogram MI objective it
-and solves a robust (trimmed) least-squares for the rigid transform — a
+matches 4-voxel image blocks by **normalised cross-correlation** at each level
-genuinely different algorithmic family, and cross-modality tolerant because
+of a 3-level Gaussian pyramid and solves a robust (trimmed) least-squares for
-NCC is invariant to intensity scaling.
+the rigid transform — a genuinely different algorithmic family, and
-
+cross-modality tolerant because NCC is invariant to intensity scaling.
 **Default backend: the official v2.1.1 release in Docker.** The v2.x release
 binaries need GLIBC_2.34, which this host (Ubuntu 20.04, glibc 2.31) lacks, so
 the v2.1.1 Ubuntu zip is packaged in the `niftyreg:2.1.1-ubuntu` image
 (nvidia/cuda ubuntu22.04 base for the CUDA plugin). GPU runs use
 `docker run --gpus all` + `-platf 1 -gpuid N` and are ~20× faster than CPU;
 the v2.1.1 CUDA path supports **linear interpolation only** (`-interp 1`), so
 GPU runs use linear and CPU/conda runs use cubic (`-interp 3`) — the poses
 they recover on these cases agree to ~4 decimal places. The conda-forge
 build (1.5.69, CPU) remains the fallback when the Docker image is absent
 (`--niftyreg-backend conda`).
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
 ```bash
-# default: build the v2.1.1 Docker image (one-time)
+# one-time, in the benchmark env (provides reg_aladin + reg_f3d, NiftyReg 1.5.69)
 docker build -t niftyreg:2.1.1-ubuntu niftyreg/   # from NiftyReg-Ubuntu-v2.1.1.zip
 # fallback: conda-forge build (provides reg_aladin + reg_f3d, NiftyReg 1.5.69)
 conda install -c conda-forge niftyreg
 ```
@ -438,10 +322,8 @@ conda install -c conda-forge niftyreg
   (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 1` (GPU) / `-interp 3` (CPU) internal warping, default
+   coarse), `-interp 3` cubic internal warping, default initialisation
-   initialisation (image centres). GPU: `--gpus all`, `-platf 1`,
+   (image centres).
   `-gpuid <auto by pid>`; CPU/conda: `-platf 0` (v2.1.1) or no `-platf`
   (1.5.69).
 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
@ -454,17 +336,12 @@ conda install -c conda-forge niftyreg
 The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
 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 segmentation feature extractor (**MRSegmentator** `MRSeg/MRSeg.pt`
+TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
-by default; TotalSegmentator `TS/M730.pt` and other models from
+Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
-`VBoussot/impact-torchscript-models` are selectable). It is GPU-based and
+**fixed = T1c, moving = CT** (the reverse of the others): the CT FOV
-Docker-based (alongside `synthmorph`). It follows the shared rules (full
+(~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with
-volumes, no mask, fixed = CT, moving = T1c). The metric only evaluates samples
+CT as the fixed image too many CT samples map outside the T1c buffer and the
-that map inside the moving buffer, and the CT head-neck FOV (~300×300×251 mm)
+IMPACT metric aborts ("too many samples map outside moving image buffer").
 is much larger than the T1c head crop (~220×220×150 mm), so ~14% of the samples
 are in-buffer at the auto-init. The image build's sample-validity gate
 defaults to `RequiredRatioOfValidSamples 0.25` and used to abort the M6 pairs
 on that ("too many samples map outside moving image buffer"); the parameter
 map now sets it to `0.1` (see the IMPACT specifics below).
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
@ -477,40 +354,25 @@ docker --version
 docker pull vboussot/elastix_impact
 # 3. Fetch the TorchScript model (~123 MB, one-time) into the repo's model
-#    cache (data/impact_models/). MRSeg.pt (MRSegmentator) is the default:
+#    cache (data/impact_models/). M730.pt is the recommended TotalSegmentator.
-#    it is trained on both MRI and CT, the right prior for this
+mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/TS
-#    cross-modality CT->T1c task (all models share the same 2^l+3 FOV
+curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/TS/M730.pt \
-#    architecture, so the parameter map is unchanged between them).
+    https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/TS/M730.pt
 mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg
 curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg/MRSeg.pt \
    https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/MRSeg/MRSeg.pt
 ```
-Other models (TotalSegmentator `TS/M730.pt`, `TS/M258/M291/M850`, MIND, SAM2.1,
+Other models (M258/M291/M850, MIND, SAM2.1, …) are selectable with
-Anatomix, …) are selectable with `--impact-model <HF-relative-path-or-abs>`;
+`--impact-model <HF-relative-path-or-abs>`; see the
 see the
 [impact-torchscript-models](https://huggingface.co/VBoussot/impact-torchscript-models)
-repo for the full list and FOV/layer tables. The model choice does not affect
+repo for the full list and FOV/layer tables.
 the sample-validity gate (below): the "too many samples outside moving buffer"
 check is geometric (CT head-neck FOV vs the tight T1c crop) and fires at the
 first sample evaluation, before feature extraction.
 ### How the adapter drives the image
 The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
 binary directly (`/lib/elastix-install/bin/elastix`):
-1. **Full volumes, fixed=CT, moving=T1c** — both volumes are used as-is (no
+1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box).
-   head box, no mask, no crop / pad, per the shared rules above). The metric
+   T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too
-   only evaluates spatial samples that map inside the moving buffer, and the
+   many CT samples map outside the T1c buffer and the metric aborts.
   CT FOV (~300×300×251 mm, head+neck) is much larger than the tight T1c head
   crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
   the samples are in-buffer. The image build's default validity gate
   (`RequiredRatioOfValidSamples 0.25`) aborted on that ("Too many samples map
   outside moving image buffer: 704 / 5000"); the parameter map therefore sets
   `(RequiredRatioOfValidSamples 0.1)` so the M6 pairs run on their in-buffer
   samples while a genuinely diverging transform (<10% in-buffer) still
   aborts.
 2. **Auto init (no `-t0`)** — the parameter map sets
   `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
   elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
@ -519,24 +381,9 @@ binary directly (`/lib/elastix-install/bin/elastix`):
   **Jacobian** mode (differentiable, per the recommended config), 4 semantic
   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, `(RequiredRatioOfValidSamples 0.1)`. A copy of the
+   ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid
-   exact parameter map is dumped to
+   (`--impact-gpu -1` forces CPU).
-   `data/m6reg/impact/paramdump/ParameterMap_<model>_<pid>.txt` before every
+4. **Read-back** — the result `TransformParameters.0.itk.txt` is an
   docker run (also recorded as `param_dump` in the result JSON). GPU
   auto-spread across the host GPUs by pid (`--impact-gpu -1` forces CPU).
 4. **Intensity standardization happens inside the metric** (no pipeline-side
   normalization — the model card's "standardize intensities to zero mean and
   unit variance" is already applied): at init the metric computes
   whole-volume `[min, max, mean, sigma]` for each image
   (`itk::StatisticsImageFilter`) and passes it to the model, whose
   `normalize` branch rescales every extracted patch with the image's own
   `(mean, sigma)`. It is computed over the whole buffer (CT air included) and
   identically for fixed and moving, so both modalities' features sit on the
   same normalized scale. Pre-standardizing the MHA files in the pipeline
   would be an exact no-op (z-scoring is affine-invariant: the metric
   recomputes the stats of the image it receives) and would only corrupt the
 HU scale used by MI scoring and the QA brain-window images.
 5. **Read-back** — the result `TransformParameters.0.itk.txt` is an
    `Euler3DTransform` in the world physical frame (the
    `AutomaticTransformInitialization` is folded in by elastix); the adapter
    runs it through the usual `fix_orientation` gate.
@ -544,14 +391,13 @@ HU scale used by MI scoring and the QA brain-window images.
 ## Notes / caveats
 - The benchmark compares each toolkit in its **native configuration**: all
-  register on the full volumes (no head box, no mask, no crop / pad) and use
+  register on the full volumes (no head box) and use each toolkit's native
-  each toolkit's native automatic transform initialization (no manual
+  automatic transform initialization (no manual `centroid_offset`) — elastix /
-  `centroid_offset`) — elastix / impact `AutomaticTransformInitialization`
+  impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's
-  (CenterOfGravity), antspyx's center-of-mass alignment, simpleitk identity,
+  center-of-mass alignment, simpleitk identity, synthmorph none. Production
-  synthmorph none. Production uses `--simpleitk-mode production`'s
+  uses `--simpleitk-mode production`'s multi-start tilt grid, which raises
-  multi-start tilt grid, which raises SimpleITK's acceptance at ~20× the
+  SimpleITK's acceptance at ~20× the cost; that is out of scope for the
-  cost; that is out of scope for the head-to-head but is the config the
+  head-to-head but is the config the dataset actually ships with.
  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
@ -593,16 +439,11 @@ HU scale used by MI scoring and the QA brain-window images.
    - **`-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.
-    - **GPU is the default now** (v2.1.1, `-platf 1` via
+    - **CPU-only, thread-bound** (`-omp`); ~4 min/case median (120–440 s range)
-      `docker run --gpus all`): ~20–25 s/case on the full volumes vs ~4 min
+      at full volumes with 24 threads.
-      CPU (120–440 s, `-omp` thread-bound) for the conda 1.5.69 build.
+    - **Version note.** The upstream KCL-BMEIS release zips (v2.x) need
-      `--niftyreg-gpu -1` forces the CPU platform (`-platf 0`).
+      glibc ≥ 2.33; this host (Ubuntu 20.04, glibc 2.31) uses the conda-forge
-    - **Version note.** The upstream KCL-BMEIS release binaries (v2.x) need
+      build (1.5.69) instead.
      glibc ≥ 2.34; this host (Ubuntu 20.04, glibc 2.31) runs them in the
      `niftyreg:2.1.1-ubuntu` Docker image (nvidia/cuda ubuntu22.04 base;
      the CUDA plugin needs libcusolver.so.11, provided by the image). The
      conda-forge 1.5.69 build remains as the no-Docker fallback and gives
      essentially identical poses (MI matches to ±0.001 on the 10 cases).
 - **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
@ -629,41 +470,28 @@ HU scale used by MI scoring and the QA brain-window images.
    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 segmentation
+- **Slow, by design.** Every ASGD iteration runs the TotalSegmentator
-      extractor (MRSegmentator by default) on the sampled patches with
+      extractor on the sampled patches with backpropagation. The benchmark
      backpropagation. The benchmark
      budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
      the elastix default) is ~10× the authors' 500-iter / 2000-sample
-       Jacobian budget: measured 2534 s/case (median) on the no-mask M6
+      Jacobian budget, i.e. on the order of *hours* per case on full volumes —
-       pairs, ~276 s/case at 500/2000. `--impact-iters` / `--impact-samples`
+      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.
- **Runs on M6 with a relaxed sample-validity gate.** Like every other
+   - **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT
-      pipeline it runs fixed=CT/moving=T1c on the full volumes (no mask, no
+     follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251
-      crop / pad). The CT FOV (~300×300×251 mm) is much larger than the T1c
+     mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the
-      crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
+     fixed image too many CT samples map outside the T1c buffer and the metric
-      the samples map inside the T1c moving buffer. The metric calls
+     aborts ("too many samples map outside moving image buffer"). The returned
-      `CheckNumberOfSamples()` after every evaluation and the image build
+     transform still satisfies the contract (T1c→CT).
      defaults to `RequiredRatioOfValidSamples 0.25`, so it used to abort
      ("Too many samples map outside moving image buffer: 704 / 5000"). The
      parameter map now sets `(RequiredRatioOfValidSamples 0.1)`: the gate
      still aborts a genuinely diverging transform (<10% in-buffer) while the
      FOV-mismatched M6 pairs run on their ~14% in-buffer samples — exactly
      the native "only in-buffer samples are evaluated" behaviour. In the
      2000/5000 re-run the gate still fired once, mid-registration: on
      `m6_26TYQI3C` the transform drifted until only 468/5000 samples (9.4%)
      stayed in-buffer at iteration 256 and the run aborted. The transform
      contract (T1c→CT) holds; `fix_orientation` decides the direction as
      usual (all 8 completed runs report `auto_inverted: True`).
   - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns
     the intensity centroids; no manual `-t0` translation is used.
- **History: single-init convergence.** The semantic landscape has more
+   - **Rugged single-init convergence.** The semantic landscape has more local
-      local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer
+     minima than Mattes-MI: on full volumes it no longer nails the hardest
-      nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0;
+     repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box
-      the pre-mask head-box flow got it to mid 8.3 / std 13.2). The no-mask
+     used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration
-      MRSegmentator runs show the same pattern: strong MI (0.635 on
+     budget would likely raise its acceptance; that is out of scope for the
-      `m6_22IRRKU6` at 2000/5000) but the pose still misses the concentricity
+     single-init head-to-head.
      acceptance rule there (mid 26.1 / std 28.6); it *is* accepted on
      `m6_2BBN5XA4` (mid 15.1 / std 10.7).
   - **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
--- a/niftyreg/Dockerfile
+++ b/niftyreg/Dockerfile
@ -1,35 +0,0 @@
 # NiftyReg v2.1.1 (official KCL-BMEIS Ubuntu release) in Docker.
 #
 # The host (Ubuntu 20.04, glibc 2.31) is too old for the v2.x binaries, which
 # need GLIBC_2.34; the nvidia/cuda ubuntu22.04 base (glibc 2.35) runs them
 # directly.  The CUDA runtime provides libcusolver.so.11, which the
 # libniftyreg_cuda.so plugin (dlopen'd by reg_aladin with -platf 1) links
 # against; libcuda.so.1 comes from the host driver via `docker run --gpus all`.
 # Build:
 #   docker build -t niftyreg:2.1.1-ubuntu niftyreg/
 # Run (GPU, the pipeline's default):
 #   docker run --rm --gpus all -v /path/to/work:/work -w /work \
 #       niftyreg:2.1.1-ubuntu -ref ref.nii.gz -flo flo.nii.gz \
 #       -rigOnly -aff aff.txt -platf 1
 # Run (CPU fallback):
 #   docker run --rm -v /path/to/work:/work -w /work \
 #       niftyreg:2.1.1-ubuntu -ref ref.nii.gz -flo flo.nii.gz \
 #       -rigOnly -aff aff.txt -platf 0
 # reg_aladin is the ENTRYPOINT, so the NiftyReg CLI flags follow directly.
 # -platf 0|1 = CPU|CUDA, -gpuid selects the device, -omp sets OpenMP threads.
 FROM nvidia/cuda:12.4.1-runtime-ubuntu22.04
 RUN apt-get update && apt-get install -y --no-install-recommends \
        libgomp1 libstdc++6 zlib1g \
    && rm -rf /var/lib/apt/lists/*
 COPY NiftyReg-Ubuntu-v2.1.1/ /niftyreg/
 # reg_aladin searches /usr/local/lib for the CUDA plugin
 RUN chmod +x /niftyreg/* \
    && cp /niftyreg/libniftyreg_cuda.so /usr/local/lib/ \
    && ldconfig
 ENV PATH=/niftyreg:$PATH \
    LD_LIBRARY_PATH=/niftyreg:/usr/local/lib
 ENTRYPOINT ["reg_aladin"]
--- a/scripts/02_build_m6_dataset.py
+++ b/scripts/02_build_m6_dataset.py
@ -1,9 +1,8 @@
 """Build M6-2025 manifests: T1c volumes, labeled when GTV exists.
-GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on the full
+GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on a head box
-volumes -- no mask, no crop / pad; multi-start initial tilts because
+located in CT space; multi-start initial tilts because inter-session head tilt
-inter-session head tilt reaches ~30 deg and tilted MI local minima can
+reaches ~30 deg and tilted MI local minima can outscore the true alignment).
 outscore the true alignment).
 The inverse transform then warps CT and GTV (nearest) into T1c native space.
 Acceptance requires the warped CT skull ring to be concentric with the MR head
 silhouette (crop-independent check).
@ -106,10 +105,6 @@ def _skull_ring_gap(warped_ct, ha):
    gap between the MR head boundary and the CT skull ring (p97 of positive
    CT) at two levels. Returns {"mid_mm", "low_mm", "std_mm"} or None when no
    skull ring is found. Aligned = small gap with low angular variance.
    The ring test is ``>=`` (not ``>``) so it also works on a CT clipped to a
    brain window (e.g. [0, 80] HU), where the p97 threshold sits exactly on
    the clipped-skull atom and a strict comparison would give an empty ring.
    """
    posw = warped_ct[warped_ct > 0]
    if not posw.size:
@ -124,7 +119,7 @@ def _skull_ring_gap(warped_ct, ha):
        sl = sil[iz]
        if sl.sum() < 500:
            return None, None
-        sk = warped_ct[iz] >= bthr
+        sk = warped_ct[iz] > bthr
        if not sk.any():
            return None, None
        hy, hx = np.nonzero(sl)
@ -176,11 +171,14 @@ def _valid_pyramid(size, shrink, smooth, min_dim=6):
 def register_ct_to_t1c(ct, t1c, threads=12):
-    """Rigid registration with CT as fixed and T1c as moving, Mattes MI on the
+    """Rigid registration with CT as fixed and T1c as moving, Mattes MI on
-    full volumes (no mask, no crop / pad).
+    head-box images.
-    Both inputs share physical space (oriented). The transform is computed on
+    Both inputs share physical space (oriented). The registration grid is a head
-    the full CT / T1c volumes, so it applies directly to any resampling.
+    box located in CT space from the top of the head (the CTs are head-neck
    exams, so an MR-head crop would show mostly neck/air on the CT); the box
    excludes most air so MI is computed on tissue. The box is in physical
    space, so the resulting transform applies directly to the full volumes.
    Returns (transform, info): transform is the registration transform, mapping
    T1c (moving) coordinates -> CT (fixed) coordinates; resampling CT/TV onto
    the T1c grid with it applies the inverse transform to CT-space points,
@ -200,9 +198,38 @@ def register_ct_to_t1c(ct, t1c, threads=12):
    from src.common import head_mask_from_image
    ha = sitk.GetArrayFromImage(t1c)
    hm = head_mask_from_image(ha)
-    # full volumes, no mask / crop / pad: fixed = CT, moving = T1c
+    # Registration grid = a head box located in CT space. The CTs are head-neck
-    fixed, moving = ct, t1c
+    # (often head-neck-chest) exams, so cropping to the MR head box would leave
    # mostly neck/chest/air on the fixed image and MI locks onto the wrong
    # structures. Vertex = topmost tissue slice with 10mm of continuity below
    # (skips small artifacts floating above the head; bone HU is unreliable
    # here, so tissue > -150 HU is used). x/y from the tissue extent of the
    # 150mm band below the vertex. Box covers the whole head (x/y: extent +
    # 15mm margin, z: +10 to -215 from vertex).
    ca = sitk.GetArrayFromImage(ct)
    tissue = ca > -150
    csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
    vz = _ct_vertex_slice(ct)
    if vz is not 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)
    else:  # no head found: fall back to the whole CT volume
        gs = ct.GetSize()
        box_lo = tuple(float(x) for x in co)
        box_hi = (co[0] + gs[0] * csp[0], co[1] + gs[1] * csp[1], co[2] + gs[2] * csp[2])
    tsp = np.array(t1c.GetSpacing())
    shape = tuple(int(np.ceil((box_hi[d] - box_lo[d]) / tsp[d])) + 1 for d 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 box_lo))
    fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
    moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
    stages = []
@ -234,21 +261,8 @@ def register_ct_to_t1c(ct, t1c, threads=12):
    # rank candidates. Run the full pyramid from the centered start; if its
    # full-resolution geometric concentricity (skull-ring gap) is poor, repeat
    # from a grid of initial tilts and keep the best by gap.
    # Euler rotation center = the CT head centroid (the 150 mm band below the
    # vertex, same head definition as before); off maps it onto the MR head
    # centroid.  Both are init-only quantities; the registration itself runs on
    # the full volumes.
    ca = sitk.GetArrayFromImage(ct)
    csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
    tissue = ca > -150
    vz = _ct_vertex_slice(ct)
    if vz is not None:
        z0 = max(0, vz - int(round(150 / csp[2])))
        idx = np.array([z0, 0, 0]) + np.argwhere(tissue[z0:vz + 1]).mean(0)
    else:
        idx = np.argwhere(tissue).mean(0)
    ct_c = np.array(co) + idx * csp
    mr_c = np.array(t1c.GetOrigin()) + np.argwhere(hm).mean(0) * tsp
    ct_c = np.array([(box_lo[d] + box_hi[d]) / 2.0 for d in range(3)])
    off = (float(mr_c[0] - ct_c[0]), float(mr_c[1] - ct_c[1]), float(mr_c[2] - ct_c[2]))
    rad = np.pi / 180.0
@ -262,9 +276,7 @@ def register_ct_to_t1c(ct, t1c, threads=12):
        return tr0, gap
    def is_good(g):
-        # |mid| (not 0 <= mid): full-volume registration can land the CT skull
+        return bool(g) and g["std_mm"] <= 12.0 and 0.0 <= g["mid_mm"] <= 20.0
        # ring slightly inside the MR silhouette, which is equally concentric
        return bool(g) and g["std_mm"] <= 12.0 and abs(g["mid_mm"]) <= 20.0
    tr_reg, gap = full_run(0, 0, 0)
    if not is_good(gap):
@ -313,6 +325,7 @@ def register_ct_to_t1c(ct, t1c, threads=12):
                             "offset": [round(float(x), 3) for x in t0], "scale": scale},
            "coverage_head_box": round(coverage, 3),
            "skull_ring_gap_mm": skull_gap,
            "ct_head_box_mm": [[round(float(x), 1) for x in box_lo], [round(float(x), 1) for x in box_hi]],
            "seconds_total": round(time.time() - t_start, 1)}
    print("reg done %.0fs (rigid %.1fs) mi %.3f -> %.3f cov %.2f%s" % (
        info["seconds_total"], sum(s["seconds"] for s in stages),
@ -360,7 +373,7 @@ def reg_worker(row):
        tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1))
    tv_in_ct = tv_in_ct_frac >= 0.5  # slab must straddle the lesion
    gap = info["skull_ring_gap_mm"]
-    if gap and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
+    if gap and 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
        align = "concentric"
    else:
        align = "not_concentric"
--- a/scripts/bench_registration.py
+++ b/scripts/bench_registration.py
@ -20,10 +20,6 @@ higher = better in-sync = better aligned.  The production geometric QA
  - accepted           : production rule (nvox>0 and tv_in_ct and concentric).
  - reg_wall_s         : wall-clock registration time (excludes I/O + scoring).
 `--clip-ct LO HI` (e.g. `--clip-ct 0 80`, brain window) clips the CT to
 [LO, HI] before registration AND before scoring, so a clipped vs unclipped
 A/B compares the same input to both the toolkit and the MI score.
 Because some toolkits are not pip-installable into every
 `longitudinal-313` env (antspyx needs a cp310-cp313 wheel, see
 REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under
@ -97,14 +93,6 @@ def _m6():
    return RP._load_m6()
 def _clip_image(img, lo, hi):
    """Clip image intensities to [lo, hi] (geometry preserved)."""
    a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
    out = sitk.GetImageFromArray(a)
    out.CopyInformation(img)
    return out
 # --------------------------------------------------------------------------- #
 # case selection
 # --------------------------------------------------------------------------- #
@ -164,7 +152,7 @@ def evaluate(ct, t1c, tv, T):
    if nvox > 0:
        tvd = ndimage.binary_dilation(warped_tv > 0.5, iterations=3)
        tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1))
-    ok_gap = bool(gap) and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0
+    ok_gap = bool(gap) and 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0
    accepted = nvox > 0 and tv_in_ct_frac >= 0.5 and ok_gap
    return {
        "skull_ring_gap_mm": gap,
@ -249,9 +237,7 @@ def _call(pipeline, ct, t1c, tv, args):
        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,
+        return RP.register_niftyreg(ct, t1c, threads=args.threads)
                                    backend=args.niftyreg_backend,
                                    gpu=args.niftyreg_gpu)
    raise KeyError(pipeline)
@ -269,11 +255,6 @@ def _run_case(pipeline, case, args):
                    "alignment": "ct_incomplete", "accepted": False,
                    "note": "CT head incomplete (data issue, not a pipeline failure)"})
        return out
    if args.clip_ct is not None:
        # clip BEFORE registration AND scoring so the A/B compares like for like
        lo, hi = args.clip_ct
        ct = _clip_image(ct, lo, hi)
        out["clip_ct"] = [lo, hi]
    t0 = time.time()
    T, meta = _call(pipeline, ct, t1c, tv, args)
    out["reg_wall_s"] = round(time.time() - t0, 1)
@ -399,8 +380,8 @@ def write_plot(rows, out_png):
        print("no available pipelines to plot")
        return
    labels = list(ok_rows)
-    regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") or 0 for p in labels]
+    regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels]
-    mi = [(ok_rows[p].get("mi_ct_mr") or {}).get("median") 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))
    colors = ["#4878CF", "#6ACC65"]
    for ax, (title, vals, fmt, ylab) in zip(axes, [
@ -506,22 +487,6 @@ def parse_args():
    ap.add_argument("--impact-samples", type=int, default=5000,
                    help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
                         "elastix default)")
    ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"],
                    default="auto",
                    help="niftyreg reg_aladin build: docker = official v2.1.1 in the "
                         "niftyreg:2.1.1-ubuntu image (host glibc is too old for the "
                         "v2.x binaries), conda = conda-forge 1.5.69 on PATH, "
                         "auto (default) = docker when the image is present else conda")
    ap.add_argument("--niftyreg-gpu", type=int, default=None,
                    help="niftyreg CUDA device index (docker backend; default: "
                         "auto-spread over the host GPUs by pid, -1 forces CPU "
                         "-platf 0; on the GPU the run uses 'docker run --gpus all "
                         "-platf 1' with linear interpolation, the v2.1.1 CUDA "
                         "limit)")
    ap.add_argument("--clip-ct", nargs=2, type=float, default=None, metavar=("LO", "HI"),
                    help="clip the CT (HU) to [LO, HI] BEFORE registration and scoring "
                         "(e.g. --clip-ct 0 80 = brain window, skull/air removed; "
                         "default: no clipping)")
    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
@ -23,26 +23,15 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
 Pipelines (listed alphabetically)
 ---------------------------------
-Shared rules for every pipeline:
+All pipelines register on the **full volumes** (no head box) and use each
-
+toolkit's **native automatic transform initialization** (no manual
- **No mask, no crop / pad**: registration always runs on the full input
+``centroid_offset``): elastix / impact use elastix's
  volumes exactly as given (no ``FixedMask`` / ``MovingMask``, no head-box
  resampling, no padding).
 - **CT is always the fixed image** (moving = T1c).  When both inputs are MR
  (no CT present) the first argument is kept as fixed (see
  :func:`fixed_moving`).
 - **CT brain-window clip for the pipelines that benefit from it**: per the
  CT-clipping A/B (see REGISTRATION_BENCHMARK.md), ``niftyreg`` and
  ``synthmorph`` clip the CT to [0, 80] HU before registration
  (``clip_ct`` parameter, default on, MR inputs never clipped); the
  Mattes-MI pipelines (antspyx / elastix / simpleitk) register the CT as-is.
 Each pipeline uses its toolkit's **native automatic transform initialization**
 (no manual ``centroid_offset``): elastix / impact use elastix's
 ``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its
 default center-of-mass alignment, simpleitk starts from identity, niftyreg
 (reg_aladin) initializes from the image centres, and synthmorph needs none
-(global deep method).
+(global deep method).  All use **fixed = CT, moving = T1c** except IMPACT,
 which keeps **fixed = T1c, moving = CT** (swapping breaks its
 TotalSegmentator metric -- see ``register_impact``).
 - ``register_antspyx``    : ANTs / ANTsPy ``ants.registration`` (Rigid).
  Default ``initial_transform=None`` computes a translation aligning the image
@ -53,9 +42,9 @@ default center-of-mass alignment, simpleitk starts from identity, niftyreg
 - ``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).  full volumes, fixed=CT/moving=T1c (no mask -- shared
+  GPU (Jacobian mode).  Kept as fixed=T1c/moving=CT (see above); elastix
-  rules above); elastix ``AutomaticTransformInitialization`` (CenterOfGravity)
+  ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
-  provides the start.  Needs Docker + the image + a model asset.
+  Needs Docker + the image + a model asset.
 - ``register_niftyreg``   : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
  symmetric block-matching + trimmed least squares (Modat et al., J. Med.
  Img. 2014), NCC-based block similarity (cross-modality tolerant).  Rigid
@ -134,50 +123,6 @@ def _warp_ct(ct, t1c, T):
    return sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
 def _is_ct(img):
    """Modality heuristic: CT is in Hounsfield units and carries air at large
    negative values; MR intensities are non-negative."""
    return float(sitk.GetArrayFromImage(img).min()) < -50.0
 def _modality_name(img):
    return "CT" if _is_ct(img) else "MR"
 def fixed_moving(img_a, img_b):
    """Shared fixed-image rule: the CT is always the fixed image.  When both
    inputs are MR (or both CT) the given order is kept, i.e. ``img_a`` is the
    fixed image.  Returns ``(fixed, moving)``.
    Note: a CT that has been intensity-clipped to a non-negative window (e.g.
    [0, 80] HU) is no longer detectable as CT, so the given order is kept --
    callers should pass the CT first (as all pipelines do)."""
    if _is_ct(img_b) and not _is_ct(img_a):
        return img_b, img_a
    return img_a, img_b
 def _clip_intensity(img, lo, hi):
    """Clip image intensities to [lo, hi] (geometry preserved)."""
    a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
    out = sitk.GetImageFromArray(a)
    out.CopyInformation(img)
    return out
 def _maybe_clip_ct(fixed, clip_ct):
    """Apply the pipeline's CT brain-window clip to the fixed image.
    The A/B benchmark (``--clip-ct 0 80``) showed clipping the CT to a brain
    window helps the block-matching / deep pipelines (niftyreg, synthmorph)
    and hurts the Mattes-MI ones (antspyx, elastix, simpleitk), so those two
    pipelines clip the CT by default.  Returns ``(fixed, clipped: bool)``;
    MR fixed images (both-MR pairs) are never clipped."""
    if clip_ct and _is_ct(fixed):
        return _clip_intensity(fixed, *clip_ct), True
    return fixed, False
 def fix_orientation(ct, t1c, cand):
    """Return (best_T, inverted) choosing between ``cand`` and ``cand^-1``.
@ -258,9 +203,8 @@ 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, no mask / crop / pad; CT is the fixed image (or, if both
+    # full volumes, fixed = CT, moving = T1c (no head box)
-    # inputs are MR, the first argument is)
+    fixed, moving = ct, t1c
    fixed, moving = fixed_moving(ct, t1c)
    shape = np.array(ct.GetSize())
    # Identity start (no centroid init); the multi-resolution Mattes-MI pyramid
    # recovers the gross inter-session offset.  The Euler center is the fixed
@ -294,7 +238,6 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
        "transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
                      f"{iters_full} iters/level, identity start)"),
        "convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
        "fixed": _modality_name(fixed), "moving": _modality_name(moving),
        "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
        "skull_ring_gap_mm": gap,
    }
@ -372,8 +315,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
    after a successful registration.
    """
    itk = _elastix_api()
-    # full volumes, no mask / crop / pad; CT is the fixed image
+    # full volumes, fixed = CT, moving = T1c (no head box)
-    fixed, moving = fixed_moving(ct, t1c)
+    fixed = ct
    moving = t1c
    params = {
        "ImageDimension": 3,
        "FixedInternalImagePixelType": "float",
@ -429,8 +373,7 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
        "metric": "MattesMutualInformation (full volumes)",
        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)",
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
-        "fixed": _modality_name(fixed), "moving": _modality_name(moving),
+        "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
        "seconds_reg": round(dt, 1),
    }
    return T, meta
@ -465,8 +408,8 @@ 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, no mask / crop / pad; CT is the fixed image
+    # full volumes, fixed = CT, moving = T1c (no head box)
-    fx, mv = fixed_moving(ct, t1c)
+    fx, mv = ct, t1c
    import time
    t0 = time.time()
    # No manual centroid init: antspyx's default initial_transform=None computes
@ -488,8 +431,7 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
        "pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
        "metric": f"ANTs {aff_metric} (multi-resolution, full volumes)",
        "transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)",
-        "auto_inverted": inverted, "fixed": _modality_name(fx), "moving": _modality_name(mv),
+        "auto_inverted": inverted, "aff_iterations": list(aff_iterations),
        "aff_iterations": list(aff_iterations),
        "aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
    }
    return T, meta
@ -542,8 +484,7 @@ def _lta_to_sitk(path):
 def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
-                        extent=256, keep_scratch=False,
+                        extent=256, keep_scratch=False):
                        clip_ct=(0.0, 80.0)):
    """SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image
    ``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``,
    pretrained ``synthmorph.rigid.1.h5`` weights).
@ -551,10 +492,7 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
    SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware,
    acquisition-agnostic deep registration: each volume is normalized to an
    isotropic 1-mm LIA network space and the (pretrained) network estimates the
-    transform, so no skull-stripping or intensity preprocessing is needed --
+    transform, so no skull-stripping or intensity preprocessing is needed.  Being
    except the one it benefits from: the CT is clipped to a brain window
    (``clip_ct=(0.0, 80.0)`` by default, per the CT-clipping A/B) before
    registration; pass ``clip_ct=None`` to register the CT as-is.  Being
    a global deep method it recovers the gross inter-session offset on its own,
    so **no centroid init is used**.  The pipeline
      - runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible
@ -586,15 +524,11 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
    # per-case scratch on a local (non-NFS) tmp fs
    tmp = tempfile.mkdtemp(prefix="synthmorph_reg_")
    try:
-        # full volumes, no mask / crop / pad; CT is the fixed image
+        # fixed = CT, moving = T1c (shared convention for all pipelines)
-        fixed, moving = fixed_moving(ct, t1c)
+        sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
-        fixed_name = _modality_name(fixed)  # before clipping (a clipped CT
+                        os.path.join(tmp, "fixed.nii.gz"), True)      # CT
-        # is no longer intensity-detectable)
+        sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
-        fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
+                        os.path.join(tmp, "moving.nii.gz"), True)     # T1c
        sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
                        os.path.join(tmp, "fixed.nii.gz"), True)
        sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
                        os.path.join(tmp, "moving.nii.gz"), True)
        # run as the host uid/gid so the container-written LTA is host-owned
        # (clean scratch removal) and TF still reaches the GPU.
@ -636,11 +570,10 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
                      f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, "
                      f"extent {extent})"),
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
-        "fixed": fixed_name, "moving": _modality_name(moving),
+        "fixed": "CT", "moving": "T1c",
        "seconds_reg": round(dt, 1),
        "docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu,
        "extent": extent,
        "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
    }
    if keep_scratch:
        meta["scratch"] = tmp
@ -654,7 +587,7 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
 IMPACT_IMAGE = "vboussot/elastix_impact"
 IMPACT_ELASTIX_BIN = "/lib/elastix-install/bin/elastix"
 IMPACT_MODEL_CACHE = "data/impact_models"
-IMPACT_DEFAULT_MODEL = "MRSeg/MRSeg.pt"
+IMPACT_DEFAULT_MODEL = "TS/M730.pt"
 IMPACT_HF_MODEL_URL = (
    "https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/{rel}"
 )
@ -757,14 +690,6 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
        f'(ImpactUseMixedPrecision "{mixed}")',
        "(ImpactFeaturesMapUpdateInterval -1)",
        '(ImpactWriteFeatureMaps "false")',
        # safety gate: abort if < ratio of the spatial samples map inside the
        # moving buffer.  The build default is 0.25, but the no-mask M6 pairs
        # (CT head-neck FOV vs the tight T1c crop) only have ~14% in-buffer at
        # the CenterOfGravity auto-init, so elastix aborts at the first metric
        # evaluation.  0.1 keeps the gate (a diverging transform still aborts)
        # while letting the FOV-mismatched runs proceed -- the metric natively
        # evaluates only the in-buffer samples.
        "(RequiredRatioOfValidSamples 0.1)",
        "",
        "// rigid 6-DOF to be optimized (auto center init, no -t0)",
        '(Transform "EulerTransform")',
@ -783,8 +708,7 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
        '(FixedInternalImagePixelType "float")',
        '(MovingInternalImagePixelType "float")',
        '(UseDirectionCosines "true")',
-    ]
+        "",
    L += [
        '(Registration "MultiMetricMultiResolutionRegistration")',
        f"(NumberOfResolutions {n})",
        f"(MaximumNumberOfIterations {iters})",
@ -824,25 +748,16 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
 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=CT, moving=T1c) in the Docker image
+    """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image
    ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
    IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
-    one: a TorchScript segmentation model (``MRSeg/MRSeg.pt`` by default,
+    one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run
-    MRSegmentator -- robust multi-modality segmentation on MRI *and* CT, the
+    on randomly sampled patches of the fixed and moving images and the
    right prior for this cross-modality CT->T1c task; TotalSegmentator
    ``TS/M730.pt`` and the other models in the HF repo
    ``VBoussot/impact-torchscript-models`` work too via ``model=``) is run on
    randomly sampled patches of the fixed and moving images and the
    registration minimizes the distance between the extracted feature maps
-    (L2, per the recommended config).  Intensity standardization (zero mean /
+    (L2, per the recommended config).  This is a genuinely different metric from
-    unit variance per the model card) is applied *inside the metric*: at init
+    the Mattes-MI pipelines above -- it is the cross-modality tool the image was
-    it computes whole-volume [min, max, mean, sigma] per image and the model's
+    built for.
    normalize branch rescales every patch with that image's own (mean, sigma),
    so no pipeline-side pre-normalization is needed (and would be an exact
    no-op -- z-scoring is affine-invariant).  This is a genuinely different
    metric from 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
    manual centroid init: elastix's ``AutomaticTransformInitialization``
@ -850,17 +765,12 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
    ``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level
    scale pyramid (no image resampling).
-    Like every other pipeline this runs with **fixed=CT, moving=T1c** on the
+    NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other
-    full volumes: no mask, no crop / pad (shared rules above).  The IMPACT
+    pipelines).  The CT FOV (~300x300x251 mm) is much larger than the T1c crop
-    metric draws random spatial samples from the fixed (CT) buffer and
+    (~220x220x150 mm); with CT as the fixed image, too many CT samples map
-    evaluates them only where they map inside the moving (T1c) buffer; when
+    outside the T1c buffer and the IMPACT metric aborts ("too many samples map
-    too few samples are in-buffer, elastix aborts ("too many samples map
+    outside moving image buffer").  The returned transform still satisfies the
-    outside moving image buffer").  On the M6 head-neck CTs (FOV
+    unified contract (T1c->CT).
    ~300x300x251 mm) vs the tight T1c head crop (~220x220x150 mm) only ~14%
    of the samples are in-buffer at the CenterOfGravity auto-init, so the
    M6 runs abort -- IMPACT needs a fixed volume whose FOV is comparable to
    the moving crop (e.g. a head-only CT).  The returned transform satisfies
    the 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'
@ -882,10 +792,13 @@ 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, no mask / crop / pad; CT is the fixed image (or, if both
+    # full volumes, fixed = T1c, moving = CT (no head box, no manual init --
-    # inputs are MR, the first argument).  No manual init -- elastix's
+    # elastix's AutomaticTransformInitialization aligns the image centers).
-    # AutomaticTransformInitialization aligns the image centers.
+    # T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger
-    fixed, moving = fixed_moving(ct, t1c)
+    # than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples
    # map outside the T1c buffer and the IMPACT metric aborts.
    fixed = t1c
    moving = ct
    # GPU: auto-spread across the host GPUs by pid, or forced value
    if gpu is None:
@ -893,8 +806,8 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
        gpu = (os.getpid() % ngpu) if ngpu else -1
    # per-level feature-extraction scales (recommended 6/3/1.5/1 mm), scaled so
-    # the finest level matches the moving image's base spacing
+    # the finest level matches the box's (T1c) base spacing
-    base = float(np.median(moving.GetSpacing()))
+    base = float(np.median(t1c.GetSpacing()))
    n = len(scales)
    voxel_sizes = [base * s for s in scales]
    # image pyramid: IMPACT's docs recommend NO image resampling (all-1 schedules)
@ -927,22 +840,10 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
                        os.path.join(datadir, "Fixed_image.mha"), True)
        sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
                        os.path.join(datadir, "Moving_image.mha"), True)
        param_text = _impact_param_map(os.path.join("/Data/Models", model_rel),
                                       gpu, iters, samples, voxel_sizes,
                                       rescale=schedule, smoothing=smoothe)
        with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
-            f.write(param_text)
+            f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
-
+                                       gpu, iters, samples, voxel_sizes,
-        # debug aid: keep a persistent copy of the exact parameter map elastix
+                                       rescale=schedule, smoothing=smoothe))
        # is about to run, written BEFORE the docker call so an abort leaves it
        # behind for inspection (per-pid name: parallel workers do not clobber
        # each other's dumps)
        dump_dir = d("data/m6reg/impact/paramdump")
        os.makedirs(dump_dir, exist_ok=True)
        family = os.path.basename(os.path.dirname(model_rel)) or "model"
        param_dump = os.path.join(dump_dir, f"ParameterMap_{family}_{os.getpid()}.txt")
        with open(param_dump, "w") as f:
            f.write(param_text)
        cmd = ["docker", "run", "--rm"]
        if gpu >= 0:
@ -969,28 +870,26 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
            raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
        # the optimized Euler (world frame) is the full combination transform
-        # (AutomaticTransformInitialization is folded in by elastix); elastix
+        # (AutomaticTransformInitialization is folded in by elastix); it maps
-        # maps moving(T1c) -> fixed(CT), which is the contract direction;
+        # fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT);
-        # fix_orientation verifies it against the skull ring.
+        # fix_orientation resolves the direction.
        T_main = sitk.ReadTransform(itk_tf)
        T, inverted = fix_orientation(ct, t1c, T_main)
    finally:
        if not keep_scratch:
            shutil.rmtree(tmp, ignore_errors=True)
    family = os.path.basename(os.path.dirname(model_rel)) \
        or os.path.splitext(os.path.basename(model_rel))[0]
    meta = {
        "pipeline": "impact",
        "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
-        "metric": f"IMPACT semantic ({family} {os.path.basename(model_rel)}, L2, Jacobian)",
+        "metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)",
        "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": _modality_name(fixed), "moving": _modality_name(moving),
+        "fixed": "T1c", "moving": "CT",
        "seconds_reg": round(dt, 1),
        "docker_image": IMPACT_IMAGE, "model": model_rel, "gpu": gpu,
-        "iters": iters, "samples": samples, "param_dump": param_dump,
+        "iters": iters, "samples": samples,
    }
    if keep_scratch:
        meta["scratch"] = tmp
@ -1002,9 +901,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
 # --------------------------------------------------------------------------- #
 NIFTYREG_DOCKER_IMAGE = "niftyreg:2.1.1-ubuntu"
 def _niftyreg_bin():
    """Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable."""
    b = shutil.which("reg_aladin")
@ -1017,26 +913,6 @@ def _niftyreg_bin():
        "re-run in the activated env.")
 def _niftyreg_docker_ok():
    """Return (ok, note): is the NiftyReg v2.1.1 Docker image available?
    The official KCL-BMEIS v2.x release binaries need GLIBC_2.34, which this
    host (Ubuntu 20.04, glibc 2.31) lacks; the ``niftyreg:2.1.1-ubuntu`` image
    (built from the release zip, see niftyreg/Dockerfile) runs them instead.
    """
    if shutil.which("docker") is None:
        return False, "docker CLI not found on PATH"
    try:
        r = subprocess.run(["docker", "image", "inspect", NIFTYREG_DOCKER_IMAGE],
                           capture_output=True, text=True, timeout=30)
        if r.returncode != 0:
            return False, (f"docker image {NIFTYREG_DOCKER_IMAGE} not built "
                           f"(run: docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/)")
    except Exception as e:  # noqa
        return False, f"docker image check failed ({type(e).__name__}: {e})"
    return True, NIFTYREG_DOCKER_IMAGE
 def _niftyreg_version(binpath):
    try:
        r = subprocess.run([binpath, "--version"], capture_output=True,
@ -1062,8 +938,7 @@ def _parse_niftyreg_affine(path):
    return T
-def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False,
+def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False):
                      backend="auto", gpu=None, clip_ct=(0.0, 80.0)):
    """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
@ -1076,96 +951,32 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
    transforms concurrently; ``-maxit 5`` (the tool default) is too coarse for
    these repositioned pairs -- the benchmark uses 15.
-    The CT is clipped to a brain window (``clip_ct=(0.0, 80.0)`` by default,
+    Full volumes (no head box); default initialisation (image centres);
-    per the CT-clipping A/B) before registration; pass ``clip_ct=None`` to
+    ``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3`` (cubic) internal
-    register the CT as-is.  MR inputs are never clipped.
+    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.
-    Full volumes (no head box, no mask / crop / pad); default initialisation
+    Preconditions (else :class:`PipelineUnavailable`): ``reg_aladin`` on PATH
-    (image centres); ``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3``
+    (``conda install -c conda-forge niftyreg``).
    (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.
    ``backend`` selects the reg_aladin build: ``"docker"`` runs the official
    KCL-BMEIS **v2.1.1** release in the ``niftyreg:2.1.1-ubuntu`` Docker image
    (built from NiftyReg-Ubuntu-v2.1.1.zip, see niftyreg/Dockerfile -- the
    v2.x binaries need GLIBC_2.34, which this host lacks); ``"conda"`` uses
    the conda-forge ``reg_aladin`` on PATH (NiftyReg 1.5.69); ``"auto"``
    (default) prefers Docker when the image is present and falls back to
    conda.
    ``gpu`` selects the CUDA device for the Docker backend: ``None`` (default)
    auto-spreads over the host GPUs by pid, ``-1`` forces CPU.  On the GPU the
    run uses ``docker run --gpus all`` + ``-platf 1 -gpuid N`` and
    ``-interp 1`` (the v2.1.1 CUDA path supports linear interpolation only);
    on the CPU it uses ``-platf 0`` with ``-interp 3`` (cubic), as does the
    conda backend.  The GPU is ~20x faster than CPU on the full volumes.
    Preconditions (else :class:`PipelineUnavailable`): the Docker image for
    backend docker/auto, or ``reg_aladin`` on PATH (``conda install -c
    conda-forge niftyreg``).
    """
    import time
    if backend not in ("auto", "conda", "docker"):
        raise PipelineUnavailable(
            f"niftyreg: unknown backend {backend!r} (use auto|conda|docker)")
    dok, dnote = _niftyreg_docker_ok()
    use_docker = backend == "docker" or (backend == "auto" and dok)
    if use_docker:
        if not dok:
            raise PipelineUnavailable(f"niftyreg: {dnote}")
        binpath = None
        if gpu is None:
            ngpu = _n_gpus()
            gpu = (os.getpid() % ngpu) if ngpu else -1
        use_gpu = gpu >= 0
        platf = "1" if use_gpu else "0"
        # the v2.1.1 CUDA path supports linear interpolation only
        interp = "1" if use_gpu else "3"
        toolkit = (f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE}, "
                   f"CUDA)" if use_gpu
                   else f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE})")
    else:
    binpath = _niftyreg_bin()
-        gpu, use_gpu, platf = -1, False, None
+    # full volumes, reference = CT, floating = T1c (no head box)
        interp = "3"
        toolkit = _niftyreg_version(binpath)
    # full volumes, no mask / crop / pad; reference (fixed) = CT, floating = T1c
    fixed, moving = fixed_moving(ct, t1c)
    fixed_name = _modality_name(fixed)  # before clipping (a clipped CT is no
    # longer intensity-detectable)
    fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
    tmp = tempfile.mkdtemp(prefix="niftyreg_reg_")
    dt = 0.0
    try:
-        sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
+        sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
                        os.path.join(tmp, "ref.nii.gz"), True)
-        sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
+        sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
                        os.path.join(tmp, "flo.nii.gz"), True)
-        flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
+        cmd = [binpath, "-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
               "-rigOnly", "-aff", "aff.txt",
               "-ln", str(levels), "-maxit", str(maxit),
-                 "-interp", interp, "-omp", str(threads)]
+               "-interp", "3", "-omp", str(threads), "-voff"]
        if platf is not None:
            flags += ["-platf", platf]
            if use_gpu:
                flags += ["-gpuid", str(gpu)]
        flags += ["-voff"]
        if use_docker:
            # reg_aladin is the image ENTRYPOINT; the work dir is /work.
            # --gpus all: the v2.1.1 Ubuntu build ships libniftyreg_cuda.so
            # (dlopen'd by reg_aladin with -platf 1); the NVIDIA container
            # toolkit provides libcuda.so.1, and -gpuid selects the device.
            cmd = ["docker", "run", "--rm"]
            if use_gpu:
                cmd += ["--gpus", "all"]
            cmd += ["-v", f"{tmp}:/work", "-w", "/work",
                    NIFTYREG_DOCKER_IMAGE] + flags
        else:
            cmd = [binpath] + flags
        t0 = time.time()
        proc = subprocess.run(cmd, capture_output=True, text=True,
                              cwd=tmp, timeout=7200)
@ -1188,16 +999,14 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
    meta = {
        "pipeline": "niftyreg",
-        "toolkit": toolkit,
+        "toolkit": _niftyreg_version(binpath),
        "backend": "docker" if use_docker else "conda",
        "metric": "NCC block matching + trimmed least squares (symmetric, full volumes)",
        "transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, "
                      f"{maxit} TLS iters/level, cubic interp)"),
        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
-        "fixed": fixed_name, "moving": _modality_name(moving),
+        "fixed": "CT", "moving": "T1c",
        "seconds_reg": round(dt, 1),
-        "levels": levels, "maxit": maxit, "interp": interp, "gpu": gpu,
+        "levels": levels, "maxit": maxit,
        "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
    }
    if keep_scratch:
        meta["scratch"] = tmp
@ -1238,20 +1047,11 @@ def availability():
                    out[name] = (False, f"ants(antspyx) not importable ({type(e).__name__})")
            elif name == "niftyreg":
                b = shutil.which("reg_aladin")
                dok, dnote = _niftyreg_docker_ok()
                if dok:
                    ver = f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE})"
                if b:
-                        ver += f" | conda fallback: {_niftyreg_version(b)}"
+                    out[name] = (True, _niftyreg_version(b))
                    out[name] = (True, ver)
                elif b:
                    out[name] = (True, _niftyreg_version(b)
                                 + " | docker v2.1.1 image missing (build: "
                                 f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/)")
                else:
-                    out[name] = (False, "reg_aladin not on PATH and no docker image; "
+                    out[name] = (False, "reg_aladin not on PATH "
-                                        "conda install -c conda-forge niftyreg or "
+                                        "(conda install -c conda-forge niftyreg)")
                                        f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/")
            elif name == "synthmorph":
                ok, note = _synthmorph_docker_ok()
                if not ok: