xfr/longitudinal
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

feat(pipelines): integrate niftyreg and adopt Mattes MI as the benchmark headline score

Browse source 
Add `register_niftyreg` as the sixth rigid CT→T1c pipeline using
`reg_aladin` (symmetric block matching + trimmed least squares, NCC
similarity). Runs natively via `conda install -c conda-forge niftyreg`;
no Docker required.

Change the benchmark's primary reported metric from geometric skull-ring
QA to the achieved Mattes MI (`mi_ct_mr`, nats) between the warped CT
and the original MR over the MR head. Geometric QA fields are still
recorded per case in the JSON output.

Raise the SimpleITK benchmark-mode optimizer budget from 150 to 2000
iterations per pyramid level so all toolkits are compared on an equal
iteration footing.

Update `REGISTRATION_BENCHMARK.md` (five→six pipelines, niftyreg row,
elastix bin count 10→32), `AGENTS.md` (toolkit list), and
`bench_registration.py` docstrings/usage accordingly.
This commit is contained in:
Furen Xiao 2026-09-29 09:47:28 +08:00
parent 8c083abe92
commit a65189cd33
4 changed files with 408 additions and 218 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
AGENTS.md
View file

@ -10,7 +10,7 @@ Activate with:
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 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. 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`), plus Docker images `freesurfer/synthmorph` and `vboussot/elastix_impact`.
## Project ## Project

338
REGISTRATION_BENCHMARK.md
View file

@ -1,6 +1,6 @@
# Registration Pipeline Benchmark — CT → T1c (M6) # Registration Pipeline Benchmark — CT → T1c (M6)
Benchmarks **five rigid, cross-modality registration pipelines** that align a 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 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 T1c (MR) native grid. This is the registration step of
`scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT `scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT
@ -27,30 +27,39 @@ register on the **full CT / T1c volumes** (no head box) and use each toolkit's
**native automatic transform initialization** (no manual centroid offset): **native automatic transform initialization** (no manual centroid offset):
elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity), elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
antspyx its default center-of-mass alignment, simpleitk starts from identity, antspyx its default center-of-mass alignment, simpleitk starts from identity,
synthmorph needs none. All use **fixed = CT, moving = T1c** except `impact`, niftyreg (reg_aladin) initializes from the image centres, synthmorph needs
which keeps **fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much none. All use **fixed = CT, moving = T1c** except `impact`, which keeps
larger than the T1c crop (~220×220×150 mm), so with CT as the fixed image too **fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than
many CT samples map outside the T1c buffer and the IMPACT metric aborts the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
("too many samples map outside moving image buffer"). The comparison is on samples map outside the T1c buffer and the IMPACT metric aborts ("too many
each toolkit's native optimizer / metric / global-search capability. samples map outside moving image buffer"). The comparison is on each
toolkit's native optimizer / metric / global-search capability.
## The five pipelines ## The six pipelines
Listed alphabetically. All register on the **full volumes** (no head box) and Listed alphabetically. All register on the **full volumes** (no head box) and
use each toolkit's **native automatic transform initialization** (no manual use each toolkit's **native automatic transform initialization** (no manual
centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization` centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk (CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
starts from identity, and synthmorph needs none (global deep method). All use starts from identity, niftyreg initializes from the image centres, and
**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c, moving synthmorph needs none (global deep method). All use **fixed = CT, moving =
= CT** (swapping breaks its TotalSegmentator metric — see below). T1c** except `impact`, which keeps **fixed = T1c, moving = CT** (swapping
breaks its TotalSegmentator metric — see below).
| # | pipeline | fixed / moving | init | metric | transform | status in this env | | # | 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) | | 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 (10 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) | | 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (32 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
| 3 | `impact` | **T1c / CT** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) | | 3 | `impact` | **T1c / CT** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) |
| 4 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | | 4 | `niftyreg` | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (main env, conda-forge `reg_aladin` 1.5.69) |
| 5 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) | | 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, `simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid,
centered start — the default, comparable to the other tools) and centered start — the default, comparable to the other tools) and
@ -74,128 +83,123 @@ m6_2BBN5XA4_20250411 m6_2C4TKI6H_20240506 m6_2CDNN7BH_20240102
## Metrics (identical for every pipeline) ## Metrics (identical for every pipeline)
Production QA, crop-independent: The headline alignment score is **`mi_ct_mr`** — the registration's *achieved*
Mattes MI, read off as a standalone score:
- **`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 - **`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**, 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 restricted to the MR head and each image percentile-binned independently
(SimpleITK Mattes-MI style). This is the registration's *achieved* metric read (SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are
off as a standalone alignment score: higher = the warped-CT and MR intensities better in-sync = better aligned; it collapses toward 0 when the CT is
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
rotated/translated away from the MR. 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). - **`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 ⇔ `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`).
## Results (10 cases) ## 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) | Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
|---|---|---|---|---|---|---|---|---| the MR head; higher = better).
| **antspyx** | ANTsPy 0.6.3 | 9 | **4/10 (0.40)** | 12.8 / 15.6 | 15.6 / 16.5 | 36.7 | 0.93 | **0.652 / 0.644** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 13.6 / 15.7 | 14.8 / 17.0 | 2.6 | 0.93 | **0.641 / 0.640** |
| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **4/10 (0.40)** | 11.6 / 11.2 | 16.2 / 17.3 | **833.5** | 0.93 | **0.511 / 0.493** |
| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 14.4 / 16.6 | 17.9 / 17.5 | 9.7 | 0.93 | **0.651 / 0.641** |
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 21.8 / 21.5 | 25.0 | 0.93 | **0.238 / 0.308** |
*(antspyx and impact are stochastic (random metric/patch subsampling), so their | pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
accepted counts are one draw each — see Interpretation. impact's s/reg is ~14
min/case on GPU (full volumes), orders of magnitude above the intensity-metric
tools. elastix, simpleitk and synthmorph are deterministic at fixed threads.)*
Per-case (mid / std, `OK` = accepted):
| case | antspyx | elastix | impact | simpleitk | synthmorph |
|---|---|---|---|---|---| |---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** |
| m6_22IRRKU6_20240115 | 24.6 / 28.6 | 24.0 / 28.3 | 22.7 / 27.0 | 26.3 / 29.6 | 10.4 / 23.7 | | **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** |
| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 22.3 / 23.5 | 9.1 / 16.2 | 21.9 / 23.2 | 20.3 / 23.3 | | **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** |
| m6_24V62QQM_20230223 | 9.5 / 15.6 | 11.6 / 14.8 | 8.6 / 7.2 **OK** | 14.4 / 23.2 | 11.1 / 21.8 | | **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** |
| m6_24V62QQM_20231017 | 11.7 / 7.3 **OK** | 12.2 / 7.4 **OK** | 4.6 / 6.6 **OK** | 11.7 / 7.5 **OK** | 11.5 / 8.4 **OK** | | **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** |
| m6_24V62QQM_20240820 | 10.8 / 9.1 **OK** | 11.4 / 9.0 **OK** | 11.6 / 4.9 **OK** | 11.5 / 9.1 **OK** | 43.4 / 15.8 | | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** |
| m6_26TYQI3C_20221027 | 12.3 / 4.4 **OK** | 10.6 / 3.1 **OK** | 11.9 / 4.5 **OK** | 11.3 / 4.6 **OK** | 11.8 / 2.9 **OK** |
| m6_2BBN5XA4_20250411 | 19.1 / 7.7 **OK** | 18.2 / 14.1 | −1.7 / 29.3 | 18.8 / 7.5 **OK** | 52.4 / 53.0 | \* `impact` at its **previous** 500-iter / 2000-sample budget; the new default
| m6_2C4TKI6H_20240506 | 12.8 / 17.2 | 13.6 / 20.0 | 15.4 / 23.2 | 14.0 / 17.9 | 16.7 / 28.2 | (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
| m6_2CDNN7BH_20240102 | 17.5 / 34.6 | 17.6 / 32.8 | 18.9 / 37.1 | 19.3 / 35.0 | 18.9 / 16.8 | (hours/case) and was not re-run.
Per-case MI (nats; **bold** = best for the case):
| case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
|---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 0.649 | 0.472 | **0.624** | 0.296 | 0.642 | 0.235 |
| m6_23XIMO4M_20230222 | 0.655 | 0.416 | 0.474 | **0.634** | 0.652 | 0.510 |
| m6_24V62QQM_20230223 | 0.703 | **0.701** | 0.554 | 0.260 | 0.663 | 0.238 |
| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.419 | 0.493 | **0.602** | 0.516 |
| m6_24V62QQM_20240820 | 0.652 | 0.643 | 0.458 | 0.118 | **0.652** | 0.096 |
| m6_26TYQI3C_20221027 | 0.657 | 0.653 | 0.511 | 0.608 | **0.695** | 0.604 |
| m6_2BBN5XA4_20250411 | 0.576 | 0.287 | 0.155 | 0.069 | **0.577** | 0.060 |
| m6_2C4TKI6H_20240506 | 0.645 | 0.401 | **0.626** | 0.231 | 0.634 | 0.209 |
| m6_2CDNN7BH_20240102 | 0.655 | 0.516 | **0.618** | 0.344 | 0.651 | 0.301 |
### Interpretation ### Interpretation
- **Translation is recovered well by all five** (mid gaps mostly 4–16 mm, - **SimpleITK and ANTsPy tie for best** (median MI 0.651 / 0.652): on the
within the `≤ 20 mm` band, except the grossly-drifted cases). well-positioned cases they cluster at 0.60–0.70 and both survive the harder
`tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for repositioned cases (0.58–0.65 on `m6_2BBN5XA4` / `m6_22IRRKU6`). SimpleITK's
every pipeline, so the target is in the right place. per-case MI is **unchanged** from the previous 150-iteration budget — RSGD
- **The binding failure is rotation, and removing the head box makes it converges well before the 2000-iteration cap, so the cap is now headroom
harder.** On the repositioned cases `std > 12 mm` — the ~30° inter-session rather than a constraint.
tilt produces local minima a single auto-start locks onto, and without the - **The grossly repositioned cases define the ranking.** On
head-box crop the full-volume metric sees more neck/air confounders. Notably, `m6_2BBN5XA4` and `m6_24V62QQM_20240820` (~30° inter-session tilt) the tools
the hardest tilt case `m6_22IRRKU6` is now **rejected by all five** (std split: ANTs/SimpleITK hold 0.58–0.65, elastix drops to 0.29–0.64, niftyreg
23–30 mm), which the head-box config got through with IMPACT. This is exactly 0.07–0.12, synthmorph 0.06–0.10 — single-start optimizers locking onto local
why production `02_build_m6_dataset.py` uses a **multi-start tilt grid**. minima. This is exactly why production `02_build_m6_dataset.py` uses a
- **The easy cases are now accepted by nearly all pipelines.** **multi-start tilt grid**.
`m6_24V62QQM_20231017`, `m6_24V62QQM_20240820` and `m6_26TYQI3C` are accepted - **Elastix is second at the new budget** (2000 iters per resolution, 32 bins,
by ANTs/Elastix/IMPACT/SimpleITK (SynthMorph accepts the two cleanest). The 5000 samples): median 0.516, below the 0.641 recorded at the old 500/500/250
native automatic inits (center-of-gravity / center-of-mass) place the heads + 10-bin setting. The long ASGD run nails the clean cases (0.64–0.70) but its
close enough that the optimizers converge on the well-positioned cases. single CenterOfGravity start still locks onto local minima on the
- **SimpleITK, ANTs and IMPACT tie at 4/10; Elastix 3/10; SynthMorph 2/10.** repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
ANTs is no longer the weakest — on the full volumes its native center-of-mass 0.416).
init works well (median std 15.6, up from the head-box+manual-pre-shift - **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI
31.7). It is stochastic (random MI subsampling), so the count varies by draw. on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618),
- **IMPACT (semantic metric) is 4/10 but the slowest (~14 min/case on full where anatomical-segmentation similarity out-resolves intensity MI — but its
volumes) and lost its standout on the hardest tilt** (`m6_22IRRKU6` now median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is
mid 22.7 / std 27.0 vs mid 8.3 / std 13.2 with the head box). Matching also by far the slowest (~14 min/case at the shown budget; the new 2000/5000
*anatomical semantics* still gives it the best median mid (11.6) and a budget would be on the order of hours/case).
competitive MI (0.511), but the full-volume registration is a larger problem - **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
for the local ASGD. It also drifts grossly on `m6_2BBN5XA4` (mid −1.7, the is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class;
transform flips to the far side). `m6_26TYQI3C` 0.608) but falls into the same tilt local minima
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.6 s (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower
(warm workers), SimpleITK ~9.7 s, ANTs ~36.7 s (CPU-bound, full volumes), than elastix at full volumes.
**SynthMorph ~25 s** (GPU: a single pretrained-model forward pass per case), - **SynthMorph is weakest of the six on these cross-modality CT pairs**
**IMPACT ~834 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator (median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30°
extractor on the sampled patches with backprop). tilt defeats it on the repositioned cases (0.06–0.10) even though it matches
- **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**: the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no
medians 0.651 (SimpleITK) / 0.641 (Elastix) / 0.511 (IMPACT) / 0.652 (ANTs) / centroid init** (global deep method) and is deterministic.
0.238 (SynthMorph). - **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the
The intensity-metric tools cluster at ~0.64 (a genuine alignment signal), CT head everywhere, so coverage is not the discriminator — MI is.
while SynthMorph's 0.238 reflects its tilt failures. A high MI means the CT - **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s,
skull/tissue intensities fall on the matching MR structures after the inverse SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s,
transform. **NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous
- **SynthMorph (deep, pretrained) is near-weakest of the five on these budget; GPU feature extraction with backprop per iteration).
cross-modality CT pairs.** Its rigid model is trained on brain *MRI*
(T1/T2, anatomy-aware but intensity-domain MRI), so it recovers gross
translation well (mid 10–20 mm on most cases) but is repeatedly defeated by
the ~30° inter-session tilt: 7 of 9 cases have `std > 15 mm` (rotation error)
and two drift grossly (`m6_2BBN5XA4` mid 52 / std 53, `m6_24V62QQM_20240820`
mid 43). It accepts the clean `m6_26TYQI3C` (mid 11.8 / std 2.9, MI 0.60) and
`m6_24V62QQM_20231017` (mid 11.5 / std 8.4) — where it is competitive with the
intensity tools. Notably it needs **no centroid init** (it is a global deep
method that recovers the gross offset on its own). Its low median MI (0.238)
reflects the tilt failures. On same-modality MRI it would likely rank much
higher; this CT benchmark is an out-of-distribution stress test for it. It is
deterministic (one fixed pretrained model, no random subsampling), so unlike
ANTs/IMPACT the per-case numbers are stable run to run.
- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
per-case numbers carry run-to-run variance (this table is one draw); per-case numbers carry run-to-run variance (this table is one draw);
SimpleITK's MI (and registration) is deterministic at fixed thread count SimpleITK's MI (and registration) is deterministic at fixed thread count;
(RSGD wall time is thread-sensitive, convergence is not); the in-process the in-process Elastix ASGD is deterministic given its fixed seed;
Elastix ASGD is deterministic given its fixed seed. **IMPACT is the most SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**:
stochastic**: its Jacobian mode draws random patches every iteration, so the its Jacobian mode draws random patches every iteration, so the final pose
final pose — and even *which* cases get accepted — changes run to run (two changes run to run. Treat its per-case numbers as one sample from a
independent 10-case draws both accepted 3/10 but on different cases). Treat distribution.
its per-case numbers as one sample from a distribution.
## Reproducing ## Reproducing
```bash ```bash
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here) # main env — SimpleITK + Elastix + ANTsPy + NiftyReg (pip/conda-installed here)
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx \ python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg \
--n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark --n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT # Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
@ -211,10 +215,12 @@ python scripts/bench_registration.py --merge
Other flags: `--sample random --seed 0` (reproducible random 10), Other flags: `--sample random --seed 0` (reproducible random 10),
`--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`, `--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`,
`--threads T`, `--simpleitk-mode {benchmark,production}`, `--threads T` (also the OpenMP thread count for `reg_aladin`),
`--simpleitk-mode {benchmark,production}`,
`--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`, `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`,
`--synthmorph-extent <int>`, `--impact-model <rel|abs>`, `--synthmorph-extent <int>`, `--impact-model <rel|abs>`,
`--impact-gpu <int|-1>`, `--impact-iters N`, `--impact-samples N`, `--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). `--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
## Visual QA images (checkerboard + color overlay) ## Visual QA images (checkerboard + color overlay)
@ -236,9 +242,9 @@ 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 white); the MR uses a 2–98 % head window. The coronal and sagittal views are
oriented superior-up to match the axial view. oriented superior-up to match the axial view.
The title carries the pipeline, case key, Mattes MI, and the skull-ring gap The title carries the pipeline, case key and the achieved Mattes MI, so a
(mid/std), so a glance tells you both the achieved-metric score and whether the glance tells you whether the visuals agree with the score. Generate for the
visuals agree. Generate for the available pipelines with: available pipelines with:
```bash ```bash
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
@ -290,6 +296,41 @@ No extra model assets: the rigid weights ship inside the image
removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
concurrent benchmark workers can share a GPU. 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 — the release zips for
v2.x need a newer glibc, so this env uses the conda-forge build):
**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.
Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command):
```bash
# one-time, in the benchmark env (provides reg_aladin + reg_f3d, NiftyReg 1.5.69)
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 3` cubic internal warping, default initialisation
(image centres).
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 ## IMPACT pipeline (Docker) — one-time setup
The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact) The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
@ -359,15 +400,17 @@ binary directly (`/lib/elastix-install/bin/elastix`):
head-to-head but is the config the dataset actually ships with. head-to-head but is the config the dataset actually ships with.
- `synthmorph` is deterministic (one fixed pretrained model, no random - `synthmorph` is deterministic (one fixed pretrained model, no random
subsampling), so its per-case numbers are stable run to run — unlike ANTs / subsampling), so its per-case numbers are stable run to run — unlike ANTs /
IMPACT. Its low acceptance (2/10) is a genuine cross-modality / tilt IMPACT. Its low MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
weakness, not variance. genuine cross-modality / tilt weakness, not variance.
- `fix_orientation` (elastix/ants/synthmorph/impact) chooses between a transform - `fix_orientation` (elastix/ants/niftyreg/synthmorph/impact) chooses between a
and its inverse by skull-ring concentricity — a correctness gate that makes transform and its inverse by skull-ring concentricity — a correctness gate
the result independent of each toolkit's transform-direction convention. With that makes the result independent of each toolkit's transform-direction
the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the convention. With the fixed=CT/moving=T1c swap, `elastix`'s combination
inverse of the contract, so it reports `auto_inverted: True` (the gate transform is now the inverse of the contract, so it reports
corrects it); antspyx/synthmorph/impact return the contract direction `auto_inverted: True` (the gate corrects it); antspyx/synthmorph/impact
directly (`auto_inverted: False`). 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): - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU):
- **Cross-modality / tilt weakness.** The `rigid` model is trained on brain - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain
*MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly
@ -382,6 +425,25 @@ binary directly (`/lib/elastix-install/bin/elastix`):
`TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU. `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
- **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK
(~10 s) / ANTs (~37 s) and IMPACT (~834 s). (~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.
- **CPU-only, thread-bound** (`-omp`); ~4 min/case median (120–440 s range)
at full volumes with 24 threads.
- **Version note.** The upstream KCL-BMEIS release zips (v2.x) need
glibc ≥ 2.33; this host (Ubuntu 20.04, glibc 2.31) uses the conda-forge
build (1.5.69) instead.
- **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python - **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python
API — no CLI binary involved): API — no CLI binary involved):
- Elastix 5.x renamed components vs the classic parameter files: the rigid - Elastix 5.x renamed components vs the classic parameter files: the rigid
@ -403,15 +465,19 @@ binary directly (`/lib/elastix-install/bin/elastix`):
(`GetCombinationTransform().TransformPoint`) into an affine — independent (`GetCombinationTransform().TransformPoint`) into an affine — independent
of elastix's Euler-angle parameterization — and wrapped in a SimpleITK of elastix's Euler-angle parameterization — and wrapped in a SimpleITK
`AffineTransform(3)`. `AffineTransform(3)`.
- Registration wall time is ~2.4 s/case median at full budget (ASGD - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
500/500/250 over 3 resolutions). (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, - **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0,
GPU, Jacobian mode): GPU, Jacobian mode):
- **Slow, by design.** Every ASGD iteration runs the TotalSegmentator - **Slow, by design.** Every ASGD iteration runs the TotalSegmentator
extractor on the sampled patches with backpropagation, so the recommended extractor on the sampled patches with backpropagation. The benchmark
500-iter / 2000-sample budget is ~14 min/case on full volumes (vs 2–37 s budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
for the other tools). `--impact-iters` / `--impact-samples` trade accuracy the elastix default) is ~10× the authors' 500-iter / 2000-sample
for speed. Jacobian budget, i.e. on the order of *hours* per case on full volumes —
the results shown for `impact` are from the previous 500/2000 budget and
were not re-run at the new budget. `--impact-iters` / `--impact-samples`
trade accuracy for speed.
- **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT - **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT
follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251 follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251
mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the

109
scripts/bench_registration.py
View file

@ -1,37 +1,41 @@
"""Benchmark the 5 rigid CT->T1c registration pipelines on N M6 test cases. """Benchmark the 6 rigid CT->T1c registration pipelines on N M6 test cases.
Pipelines (see src/registration_pipelines.py for the shared transform contract): Pipelines (see src/registration_pipelines.py for the shared transform contract):
simpleitk | elastix | antspyx | synthmorph | impact simpleitk | elastix | antspyx | synthmorph | impact | niftyreg
Each pipeline registers a head-and-neck CT to a tightly-cropped T1c MRI 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 (cross-modality, inter-session tilt up to ~30 deg) and the GTV is warped onto
the T1c native grid. Every pipeline is scored with the SAME crop-independent the T1c native grid. The headline alignment score is the achieved Mattes MI
geometric QA used in production (scripts/02_build_m6_dataset.py): (mi_ct_mr, nats) between the warped-CT and the original MR over the MR head --
higher = better in-sync = better aligned. The production geometric QA
(scripts/02_build_m6_dataset.py) is still recorded per case in the JSON:
- skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR - mi_ct_mr : Mattes MI (nats), warped-CT ~ MR over the MR head
head silhouette (mid/low level + angular std). (the reported alignment score).
Aligned = 0 <= mid <= 20 mm and std <= 12 mm.
- tv_in_ct_frac : fraction of the (dilated) warped GTV that falls inside - tv_in_ct_frac : fraction of the (dilated) warped GTV that falls inside
the warped-CT head tissue (the slab must straddle it). the warped-CT head tissue (the slab must straddle it).
- tv_nvox : warped GTV volume (>0.5 voxels). - 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).
- accepted : production rule (nvox>0 and tv_in_ct and concentric). - accepted : production rule (nvox>0 and tv_in_ct and concentric).
- reg_wall_s : wall-clock registration time (excludes I/O + scoring). - reg_wall_s : wall-clock registration time (excludes I/O + scoring).
Because some toolkits are not pip-installable into the main Because some toolkits are not pip-installable into every
`longitudinal-313` env (antspyx is kept in an isolated env, see `longitudinal-313` env (antspyx needs a cp310-cp313 wheel, see
REGISTRATION_BENCHMARK.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 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 pipelines that are importable there. elastix runs in the main env through the
ITK-Elastix python API (pip ``itk`` + ``itk-elastix``). ``synthmorph`` and ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively
``impact`` run in Docker images (``freesurfer/synthmorph`` and (``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see ``synthmorph`` and ``impact`` run in Docker images (``freesurfer/synthmorph``
REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the and ``vboussot/elastix_impact`` respectively; needs Docker + the pulled image,
see REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the
per-pipeline JSON files (accumulated across envs) into one summary + plot. per-pipeline JSON files (accumulated across envs) into one summary + plot.
Usage (repo root): Usage (repo root):
# main env (SimpleITK + elastix + antspyx) # main env (SimpleITK + elastix + antspyx + niftyreg)
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg
# synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md # synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
python scripts/bench_registration.py --pipelines synthmorph impact python scripts/bench_registration.py --pipelines synthmorph impact
@ -232,6 +236,8 @@ def _call(pipeline, ct, t1c, tv, args):
if pipeline == "impact": if pipeline == "impact":
return RP.register_impact(ct, t1c, model=args.impact_model, gpu=args.impact_gpu, return RP.register_impact(ct, t1c, model=args.impact_model, gpu=args.impact_gpu,
iters=args.impact_iters, samples=args.impact_samples) iters=args.impact_iters, samples=args.impact_samples)
if pipeline == "niftyreg":
return RP.register_niftyreg(ct, t1c, threads=args.threads)
raise KeyError(pipeline) raise KeyError(pipeline)
@ -255,10 +261,7 @@ def _run_case(pipeline, case, args):
out.update(meta) out.update(meta)
out.update(evaluate(ct, t1c, tv, T)) out.update(evaluate(ct, t1c, tv, T))
if getattr(args, "qa_images", False): if getattr(args, "qa_images", False):
g = out.get("skull_ring_gap_mm") title = f"{pipeline} | {key} | MI={out.get('mi_ct_mr')}"
title = (f"{pipeline} | {key} | MI={out.get('mi_ct_mr')} | "
f"gap mid {g['mid_mm']}mm std {g['std_mm']}mm" if g
else f"{pipeline} | {key}")
qp = os.path.join(args.qa_out, f"{pipeline}__{key}.png") qp = os.path.join(args.qa_out, f"{pipeline}__{key}.png")
try: try:
wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear)) wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
@ -287,12 +290,8 @@ def run_pipeline(pipeline, cases, args):
try: try:
results.append(fu.result()) results.append(fu.result())
r = results[-1] r = results[-1]
g = r.get("skull_ring_gap_mm")
reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT" 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']}: accepted={r['accepted']} " print(f" {c['key']}: mi={r.get('mi_ct_mr')} "
f"gap={g['mid_mm'] if g else 'None'}mm "
f"std={g['std_mm'] if g else 'None'}mm "
f"mi={r.get('mi_ct_mr')} "
f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True) f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True)
except Exception as e: # noqa except Exception as e: # noqa
errors.append({"key": c["key"], "error": repr(e)}) errors.append({"key": c["key"], "error": repr(e)})
@ -348,27 +347,24 @@ def _fmt(v, spec="{}", default="-"):
def print_table(rows): def print_table(rows):
order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)] order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)]
print("\n" + "=" * 100) print("\n" + "=" * 92)
hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} " hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} "
f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}") f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
print(hdr) print(hdr)
print("-" * 100) print("-" * 92)
for p in order: for p in order:
r = next(x for x in rows if x["pipeline"] == p) r = next(x for x in rows if x["pipeline"] == p)
if not r.get("available"): if not r.get("available"):
print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} " print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} "
f"{'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]") f"{'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]")
continue continue
gm = (r.get("skull_ring_gap_mid_mm") or {}).get("median")
gs = (r.get("skull_ring_gap_std_mm") or {}).get("median")
rt = (r.get("reg_wall_s") or {}).get("median") rt = (r.get("reg_wall_s") or {}).get("median")
tv = (r.get("tv_in_ct_frac") or {}).get("mean") tv = (r.get("tv_in_ct_frac") or {}).get("mean")
mi = (r.get("mi_ct_mr") or {}).get("median") mi = (r.get("mi_ct_mr") or {}).get("median")
print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} " 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} {r['accepted_rate']:>7.2f} " f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} "
f"{_fmt(gm, '{:.1f}'):>7s} {_fmt(gs, '{:.1f}'):>7s} {_fmt(rt, '{:.1f}'):>7s} "
f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}") f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
print("=" * 100) print("=" * 92)
def write_plot(rows, out_png): def write_plot(rows, out_png):
@ -384,19 +380,13 @@ def write_plot(rows, out_png):
print("no available pipelines to plot") print("no available pipelines to plot")
return return
labels = list(ok_rows) labels = list(ok_rows)
gap_mid = [ok_rows[p].get("skull_ring_gap_mid_mm", {}).get("median") or 0 for p in labels]
gap_std = [ok_rows[p].get("skull_ring_gap_std_mm", {}).get("median") or 0 for p in labels]
regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels] regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels]
acc = [ok_rows[p].get("accepted_rate") or 0 for p in labels]
mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels] mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels]
fig, axes = plt.subplots(1, 5, figsize=(20, 4)) fig, axes = plt.subplots(1, 2, figsize=(10, 4))
colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD", "#E08E45"] colors = ["#4878CF", "#6ACC65"]
for ax, (title, vals, fmt, ylab, inv) in zip(axes, [ for ax, (title, vals, fmt, ylab) in zip(axes, [
("Skull-ring gap mid (mm, median)", gap_mid, "{:.1f}", "mm", True), ("Registration time (s, median)", regt, "{:.0f}", "s"),
("Skull-ring gap std (mm, median)", gap_std, "{:.1f}", "mm", True), ("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]):
("Registration time (s, median)", regt, "{:.0f}", "s", True),
("Acceptance rate", acc, "{:.2f}", "frac", False),
("Mean Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats", False)]):
bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)]) bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)])
for b, v in zip(bars, vals): for b, v in zip(bars, vals):
ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()), ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()),
@ -426,10 +416,11 @@ def merge(args):
print_table(rows) print_table(rows)
summary = { summary = {
"timestamp": datetime.now().isoformat(timespec="seconds"), "timestamp": datetime.now().isoformat(timespec="seconds"),
"protocol": ("rigid CT->T1c; cross-modality; scored by production skull-ring " "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI "
"concentricity (aligned: 0<=mid<=20mm, std<=12mm) + GTV-in-CT fraction " "(nats, warped-CT ~ MR over the MR head, higher = better); reg time "
"(>=0.5); reg time excludes I/O+scoring; per-pipeline best-of-T vs T^-1 " "excludes I/O+scoring; per-case JSON retains the production skull-ring "
"orientation gate for elastix/antspyx"), "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), "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", [])}), "case_keys": sorted({c["key"] for r in rows if r.get("available") for c in r.get("cases", [])}),
"results": rows, "results": rows,
@ -440,18 +431,16 @@ def merge(args):
print("summary:", sf) print("summary:", sf)
write_plot(rows, os.path.join(out_dir, "bench_registration.png")) write_plot(rows, os.path.join(out_dir, "bench_registration.png"))
# markdown table for the report # markdown table for the report
md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT | MI (nats) |", 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]: for p in [x["pipeline"] for x in rows]:
r = next(x for x in rows if x["pipeline"] == p) r = next(x for x in rows if x["pipeline"] == p)
if not r.get("available"): if not r.get("available"):
md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | - | - | - | ({r['reason']}) |") md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | ({r['reason']}) |")
continue continue
md.append( md.append(
f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | " f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | "
f"{r['n_error']} | {r['accepted_rate']:.2f} | " f"{r['n_error']} | "
f"{(r.get('skull_ring_gap_mid_mm') or {}).get('median')} | "
f"{(r.get('skull_ring_gap_std_mm') or {}).get('median')} | "
f"{(r.get('reg_wall_s') or {}).get('median')} | " f"{(r.get('reg_wall_s') or {}).get('median')} | "
f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | " f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | "
f"{(r.get('mi_ct_mr') or {}).get('median')} |") f"{(r.get('mi_ct_mr') or {}).get('median')} |")
@ -492,12 +481,12 @@ def parse_args():
ap.add_argument("--impact-gpu", type=int, default=None, ap.add_argument("--impact-gpu", type=int, default=None,
help="GPU device index for IMPACT (default: auto-spread over " help="GPU device index for IMPACT (default: auto-spread over "
"the host GPUs by pid; -1 forces CPU)") "the host GPUs by pid; -1 forces CPU)")
ap.add_argument("--impact-iters", type=int, default=500, ap.add_argument("--impact-iters", type=int, default=2000,
help="IMPACT per-resolution max iterations (default 500, the authors' " help="IMPACT MaximumNumberOfIterations (default 2000, the benchmark "
"recommended Jacobian budget; lower to trade accuracy for speed)") "budget; ~10x the authors' 500-iter Jacobian budget, much slower)")
ap.add_argument("--impact-samples", type=int, default=2000, ap.add_argument("--impact-samples", type=int, default=5000,
help="IMPACT random spatial samples per iteration (default 2000, " help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
"the recommended budget)") "elastix default)")
ap.add_argument("--qa-images", action="store_true", ap.add_argument("--qa-images", action="store_true",
help="also write a checkreg-style PNG (checkerboard blend + color " help="also write a checkreg-style PNG (checkerboard blend + color "
"overlay of the inverse-warped CT vs the original MR) per case") "overlay of the inverse-warped CT vs the original MR) per case")

169
src/registration_pipelines.py
View file

@ -1,4 +1,4 @@
"""Five rigid CT -> T1c registration pipelines behind one unified contract. """Six rigid CT -> T1c registration pipelines behind one unified contract.
Problem (shared with ``scripts/02_build_m6_dataset.py``) Problem (shared with ``scripts/02_build_m6_dataset.py``)
-------------------------------------------------------- --------------------------------------------------------
@ -27,10 +27,11 @@ All pipelines register on the **full volumes** (no head box) and use each
toolkit's **native automatic transform initialization** (no manual toolkit's **native automatic transform initialization** (no manual
``centroid_offset``): elastix / impact use elastix's ``centroid_offset``): elastix / impact use elastix's
``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its ``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its
default center-of-mass alignment, simpleitk starts from identity, and default center-of-mass alignment, simpleitk starts from identity, niftyreg
synthmorph needs none (global deep method). All use **fixed = CT, moving = (reg_aladin) initializes from the image centres, and synthmorph needs none
T1c** except IMPACT, which keeps **fixed = T1c, moving = CT** (swapping breaks (global deep method). All use **fixed = CT, moving = T1c** except IMPACT,
its TotalSegmentator metric -- see ``register_impact``). which keeps **fixed = T1c, moving = CT** (swapping breaks its
TotalSegmentator metric -- see ``register_impact``).
- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). - ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
Default ``initial_transform=None`` computes a translation aligning the image Default ``initial_transform=None`` computes a translation aligning the image
@ -44,6 +45,11 @@ its TotalSegmentator metric -- see ``register_impact``).
GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix
``AutomaticTransformInitialization`` (CenterOfGravity) provides the start. ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
Needs Docker + the image + a model asset. Needs Docker + the image + a model asset.
- ``register_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
symmetric block-matching + trimmed least squares (Modat et al., J. Med.
Img. 2014), NCC-based block similarity (cross-modality tolerant). Rigid
6-DOF (``-rigOnly``), 3-level pyramid, 15 TLS iterations/level. Needs
``conda install -c conda-forge niftyreg`` (provides ``reg_aladin``).
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on the - ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on the
full volumes. Reuses the *proven* production ``register_ct_to_t1c`` in full volumes. Reuses the *proven* production ``register_ct_to_t1c`` in
``mode="production"``; the benchmark mode is a single identity-start ``mode="production"``; the benchmark mode is a single identity-start
@ -71,6 +77,7 @@ applied to elastix / antspyx / synthmorph / impact, not to the trusted
production SimpleITK path. production SimpleITK path.
""" """
import os import os
import re
import shutil import shutil
import subprocess import subprocess
import importlib.util import importlib.util
@ -162,13 +169,13 @@ def centroid_offset(ct, t1c, ct_box=None):
# 2. SimpleITK (production baseline -- reused verbatim) # 2. SimpleITK (production baseline -- reused verbatim)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150): def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
"""SimpleITK rigid CT->T1c, Mattes MI, full volumes (no head box). """SimpleITK rigid CT->T1c, Mattes MI, full volumes (no head box).
mode="benchmark" (default): a single multi-resolution pyramid run from the mode="benchmark" (default): a single multi-resolution pyramid run from the
identity start, capped at ``iters_full`` iterations at full resolution. identity start, with ``iters_full`` (=2000) optimizer iterations at every
This is the *standard* SimpleITK algorithm and matches the budget the pyramid level. This is the *standard* SimpleITK algorithm and matches the
other toolkits get in the benchmark (one rigid MI registration), so the 2000-iteration budget the other toolkits get in the benchmark, so the
pipelines are compared on equal footing. No centroid init is used. pipelines are compared on equal footing. No centroid init is used.
mode="production": the exact multi-start tilt-grid ``register_ct_to_t1c`` mode="production": the exact multi-start tilt-grid ``register_ct_to_t1c``
from scripts/02_build_m6_dataset.py (robust to the ~30 deg inter-session from scripts/02_build_m6_dataset.py (robust to the ~30 deg inter-session
@ -297,6 +304,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
CT / T1c volumes (no head box). No manual centroid init: elastix's CT / T1c volumes (no head box). No manual centroid init: elastix's
``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity
centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales. centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales.
Budget: ``MaximumNumberOfIterations`` 2000 per resolution; the sampler /
metric run at their defaults (``NumberOfSpatialSamples`` 5000,
``NumberOfHistogramBins`` 32).
Elastix 5.x Elastix 5.x
component names: sampler ``Random`` (not ``RandomSampler``), metric key component names: sampler ``Random`` (not ``RandomSampler``), metric key
``Metric`` (not ``ImageMetric``). An output directory is required for the ``Metric`` (not ``ImageMetric``). An output directory is required for the
@ -321,9 +331,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
"AutomaticTransformInitializationMethod": "CenterOfGravity", "AutomaticTransformInitializationMethod": "CenterOfGravity",
"AutomaticScalesEstimation": True, "AutomaticScalesEstimation": True,
"Optimizer": "AdaptiveStochasticGradientDescent", "Optimizer": "AdaptiveStochasticGradientDescent",
"MaximumNumberOfIterations": [500, 500, 250], "MaximumNumberOfIterations": [2000, 2000, 2000],
"Metric": "AdvancedMattesMutualInformation", "Metric": "AdvancedMattesMutualInformation",
"NumberOfHistogramBins": 10, "NumberOfHistogramBins": 32,
"ImageSampler": "Random", "ImageSampler": "Random",
"ImageSamplerPatchSize": [25, 25, 25], "ImageSamplerPatchSize": [25, 25, 25],
"WriteResultImage": False, "WriteResultImage": False,
@ -735,7 +745,7 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
return "\n".join(L) + "\n" return "\n".join(L) + "\n"
def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
threads=12, keep_scratch=False, scales=(6.0, 3.0, 1.5, 1.0), threads=12, keep_scratch=False, scales=(6.0, 3.0, 1.5, 1.0),
downsample=False): downsample=False):
"""IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image
@ -761,11 +771,11 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
outside the T1c buffer and the IMPACT metric aborts ("too many samples map outside the T1c buffer and the IMPACT metric aborts ("too many samples map
outside moving image buffer"). The returned transform still satisfies the outside moving image buffer"). The returned transform still satisfies the
unified contract (T1c->CT). unified contract (T1c->CT).
Budget: the default ``iters=500`` / ``samples=2000`` is the authors' Budget: the default ``iters=2000`` (MaximumNumberOfIterations) /
recommended Jacobian budget. It is *much* slower than the intensity-metric ``samples=5000`` (NumberOfSpatialSamples, the elastix default) runs the
pipelines (minutes per case on GPU vs seconds) because every iteration runs full benchmark budget for best results -- roughly 10x the authors'
the TotalSegmentator feature extractor on the sampled patches with 500/2000 Jacobian budget, i.e. hours per case on full volumes. Lower
backpropagation. Lower ``iters``/``samples`` to trade accuracy for speed. ``iters``/``samples`` to trade accuracy for speed.
The result is read back from ``TransformParameters.0.itk.txt`` The result is read back from ``TransformParameters.0.itk.txt``
(an ``Euler3DTransform`` in the world physical frame) and checked by (an ``Euler3DTransform`` in the world physical frame) and checked by
:func:`fix_orientation`. :func:`fix_orientation`.
@ -886,10 +896,128 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
return T, meta return T, meta
# --------------------------------------------------------------------------- #
# 6. NiftyReg (KCL-BMEIS suite, reg_aladin)
# --------------------------------------------------------------------------- #
def _niftyreg_bin():
"""Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable."""
b = shutil.which("reg_aladin")
if b:
return b
raise PipelineUnavailable(
"niftyreg: reg_aladin not found on PATH. Install NiftyReg with "
"`conda install -c conda-forge niftyreg` (NiftyReg 1.5.69, the "
"KCL-BMEIS suite: reg_aladin rigid/affine + reg_f3d non-linear) and "
"re-run in the activated env.")
def _niftyreg_version(binpath):
try:
r = subprocess.run([binpath, "--version"], capture_output=True,
text=True, timeout=60)
m = re.search(r"(\d+\.\d+\.\d+)", r.stdout + r.stderr)
if m:
return f"NiftyReg {m.group(1)} (reg_aladin)"
except Exception: # noqa
pass
return "NiftyReg (reg_aladin)"
def _parse_niftyreg_affine(path):
"""Parse reg_aladin's output affine (4x4, row-major) into a SimpleITK
AffineTransform. The file is *documented* as
``Affine * Reference = Floating``; the direction that actually satisfies
the contract is verified empirically by :func:`fix_orientation` (on the
M6 cases the gate keeps the file's affine as-is)."""
M = np.array([float(x) for x in open(path).read().split()]).reshape(4, 4)
T = sitk.AffineTransform(3)
T.SetMatrix(M[:3, :3].ravel())
T.SetTranslation(M[:3, 3])
return T
def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False):
"""NiftyReg rigid CT->T1c (reference=CT, floating=T1c) via ``reg_aladin``.
reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global
block-matching registration: at each level of a ``levels``-level Gaussian
pyramid it matches 4-voxel blocks by normalised cross-correlation
(NCC -- a cross-modality-tolerant block similarity, unlike the histogram
MI objectives of the other pipelines), then solves a trimmed least squares
(TLS) problem for the rigid transform, iterating ``maxit`` times per
level. The symmetric variant (default) optimises forward and backward
transforms concurrently; ``-maxit 5`` (the tool default) is too coarse for
these repositioned pairs -- the benchmark uses 15.
Full volumes (no head box); default initialisation (image centres);
``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3`` (cubic) internal
warping. The output affine is documented as ``Affine * Reference =
Floating``; :func:`fix_orientation` resolves the direction that actually
satisfies the contract against the skull ring (on the M6 cases it keeps
the file's affine as-is, ``auto_inverted: False``) and returns the T1c ->
CT transform.
Preconditions (else :class:`PipelineUnavailable`): ``reg_aladin`` on PATH
(``conda install -c conda-forge niftyreg``).
"""
import time
binpath = _niftyreg_bin()
# full volumes, reference = CT, floating = T1c (no head box)
tmp = tempfile.mkdtemp(prefix="niftyreg_reg_")
dt = 0.0
try:
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
os.path.join(tmp, "ref.nii.gz"), True)
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
os.path.join(tmp, "flo.nii.gz"), True)
cmd = [binpath, "-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
"-rigOnly", "-aff", "aff.txt",
"-ln", str(levels), "-maxit", str(maxit),
"-interp", "3", "-omp", str(threads), "-voff"]
t0 = time.time()
proc = subprocess.run(cmd, capture_output=True, text=True,
cwd=tmp, timeout=7200)
dt = time.time() - t0
if proc.returncode != 0:
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
raise RuntimeError(f"niftyreg: reg_aladin failed (rc={proc.returncode}):\n{tail}")
aff = os.path.join(tmp, "aff.txt")
if not os.path.exists(aff):
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
raise RuntimeError(f"niftyreg: no affine output produced:\n{tail}")
T = _parse_niftyreg_affine(aff)
# aff.txt is documented as CT -> T1c (Affine*Reference=Floating);
# the contract is T1c -> CT. fix_orientation resolves the direction
# empirically (on these cases it keeps the file's affine as-is).
T, inverted = fix_orientation(ct, t1c, T)
finally:
if not keep_scratch:
shutil.rmtree(tmp, ignore_errors=True)
meta = {
"pipeline": "niftyreg",
"toolkit": _niftyreg_version(binpath),
"metric": "NCC block matching + trimmed least squares (symmetric, full volumes)",
"transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, "
f"{maxit} TLS iters/level, cubic interp)"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "CT", "moving": "T1c",
"seconds_reg": round(dt, 1),
"levels": levels, "maxit": maxit,
}
if keep_scratch:
meta["scratch"] = tmp
return T, meta
REGISTERS = { REGISTERS = {
"antspyx": register_antspyx, "antspyx": register_antspyx,
"elastix": register_elastix, "elastix": register_elastix,
"impact": register_impact, "impact": register_impact,
"niftyreg": register_niftyreg,
"simpleitk": register_simpleitk, "simpleitk": register_simpleitk,
"synthmorph": register_synthmorph, "synthmorph": register_synthmorph,
} }
@ -917,6 +1045,13 @@ def availability():
out[name] = (True, f"ants importable v{getattr(ants, '__version__', '?')}") out[name] = (True, f"ants importable v{getattr(ants, '__version__', '?')}")
except Exception as e: # noqa except Exception as e: # noqa
out[name] = (False, f"ants(antspyx) not importable ({type(e).__name__})") out[name] = (False, f"ants(antspyx) not importable ({type(e).__name__})")
elif name == "niftyreg":
b = shutil.which("reg_aladin")
if b:
out[name] = (True, _niftyreg_version(b))
else:
out[name] = (False, "reg_aladin not on PATH "
"(conda install -c conda-forge niftyreg)")
elif name == "synthmorph": elif name == "synthmorph":
ok, note = _synthmorph_docker_ok() ok, note = _synthmorph_docker_ok()
if not ok: if not ok: