fix(benchmark): lower IMPACT sample-validity gate to unblock M6 runs

The no-mask CT→T1c FOV mismatch left only ~14% of spatial samples
inside the moving buffer at auto-init, below Elastix's default
RequiredRatioOfValidSamples of 0.25, so the IMPACT pipeline aborted
on every M6 case. Set the gate to 0.1 so the metric's native
in-buffer evaluation proceeds while still catching diverging
transforms.

Switch the default semantic model from TotalSegmentator (TS/M730.pt)
to MRSegmentator (MRSeg/MRSeg.pt), a multi-modality MRI+CT model
better suited to this cross-modality task. Update docstrings and
REGISTRATION_BENCHMARK.md with the new per-case MI results.

.gitignore

@@ -8,4 +8,6 @@ runs_nnu/
 runs_monai/
 __pycache__/
 
-elastix/
+elastix/
+niftyreg/NiftyReg-Ubuntu-v2.1.1.zip
+niftyreg/NiftyReg-Ubuntu-v2.1.1/

REGISTRATION_BENCHMARK.md

@@ -22,36 +22,50 @@ 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**. Shared rules for
-register on the **full CT / T1c volumes** (no head box) and use each toolkit's
+every pipeline:
-**native automatic transform initialization** (no manual centroid offset):
+
-elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
+- **No mask, no crop / pad**: registration runs on the full input volumes
-antspyx its default center-of-mass alignment, simpleitk starts from identity,
+  exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
-niftyreg (reg_aladin) initializes from the image centres, synthmorph needs
+  padding).
-none. All use **fixed = CT, moving = T1c** except `impact`, which keeps
+- **CT is always the fixed image** (moving = T1c); when both inputs are MR
-**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than
+  the first argument is kept as fixed.
-the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
+
-samples map outside the T1c buffer and the IMPACT metric aborts ("too many
+All pipelines also use each toolkit's **native automatic transform
-samples map outside moving image buffer"). The comparison is on each
+initialization** (no manual centroid offset): elastix / impact use
-toolkit's native optimizer / metric / global-search capability.
+`AutomaticTransformInitialization` (CenterOfGravity), antspyx its default
+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) and
+Listed alphabetically. All register on the **full volumes** (no head box, no
-use each toolkit's **native automatic transform initialization** (no manual
+mask, no crop / pad) and use each toolkit's **native automatic transform
-centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
+initialization** (no manual centroid offset): elastix / impact use elastix's
-(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
+`AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default
-starts from identity, niftyreg initializes from the image centres, and
+center-of-mass alignment, simpleitk starts from identity, niftyreg
-synthmorph needs none (global deep method). All use **fixed = CT, moving =
+initializes from the image centres, and synthmorph needs none (global deep
-T1c** except `impact`, which keeps **fixed = T1c, moving = CT** (swapping
+method). **All use fixed = CT, moving = T1c** (shared rule above; when both
-breaks its TotalSegmentator metric — see below).
+inputs are MR the first argument is fixed). `niftyreg` and `synthmorph`
+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`     | **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`     | **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) |
-| 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 | `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`) |
 | 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) |
 
@@ -62,10 +76,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,
-centered start — the default, comparable to the other tools) and
+identity start — the default, comparable to the other tools) and
 `--simpleitk-mode production` (the full multi-start tilt grid from
-`02_build_m6_dataset.py`, ~20× slower, the most robust but not used for the
+`02_build_m6_dataset.py`, now also on the full volumes per the no-crop rule,
-head-to-head).
+~20× slower, the most robust but not used for the head-to-head).
 
 ## Test cases
 
@@ -108,40 +122,45 @@ 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 ⇔ `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`).
+  production acceptance rule; aligned ⇔ `|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)
+## Results (10 cases, no clip)
 
 Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
-the MR head; higher = better).
+the MR head; higher = better) with the **unclipped** CT, so all rows are
+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 | 35.0 | 0.93 | **0.652 / 0.644** |
+| **antspyx** | ANTsPy 0.6.3 | 9 | 36.6 | 0.93 | **0.653 / 0.644** |
-| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** |
+| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.517 / 0.520** |
-| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** |
+| **impact** | Docker Elastix 5.2 + IMPACT (MRSegmentator) | 8 | 2534 | 0.93 | **0.591 / 0.575** |
-| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** |
+| **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 20.8 | 0.93 | **0.292 / 0.338** |
-| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** |
+| **simpleitk** | SimpleITK 2.5.6 | 9 | 10.6 | 0.93 | **0.651 / 0.641** |
-| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** |
+| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 24.5 | 0.93 | **0.277 / 0.328** |
 
-\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default
+Per-case MI (nats, unclipped score; **bold** = best for the case):
-(MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
-(hours/case) and was not re-run.
-
-Per-case MI (nats; **bold** = best for the case):
 
 | case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
 |---|---|---|---|---|---|---|
 | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
-| m6_22IRRKU6_20240115 | 0.649 | 0.472 | **0.624** | 0.296 | 0.642 | 0.235 |
+| m6_22IRRKU6_20240115 | **0.649** | 0.472 | 0.635 | 0.292 | 0.642 | 0.277 |
-| m6_23XIMO4M_20230222 | 0.655 | 0.416 | 0.474 | **0.634** | 0.652 | 0.510 |
+| m6_23XIMO4M_20230222 | **0.655** | 0.416 | 0.478 | 0.632 | 0.652 | 0.616 |
-| m6_24V62QQM_20230223 | 0.703 | **0.701** | 0.554 | 0.260 | 0.663 | 0.238 |
+| m6_24V62QQM_20230223 | **0.704** | 0.701 | 0.697 | 0.260 | 0.663 | 0.256 |
-| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.419 | 0.493 | **0.602** | 0.516 |
+| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.556 | 0.493 | **0.602** | 0.484 |
-| m6_24V62QQM_20240820 | 0.652 | 0.643 | 0.458 | 0.118 | **0.652** | 0.096 |
+| m6_24V62QQM_20240820 | **0.653** | 0.644 | 0.531 | 0.118 | 0.652 | 0.111 |
-| m6_26TYQI3C_20221027 | 0.657 | 0.653 | 0.511 | 0.608 | **0.695** | 0.604 |
+| m6_26TYQI3C_20221027 | 0.659 | 0.653 | gate-abort | 0.609 | **0.695** | 0.590 |
-| m6_2BBN5XA4_20250411 | 0.576 | 0.287 | 0.155 | 0.069 | **0.577** | 0.060 |
+| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | 0.450 | 0.069 | **0.577** | 0.061 |
-| m6_2C4TKI6H_20240506 | 0.645 | 0.401 | **0.626** | 0.231 | 0.634 | 0.209 |
+| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | 0.629 | 0.233 | 0.634 | 0.231 |
-| m6_2CDNN7BH_20240102 | 0.655 | 0.516 | **0.618** | 0.344 | 0.651 | 0.301 |
+| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | 0.625 | 0.340 | 0.651 | 0.323 |
 
 ### Interpretation
 
@@ -163,52 +182,138 @@ Per-case MI (nats; **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) is the strongest on the hardest cases** — top MI
+- **IMPACT (semantic metric) now runs on these pairs** — the no-mask abort was
-  on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618),
+  elastix's sample-validity gate (build default `RequiredRatioOfValidSamples`
-  where anatomical-segmentation similarity out-resolves intensity MI — but its
+  0.25) rejecting the ~14% in-buffer fraction of the CT head-neck FOV vs the
-  median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is
+  tight T1c crop; the parameter map sets it to 0.1 (see the IMPACT section).
-  also by far the slowest (~14 min/case at the shown budget; the new 2000/5000
+  With MRSegmentator at the 2000/5000 budget it is top-three in MI on 7 of the
-  budget would be on the order of hours/case).
+  8 registered cases (median 0.591): 0.635 on the hardest `m6_22IRRKU6`
+  (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.634, best in class;
+  is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class;
-  `m6_26TYQI3C` 0.608) but falls into the same tilt local minima
+  `m6_26TYQI3C` 0.609) but falls into the same tilt local minima
-  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower
+  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). The v2.1.1 Docker build
-  than elastix at full volumes.
+  on the GPU (`-platf 1`) is ~21 s/case, ~4× slower than elastix (vs ~4 min
+  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.238): its `rigid` model is trained on brain *MRI*, so the ~30°
+  (median 0.277 with its default CT clip; 0.241 unclipped): its `rigid` model
-  tilt defeats it on the repositioned cases (0.06–0.10) even though it matches
+  is trained on brain *MRI*, so the ~30° tilt defeats it on the repositioned
-  the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no
+  cases (0.06–0.11) even though it matches the intensity tools on the
-  centroid init** (global deep method) and is deterministic.
+  cleanest (`m6_26TYQI3C` 0.590). The [0,80] brain-window clip of the CT is
+  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 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s,
+  SimpleITK 10.6 s, **NiftyReg 20.8 s** (v2.1.1, GPU block matching),
-  **NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous
+  SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT is by far the
-  budget; GPU feature extraction with backprop per iteration).
+  slowest — 2534 s/case median at the 2000/5000 benchmark budget, ~276 s/case
+  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 and NiftyReg are deterministic. **IMPACT is the most stochastic**:
+  SynthMorph is deterministic; NiftyReg is deterministic on the CPU
+  (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 (pip/conda-installed here)
+# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1)
 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
 ```
@@ -249,11 +354,9 @@ 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 --n-cases 10 --qa-images
+    --pipelines simpleitk elastix antspyx niftyreg --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.)
@@ -299,20 +402,33 @@ 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 — the release zips for
+suite (KCL-BMEIS, the actively-maintained NiftyReg): **symmetric block
-v2.x need a newer glibc, so this env uses the conda-forge build):
+matching + trimmed least squares** (Modat et al., J. Med. Img. 1(2) 024003,
-**symmetric block matching + trimmed least squares** (Modat et al.,
+2014). Instead of a histogram MI objective it matches 4-voxel image blocks by
-J. Med. Img. 1(2) 024003, 2014). Instead of a histogram MI objective it
+**normalised cross-correlation** at each level of a 3-level Gaussian pyramid
-matches 4-voxel image blocks by **normalised cross-correlation** at each level
+and solves a robust (trimmed) least-squares for the rigid transform — a
-of a 3-level Gaussian pyramid and solves a robust (trimmed) least-squares for
+genuinely different algorithmic family, and cross-modality tolerant because
-the rigid transform — a genuinely different algorithmic family, and
+NCC is invariant to intensity scaling.
-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
-# one-time, in the benchmark env (provides reg_aladin + reg_f3d, NiftyReg 1.5.69)
+# default: build the v2.1.1 Docker image (one-time)
+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
 ```
 
@@ -322,8 +438,10 @@ 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 3` cubic internal warping, default initialisation
+   coarse), `-interp 1` (GPU) / `-interp 3` (CPU) internal warping, default
-   (image centres).
+   initialisation (image centres). GPU: `--gpus all`, `-platf 1`,
+   `-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
@@ -336,12 +454,17 @@ 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 **TotalSegmentator** feature extractor. It is GPU-based and
+TorchScript segmentation feature extractor (**MRSegmentator** `MRSeg/MRSeg.pt`
-Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
+by default; TotalSegmentator `TS/M730.pt` and other models from
-**fixed = T1c, moving = CT** (the reverse of the others): the CT FOV
+`VBoussot/impact-torchscript-models` are selectable). It is GPU-based and
-(~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with
+Docker-based (alongside `synthmorph`). It follows the shared rules (full
-CT as the fixed image too many CT samples map outside the T1c buffer and the
+volumes, no mask, fixed = CT, moving = T1c). The metric only evaluates samples
-IMPACT metric aborts ("too many samples map outside moving image buffer").
+that map inside the moving buffer, and the CT head-neck FOV (~300×300×251 mm)
+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):
@@ -354,25 +477,40 @@ 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/). M730.pt is the recommended TotalSegmentator.
+#    cache (data/impact_models/). MRSeg.pt (MRSegmentator) is the default:
-mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/TS
+#    it is trained on both MRI and CT, the right prior for this
-curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/TS/M730.pt \
+#    cross-modality CT->T1c task (all models share the same 2^l+3 FOV
-    https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/TS/M730.pt
+#    architecture, so the parameter map is unchanged between them).
+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 (M258/M291/M850, MIND, SAM2.1, …) are selectable with
+Other models (TotalSegmentator `TS/M730.pt`, `TS/M258/M291/M850`, MIND, SAM2.1,
-`--impact-model <HF-relative-path-or-abs>`; see the
+Anatomix, …) are selectable with `--impact-model <HF-relative-path-or-abs>`;
+see the
 [impact-torchscript-models](https://huggingface.co/VBoussot/impact-torchscript-models)
-repo for the full list and FOV/layer tables.
+repo for the full list and FOV/layer tables. The model choice does not affect
+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=T1c and moving=CT are used as-is (no head box).
+1. **Full volumes, fixed=CT, moving=T1c** — both volumes are used as-is (no
-   T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too
+   head box, no mask, no crop / pad, per the shared rules above). The metric
-   many CT samples map outside the T1c buffer and the metric aborts.
+   only evaluates spatial samples that map inside the moving buffer, and the
+   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
@@ -381,9 +519,24 @@ 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. GPU auto-spread across the host GPUs by pid
+   ASGD, `Random` sampler, `(RequiredRatioOfValidSamples 0.1)`. A copy of the
-   (`--impact-gpu -1` forces CPU).
+   exact parameter map is dumped to
-4. **Read-back** — the result `TransformParameters.0.itk.txt` is an
+   `data/m6reg/impact/paramdump/ParameterMap_<model>_<pid>.txt` before every
+   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.
@@ -391,13 +544,14 @@ binary directly (`/lib/elastix-install/bin/elastix`):
 ## Notes / caveats
 
 - The benchmark compares each toolkit in its **native configuration**: all
-  register on the full volumes (no head box) and use each toolkit's native
+  register on the full volumes (no head box, no mask, no crop / pad) and use
-  automatic transform initialization (no manual `centroid_offset`) — elastix /
+  each toolkit's native automatic transform initialization (no manual
-  impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's
+  `centroid_offset`) — elastix / impact `AutomaticTransformInitialization`
-  center-of-mass alignment, simpleitk identity, synthmorph none. Production
+  (CenterOfGravity), antspyx's center-of-mass alignment, simpleitk identity,
-  uses `--simpleitk-mode production`'s multi-start tilt grid, which raises
+  synthmorph none. Production uses `--simpleitk-mode production`'s
-  SimpleITK's acceptance at ~20× the cost; that is out of scope for the
+  multi-start tilt grid, which raises SimpleITK's acceptance at ~20× the
-  head-to-head but is the config the dataset actually ships with.
+  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
@@ -439,11 +593,16 @@ binary directly (`/lib/elastix-install/bin/elastix`):
     - **`-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)
+    - **GPU is the default now** (v2.1.1, `-platf 1` via
-      at full volumes with 24 threads.
+      `docker run --gpus all`): ~20–25 s/case on the full volumes vs ~4 min
-    - **Version note.** The upstream KCL-BMEIS release zips (v2.x) need
+      CPU (120–440 s, `-omp` thread-bound) for the conda 1.5.69 build.
-      glibc ≥ 2.33; this host (Ubuntu 20.04, glibc 2.31) uses the conda-forge
+      `--niftyreg-gpu -1` forces the CPU platform (`-platf 0`).
-      build (1.5.69) instead.
+    - **Version note.** The upstream KCL-BMEIS release binaries (v2.x) need
+      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
@@ -470,28 +629,41 @@ binary directly (`/lib/elastix-install/bin/elastix`):
     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
+- **Slow, by design.** Every ASGD iteration runs the segmentation
-      extractor on the sampled patches with backpropagation. The benchmark
+      extractor (MRSegmentator by default) on the sampled patches with
-      budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
+      backpropagation. The benchmark
-      the elastix default) is ~10× the authors' 500-iter / 2000-sample
+budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
-      Jacobian budget, i.e. on the order of *hours* per case on full volumes —
+       the elastix default) is ~10× the authors' 500-iter / 2000-sample
-      the results shown for `impact` are from the previous 500/2000 budget and
+       Jacobian budget: measured 2534 s/case (median) on the no-mask M6
-      were not re-run at the new budget. `--impact-iters` / `--impact-samples`
+       pairs, ~276 s/case at 500/2000. `--impact-iters` / `--impact-samples`
-      trade accuracy for speed.
+       trade accuracy for speed.
-   - **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT
+- **Runs on M6 with a relaxed sample-validity gate.** Like every other
-     follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251
+      pipeline it runs fixed=CT/moving=T1c on the full volumes (no mask, no
-     mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the
+      crop / pad). The CT FOV (~300×300×251 mm) is much larger than the T1c
-     fixed image too many CT samples map outside the T1c buffer and the metric
+      crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
-     aborts ("too many samples map outside moving image buffer"). The returned
+      the samples map inside the T1c moving buffer. The metric calls
-     transform still satisfies the contract (T1c→CT).
+      `CheckNumberOfSamples()` after every evaluation and the image build
+      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.
-   - **Rugged single-init convergence.** The semantic landscape has more local
+- **History: single-init convergence.** The semantic landscape has more
-     minima than Mattes-MI: on full volumes it no longer nails the hardest
+      local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer
-     repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box
+      nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0;
-     used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration
+      the pre-mask head-box flow got it to mid 8.3 / std 13.2). The no-mask
-     budget would likely raise its acceptance; that is out of scope for the
+      MRSegmentator runs show the same pattern: strong MI (0.635 on
-     single-init head-to-head.
+      `m6_22IRRKU6` at 2000/5000) but the pose still misses the concentricity
+      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

niftyreg/Dockerfile

@@ -0,0 +1,35 @@
+# 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"]

scripts/02_build_m6_dataset.py

@@ -1,8 +1,9 @@
 """Build M6-2025 manifests: T1c volumes, labeled when GTV exists.
 
-GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on a head box
+GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on the full
-located in CT space; multi-start initial tilts because inter-session head tilt
+volumes -- no mask, no crop / pad; multi-start initial tilts because
-reaches ~30 deg and tilted MI local minima can outscore the true alignment).
+inter-session head tilt reaches ~30 deg and tilted MI local minima can
+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).
@@ -105,6 +106,10 @@ 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:
@@ -119,7 +124,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)
@@ -171,14 +176,11 @@ 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
+    """Rigid registration with CT as fixed and T1c as moving, Mattes MI on the
-    head-box images.
+    full volumes (no mask, no crop / pad).
 
-    Both inputs share physical space (oriented). The registration grid is a head
+    Both inputs share physical space (oriented). The transform is computed on
-    box located in CT space from the top of the head (the CTs are head-neck
+    the full CT / T1c volumes, so it applies directly to any resampling.
-    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,
@@ -198,38 +200,9 @@ 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)
-    # Registration grid = a head box located in CT space. The CTs are head-neck
+    # full volumes, no mask / crop / pad: fixed = CT, moving = T1c
-    # (often head-neck-chest) exams, so cropping to the MR head box would leave
+    fixed, moving = ct, t1c
-    # 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 = []
 
@@ -261,8 +234,21 @@ 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
 
@@ -276,7 +262,9 @@ def register_ct_to_t1c(ct, t1c, threads=12):
         return tr0, gap
 
     def is_good(g):
-        return bool(g) and g["std_mm"] <= 12.0 and 0.0 <= g["mid_mm"] <= 20.0
+        # |mid| (not 0 <= mid): full-volume registration can land the CT skull
+        # 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):
@@ -325,7 +313,6 @@ 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),
@@ -373,7 +360,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 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
+    if gap and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
         align = "concentric"
     else:
         align = "not_concentric"

scripts/bench_registration.py

@@ -20,6 +20,10 @@ 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
@@ -93,6 +97,14 @@ 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
 # --------------------------------------------------------------------------- #
@@ -152,7 +164,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 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0
+    ok_gap = bool(gap) and abs(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,
@@ -237,7 +249,9 @@ 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)
 
 
@@ -255,6 +269,11 @@ 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)
@@ -380,8 +399,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", {}).get("median") or 0 for p in labels]
+    regt = [(ok_rows[p].get("reg_wall_s") 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]
+    mi = [(ok_rows[p].get("mi_ct_mr") or {}).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, [
@@ -487,6 +506,22 @@ 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")

src/registration_pipelines.py

@@ -23,15 +23,26 @@ 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
+Shared rules for every pipeline:
-toolkit's **native automatic transform initialization** (no manual
+
-``centroid_offset``): elastix / impact use elastix's
+- **No mask, no crop / pad**: registration always runs on the full input
+  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).  All use **fixed = CT, moving = T1c** except IMPACT,
+(global deep method).
-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
@@ -42,9 +53,9 @@ TotalSegmentator metric -- see ``register_impact``).
 - ``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).  full volumes, fixed=CT/moving=T1c (no mask -- shared
-  ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
+  rules above); elastix ``AutomaticTransformInitialization`` (CenterOfGravity)
-  Needs Docker + the image + a model asset.
+  provides the start.  Needs Docker + the image + a model asset.
 - ``register_niftyreg``   : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
   symmetric block-matching + trimmed least squares (Modat et al., J. Med.
   Img. 2014), NCC-based block similarity (cross-modality tolerant).  Rigid
@@ -123,6 +134,50 @@ 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``.
 
@@ -203,8 +258,9 @@ 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)
+    # full volumes, no mask / crop / pad; CT is the fixed image (or, if both
-    fixed, moving = ct, t1c
+    # inputs are MR, the first argument is)
+    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
@@ -236,8 +292,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
         "toolkit": f"SimpleITK {sitk.Version_VersionString()}",
         "metric": "MattesMutualInformation (full volumes)",
         "transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
-                      f"{iters_full} iters/level, identity start)"),
+                       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,
     }
@@ -315,9 +372,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
     after a successful registration.
     """
     itk = _elastix_api()
-    # full volumes, fixed = CT, moving = T1c (no head box)
+    # full volumes, no mask / crop / pad; CT is the fixed image
-    fixed = ct
+    fixed, moving = fixed_moving(ct, t1c)
-    moving = t1c
     params = {
         "ImageDimension": 3,
         "FixedInternalImagePixelType": "float",
@@ -373,7 +429,8 @@ 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": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
+        "fixed": _modality_name(fixed), "moving": _modality_name(moving),
+        "seconds_reg": round(dt, 1),
     }
     return T, meta
 
@@ -408,8 +465,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, fixed = CT, moving = T1c (no head box)
+    # full volumes, no mask / crop / pad; CT is the fixed image
-    fx, mv = ct, t1c
+    fx, mv = fixed_moving(ct, t1c)
     import time
     t0 = time.time()
     # No manual centroid init: antspyx's default initial_transform=None computes
@@ -431,7 +488,8 @@ 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, "aff_iterations": list(aff_iterations),
+        "auto_inverted": inverted, "fixed": _modality_name(fx), "moving": _modality_name(mv),
+        "aff_iterations": list(aff_iterations),
         "aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
     }
     return T, meta
@@ -484,7 +542,8 @@ 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).
@@ -492,7 +551,10 @@ 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.  Being
+    transform, so no skull-stripping or intensity preprocessing is needed --
+    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
@@ -524,11 +586,15 @@ 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:
-        # fixed = CT, moving = T1c (shared convention for all pipelines)
+        # full volumes, no mask / crop / pad; CT is the fixed image
-        sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
+        fixed, moving = fixed_moving(ct, t1c)
-                        os.path.join(tmp, "fixed.nii.gz"), True)      # CT
+        fixed_name = _modality_name(fixed)  # before clipping (a clipped CT
-        sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
+        # is no longer intensity-detectable)
-                        os.path.join(tmp, "moving.nii.gz"), True)     # T1c
+        fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
+        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.
@@ -570,10 +636,11 @@ 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": "CT", "moving": "T1c",
+        "fixed": fixed_name, "moving": _modality_name(moving),
         "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
@@ -587,7 +654,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 = "TS/M730.pt"
+IMPACT_DEFAULT_MODEL = "MRSeg/MRSeg.pt"
 IMPACT_HF_MODEL_URL = (
     "https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/{rel}"
 )
@@ -690,6 +757,14 @@ 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")',
@@ -708,7 +783,8 @@ 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})",
@@ -748,16 +824,25 @@ 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=T1c, moving=CT) in the Docker image
+    """IMPACT rigid CT->T1c (fixed=CT, moving=T1c) 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 TotalSegmentator model (``TS/M730.pt`` by default) is run
+    one: a TorchScript segmentation model (``MRSeg/MRSeg.pt`` by default,
-    on randomly sampled patches of the fixed and moving images and the
+    MRSegmentator -- robust multi-modality segmentation on MRI *and* CT, 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).  This is a genuinely different metric from
+    (L2, per the recommended config).  Intensity standardization (zero mean /
-    the Mattes-MI pipelines above -- it is the cross-modality tool the image was
+    unit variance per the model card) is applied *inside the metric*: at init
-    built for.
+    it computes whole-volume [min, max, mean, sigma] per image and the model's
+    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``
@@ -765,12 +850,17 @@ 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).
 
-    NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other
+    Like every other pipeline this runs with **fixed=CT, moving=T1c** on the
-    pipelines).  The CT FOV (~300x300x251 mm) is much larger than the T1c crop
+    full volumes: no mask, no crop / pad (shared rules above).  The IMPACT
-    (~220x220x150 mm); with CT as the fixed image, too many CT samples map
+    metric draws random spatial samples from the fixed (CT) buffer and
-    outside the T1c buffer and the IMPACT metric aborts ("too many samples map
+    evaluates them only where they map inside the moving (T1c) buffer; when
-    outside moving image buffer").  The returned transform still satisfies the
+    too few samples are in-buffer, elastix aborts ("too many samples map
-    unified contract (T1c->CT).
+    outside moving image buffer").  On the M6 head-neck CTs (FOV
+    ~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'
@@ -792,13 +882,10 @@ 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 --
+    # full volumes, no mask / crop / pad; CT is the fixed image (or, if both
-    # elastix's AutomaticTransformInitialization aligns the image centers).
+    # inputs are MR, the first argument).  No manual init -- elastix's
-    # T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger
+    # AutomaticTransformInitialization aligns the image centers.
-    # than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples
+    fixed, moving = fixed_moving(ct, t1c)
-    # 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:
@@ -806,8 +893,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 box's (T1c) base spacing
+    # the finest level matches the moving image's base spacing
-    base = float(np.median(t1c.GetSpacing()))
+    base = float(np.median(moving.GetSpacing()))
     n = len(scales)
     voxel_sizes = [base * s for s in scales]
     # image pyramid: IMPACT's docs recommend NO image resampling (all-1 schedules)
@@ -840,10 +927,22 @@ 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)
-        with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
+        param_text = _impact_param_map(os.path.join("/Data/Models", model_rel),
-            f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
                                        gpu, iters, samples, voxel_sizes,
-                                       rescale=schedule, smoothing=smoothe))
+                                       rescale=schedule, smoothing=smoothe)
+        with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
+            f.write(param_text)
+
+        # debug aid: keep a persistent copy of the exact parameter map elastix
+        # 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:
@@ -870,26 +969,28 @@ 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); it maps
+        # (AutomaticTransformInitialization is folded in by elastix); elastix
-        # fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT);
+        # maps moving(T1c) -> fixed(CT), which is the contract direction;
-        # fix_orientation resolves the direction.
+        # fix_orientation verifies it against the skull ring.
         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 (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)",
+        "metric": f"IMPACT semantic ({family} {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": "T1c", "moving": "CT",
+        "fixed": _modality_name(fixed), "moving": _modality_name(moving),
         "seconds_reg": round(dt, 1),
         "docker_image": IMPACT_IMAGE, "model": model_rel, "gpu": gpu,
-        "iters": iters, "samples": samples,
+        "iters": iters, "samples": samples, "param_dump": param_dump,
     }
     if keep_scratch:
         meta["scratch"] = tmp
@@ -901,6 +1002,9 @@ 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")
@@ -913,6 +1017,26 @@ 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,
@@ -938,7 +1062,8 @@ 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
@@ -951,32 +1076,96 @@ 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.
 
-    Full volumes (no head box); default initialisation (image centres);
+    The CT is clipped to a brain window (``clip_ct=(0.0, 80.0)`` by default,
-    ``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3`` (cubic) internal
+    per the CT-clipping A/B) before registration; pass ``clip_ct=None`` to
-    warping.  The output affine is documented as ``Affine * Reference =
+    register the CT as-is.  MR inputs are never clipped.
-    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
+    Full volumes (no head box, no mask / crop / pad); default initialisation
-    (``conda install -c conda-forge niftyreg``).
+    (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.
+
+    ``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
 
-    binpath = _niftyreg_bin()
+    if backend not in ("auto", "conda", "docker"):
-    # full volumes, reference = CT, floating = T1c (no head box)
+        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
+        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(ct, sitk.sitkFloat32),
+        sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
                         os.path.join(tmp, "ref.nii.gz"), True)
-        sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
+        sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
                         os.path.join(tmp, "flo.nii.gz"), True)
-        cmd = [binpath, "-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
+        flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
-               "-rigOnly", "-aff", "aff.txt",
+                 "-rigOnly", "-aff", "aff.txt",
-               "-ln", str(levels), "-maxit", str(maxit),
+                 "-ln", str(levels), "-maxit", str(maxit),
-               "-interp", "3", "-omp", str(threads), "-voff"]
+                 "-interp", interp, "-omp", str(threads)]
+        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)
@@ -999,14 +1188,16 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
 
     meta = {
         "pipeline": "niftyreg",
-        "toolkit": _niftyreg_version(binpath),
+        "toolkit": toolkit,
+        "backend": "docker" if use_docker else "conda",
         "metric": "NCC block matching + trimmed least squares (symmetric, full volumes)",
         "transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, "
                       f"{maxit} TLS iters/level, cubic interp)"),
         "convention": "T1c->CT (verified)", "auto_inverted": inverted,
-        "fixed": "CT", "moving": "T1c",
+        "fixed": fixed_name, "moving": _modality_name(moving),
         "seconds_reg": round(dt, 1),
-        "levels": levels, "maxit": maxit,
+        "levels": levels, "maxit": maxit, "interp": interp, "gpu": gpu,
+        "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
     }
     if keep_scratch:
         meta["scratch"] = tmp
@@ -1047,11 +1238,20 @@ def availability():
                     out[name] = (False, f"ants(antspyx) not importable ({type(e).__name__})")
             elif name == "niftyreg":
                 b = shutil.which("reg_aladin")
-                if b:
+                dok, dnote = _niftyreg_docker_ok()
-                    out[name] = (True, _niftyreg_version(b))
+                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, 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 "
+                    out[name] = (False, "reg_aladin not on PATH and no docker image; "
-                                        "(conda install -c conda-forge niftyreg)")
+                                        "conda install -c conda-forge niftyreg or "
+                                        f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/")
             elif name == "synthmorph":
                 ok, note = _synthmorph_docker_ok()
                 if not ok: