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_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).
+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`.
 
 ## Project
 

REGISTRATION_BENCHMARK.md

@@ -0,0 +1,726 @@
+# 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).

REGISTRATION_BENCHMARK_DEFORMABLE.md

@@ -1,456 +0,0 @@
-# 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.

REGISTRATION_BENCHMARK_RIGID.md

@@ -1,833 +0,0 @@
-# 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).

scripts/bench_registration.py

@@ -2,9 +2,6 @@
 
 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
@@ -14,10 +11,9 @@ 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, higher = better).
-  - impact_loss        : IMPACT semantic loss (vboussot/ImpactLoss; TotalSeg
-                         M730 deep features, weighted L1), warped-CT ~ MR over
-                         the full co-registered volumes (lower = better).
+                         (the reported alignment score).
+  - tv_in_ct_frac      : fraction of the (dilated) warped GTV that falls inside
+                         the warped-CT head tissue (the slab must straddle it).
   - tv_nvox            : warped GTV volume (>0.5 voxels).
   - skull_ring_gap_mm  : concentricity of the warped-CT skull ring vs the MR
                          head silhouette (mid/low level + angular std).
@@ -30,16 +26,15 @@ 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_RIGID.md), the benchmark is env-portable: run it under
+REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under
 whichever conda env has a given toolkit installed, and it records only the
 pipelines that are importable there.  elastix runs in the main env through the
 ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively
 (``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
-``impact`` runs the ``vboussot/elastix_impact`` Docker image (GPU; the rigid
-MIND + TotalSegmentator static config, needs Docker + the pulled image + the
-model assets — see REGISTRATION_BENCHMARK_RIGID.md); ``synthmorph`` runs in a Docker
-image (``freesurfer/synthmorph``; needs Docker + the pulled image).  The
-`--merge`
+``impact`` runs the ``impact-reg-konfai`` package (pip; GPU; one-time elastix
+shim setup on hosts too old for the stock binary — see
+REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker 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.
 
@@ -48,7 +43,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_RIGID.md
+  # synthmorph (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
   python scripts/bench_registration.py --pipelines synthmorph
 
   # any pipeline that is not importable in the current env is recorded as
@@ -104,54 +99,6 @@ 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)
@@ -206,12 +153,10 @@ 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 scores: Mattes MI between the CT
-    # inverse-transformed into MR coordinates (warped_ct) and the original MR
-    # (ha), over the MR head; IMPACT semantic loss over the full volumes.
+    # achieved-metric alignment score: Mattes MI between the CT inverse-transformed
+    # into MR coordinates (warped_ct) and the original MR (ha), over the MR head.
     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)
@@ -226,7 +171,6 @@ 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",
@@ -293,32 +237,6 @@ 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":
@@ -330,9 +248,8 @@ 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, gpu=args.impact_gpu,
-                                  iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS,
-                                  samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES)
+        return RP.register_impact(ct, t1c, preset=args.impact_preset, gpu=args.impact_gpu,
+                                  iters=args.impact_iters, samples=args.impact_samples)
     if pipeline == "niftyreg":
         return RP.register_niftyreg(ct, t1c, threads=args.threads,
                                     backend=args.niftyreg_backend,
@@ -351,7 +268,6 @@ 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
@@ -397,7 +313,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"impact_loss={r.get('impact_loss')} {reg}", flush=True)
+                      f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True)
             except Exception as e:  # noqa
                 errors.append({"key": c["key"], "error": repr(e)})
                 print(f"   {c['key']}: ERROR {e!r}", flush=True)
@@ -409,7 +325,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]
-    ilv = [r["impact_loss"] for r in registered if r.get("impact_loss") is not None]
+    tvfr = [r["tv_in_ct_frac"] for r in registered if r.get("tv_in_ct_frac") 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 = {
@@ -425,8 +341,7 @@ 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,
-        "impact_loss": {"mean": round(float(np.mean(ilv)), 4),
-                        "median": round(float(np.median(ilv)), 4)} if ilv else None,
+        "tv_in_ct_frac": {"mean": round(float(np.mean(tvfr)), 3)} if tvfr 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,
@@ -452,28 +367,28 @@ def _fmt(v, spec="{}", default="-"):
 
 
 def print_table(rows):
-    order = [p for p in RP.ALL_PIPELINE_NAMES if any(r["pipeline"] == p for r in rows)]
-    print("\n" + "=" * 96)
+    order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)]
+    print("\n" + "=" * 92)
     hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} "
-           f"{'s/reg':>7s} {'ImpactLoss':>10s} {'MI':>7s}")
+           f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
     print(hdr)
-    print("-" * 96)
+    print("-" * 92)
     for p in order:
         r = next(x for x in rows if x["pipeline"] == p)
         if not r.get("available"):
             print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} "
-                  f"{'-':>7s} {'-':>10s} {'-':>7s}   [{r['reason']}]")
+                  f"{'-':>7s} {'-':>7s} {'-':>7s}   [{r['reason']}]")
             continue
         rt = (r.get("reg_wall_s") or {}).get("median")
-        il = (r.get("impact_loss") or {}).get("median")
+        tv = (r.get("tv_in_ct_frac") or {}).get("mean")
         mi = (r.get("mi_ct_mr") or {}).get("median")
         print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} "
               f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} "
-              f"{_fmt(il, '{:.3f}'):>10s} {_fmt(mi, '{:.3f}'):>7s}")
-    print("=" * 96)
+              f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
+    print("=" * 92)
 
 
-def write_plot(rows, out_png, label="rigid"):
+def write_plot(rows, out_png):
     try:
         import matplotlib
         matplotlib.use("Agg")
@@ -487,19 +402,13 @@ def write_plot(rows, out_png, label="rigid"):
         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, 3, figsize=(14, 4))
-    colors = ["#6ACC65", "#F5A623", "#4878CF"]
+    fig, axes = plt.subplots(1, 2, figsize=(10, 4))
+    colors = ["#4878CF", "#6ACC65"]
     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")]):
-        # 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])
+        bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)])
         for b, v in zip(bars, vals):
             ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()),
                         ha="center", va="bottom", fontsize=8)
@@ -507,7 +416,7 @@ def write_plot(rows, out_png, label="rigid"):
         ax.set_ylabel(ylab)
         ax.set_xticks(range(len(vals)))
         ax.set_xticklabels(labels, fontsize=9)
-    plt.suptitle(f"CT->T1c registration benchmark (M6, cross-modality {label})", fontsize=11)
+    plt.suptitle("CT->T1c registration benchmark (M6, cross-modality rigid)", fontsize=11)
     fig.tight_layout()
     fig.savefig(out_png, dpi=120)
     plt.close(fig)
@@ -517,7 +426,7 @@ def write_plot(rows, out_png, label="rigid"):
 def merge(args):
     out_dir = os.path.join(args.out_dir)
     rows = []
-    for p in RP.ALL_PIPELINE_NAMES:
+    for p in RP.ALL_PIPELINES:
         fp = os.path.join(out_dir, p + ".json")
         if os.path.exists(fp):
             with open(fp) as f:
@@ -526,30 +435,13 @@ 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": protocol,
+        "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI "
+                      "(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,
@@ -558,10 +450,9 @@ 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"),
-               label="deformable" if deformable else "rigid")
+    write_plot(rows, os.path.join(out_dir, "bench_registration.png"))
     # markdown table for the report
-    md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |",
+    md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | tv-in-CT | MI (nats, median) |",
           "|---|---|---|---|---|---|---|---|---|"]
     for p in [x["pipeline"] for x in rows]:
         r = next(x for x in rows if x["pipeline"] == p)
@@ -572,7 +463,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('impact_loss') or {}).get('median')} | "
+            f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | "
             f"{(r.get('mi_ct_mr') or {}).get('median')} |")
     mf = os.path.join(out_dir, "summary.md")
     with open(mf, "w") as f:
@@ -582,10 +473,6 @@ 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)")
@@ -609,16 +496,24 @@ 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="IMPACT MaximumNumberOfIterations, one value for all "
-                         "levels (default: the reference map's per-level budget "
-                         f"{RP.IMPACT_DEFAULT_ITERS})")
+                    help="preset max_iterations override, one value for all "
+                         "resolutions (only effective for presets without a "
+                         "per-resolution model matrix; default: the preset's "
+                         "own budget, 400/300/200/200 for Elastix_IMPACT_Static)")
     ap.add_argument("--impact-samples", type=int, default=None,
-                    help="IMPACT NumberOfSpatialSamples per iteration (default: "
-                         f"{RP.IMPACT_DEFAULT_SAMPLES})")
+                    help="preset spatial_samples override (NumberOfSpatialSamples; "
+                         "default: the preset's own value, 2000)")
     ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"],
                     default="auto",
                     help="niftyreg reg_aladin build: docker = official v2.1.1 in the "
@@ -647,9 +542,6 @@ 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
@@ -658,7 +550,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_PIPELINE_NAMES)
+    unknown = set(pipelines) - set(RP.ALL_PIPELINES)
     if unknown:
         raise SystemExit(f"unknown pipeline(s): {unknown}")
 

scripts/investigate_synthmorph.py

@@ -1,601 +0,0 @@
-"""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()

scripts/investigate_synthmorph_mm.py

@@ -1,101 +0,0 @@
-"""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()

scripts/setup_impact_konfai.sh

@@ -0,0 +1,110 @@
+#!/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"

src/impact_loss.py

@@ -1,152 +0,0 @@
-# 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