longitudinal/REGISTRATION_BENCHMARK.md

# Registration Pipeline Benchmark — CT → T1c (M6)

Benchmarks **five 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**. All pipelines register
on the same CT-space head box (MI computed on head tissue, not neck/air) and,
for the single-init tools, from the **same centroid-aligned start**
(`centroid_offset`), so the comparison is on the optimizer / metric / transform
handling — not on who finds the ~30–60 mm gross offset first.

## The four pipelines

| # | pipeline | toolkit | metric | transform | status in this env |
|---|----------|---------|--------|-----------|--------------------|
| 1 | `elastix`   | ITK-Elastix 0.25.4 (pip `itk` + `itk-elastix`, elastix 5.3.1) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x replacement for `Rigid3DTransform`), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
| 2 | `simpleitk` | SimpleITK 2.5.6 | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
| 3 | `antspyx`   | ANTsPy / ANTs 0.6.3 | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
| 4 | `deepreg`   | Oxford DeepReg (CycleGAN) | CycleGAN + reg. loss (pretrained) | 3D displacement field | **not installed** (see below) |
| 5 | `impact`    | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT) | **IMPACT semantic** (TotalSegmentator `TS/M730.pt` features, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales (no image resampling) | ✅ ran (main env, Docker + GPU) |

`simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid,
centered start — the default, comparable to the other tools) and
`--simpleitk-mode production` (the full multi-start tilt grid from
`02_build_m6_dataset.py`, ~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)

Production QA, crop-independent:

- **`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. **Aligned ⇔
  `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`.** `mid` measures translational
  concentricity; `std` measures rotational concentricity.
- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT
  head tissue (the slab must straddle the lesion). Acceptance needs `≥ 0.5`.
- **`tv_nvox`** — warped GTV volume.
- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the
  production acceptance rule).
- **`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). This is the registration's *achieved* metric read
  off as a standalone alignment score: 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.
- **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring).

## Results (10 cases)

| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|---|---|---|---|---|---|---|---|---|
| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.6 | 0.93 | **0.643 / 0.630** |
| **antspyx** | ANTsPy 0.6.3 | 9 | **0/10 (0.00)** | 14.9 / 16.5 | 18.4 / 22.7 | 14.5 | 0.93 | **0.598 / 0.529** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | **5/10 (0.50)** | 15.3 / 16.1 | 9.9 / 13.8 | 3.2 | 0.93 | **0.562 / 0.549** |
| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 12.4 / 19.0 | 15.1 / 16.7 | **562.3** | 0.93 | **0.475 / 0.448** |
| deepreg | — | — | — | — | — | — | — | — |

*(antspyx and impact are stochastic (random metric/patch subsampling), so their
accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.4
min/case on GPU, orders of magnitude above the intensity-metric tools.)*

Per-case (mid / std, `OK` = accepted):

| case | simpleitk | antspyx | elastix | impact |
|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 24.4 / 27.7 | 26.4 / 30.3 | **8.4 / 12.6** |
| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.4 / 24.8 | 21.1 / 23.3 | 21.4 / 24.1 |
| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.5 / 18.4 | 10.6 / 8.2 **OK** | 10.0 / 8.3 **OK** |
| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 14.9 / 15.5 | 11.0 / 6.1 **OK** | 11.1 / 8.8 **OK** |
| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.3 / 17.5 | 15.3 / 14.0 | 54.4 / 15.1 |
| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 17.7 / 37.6 | 11.4 / 3.2 **OK** | 23.9 / 22.6 |
| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** | 10.1 / 5.9 **OK** |
| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.1 / 14.3 | 12.0 / 6.8 **OK** | 12.4 / 22.3 |
| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 17.2 / 35.4 | 18.9 / 22.7 | 19.6 / 30.3 |

### Interpretation

- **Translation is recovered well by all four** (mid gaps mostly 6–20 mm,
  within the `≤ 20 mm` band). `tv-in-CT` is identical (~0.93): the warped GTV
  lands inside the CT head for every pipeline, so the target is in the right
  place.
- **The binding failure is rotation** — `std > 12 mm` on the repositioned cases
  (the ~30° inter-session tilt produces local minima that a single centered
  start can lock onto). This is exactly why production `02_build_m6_dataset.py`
  uses a **multi-start tilt grid**.
- **Elastix is the most accurate of the standard single-init configs**
  (5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 18.4) and
  the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s
  SimpleITK / 14.5 s ANTs): ASGD's adaptive step sizes climb out of small
  misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223`
  (std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on
  (std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both
  others accept — different optimizers settle in different local minima on the
  hard tilted cases.
- **ANTs' Rigid is the weakest of the intensity-metric tools** (0/10 in this
  draw; it is stochastic, so a prior draw gave 1/10): it frequently locks onto
  the tilted local minimum. Note `m6_26TYQI3C`: ANTs nailed the position (mid
  17.7) but missed the orientation (std 37.6) — a textbook tilted local minimum.
- **IMPACT (semantic metric) is the standout on the hardest repositioned case**
  and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every
  intensity tool fails, std 28–30 mm) it gets **mid 8.4 / std 12.6** — the best
  translation of any pipeline and within 0.6 mm of the acceptance threshold —
  because matching *anatomical semantics* (TotalSegmentator features) is far
  more discriminating for gross pose than raw cross-modality intensity. But the
  semantic landscape is rugged: single-init ASGD lands in bad local minima on
  4 of 9 cases (std 15–30 mm), and one (`m6_24V62QQM_20240820`) even drifts the
  translation (mid 54 mm). Net 3/10 accepted — between ANTs and SimpleITK — at
  ~9.4 min/case, two orders of magnitude slower than the intensity tools.
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~3.2 s
  (warm workers; a worker's first case pays a ~10 s warm-up), SimpleITK ~8.6 s,
  ANTs ~14.5 s (CPU-bound), **IMPACT ~562 s** (GPU-bound: every ASGD iteration
  runs the TotalSegmentator extractor on the sampled patches with backprop).
- **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**:
  medians 0.643 (SimpleITK) / 0.598 (ANTs) / 0.562 (Elastix) / 0.475 (IMPACT),
  and it diverges sharply on the tilted cases — ANTs `m6_26TYQI3C` drops to
  ~0.1 nats (rotation ~38 mm off) while SimpleITK holds 0.670. A high MI means
  the CT skull/tissue intensities fall on the matching MR structures after the
  inverse transform. IMPACT's lower median MI is consistent with its less
  reliable single-init convergence.
- **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
  (RSGD wall time is thread-sensitive, convergence is not); the in-process
  Elastix ASGD is deterministic given its fixed seed. **IMPACT is the most
  stochastic**: its Jacobian mode draws random patches every iteration, so the
  final pose — and even *which* cases get accepted — changes run to run (two
  independent 10-case draws both accepted 3/10 but on different cases). Treat
  its per-case numbers as one sample from a distribution.

## Reproducing

```bash
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here); deepreg
# reports unavailable and is recorded as such
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
python scripts/bench_registration.py --n-cases 10 \
    --workers 4 --threads 24 --simpleitk-mode benchmark

# IMPACT (needs Docker + the pulled image + a model; see the Impact section
# below for the one-time setup)
python scripts/bench_registration.py --pipelines impact --n-cases 10 \
    --workers 4 --threads 24

# 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`, `--simpleitk-mode {benchmark,production}`,
`--deepreg-ckpt <.ckpt>` (DeepReg), `--impact-model <rel|abs>`,
`--impact-gpu <int|-1>`, `--impact-iters N`, `--impact-samples N`,
`--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, Mattes MI, and the skull-ring gap
(mid/std), so a glance tells you both the achieved-metric score and whether the
visuals agree. 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 --n-cases 10 --qa-images
python scripts/bench_registration.py \
    --pipelines impact --n-cases 10 --qa-images
```

(9 images each; the incomplete-CT case is skipped.)

## Enabling the blocked pipeline

**DeepReg (Oxford)** — the PyPI `deepreg` (0.1.2) is an unrelated stub that pins
`numpy==1.18.5`; do not install it. DeepReg is not a one-shot rigid tool: it is a
CycleGAN image-to-image registration that must be **trained or run from a
pretrained checkpoint**. To benchmark it:

```bash
git clone https://github.com/OB-CVHU/DeepReg && cd DeepReg && pip install -e .
python scripts/bench_registration.py --pipelines deepreg --n-cases 10 \
    --deepreg-ckpt /path/to/pretrained_3d.ckpt
```

The adapter runs `deepreg test` to predict a voxel displacement field, converts
it to a SimpleITK `DisplacementFieldTransform`, and scores it like the others.

## 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 **TotalSegmentator** feature extractor. It is the only pipeline that
is GPU-based and the only one that is Docker-based.

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/). M730.pt is the recommended TotalSegmentator.
mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/TS
curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/TS/M730.pt \
    https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/TS/M730.pt
```

Other models (M258/M291/M850, MIND, SAM2.1, …) 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.

### How the adapter drives the image

The image's `CMD.sh` only runs `elastix -f -m -p -out` (no `-t0`), so the
adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`)
to add the shared centered start:

1. **Head box** — resamples fixed=T1c and moving=CT into the same CT-space head
   box (physical frame) the other pipelines use, so MI/IMPACT is computed on
   head tissue, not neck/air.
2. **Shared init via `-t0`** — the `centroid_offset` translation is written to a
   `TranslationTransform` parameter file and passed with `-t0`; with
   `HowToCombineTransforms "Compose"` the final transform is
   `T = T_euler(T_init(p))`. (Elastix 5.x removed in-file multi-`Transform`
   chains and `UseScalingFactors`, so the initial transform must be a separate
   `-t0` file.)
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 box's base spacing) with
   **no image resampling** (all-1 pyramid schedules, as the IMPACT docs require),
   ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid
   (`--impact-gpu -1` forces CPU).
4. **Read-back** — the result `TransformParameters.0.itk.txt` is an
   `Euler3DTransform` in the box/world physical frame; the adapter composes it
   with the `-t0` translation and runs the usual `fix_orientation` gate.

## Notes / caveats

- The benchmark compares **standard single-init** configurations (centered
  start). 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.
- `deepreg` could not be executed here, so its row is "not installed", not
  measured; its code path is implemented and unit-verified for import/convention.
- `fix_orientation` (elastix/ants/deepreg/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. In
  practice it never fired (all `auto_inverted: False`): elastix/ants/impact
  returned the T1c→CT direction as expected.
- **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.
  - The initial transform (the shared `centroid_offset` start) is passed as an
    `InitialTransformParameterObject` (a `TranslationTransform` with
    `NumberOfParameters 3`); the C++ filter only wires it up when an
    **output directory is set** (`SetOutputDirectory`), so 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.
  - The result transform is recovered by point-mapping
    (`GetCombinationTransform().TransformPoint`) into an affine — independent
    of elastix's Euler-angle parameterization — and wrapped in a SimpleITK
    `AffineTransform(3)`.
  - Registration wall time is ~1.8 s/case at full budget (ASGD 500/500/250
    over 3 resolutions); each worker's first case pays a ~10 s warm-up, so the
    median wall is 3.2 s and the mean is inflated by warm-up.
- **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0,
  GPU, Jacobian mode):
  - **Slow, by design.** Every ASGD iteration runs the TotalSegmentator
    extractor on the sampled patches with backpropagation, so the recommended
    500-iter / 2000-sample budget is ~9 min/case (vs 3–15 s for the intensity
    tools). `--impact-iters` / `--impact-samples` trade accuracy for speed.
  - **Rugged single-init convergence.** The semantic landscape has more local
    minima than Mattes-MI: the pipeline nails the hardest repositioned case
    (`m6_22IRRKU6`, mid 8.4 mm — best of all) but lands in bad minima on some
    others. A multi-start or higher-iteration budget would likely raise its
    acceptance; that is out of scope for the single-init head-to-head.
  - **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 shared centered start is a separate `-t0` translation file,
    and the result is read back as an `Euler3DTransform` (`.itk.txt`) and
    composed with `-t0`. The image's own `CMD.sh` doesn't pass `-t0`, 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).