longitudinal/REGISTRATION_BENCHMARK.md
Furen Xiao 73234db39e fix(benchmark): lower IMPACT sample-validity gate to unblock M6 runs
The no-mask CT→T1c FOV mismatch left only ~14% of spatial samples
inside the moving buffer at auto-init, below Elastix's default
RequiredRatioOfValidSamples of 0.25, so the IMPACT pipeline aborted
on every M6 case. Set the gate to 0.1 so the metric's native
in-buffer evaluation proceeds while still catching diverging
transforms.

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

678 lines
No EOL
40 KiB
Markdown
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

# Registration Pipeline Benchmark — CT → T1c (M6)
Benchmarks **six rigid, cross-modality registration pipelines** that align a
head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be warped onto the
T1c (MR) native grid. This is the registration step of
`scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT
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).
- **CT is always the fixed image** (moving = T1c); when both inputs are MR
the first argument is kept as fixed.
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, so with the CT (head-neck, ~300×300×251 mm) as
the fixed image and the tight T1c head crop (~220×220×150 mm) as moving, most
of the CT buffer maps outside the T1c buffer (~14% of samples in-buffer at the
auto-init). Elastix's sample-validity gate (build default
`RequiredRatioOfValidSamples 0.25`) used to abort the M6 pairs on that ("too
many samples map outside moving image buffer"); the parameter map now sets
`RequiredRatioOfValidSamples 0.1` so the FOV-mismatched runs proceed — see the
IMPACT section. The comparison is on each toolkit's native optimizer / metric
/ global-search capability.
## The six pipelines
Listed alphabetically. All register on the **full volumes** (no head box, no
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 = CT, moving = T1c** (shared rule above; when both
inputs are MR the first argument is fixed). `niftyreg` and `synthmorph`
clip the CT to [0, 80] HU before registration (per the CT-clipping A/B
below; `clip_ct=None` disables it).
| # | pipeline | fixed / moving | init | metric | transform | status in this env |
|---|----------|----------------|------|--------|-----------|--------------------|
| 1 | `antspyx` | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
| 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (32 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
| 3 | `impact` | **CT / T1c** | auto CenterOfGravity | **IMPACT semantic** (MRSegmentator `MRSeg/MRSeg.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (Docker + GPU; sample gate 0.1 — see IMPACT section) |
| 4 | `niftyreg` | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (Docker `niftyreg:2.1.1-ubuntu`, GPU `-platf 1`) |
| 5 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
| 6 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
Benchmark budget (best-results settings): `elastix`, `impact` and `simpleitk`
run with **MaximumNumberOfIterations = 2000** (per resolution / per pyramid
level), with the sampler and metric left at their defaults
(**NumberOfSpatialSamples 5000**, **NumberOfHistogramBins 32**). `antspyx`,
`niftyreg` and `synthmorph` keep their native defaults.
`simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid,
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 = CT. `niftyreg` and `synthmorph` clip the CT to [0, 80] HU *before
registration* (per the CT-clipping A/B below; `clip_ct=None` disables it);
`antspyx`, `elastix` and `simpleitk` register the CT as-is. `niftyreg` is the
v2.1.1 Docker build on the GPU (`-platf 1`); `impact` cannot run on these
pairs under the no-mask rule (see the IMPACT section).
| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|---|---|---|---|---|---|
| **antspyx** | ANTsPy 0.6.3 | 9 | 36.6 | 0.93 | **0.653 / 0.644** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.517 / 0.520** |
| **impact** | Docker Elastix 5.2 + IMPACT (MRSegmentator) | 8 | 2534 | 0.93 | **0.591 / 0.575** |
| **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 20.8 | 0.93 | **0.292 / 0.338** |
| **simpleitk** | SimpleITK 2.5.6 | 9 | 10.6 | 0.93 | **0.651 / 0.641** |
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 24.5 | 0.93 | **0.277 / 0.328** |
Per-case MI (nats, unclipped score; **bold** = best for the case):
| case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
|---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | **0.649** | 0.472 | 0.635 | 0.292 | 0.642 | 0.277 |
| m6_23XIMO4M_20230222 | **0.655** | 0.416 | 0.478 | 0.632 | 0.652 | 0.616 |
| m6_24V62QQM_20230223 | **0.704** | 0.701 | 0.697 | 0.260 | 0.663 | 0.256 |
| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.556 | 0.493 | **0.602** | 0.484 |
| m6_24V62QQM_20240820 | **0.653** | 0.644 | 0.531 | 0.118 | 0.652 | 0.111 |
| m6_26TYQI3C_20221027 | 0.659 | 0.653 | gate-abort | 0.609 | **0.695** | 0.590 |
| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | 0.450 | 0.069 | **0.577** | 0.061 |
| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | 0.629 | 0.233 | 0.634 | 0.231 |
| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | 0.625 | 0.340 | 0.651 | 0.323 |
### Interpretation
- **SimpleITK and ANTsPy tie for best** (median MI 0.651 / 0.652): on the
well-positioned cases they cluster at 0.60–0.70 and both survive the harder
repositioned cases (0.58–0.65 on `m6_2BBN5XA4` / `m6_22IRRKU6`). SimpleITK's
per-case MI is **unchanged** from the previous 150-iteration budget — RSGD
converges well before the 2000-iteration cap, so the cap is now headroom
rather than a constraint.
- **The grossly repositioned cases define the ranking.** On
`m6_2BBN5XA4` and `m6_24V62QQM_20240820` (~30° inter-session tilt) the tools
split: ANTs/SimpleITK hold 0.58–0.65, elastix drops to 0.29–0.64, niftyreg
0.07–0.12, synthmorph 0.06–0.10 — single-start optimizers locking onto local
minima. This is exactly why production `02_build_m6_dataset.py` uses a
**multi-start tilt grid**.
- **Elastix is second at the new budget** (2000 iters per resolution, 32 bins,
5000 samples): median 0.516, below the 0.641 recorded at the old 500/500/250
+ 10-bin setting. The long ASGD run nails the clean cases (0.64–0.70) but its
single CenterOfGravity start still locks onto local minima on the
repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
0.416).
- **IMPACT (semantic metric) now runs on these pairs** — the no-mask abort was
elastix's sample-validity gate (build default `RequiredRatioOfValidSamples`
0.25) rejecting the ~14% in-buffer fraction of the CT head-neck FOV vs the
tight T1c crop; the parameter map sets it to 0.1 (see the IMPACT section).
With MRSegmentator at the 2000/5000 budget it is top-three in MI on 7 of the
8 registered cases (median 0.591): 0.635 on the hardest `m6_22IRRKU6`
(masked-era 0.624), 0.697 on `m6_24V62QQM_20230223` (near-best), 0.629 on
`m6_2C4TKI6H` — and it shares the production-acceptance of the grossly
repositioned `m6_2BBN5XA4` (std 10.7 mm) with antspyx/simpleitk, the only
non-intensity tool to do so. One run diverges mid-registration
(`m6_26TYQI3C`: the in-buffer fraction falls below 10% at iteration 256 and
the gate aborts — a stochastic-draw outcome). It remains by far the slowest
(~42 min/case at 2000/5000).
- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
is strong on moderate cases (`m6_23XIMO4M` 0.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): elastix 6.0 s,
SimpleITK 10.6 s, **NiftyReg 20.8 s** (v2.1.1, GPU block matching),
SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT is by far the
slowest — 2534 s/case median at the 2000/5000 benchmark budget, ~276 s/case
at the authors' 500/2000; GPU feature extraction with backprop per
iteration. No-mask runs are faster per iteration than the masked-era
~833 s/case @500/2000 because only the ~14% in-buffer samples are scored.)
- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
per-case numbers carry run-to-run variance (this table is one draw);
SimpleITK's MI (and registration) is deterministic at fixed thread count;
the in-process Elastix ASGD is deterministic given its fixed seed;
SynthMorph is deterministic; NiftyReg is deterministic on the CPU
(conda 1.5.69) but shows small run-to-run variation on the GPU (v2.1.1,
~±0.01–0.05 MI on the weak cases). **IMPACT is the most stochastic**:
its Jacobian mode draws random patches every iteration, so the final pose
changes run to run. Treat its per-case numbers as one sample from a
distribution.
## CT clipping A/B (`--clip-ct 0 80`)
Question: does clipping the CT to a brain window ([0, 80] HU — skull and air
removed) **before registration** improve the result? The harness clips the CT
before registration AND before scoring (percentile-binned MI is scale-free,
so the A/B compares like for like). Same 10 cases, same budgets.
Median achieved MI (nats) and wall-clock (9 complete-CT cases):
| pipeline | no clip | clip [0,80] | Δ (median MI) | clip wins (cases by MI) |
|---|---|---|---|---|
| antspyx | 0.653 / 36.6 s | 0.466 / 34.7 s | **−0.186** | 0/9 |
| elastix | 0.517 / 6.0 s | 0.425 / 5.9 s | **−0.092** | 1/9 |
| niftyreg | 0.296 / 23.8 s | 0.320 / 20.2 s | +0.024 | 5/9 |
| simpleitk | 0.651 / 10.6 s | 0.436 / 11.1 s | **−0.215** | 0/9 |
| synthmorph | 0.241 / 26.0 s | 0.318 / 25.6 s | +0.077 | 5/9 |
| impact | abort | abort | - | - (historical: at the time both aborted on the old 0.25 sample-validity gate; the gate is now 0.1, A/B to be re-done) |
**Verdict: no — clipping the CT to [0, 80] HU does not improve the result; do
not clip.**
- The headline alignment score drops for **all three leading intensity
pipelines**: antspyx 0.653→0.466, simpleitk 0.651→0.436, elastix
0.517→0.425 (median nats; simpleitk and antspyx lose on **all 9** cases,
elastix on 8/9).
- The two weaker pipelines show small median-MI gains (niftyreg +0.024,
synthmorph +0.077) but only mixed case-by-case (5/9), and their geometric
QA does not support a real improvement.
- **Geometric concentricity (skull-ring gap) is a wash**: mean mid/std move by
~1–5 mm in either direction (niftyreg mid 23.6→18.1 the only clear gain;
antspyx std 16.4→18.0 the clear loss). Production acceptance even falls
(no-clip 4/4/2/4/2 → clip 1/0/1/0/2).
- **Mechanism.** Clipping is a saturating (non-invertible) intensity map: the
skull atom (800+ HU → 80) removes the strongest intensity anchor
(skull ↔ MR scalp) and compresses the CT rank distribution, so the
percentile-binned MI between the warped CT and the MR is inherently lower
*for the same pose*. The pose itself does not systematically degrade — the
head is rigid, so skull-ring concentricity is its fingerprint, and that is
unchanged. The MI drop is therefore mostly a score-scale effect, but it is
exactly the score the benchmark reports, and there is no geometric benefit
to pay it for.
- Scoring note: `_skull_ring_gap` tests the ring with `>= p97` (not `>`) so
the orientation gate and the gap QA also work on a clipped CT (where p97
sits exactly on the clipped-skull atom).
- **Decision / adopted default.** The pipelines were updated: `niftyreg` and
`synthmorph` clip the CT to [0, 80] HU before registration
(`clip_ct=(0.0, 80.0)` default, `clip_ct=None` disables it; MR inputs are
never clipped); `antspyx` / `elastix` / `simpleitk` register the CT as-is.
Re-checked on the **unclipped score scale** (the protocol used by the table
above): synthmorph keeps a genuine gain (0.241→0.277 median, 7/9 cases),
niftyreg is neutral (0.296→0.292; its apparent +0.024 above is mostly the
clip-scoring effect). All rows of the main table are scored on the
unclipped warped CT, so pipelines remain comparable.
## Reproducing
```bash
# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1)
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg \
--n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
# force the conda 1.5.69 CPU build: --niftyreg-backend conda
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
# sections below for the one-time setup)
python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \
--workers 4 --threads 24
# NOTE: impact now runs on the M6 pairs (the parameter map relaxes the
# sample-validity gate to RequiredRatioOfValidSamples 0.1) but is by far the
# slowest: ~5 min/case at 500/2000, ~42 min/case (2534 s median) at 2000/5000
# (--impact-iters / --impact-samples) — see the IMPACT section
python scripts/bench_registration.py --pipelines impact --n-cases 10 \
--workers 4 --threads 24
# A/B: does clipping the CT to a brain window help? (--clip-ct clips the CT
# for ALL pipelines, before registration and scoring; niftyreg/synthmorph
# already clip by default, so this only affects antspyx/elastix/simpleitk)
python scripts/bench_registration.py --n-cases 10 --workers 4 --threads 24 \
--out-dir results/bench_reg_noclip
python scripts/bench_registration.py --n-cases 10 --workers 4 --threads 24 \
--clip-ct 0 80 --out-dir results/bench_reg_clip0_80
python scripts/bench_registration.py --merge --out-dir results/bench_reg_noclip
python scripts/bench_registration.py --merge --out-dir results/bench_reg_clip0_80
# combine the per-pipeline JSONs accumulated across envs into one table + plot
python scripts/bench_registration.py --merge
```
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-model <rel|abs>`,
`--impact-gpu <int|-1>`, `--impact-iters N` (default 2000),
`--impact-samples N` (default 5000),
`--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 `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in shared world
space; with fixed=CT/moving=T1c the `-T` matrix is the T1c→CT contract
transform) 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 CT -flo T1c -rigOnly`
(symmetric variant on by default, which matters: the non-symmetric mode
drifts grossly on these repositioned pairs). `-ln 3` pyramid,
`-maxit 15` TLS iterations per level (the tool's default of 5 is too
coarse), `-interp 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` (CT→T1c); the adapter parses
it into a SimpleITK `AffineTransform` and runs it through the usual
`fix_orientation` gate (skull-ring concentricity). On the M6 cases the gate
keeps the file's affine as-is (`auto_inverted: False`), i.e. it is already
the contract transform (T1c→CT).
## IMPACT pipeline (Docker) — one-time setup
The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121)
*semantic* similarity metric compiled in as a native Elastix metric, plus a
TorchScript segmentation feature extractor (**MRSegmentator** `MRSeg/MRSeg.pt`
by default; TotalSegmentator `TS/M730.pt` and other models from
`VBoussot/impact-torchscript-models` are selectable). It is GPU-based and
Docker-based (alongside `synthmorph`). It follows the shared rules (full
volumes, no mask, fixed = CT, moving = T1c). The metric only evaluates samples
that map inside the moving buffer, and the CT head-neck FOV (~300×300×251 mm)
is much larger than the T1c head crop (~220×220×150 mm), so ~14% of the samples
are in-buffer at the auto-init. The image build's sample-validity gate
defaults to `RequiredRatioOfValidSamples 0.25` and used to abort the M6 pairs
on that ("too many samples map outside moving image buffer"); the parameter
map now sets it to `0.1` (see the IMPACT specifics below).
Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command):
```bash
# 1. Docker (with NVIDIA container toolkit so --gpus works)
docker --version
# 2. Pull the image (~23.5 GB, one-time)
docker pull vboussot/elastix_impact
# 3. Fetch the TorchScript model (~123 MB, one-time) into the repo's model
# cache (data/impact_models/). MRSeg.pt (MRSegmentator) is the default:
# it is trained on both MRI and CT, the right prior for this
# cross-modality CT->T1c task (all models share the same 2^l+3 FOV
# architecture, so the parameter map is unchanged between them).
mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg
curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg/MRSeg.pt \
https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/MRSeg/MRSeg.pt
```
Other models (TotalSegmentator `TS/M730.pt`, `TS/M258/M291/M850`, MIND, SAM2.1,
Anatomix, …) are selectable with `--impact-model <HF-relative-path-or-abs>`;
see the
[impact-torchscript-models](https://huggingface.co/VBoussot/impact-torchscript-models)
repo for the full list and FOV/layer tables. The model choice does not affect
the sample-validity gate (below): the "too many samples outside moving buffer"
check is geometric (CT head-neck FOV vs the tight T1c crop) and fires at the
first sample evaluation, before feature extraction.
### How the adapter drives the image
The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
binary directly (`/lib/elastix-install/bin/elastix`):
1. **Full volumes, fixed=CT, moving=T1c** — both volumes are used as-is (no
head box, no mask, no crop / pad, per the shared rules above). The metric
only evaluates spatial samples that map inside the moving buffer, and the
CT FOV (~300×300×251 mm, head+neck) is much larger than the tight T1c head
crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
the samples are in-buffer. The image build's default validity gate
(`RequiredRatioOfValidSamples 0.25`) aborted on that ("Too many samples map
outside moving image buffer: 704 / 5000"); the parameter map therefore sets
`(RequiredRatioOfValidSamples 0.1)` so the M6 pairs run on their in-buffer
samples while a genuinely diverging transform (<10% in-buffer) still
aborts.
2. **Auto init (no `-t0`)** — the parameter map sets
`AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
optimizer scales.
3. **Parameter map** — rigid `EulerTransform` (6-DOF), `Impact` metric in
**Jacobian** mode (differentiable, per the recommended config), 4 semantic
scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the T1c base spacing) with
**no image resampling** (all-1 pyramid schedules, as the IMPACT docs require),
ASGD, `Random` sampler, `(RequiredRatioOfValidSamples 0.1)`. A copy of the
exact parameter map is dumped to
`data/m6reg/impact/paramdump/ParameterMap_<model>_<pid>.txt` before every
docker run (also recorded as `param_dump` in the result JSON). GPU
auto-spread across the host GPUs by pid (`--impact-gpu -1` forces CPU).
4. **Intensity standardization happens inside the metric** (no pipeline-side
normalization — the model card's "standardize intensities to zero mean and
unit variance" is already applied): at init the metric computes
whole-volume `[min, max, mean, sigma]` for each image
(`itk::StatisticsImageFilter`) and passes it to the model, whose
`normalize` branch rescales every extracted patch with the image's own
`(mean, sigma)`. It is computed over the whole buffer (CT air included) and
identically for fixed and moving, so both modalities' features sit on the
same normalized scale. Pre-standardizing the MHA files in the pipeline
would be an exact no-op (z-scoring is affine-invariant: the metric
recomputes the stats of the image it receives) and would only corrupt the
HU scale used by MI scoring and the QA brain-window images.
5. **Read-back** — the result `TransformParameters.0.itk.txt` is an
`Euler3DTransform` in the world physical frame (the
`AutomaticTransformInitialization` is folded in by elastix); the adapter
runs it through the usual `fix_orientation` gate.
## 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/impact) 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. With the fixed=CT/moving=T1c swap, `elastix`'s combination
transform is now the inverse of the contract, so it reports
`auto_inverted: True` (the gate corrects it); antspyx/synthmorph/impact
return the contract direction directly (`auto_inverted: False`). `niftyreg`'s
`-aff` file is documented as CT→T1c but the gate keeps it as-is on these
cases (`auto_inverted: False`) — the empirical gate, not the docs, decides.
- **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU):
- **Cross-modality / tilt weakness.** The `rigid` model is trained on brain
*MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly
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 (~834 s).
- **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`,
symmetric block matching + TLS, CPU/OpenMP):
- **Different algorithmic family.** No histogram MI objective: 4-voxel
blocks are matched by normalised cross-correlation and the rigid
transform is solved by a trimmed least squares, iterated per pyramid
level (`-ln 3 -maxit 15`). NCC's scale invariance is what makes it
cross-modality tolerant.
- **Symmetric mode matters.** The default symmetric variant (optimising
forward and backward transforms concurrently) is the robust one on these
repositioned pairs; `-noSym` drifts grossly (gap std 45 mm on a clean
case every other tool passes).
- **`-maxit 5` (the tool default) is too coarse** — on the clean case
`m6_26TYQI3C` it reaches MI 0.41 vs 0.61 at 15 TLS iterations/level;
15 is the benchmark setting.
- **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=CT/moving=T1c the combination transform is CT→T1c (the inverse
of the contract), so `fix_orientation` reports `auto_inverted: True` and
inverts it. 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.0,
GPU, Jacobian mode):
- **Slow, by design.** Every ASGD iteration runs the segmentation
extractor (MRSegmentator by default) on the sampled patches with
backpropagation. The benchmark
budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
the elastix default) is ~10× the authors' 500-iter / 2000-sample
Jacobian budget: measured 2534 s/case (median) on the no-mask M6
pairs, ~276 s/case at 500/2000. `--impact-iters` / `--impact-samples`
trade accuracy for speed.
- **Runs on M6 with a relaxed sample-validity gate.** Like every other
pipeline it runs fixed=CT/moving=T1c on the full volumes (no mask, no
crop / pad). The CT FOV (~300×300×251 mm) is much larger than the T1c
crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
the samples map inside the T1c moving buffer. The metric calls
`CheckNumberOfSamples()` after every evaluation and the image build
defaults to `RequiredRatioOfValidSamples 0.25`, so it used to abort
("Too many samples map outside moving image buffer: 704 / 5000"). The
parameter map now sets `(RequiredRatioOfValidSamples 0.1)`: the gate
still aborts a genuinely diverging transform (<10% in-buffer) while the
FOV-mismatched M6 pairs run on their ~14% in-buffer samples — exactly
the native "only in-buffer samples are evaluated" behaviour. In the
2000/5000 re-run the gate still fired once, mid-registration: on
`m6_26TYQI3C` the transform drifted until only 468/5000 samples (9.4%)
stayed in-buffer at iteration 256 and the run aborted. The transform
contract (T1c→CT) holds; `fix_orientation` decides the direction as
usual (all 8 completed runs report `auto_inverted: True`).
- **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns
the intensity centroids; no manual `-t0` translation is used.
- **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 show the same pattern: strong MI (0.635 on
`m6_22IRRKU6` at 2000/5000) but the pose still misses the concentricity
acceptance rule there (mid 26.1 / std 28.6); it *is* accepted on
`m6_2BBN5XA4` (mid 15.1 / std 10.7).
- **Stochastic.** Jacobian mode resamples patches each iteration, so results
vary run to run (see Stochasticity).
- **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image
bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors`
are gone. The result is read back as an `Euler3DTransform` (`.itk.txt`).
The image's own `CMD.sh` is a generic launcher, so the adapter invokes
`/lib/elastix-install/bin/elastix` directly.
- **Models are not in the image.** The TorchScript model(s) must be present
under the mounted `/Data/Models` (the repo caches them in
`data/impact_models/`, bind-mounted read-only).
- **CPU fallback:** `--impact-gpu -1` runs without `--gpus` and with
`ImpactGPU -1` / mixed precision off (much slower).