feat(benchmark): add ANTs deformable variant pipelines

Extend registration_pipelines and bench_registration with stock antspyx ANTs presets (SyNRA, SyN, Repro sr/br, QuickRepro sr/br) plus the legacy Rigid -> SyNOnly workflow, and update the deformable benchmark documentation with follow-up ANTs-variant runs and results.
Refactor code structure for improved readability and maintainability
2026-10-03 13:26:12 +08:00 · 2026-10-03 13:25:43 +08:00 · 2026-10-02 19:00:07 +08:00 · 2026-10-02 18:59:08 +08:00 · 2026-09-30 21:23:18 +08:00
10 changed files with 3691 additions and 1189 deletions
--- a/AGENTS.md
+++ b/AGENTS.md
@ -10,7 +10,7 @@ Activate with:
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 ```
-Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), `impact-reg-konfai` (pip; IMPACT — on this host the elastix binary needs the one-time shim from `scripts/setup_impact_konfai.sh`), plus the Docker image `freesurfer/synthmorph`.
+Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK_RIGID.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), plus the Docker images `freesurfer/synthmorph` and `vboussot/elastix_impact` (IMPACT; rigid MIND + TotalSegmentator config, feature models under `data/impact_models`). The benchmark's ImpactLoss scoring metric uses `src/impact_loss.py` (vendored from vboussot/ImpactLoss; runs on the env's torch).
 ## Project
--- a/REGISTRATION_BENCHMARK.md
+++ b/REGISTRATION_BENCHMARK.md
@ -1,726 +0,0 @@
 # Registration Pipeline Benchmark — CT → T1c (M6)
 Benchmarks **six cross-modality registration pipelines** (five rigid;
 `impact` is the native deformable B-spline preset of the ImpactReg package)
 that align a head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be
 warped onto the T1c (MR) native grid. This is the registration step of
 `scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT
 Hounsfield vs MR intensity) and the patient is repositioned between the CT and
 MR sessions, so the inter-session head tilt can reach ~30°.
 - Pipelines: [`src/registration_pipelines.py`](src/registration_pipelines.py)
 - Benchmark harness: [`scripts/bench_registration.py`](scripts/bench_registration.py)
 - Results: `results/bench_registration/` (per-pipeline JSON, `summary.json`,
  `summary.md`, `bench_registration.png`)
 ## Unified contract
 Every pipeline returns a SimpleITK transform `T` such that
 ```python
 sitk.Resample(ct, t1c, T, sitk.sitkLinear)              # CT  onto the MR grid
 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
 every pipeline:
 - **No mask, no crop / pad**: registration runs on the full input volumes
  exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
  padding).
 - **The T1c (MR) is always the fixed image** (moving = CT) for every
  pipeline.  `impact`'s ImpactReg **MR_CT** convention coincides with this
  rule: its output displacement field (defined on the fixed grid) is the
  contract transform T1c→CT directly — see the IMPACT section.  (A CT clipped
  to a non-negative window reads as MR, so the T1c slot decides; both-CT
  pairs keep the given order.)
 All pipelines also use each toolkit's **native automatic transform
 initialization** (no manual centroid offset): elastix / impact use
 `AutomaticTransformInitialization` (CenterOfGravity), antspyx its default
 center-of-mass alignment, simpleitk starts from identity, 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. With the ImpactReg MR_CT direction the samples
 are drawn from the tight T1c head crop (fixed, ~220×220×150 mm) and map mostly
 inside the large CT head-neck FOV (moving, ~300×300×251 mm), so the gate
 rarely fires; the pipeline still sets `RequiredRatioOfValidSamples 0.1` (vs
 the build default 0.25) as a safety override for the FOV mismatch — 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
 mask, no crop / pad) and use each toolkit's **native automatic transform
 initialization** (no manual centroid offset): elastix / impact use elastix's
 `AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default
 center-of-mass alignment, simpleitk starts from identity, niftyreg
 initializes from the image centres, and synthmorph needs none (global deep
 method). **All use fixed = T1c, moving = CT** (shared rule above;
 `impact`'s ImpactReg MR_CT convention coincides with it). `niftyreg` and
 `synthmorph` clip the CT (the moving image) 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`    | T1c / CT | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
 | 2 | `elastix`    | T1c / CT | 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** (MR_CT) | auto CenterOfGravity (override) | **IMPACT semantic** (MIND `R1D2_3D` L1 + TotalSegmentator `M730` deepest-layer Dice, per-level weighted, Static) + Mattes MI (w 0.4) + bending energy | B-spline DVF (preset native), ASGD, 4 scales | ✅ ran (impact-reg-konfai + GPU; `Elastix_IMPACT_Static` preset — see IMPACT section) |
 | 4 | `niftyreg`   | T1c / CT | 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`  | T1c / CT | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
 | 6 | `synthmorph` | T1c / CT | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
 Benchmark budget (best-results settings): `elastix` and `simpleitk` run with
 **MaximumNumberOfIterations = 2000** (per resolution / per pyramid level),
 with the sampler and metric left at their defaults (**NumberOfSpatialSamples
 5000**, **NumberOfHistogramBins 32**). `impact` runs the published ImpactReg
 preset **`Elastix_IMPACT_Static`** (the `impact-reg-konfai` package) — an
 elastix B-spline (4 scales, FinalGridSpacing 14 mm, per-level
 **MaximumNumberOfIterations 400/300/200/200**, **NumberOfSpatialSamples
 2000**) driven by the IMPACT metric: a per-level **MIND (L1) +
 TotalSegmentator M730 (deepest layer, Dice)** combination, layer weights
 0.2/0.8 → 0.7/0.3, Static mode, + Mattes MI weight 0.4 + bending energy —
 with CenterOfGravity auto-init and `RequiredRatioOfValidSamples 0.1` applied
 as raw-elastix overrides for the M6 pairs. `antspyx`, `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
 `--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,
 ~20× slower, the most robust but not used for the head-to-head).
 ## Test cases
 First 10 labeled M6 cases (`--n-cases 10 --sample first`), reproducible. 9 have
 a complete CT; 1 (`m6_22GD5VL2_20230724`) has an incomplete CT head and is
 reported as `incomplete-CT` for every pipeline (a data issue, not a pipeline
 failure).
 ```
 m6_22GD5VL2_20230724  (incomplete CT)
 m6_22IRRKU6_20240115  m6_23XIMO4M_20230222  m6_24V62QQM_20230223
 m6_24V62QQM_20231017  m6_24V62QQM_20240820  m6_26TYQI3C_20221027
 m6_2BBN5XA4_20250411  m6_2C4TKI6H_20240506  m6_2CDNN7BH_20240102
 ```
 ## Metrics (identical for every pipeline)
 The headline alignment score is **`mi_ct_mr`** — the registration's *achieved*
 Mattes MI, read off as a standalone score:
 - **`mi_ct_mr`** — mean Mattes MI (nats) between the **CT inverse-transformed
  into MR coordinates** (the warped CT on the T1c grid) and the **original MR**,
  restricted to the MR head and each image percentile-binned independently
  (SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are
  better in-sync = better aligned; it collapses toward 0 when the CT is
  rotated/translated away from the MR. Tables and plots report this (plus
  registration time and GTV-in-CT coverage).
 - **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT
  head tissue (the slab must straddle the lesion).
 - **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring).
 The production geometric QA is still computed per case and kept in the JSON
 (it drives dataset acceptance in `02_build_m6_dataset.py`, but its summary
 numbers were found to correlate poorly with perceived alignment, so they are no
 longer reported in the tables/plots):
 - **`skull_ring_gap_mm`** — march rays from the MR head center; the gap between
  the MR head silhouette boundary and the warped-CT skull ring (p97 of positive
  CT) at mid and low levels, plus the angular std. `mid` measures
  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` —
  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)
 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
 comparable on the same score scale. Full volumes, no mask / crop / pad,
 fixed = T1c / moving = CT (shared rule; `impact`'s MR_CT convention
 coincides with it). `niftyreg` and `synthmorph` clip the CT (the moving
 image) 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` runs the `impact-reg-konfai` package on the GPU (see the IMPACT
 section).
 | pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
 |---|---|---|---|---|---|
 | **antspyx** | ANTsPy 0.63 | 9 | 15.5 | 0.93 | **0.652 / 0.650** |
 | **elastix** | ITK-Elastix 0.25.4 | 9 | 9.8 | 0.93 | **0.646 / 0.638** |
 | **impact** | impact-reg-konfai (Elastix 5.2 + IMPACT, `Elastix_IMPACT_Static` preset, MIND+M730, deformable) | 9 | 216.4 | 0.93 | **0.578 / 0.558** |
 | **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 16.2 | 0.93 | **0.292 / 0.339** |
 | **simpleitk** | SimpleITK 2.5.6 | 9 | 9.6 | 0.93 | **0.647 / 0.642** |
 | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 27.9 | 0.93 | **0.277 / 0.328** |
 Per-case MI (nats, unclipped score; **bold** = best for the case). All rows
 are re-runs under the **fixed = T1c / moving = CT** convention (the shared
 rule); the CT-fixed-era numbers are in the git history. The `impact` column
 is the `impact-reg-konfai` deformable `Elastix_IMPACT_Static` preset (the
 retired rigid Docker config is archived in
 `results/bench_registration/impact_docker_legacy.json`).
 | 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.648** | 0.644 | 0.542 | 0.292 | 0.647 | 0.277 |
 | m6_23XIMO4M_20230222 | **0.654** | 0.652 | 0.543 | 0.632 | **0.654** | 0.616 |
 | m6_24V62QQM_20230223 | 0.703 | 0.698 | 0.667 | 0.260 | **0.704** | 0.256 |
 | m6_24V62QQM_20231017 | **0.601** | 0.590 | 0.587 | 0.493 | 0.597 | 0.484 |
 | m6_24V62QQM_20240820 | **0.652** | 0.646 | 0.612 | 0.118 | 0.638 | 0.111 |
 | m6_26TYQI3C_20221027 | **0.718** | 0.650 | 0.717 | 0.609 | 0.670 | 0.590 |
 | m6_2BBN5XA4_20250411 | 0.577 | 0.577 | 0.258 | 0.069 | **0.579** | 0.061 |
 | m6_2C4TKI6H_20240506 | **0.645** | 0.638 | 0.578 | 0.233 | 0.636 | 0.231 |
 | m6_2CDNN7BH_20240102 | **0.655** | 0.650 | 0.523 | 0.340 | 0.652 | 0.323 |
 ### Interpretation
 - **ANTsPy leads, with SimpleITK and Elastix within 0.006 of it** (median MI
  0.652 / 0.647 / 0.646): on the well-positioned cases the three cluster at
  0.60–0.70 and all now survive the harder repositioned cases (0.58–0.65 on
  `m6_2BBN5XA4` / `m6_22IRRKU6`).
 - **Switching to fixed = T1c / moving = CT lifted Elastix sharply**
  (0.517 → 0.646 median, and it is now #3 rather than #4). The Mattes-MI
  pipelines sample spatial pairs from the fixed buffer and discard the ones
  that map outside the moving image; with the small T1c head crop as fixed and
  the large CT head+neck FOV as moving, (nearly) all fixed samples are
  in-buffer, versus only a minority when the big CT was fixed. Elastix's
  repositioned cases jumped (`m6_2BBN5XA4` 0.287→0.577, `m6_2C4TKI6H`
  0.401→0.638, `m6_23XIMO4M` 0.416→0.652). The other pipelines are essentially
  unchanged by the swap (antspyx 0.653→0.652, simpleitk 0.651→0.647, niftyreg
  0.292→0.292, synthmorph 0.277→0.277) — they were not sample-coverage
  limited.
 - **The grossly repositioned cases still separate the leaders from the
  tail.** On `m6_2BBN5XA4` / `m6_24V62QQM_20240820` (~30° inter-session tilt)
  ANTs/SimpleITK/Elastix hold 0.58–0.65, while niftyreg (0.07–0.12) and
  synthmorph (0.06–0.11) — the single-start block-matching and pretrained-MRI
  methods — lock onto local minima. This is exactly why production
  `02_build_m6_dataset.py` uses a **multi-start tilt grid**.
 - **IMPACT now runs the published `Elastix_IMPACT_Static` preset**
  (`impact-reg-konfai` package; the old Docker backend is retired): per-level
  **MIND + TotalSegmentator M730** combination (L1 + deepest-layer Dice,
  weights 0.2/0.8 → 0.7/0.3, Static mode, whole-image features), the preset's
  native deformable B-spline (400/300/200/200 iterations, 2000 samples,
  + Mattes MI 0.4 + bending energy), registered in the MR_CT direction
  (fixed=T1c) with CenterOfGravity auto-init and the 0.1 sample-validity gate
  as overrides. It is the best of the IMPACT variants on M6 and now **second
  overall at median MI 0.578** (behind SimpleITK/ANTsPy 0.65; ahead of
  elastix 0.52), a clear step up from the old rigid Docker config (median
  0.467, archived in `impact_docker_legacy.json`): the deformable refinement
  + MR_CT direction lifts the repositioned cases hardest (`m6_2C4TKI6H`
  0.334→0.578, `m6_23XIMO4M` 0.362→0.543, `m6_22IRRKU6` 0.467→0.542,
  `m6_2CDNN7BH` 0.332→0.523; `m6_26TYQI3C` 0.717 — now the best tool on that
  case) at a small cost on two well-aligned ones (`m6_24V62QQM_20230223`
  0.700→0.667, `m6_24V62QQM_20240820` 0.646→0.612, still accepted).
  Concentricity also improved (median mid 12.0 / std 21.9 mm vs 21.1 / 28.0;
  3/9 accepted). It is the slowest of the six (~216 s/case, the deformable
  optimization + KonfAI wrapper; the old rigid config was ~78 s). The
  `m6_2BBN5XA4` ~30° tilt case is still its weakest (0.258) — the single
  CenterOfGravity start remains the limiting factor. Earlier Docker-era
  variants: MIND+MRSegmentator median 0.417 (~72 s), single-model MIND
  median 0.406 (~37 s), MRSegmentator-only Jacobian configs 0.631
  (500/2000) and 0.591 (2000/5000) at ~5×/~35× the cost (archived:
  `results/bench_registration/impact_marseg_*.json`, `impact_2000_5000.json`,
  `impact_500_2000_explicit_smoothing.json`, `impact_mind_500_ref.json`).
 - **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
  is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class;
  `m6_26TYQI3C` 0.609) but falls into the same tilt local minima
  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). The v2.1.1 Docker build
  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.277 with its default CT clip; 0.241 unclipped): its `rigid` model
  is trained on brain *MRI*, so the ~30° tilt defeats it on the repositioned
  cases (0.06–0.11) even though it matches the intensity tools on the
  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, T1c-fixed re-run):
  SimpleITK 9.6 s, elastix 9.8 s, **ANTs 15.5 s** (down from 36.6 s — the
  smaller fixed image also shaves ANTs' per-level work), **NiftyReg 16.2 s**
  (v2.1.1, GPU block matching), SynthMorph 27.9 s (GPU forward pass),
  **IMPACT 216.4 s** (`Elastix_IMPACT_Static` preset — Static feature mode, no
  backprop, but the deformable B-spline optimization at 1100 total iterations
  + the KonfAI wrapper; the retired rigid Docker config was 78.0 s). The
  earlier MRSegmentator Jacobian configs were far slower still: ~354 s/case at
  500/2000, 2534 s/case at 2000/5000 (GPU feature extraction with backprop
  per iteration).
 - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
  per-case numbers carry run-to-run variance (this table is one draw);
  SimpleITK's MI (and registration) is deterministic at fixed thread count;
  the in-process Elastix ASGD is deterministic given its fixed seed;
  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 in the Jacobian preset
  (`Elastix_IMPACT_Jacobian`) is the most stochastic**: it draws random
  patches through the network every iteration, so the final pose changes run
  to run — treat such numbers as one sample from a distribution. The default
  Static preset re-samples coordinates every iteration as well, but the
  feature maps are deterministic given the pose.
 ## 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)
 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
 # impact (pip; GPU; ~120-220 s/case; one-time elastix shim on hosts too old
 # for the stock binary — see the IMPACT section). --impact-preset selects
 # another published preset (e.g. Elastix_IMPACT_Jacobian, Generic_Rigid);
 # --impact-iters / --impact-samples are no-ops for the default preset's
 # per-resolution matrix
 python scripts/bench_registration.py --pipelines impact --n-cases 10 \
    --workers 3 --threads 24
 # synthmorph (Docker, needs the pulled image; see the SynthMorph section)
 python scripts/bench_registration.py --pipelines synthmorph --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
 ```
 Other flags: `--sample random --seed 0` (reproducible random 10),
 `--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`,
 `--threads T` (also the OpenMP thread count for `reg_aladin`),
 `--simpleitk-mode {benchmark,production}`,
 `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`,
 `--synthmorph-extent <int>`, `--impact-preset <name>` (default
 `Elastix_IMPACT_Static`), `--impact-gpu <int|-1>`, `--impact-iters N`
 (preset max_iterations override; no-op for matrix presets),
 `--impact-samples N` (preset spatial_samples override),
 `--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
 ## Visual QA images (checkerboard + color overlay)
 `--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to
 `<out-dir>/qa/<pipeline>__<case>.png` for eyeballing the alignment.  Each image
 is a 3×2 grid: **rows** are the three orthogonal mid-slices through the head
 center (**mid-axial, mid-coronal, mid-sagittal**), and **columns** are:
 - **left — checkerboard blend:** alternating blocks of the CT inverse-transformed
  into MR coordinates and the original MR; a clean registration gives a seamless
  mosaic, a misregistration shows hard edges at the block boundaries.
 - **right — additive color overlay:** CT in red, MR in green. Structures that
  coincide blend to yellow; rotation/translation error appears as red/green
  fringing at the skull and brain edges.
 The CT is shown in a **brain window (L40/W80 = 0–80 HU)** so the brain
 parenchyma is visible and comparable to the MR (the skull, ~800+ HU, clips to
 white); the MR uses a 2–98 % head window. The coronal and sagittal views are
 oriented superior-up to match the axial view.
 The title carries the pipeline, case key and the achieved Mattes MI, so a
 glance tells you whether the visuals agree with the score. Generate for the
 available pipelines with:
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 python scripts/bench_registration.py \
    --pipelines simpleitk elastix antspyx niftyreg --n-cases 10 --qa-images
 python scripts/bench_registration.py \
    --pipelines synthmorph --n-cases 10 --qa-images
 ```
 (9 images each; the incomplete-CT case is skipped.)
 ## SynthMorph pipeline (Docker) — one-time setup
 The `synthmorph` pipeline runs the [`freesurfer/synthmorph`](https://hub.docker.com/r/freesurfer/synthmorph)
 image: [SynthMorph](https://synthmorph.io) (Hoffmann et al., IEEE TMI 2022), an
 anatomy-aware, acquisition-agnostic **deep** registration (TensorFlow
 `VxmAffineFeatureDetector` with the pretrained `synthmorph.rigid.1.h5` weights).
 No skull-stripping or intensity preprocessing is needed: each volume is
 normalized to an isotropic 1-mm LIA network space internally. It is the
 second Docker-based pipeline and the only deep-learning one.
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
 ```bash
 # 1. Docker (with NVIDIA container toolkit so --gpus works)
 docker --version
 # 2. Pull the image (~11 GB, one-time)
 docker pull freesurfer/synthmorph
 ```
 No extra model assets: the rigid weights ship inside the image
 (`/freesurfer/models`).
 ### How the adapter drives the image
 1. **Rigid matrix model, no init** — `register -m rigid` on GPU when one is
   visible (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces
   CPU). SynthMorph is a global deep method, so it recovers the gross
   inter-session offset on its own and needs **no** initial transform.
 2. **Read-back** — the model estimates first→second, so the pair is passed as
   `(fixed=T1c, moving=CT)` and the `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in
   shared world space) is the T1c→CT contract transform directly; it is
   parsed into a SimpleITK `AffineTransform` and run through the usual
   `fix_orientation` gate.
 3. **GPU sharing** — the container is run as the host uid (clean scratch
   removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
   concurrent benchmark workers can share a GPU.
 ## NiftyReg pipeline — one-time setup
 The `niftyreg` pipeline runs `reg_aladin` from the [NiftyReg](https://github.com/KCL-BMEIS/niftyreg)
 suite (KCL-BMEIS, the actively-maintained NiftyReg): **symmetric block
 matching + trimmed least squares** (Modat et al., J. Med. Img. 1(2) 024003,
 2014). Instead of a histogram MI objective it matches 4-voxel image blocks by
 **normalised cross-correlation** at each level of a 3-level Gaussian pyramid
 and solves a robust (trimmed) least-squares for the rigid transform — a
 genuinely different algorithmic family, and cross-modality tolerant because
 NCC is invariant to intensity scaling.
 **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)
 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
 ```
 ### How the adapter drives the tool
 1. **Rigid 6-DOF, full volumes** — `reg_aladin -ref T1c -flo CT -rigOnly`
   (the T1c is the reference/fixed image, the CT the floating/moving one;
   the CT is clipped to [0,80] HU first). 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 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`; the adapter parses it into
   a SimpleITK `AffineTransform` and runs it through the usual
   `fix_orientation` gate (skull-ring concentricity), which resolves the
   direction empirically rather than trusting the docs. Under
   `-ref T1c -flo CT` the file's affine comes out the other way, so the gate
   inverts it on the M6 cases (8/9 `auto_inverted: True`) to the contract
   transform (T1c→CT).
 ## IMPACT pipeline (impact-reg-konfai) — one-time setup
 The `impact` pipeline runs the
 [impact-reg-konfai](https://pypi.org/project/impact-reg-konfai/) package:
 [IMPACT](https://arxiv.org/abs/2503.24121) registration presets from the
 [`VBoussot/ImpactReg`](https://huggingface.co/VBoussot/ImpactReg) Hugging Face
 repo, driven through the [KonfAI](https://github.com/fideus-labs/KonfAI)
 runtime — no Docker. The default preset is **`Elastix_IMPACT_Static`**:
 Elastix 5.2 with the IMPACT *semantic* similarity metric — a per-level
 combination of the **MIND** handcrafted descriptor (`MIND/R1D2_3D.pt`, L1,
 subset 32) and **TotalSegmentator** `TS/M730.pt` (deepest layer, Dice, subset
 64), layer weights 0.2/0.8 → 0.7/0.3 from coarse to fine, **Static** mode
 (whole-image features, no backprop), a 4-scale B-spline (FinalGridSpacing
 14 mm), 400/300/200/200 iterations per level, 2000 spatial samples, + Mattes
 MI (w 0.4) + bending energy — the published successor of the `MR_CT_TS`
 reference parameter map this repo used to hand-build (and replace the
 `vboussot/elastix_impact` Docker image). Other presets are selectable with
 `--impact-preset` (`Elastix_IMPACT_Jacobian`, `Generic_Rigid`,
 `Generic_Rigid_BSpline`, `FireANTs_Anatomix/IMPACT/SyN`,
 `ConvexAdam_Coarse/Composite`; the live database is listed by
 `impact-reg-konfai register --help`). It is GPU-based (like `synthmorph`),
 but pip-installed.
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
 ```bash
 # in the longitudinal-313 env (Python >= 3.11):
 pip install impact-reg-konfai          # CLI + KonfAI runtime (pulls h5py etc.)
 bash scripts/setup_impact_konfai.sh    # one-time, this host only: the elastix shim
 ```
 The CLI then auto-downloads the one-time assets on first `register`: the
 preset app (HF `VBoussot/ImpactReg`), the TorchScript feature models (HF
 `VBoussot/impact-torchscript-models`: MIND ~15 KB + TS/M730 ~200 MB) and the
 **elastix-impact 1.0.0** binary (GitHub `vboussot/ImpactElastix`, ~40 MB)
 into `~/.cache` (HF cache + `~/.cache/konfai`).
 ### Why the elastix shim (this host only)
 The stock binary needs **glibc ≥ 2.35** and links **LibTorch 2.8.0**; this
 host (Ubuntu 20.04, glibc 2.31, env torch 2.14) satisfies neither, so
 `scripts/setup_impact_konfai.sh` builds a self-contained shim at
 `~/.cache/konfai/elastix-impact-shim`: the glibc 2.35 loader + libs extracted
 from the `niftyreg:2.1.1-ubuntu` image (docker is needed once for that copy
 only), the self-contained LibTorch 2.8.0 cu128 from download.pytorch.org
 (~3.8 GB zip), and a wrapper `bin/elastix` that runs the stock binary through
 the shim loader with the LibTorch 2.8.0 library path (stripping the inherited
 `LD_LIBRARY_PATH`, so the env's torch 2.14 cannot leak in). The pipeline
 selects the shim via `KONFAI_ELASTIX_DIR` (`_impact_elastix_root`); a host
 new enough to run the stock binary directly needs no shim at all.
 ### How the adapter drives the CLI
 1. **Full volumes, fixed=T1c (MR), moving=CT** — the ImpactReg **MR_CT**
   convention, which coincides with the shared rule (the T1c is the fixed
   image): the preset's output is a displacement field on the **fixed** grid
   mapping fixed→moving, so the field is the contract transform T1c→CT
   **directly** (displacement fields cannot be inverted with `GetInverse`,
   which a fixed=CT direction would require). It also keeps
   the metric's sample-validity gate healthy: samples drawn from the tight
   T1c head crop (~220×220×150 mm) map mostly inside the large CT head-neck
   FOV (~300×300×251 mm), unlike fixed=CT where only ~14% were in-buffer at
   the auto-init.
 2. **M6 robustness overrides** — the preset's map is tuned for well-aligned
   pairs, so the adapter forwards raw elastix lines through the package's
   `parameter_overrides` knob: `AutomaticTransformInitialization=true` +
   `AutomaticTransformInitializationMethod=CenterOfGravity` (the preset's
   GeometricalCenter start leaves the intensity centroids ~50-100 mm apart on
   the head+neck FOV mismatch) and `RequiredRatioOfValidSamples 0.1` (build
   default 0.25; keeps the divergence gate while letting the FOV-mismatched
   runs proceed).
 3. **`--fields-only`** — the adapter derives the warp itself (the package's
   derived "Moved" image would be an extra full-size resample of the same
   content).
 4. **GPU auto-spread** across the host GPUs by pid (`--impact-gpu -1` forces
   CPU; the engine then turns mixed precision off).
 5. **Read-back** — the transform comes from `<out>/P000/Transform.h5` (an
   ITK-form HDF5 `DisplacementFieldTransform` on the T1c grid; the SimpleITK
   build here has no ITK HDF5 transform reader, so it is parsed with h5py)
   and direction-checked by `_orient_impact`: the skull-ring gap of the field
   vs a sign-flipped (first-order reverse) proxy — the `fix_orientation`
   analogue for a field with no `GetInverse`.
 6. **Intensity normalization happens inside the metric / model** (no
   pipeline-side normalization): the MIND descriptor normalizes each 5-voxel
   neighbourhood to [0, 1] locally (modality-agnostic by construction); the
   segmentation models (TS/MRSeg) are standardized to zero mean / unit
   variance per the model card. Pre-normalizing the volumes would be an
   exact no-op for the z-scored models and would only corrupt the 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
  each toolkit's native automatic transform initialization (no manual
  `centroid_offset`) — elastix / impact `AutomaticTransformInitialization`
  (CenterOfGravity), antspyx's center-of-mass alignment, simpleitk identity,
  synthmorph none. Production uses `--simpleitk-mode production`'s
  multi-start tilt grid, which raises SimpleITK's acceptance at ~20× the
  cost; that is out of scope for the head-to-head but is the config the
  dataset actually ships with.
 - `synthmorph` is deterministic (one fixed pretrained model, no random
  subsampling), so its per-case numbers are stable run to run — unlike ANTs /
  IMPACT. Its low MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
  genuine cross-modality / tilt weakness, not variance.
 - `fix_orientation` (elastix/ants/niftyreg/synthmorph) chooses between a
  transform and its inverse by skull-ring concentricity — a correctness gate
  that makes the result independent of each toolkit's transform-direction
  convention, so the fixed/moving swap is absorbed by the gate rather than
  re-derived per toolkit. Under the shared fixed=T1c / moving=CT rule the raw
  transforms come out in the contract direction for elastix (combination
  transform = fixed→moving), antspyx (ANTs forward transform) and synthmorph
  (the `-T` LTA, model estimates first→second) — completed runs report
  `auto_inverted: False` for those; `niftyreg`'s `-aff` file comes out the
  other way and the gate inverts it (8/9 `auto_inverted: True`). `impact`
  returns a displacement field (no `GetInverse` in ITK), so it uses
  `_orient_impact` — the same skull-ring decision against a sign-flipped
  (first-order reverse) proxy field; its MR_CT field is in the contract
  direction (`auto_inverted: False`).
 - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU):
    - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain
      *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly
      defeated by the ~30° inter-session tilt (7/9 cases have `std > 15 mm`).
      It is a legitimate data point for how far a pretrained MRI deep model
      generalizes out-of-domain to CT — not a bug.
    - **No centroid init needed.** Being a global deep method it recovers the
      gross offset on its own; the `-i` initial-transform path is not used
      (which also sidesteps this image's float64/float32 bug on the `-i` code
      path).
    - **GPU sharing.** The container runs as the host uid with
      `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
 - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK
       (~10 s) / ANTs (~37 s) and IMPACT (~216 s, deformable preset).
 - **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`,
   symmetric block matching + TLS, CPU/OpenMP):
    - **Different algorithmic family.** No histogram MI objective: 4-voxel
      blocks are matched by normalised cross-correlation and the rigid
      transform is solved by a trimmed least squares, iterated per pyramid
      level (`-ln 3 -maxit 15`). NCC's scale invariance is what makes it
      cross-modality tolerant.
    - **Symmetric mode matters.** The default symmetric variant (optimising
      forward and backward transforms concurrently) is the robust one on these
      repositioned pairs; `-noSym` drifts grossly (gap std 45 mm on a clean
      case every other tool passes).
    - **`-maxit 5` (the tool default) is too coarse** — on the clean case
      `m6_26TYQI3C` it reaches MI 0.41 vs 0.61 at 15 TLS iterations/level;
      15 is the benchmark setting.
    - **GPU is the default now** (v2.1.1, `-platf 1` via
      `docker run --gpus all`): ~20–25 s/case on the full volumes vs ~4 min
      CPU (120–440 s, `-omp` thread-bound) for the conda 1.5.69 build.
      `--niftyreg-gpu -1` forces the CPU platform (`-platf 0`).
    - **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
    3D transform is `EulerTransform` (6-DOF, `Rigid3DTransform` is gone from
    the component set), the metric key is `Metric` (not `ImageMetric`), the
    random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x
    names and is not compatible with elastix ≤ 4.x parameter files.
  - `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the
    start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the
    optimizer scales. An **output directory** is still set
    (`SetOutputDirectory`) for the C++ filter to run; the adapter uses a
    scratch dir under `data/m6reg/elastix/itkelastix/`.
  - `WriteResultImage` must be the string `"false"` (Python `False`); the
    default (absent) makes the filter raise "Could not read result image"
    after a *successful* registration.
 - With fixed=T1c/moving=CT the combination transform is T1c→CT — the
     contract direction — so `fix_orientation` keeps it as-is
     (`auto_inverted: False`). The result is recovered by point-mapping
     (`GetCombinationTransform().TransformPoint`) into an affine — independent
     of elastix's Euler-angle parameterization — and wrapped in a SimpleITK
     `AffineTransform(3)`.
  - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
    (ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their
    defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32).
 - **IMPACT specifics** (`impact-reg-konfai` 1.8.6, KonfAI runtime,
  elastix-impact 5.2, GPU; preset `Elastix_IMPACT_Static`):
    - **Runs on M6 with the T1c-fixed direction + relaxed sample-validity
      gate.** Full volumes, no mask, no crop / pad, fixed=T1c/moving=CT (the
      shared rule; the ImpactReg MR_CT convention coincides with it): the
      preset's output field is on the fixed grid, so that direction is the
      contract transform (T1c→CT) directly, and the metric's samples (drawn
      from the tight T1c crop) map mostly inside the large CT FOV. The package's generated map still defaults to
      `RequiredRatioOfValidSamples 0.25`, so the adapter overrides it to 0.1
      (a genuinely diverging transform, <10% in-buffer, still aborts). In the
      old Docker era (fixed=CT) the gate aborted the no-mask M6 pairs until
      the same 0.1 relaxation. The transform contract holds; `_orient_impact`
      decides the direction as usual (completed runs report
      `auto_inverted: False` — the MR_CT field is already the contract).
    - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) is
      applied as an override (the preset's GeometricalCenter start is too
      crude for the head+neck FOV mismatch); no manual translation is used.
    - **Deformable, not rigid.** The preset's native transform is a
      RecursiveBSpline (FinalGridSpacing 14 mm, 4 scales) regularized by the
      bending-energy metric; the old pipeline forced a rigid `EulerTransform`
      in place of it. The benchmark contract (warp the CT concentric with the
      MR head) is transform-agnostic, so this is a native-config comparison
      point like every other row.
    - **History: single-init convergence.** The semantic landscape has more
      local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer
      nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0;
      the pre-mask head-box flow got it to mid 8.3 / std 13.2). The no-mask
      MRSegmentator runs showed the same pattern (0.635 MI on `m6_22IRRKU6`
      at 2000/5000 but mid 26.1 / std 28.6 there; accepted on `m6_2BBN5XA4`,
      mid 15.1 / std 10.7). The single CenterOfGravity start on the ~30°
      tilt pairs remains IMPACT's limiting factor — a multi-start (or
      `--tta` flips through the preset app) would be the next step.
    - **Stochastic.** The Jacobian preset (`Elastix_IMPACT_Jacobian`)
      resamples patches through the network each iteration, so results vary
      run to run; the default Static preset re-samples coordinates each
      iteration but the feature maps are deterministic given the pose (see
      Stochasticity).
    - **The elastix shim (this host).** The stock binary needs glibc ≥ 2.35
      + LibTorch 2.8.0; the host (glibc 2.31, torch 2.14) runs it through
      `~/.cache/konfai/elastix-impact-shim` (built by
      `scripts/setup_impact_konfai.sh`, selected via `KONFAI_ELASTIX_DIR`).
      Same glibc story as the niftyreg v2.x binaries.
    - **Models / presets are not in the repo.** They auto-download from HF
      (`VBoussot/ImpactReg`, `VBoussot/impact-torchscript-models`) into
      `~/.cache/huggingface` on first use; the old `data/impact_models/`
      cache is obsolete.
    - **CPU fallback:** `--impact-gpu -1` (no `--gpu` flag; the engine runs
      the feature models on CPU with mixed precision off — much slower).
--- a/REGISTRATION_BENCHMARK_DEFORMABLE.md
+++ b/REGISTRATION_BENCHMARK_DEFORMABLE.md
@ -0,0 +1,456 @@
 # Registration Pipeline Benchmark — CT → T1c (M6) — Deformable
 Companion to [REGISTRATION_BENCHMARK_RIGID.md](REGISTRATION_BENCHMARK_RIGID.md):
 the **deformable** variant of the same six-pipeline benchmark. The rigid
 benchmark asked "which toolkit aligns the head best with 6 DOF?"; this one
 asks "which toolkit's native deformable workflow aligns the head (and the
 tissue) best?" on the same 10 M6 CT→T1c test cases, the same unified
 contract, and the same metrics.
 Run date: 2026-10-02 (main run 05:43–07:52 local; serial `synthmorph`
 re-run 09:58–10:07 local; accelerated `simpleitk` re-run 15:31–16:26
 local). Follow-up ANTs-variant runs 2026-10-03 (08:44–09:12 local: five
 stock antspyx presets — direct `SyNRA` / `SyN`, `Repro[sr|br]`,
 `QuickRepro[sr|br]`; 11:00–11:26 local: the legacy two-call
 `Rigid → SyNOnly` workflow; 12:44–12:49 local: re-run of the stock `SyN`
 preset as `ants_syn`, replacing the initial `SyNOnly` row) on the same
 10 cases (see the *ANTs variants* section). 10 M6 test cases, 3× NVIDIA A100 80GB PCIe, 96-core host, conda
 env `longitudinal-313`. Main run: 6 workers (24 threads each, 2 workers
 per GPU). `synthmorph` was re-run serially (1 worker, one container per
 GPU at a time) because the concurrent run lost 6/9 cases to cuDNN
 scratch-space allocation failures (see the SynthMorph note under
 Interpretation); the re-run registered 9/9, matching the concurrent
 run's scores to 4 decimals on the 3 cases it completed. `simpleitk`
 was re-run after its spline grid was accelerated from 16 mm to 24 mm
 (see the simpleitk row of the pipeline table; the original 16 mm results
 are archived as `simpleitk_16mm_legacy.json`).
 ## Unified contract
 Identical to the rigid benchmark: the returned transform `T` must satisfy
 `sitk.Resample(ct, t1c, T, ...)` = CT on the T1c grid, i.e. **T maps T1c →
 CT coordinates**. fixed = T1c, moving = CT, full volumes (no head box, no
 mask, no crop, no pad), each toolkit's native automatic transform
 initialization, CT clipped to [0, 80] HU only for the pipelines where the
 rigid A/B showed clipping helps (niftyreg, synthmorph).
 One deformable-specific change: every pipeline's output (elastix BSpline
 combination, ANTs warp + affine composite, NiftyReg affine + F3D spline,
 SimpleITK Euler + BSpline, SynthMorph joint field, IMPACT RecursiveBSpline)
 is **rasterised to a SimpleITK `DisplacementFieldTransform` on the T1c
 grid** before scoring, so all six pipelines return the same transform type
 and the shared scoring / GTV warping / QA code is used verbatim.
 ## The six pipelines (deformable configurations)
 Each pipeline runs its toolkit's **native deformable workflow** — a linear
 init stage + a deformable stage, or a single deformable stage where that is
 the toolkit's canonical setup:
 | pipeline | deformable configuration | metric |
 |---|---|---|
 | **antspyx** | `type_of_transform='SyN'` (ANTs' native deformable): Affine stage (2100×1200×1200×0, shrink 4×2×2×1) + SyN[0.2,3,0] deformable stage (40×20×0, shrink 4×2×1), center-of-mass auto init | Mattes MI (32 bins), both stages |
 | **elastix** | two-stage (the toolkit's canonical deformable workflow): EulerTransform rigid stage (2000 iters × 3 resolutions, CenterOfGravity auto init) → BSplineTransform stage (16 mm grid, spline order 3, 2000 iters × 3 resolutions) initialized from stage 1 via `-t0` | Mattes MI (32 bins), 5000 samples/level |
 | **impact** | the ImpactReg `MR_CT_TS` **reference map in its native deformable form**: RecursiveBSplineTransform (14 mm final grid, order 3) driven by the MIND + TotalSegmentator static feature combination (4 per-level scales 400/300/200/200 iters, 2000 samples) with the reference's three-metric combination: IMPACT (w 1.0) + Mattes MI (w 0.4) + TransformBendingEnergyPenalty (w 10). Two M6 adaptations, same as the rigid benchmark: CenterOfGravity auto init (the reference's GeometricalCenter start leaves the heads ~100 mm apart and aborts) and RequiredRatioOfValidSamples 0.1 | IMPACT semantic (L1 + Dice) + Mattes MI + bending energy |
 | **niftyreg** | reg_aladin rigid init (3-level pyramid, 15 TLS iters/level) → reg_f3d non-rigid (F3D B-spline Free-Form deformation, Modat et al. 2010; NMI metric, symmetric bending/linear regularisation, 3-level pyramid, 150 iters at the final level, 8 mm final grid) | NCC block matching (init) + NMI (F3D) |
 | **simpleitk** | two-stage (the toolkit's canonical deformable demo): Euler3D stage (2000 iters/level, pyramid [4,2,1]) → BSplineTransform3D stage (24 mm grid, 50 iters/level on pyramid [4,2]), the BSpline's control points initialised to reproduce the Euler stage exactly. The 24 mm grid (5577 parameters) is the accelerated setting (the 2026-10-02 original used 16 mm / 13005 parameters, ~61 min/case under the harness): this build's spline per-iteration cost scales ~linearly with the control-point count (~2.8 s/iter at the shrink-4 level and ~20 s/iter at shrink-2 for 24 mm vs ~6.4 / ~45.8 s for 16 mm — independent of the metric sample count, which is fixed in this SimpleITK build; LBFGS2 was tested and fails to move the transform on this noisy MI), and the full re-run keeps the 16 mm quality within 0.003–0.029 MI per case (median 0.698 vs 0.705, ImpactLoss median unchanged) while retaining simpleitk's `m6_26TYQI3C` MI/ImpactLoss win — see the acceleration note under Interpretation. A full-resolution polish pass (10 iters at shrink-1) was tested and dropped the achieved MI (0.6891 at 32 mm), so the pyramid stays [4,2] | Mattes MI (32 bins), both stages |
 | **synthmorph** | `joint` model (affine + deformable in one symmetric network pass — the toolkit's canonical deformable; the `deformable` model requires a pre-existing field and would bias the comparison) | anatomy-aware contrast-invariant network (pretrained) |
 ## Test cases
 First 10 M6 labelled cases of `m6.discover()` (all with a GTV label):
 `m6_22GD5VL2_20230724`, `m6_22IRRKU6_20240115`, `m6_23XIMO4M_20230222`,
 `m6_24V62QQM_20230223`, `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820`,
 `m6_26TYQI3C_20221027`, `m6_2BBN5XA4_20250411`, `m6_2C4TKI6H_20240506`,
 `m6_2CDNN7BH_20240102`.
 (`m6_22GD5VL2_20230724` has an incomplete CT head — a known data issue,
 not a pipeline failure; it is scored as ct_incomplete in every pipeline,
 as in the rigid benchmark.)
 (Ground-truth reference transforms are **not** part of this benchmark —
 scoring is by achieved image similarity only, as in the rigid benchmark.)
 ## Metrics (identical for every pipeline)
 - **Mattes MI** (nats) between the warped CT and the T1c over the MR head
  (largest connected component of the T1c mask) — achieved alignment, higher
  = better.
 - **IMPACT loss** (vboussot/ImpactLoss; pretrained TotalSegmentator
  `TS/M730` deep features — weighted L1 over the two deepest layers, full
  volumes; the same scorer as the rigid benchmark) between the warped CT
  and the T1c — semantic-anatomy similarity, lower = better.
 - **Registration wall time** (seconds, excludes I/O + scoring).
 Per-case JSONs additionally retain the production skull-ring QA (gap +
 concentricity, and the acceptance verdict) for reference.
 ### Result
 Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ T1c
 over the MR head; higher = better) **and** the **ImpactLoss** semantic score
 (warped-CT ~ T1c over the full volumes; lower = better), on the same
 unclipped scoring volumes as the rigid benchmark, so all rows are comparable
 on the same scales and against the rigid numbers. Full volumes, no mask /
 crop / pad, fixed = T1c / moving = CT (shared rule). `niftyreg` and
 `synthmorph` clip the CT (the moving image) to [0, 80] HU *before
 registration* (per the rigid benchmark's clipping A/B); `antspyx`,
 `elastix`, `impact` and `simpleitk` register the CT as-is. All six
 transforms are deformable (each toolkit's native workflow) and are
 rasterised to T1c-grid displacement fields before scoring.
 All six pipelines registered 9/9 cases (the 10th is the known
 incomplete-CT data issue, unscored in every pipeline). **Bold** = best
 for the column (MI higher is better, ImpactLoss lower is better).
 | pipeline | toolkit | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) |
 |---|---|---|---|---|---|
 | **antspyx** | ANTsPy 0.6.3 | 9 | 31.7 | 0.158 / 0.156 | **0.723 / 0.729** |
 | **elastix** | ITK-Elastix 0.25.4 | 9 | 45.8 | 0.164 / 0.164 | 0.622 / 0.627 |
 | **impact** | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT, deformable MIND+TS reference map) | 9 | 273.2 | 0.159 / 0.159 | 0.553 / 0.547 |
 | **niftyreg** | NiftyReg 2.1.1 (reg_aladin + reg_f3d, Docker + CUDA) | 9 | 38.1 | 0.161 / 0.161 | 0.666 / 0.686 |
 | **simpleitk** | SimpleITK 2.5.6 | 9 | 1369.9 | 0.158 / 0.159 | 0.698 / 0.704 |
 | **synthmorph** | Docker `freesurfer/synthmorph` (TensorFlow, pretrained) | 9 | 49.1 | **0.156 / 0.156** | 0.627 / 0.640 |
 Per-case MI (nats, unclipped score; **bold** = best for the case). The
 `synthmorph` column is from the serial re-run (see the SynthMorph note
 under Interpretation):
 | 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.679** | 0.580 | 0.548 | 0.637 | 0.654 | 0.620 |
 | m6_23XIMO4M_20230222 | **0.723** | 0.622 | 0.573 | 0.711 | 0.693 | 0.618 |
 | m6_24V62QQM_20230223 | **0.789** | 0.570 | 0.635 | 0.729 | 0.738 | 0.701 |
 | m6_24V62QQM_20231017 | **0.720** | 0.622 | 0.580 | 0.665 | 0.692 | 0.628 |
 | m6_24V62QQM_20240820 | **0.738** | 0.674 | 0.462 | 0.707 | 0.700 | 0.666 |
 | m6_26TYQI3C_20221027 | 0.799 | 0.694 | 0.537 | 0.793 | **0.801** | 0.709 |
 | m6_2BBN5XA4_20250411 | **0.666** | 0.585 | 0.436 | 0.612 | 0.641 | 0.569 |
 | m6_2C4TKI6H_20240506 | **0.719** | 0.635 | 0.553 | 0.656 | 0.698 | 0.627 |
 | m6_2CDNN7BH_20240102 | **0.726** | 0.660 | 0.602 | 0.666 | 0.718 | 0.624 |
 Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest)
 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.1541 | 0.1595 | 0.1532 | 0.1530 | 0.1572 | **0.1512** |
 | m6_23XIMO4M_20230222 | 0.1581 | 0.1657 | **0.1567** | 0.1605 | 0.1603 | 0.1573 |
 | m6_24V62QQM_20230223 | 0.1473 | 0.1634 | 0.1497 | 0.1541 | 0.1502 | **0.1470** |
 | m6_24V62QQM_20231017 | **0.1518** | 0.1591 | 0.1547 | 0.1620 | 0.1548 | **0.1518** |
 | m6_24V62QQM_20240820 | **0.1545** | 0.1646 | 0.1694 | 0.1706 | 0.1584 | 0.1562 |
 | m6_26TYQI3C_20221027 | 0.1580 | 0.1638 | 0.1656 | 0.1672 | **0.1562** | 0.1602 |
 | m6_2BBN5XA4_20250411 | 0.1624 | 0.1663 | 0.1669 | 0.1648 | 0.1631 | **0.1616** |
 | m6_2C4TKI6H_20240506 | 0.1594 | 0.1640 | 0.1585 | **0.1552** | 0.1626 | 0.1564 |
 | m6_2CDNN7BH_20240102 | 0.1590 | 0.1662 | 0.1597 | 0.1605 | 0.1633 | **0.1583** |
 ### Interpretation
 - **`antspyx` (Affine + SyN) leads on achieved MI** (median 0.723; best on
  8 of 9 cases and within 0.002 of best on the other —
  `m6_26TYQI3C`, 0.799 vs `simpleitk`'s 0.801) at the best median
  ImpactLoss of the gradient-based tools (0.158, within 0.002 of
  `synthmorph`'s 0.156), and it is the fastest full-workflow pipeline
  (31.7 s/case) — the best time/quality trade of the six.
 - **`simpleitk` is second on MI** (0.698; it keeps its per-case win on
  `m6_26TYQI3C` — 0.801, still the highest MI of the whole benchmark,
  though now only 0.001 ahead of `antspyx` — and cedes its
  `m6_2CDNN7BH` win to `antspyx`, 0.718 vs 0.726) at a cost of ~23
  min/case under the harness — ~43× the fastest pipeline after the
  spline-grid acceleration below (this build's BSpline per-iteration
  cost scales ~linearly with the control-point count, independent of the
  metric sample count — see the pipeline table).
 - **`simpleitk` acceleration (2026-10-02):** the original setting (16 mm
  spline grid, 13005 parameters) ran 3658.9 s/case (median) under the
  harness. The spline stage's per-iteration cost in this build scales
  ~linearly with the control-point count (16 mm: ~6.4 s/iter at the
  shrink-4 level, ~45.8 s at shrink-2; 24 mm: ~2.8 / ~20 s; 32 mm:
  ~1.2 / ~9.6 s), and the metric sample count is fixed in this SimpleITK
  build (`SetMetricAsMattesMutualInformation` exposes only the histogram
  bins); LBFGS2 was tested and its line search rejects every step on
  this noisy MI (the transform does not move). Switching to the **24 mm
  grid (5577 parameters) cut the wall to 1369.9 s median (2.7×)** at a
  cost of −0.003…−0.029 MI on 8 of 9 cases (+0.017 on `m6_2C4TKI6H`;
  median 0.698 vs 0.705), with the median ImpactLoss unchanged (0.158)
  and `simpleitk`'s `m6_26TYQI3C` MI/ImpactLoss wins intact. A 32 mm
  grid (a further ~1.9×) was rejected: it drops below `antspyx`'s 0.7994
  on `m6_26TYQI3C` (0.7902, and 0.7963 even with 100 fine-level iters).
  A full-resolution polish pass (10 iters at shrink-1) was also tested
  and dropped the achieved MI (0.6891 at 32 mm). The original 16 mm
  results are archived as `simpleitk_16mm_legacy.json`.
 - **`niftyreg`'s `reg_f3d` stage delivers the largest deformable gain of
  any pipeline** (0.292 → 0.666 median MI, +0.374): it recovers the two
  gross-offset cases where the rigid block-matching locked onto local
  minima (`m6_24V62QQM_20240820` 0.118 → 0.707, `m6_2BBN5XA4` 0.069 →
  0.612), lifts the toolkit to third on MI, and it is the most
  concentric of the six (median skull-ring gap 5.5 mm mid / 3.8 mm std,
  6/9 accepted).
 - **The deformable stage does not uniformly pay off.** Against the rigid
  benchmark's per-case MI it improves `antspyx` (0.652 → 0.723 median),
  `simpleitk` (0.647 → 0.698, accelerated 24 mm setting) and `niftyreg`
  (0.292 → 0.666), but
  **degrades** `elastix` (0.646 → 0.622; mixed per case, worst regression
  `m6_24V62QQM_20230223` 0.698 → 0.570) and `impact` (0.652 → 0.553 —
  lower than its own rigid init on **all 9 cases**). The two degradations
  are the coarse-grid BSpline configurations (14–16 mm) with heavy
  bending regularisation — for `impact` the reference map's
  TransformBendingEnergyPenalty at weight 10 — so the deformable stage
  smooths rather than refines, and (for `impact`) its IMPACT-dominant
  objective is simply not the achieved-MI objective.
 - **`impact`'s deformable reference map is the weakest on the intensity
  scale** (MI 0.553) even though its *semantic* score is competitive
  (ImpactLoss 0.159, within 0.003 of the best median) — a
  metric-objective/achieved-MI divergence the rigid benchmark did not
  show (there the two scores agreed, with `impact` in the top tier on
  both). At 273 s/case it is the worst time/quality trade of the six.
 - **`synthmorph` (serial re-run, 9/9) splits the leaderboards: best on
  the semantic scale, mid-pack on the intensity scale.** Median
  ImpactLoss 0.156 — the best of the six — with the per-case best on 4 of
  9 (plus a 4-decimal tie with `antspyx` on `m6_24V62QQM_20231017`),
  while its MI (0.627) is fourth, statistically tied with `elastix`
  (0.622). The pretrained contrast-invariant network optimises learned
  anatomy rather than raw CT/MR intensity synchrony, so it wins where
  anatomy matters and loses where the intensity histogram matters. Fast
  (49.1 s/case) and the read-back cross-check against the tool's own
  resampling is NCC 1.0 on all 9 cases.
 - **The two score families discriminate differently here than in the
  rigid benchmark.** The ImpactLoss scale is nearly saturated (medians
  0.156–0.164, ~5 % spread; per-case bests spread across five tools,
  `elastix` with 0/9), whereas the MI scale spans 0.553–0.723
  (0.17 nats absolute) and the per-case bests go 8/9 `antspyx`, 1/9
 `simpleitk`, 0/9 the rest. The deformable stage's main effect is fine
   intensity alignment on top of the gross anatomy the rigid stages
   already place — exactly what MI measures and ImpactLoss is less
   sensitive to. For ranking the deformable workflows, achieved MI is
   the more discriminative score.
 - **`synthmorph` infrastructure note.** In the main concurrent run (6
  workers, 2 containers per GPU) 6 of 9 cases failed with `CUDNN failed
  to allocate the scratch space ... [Op:Conv3D]` inside the joint model's
  256-channel 128³ hyperconvolution. Each case is a fresh
  `docker run --rm` container, so the failures are cuDNN scratch-space
  contention when two SynthMorph containers share one A100, not a
  pipeline fault: the serial re-run registered 9/9, reproduced the
  concurrent run's MI/ImpactLoss to 4 decimals on its 3 completed cases
  (deterministic model), and ran ~17 s/case faster without the
  contention (49.1 s vs 65.7 s median). The concurrent-run JSON is archived as
  `synthmorph_concurrent_run.json`.
 - **Production skull-ring QA (reference only, kept in the per-case
  JSONs):** acceptance 6/9 `niftyreg`, 5/9 `elastix` / `synthmorph`,
  4/9 `antspyx` / `simpleitk`, 2/9 `impact`. Concentricity: `niftyreg`
  best (median gap 5.5 mm mid / 3.8 mm std), `simpleitk` (19.1 mm std,
  accelerated setting) and `impact` (18.9 mm) worst on the rotational
  term. The same weak correlation with perceived alignment noted in the
  rigid benchmark, so it is reported for reference only.
 ## ANTs variants: direct SyNRA / SyN and the Repro / QuickRepro workflows
 The current `antspyx` row above is antspyx's stock **`SyN` preset** (Affine
 2100×1200×1200×0 + SyN[0.2,3,0] 40×20×0, Mattes MI 32, center-of-mass init).
 As a follow-up (2026-10-03, 6 workers as the main run) seven further
 ANTs configurations were benchmarked on the same 10 cases and the same
 contract — the question: is there a better stock ANTs deformable workflow
 for CT→T1c than the `SyN` preset?
 | pipeline | antspyx preset | stages (all COM init unless noted) |
 |---|---|---|
 | **ants_rigid_synonly** | legacy two-call (`Rigid` + `SyNOnly`) | Rigid[0.25] (Mattes 32, 500×500×200, 4×2×1, 3×2×0) → SyNOnly SyN[0.2,3,0] (Mattes 32, 100×100×0, 4×2×1, 2×1×0) initialized from the Rigid stage; `Composite.h5` read-back |
 | **ants_synra** | `SyNRA` | Rigid[0.25] + Affine[0.25] (each Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1) |
 | **ants_syn** | `SyN` | Affine[0.25] (Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1); COM init — the stock preset, identical configuration to the existing `antspyx` row |
 | **ants_synreprosr** | `antsRegistrationSyNRepro[sr]` | Rigid[0.1] (GC 1/1/Regular/0.25, 1000×500×250×100, 8×4×2×1, 3×2×1×0vox) + SyN[0.2,3,0] (CC 1/2, 50×20, 2×1, 1×0vox) |
 | **ants_synreprobr** | `antsRegistrationSyNRepro[br]` | the `Repro` Rigid[0.1] stage + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×20, 2×1, 1×0vox) |
 | **ants_synquickreprosr** | `antsRegistrationSyNQuickRepro[sr]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + SyN[0.2,3,0] (CC 1/2, 50×0, 2×1) |
 | **ants_synquickreprobr** | `antsRegistrationSyNQuickRepro[br]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×0, 2×1) |
 `sr`/`br` select the deformable stage (SyN vs BSplineSyN, antspyx preset
 naming — not a brain mask); these are the same stage configurations the
 `antsRegistrationSyNRepro*.sh` / `antsRegistrationSyNQuickRepro*.sh`
 workflows use. `ants_rigid_synonly` is the two-call workflow that an early
 version of this benchmark used: a Rigid stage, then a
 *deeper* SyNOnly stage (100 iterations at every level incl. full
 resolution, vs the `SyN` preset's 40×20×0) initialized from the Rigid
 stage's transform. Full volumes, no CT clip, each toolkit preset's own
 metrics (GC/CC/Mattes as above). As for the other rows, the stage files
 are rasterised to the T1c-grid displacement field and scored with the
 shared code; the composite read-back was cross-checked per case against
 ANTs' own resampling (**NCC 1.0 on all 63 scored cases**,
 `own_warp_ncc_head` in the JSONs). One read-back subtlety: when ANTs runs Rigid + Affine as separate
 linear stages (`SyNRA`), it writes the merged linear stage as a *binary*
 `0GenericAffine.mat` that SimpleITK cannot parse; the pipeline recovers it
 by point-probing the transform (exact for affine maps, verified to
 2×10⁻⁴ mm).
 All six variants registered 9/9 (the 10th is the known incomplete-CT data
 issue). **Bold** = best for the column; `antspyx` is the existing row,
 re-listed for reference.
 | pipeline | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) |
 |---|---|---|---|---|
 | **ants_rigid_synonly** | 9 | 208.4 | 0.156 / 0.155 | **0.753 / 0.763** |
 | **antspyx** (existing) | 9 | 31.7 | 0.158 / 0.156 | 0.723 / 0.729 |
 | **ants_synra** | 9 | 43.6 | 0.158 / 0.156 | 0.722 / 0.729 |
 | **ants_syn** | 9 | 33.8 | 0.158 / 0.156 | 0.722 / 0.729 |
 | **ants_synreprosr** | 9 | 294.7 | **0.151 / 0.153** | 0.530 / 0.529 |
 | **ants_synreprobr** | 9 | 370.3 | 0.156 / 0.157 | 0.541 / 0.537 |
 | **ants_synquickreprosr** | 9 | 83.7 | 0.154 / 0.156 | 0.542 / 0.534 |
 | **ants_synquickreprobr** | 9 | 96.0 | 0.158 / 0.158 | 0.549 / 0.538 |
 Per-case MI (nats, unclipped score; **bold** = best for the case):
 | case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr |
 |---|---|---|---|---|---|---|---|---|
 | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
 | m6_22IRRKU6_20240115 | 0.679 | **0.741** | 0.681 | 0.681 | 0.623 | 0.646 | 0.630 | 0.627 |
 | m6_23XIMO4M_20230222 | 0.723 | **0.760** | 0.722 | 0.722 | 0.530 | 0.544 | 0.529 | 0.530 |
 | m6_24V62QQM_20230223 | 0.789 | **0.824** | 0.789 | 0.787 | 0.520 | 0.528 | 0.542 | 0.549 |
 | m6_24V62QQM_20231017 | **0.720** | 0.716 | **0.720** | **0.720** | 0.473 | 0.475 | 0.469 | 0.462 |
 | m6_24V62QQM_20240820 | 0.738 | **0.753** | 0.738 | 0.738 | 0.524 | 0.525 | 0.517 | 0.524 |
 | m6_26TYQI3C_20221027 | 0.799 | **0.857** | 0.798 | 0.799 | 0.557 | 0.541 | 0.599 | 0.632 |
 | m6_2BBN5XA4_20250411 | 0.665 | **0.696** | 0.665 | 0.666 | 0.230 | 0.303 | 0.235 | 0.297 |
 | m6_2C4TKI6H_20240506 | 0.719 | **0.751** | 0.719 | 0.719 | 0.672 | 0.649 | 0.654 | 0.621 |
 | m6_2CDNN7BH_20240102 | 0.726 | **0.765** | 0.726 | 0.725 | 0.629 | 0.622 | 0.631 | 0.599 |
 Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) for
 the case):
 | case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr |
 |---|---|---|---|---|---|---|---|---|
 | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
 | m6_22IRRKU6_20240115 | 0.1541 | 0.1530 | 0.1542 | 0.1540 | **0.1442** | 0.1482 | 0.1474 | 0.1493 |
 | m6_23XIMO4M_20230222 | 0.1581 | 0.1565 | 0.1581 | 0.1581 | **0.1505** | 0.1549 | 0.1535 | 0.1585 |
 | m6_24V62QQM_20230223 | 0.1473 | **0.1464** | 0.1472 | 0.1472 | 0.1515 | 0.1553 | 0.1535 | 0.1560 |
 | m6_24V62QQM_20231017 | 0.1518 | **0.1512** | 0.1516 | 0.1519 | 0.1559 | 0.1575 | 0.1577 | 0.1581 |
 | m6_24V62QQM_20240820 | **0.1545** | **0.1545** | 0.1547 | 0.1546 | 0.1553 | 0.1578 | 0.1586 | 0.1586 |
 | m6_26TYQI3C_20221027 | 0.1580 | 0.1562 | 0.1578 | 0.1579 | **0.1561** | 0.1574 | 0.1602 | 0.1610 |
 | m6_2BBN5XA4_20250411 | 0.1624 | **0.1611** | 0.1624 | 0.1624 | 0.1680 | 0.1695 | 0.1690 | 0.1707 |
 | m6_2C4TKI6H_20240506 | 0.1594 | 0.1582 | 0.1594 | 0.1593 | **0.1490** | 0.1526 | 0.1517 | 0.1545 |
 | m6_2CDNN7BH_20240102 | 0.1590 | 0.1582 | 0.1589 | 0.1589 | **0.1493** | 0.1560 | 0.1538 | 0.1580 |
 ### Interpretation (ANTs variants)
 - **`ants_rigid_synonly` is the best ANTs configuration of the eight** —
  the only one that beats the current `antspyx` row: median MI **0.753**
  (per-case best on 8 of 9; its 0.857 on `m6_26TYQI3C` is the highest MI
  of the whole benchmark, ahead of the previous 0.801) at ImpactLoss
  0.156 (second-best median after `synthmorph`'s 0.156). The gain over
  the `SyN` preset comes from the *deeper* SyN stage: 100 iterations at
  every level **including full resolution** (100×100×0 vs 40×20×0)
  refining the Rigid stage's pose. The cost is 208.4 s/case (~6.6× the
  preset) — the full-resolution SyN iterations dominate.
 - **`ants_synra` is indistinguishable from the current `antspyx` row** on
  quality (median MI 0.722 vs 0.723, ImpactLoss 0.158 vs 0.158; per-case
  deltas within 0.002 nats, ties on 4 of 9) at **+38 % wall time**
  (43.6 s vs 31.7 s): the explicit Rigid stage ahead of the Affine stage
  buys nothing — on these pairs the 12-DOF stage from the center-of-mass
  init already finds the gross pose, so the extra 6-DOF stage just
  re-solves it.
 - **`ants_syn` reproduces the existing `antspyx` row** (median MI 0.722
  vs 0.723, ImpactLoss 0.158, s/reg 33.8 vs 31.7; per-case deltas within
  0.002 nats) — as expected, since it is the *same stock `SyN` preset*
  (Affine + SyN[0.2,3,0], COM init) re-run in this follow-up batch: it
  serves as a consistency check that the shared read-back / scoring path
  is configuration-stable.
 - **The Repro / QuickRepro family is semantically competitive but weaker on
  intensity:** median MI 0.530–0.549 (vs 0.723 for the `SyN` preset),
  driven by the GC-based Rigid stage and the short 2-level SyN/BSplineSyN
  refinement (50×20 or 50×0 iters, 2×1 pyramid, CC 2-bin metric) — a
  workflow tuned for intra-modality T1w-T1w. `ants_synreprosr` still takes
  the **best ImpactLoss median of all thirteen pipelines here** (0.151;
  per-case wins on 4 of 9), so the semantic/intensity split of the main
  table carries over to the stock ANTs workflows.
 - **Time tracks budget, not quality:** the full Repro Rigid stage
  (1000×500×250×100, 8×4×2×1, GC) dominates — `ants_synreprobr` is the
  slowest ANTs configuration (370 s/case) at the same (low) MI as the
  84–96 s `QuickRepro` variants (which drop the 4th rigid level and the
  2nd SyN level, changing the MI by < 0.02). Only `ants_rigid_synonly`
  challenges the `SyN` preset's quality (better MI at 6.6× the time);
  every other variant is worse for more or for the same.
 - **Production skull-ring QA (reference only):** acceptance 4/9
  `ants_syn` / `ants_synra` / `ants_rigid_synonly` (= baseline), 2/9
  `ants_synreprobr`, 1/9 the rest; rotational dispersion worst in the
  `Repro` / `QuickRepro` family (24–29 mm std) vs 14–15 mm for
  `ants_syn` / `ants_synra` / `ants_rigid_synonly` — consistent with
  their alignment.
 ## Transform conventions (deformable-specific)
 Unlike the rigid benchmark (where the elastix-style toolkits came out in the
 opposite direction 8/9 cases and needed the best-of-T vs T⁻¹ orientation
 gate), every deformable output here is consumed in the documented
 fixed→moving = T1c→CT direction; the gate is replaced by an analytic
 read-back plus a skull-ring sanity check. Two non-obvious conventions were
 established empirically (validated against each tool's own resampling, NCC
 ≥ 0.998 on the head):
 - **NiftyReg and surfa (FreeSurfer's NIfTI stack) decode RAS+ NIfTI qforms
  naively, i.e. they work in an x/y-mirrored world.** A displacement delta
  (dx, dy, dz) stored by such a tool is (−dx, −dy, dz) in the true RAS
  world; the field grid (file indices) is unaffected. The read-back flips
  the x/y components of the dense field. (ANTs writes correct RAS fields —
  no flip needed.)
 - **NiftyReg reg_f3d writes only the control-point grid** (outputCPP.nii);
  `reg_transform -disp` evaluates it to a dense field on the reference grid.
  The reg_aladin affine is folded into the control points before
  optimization (reg_f3d initialises the spline with the affine's
  deformation field), so the dense field **is** the full transform — no
  separate composition step. (The `-comp` subcommand was tested and does
  not reproduce reg_f3d's internal transform; it is not used.)
 - **SynthMorph joint**: driven as `register -m joint moving fixed`
  (moving = clipped CT first, fixed = T1c second); the `-t` output is the
  displacement field on the fixed (T1c) grid — the contract transform
  directly. The tool's own `-o` resampling is written as a cross-check
  (per-case `own_resample_ncc_head` in the JSONs).
 - **ANTs SyN**: the two stage outputs (001Warp.nii.gz + 000GenericAffine.mat)
  are composed as warp-then-affine (verified to reproduce ANTs' own
  `warpedmovout` to < 0.001 MI) and rasterised.
 - **Elastix / IMPACT**: the combination transform (auto-init folded in) is
  rasterised — elastix in-process via `itk.TransformToDisplacementFieldFilter`
  (exact, machine precision on affine fields), IMPACT by evaluating the
  BSpline control-point file directly (SimpleITK cannot read elastix's
  BSpline .txt format; scipy's order-3 B-spline = ITK's kernel).
 ## Reproducing
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 cd /mnt/b4/xfr/git26/longitudinal
 # all six deformable pipelines (results go to results/bench_registration_deformable/)
 # the 2026-10-02 main run used --workers 6 --threads 24 (~2 h 9 min wall,
 # dominated by simpleitk's ~111 min phase)
 python scripts/bench_registration.py --deformable --workers 6 --threads 24 --qa-images
 # single pipeline / case subset
 python scripts/bench_registration.py --deformable --pipelines simpleitk impact --n-cases 3
 # the ANTs-variant follow-up (direct SyNRA / SyN + Repro / QuickRepro
 # sr/br; results go to the same results/bench_registration_deformable/ dir,
 # 2026-10-03 run, 6 workers, ~28 min wall; `ants_syn` is the re-run of the
 # stock SyN preset that replaced the initial SyNOnly row)
 python scripts/bench_registration.py --deformable --pipelines "ants_synra ants_syn ants_synreprosr ants_synreprobr ants_synquickreprosr ants_synquickreprobr" --workers 6 --threads 24 --qa-images
 # the legacy two-call Rigid -> SyNOnly(100,100,100) workflow (2026-10-03,
 # 6 workers, ~25 min wall; dominated by the full-resolution SyN iters)
 python scripts/bench_registration.py --deformable --pipelines "ants_rigid_synonly" --workers 6 --threads 24 --qa-images
 # serial SynthMorph re-run (the 2026-10-02 final numbers): one container per
 # GPU at a time; two concurrent containers on one A100 hit cuDNN
 # scratch-space OOM in the joint model's 256-ch 128^3 hyperconvolution
 python scripts/bench_registration.py --deformable --pipelines synthmorph --workers 1 --threads 24 --qa-images
 # combine the per-pipeline JSONs into summary.json / summary.md / plot
 python scripts/bench_registration.py --deformable --merge
 ```
 ## Visual QA images
 `--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to
 `results/bench_registration_deformable/qa/<pipeline>__<case>.png`
 (9 per pipeline; the incomplete-CT case is skipped) — same format as the
 rigid benchmark: a 3×2 grid of the mid-axial / mid-coronal / mid-sagittal
 slices through the head center, **left column** checkerboard blend of the
 warped CT vs the T1c (clean registration = seamless mosaic), **right
 column** additive color overlay (CT red, T1c green; agreement = yellow,
 misregistration = red/green fringing). CT is brain-windowed 0–80 HU, T1c
 2–98 % head window; the title carries the pipeline, case key and the
 achieved Mattes MI.
--- a/REGISTRATION_BENCHMARK_RIGID.md
+++ b/REGISTRATION_BENCHMARK_RIGID.md
@ -0,0 +1,833 @@
 # Registration Pipeline Benchmark — CT → T1c (M6)
 Benchmarks **six cross-modality registration pipelines — all rigid** (6-DOF)
 — that align a head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be
 warped onto the T1c (MR) native grid. This is the registration step of
 `scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT
 Hounsfield vs MR intensity) and the patient is repositioned between the CT and
 MR sessions, so the inter-session head tilt can reach ~30°.
 - Pipelines: [`src/registration_pipelines.py`](src/registration_pipelines.py)
 - Benchmark harness: [`scripts/bench_registration.py`](scripts/bench_registration.py)
 - Results: `results/bench_registration/` (per-pipeline JSON, `summary.json`,
  `summary.md`, `bench_registration.png`)
 ## Unified contract
 Every pipeline returns a SimpleITK transform `T` such that
 ```python
 sitk.Resample(ct, t1c, T, sitk.sitkLinear)              # CT  onto the MR grid
 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
 every pipeline:
 - **No mask, no crop / pad**: registration runs on the full input volumes
  exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
  padding).
 - **The T1c (MR) is always the fixed image** (moving = CT) for every
  pipeline.  (A CT clipped to a non-negative window reads as MR, so the T1c
  slot decides; both-CT pairs keep the given order.)
 - **Every transform is rigid (6-DOF)** — Euler / rigid-matrix parameterizations
  per toolkit.  No affine, no deformable.
 All pipelines also use each toolkit's **native automatic transform
 initialization** (no manual centroid offset): elastix / impact use
 `AutomaticTransformInitialization` (CenterOfGravity), antspyx its default
 center-of-mass alignment, simpleitk starts from identity, 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. With fixed=T1c the samples are drawn from the
 tight T1c head crop (fixed, ~220×220×150 mm) and map mostly inside the large
 CT head-neck FOV (moving, ~300×300×251 mm), so the gate rarely fires; the
 pipeline still sets `RequiredRatioOfValidSamples 0.1` (vs the build default
 0.25) as a safety override for the FOV mismatch — 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
 mask, no crop / pad) and use each toolkit's **native automatic transform
 initialization** (no manual centroid offset): elastix / impact use elastix's
 `AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default
 center-of-mass alignment, simpleitk starts from identity, niftyreg
 initializes from the image centres, and synthmorph needs none (global deep
 method). **All use fixed = T1c, moving = CT** (shared rule above). `niftyreg`
 and `synthmorph` clip the CT (the moving image) 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`    | T1c / CT | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
 | 2 | `elastix`    | T1c / CT | 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** (MIND `R1D2_3D` L1 + TotalSegmentator `M852`/`M850` deepest-layer Dice, per-level weighted, Static) + Mattes MI (w 0.4) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (Docker `vboussot/elastix_impact` + GPU — see IMPACT section) |
 | 4 | `niftyreg`   | T1c / CT | 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`  | T1c / CT | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
 | 6 | `synthmorph` | T1c / CT | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
 Benchmark budget (best-results settings): `elastix` and `simpleitk` run with
 **MaximumNumberOfIterations = 2000** (per resolution / per pyramid level),
 with the sampler and metric left at their defaults (**NumberOfSpatialSamples
 5000**, **NumberOfHistogramBins 32**). `impact` runs the MIND + TotalSegmentator
 **static** feature combination of the ImpactReg `MR_CT_TS` reference parameter
 map (per `VBoussot/ImpactReg`) as a **rigid** `EulerTransform`: 4 levels,
 per-level voxel sizes 6/3/2/2 mm (no image-resampling pyramid), per-level
 **MaximumNumberOfIterations 400/300/200/200**, **NumberOfSpatialSamples 2000**,
 + Mattes MI weight 0.4, CenterOfGravity auto-init,
 `RequiredRatioOfValidSamples 0.1`. `antspyx`, `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
 `--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,
 ~20× slower, the most robust but not used for the head-to-head).
 ## Test cases
 First 10 labeled M6 cases (`--n-cases 10 --sample first`), reproducible. 9 have
 a complete CT; 1 (`m6_22GD5VL2_20230724`) has an incomplete CT head and is
 reported as `incomplete-CT` for every pipeline (a data issue, not a pipeline
 failure).
 ```
 m6_22GD5VL2_20230724  (incomplete CT)
 m6_22IRRKU6_20240115  m6_23XIMO4M_20230222  m6_24V62QQM_20230223
 m6_24V62QQM_20231017  m6_24V62QQM_20240820  m6_26TYQI3C_20221027
 m6_2BBN5XA4_20250411  m6_2C4TKI6H_20240506  m6_2CDNN7BH_20240102
 ```
 ## Metrics (identical for every pipeline)
 Alignment is scored two ways, both as the registration's *achieved* metric read
 off as a standalone score: the intensity-based **Mattes MI** (`mi_ct_mr`) and
 the semantic **ImpactLoss** (`impact_loss`):
 - **`mi_ct_mr`** — mean Mattes MI (nats) between the **CT inverse-transformed
  into MR coordinates** (the warped CT on the T1c grid) and the **original MR**,
  restricted to the MR head and each image percentile-binned independently
  (SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are
  better in-sync = better aligned; it collapses toward 0 when the CT is
  rotated/translated away from the MR.
 - **`impact_loss`** — the IMPACT semantic loss from
  [vboussot/ImpactLoss](https://github.com/vboussot/ImpactLoss) (the metric
  family the `impact` pipeline registers with): a weighted L1 distance between
  the deep feature maps of the **warped CT** and the **original MR** extracted
  by the pretrained TotalSegmentator `TS/M730.pt` TorchScript model (two
  deepest layers, weights 1/1, no resampling), over the **full co-registered
  volumes**. The model standardizes intensity internally (per-image z-score),
  so both volumes go in raw and the score is invariant to affine intensity
  scaling. Lower = the two volumes carry more in-sync anatomy = better
  aligned; it rises when the CT is rotated/translated away from the MR.
 - **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring).
 Tables and plots report `mi_ct_mr` and `impact_loss` (plus registration time).
 The production geometric QA is still computed per case and kept in the JSON
 (it drives dataset acceptance in `02_build_m6_dataset.py`, but its summary
 numbers were found to correlate poorly with perceived alignment, so they are no
 longer reported in the tables/plots):
 - **`skull_ring_gap_mm`** — march rays from the MR head center; the gap between
  the MR head silhouette boundary and the warped-CT skull ring (p97 of positive
  CT) at mid and low levels, plus the angular std. `mid` measures
  translational concentricity; `std` measures rotational concentricity.
 - **`tv_nvox`** — warped GTV volume.
 - **`accepted`** — the production acceptance rule of `02_build_m6_dataset.py`
  (warped GTV non-empty and concentric with the head; 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, no clip)
 Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
 the MR head; higher = better) **and** the **ImpactLoss** semantic score
 (warped-CT ~ MR over the full volumes; lower = better), with the
 **unclipped** CT, so all rows are comparable on the same score scales. Full
 volumes, no mask / crop / pad, fixed = T1c / moving = CT (shared rule).
 `niftyreg` and `synthmorph` clip the CT (the moving image) to [0, 80] HU
 *before registration* (per the CT-clipping A/B below; `clip_ct=None`
 disables it); `antspyx`, `elastix`, `impact` and `simpleitk` register the CT
 as-is. `niftyreg` is the v2.1.1 Docker build on the GPU (`-platf 1`);
 `impact` runs the `vboussot/elastix_impact` Docker image on the GPU (see the
 IMPACT section). All six transforms are rigid (6-DOF).
 | pipeline | toolkit | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) |
 |---|---|---|---|---|---|
 | **antspyx** | ANTsPy 0.63 | 9 | 16.7 | 0.157 / 0.155 | **0.652 / 0.650** |
 | **elastix** | ITK-Elastix 0.25.4 | 9 | 11.7 | 0.157 / 0.156 | 0.646 / 0.638 |
 | **impact** | Docker `vboussot/elastix_impact` (Elastix 5.2 + IMPACT, rigid MIND+TS static config) | 9 | 111.2 | 0.157 / 0.156 | **0.652 / 0.641** |
 | **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 15.4 | 0.166 / 0.170 | 0.292 / 0.338 |
 | **simpleitk** | SimpleITK 2.5.6 | 9 | 10.8 | **0.156 / 0.156** | 0.647 / 0.642 |
 | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 22.2 | 0.167 / 0.170 | 0.277 / 0.328 |
 Per-case MI (nats, unclipped score; **bold** = best for the case). All rows
 are re-runs under the **fixed = T1c / moving = CT** convention (the shared
 rule); the CT-fixed-era numbers and the retired konfai-preset run are in the
 git history / archived JSONs. The `impact` column is the Docker rigid
 MIND+TS static config (the CT-fixed-era Docker run is archived in
 `results/bench_registration/impact_docker_legacy.json`).
 | 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.648** | 0.644 | 0.645 | 0.292 | 0.647 | 0.277 |
 | m6_23XIMO4M_20230222 | **0.655** | 0.652 | 0.654 | 0.632 | 0.654 | 0.616 |
 | m6_24V62QQM_20230223 | 0.703 | 0.698 | 0.700 | 0.260 | **0.704** | 0.256 |
 | m6_24V62QQM_20231017 | **0.602** | 0.590 | 0.596 | 0.493 | 0.597 | 0.484 |
 | m6_24V62QQM_20240820 | **0.652** | 0.646 | 0.652 | 0.118 | 0.638 | 0.111 |
 | m6_26TYQI3C_20221027 | **0.718** | 0.650 | 0.658 | 0.609 | 0.670 | 0.590 |
 | m6_2BBN5XA4_20250411 | 0.577 | 0.577 | 0.572 | 0.069 | **0.579** | 0.061 |
 | m6_2C4TKI6H_20240506 | **0.645** | 0.638 | 0.639 | 0.233 | 0.636 | 0.231 |
 | m6_2CDNN7BH_20240102 | **0.655** | 0.650 | 0.655 | 0.340 | 0.652 | 0.323 |
 Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) 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.1507 | 0.1508 | 0.1510 | 0.1656 | **0.1504** | 0.1664 |
 | m6_23XIMO4M_20230222 | 0.1565 | **0.1563** | 0.1564 | 0.1570 | 0.1564 | 0.1575 |
 | m6_24V62QQM_20230223 | **0.1466** | 0.1470 | 0.1469 | 0.1664 | 0.1467 | 0.1665 |
 | m6_24V62QQM_20231017 | **0.1520** | 0.1539 | 0.1534 | 0.1577 | 0.1532 | 0.1576 |
 | m6_24V62QQM_20240820 | **0.1566** | 0.1571 | 0.1567 | 0.1840 | 0.1585 | 0.1842 |
 | m6_26TYQI3C_20221027 | **0.1554** | 0.1620 | 0.1622 | 0.1623 | 0.1604 | 0.1620 |
 | m6_2BBN5XA4_20250411 | 0.1617 | **0.1612** | 0.1615 | 0.1847 | 0.1616 | 0.1856 |
 | m6_2C4TKI6H_20240506 | 0.1573 | 0.1579 | 0.1580 | 0.1754 | **0.1547** | 0.1757 |
 | m6_2CDNN7BH_20240102 | 0.1591 | **0.1590** | 0.1591 | 0.1734 | 0.1587 | 0.1741 |
 ### Interpretation
 - **The top four are a statistical tie** (median MI 0.646–0.652, median
  ImpactLoss 0.156–0.157): `antspyx`, `impact`, `simpleitk` and `elastix`
  all recover the gross offset on every case, including the two ~30°
  repositioned ones (`m6_2BBN5XA4` 0.572–0.579, `m6_24V62QQM_20240820`
  0.638–0.652). The two alignment scores — an intensity histogram and a deep
  semantic one — agree on the ranking.
 - **IMPACT (Docker, rigid, T1c-fixed) is in the top tier** at median MI
  0.652 (0.6523) — best or within 0.006 of the best on 8 of 9 cases
  (`m6_24V62QQM_20230223` 0.700, `m6_26TYQI3C` 0.658 where it trails
  ANTsPy's 0.718). This is a large step
  up from the legacy CT-fixed Docker config (median 0.467, archived in
  `impact_docker_legacy.json`): the T1c-fixed direction fixes the
  sample-coverage problem the same way it fixed Elastix (the metric draws
  samples from the fixed T1c crop and evaluates them where they map into the
  large CT FOV, instead of the other way round; previously only ~14% of the
  samples were in-buffer at the auto-init). The `m6_2BBN5XA4` ~30° tilt case
  that defeated every earlier IMPACT variant (0.223 CT-fixed Docker, 0.258
  konfai preset) now registers at 0.572. Cost: the slowest of the six
  (~111 s/case, ~10× the intensity tools). Concentricity is comparable to
  the MI tools (4/9 accepted by the production rule).
 - **ImpactLoss corroborates the MI leaderboard from a different metric
  family.** The top four cluster within 0.001 of each other on the median
  (0.1564–0.1571), and on every case the best (lowest) ImpactLoss belongs to
  one of them (antspyx 4/9, elastix 3/9, simpleitk 2/9, impact 0/9);
  `niftyreg` (0.1664) and `synthmorph` (0.1665) are clearly higher, and the
  gap widens on the hard cases (0.184–0.186 for the two tail tools vs
  0.155–0.162 for the top four). A semantic score computed by pretrained
  segmentation features (not the registration's own objective) ranking the
  pipelines the same way the intensity score does is a strong check that the
  leaderboard is not an artifact of the MI scale.
 - **The grossly repositioned cases still separate the leaders from the
   tail.** On `m6_2BBN5XA4` / `m6_24V62QQM_20240820` (the two most
   repositioned pairs — large inter-session translation; the required
   rotations are only ~10° per the references, see the SynthMorph
   investigation below) the top four hold 0.57–0.65, while niftyreg
   (0.069 / 0.118) and synthmorph (0.061 / 0.111) — the single-start
   block-matching and pretrained-MRI methods — lock onto local minima. This
   is exactly why production `02_build_m6_dataset.py` uses a **multi-start
   tilt grid**.
 - **NiftyReg (reg_aladin) is second-from-bottom**: single-start symmetric
  block matching is strong on moderate cases (`m6_23XIMO4M` 0.632;
  `m6_26TYQI3C` 0.609) but falls into the same tilt local minima
  (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). The v2.1.1 Docker build
  on the GPU (`-platf 1`) is ~15 s/case (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.277 with its default CT clip; 0.241 unclipped). 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). The low score is *not* a tilt effect: the required rotations
  are ≤10.6° (the reference tools' own). The investigation (below) shows the
  `rigid` matrix model fails gross-pose estimation case-by-case on the
  repositioned pairs (19–66 mm translation errors; on `m6_24V62QQM_20240820`
  its transform is worse than no registration at all), fails equally on an
  in-domain MR→MR pair with a ~55° relative pose, and is not fixed by the
  `affine` model — while the `joint` (affine + deformable) model recovers
  all three tested failure cases to the best ImpactLoss of the six tools.
 - **Speed** (median wall, 3–4 cases in parallel, 24 threads, T1c-fixed):
  SimpleITK 10.8 s, elastix 11.7 s, NiftyReg 15.4 s (v2.1.1, GPU block
  matching), ANTs 16.7 s, SynthMorph 22.2 s (GPU forward pass), **IMPACT
  111.2 s** (Docker, static feature mode, no backprop, 400/300/200/200
  rigid ASGD iterations).
 - **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
  (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 in the **Static** config
  re-samples spatial coordinates every iteration, but the feature maps are
  deterministic given the pose. The earlier Jacobian-mode Docker experiments
  (archived in `results/bench_registration/impact_*.json`) resampled patches
  through the network each iteration, so their per-case numbers are one
  sample from a distribution — e.g. the MRSegmentator-only Jacobian configs
  reached 0.631 (500/2000) and 0.591 (2000/5000) at ~5×/~35× the Static
  cost, and the MIND+TS/MRSeg static variants 0.417/0.406 (~72/~37 s).
 ## 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.
 ## SynthMorph low-performance investigation (2026-09-30)
 Question: *why* is `synthmorph` weakest (median MI 0.277; 0.06–0.11 on the
 repositioned cases)? The earlier explanation ("trained on brain MRI, defeated
 by the ~30° tilt") was checked and turned out to be wrong in its specifics.
 Everything below was scored with the benchmark's own scoring
 (`bench_registration.py::evaluate`); results in
 `results/synthmorph_investigation/`, reproducible with
 `scripts/investigate_synthmorph.py` (parts A–D, F) and
 `scripts/investigate_synthmorph_mm.py` (part E).
 **0. It is not a pipeline bug.** The `rigid` run is deterministic and
 reproduces the benchmark JSON exactly (e.g. `m6_26TYQI3C` MI 0.5904, gap
 13.1/4.9). The read-back direction is right (8/9 `auto_inverted: false`,
 verified by the `fix_orientation` skull-ring gate), and the matrix is rigid
 (`det 0.999–1.001`, `‖RᵀR−I‖ ≤ 0.0015`). The lone `auto_inverted: true`
 (`m6_2BBN5XA4`) is the gate choosing the *less bad* of two wrong directions —
 a symptom of the failure, not its cause.
 **1. The failure is a gross-pose estimation failure, not a tilt failure.**
 Each case's synthmorph transform was decomposed against the `simpleitk`
 benchmark reference (MI 0.58–0.70 on all 9 cases; its rotations are
 **1.4–10.6°**, not ~30°):
 | case | MI ref / identity / rigid | head-centre err (mm) | rot err (deg) |
 |---|---|---|---|
 | m6_22IRRKU6_20240115 | 0.647 / 0.250 / **0.277** | 28.2 | 6.2 |
 | m6_23XIMO4M_20230222 | 0.654 / 0.246 / **0.616** | 1.9 | 2.2 |
 | m6_24V62QQM_20230223 | 0.704 / 0.276 / **0.256** | 25.6 | 3.7 |
 | m6_24V62QQM_20231017 | 0.597 / 0.331 / **0.484** | 3.8 | 17.7 |
 | m6_24V62QQM_20240820 | 0.638 / 0.274 / **0.111** | 66.3 | 5.9 |
 | m6_26TYQI3C_20221027 | 0.670 / 0.397 / **0.590** | 6.5 | 10.2 |
 | m6_2BBN5XA4_20250411 | 0.579 / 0.046 / **0.061** | 61.1 | 22.1 |
 | m6_2C4TKI6H_20240506 | 0.636 / 0.231 / **0.231** | 40.6 | 3.4 |
 | m6_2CDNN7BH_20240102 | 0.652 / 0.197 / **0.323** | 19.4 | 3.0 |
 - 6/9 cases carry a 19–66 mm head-centre error (translation), and 3/9 a
  10–22° rotation error; the two "clean" failures are translation-only.
 - On `m6_24V62QQM_20240820` the model's transform (MI 0.111) is **worse than
  registering nothing at all** (identity MI 0.274); on `m6_2C4TKI6H` it is
 exactly identity-level (0.231 = 0.231). The net gain over no-registration
   is ≤0.03 on 5/9 cases.
 - The error *vectors* are case-specific, not a systematic bias: the fraction
  of the required network-space translation that goes uncorrected is ≈0 on
  the successes (23XIMO4M 0.07, 20231017 −0.01), ≈1.0 on `m6_24V62QQM_20240820`
  (the model outputs ≈ its network-space identity — the null that aligns the
  two *volume* centres), and −1.7 on `m6_2BBN5XA4` (it moves the head 1.7× the
  required amount, in the wrong direction).
 **2. Cross-modality is not the (only) cause.** An in-domain control — two
 T1c scans of the *same* subject 8 months apart (`m6_24V62QQM`, both tight
 head crops, no CT, no HU issue) — has a true relative pose of **~55–57°**
 (benchmark `simpleitk` mode 55.3°, production multi-start tilt grid 56.6°,
 head centres agreeing to 1.2 mm) and SynthMorph `rigid` fails it equally:
 18.5° rotation error, 18 mm head-centre error, MI 0.361 vs the reference's
 0.689 (ImpactLoss 0.140 vs 0.114).
 **3. The `affine` matrix model does not help.** Same 9 cases, same budget:
 MI within ±0.025 of `rigid` on 8/9 cases (`m6_23XIMO4M`: −0.096, the one case
 where rigid succeeded); the linear part is a near-pure rotation
 (scale 0.98–1.03, anisotropy ≤2%) — the extra 6 DOF are essentially unused.
 **4. The `joint` (affine + deformable) model recovers the failures.** Run
 directly in Docker (`-m joint`, outside the rigid contract; the tool's own
 `-o` resampling warps the CT onto the T1c grid, scored as-is) on the three
 worst cases:
 | case | MI rigid → joint | ImpactLoss joint (vs best intensity tool) |
 |---|---|---|
 | m6_24V62QQM_20230223 | 0.256 → **0.529** | **0.133** (simpleitk 0.147) |
 | m6_24V62QQM_20240820 | 0.111 → **0.481** | **0.132** (antspyx 0.157) |
 | m6_2BBN5XA4_20250411 | 0.061 → **0.433** | **0.129** (antspyx 0.162) |
 The model family's anatomy is fine — its *semantic* alignment on the two
 grossly repositioned cases is the **best of the six tools** — but the matrix
 head (6- or 12-DOF) cannot estimate the gross inter-session pose, and the
 deformable field absorbs what it misses (~40 s/case on the GPU).
 **5. Pipeline-side suspects were ruled out.**
 - **Network-space extent truncation: minor.** Each volume is resampled into
  its own 256 mm isotropic box (centred on its own volume centre, per the
  image code); with the head placed in CT space via the reference transform,
  only 0–6.5% of head voxels fall outside the CT box (0% T1c). Re-centring
  the CT volume on the head (crop to head bbox + 15 mm) changes *nothing*
  (ΔMI ≤ 0.012, Δhead-centre-err ≤ 1.2 mm on all 9 cases) — the model is not
  anchored to the volume centre.
 - **Min-max normalization: a real, quantified cross-modality penalty.** The
  network min-max-normalizes each volume; over the raw CT's full HU range
  (−3024…+3071, air/metal) the brain (0–80 HU) collapses to a near-flat
  **0.496–0.509** band (air sits at ~0.33, the skull at 0.63–1.0), so the
  parenchyma carries almost no internal contrast; the T1c brain spans
  ~0–0.6 (p95 0.26–0.58), so the raw CT reads as a flat-gray blob with a
  bright ring. The [0,80] clip restores full-range brain contrast (this is
  the mechanism behind the CT-clip A/B win, +0.036–0.077 MI). Clipping does
  not make the model domain-typical, which is why it only partially helps.
 **Conclusion.** SynthMorph's matrix models are a *local pose prior*: on the
 in-distribution cases (small relative pose) they estimate the M6 alignment
 well (23XIMO4M 1.9 mm / 2.2°, 20231017 3.8 mm / 17.7°, 26TYQI3C 6.5 mm /
 10.2° residual), but the M6 inter-session repositioning (tens of mm
 translation and, in the same-patient control, a 55° world-frame rotation) is
 outside the range the pretrained model (trained on same-subject MRI pairs)
 reliably inverts, and the output is
 then case-specifically wrong — do-nothing, half-done, or opposite. It is a
 valid data point for how far a pretrained-MRI deep matrix model generalizes
 out-of-domain, and it should not be expected to compete with the intensity
 tools on repositioned CT pairs in rigid mode. If a SynthMorph variant is
 wanted for these data, the `joint` model is the one that works (at the cost
 of leaving the rigid-transform contract).
 Reproducing (needs the `freesurfer/synthmorph` image + GPU; ~15 min):
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 python scripts/investigate_synthmorph.py --parts ABCD --workers 3  # A ref+rigid, B affine, C joint, D network-space
 python scripts/investigate_synthmorph.py --parts F --workers 3     # re-centred CT probe (needs ABCD first)
 python scripts/investigate_synthmorph_mm.py --gpu 0                # in-domain MR->MR control
 ```
 ## Reproducing
 ```bash
 # 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
 # impact (Docker vboussot/elastix_impact; GPU; ~110-120 s/case; see the IMPACT
 # section). --impact-iters / --impact-samples override the reference map's
 # per-level budget (one iters value for all levels)
 python scripts/bench_registration.py --pipelines impact --n-cases 10 \
    --workers 3 --threads 24
 # synthmorph (Docker, needs the pulled image; see the SynthMorph section)
 python scripts/bench_registration.py --pipelines synthmorph --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
 ```
 Other flags: `--sample random --seed 0` (reproducible random 10),
 `--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`,
 `--threads T` (also the OpenMP thread count for `reg_aladin`),
 `--simpleitk-mode {benchmark,production}`,
 `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`,
 `--synthmorph-extent <int>`, `--impact-gpu <int|-1>`, `--impact-iters N`
 (one MaximumNumberOfIterations value for all levels; default the reference
 map's 400/300/200/200), `--impact-samples N` (NumberOfSpatialSamples; default
 2000),
 `--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
 ## Visual QA images (checkerboard + color overlay)
 `--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to
 `<out-dir>/qa/<pipeline>__<case>.png` for eyeballing the alignment.  Each image
 is a 3×2 grid: **rows** are the three orthogonal mid-slices through the head
 center (**mid-axial, mid-coronal, mid-sagittal**), and **columns** are:
 - **left — checkerboard blend:** alternating blocks of the CT inverse-transformed
  into MR coordinates and the original MR; a clean registration gives a seamless
  mosaic, a misregistration shows hard edges at the block boundaries.
 - **right — additive color overlay:** CT in red, MR in green. Structures that
  coincide blend to yellow; rotation/translation error appears as red/green
  fringing at the skull and brain edges.
 The CT is shown in a **brain window (L40/W80 = 0–80 HU)** so the brain
 parenchyma is visible and comparable to the MR (the skull, ~800+ HU, clips to
 white); the MR uses a 2–98 % head window. The coronal and sagittal views are
 oriented superior-up to match the axial view.
 The title carries the pipeline, case key and the achieved Mattes MI, so a
 glance tells you whether the visuals agree with the score. Generate for the
 available pipelines with:
 ```bash
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 python scripts/bench_registration.py \
    --pipelines "simpleitk elastix antspyx niftyreg" --n-cases 10 --qa-images
 python scripts/bench_registration.py \
    --pipelines "synthmorph impact" --n-cases 10 --qa-images
 ```
 (9 images each; the incomplete-CT case is skipped.)
 ## SynthMorph pipeline (Docker) — one-time setup
 The `synthmorph` pipeline runs the [`freesurfer/synthmorph`](https://hub.docker.com/r/freesurfer/synthmorph)
 image: [SynthMorph](https://synthmorph.io) (Hoffmann et al., IEEE TMI 2022), an
 anatomy-aware, acquisition-agnostic **deep** registration (TensorFlow
 `VxmAffineFeatureDetector` with the pretrained `synthmorph.rigid.1.h5` weights).
 No skull-stripping or intensity preprocessing is needed: each volume is
 normalized to an isotropic 1-mm LIA network space internally. It is the
 second Docker-based pipeline and the only deep-learning one.
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
 ```bash
 # 1. Docker (with NVIDIA container toolkit so --gpus works)
 docker --version
 # 2. Pull the image (~11 GB, one-time)
 docker pull freesurfer/synthmorph
 ```
 No extra model assets: the rigid weights ship inside the image
 (`/freesurfer/models`).
 ### How the adapter drives the image
 1. **Rigid matrix model, no init** — `register -m rigid` on GPU when one is
   visible (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces
   CPU). SynthMorph is a global deep method, so it recovers the gross
   inter-session offset on its own and needs **no** initial transform.
 2. **Read-back** — the model estimates first→second, so the pair is passed as
   `(fixed=T1c, moving=CT)` and the `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in
   shared world space) is the T1c→CT contract transform directly; it is
   parsed into a SimpleITK `AffineTransform` and run through the usual
   `fix_orientation` gate.
 3. **GPU sharing** — the container is run as the host uid (clean scratch
   removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
   concurrent benchmark workers can share a GPU.
 ## NiftyReg pipeline — one-time setup
 The `niftyreg` pipeline runs `reg_aladin` from the [NiftyReg](https://github.com/KCL-BMEIS/niftyreg)
 suite (KCL-BMEIS, the actively-maintained NiftyReg): **symmetric block
 matching + trimmed least squares** (Modat et al., J. Med. Img. 1(2) 024003,
 2014). Instead of a histogram MI objective it matches 4-voxel image blocks by
 **normalised cross-correlation** at each level of a 3-level Gaussian pyramid
 and solves a robust (trimmed) least-squares for the rigid transform — a
 genuinely different algorithmic family, and cross-modality tolerant because
 NCC is invariant to intensity scaling.
 **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)
 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
 ```
 ### How the adapter drives the tool
 1. **Rigid 6-DOF, full volumes** — `reg_aladin -ref T1c -flo CT -rigOnly`
   (the T1c is the reference/fixed image, the CT the floating/moving one;
   the CT is clipped to [0,80] HU first). 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 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`; the adapter parses it into
   a SimpleITK `AffineTransform` and runs it through the usual
   `fix_orientation` gate (skull-ring concentricity), which resolves the
   direction empirically rather than trusting the docs. Under
   `-ref T1c -flo CT` the file's affine comes out the other way, so the gate
   inverts it on the M6 cases (8/9 `auto_inverted: True`) to the contract
   transform (T1c→CT).
 ## IMPACT pipeline (Docker `vboussot/elastix_impact`) — one-time setup
 The `impact` pipeline runs the
 [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
 image: Elastix 5.2 with the [IMPACT](https://arxiv.org/abs/2503.24121)
 *semantic* similarity metric ([vboussot/ImpactLoss](https://github.com/vboussot/ImpactLoss))
 as a native Elastix metric. The parameter map is the MIND + TotalSegmentator
 **static** feature combination of the `MR_CT_TS` reference map in the
 [`VBoussot/ImpactReg`](https://huggingface.co/VBoussot/ImpactReg) repo, with a
 **rigid 6-DOF `EulerTransform`** in place of the reference's RecursiveBSpline
 (its bending-energy penalty is zero for a rigid transform, so it is dropped) —
 i.e. the reference metric/budget, rigid-only, on the GPU.
 Prerequisites (check with `RP.availability()`; a missing piece raises
 `PipelineUnavailable` with the exact command):
 ```bash
 # 1. Docker (with NVIDIA container toolkit so --gpus works)
 docker --version
 # 2. Pull the image (~38 GB, one-time)
 docker pull vboussot/elastix_impact
 # 3. Feature models under the cache (the pipeline prints the exact
 #    curl commands if any are missing):
 #    data/impact_models/MIND/R1D2_3D.pt   (~15 KB)
 #    data/impact_models/TS/M852.pt        (~200 MB)
 #    data/impact_models/TS/M850.pt        (~200 MB)
 #    (from HF VBoussot/impact-torchscript-models, or `sh
 #     Data/Models/download_models.sh` from the ImpactLoss repo)
 ```
 ### How the adapter drives the image
 1. **Full volumes, fixed=T1c (MR), moving=CT** — the shared fixed-image rule:
   the metric draws its random spatial samples from the fixed buffer and
   evaluates them where they map inside the moving buffer. With the tight T1c
   head crop (~220×220×150 mm) as fixed, the samples map mostly inside the
   large CT head-neck FOV (~300×300×251 mm), so the
   `RequiredRatioOfValidSamples 0.1` gate (build default 0.25) rarely fires.
 2. **Parameter map** (the rigid `MR_CT_TS` reference config):
   `(Metric "Impact" "AdvancedMattesMutualInformation")` with weights
   1.0 / 0.4; per level (coarse → fine) the **MIND** handcrafted descriptor
   (`MIND/R1D2_3D.pt`, L1, subset 32, layer mask `"1"`) combined with
   **TotalSegmentator** `TS/M852.pt` (levels 0-1) / `TS/M850.pt` (levels
   2-3, Dice, subset 64, deepest layers), layer weights 0.2/0.8, 0.3/0.7,
   0.6/0.4, 0.7/0.3, at ImpactVoxelSize 6/3/2/2 mm (scaled to the moving
   spacing); **Static** mode (whole-image features, no backprop); no
   image-resampling pyramid (all-1 schedules — the scale is carried by
   ImpactVoxelSize, per the IMPACT docs); `RandomCoordinate` sampler; ASGD
   with per-level MaximumNumberOfIterations 400/300/200/200 and
   NumberOfSpatialSamples 2000; `EulerTransform` +
   `AutomaticTransformInitialization` (CenterOfGravity) +
   `AutomaticScalesEstimation`. A copy of the exact map is dumped to
   `data/m6reg/impact/paramdump/` per run for inspection.
 3. **GPU auto-spread** across the host GPUs by pid (`docker run --gpus all`,
   `ImpactGPU` = the worker's device; `--impact-gpu -1` forces CPU and turns
   mixed precision off).
 4. **Read-back** — the transform comes from
   `Out/TransformParameters.0.itk.txt` (an `Euler3DTransform` in the world
   physical frame; the full combination transform, with the auto-init folded
   in by elastix) and is direction-checked by `fix_orientation` against the
   skull ring (with fixed=T1c / moving=CT the gate keeps the file's transform,
   `auto_inverted: False`).
 5. **Intensity normalization happens inside the metric / models** (no
   pipeline-side normalization): the MIND descriptor normalizes each 5-voxel
   neighbourhood to [0, 1] locally (modality-agnostic by construction); the
   segmentation models are standardized to zero mean / unit variance per the
   model card. Pre-normalizing the volumes would be an exact no-op for the
   z-scored models and would only corrupt the 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
  each toolkit's native automatic transform initialization (no manual
  `centroid_offset`) — elastix / impact `AutomaticTransformInitialization`
  (CenterOfGravity), antspyx's center-of-mass alignment, simpleitk identity,
  synthmorph none. Production uses `--simpleitk-mode production`'s
  multi-start tilt grid, which raises SimpleITK's acceptance at ~20× the
  cost; that is out of scope for the head-to-head but is the config the
  dataset actually ships with.
 - `synthmorph` is deterministic (one fixed pretrained model, no random
   subsampling), so its per-case numbers are stable run to run — unlike ANTs /
   IMPACT. Its low MI (median 0.277 with the default CT clip; 0.06–0.11 on the
   repositioned cases) is a genuine gross-pose-estimation weakness of the
   pretrained matrix model (case-by-case, not variance) — see the
   "SynthMorph low-performance investigation" section above.
 - `fix_orientation` (elastix/ants/impact/niftyreg/synthmorph) chooses between
  a transform and its inverse by skull-ring concentricity — a correctness gate
  that makes the result independent of each toolkit's transform-direction
  convention, so the fixed/moving swap is absorbed by the gate rather than
 re-derived per toolkit. Under the shared fixed=T1c / moving=CT rule the raw
   transforms come out in the contract direction for elastix (combination
   transform = fixed→moving), impact (its `TransformParameters.0.itk.txt`,
   the full combination Euler transform) and antspyx (ANTs forward transform)
   and synthmorph (the `-T` LTA, model estimates first→second) — completed
   runs report `auto_inverted: False` for those; `niftyreg`'s `-aff` file
   comes out the other way and the gate inverts it (8/9 `auto_inverted: True`).
 - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU):
    - **Gross-pose estimation weakness (measured, not a bug).** The `rigid`
      matrix model is trained on brain *MRI* (same-subject pairs, small
      relative poses). On these repositioned CT↔T1c pairs it estimates the
      translation wrongly on 6/9 cases (19–66 mm head-centre error; the two
      worst, `m6_2BBN5XA4` / `m6_24V62QQM_20240820`, at 61 / 66 mm) and the
      rotation on 3/9 (10–22°) — and the required rotations are only 1.4–10.6°,
      so it is *not* a "defeated by tilt" story; on `m6_24V62QQM_20240820`
      its transform scores below the unregistered (identity) one. An in-domain
      MR→MR control (same subject, ~55° relative pose) fails equally, and the
      `affine` model changes nothing, while the `joint` (affine + deformable)
      model recovers all three tested failures to the best ImpactLoss of the
      six tools. Full evidence: the "SynthMorph low-performance investigation"
      section.
    - **No centroid init needed.** Being a global deep method it recovers the
      gross offset on its own; the `-i` initial-transform path is not used
      (which also sidesteps this image's float64/float32 bug on the `-i` code
      path).
    - **GPU sharing.** The container runs as the host uid with
      `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
 - **~28 s/case** (GPU: one pretrained forward pass), between SimpleITK
       (~11 s) / ANTs (~17 s) and IMPACT (~115 s, rigid docker config).
 - **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`,
   symmetric block matching + TLS, CPU/OpenMP):
    - **Different algorithmic family.** No histogram MI objective: 4-voxel
      blocks are matched by normalised cross-correlation and the rigid
      transform is solved by a trimmed least squares, iterated per pyramid
      level (`-ln 3 -maxit 15`). NCC's scale invariance is what makes it
      cross-modality tolerant.
    - **Symmetric mode matters.** The default symmetric variant (optimising
      forward and backward transforms concurrently) is the robust one on these
      repositioned pairs; `-noSym` drifts grossly (gap std 45 mm on a clean
      case every other tool passes).
    - **`-maxit 5` (the tool default) is too coarse** — on the clean case
      `m6_26TYQI3C` it reaches MI 0.41 vs 0.61 at 15 TLS iterations/level;
      15 is the benchmark setting.
    - **GPU is the default now** (v2.1.1, `-platf 1` via
      `docker run --gpus all`): ~20–25 s/case on the full volumes vs ~4 min
      CPU (120–440 s, `-omp` thread-bound) for the conda 1.5.69 build.
      `--niftyreg-gpu -1` forces the CPU platform (`-platf 0`).
    - **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
    3D transform is `EulerTransform` (6-DOF, `Rigid3DTransform` is gone from
    the component set), the metric key is `Metric` (not `ImageMetric`), the
    random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x
    names and is not compatible with elastix ≤ 4.x parameter files.
  - `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the
    start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the
    optimizer scales. An **output directory** is still set
    (`SetOutputDirectory`) for the C++ filter to run; the adapter uses a
    scratch dir under `data/m6reg/elastix/itkelastix/`.
  - `WriteResultImage` must be the string `"false"` (Python `False`); the
    default (absent) makes the filter raise "Could not read result image"
    after a *successful* registration.
 - With fixed=T1c/moving=CT the combination transform is T1c→CT — the
     contract direction — so `fix_orientation` keeps it as-is
     (`auto_inverted: False`). The result is recovered by point-mapping
     (`GetCombinationTransform().TransformPoint`) into an affine — independent
     of elastix's Euler-angle parameterization — and wrapped in a SimpleITK
     `AffineTransform(3)`.
  - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
    (ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their
    defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32).
 - **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2 +
  IMPACT metric, GPU; rigid MIND+TS static config of the `MR_CT_TS` reference
  map):
    - **Runs on M6 with the T1c-fixed direction + relaxed sample-validity
      gate.** Full volumes, no mask, no crop / pad, fixed=T1c/moving=CT (the
      shared rule): the metric's samples (drawn from the tight T1c head crop)
      map mostly inside the large CT head-neck FOV, so the gate rarely fires.
      The map still sets `RequiredRatioOfValidSamples 0.1` (build default
      0.25) as a safety override — a genuinely diverging transform, <10%
      in-buffer, still aborts. The combination Euler transform comes out in
      the contract direction (completed runs report `auto_inverted: False`).
    - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity)
      aligns the intensity centroids (no manual translation); the reference
      map's GeometricalCenter start is too crude for the head+neck FOV
      mismatch.
    - **Rigid, like every other pipeline.** The reference map's
      RecursiveBSpline + bending-energy regularizer is replaced by a 6-DOF
      `EulerTransform` (the bending penalty is zero for a rigid transform).
    - **History: single-init convergence.** The semantic landscape has more
      local minima than Mattes-MI. The T1c-fixed direction fixes most of the
      hard cases (the repositioned `m6_22IRRKU6` / `m6_23XIMO4M` now register
      as well as the MI tools do), but the single CenterOfGravity start on the
      ~30° tilt pair `m6_2BBN5XA4` remains a weakness — a multi-start tilt
      grid (as production SimpleITK uses) would be the next step.
    - **Stochastic.** The config is **Static** mode: it re-samples spatial
      coordinates every iteration, but the feature maps are deterministic
      given the pose (see Stochasticity). The older Jacobian-mode configs
      (docker-era experiments, archived in `results/bench_registration/impact_*.json`)
      resampled patches through the network each iteration, so their
      per-case numbers are one sample from a distribution.
    - **Models are under the data root, not the repo.** The TorchScript
      feature models live in `data/impact_models/` (MIND ~15 KB + TS ~200 MB
      each, from HF `VBoussot/impact-torchscript-models`); the pipeline
      raises with the exact `curl` commands if one is missing.
    - **CPU fallback:** `--impact-gpu -1` (`ImpactGPU -1`, mixed precision
      off — much slower).
--- a/scripts/bench_registration.py
+++ b/scripts/bench_registration.py
@ -2,6 +2,9 @@
 Pipelines (see src/registration_pipelines.py for the shared transform contract):
  simpleitk | elastix | antspyx | synthmorph | impact | niftyreg
  + the ANTs variants (deformable only): ants_synra | ants_syn |
    ants_synreprosr | ants_synreprobr | ants_synquickreprosr |
    ants_synquickreprobr
 Each pipeline registers a head-and-neck CT to a tightly-cropped T1c MRI
 (cross-modality, inter-session tilt up to ~30 deg) and the GTV is warped onto
@ -11,9 +14,10 @@ higher = better in-sync = better aligned.  The production geometric QA
 (scripts/02_build_m6_dataset.py) is still recorded per case in the JSON:
  - mi_ct_mr           : Mattes MI (nats), warped-CT ~ MR over the MR head
-                         (the reported alignment score).
+                         (the reported alignment score, higher = better).
-  - tv_in_ct_frac      : fraction of the (dilated) warped GTV that falls inside
+  - impact_loss        : IMPACT semantic loss (vboussot/ImpactLoss; TotalSeg
-                         the warped-CT head tissue (the slab must straddle it).
+                         M730 deep features, weighted L1), warped-CT ~ MR over
                         the full co-registered volumes (lower = better).
  - tv_nvox            : warped GTV volume (>0.5 voxels).
  - skull_ring_gap_mm  : concentricity of the warped-CT skull ring vs the MR
                         head silhouette (mid/low level + angular std).
@ -26,15 +30,16 @@ 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
+REGISTRATION_BENCHMARK_RIGID.md), the benchmark is env-portable: run it under
 whichever conda env has a given toolkit installed, and it records only the
 pipelines that are importable there.  elastix runs in the main env through the
 ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively
 (``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
-``impact`` runs the ``impact-reg-konfai`` package (pip; GPU; one-time elastix
+``impact`` runs the ``vboussot/elastix_impact`` Docker image (GPU; the rigid
-shim setup on hosts too old for the stock binary — see
+MIND + TotalSegmentator static config, needs Docker + the pulled image + the
-REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker image
+model assets — see REGISTRATION_BENCHMARK_RIGID.md); ``synthmorph`` runs in a Docker
-(``freesurfer/synthmorph``; needs Docker + the pulled image).  The `--merge`
+image (``freesurfer/synthmorph``; needs Docker + the pulled image).  The
 `--merge`
 mode then combines the per-pipeline JSON files (accumulated across envs) into
 one summary + plot.
@ -43,7 +48,7 @@ Usage (repo root):
  source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
  python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg impact
-  # synthmorph (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
+  # synthmorph (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK_RIGID.md
  python scripts/bench_registration.py --pipelines synthmorph
  # any pipeline that is not importable in the current env is recorded as
@ -99,6 +104,54 @@ def _m6():
    return RP._load_m6()
 # the ImpactLoss scoring model (the repo's default baseline for MR/CT)
 IMPACT_SCORE_MODEL = "TS/M730.pt"
 _IMPACT_LOSS_STATE = None
 def _impact_loss_scorer():
    """Lazily build the ImpactLoss semantic-loss scorer (once per process) on
    the GPU auto-picked for this worker.  Returns ``(fn, device)`` or None if
    torch / the model are unavailable (the score is then omitted)."""
    global _IMPACT_LOSS_STATE
    if _IMPACT_LOSS_STATE is None:
        try:
            import warnings
            import torch
            from src.impact_loss import IMPACTReg
            ngpu = torch.cuda.device_count() if torch.cuda.is_available() else 0
            if ngpu:
                torch.cuda.set_device(os.getpid() % ngpu)
            with warnings.catch_warnings():
                # torch.jit.load is marked deprecated by upstream torch 2.14;
                # the IMPACT models are TorchScript and jit.load is required
                warnings.simplefilter("ignore", FutureWarning)
                fn = IMPACTReg(IMPACT_SCORE_MODEL, shape=[0, 0, 0], in_channels=1,
                               weights=[1.0, 1.0], model_cache=d("data/impact_models"))
            _IMPACT_LOSS_STATE = (fn, torch.device("cuda" if ngpu else "cpu"))
        except Exception:  # noqa
            _IMPACT_LOSS_STATE = False
    return _IMPACT_LOSS_STATE or None
 def _impact_loss_value(warped_ct, mr):
    """IMPACT semantic loss (lower = more in-sync) between the co-registered
    pair: the CT inverse-transformed into MR coordinates vs the original MR,
    over the full volumes (the IMPACT models standardize intensity inside
    the metric, so both go in raw)."""
    st = _impact_loss_scorer()
    if st is None:
        return None
    fn, dev = st
    import torch
    a = torch.as_tensor(np.ascontiguousarray(warped_ct, dtype="float32"),
                        device=dev).unsqueeze(0).unsqueeze(0)
    b = torch.as_tensor(np.ascontiguousarray(mr, dtype="float32"),
                        device=dev).unsqueeze(0).unsqueeze(0)
    with torch.no_grad():
        return float(fn(a, b).item())
 def _clip_image(img, lo, hi):
    """Clip image intensities to [lo, hi] (geometry preserved)."""
    a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
@ -153,10 +206,12 @@ def evaluate(ct, t1c, tv, T):
    warped_ct = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
    warped_tv = sitk.GetArrayFromImage(sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0.0))
    gap = m6._skull_ring_gap(warped_ct, ha)
-    # achieved-metric alignment score: Mattes MI between the CT inverse-transformed
+    # achieved-metric alignment scores: Mattes MI between the CT
-    # into MR coordinates (warped_ct) and the original MR (ha), over the MR head.
+    # inverse-transformed into MR coordinates (warped_ct) and the original MR
    # (ha), over the MR head; IMPACT semantic loss over the full volumes.
    ha_f = ha.astype("float32", copy=False)
    mi = _mi_masked(warped_ct, ha_f, head_mask_from_image(ha_f))
    il = _impact_loss_value(warped_ct, ha_f)
    posw = warped_ct[warped_ct > 0]
    tissue_w = (warped_ct > max(2.0, float(np.percentile(posw, 2)))) if posw.size \
        else np.zeros_like(warped_ct, bool)
@ -171,6 +226,7 @@ def evaluate(ct, t1c, tv, T):
    return {
        "skull_ring_gap_mm": gap,
        "mi_ct_mr": round(mi, 4) if mi is not None else None,
        "impact_loss": round(il, 4) if il is not None else None,
        "tv_nvox": nvox,
        "tv_in_ct_frac": round(tv_in_ct_frac, 3),
        "alignment": "concentric" if ok_gap else "not_concentric",
@ -237,6 +293,32 @@ def make_qa_png(warped_ct, mr, out_png, title, block=24, ct_brain_window=(0.0, 8
 def _call(pipeline, ct, t1c, tv, args):
    if args.deformable:
        if pipeline == "simpleitk":
            return RP.register_simpleitk_deformable(ct, t1c, threads=args.threads)
        if pipeline == "elastix":
            return RP.register_elastix_deformable(ct, t1c, workdir=args.elastix_dir,
                                                  threads=args.threads)
        if pipeline == "antspyx":
            return RP.register_antspyx_deformable(ct, t1c)
        if pipeline in RP.ANT_VARIANTS:
            return RP.register_ants_variant_deformable(ct, t1c, variant=pipeline)
        if pipeline == "ants_rigid_synonly":
            return RP.register_ants_rigid_synonly_deformable(ct, t1c)
        if pipeline == "synthmorph":
            return RP.register_synthmorph_deformable(ct, t1c,
                                                     gpu=args.synthmorph_gpu,
                                                     threads=args.threads,
                                                     extent=args.synthmorph_extent)
        if pipeline == "impact":
            return RP.register_impact_deformable(ct, t1c, gpu=args.impact_gpu,
                                                 iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS,
                                                 samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES)
        if pipeline == "niftyreg":
            return RP.register_niftyreg_deformable(ct, t1c, threads=args.threads,
                                                   backend=args.niftyreg_backend,
                                                   gpu=args.niftyreg_gpu)
        raise KeyError(pipeline)
    if pipeline == "simpleitk":
        return RP.register_simpleitk(ct, t1c, threads=args.threads, mode=args.simpleitk_mode)
    if pipeline == "elastix":
@ -248,8 +330,9 @@ def _call(pipeline, ct, t1c, tv, args):
                                      gpu=args.synthmorph_gpu, threads=args.threads,
                                      extent=args.synthmorph_extent)
    if pipeline == "impact":
-        return RP.register_impact(ct, t1c, preset=args.impact_preset, gpu=args.impact_gpu,
+        return RP.register_impact(ct, t1c, gpu=args.impact_gpu,
-                                  iters=args.impact_iters, samples=args.impact_samples)
+                                  iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS,
                                  samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES)
    if pipeline == "niftyreg":
        return RP.register_niftyreg(ct, t1c, threads=args.threads,
                                    backend=args.niftyreg_backend,
@ -268,6 +351,7 @@ def _run_case(pipeline, case, args):
    out = {"key": key, "subject": case["subject"], "date": case["date"]}
    if not ct_complete(ct):
        out.update({"skull_ring_gap_mm": None, "tv_nvox": 0, "tv_in_ct_frac": None,
                    "impact_loss": None,
                    "alignment": "ct_incomplete", "accepted": False,
                    "note": "CT head incomplete (data issue, not a pipeline failure)"})
        return out
@ -313,7 +397,7 @@ def run_pipeline(pipeline, cases, args):
                r = results[-1]
                reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT"
                print(f"   {c['key']}: mi={r.get('mi_ct_mr')} "
-                      f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True)
+                      f"impact_loss={r.get('impact_loss')} {reg}", flush=True)
            except Exception as e:  # noqa
                errors.append({"key": c["key"], "error": repr(e)})
                print(f"   {c['key']}: ERROR {e!r}", flush=True)
@ -325,7 +409,7 @@ def run_pipeline(pipeline, cases, args):
    gaps_mid = [r["skull_ring_gap_mm"]["mid_mm"] for r in registered if r.get("skull_ring_gap_mm")]
    gaps_std = [r["skull_ring_gap_mm"]["std_mm"] for r in registered if r.get("skull_ring_gap_mm")]
    regt = [r["reg_wall_s"] for r in registered]
-    tvfr = [r["tv_in_ct_frac"] for r in registered if r.get("tv_in_ct_frac") is not None]
+    ilv = [r["impact_loss"] for r in registered if r.get("impact_loss") is not None]
    mifr = [r["mi_ct_mr"] for r in registered if r.get("mi_ct_mr") is not None]
    first = (registered or results)[0] if (registered or results) else {}
    row = {
@ -341,7 +425,8 @@ def run_pipeline(pipeline, cases, args):
                                  "median": round(float(np.median(gaps_mid)), 1)} if gaps_mid else None,
        "skull_ring_gap_std_mm": {"mean": round(float(np.mean(gaps_std)), 1),
                                  "median": round(float(np.median(gaps_std)), 1)} if gaps_std else None,
-        "tv_in_ct_frac": {"mean": round(float(np.mean(tvfr)), 3)} if tvfr else None,
+        "impact_loss": {"mean": round(float(np.mean(ilv)), 4),
                        "median": round(float(np.median(ilv)), 4)} if ilv else None,
        "mi_ct_mr": {"mean": round(float(np.mean(mifr)), 4),
                     "median": round(float(np.median(mifr)), 4)} if mifr else None,
        "cases": results, "errors": errors,
@ -367,28 +452,28 @@ def _fmt(v, spec="{}", default="-"):
 def print_table(rows):
-    order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)]
+    order = [p for p in RP.ALL_PIPELINE_NAMES if any(r["pipeline"] == p for r in rows)]
-    print("\n" + "=" * 92)
+    print("\n" + "=" * 96)
    hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} "
-           f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
+           f"{'s/reg':>7s} {'ImpactLoss':>10s} {'MI':>7s}")
    print(hdr)
-    print("-" * 92)
+    print("-" * 96)
    for p in order:
        r = next(x for x in rows if x["pipeline"] == p)
        if not r.get("available"):
            print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} "
-                  f"{'-':>7s} {'-':>7s} {'-':>7s}   [{r['reason']}]")
+                  f"{'-':>7s} {'-':>10s} {'-':>7s}   [{r['reason']}]")
            continue
        rt = (r.get("reg_wall_s") or {}).get("median")
-        tv = (r.get("tv_in_ct_frac") or {}).get("mean")
+        il = (r.get("impact_loss") or {}).get("median")
        mi = (r.get("mi_ct_mr") or {}).get("median")
        print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} "
              f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} "
-              f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
+              f"{_fmt(il, '{:.3f}'):>10s} {_fmt(mi, '{:.3f}'):>7s}")
-    print("=" * 92)
+    print("=" * 96)
-def write_plot(rows, out_png):
+def write_plot(rows, out_png, label="rigid"):
    try:
        import matplotlib
        matplotlib.use("Agg")
@ -402,13 +487,19 @@ def write_plot(rows, out_png):
        return
    labels = list(ok_rows)
    regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") or 0 for p in labels]
    il = [(ok_rows[p].get("impact_loss") or {}).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))
+    fig, axes = plt.subplots(1, 3, figsize=(14, 4))
-    colors = ["#4878CF", "#6ACC65"]
+    colors = ["#6ACC65", "#F5A623", "#4878CF"]
    for ax, (title, vals, fmt, ylab) in zip(axes, [
        ("Registration time (s, median)", regt, "{:.0f}", "s"),
        ("ImpactLoss: warped-CT ~ MR (semantic, median; lower = better)", il, "{:.3f}", "loss"),
        ("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]):
-        bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)])
+        # wide time ranges (deformable: simpleitk ~100x the others) -> log axis
        nz = [v for v in vals if v > 0]
        if len(nz) >= 2 and max(nz) / min(nz) > 30:
            ax.set_yscale("log")
        bars = ax.bar(range(len(vals)), vals, color=colors[0])
        for b, v in zip(bars, vals):
            ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()),
                        ha="center", va="bottom", fontsize=8)
@ -416,7 +507,7 @@ def write_plot(rows, out_png):
        ax.set_ylabel(ylab)
        ax.set_xticks(range(len(vals)))
        ax.set_xticklabels(labels, fontsize=9)
-    plt.suptitle("CT->T1c registration benchmark (M6, cross-modality rigid)", fontsize=11)
+    plt.suptitle(f"CT->T1c registration benchmark (M6, cross-modality {label})", fontsize=11)
    fig.tight_layout()
    fig.savefig(out_png, dpi=120)
    plt.close(fig)
@ -426,7 +517,7 @@ def write_plot(rows, out_png):
 def merge(args):
    out_dir = os.path.join(args.out_dir)
    rows = []
-    for p in RP.ALL_PIPELINES:
+    for p in RP.ALL_PIPELINE_NAMES:
        fp = os.path.join(out_dir, p + ".json")
        if os.path.exists(fp):
            with open(fp) as f:
@ -435,13 +526,30 @@ def merge(args):
        raise SystemExit(f"no pipeline results found in {out_dir}")
    print(f"[merge] {len(rows)} pipeline result files in {out_dir}")
    print_table(rows)
    deformable = (os.path.basename(os.path.normpath(args.out_dir))
              == "bench_registration_deformable")
    if deformable:
        protocol = ("deformable CT->T1c; cross-modality; each pipeline runs its "
                    "toolkit's native deformable workflow (linear init + deformable "
                    "stage, or a single deformable stage where that is native); the "
                    "result is rasterised to the T1c-grid displacement field of the "
                    "unified contract (T1c->CT); alignment scored by achieved Mattes MI "
                    "(nats, warped-CT ~ MR over the MR head, higher = better) and the "
                    "IMPACT semantic loss (vboussot/ImpactLoss, TS/M730 features, "
                    "warped-CT ~ MR over the full volumes, lower = better); reg time "
                    "excludes I/O+scoring; per-case JSON retains the production "
                    "skull-ring gap / acceptance QA")
    else:
        protocol = ("rigid CT->T1c; cross-modality; alignment scored by achieved "
                    "Mattes MI (nats, warped-CT ~ MR over the MR head, higher = "
                    "better) and the IMPACT semantic loss (vboussot/ImpactLoss, "
                    "TS/M730 features, warped-CT ~ MR over the full volumes, lower "
                    "= better); reg time excludes I/O+scoring; per-case JSON retains "
                    "the production skull-ring gap / acceptance QA; best-of-T vs T^-1 "
                    "orientation gate for the elastix-style toolkits")
    summary = {
        "timestamp": datetime.now().isoformat(timespec="seconds"),
-        "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI "
+        "protocol": protocol,
                      "(nats, warped-CT ~ MR over the MR head, higher = better); reg time "
                      "excludes I/O+scoring; per-case JSON retains the production skull-ring "
                      "gap / acceptance QA; best-of-T vs T^-1 orientation gate for the "
                      "elastix-style toolkits"),
        "n_cases": next((r["n_cases"] for r in rows if r.get("available")), None),
        "case_keys": sorted({c["key"] for r in rows if r.get("available") for c in r.get("cases", [])}),
        "results": rows,
@ -450,9 +558,10 @@ def merge(args):
    with open(sf, "w") as f:
        json.dump(summary, f, indent=1)
    print("summary:", sf)
-    write_plot(rows, os.path.join(out_dir, "bench_registration.png"))
+    write_plot(rows, os.path.join(out_dir, "bench_registration.png"),
               label="deformable" if deformable else "rigid")
    # markdown table for the report
-    md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | tv-in-CT | MI (nats, median) |",
+    md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |",
          "|---|---|---|---|---|---|---|---|---|"]
    for p in [x["pipeline"] for x in rows]:
        r = next(x for x in rows if x["pipeline"] == p)
@ -463,7 +572,7 @@ def merge(args):
            f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | "
            f"{r['n_error']} | "
            f"{(r.get('reg_wall_s') or {}).get('median')} | "
-            f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | "
+            f"{(r.get('impact_loss') or {}).get('median')} | "
            f"{(r.get('mi_ct_mr') or {}).get('median')} |")
    mf = os.path.join(out_dir, "summary.md")
    with open(mf, "w") as f:
@ -473,6 +582,10 @@ def merge(args):
 def parse_args():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--deformable", action="store_true",
                    help="run the deformable variants of the pipelines "
                         "(each toolkit's native deformable workflow; results go "
                         "to results/bench_registration_deformable by default)")
    ap.add_argument("--pipelines", default=" ".join(RP.ALL_PIPELINES),
                    help="space-separated subset of: " + " ".join(RP.ALL_PIPELINES))
    ap.add_argument("--n-cases", type=int, default=10, help="number of M6 test cases (default 10)")
@ -496,24 +609,16 @@ def parse_args():
    ap.add_argument("--synthmorph-extent", type=int, default=256,
                    help="SynthMorph isotropic registration-space extent in mm "
                         "(192 or 256; default 256)")
    ap.add_argument("--impact-preset", default=RP.IMPACT_DEFAULT_PRESET,
                    help="ImpactReg preset from the published database "
                         f"(default: {RP.IMPACT_DEFAULT_PRESET}; e.g. "
                         "Elastix_IMPACT_Jacobian, Generic_Rigid, "
                         "Generic_Rigid_BSpline, FireANTs_SyN, FireANTs_IMPACT, "
                         "ConvexAdam_Coarse/Composite; list live with: "
                         "impact-reg-konfai register --help)")
    ap.add_argument("--impact-gpu", type=int, default=None,
                    help="GPU device index for IMPACT (default: auto-spread over "
                         "the host GPUs by pid; -1 forces CPU)")
    ap.add_argument("--impact-iters", type=int, default=None,
-                    help="preset max_iterations override, one value for all "
+                    help="IMPACT MaximumNumberOfIterations, one value for all "
-                         "resolutions (only effective for presets without a "
+                         "levels (default: the reference map's per-level budget "
-                         "per-resolution model matrix; default: the preset's "
+                         f"{RP.IMPACT_DEFAULT_ITERS})")
                         "own budget, 400/300/200/200 for Elastix_IMPACT_Static)")
    ap.add_argument("--impact-samples", type=int, default=None,
-                    help="preset spatial_samples override (NumberOfSpatialSamples; "
+                    help="IMPACT NumberOfSpatialSamples per iteration (default: "
-                         "default: the preset's own value, 2000)")
+                         f"{RP.IMPACT_DEFAULT_SAMPLES})")
    ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"],
                    default="auto",
                    help="niftyreg reg_aladin build: docker = official v2.1.1 in the "
@ -542,6 +647,9 @@ def parse_args():
 def main():
    args = parse_args()
    if args.deformable:
        if args.out_dir == os.path.join(ROOT, "results", "bench_registration"):
            args.out_dir = os.path.join(ROOT, "results", "bench_registration_deformable")
    if args.merge:
        merge(args)
        return
@ -550,7 +658,7 @@ def main():
        args.qa_out = os.path.abspath(args.qa_out or os.path.join(args.out_dir, "qa"))
        os.makedirs(args.qa_out, exist_ok=True)
    pipelines = [p.lower() for p in args.pipelines.split()]
-    unknown = set(pipelines) - set(RP.ALL_PIPELINES)
+    unknown = set(pipelines) - set(RP.ALL_PIPELINE_NAMES)
    if unknown:
        raise SystemExit(f"unknown pipeline(s): {unknown}")
--- a/scripts/investigate_synthmorph.py
+++ b/scripts/investigate_synthmorph.py
@ -0,0 +1,601 @@
 """Investigate the low performance of the ``synthmorph`` pipeline on the M6
 CT->T1c rigid benchmark (see REGISTRATION_BENCHMARK_RIGID.md).
 The benchmark median MI is 0.277 for synthmorph vs ~0.65 for the intensity
 tools, with the ~30 deg repositioned cases collapsing to 0.06-0.11.  This
 script decomposes *why*, scoring everything with the exact benchmark scoring
 (``scripts/bench_registration.py::evaluate``):
 A. Transform decomposition -- for each case, run synthmorph ``rigid`` (the
   benchmark config; deterministic, so it reproduces the bench JSON) and a
   ``simpleitk`` benchmark reference, and report per case:
   - achieved MI / ImpactLoss / skull-ring gap for both,
   - rotation angle + axis of each transform (world RAS frame),
   - the *relative* rotation error  angle(R_ref^T R_sm),
   - the head-centre placement error  |T_sm(c) - T_ref(c)|  (c = T1c head
     centroid),
   - the identity-transform score (the "do nothing" baseline),
   - a numeric rigidity check of the synthmorph matrix (R^T R ~ I, det = 1).
 B. Affine model -- run synthmorph ``affine`` (12-DOF matrix model) on the
   same cases and score it: does the bigger matrix family recover the tilt
   that rigid does not?  The linear part is decomposed into rotation /
   scale / shear (polar).
 C. Joint (affine + deformable) upper bound -- for the worst cases, run
   ``-m joint`` directly in Docker (outside the rigid contract: the output is
   a displacement field, applied with the tool's own ``-o`` resampling) and
   score the warped CT image.  Answers: is the low score a matrix-model
   limitation or a model-family (pretrained-on-MRI) limitation?
 D. Network-space mechanics -- static (no registration) check of what the
   network actually sees, emulating ``network_space`` + ``normalize`` from
   the image code (isotropic 1-mm LIA grid, extent 256, centred on each
   volume's own centre, fill 0, min-max into [0,1]):
   - head truncation: fraction of head-mask voxels (the T1c head mapped into
     CT space via the reference transform) outside the 256 mm box, for the CT
     and T1c network spaces,
   - intensity band: the min-max-normalized values over the head region for
     the T1c, the [0,80]-clipped CT (the benchmark default) and the raw CT.
 Usage (repo root, conda env longitudinal-313):
    python scripts/investigate_synthmorph.py --workers 3
    python scripts/investigate_synthmorph.py --parts D          # static only
    python scripts/investigate_synthmorph.py --keys m6_2BBN5XA4_20250411
 """
 import argparse
 import importlib.util
 import json
 import math
 import os
 import shutil
 import subprocess
 import sys
 import tempfile
 import time
 from concurrent.futures import ProcessPoolExecutor, as_completed
 ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 sys.path.insert(0, ROOT)
 import numpy as np
 import SimpleITK as sitk
 from scipy.ndimage import center_of_mass
 from src.common import d, read_nii_img, head_mask_from_image
 from src import registration_pipelines as RP
 def _load_module(name, rel):
    p = os.path.join(ROOT, rel)
    spec = importlib.util.spec_from_file_location(name, p)
    m = importlib.util.module_from_spec(spec)
    spec.loader.exec_module(m)
    return m
 BR = _load_module("bench_registration", "scripts/bench_registration.py")
 M6 = RP._load_m6()
 # the worst 3 bench cases (rigid MI): ~30 deg tilt x2 + a mid-tier case
 JOINT_KEYS = ("m6_2BBN5XA4_20250411", "m6_24V62QQM_20240820",
              "m6_24V62QQM_20230223")
 # --------------------------------------------------------------------------- #
 # transform decomposition
 # ---------------------------------------------------------------------------
 def _rot_angle(R):
    c = (np.trace(R) - 1.0) / 2.0
    return math.degrees(math.acos(float(np.clip(c, -1.0, 1.0))))
 def _rot_axis(R, theta_deg):
    if theta_deg <= 0.5:
        return np.zeros(3)
    axis = np.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]])
    n = np.linalg.norm(axis)
    return axis / n if n > 1e-9 else np.zeros(3)
 def _to_affine(T):
    """Recover (R, t) of any SimpleITK transform by point mapping (world RAS)."""
    o = np.array(T.TransformPoint((0.0, 0.0, 0.0)), dtype=float)
    R = np.empty((3, 3), dtype=float)
    for i in range(3):
        p = [0.0, 0.0, 0.0]
        p[i] = 1.0
        R[:, i] = np.array(T.TransformPoint(tuple(p)), dtype=float) - o
    return R, o
 def _decompose_affine(T):
    """(R, t, theta_deg, axis, rigidity) of a SimpleITK transform in world RAS."""
    R, t = _to_affine(T)
    theta = _rot_angle(R)
    rig = {"det": round(float(np.linalg.det(R)), 6),
           "max_devi_from_I": round(float(np.linalg.norm(R.T @ R - np.eye(3))), 6)}
    return R, t, theta, _rot_axis(R, theta), rig
 def _angle_between_rotations(A, B):
    """angle of A^T B in degrees (rotation B needs on top of A)."""
    return _rot_angle(A.T @ B)
 def _polar_rotation(M):
    """rotation part of a general 3x3 linear map (SVD polar)."""
    U, _s, Vt = np.linalg.svd(M)
    R_p = U @ Vt
    if np.linalg.det(R_p) < 0:
        U[:, -1] *= -1
        R_p = U @ Vt
    return R_p
 def head_centroid_world(t1c):
    ha = sitk.GetArrayFromImage(t1c).astype("float32")
    idx = center_of_mass(head_mask_from_image(ha))
    sp = np.array(t1c.GetSpacing())
    o = np.array(t1c.GetOrigin())
    D = np.array(t1c.GetDirection()).reshape(3, 3)
    return o + D @ (np.array(idx[::-1]) * sp)
 def score_identity(ct, t1c):
    """The "do nothing" baseline: CT resampled onto the T1c grid with the
    identity transform (i.e. no registration at all)."""
    ha = sitk.GetArrayFromImage(t1c).astype("float32")
    wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, sitk.Transform(),
                                              sitk.sitkLinear))
    gap = M6._skull_ring_gap(wc, ha)
    mi = BR._mi_masked(wc, ha, head_mask_from_image(ha))
    return {"mi": mi, "gap": gap}
 # --------------------------------------------------------------------------- #
 # part A/B/D worker
 # ---------------------------------------------------------------------------
 def _case_task(case, args):
    out = {"key": case["key"], "subject": case["subject"], "date": case["date"]}
    ct = read_nii_img(case["ct"])
    t1c = read_nii_img(case["img"])
    tv = read_nii_img(case["label"])
    if not BR.ct_complete(ct):
        out["note"] = "ct_incomplete (skipped, same as the benchmark)"
        return out
    gpu = (args.gpu_base + abs(hash(out["key"])) % max(1, args.n_gpus)) % max(1, args.n_gpus)
    parts = set(args.parts)
    # geometry facts
    for nm, img in (("ct", ct), ("t1c", t1c)):
        sz = img.GetSize()
        sp = img.GetSpacing()
        c = np.array(img.TransformIndexToPhysicalPoint(
            ((sz[0] - 1) // 2, (sz[1] - 1) // 2, (sz[2] - 1) // 2)))
        out[nm] = {"shape_xyz": list(sz),
                   "spacing_xyz": [round(float(s), 3) for s in sp],
                   "fov_mm_xyz": [round(float(s0 * s1), 1) for s0, s1 in zip(sz, sp)],
                   "center_world": [round(float(x), 1) for x in c]}
    c_t1c = head_centroid_world(t1c)
    out["t1c_head_centroid_world"] = [round(float(x), 1) for x in c_t1c]
    R_ref = t_ref = None
    if parts & {"A", "B", "D"}:
        # reference: simpleitk benchmark mode (deterministic at fixed threads)
        T_ref, meta_ref = RP.register_simpleitk(ct, t1c, threads=16,
                                                mode="benchmark")
        sref = BR.evaluate(ct, t1c, tv, T_ref)
        R_ref, t_ref, th_ref, ax_ref, _ = _decompose_affine(T_ref)
        out["reference"] = {
            "tool": meta_ref["toolkit"],
            "R": R_ref.tolist(), "t": t_ref.tolist(),
            "mi": sref["mi_ct_mr"], "impact_loss": sref["impact_loss"],
            "gap": sref["skull_ring_gap_mm"],
            "rotation_deg": round(th_ref, 2),
            "axis": [round(float(a), 3) for a in ax_ref],
            "translation": [round(float(x), 1) for x in t_ref],
            "head_centre_ct": [round(float(x), 1) for x in (R_ref @ c_t1c + t_ref)],
        }
        out["identity"] = score_identity(ct, t1c)
    if "A" in parts:
        t0 = time.time()
        T_sm, meta_sm = RP.register_synthmorph(ct, t1c, model="rigid",
                                               gpu=int(gpu), extent=256)
        ssm = BR.evaluate(ct, t1c, tv, T_sm)
        R_sm, t_sm, th_sm, ax_sm, rig = _decompose_affine(T_sm)
        p_ref = R_ref @ c_t1c + t_ref
        p_sm = R_sm @ c_t1c + t_sm
        out["rigid"] = {
            "seconds": round(time.time() - t0, 1),
            "auto_inverted": meta_sm["auto_inverted"],
            "clip_ct": meta_sm["clip_ct"],
            "mi": ssm["mi_ct_mr"], "impact_loss": ssm["impact_loss"],
            "gap": ssm["skull_ring_gap_mm"],
            "rotation_deg": round(th_sm, 2),
            "axis": [round(float(a), 3) for a in ax_sm],
            "translation": [round(float(x), 1) for x in t_sm],
            "matrix_rigidity": rig,
            "rot_err_vs_ref_deg": round(_angle_between_rotations(R_ref, R_sm), 2),
            "head_centre_err_mm": round(float(np.linalg.norm(p_sm - p_ref)), 1),
        }
    if "B" in parts:
        t0 = time.time()
        T_af, meta_af = RP.register_synthmorph(ct, t1c, model="affine",
                                               gpu=int(gpu), extent=256)
        saf = BR.evaluate(ct, t1c, tv, T_af)
        M_af, t_af = _to_affine(T_af)
        svals = np.linalg.svd(M_af, compute_uv=False)
        lin = {"scale_mean": round(float(np.mean(svals)), 4),
               "scale_iso_dev": round(float(np.std(svals)) / max(float(np.mean(svals)), 1e-12), 4),
               "rotation_deg": round(_rot_angle(_polar_rotation(M_af)), 2),
               "det": round(float(np.linalg.det(M_af)), 4)}
        row = {"seconds": round(time.time() - t0, 1),
               "auto_inverted": meta_af["auto_inverted"],
               "mi": saf["mi_ct_mr"], "impact_loss": saf["impact_loss"],
               "gap": saf["skull_ring_gap_mm"], "linear": lin}
        if R_ref is not None:
            row["rot_err_vs_ref_deg"] = round(
                _angle_between_rotations(R_ref, _polar_rotation(M_af)), 2)
            row["head_centre_err_mm"] = round(
                float(np.linalg.norm(M_af @ c_t1c + t_af - (R_ref @ c_t1c + t_ref))), 1)
        out["affine"] = row
    if "D" in parts:
        out["network_space"] = _network_space_analysis(ct, t1c, out,
                                                       R_ref=R_ref, t_ref=t_ref)
    return out
 # --------------------------------------------------------------------------- #
 # part D: what the network sees
 # ---------------------------------------------------------------------------
 def _network_grid(img, extent=256):
    """The matrix-model network space of ``img``: isotropic 1-mm LIA grid,
    ``extent``^3, centred on the volume's own centre (mirrors
    ``network_space`` in the container's registration.py)."""
    sz = img.GetSize()
    c = np.array(img.TransformIndexToPhysicalPoint(
        ((sz[0] - 1) // 2, (sz[1] - 1) // 2, (sz[2] - 1) // 2)))
    # LIA axes in the RAS world: x_L = -x_R, y_I = -y_A, z_A = +z_R
    origin = (c[0] + extent / 2.0 - 0.5, c[1] + extent / 2.0 - 0.5,
              c[2] - extent / 2.0 + 0.5)
    net = sitk.Image((extent, extent, extent), sitk.sitkFloat32)
    net.SetSpacing((1.0, 1.0, 1.0))
    net.SetOrigin(tuple(float(x) for x in origin))
    net.SetDirection((-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0))
    return net, c
 def _to_network(img, net):
    a = sitk.GetArrayFromImage(
        sitk.Resample(img, net, sitk.Transform(), sitk.sitkLinear, 0.0))
    lo, hi = float(a.min()), float(a.max())
    if hi - lo < 1e-9:
        return np.zeros_like(a), lo, hi
    return (a - lo) / (hi - lo), lo, hi
 def _network_space_analysis(ct, t1c, out, extent=256, R_ref=None, t_ref=None):
    """Extent truncation + min-max intensity band (what the model sees)."""
    res = {"extent": extent}
    ha = sitk.GetArrayFromImage(t1c).astype("float32")
    head_t1c = sitk.GetImageFromArray(head_mask_from_image(ha).astype("float32"))
    head_t1c.CopyInformation(t1c)
    if R_ref is not None:
        T_ref = sitk.AffineTransform(3)
        T_ref.SetMatrix(list(np.array(R_ref).ravel()))
        T_ref.SetTranslation(list(np.array(t_ref)))
        # T_ref maps T1c -> CT, so rasterising the T1c head onto the CT grid
        # needs the inverse as the destination -> source map
        head_ct = sitk.Resample(head_t1c, ct, T_ref.GetInverse(),
                                sitk.sitkNearestNeighbor, 0.0)
    else:
        head_ct = head_t1c
    # truncation: head voxels outside each volume's OWN extent box
    for nm, img, head in (("ct", ct, head_ct), ("t1c", t1c, head_t1c)):
        _net, c = _network_grid(img, extent)
        idx = np.argwhere(sitk.GetArrayFromImage(head) > 0.5)
        if idx.size == 0:
            res[f"{nm}_head_trunc_frac"] = None
            continue
        sp = np.array(img.GetSpacing())
        o = np.array(img.GetOrigin())
        D = np.array(img.GetDirection()).reshape(3, 3)
        world = o[None, :] + (idx[:, ::-1] * sp) @ D.T  # idx is (z,y,x)
        half = extent / 2.0
        inside = (np.abs(world - c[None, :]) <= half).all(axis=1)
        head_c = world.mean(axis=0)
        res[f"{nm}_head_voxels"] = int(idx.shape[0])
        res[f"{nm}_head_trunc_frac"] = round(float(1.0 - inside.mean()), 4)
        res[f"{nm}_head_off_volume_centre_mm"] = \
            [round(float(x), 1) for x in (head_c - c)]
    # min-max normalized band over the head region (emulating normalize=True)
    bands = {}
    for nm, img, head in (("t1c", t1c, head_t1c),
                          ("ct_clip", RP._clip_intensity(ct, 0.0, 80.0), head_ct),
                          ("ct_raw", ct, head_ct)):
        net, _c = _network_grid(img, extent)
        a, lo, hi = _to_network(img, net)
        hnet = sitk.GetArrayFromImage(
            sitk.Resample(head, net, sitk.Transform(),
                          sitk.sitkNearestNeighbor, 0.0)) > 0.5
        v = a[hnet]
        bands[nm] = {"src_min": round(lo, 1), "src_max": round(hi, 1),
                     "head_nvox_net": int(v.size),
                     "p5": round(float(np.percentile(v, 5)), 4) if v.size else None,
                     "p50": round(float(np.percentile(v, 50)), 4) if v.size else None,
                     "p95": round(float(np.percentile(v, 95)), 4) if v.size else None,
                     "band_width_p95_p5": (round(float(np.percentile(v, 95)
                                                       - np.percentile(v, 5)), 4)
                                            if v.size else None)}
    res["normalized_band"] = bands
    return res
 # --------------------------------------------------------------------------- #
 # part C: joint (affine + deformable) upper bound
 # ---------------------------------------------------------------------------
 def _joint_task(case, args):
    out = {"key": case["key"], "subject": case["subject"], "date": case["date"],
           "model": "joint (affine + deformable)"}
    ct = read_nii_img(case["ct"])
    t1c = read_nii_img(case["img"])
    if not BR.ct_complete(ct):
        out["note"] = "ct_incomplete (skipped)"
        return out
    gpu = (args.gpu_base + abs(hash(out["key"])) % max(1, args.n_gpus)) % max(1, args.n_gpus)
    tmp = tempfile.mkdtemp(prefix="sm_joint_")
    try:
        ct_clip = RP._clip_intensity(ct, 0.0, 80.0)
        sitk.WriteImage(sitk.Cast(ct_clip, sitk.sitkFloat32),
                        os.path.join(tmp, "moving.nii.gz"), True)
        sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
                        os.path.join(tmp, "fixed.nii.gz"), True)
        # CLI order is (moving, fixed); -o saves moving registered to fixed
        cmd = ["docker", "run", "--rm",
               "--user", f"{os.getuid()}:{os.getgid()}",
               "-e", "TF_FORCE_GPU_ALLOW_GROWTH=true",
               "--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}",
               "-v", f"{tmp}:/data", RP.SYNTHMORPH_IMAGE,
               "register", "-m", "joint", "-e", "256", "-g",
               "-o", "/data/warped.nii.gz",
               "/data/moving.nii.gz", "/data/fixed.nii.gz"]
        t0 = time.time()
        proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200)
        dt = time.time() - t0
        if proc.returncode != 0:
            out["error"] = (f"docker failed rc={proc.returncode}:\n"
                            + "\n".join((proc.stdout + proc.stderr).splitlines()[-15:]))
            return out
        warped = sitk.GetArrayFromImage(
            sitk.ReadImage(os.path.join(tmp, "warped.nii.gz"))).astype("float32")
        ha = sitk.GetArrayFromImage(t1c).astype("float32")
        gap = M6._skull_ring_gap(warped, ha)
        mi = BR._mi_masked(warped, ha, head_mask_from_image(ha))
        il = BR._impact_loss_value(warped, ha)
        out.update({"seconds": round(dt, 1), "mi": mi, "impact_loss": il,
                    "gap": gap, "gpu": int(gpu)})
    finally:
        if not args.keep_scratch:
            shutil.rmtree(tmp, ignore_errors=True)
    return out
 # --------------------------------------------------------------------------- #
 # part F: re-centred CT probe (volume-centre anchoring hypothesis)
 # ---------------------------------------------------------------------------
 def _recenter_task(case, args, ref_map):
    """H5 probe: the matrix model's network space is centred on each volume's
    OWN centre.  The T1c is a tight head crop (centre ~ head centre), but the
    head+neck CT's volume centre sits 15-80 mm below the head.  If the model
    under-corrects that offset, cropping the CT to the head box (re-centring
    the volume on the head) should remove the translation miss.  The crop is
    world-preserving, so the returned transform scores against the FULL CT
    exactly like the benchmark."""
    out = {"key": case["key"], "subject": case["subject"], "date": case["date"]}
    ct = read_nii_img(case["ct"])
    t1c = read_nii_img(case["img"])
    tv = read_nii_img(case["label"])
    if not BR.ct_complete(ct):
        out["note"] = "ct_incomplete (skipped)"
        return out
    ref = ref_map.get(case["key"])
    if not ref:
        out["error"] = "no reference transform in partABD.json"
        return out
    R = np.array(ref["R"])
    t = np.array(ref["t"])
    T_ref = sitk.AffineTransform(3)
    T_ref.SetMatrix(list(R.ravel()))
    T_ref.SetTranslation(list(t))
    gpu = (args.gpu_base + abs(hash(out["key"])) % max(1, args.n_gpus)) % max(1, args.n_gpus)
    # head in CT space (T_ref maps T1c -> CT, so rasterise with its inverse)
    ha = sitk.GetArrayFromImage(t1c).astype("float32")
    head_t1c = sitk.GetImageFromArray(head_mask_from_image(ha).astype("float32"))
    head_t1c.CopyInformation(t1c)
    head_ct = sitk.Resample(head_t1c, ct, T_ref.GetInverse(),
                            sitk.sitkNearestNeighbor, 0.0)
    m = sitk.GetArrayFromImage(head_ct) > 0.5
    idx = np.argwhere(m)  # (z, y, x)
    lo, hi = idx.min(0), idx.max(0)
    cen = (lo + hi + 1) / 2.0
    half = (hi - lo) / 2.0 + 15.0 / np.array(ct.GetSpacing())[::-1]  # (z,y,x)
    sz = np.array(ct.GetSize())[::-1]  # GetSize is (x,y,z); work in (z,y,x)
    b0 = np.clip(np.floor(cen - half).astype(int), 0, sz - 1)
    b1 = np.clip(np.ceil(cen + half).astype(int), 0, sz)
    def _volcentre(img):
        s = img.GetSize()
        return np.array(img.TransformIndexToPhysicalPoint(
            ((s[0] - 1) // 2, (s[1] - 1) // 2, (s[2] - 1) // 2)))
    hc = T_ref.TransformPoint(tuple(float(x) for x in head_centroid_world(t1c)))
    out["before"] = {"head_off_volume_centre_mm": [
        round(float(x), 1) for x in (np.array(hc) - _volcentre(ct))]}
    # SimpleITK image slicing is in (x, y, z) order; b0/b1 are (z, y, x)
    ct_crop = ct[int(b0[2]):int(b1[2]), int(b0[1]):int(b1[1]),
                 int(b0[0]):int(b1[0])]
    out["crop"] = {"shape": list(ct_crop.GetSize()),
                   "head_off_volume_centre_mm": [
                       round(float(x), 1)
                       for x in (np.array(hc) - _volcentre(ct_crop))]}
    t0 = time.time()
    T_sm, meta_sm = RP.register_synthmorph(ct_crop, t1c, model="rigid",
                                           gpu=int(gpu), extent=256)
    ssm = BR.evaluate(ct, t1c, tv, T_sm)  # full CT, as in the benchmark
    R_sm, t_sm = _to_affine(T_sm)
    c_t1c = head_centroid_world(t1c)
    out["rigid_recentered"] = {
        "seconds": round(time.time() - t0, 1),
        "auto_inverted": meta_sm["auto_inverted"],
        "mi": ssm["mi_ct_mr"], "impact_loss": ssm["impact_loss"],
        "gap": ssm["skull_ring_gap_mm"],
        "rot_err_vs_ref_deg": round(_angle_between_rotations(R, R_sm), 2),
        "head_centre_err_mm": round(
            float(np.linalg.norm((R_sm @ c_t1c + t_sm) - (R @ c_t1c + t))), 1),
    }
    return out
 # --------------------------------------------------------------------------- #
 # main
 # ---------------------------------------------------------------------------
 def parse_args():
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--parts", default="ABCD",
                    help="which parts to run (letters A/B/C/D/F; F = re-centred "
                         "CT probe, needs the partABD.json reference)")
    ap.add_argument("--n-cases", type=int, default=10)
    ap.add_argument("--keys", default=None, help="comma-separated case keys")
    ap.add_argument("--workers", type=int, default=3)
    ap.add_argument("--gpu-base", type=int, default=0)
    ap.add_argument("--n-gpus", type=int, default=3)
    ap.add_argument("--out-dir", default=os.path.join(ROOT, "results",
                                                      "synthmorph_investigation"))
    ap.add_argument("--keep-scratch", action="store_true")
    return ap.parse_args()
 def main():
    args = parse_args()
    os.makedirs(args.out_dir, exist_ok=True)
    labeled, _ = M6.discover()
    pool = list(labeled)
    if args.keys:
        bykey = {r["key"]: r for r in pool}
        missing = [k for k in args.keys.split(",") if k not in bykey]
        if missing:
            raise SystemExit(f"keys not found: {missing}")
        cases = [bykey[k.strip()] for k in args.keys.split(",")]
    else:
        cases = pool[:args.n_cases]
    print(f"[investigate] {len(cases)} cases, parts={args.parts}")
    rows, joint_rows, recenter_rows = [], [], []
    do_cases = bool(set(args.parts) & {"A", "B", "D"})
    do_joint = "C" in args.parts
    do_recenter = "F" in args.parts
    ref_map = {}
    if do_recenter:
        fp0 = os.path.join(args.out_dir, "partABD.json")
        if os.path.exists(fp0):
            with open(fp0) as f:
                ref_map = {r["key"]: r["reference"]
                           for r in json.load(f) if "reference" in r}
    with ProcessPoolExecutor(max_workers=args.workers) as ex:
        futs = {}
        if do_cases:
            for c in cases:
                futs[ex.submit(_case_task, c, args)] = ("case", c)
        if do_joint:
            for c in cases:
                if c["key"] in JOINT_KEYS:
                    futs[ex.submit(_joint_task, c, args)] = ("joint", c)
        if do_recenter:
            for c in cases:
                futs[ex.submit(_recenter_task, c, args, ref_map)] = ("recenter", c)
        for fu in as_completed(futs):
            kind, c = futs[fu]
            try:
                r = fu.result()
                if kind == "joint":
                    joint_rows.append(r)
                    print(f"   [joint] {c['key']}: mi={r.get('mi')} "
                          f"il={r.get('impact_loss')} {r.get('seconds')}s "
                          f"err={(r.get('error') or '')[:100]}", flush=True)
                elif kind == "recenter":
                    recenter_rows.append(r)
                    rr = r.get("rigid_recentered", {})
                    print(f"   [recenter] {c['key']}: mi={rr.get('mi')} "
                          f"ctr_err={rr.get('head_centre_err_mm')} "
                          f"rot_err={rr.get('rot_err_vs_ref_deg')} "
                          f"off_before={r.get('before', {}).get('head_off_volume_centre_mm')} "
                          f"off_after={r.get('crop', {}).get('head_off_volume_centre_mm')} "
                          f"note={(r.get('note') or r.get('error') or '')[:60]}", flush=True)
                else:
                    rows.append(r)
                    print(f"   [case] {c['key']}: ref.mi={r.get('reference', {}).get('mi')} "
                          f"rigid.mi={r.get('rigid', {}).get('mi')} "
                          f"affine.mi={r.get('affine', {}).get('mi')} "
                          f"note={r.get('note', r.get('error', ''))[:60]}", flush=True)
            except Exception as e:  # noqa
                print(f"   {c['key']}: ERROR {e!r}", flush=True)
                (joint_rows if kind == "joint" else
                 recenter_rows if kind == "recenter" else rows).append(
                    {"key": c["key"], "error": repr(e)})
    rows.sort(key=lambda x: x["key"])
    joint_rows.sort(key=lambda x: x["key"])
    recenter_rows.sort(key=lambda x: x["key"])
    if do_cases:
        fp = os.path.join(args.out_dir, "partABD.json")
        merged = {r["key"]: r for r in rows}
        if os.path.exists(fp):
            with open(fp) as f:
                for r in json.load(f):
                    # merge missing sections (a partial run must not wipe
                    # data from a fuller earlier run)
                    if r["key"] in merged:
                        for k, v in r.items():
                            merged[r["key"]].setdefault(k, v)
                    else:
                        merged[r["key"]] = r
        with open(fp, "w") as f:
            json.dump(sorted(merged.values(), key=lambda x: x["key"]),
                      f, indent=1)
        print("wrote", fp)
    if do_joint:
        fp = os.path.join(args.out_dir, "partC_joint.json")
        with open(fp, "w") as f:
            json.dump(joint_rows, f, indent=1)
        print("wrote", fp)
    if do_recenter:
        fp = os.path.join(args.out_dir, "partF_recenter.json")
        with open(fp, "w") as f:
            json.dump(recenter_rows, f, indent=1)
        print("wrote", fp)
 if __name__ == "__main__":
    main()
--- a/scripts/investigate_synthmorph_mm.py
+++ b/scripts/investigate_synthmorph_mm.py
@ -0,0 +1,101 @@
 """MR->MR control for the synthmorph investigation.
 The benchmark pairs are CT->T1c (out-of-domain for the MRI-trained
 SynthMorph model).  To separate the cross-modality effect from the other
 suspects (extent truncation of the CT FOV, min-max normalization), this
 registers two T1c scans of the SAME subject (m6_24V62QQM, 2023-02-23 and
 2023-10-17 -- in-domain, MRI->MRI, same tightly-cropped geometry as the
 benchmark's fixed image) and scores the result with the benchmark scoring
 (MI over the head, ImpactLoss, skull-ring gap) plus the transform
 decomposition against a SimpleITK reference.
 Usage:
    python scripts/investigate_synthmorph_mm.py [--out results/synthmorph_investigation]
 """
 import argparse
 import json
 import os
 import sys
 import time
 ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 sys.path.insert(0, ROOT)
 import numpy as np
 import SimpleITK as sitk
 from src.common import read_nii_img, head_mask_from_image
 from src import registration_pipelines as RP
 from scripts.investigate_synthmorph import (
    _to_affine, _angle_between_rotations, head_centroid_world, BR, M6)
 KEY_A = "m6_24V62QQM_20230223"
 KEY_B = "m6_24V62QQM_20231017"
 def score_pair(img_moving, img_fixed, T):
    """Benchmark scoring of the warped moving (in the fixed grid) ~ fixed."""
    ha = sitk.GetArrayFromImage(img_fixed).astype("float32")
    wc = sitk.GetArrayFromImage(sitk.Resample(img_moving, img_fixed, T,
                                              sitk.sitkLinear)).astype("float32")
    gap = M6._skull_ring_gap(wc, ha)
    mi = BR._mi_masked(wc, ha, head_mask_from_image(ha))
    il = BR._impact_loss_value(wc, ha)
    return {"mi": mi, "impact_loss": il, "gap": gap}
 def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--out", default=os.path.join(ROOT, "results",
                                                  "synthmorph_investigation"))
    ap.add_argument("--gpu", type=int, default=0)
    args = ap.parse_args()
    os.makedirs(args.out, exist_ok=True)
    labeled, _ = M6.discover()
    bykey = {r["key"]: r for r in labeled}
    a, b = bykey[KEY_A], bykey[KEY_B]
    img_a = read_nii_img(a["img"])
    img_b = read_nii_img(b["img"])
    out = {"pair": [KEY_A, KEY_B], "note": ("same-subject T1c->T1c, 8 months "
                                            "apart; in-domain control for the "
                                            "CT->T1c benchmark pairs")}
    import math
    c_fix = head_centroid_world(img_b)
    T_ref, meta_ref = RP.register_simpleitk(img_a, img_b, threads=16,
                                            mode="benchmark")
    sref = score_pair(img_a, img_b, T_ref)
    R_ref, t_ref = _to_affine(T_ref)
    th_ref = math.degrees(float(np.clip((np.trace(R_ref) - 1) / 2, -1, 1)))
    out["reference"] = {
        "tool": meta_ref["toolkit"], "R": R_ref.tolist(), "t": t_ref.tolist(),
        "mi": sref["mi"], "impact_loss": sref["impact_loss"], "gap": sref["gap"],
        "rotation_deg": round(th_ref, 2),
        "head_centre_fix": [round(float(x), 1)
                            for x in (R_ref @ c_fix + t_ref)],
    }
    t0 = time.time()
    T_sm, meta_sm = RP.register_synthmorph(img_a, img_b, model="rigid",
                                           gpu=args.gpu, extent=256)
    ssm = score_pair(img_a, img_b, T_sm)
    R_sm, t_sm = _to_affine(T_sm)
    out["synthmorph_rigid"] = {
        "seconds": round(time.time() - t0, 1),
        "auto_inverted": meta_sm["auto_inverted"],
        "mi": ssm["mi"], "impact_loss": ssm["impact_loss"], "gap": ssm["gap"],
        "rot_err_vs_ref_deg": round(_angle_between_rotations(R_ref, R_sm), 2),
        "head_centre_err_mm": round(
            float(np.linalg.norm((R_sm @ c_fix + t_sm) - (R_ref @ c_fix + t_ref))), 1),
    }
    # no-clip sanity: MR inputs are never clipped by the adapter
    fp = os.path.join(args.out, "partE_mr_mr.json")
    with open(fp, "w") as f:
        json.dump(out, f, indent=1)
    print(json.dumps(out, indent=1))
    print("wrote", fp)
 if __name__ == "__main__":
    main()
--- a/scripts/setup_impact_konfai.sh
+++ b/scripts/setup_impact_konfai.sh
@ -1,110 +0,0 @@
 #!/usr/bin/env bash
 # One-time host setup for the IMPACT registration pipeline (impact-reg-konfai).
 #
 # The stock elastix-impact 1.0.0 binary (vboussot/ImpactElastix) needs
 # glibc >= 2.35 and links LibTorch 2.8.0. This host (Ubuntu 20.04, glibc 2.31,
 # env torch 2.14) satisfies neither, so this script builds a self-contained
 # shim at ~/.cache/konfai/elastix-impact-shim:
 #   * glibc 2.35 loader + libs, extracted from the niftyreg:2.1.1-ubuntu image
 #     (ubuntu 22.04 base) -- docker is needed ONLY for this one-time copy;
 #   * LibTorch 2.8.0 cu128, self-contained (its own CUDA runtime), from
 #     download.pytorch.org;
 #   * a wrapper bin/elastix that runs the stock binary through the glibc 2.35
 #     loader with the LibTorch 2.8.0 library path (ignoring the inherited
 #     LD_LIBRARY_PATH, so the env's torch 2.14 cannot leak in).
 # The pipeline selects the shim via KONFAI_ELASTIX_DIR
 # (src/registration_pipelines.py::_impact_elastix_root).
 #
 # Run from the repo root, inside the conda env:
 #   source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
 #   pip install impact-reg-konfai        # the CLI itself
 #   bash scripts/setup_impact_konfai.sh  # this script
 #
 # First `register` afterwards auto-downloads from HF/GitHub: the preset app
 # (VBoussot/ImpactReg), the feature models (VBoussot/impact-torchscript-models)
 # and, on hosts that CAN run the stock binary directly, the binary itself.
 set -euo pipefail
 KONFAI_CACHE="${KONFAI_CACHE:-$HOME/.cache/konfai}"
 SHIM="$KONFAI_CACHE/elastix-impact-shim"
 STOCK="$KONFAI_CACHE/elastix-impact"
 GLIBC_SRC_DIR="${GLIBC_SRC_DIR:-glibc-2.35}"
 GLIBC_IMG="${GLIBC_IMG:-niftyreg:2.1.1-ubuntu}"
 LT_VER="2.8.0"
 LT_DIR="libtorch-${LT_VER}-cu128"
 LT_URL="https://download.pytorch.org/libtorch/cu128/libtorch-shared-with-deps-${LT_VER}%2Bcu128.zip"
 echo "[1/5] stock elastix-impact binary (package auto-install)"
 if [ ! -x "$STOCK/bin/elastix" ]; then
    python - <<'PY'
 from pathlib import Path
 from impact_reg_konfai.models.elastix_install import (
    ELASTIX_CACHE, get_elastix_bin, install_elastix_impact)
 if not get_elastix_bin(ELASTIX_CACHE).exists():
    # flavor auto-detected (cu128 on this host); downloads + extracts the zip
    install_elastix_impact(ELASTIX_CACHE, force_cuda=False, force_cpu=False)
 PY
 fi
 [ -x "$STOCK/bin/elastix" ] || { echo "stock binary missing: $STOCK/bin/elastix"; exit 1; }
 echo "[2/5] glibc 2.35 from $GLIBC_IMG (docker, one-time)"
 if [ ! -x "$KONFAI_CACHE/$GLIBC_SRC_DIR/ld-linux-x86-64.so.2" ]; then
    docker run --rm --entrypoint sh -v "$KONFAI_CACHE":/dest "$GLIBC_IMG" -c \
        "rm -rf /dest/$GLIBC_SRC_DIR && mkdir -p /dest/$GLIBC_SRC_DIR && \
         cp -a /lib/x86_64-linux-gnu/. /dest/$GLIBC_SRC_DIR/ && \
         cp -a /usr/lib/x86_64-linux-gnu/libstdc++.so.6* /dest/$GLIBC_SRC_DIR/ 2>/dev/null || true"
 fi
 echo "[3/5] LibTorch ${LT_VER} cu128 (self-contained)"
 if [ ! -f "$SHIM/lib/$LT_DIR/libtorch/lib/libc10.so" ]; then
    mkdir -p "$SHIM/lib"
    curl -L --fail -o "$KONFAI_CACHE/$LT_DIR.zip" "$LT_URL"
    unzip -q "$KONFAI_CACHE/$LT_DIR.zip" -d "$SHIM/lib/$LT_DIR"
    rm "$KONFAI_CACHE/$LT_DIR.zip"
 fi
 # the with-deps zip carries the CUDA libs under hashed names; expose the public
 # soname elastix links directly
 LT_LIB="$SHIM/lib/$LT_DIR/libtorch/lib"
 CUDART="$(ls "$LT_LIB"/libcudart-*.so.12 | head -1)"
 ln -sf "$(basename "$CUDART")" "$LT_LIB/libcudart.so.12"
 echo "[4/5] shim layout + wrapper"
 mkdir -p "$SHIM/bin"
 cp "$STOCK/bin/elastix" "$SHIM/bin/elastix-real"
 [ -d "$STOCK/lib" ] && cp -a "$STOCK/lib/." "$SHIM/lib/"
 [ -d "$KONFAI_CACHE/$GLIBC_SRC_DIR" ] && cp -a "$KONFAI_CACHE/$GLIBC_SRC_DIR" "$SHIM/lib/$GLIBC_SRC_DIR"
 cat > "$SHIM/bin/elastix" <<'WRAP'
 #!/bin/sh
 # elastix-impact 1.0.0 shim: the host glibc (2.31) is too old for the binary
 # (built against glibc 2.35) and the binary expects LibTorch 2.8.0, which the
 # conda env's torch (2.14) does not provide. Run it through the bundled
 # glibc 2.35 loader with the bundled self-contained LibTorch 2.8.0 (cu128).
 ROOT="$(cd "$(dirname "$0")/.." && pwd)"
 LT="$ROOT/lib/libtorch-2.8.0-cu128/libtorch/lib"
 GLIBC="$ROOT/lib/glibc-2.35"
 exec env -u LD_LIBRARY_PATH -u LD_PRELOAD \
  "$GLIBC/ld-linux-x86-64.so.2" \
  --library-path "$LT:$ROOT/lib:$GLIBC" \
  "$ROOT/bin/elastix-real" "$@"
 WRAP
 chmod +x "$SHIM/bin/elastix"
 echo "[5/5] verify"
 KONFAI_ELASTIX_DIR="$SHIM" python - <<'PY'
 import os
 from pathlib import Path
 from impact_reg_konfai.models.elastix_install import get_elastix_bin, loader_env, try_elastix
 import subprocess
 root = Path(os.environ["KONFAI_ELASTIX_DIR"])
 try_elastix(root)
 r = subprocess.run([str(get_elastix_bin(root).resolve()), "--extended-version"],
                   capture_output=True, text=True, env=loader_env(root))
 assert r.returncode == 0, r.stderr
 print(r.stdout.splitlines()[0] if r.stdout else "elastix OK")
 print("elastix-impact shim OK:", root)
 PY
 echo
 echo "Done. The impact pipeline picks the shim up automatically"
 echo "(src/registration_pipelines.py::_impact_elastix_root), e.g.:"
 echo "  python scripts/bench_registration.py --pipelines impact --n-cases 2"
--- a/src/impact_loss.py
+++ b/src/impact_loss.py
@ -0,0 +1,152 @@
 # Vendored from https://github.com/vboussot/ImpactLoss (IMPACT.py)
 # Copyright (c) 2025 Valentin Boussot.  Apache-2.0.
 #
 # IMPACT Reg: A Generic Semantic Loss for Multimodal Image Registration
 # (Boussot et al., arXiv:2503.24121).  Modification relative to upstream:
 # model files are resolved from a local cache dir first (``model_cache`` /
 # the ``IMPACT_MODELS_DIR`` env var) before falling back to the Hugging
 # Face download, so offline hosts can use the already-cached
 # ``data/impact_models`` layout.
 import os
 import torch
 from huggingface_hub import hf_hub_download
 IMPACT_HF_REPO = "VBoussot/impact-torchscript-models"
 def _resolve_model(model_name, model_cache):
    """Local cache path if present, else the Hugging Face download."""
    if model_cache:
        p = os.path.join(model_cache, model_name)
        if os.path.exists(p):
            return p
    return hf_hub_download(repo_id=IMPACT_HF_REPO, filename=model_name,
                           repo_type="model", revision=None)  # nosec B615
 class IMPACTReg(torch.nn.Module):
    """
    IMPACT-Reg Loss Module
    ------------------
    A similarity metric based on pretrained TorchScript feature extractors.
    Computes a weighted L1 distance between feature maps of two images.
    Parameters
    ----------
    model_name : str
        Name of the TorchScript model (e.g. "TS/M730").
    shape : list[int]
        Expected spatial size. Use [0, 0, 0] for no resample.
    in_channels : int
        Number of input channels (before adapting to the model).
    weights : list[float]
        Per-layer feature weights.
    model_cache : str, optional
        Directory the model files are looked up in before downloading.
    """
    def __init__(
        self, model_name: str, shape: list[int] = [0, 0], in_channels: int = 1,
        weights: list[float] = [1], model_cache: str = None,
    ) -> None:
        super().__init__()
        if model_name is None:
            return
        self.in_channels = in_channels
        self.loss = torch.nn.L1Loss()
        self.weights = weights
        self.nb_layer = len(weights)
        # Load the TorchScript model (local cache first, HF download after)
        model_cache = model_cache or os.environ.get("IMPACT_MODELS_DIR", "")
        self.model_path = _resolve_model(model_name, model_cache)
        self.model: torch.nn.Module = torch.jit.load(self.model_path, map_location=torch.device("cpu"))  # nosec B614
        self.dim = len(shape)
        self.shape = shape if all(s > 0 for s in shape) else None
        dummy_input = torch.zeros((1, self.in_channels, *(self.shape if self.shape else [128] * self.dim))).to(0)
        try:
            out = self.model.to(0)(dummy_input, torch.tensor([self.nb_layer]))
            if not isinstance(out, (list, tuple)):
                raise TypeError(f"Expected model output to be a list or tuple, but got {type(out)}.")
            if self.nb_layer != len(out):
                raise ValueError(
                    f"Loss 'L1': mismatch between the number of weights "
                    f"({self.nb_layer}) and the number of model outputs "
                    f"({len(out)}). Each output must have a corresponding weight."
                )
        except Exception as e:
            msg = (
                f"[Model Sanity Check Failed]\n"
                f"Input shape attempted: {dummy_input.shape}\n"
                f"Error: {type(e).__name__}: {e}"
            )
            raise RuntimeError(msg) from e
        self.model = None
    def preprocessing(self, tensor: torch.Tensor) -> torch.Tensor:
        """
        Resize and adjust the tensor channels to match model expectations.
        Note:
            This method does not perform any intensity normalization.
            Input images must be provided in their raw (non-normalized) form,
            i.e., in their original intensity space (e.g., HU for CT, native
            values for MRI). Intensity normalization / standardization is the
            responsibility of the IMPACT models.
        """
        if self.shape is not None and not all(
            tensor.shape[-i - 1] == size for i, size in enumerate(reversed(self.shape[2:]))
        ):
            if tensor.dtype == torch.uint8:
                mode = "nearest"
            elif len(tensor.shape) < 4:
                mode = "bilinear"
            else:
                mode = "trilinear"
            tensor = torch.nn.functional.interpolate(
                tensor, mode=mode, size=tuple(self.shape), align_corners=False
            ).type(torch.float32)
        # Replicate channels if needed (e.g., 1 -> 3)
        if tensor.shape[1] != self.in_channels:
            tensor = tensor.repeat(tuple([1, self.in_channels] + [1 for _ in range(self.dim)]))
        return [tensor, torch.tensor([self.nb_layer]),
                torch.tensor([tensor.min(), tensor.max(), tensor.mean(), tensor.std()])]
    def _compute(self, output: torch.Tensor, target: torch.Tensor) -> torch.Tensor:
        """Compute the feature-based L1 loss between two tensors."""
        loss = torch.zeros(1, device=output.device, dtype=torch.float32)
        output = self.preprocessing(output)
        target = self.preprocessing(target)
        for weight, output_features, target_features in zip(
            self.weights, self.model(*output), self.model(*target)
        ):
            if weight == 0:
                continue
            loss = loss + weight * self.loss(output_features, target_features)
        return loss
    def forward(self, output: torch.Tensor, target: torch.Tensor) -> torch.Tensor:
        """
        Main forward method.
        If the model is 2D and inputs are 3D, compute slice-wise average loss.
        """
        if self.model is None:
            self.model = torch.jit.load(self.model_path)  # nosec B614
        self.model.to(output.device)
        self.model.eval()
        loss = torch.zeros(1, device=output.device, dtype=torch.float32)
        # 3D -> 2D model case: iterate over slices
        if len(output.shape) == 5 and self.dim == 2:
            for i in range(output.shape[2]):
                loss = loss + self._compute(output[:, :, i, ...], target[:, :, i, ...])
            loss = loss / output.shape[2]
        else:
            loss = self._compute(output, target)
        return loss
--- a/src/registration_pipelines.py
+++ b/src/registration_pipelines.py
Author	SHA1	Message	Date
Furen Xiao	bec0e9834f	feat(benchmark): add ANTs deformable variant pipelines Extend registration_pipelines and bench_registration with stock antspyx ANTs presets (SyNRA, SyN, Repro sr/br, QuickRepro sr/br) plus the legacy Rigid -> SyNOnly workflow, and update the deformable benchmark documentation with follow-up ANTs-variant runs and results.	2026-10-03 13:26:12 +08:00
Furen Xiao	aa7ac26f4c	Refactor code structure for improved readability and maintainability	2026-10-03 13:25:43 +08:00
Furen Xiao	f87e5c51af	Refactor code structure for improved readability and maintainability	2026-10-02 19:00:07 +08:00
Furen Xiao	96733e370d	docs(benchmark): add rigid registration pipeline benchmark documentation Document the unified benchmark contract, per-pipeline settings, and evaluation details for the six rigid CT-to-T1c registration pipelines.	2026-10-02 18:59:08 +08:00
Furen Xiao	77653e4018	refactor(impact): replace KonfAI pip runtime with Docker elastix_impact image Migrate register_impact from impact-reg-konfai (pip + one-time glibc/LibTorch shim) to the vboussot/elastix_impact Docker image running a rigid 6-DOF EulerTransform with MIND + TotalSegmentator static features and a Mattes-MI secondary metric. Delete scripts/setup_impact_konfai.sh. Vendor ImpactLoss (vboussot/ImpactLoss) into src/impact_loss.py and add impact_loss as a per-case semantic scoring metric in the benchmark runner. Update docs (AGENTS.md, REGISTRATION_BENCHMARK.md) to reflect the all-rigid pipeline set and the new Docker-based IMPACT runtime.	2026-09-30 21:23:18 +08:00