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

Enhance registration pipelines with CT clipping and NiftyReg Docker support

Browse source 
- Added `--clip-ct LO HI` argument to clip CT intensities before registration and scoring, improving A/B comparisons.
- Implemented `_clip_image` function to handle intensity clipping of images.
- Updated registration functions to utilize CT clipping where beneficial, specifically in `register_niftyreg` and `register_synthmorph`.
- Enhanced `fixed_moving` function to ensure CT is always treated as the fixed image.
- Introduced Docker support for NiftyReg v2.1.1, including a Dockerfile for building the image.
- Updated documentation and help messages to reflect new features and usage instructions.
- Adjusted evaluation criteria for registration quality to account for intensity clipping.
This commit is contained in:
Furen Xiao 2026-09-29 17:31:24 +08:00
parent a65189cd33
commit 16b9828d15
6 changed files with 587 additions and 238 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
.gitignore vendored
View file

@ -8,4 +8,6 @@ runs_nnu/
runs_monai/ runs_monai/
__pycache__/ __pycache__/
elastix/ elastix/
niftyreg/NiftyReg-Ubuntu-v2.1.1.zip
niftyreg/NiftyReg-Ubuntu-v2.1.1/

349
REGISTRATION_BENCHMARK.md
View file

@ -22,36 +22,47 @@ 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 places the CT / GTV concentric with the MR head. In `Resample`'s
destination→source convention this `T` maps **T1c → CT**. All pipelines destination→source convention this `T` maps **T1c → CT**. Shared rules for
register on the **full CT / T1c volumes** (no head box) and use each toolkit's every pipeline:
**native automatic transform initialization** (no manual centroid offset):
elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity), - **No mask, no crop / pad**: registration runs on the full input volumes
antspyx its default center-of-mass alignment, simpleitk starts from identity, exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
niftyreg (reg_aladin) initializes from the image centres, synthmorph needs padding).
none. All use **fixed = CT, moving = T1c** except `impact`, which keeps - **CT is always the fixed image** (moving = T1c); when both inputs are MR
**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than the first argument is kept as fixed.
the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
samples map outside the T1c buffer and the IMPACT metric aborts ("too many All pipelines also use each toolkit's **native automatic transform
samples map outside moving image buffer"). The comparison is on each initialization** (no manual centroid offset): elastix / impact use
toolkit's native optimizer / metric / global-search capability. `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 and `impact` aborts on the M6
pairs ("too many samples map outside moving image buffer", ~14% of samples
in-buffer at the auto-init) — see the IMPACT section. The comparison is on
each toolkit's native optimizer / metric / global-search capability.
## The six 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, no
use each toolkit's **native automatic transform initialization** (no manual mask, no crop / pad) and use each toolkit's **native automatic transform
centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization` initialization** (no manual centroid offset): elastix / impact use elastix's
(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk `AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default
starts from identity, niftyreg initializes from the image centres, and center-of-mass alignment, simpleitk starts from identity, niftyreg
synthmorph needs none (global deep method). All use **fixed = CT, moving = initializes from the image centres, and synthmorph needs none (global deep
T1c** except `impact`, which keeps **fixed = T1c, moving = CT** (swapping method). **All use fixed = CT, moving = T1c** (shared rule above; when both
breaks its TotalSegmentator metric — see below). 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 | | # | 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 (32 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` | **CT / T1c** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ⛔ aborts on M6 (no-mask rule; see IMPACT section) |
| 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) | | 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 | | 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) | | 6 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
@ -62,10 +73,10 @@ level), with the sampler and metric left at their defaults
`niftyreg` and `synthmorph` keep their native defaults. `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 identity start — the default, comparable to the other tools) and
`--simpleitk-mode production` (the full multi-start tilt grid from `--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 `02_build_m6_dataset.py`, now also on the full volumes per the no-crop rule,
head-to-head). ~20× slower, the most robust but not used for the head-to-head).
## Test cases ## Test cases
@ -108,40 +119,45 @@ longer reported in the tables/plots):
translational concentricity; `std` measures rotational concentricity. translational concentricity; `std` measures rotational concentricity.
- **`tv_nvox`** — warped GTV volume. - **`tv_nvox`** — warped GTV volume.
- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the - **`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`). 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) ## Results (10 cases, no clip)
Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
the MR head; higher = better). 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) | | pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|---|---|---|---|---|---| |---|---|---|---|---|---|
| **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** | | **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.516 / 0.520** | | **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.517 / 0.520** |
| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** | | **impact** | Docker Elastix 5.2 + IMPACT | 0 | - | - | - (aborts, no mask) |
| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** | | **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 20.8 | 0.93 | **0.292 / 0.338** |
| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** | | **simpleitk** | SimpleITK 2.5.6 | 9 | 10.6 | 0.93 | **0.651 / 0.641** |
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** | | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 24.5 | 0.93 | **0.277 / 0.328** |
\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default Per-case MI (nats, unclipped score; **bold** = best for the case):
(MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
(hours/case) and was not re-run.
Per-case MI (nats; **bold** = best for the case):
| case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph | | case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
|---|---|---|---|---|---|---| |---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | 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_22IRRKU6_20240115 | **0.649** | 0.472 | abort | 0.292 | 0.642 | 0.277 |
| m6_23XIMO4M_20230222 | 0.655 | 0.416 | 0.474 | **0.634** | 0.652 | 0.510 | | m6_23XIMO4M_20230222 | **0.655** | 0.416 | abort | 0.632 | 0.652 | 0.616 |
| m6_24V62QQM_20230223 | 0.703 | **0.701** | 0.554 | 0.260 | 0.663 | 0.238 | | m6_24V62QQM_20230223 | **0.704** | 0.701 | abort | 0.260 | 0.663 | 0.256 |
| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.419 | 0.493 | **0.602** | 0.516 | | m6_24V62QQM_20231017 | 0.601 | 0.590 | abort | 0.493 | **0.602** | 0.484 |
| m6_24V62QQM_20240820 | 0.652 | 0.643 | 0.458 | 0.118 | **0.652** | 0.096 | | m6_24V62QQM_20240820 | **0.653** | 0.644 | abort | 0.118 | 0.652 | 0.111 |
| m6_26TYQI3C_20221027 | 0.657 | 0.653 | 0.511 | 0.608 | **0.695** | 0.604 | | m6_26TYQI3C_20221027 | 0.659 | 0.653 | abort | 0.609 | **0.695** | 0.590 |
| m6_2BBN5XA4_20250411 | 0.576 | 0.287 | 0.155 | 0.069 | **0.577** | 0.060 | | m6_2BBN5XA4_20250411 | **0.577** | 0.287 | abort | 0.069 | **0.577** | 0.061 |
| m6_2C4TKI6H_20240506 | 0.645 | 0.401 | **0.626** | 0.231 | 0.634 | 0.209 | | m6_2C4TKI6H_20240506 | **0.645** | 0.401 | abort | 0.233 | 0.634 | 0.231 |
| m6_2CDNN7BH_20240102 | 0.655 | 0.516 | **0.618** | 0.344 | 0.651 | 0.301 | | m6_2CDNN7BH_20240102 | **0.655** | 0.516 | abort | 0.340 | 0.651 | 0.323 |
### Interpretation ### Interpretation
@ -163,52 +179,128 @@ Per-case MI (nats; **bold** = best for the case):
single CenterOfGravity start still locks onto local minima on the single CenterOfGravity start still locks onto local minima on the
repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M` repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
0.416). 0.416).
- **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI - **IMPACT (semantic metric) does not run on these pairs under the no-mask
on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618), rule** — earlier (masked) runs showed it strongest on the hardest cases
where anatomical-segmentation similarity out-resolves intensity MI — but its (`m6_22IRRKU6` 0.624, `m6_2C4TKI6H` 0.626, `m6_2CDNN7BH` 0.618, median
median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is 0.511) with a stochastic draw per run; it also was by far the slowest
also by far the slowest (~14 min/case at the shown budget; the new 2000/5000 (~14 min/case at the 500/2000 budget, hours/case at 2000/5000). With the CT
budget would be on the order of hours/case). head-neck FOV as the fixed volume it aborts ("too many samples map outside
moving image buffer", ~10–20% in-buffer) — see the IMPACT section.
- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching - **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class; is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class;
`m6_26TYQI3C` 0.608) but falls into the same tilt local minima `m6_26TYQI3C` 0.609) but falls into the same tilt local minima
(`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower (`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). The v2.1.1 Docker build
than elastix at full volumes. 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** - **SynthMorph is weakest of the six on these cross-modality CT pairs**
(median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30° (median 0.277 with its default CT clip; 0.241 unclipped): its `rigid` model
tilt defeats it on the repositioned cases (0.06–0.10) even though it matches is trained on brain *MRI*, so the ~30° tilt defeats it on the repositioned
the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no cases (0.06–0.11) even though it matches the intensity tools on the
centroid init** (global deep method) and is deterministic. 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 - **`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. CT head everywhere, so coverage is not the discriminator — MI is.
- **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s, - **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s,
SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s, SimpleITK 10.6 s, **NiftyReg 20.8 s** (v2.1.1, GPU block matching),
**NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT does not run on
budget; GPU feature extraction with backprop per iteration). M6 under the no-mask rule; when it did — masked, 500/2000 budget — it was
~833 s/case, GPU feature extraction with backprop per iteration.)
- **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;
the in-process Elastix ASGD is deterministic given its fixed seed; the in-process Elastix ASGD is deterministic given its fixed seed;
SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**: 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 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 changes run to run. Treat its per-case numbers as one sample from a
distribution. 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 | - | - (still aborts; clipping does not fix the FOV geometry, 194–388/2000 samples in-buffer) |
**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 ## Reproducing
```bash ```bash
# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (pip/conda-installed here) # main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1)
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 niftyreg \ 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
# force the conda 1.5.69 CPU build: --niftyreg-backend conda
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT # Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
# sections below for the one-time setup) # sections below for the one-time setup)
python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \ python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \
--workers 4 --threads 24 --workers 4 --threads 24
# NOTE: impact aborts on the M6 pairs under the shared no-mask rule (the CT
# head-neck FOV maps mostly outside the T1c head crop) — see the IMPACT section
python scripts/bench_registration.py --pipelines impact --n-cases 10 \ python scripts/bench_registration.py --pipelines impact --n-cases 10 \
--workers 4 --threads 24 --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 # combine the per-pipeline JSONs accumulated across envs into one table + plot
python scripts/bench_registration.py --merge python scripts/bench_registration.py --merge
``` ```
@ -249,14 +341,13 @@ 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
python scripts/bench_registration.py \ python scripts/bench_registration.py \
--pipelines simpleitk elastix antspyx --n-cases 10 --qa-images --pipelines simpleitk elastix antspyx niftyreg --n-cases 10 --qa-images
python scripts/bench_registration.py \ python scripts/bench_registration.py \
--pipelines synthmorph --n-cases 10 --qa-images --pipelines synthmorph --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.) (9 images each; the incomplete-CT case is skipped. `impact` produces none on
M6 — it aborts under the no-mask rule, see the IMPACT section.)
## SynthMorph pipeline (Docker) — one-time setup ## SynthMorph pipeline (Docker) — one-time setup
@ -299,20 +390,33 @@ No extra model assets: the rigid weights ship inside the image
## NiftyReg pipeline — one-time setup ## NiftyReg pipeline — one-time setup
The `niftyreg` pipeline runs `reg_aladin` from the [NiftyReg](https://github.com/KCL-BMEIS/niftyreg) 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 suite (KCL-BMEIS, the actively-maintained NiftyReg): **symmetric block
v2.x need a newer glibc, so this env uses the conda-forge build): matching + trimmed least squares** (Modat et al., J. Med. Img. 1(2) 024003,
**symmetric block matching + trimmed least squares** (Modat et al., 2014). Instead of a histogram MI objective it matches 4-voxel image blocks by
J. Med. Img. 1(2) 024003, 2014). Instead of a histogram MI objective it **normalised cross-correlation** at each level of a 3-level Gaussian pyramid
matches 4-voxel image blocks by **normalised cross-correlation** at each level and solves a robust (trimmed) least-squares for the rigid transform — a
of a 3-level Gaussian pyramid and solves a robust (trimmed) least-squares for genuinely different algorithmic family, and cross-modality tolerant because
the rigid transform — a genuinely different algorithmic family, and NCC is invariant to intensity scaling.
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 Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command): `PipelineUnavailable` with the exact command):
```bash ```bash
# one-time, in the benchmark env (provides reg_aladin + reg_f3d, NiftyReg 1.5.69) # 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 conda install -c conda-forge niftyreg
``` ```
@ -322,8 +426,10 @@ conda install -c conda-forge niftyreg
(symmetric variant on by default, which matters: the non-symmetric mode (symmetric variant on by default, which matters: the non-symmetric mode
drifts grossly on these repositioned pairs). `-ln 3` pyramid, drifts grossly on these repositioned pairs). `-ln 3` pyramid,
`-maxit 15` TLS iterations per level (the tool's default of 5 is too `-maxit 15` TLS iterations per level (the tool's default of 5 is too
coarse), `-interp 3` cubic internal warping, default initialisation coarse), `-interp 1` (GPU) / `-interp 3` (CPU) internal warping, default
(image centres). 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, 2. **Read-back** — the output affine (`-aff`) is a 4×4 row-major matrix,
documented as `Affine * Reference = Floating` (CT→T1c); the adapter parses documented as `Affine * Reference = Floating` (CT→T1c); the adapter parses
it into a SimpleITK `AffineTransform` and runs it through the usual it into a SimpleITK `AffineTransform` and runs it through the usual
@ -337,11 +443,13 @@ The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.co
image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121) 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 *semantic* similarity metric compiled in as a native Elastix metric, plus a
TorchScript **TotalSegmentator** feature extractor. It is GPU-based and TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps Docker-based (alongside `synthmorph`). It follows the shared rules (full
**fixed = T1c, moving = CT** (the reverse of the others): the CT FOV volumes, no mask, fixed = CT, moving = T1c) — which means it **aborts on the
(~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with M6 pairs**: the IMPACT metric only evaluates samples that map inside the
CT as the fixed image too many CT samples map outside the T1c buffer and the moving buffer, and the CT head-neck FOV (~300×300×251 mm) is much larger than
IMPACT metric aborts ("too many samples map outside moving image buffer"). the T1c head crop (~220×220×150 mm), so ~14% of the samples are in-buffer at
the auto-init and elastix raises "too many samples map outside moving image
buffer" (see the IMPACT specifics below).
Prerequisites (check with `RP.availability()`; a missing piece raises Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command): `PipelineUnavailable` with the exact command):
@ -370,9 +478,15 @@ repo for the full list and FOV/layer tables.
The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
binary directly (`/lib/elastix-install/bin/elastix`): binary directly (`/lib/elastix-install/bin/elastix`):
1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box). 1. **Full volumes, fixed=CT, moving=T1c** — both volumes are used as-is (no
T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too head box, no mask, no crop / pad, per the shared rules above). **On the M6
many CT samples map outside the T1c buffer and the metric aborts. pairs this aborts**: the IMPACT 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 2000 samples are in-buffer and
elastix raises "Too many samples map outside moving image buffer: 277 /
2000". IMPACT therefore needs a fixed volume whose FOV is comparable to
the moving crop (e.g. a head-only CT) to work under the no-mask rule.
2. **Auto init (no `-t0`)** — the parameter map sets 2. **Auto init (no `-t0`)** — the parameter map sets
`AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
@ -391,13 +505,14 @@ binary directly (`/lib/elastix-install/bin/elastix`):
## Notes / caveats ## Notes / caveats
- The benchmark compares each toolkit in its **native configuration**: all - The benchmark compares each toolkit in its **native configuration**: all
register on the full volumes (no head box) and use each toolkit's native register on the full volumes (no head box, no mask, no crop / pad) and use
automatic transform initialization (no manual `centroid_offset`) — elastix / each toolkit's native automatic transform initialization (no manual
impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's `centroid_offset`) — elastix / impact `AutomaticTransformInitialization`
center-of-mass alignment, simpleitk identity, synthmorph none. Production (CenterOfGravity), antspyx's center-of-mass alignment, simpleitk identity,
uses `--simpleitk-mode production`'s multi-start tilt grid, which raises synthmorph none. Production uses `--simpleitk-mode production`'s
SimpleITK's acceptance at ~20× the cost; that is out of scope for the multi-start tilt grid, which raises SimpleITK's acceptance at ~20× the
head-to-head but is the config the dataset actually ships with. 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 - `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 MI (median 0.238; 0.06–0.10 on the repositioned cases) is a IMPACT. Its low MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
@ -439,11 +554,16 @@ binary directly (`/lib/elastix-install/bin/elastix`):
- **`-maxit 5` (the tool default) is too coarse** — on the clean case - **`-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; `m6_26TYQI3C` it reaches MI 0.41 vs 0.61 at 15 TLS iterations/level;
15 is the benchmark setting. 15 is the benchmark setting.
- **CPU-only, thread-bound** (`-omp`); ~4 min/case median (120–440 s range) - **GPU is the default now** (v2.1.1, `-platf 1` via
at full volumes with 24 threads. `docker run --gpus all`): ~20–25 s/case on the full volumes vs ~4 min
- **Version note.** The upstream KCL-BMEIS release zips (v2.x) need CPU (120–440 s, `-omp` thread-bound) for the conda 1.5.69 build.
glibc ≥ 2.33; this host (Ubuntu 20.04, glibc 2.31) uses the conda-forge `--niftyreg-gpu -1` forces the CPU platform (`-platf 0`).
build (1.5.69) instead. - **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 - **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
@ -478,20 +598,21 @@ binary directly (`/lib/elastix-install/bin/elastix`):
the results shown for `impact` are from the previous 500/2000 budget and 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` were not re-run at the new budget. `--impact-iters` / `--impact-samples`
trade accuracy for speed. trade accuracy for speed.
- **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT - **Aborts on M6 under the no-mask rule.** Like every other pipeline it
follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251 runs fixed=CT/moving=T1c on the full volumes (no mask, no crop / pad).
mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the The CT FOV (~300×300×251 mm) is much larger than the T1c crop
fixed image too many CT samples map outside the T1c buffer and the metric (~220×220×150 mm), so most of the CT buffer maps outside the T1c moving
aborts ("too many samples map outside moving image buffer"). The returned buffer and the metric aborts ("Too many samples map outside moving image
transform still satisfies the contract (T1c→CT). buffer: 277 / 2000" on M6 at the CenterOfGravity auto-init). The
transform contract (T1c→CT) still holds when a run completes; for M6
pairs a head-only CT as the fixed volume would be needed.
- **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns
the intensity centroids; no manual `-t0` translation is used. the intensity centroids; no manual `-t0` translation is used.
- **Rugged single-init convergence.** The semantic landscape has more local - **History: single-init convergence.** The semantic landscape has more
minima than Mattes-MI: on full volumes it no longer nails the hardest local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer
repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0;
used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration the pre-mask head-box flow got it to mid 8.3 / std 13.2). Under the
budget would likely raise its acceptance; that is out of scope for the current no-mask rule it does not run on M6 at all (see above).
single-init head-to-head.
- **Stochastic.** Jacobian mode resamples patches each iteration, so results - **Stochastic.** Jacobian mode resamples patches each iteration, so results
vary run to run (see Stochasticity). vary run to run (see Stochasticity).
- **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image

35
niftyreg/Dockerfile Normal file
View file

@ -0,0 +1,35 @@
# NiftyReg v2.1.1 (official KCL-BMEIS Ubuntu release) in Docker.
#
# The host (Ubuntu 20.04, glibc 2.31) is too old for the v2.x binaries, which
# need GLIBC_2.34; the nvidia/cuda ubuntu22.04 base (glibc 2.35) runs them
# directly. The CUDA runtime provides libcusolver.so.11, which the
# libniftyreg_cuda.so plugin (dlopen'd by reg_aladin with -platf 1) links
# against; libcuda.so.1 comes from the host driver via `docker run --gpus all`.
# Build:
# docker build -t niftyreg:2.1.1-ubuntu niftyreg/
# Run (GPU, the pipeline's default):
# docker run --rm --gpus all -v /path/to/work:/work -w /work \
# niftyreg:2.1.1-ubuntu -ref ref.nii.gz -flo flo.nii.gz \
# -rigOnly -aff aff.txt -platf 1
# Run (CPU fallback):
# docker run --rm -v /path/to/work:/work -w /work \
# niftyreg:2.1.1-ubuntu -ref ref.nii.gz -flo flo.nii.gz \
# -rigOnly -aff aff.txt -platf 0
# reg_aladin is the ENTRYPOINT, so the NiftyReg CLI flags follow directly.
# -platf 0|1 = CPU|CUDA, -gpuid selects the device, -omp sets OpenMP threads.
FROM nvidia/cuda:12.4.1-runtime-ubuntu22.04
RUN apt-get update && apt-get install -y --no-install-recommends \
libgomp1 libstdc++6 zlib1g \
&& rm -rf /var/lib/apt/lists/*
COPY NiftyReg-Ubuntu-v2.1.1/ /niftyreg/
# reg_aladin searches /usr/local/lib for the CUDA plugin
RUN chmod +x /niftyreg/* \
&& cp /niftyreg/libniftyreg_cuda.so /usr/local/lib/ \
&& ldconfig
ENV PATH=/niftyreg:$PATH \
LD_LIBRARY_PATH=/niftyreg:/usr/local/lib
ENTRYPOINT ["reg_aladin"]

79
scripts/02_build_m6_dataset.py
View file

@ -1,8 +1,9 @@
"""Build M6-2025 manifests: T1c volumes, labeled when GTV exists. """Build M6-2025 manifests: T1c volumes, labeled when GTV exists.
GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on a head box GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on the full
located in CT space; multi-start initial tilts because inter-session head tilt volumes -- no mask, no crop / pad; multi-start initial tilts because
reaches ~30 deg and tilted MI local minima can outscore the true alignment). inter-session head tilt reaches ~30 deg and tilted MI local minima can
outscore the true alignment).
The inverse transform then warps CT and GTV (nearest) into T1c native space. The inverse transform then warps CT and GTV (nearest) into T1c native space.
Acceptance requires the warped CT skull ring to be concentric with the MR head Acceptance requires the warped CT skull ring to be concentric with the MR head
silhouette (crop-independent check). silhouette (crop-independent check).
@ -105,6 +106,10 @@ def _skull_ring_gap(warped_ct, ha):
gap between the MR head boundary and the CT skull ring (p97 of positive gap between the MR head boundary and the CT skull ring (p97 of positive
CT) at two levels. Returns {"mid_mm", "low_mm", "std_mm"} or None when no CT) at two levels. Returns {"mid_mm", "low_mm", "std_mm"} or None when no
skull ring is found. Aligned = small gap with low angular variance. skull ring is found. Aligned = small gap with low angular variance.
The ring test is ``>=`` (not ``>``) so it also works on a CT clipped to a
brain window (e.g. [0, 80] HU), where the p97 threshold sits exactly on
the clipped-skull atom and a strict comparison would give an empty ring.
""" """
posw = warped_ct[warped_ct > 0] posw = warped_ct[warped_ct > 0]
if not posw.size: if not posw.size:
@ -119,7 +124,7 @@ def _skull_ring_gap(warped_ct, ha):
sl = sil[iz] sl = sil[iz]
if sl.sum() < 500: if sl.sum() < 500:
return None, None return None, None
sk = warped_ct[iz] > bthr sk = warped_ct[iz] >= bthr
if not sk.any(): if not sk.any():
return None, None return None, None
hy, hx = np.nonzero(sl) hy, hx = np.nonzero(sl)
@ -171,14 +176,11 @@ def _valid_pyramid(size, shrink, smooth, min_dim=6):
def register_ct_to_t1c(ct, t1c, threads=12): def register_ct_to_t1c(ct, t1c, threads=12):
"""Rigid registration with CT as fixed and T1c as moving, Mattes MI on """Rigid registration with CT as fixed and T1c as moving, Mattes MI on the
head-box images. full volumes (no mask, no crop / pad).
Both inputs share physical space (oriented). The registration grid is a head Both inputs share physical space (oriented). The transform is computed on
box located in CT space from the top of the head (the CTs are head-neck the full CT / T1c volumes, so it applies directly to any resampling.
exams, so an MR-head crop would show mostly neck/air on the CT); the box
excludes most air so MI is computed on tissue. The box is in physical
space, so the resulting transform applies directly to the full volumes.
Returns (transform, info): transform is the registration transform, mapping Returns (transform, info): transform is the registration transform, mapping
T1c (moving) coordinates -> CT (fixed) coordinates; resampling CT/TV onto T1c (moving) coordinates -> CT (fixed) coordinates; resampling CT/TV onto
the T1c grid with it applies the inverse transform to CT-space points, the T1c grid with it applies the inverse transform to CT-space points,
@ -198,38 +200,9 @@ def register_ct_to_t1c(ct, t1c, threads=12):
from src.common import head_mask_from_image from src.common import head_mask_from_image
ha = sitk.GetArrayFromImage(t1c) ha = sitk.GetArrayFromImage(t1c)
hm = head_mask_from_image(ha) hm = head_mask_from_image(ha)
# Registration grid = a head box located in CT space. The CTs are head-neck # full volumes, no mask / crop / pad: fixed = CT, moving = T1c
# (often head-neck-chest) exams, so cropping to the MR head box would leave fixed, moving = ct, t1c
# mostly neck/chest/air on the fixed image and MI locks onto the wrong
# structures. Vertex = topmost tissue slice with 10mm of continuity below
# (skips small artifacts floating above the head; bone HU is unreliable
# here, so tissue > -150 HU is used). x/y from the tissue extent of the
# 150mm band below the vertex. Box covers the whole head (x/y: extent +
# 15mm margin, z: +10 to -215 from vertex).
ca = sitk.GetArrayFromImage(ct)
tissue = ca > -150
csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
vz = _ct_vertex_slice(ct)
if vz is not None:
ztop = co[2] + vz * csp[2]
band = tissue[max(0, vz - int(round(150 / csp[2]))):vz + 1]
cols = np.nonzero(band.any(0))
box_lo = (co[0] + cols[1].min() * csp[0] - 15,
co[1] + cols[0].min() * csp[1] - 15, ztop - 215)
box_hi = (co[0] + cols[1].max() * csp[0] + 15,
co[1] + cols[0].max() * csp[1] + 15, ztop + 10)
else: # no head found: fall back to the whole CT volume
gs = ct.GetSize()
box_lo = tuple(float(x) for x in co)
box_hi = (co[0] + gs[0] * csp[0], co[1] + gs[1] * csp[1], co[2] + gs[2] * csp[2])
tsp = np.array(t1c.GetSpacing()) tsp = np.array(t1c.GetSpacing())
shape = tuple(int(np.ceil((box_hi[d] - box_lo[d]) / tsp[d])) + 1 for d in range(3))
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
ref.SetSpacing(t1c.GetSpacing())
ref.SetDirection(t1c.GetDirection())
ref.SetOrigin(tuple(float(x) for x in box_lo))
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
stages = [] stages = []
@ -261,8 +234,21 @@ def register_ct_to_t1c(ct, t1c, threads=12):
# rank candidates. Run the full pyramid from the centered start; if its # rank candidates. Run the full pyramid from the centered start; if its
# full-resolution geometric concentricity (skull-ring gap) is poor, repeat # full-resolution geometric concentricity (skull-ring gap) is poor, repeat
# from a grid of initial tilts and keep the best by gap. # from a grid of initial tilts and keep the best by gap.
# Euler rotation center = the CT head centroid (the 150 mm band below the
# vertex, same head definition as before); off maps it onto the MR head
# centroid. Both are init-only quantities; the registration itself runs on
# the full volumes.
ca = sitk.GetArrayFromImage(ct)
csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
tissue = ca > -150
vz = _ct_vertex_slice(ct)
if vz is not None:
z0 = max(0, vz - int(round(150 / csp[2])))
idx = np.array([z0, 0, 0]) + np.argwhere(tissue[z0:vz + 1]).mean(0)
else:
idx = np.argwhere(tissue).mean(0)
ct_c = np.array(co) + idx * csp
mr_c = np.array(t1c.GetOrigin()) + np.argwhere(hm).mean(0) * tsp mr_c = np.array(t1c.GetOrigin()) + np.argwhere(hm).mean(0) * tsp
ct_c = np.array([(box_lo[d] + box_hi[d]) / 2.0 for d in range(3)])
off = (float(mr_c[0] - ct_c[0]), float(mr_c[1] - ct_c[1]), float(mr_c[2] - ct_c[2])) off = (float(mr_c[0] - ct_c[0]), float(mr_c[1] - ct_c[1]), float(mr_c[2] - ct_c[2]))
rad = np.pi / 180.0 rad = np.pi / 180.0
@ -276,7 +262,9 @@ def register_ct_to_t1c(ct, t1c, threads=12):
return tr0, gap return tr0, gap
def is_good(g): def is_good(g):
return bool(g) and g["std_mm"] <= 12.0 and 0.0 <= g["mid_mm"] <= 20.0 # |mid| (not 0 <= mid): full-volume registration can land the CT skull
# ring slightly inside the MR silhouette, which is equally concentric
return bool(g) and g["std_mm"] <= 12.0 and abs(g["mid_mm"]) <= 20.0
tr_reg, gap = full_run(0, 0, 0) tr_reg, gap = full_run(0, 0, 0)
if not is_good(gap): if not is_good(gap):
@ -325,7 +313,6 @@ def register_ct_to_t1c(ct, t1c, threads=12):
"offset": [round(float(x), 3) for x in t0], "scale": scale}, "offset": [round(float(x), 3) for x in t0], "scale": scale},
"coverage_head_box": round(coverage, 3), "coverage_head_box": round(coverage, 3),
"skull_ring_gap_mm": skull_gap, "skull_ring_gap_mm": skull_gap,
"ct_head_box_mm": [[round(float(x), 1) for x in box_lo], [round(float(x), 1) for x in box_hi]],
"seconds_total": round(time.time() - t_start, 1)} "seconds_total": round(time.time() - t_start, 1)}
print("reg done %.0fs (rigid %.1fs) mi %.3f -> %.3f cov %.2f%s" % ( print("reg done %.0fs (rigid %.1fs) mi %.3f -> %.3f cov %.2f%s" % (
info["seconds_total"], sum(s["seconds"] for s in stages), info["seconds_total"], sum(s["seconds"] for s in stages),
@ -373,7 +360,7 @@ def reg_worker(row):
tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1)) tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1))
tv_in_ct = tv_in_ct_frac >= 0.5 # slab must straddle the lesion tv_in_ct = tv_in_ct_frac >= 0.5 # slab must straddle the lesion
gap = info["skull_ring_gap_mm"] gap = info["skull_ring_gap_mm"]
if gap and 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0: if gap and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
align = "concentric" align = "concentric"
else: else:
align = "not_concentric" align = "not_concentric"

43
scripts/bench_registration.py
View file

@ -20,6 +20,10 @@ higher = better in-sync = better aligned. The production geometric QA
- 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).
`--clip-ct LO HI` (e.g. `--clip-ct 0 80`, brain window) clips the CT to
[LO, HI] before registration AND before scoring, so a clipped vs unclipped
A/B compares the same input to both the toolkit and the MI score.
Because some toolkits are not pip-installable into every Because some toolkits are not pip-installable into every
`longitudinal-313` env (antspyx needs a cp310-cp313 wheel, 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
@ -93,6 +97,14 @@ def _m6():
return RP._load_m6() return RP._load_m6()
def _clip_image(img, lo, hi):
"""Clip image intensities to [lo, hi] (geometry preserved)."""
a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
out = sitk.GetImageFromArray(a)
out.CopyInformation(img)
return out
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# case selection # case selection
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
@ -152,7 +164,7 @@ def evaluate(ct, t1c, tv, T):
if nvox > 0: if nvox > 0:
tvd = ndimage.binary_dilation(warped_tv > 0.5, iterations=3) tvd = ndimage.binary_dilation(warped_tv > 0.5, iterations=3)
tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1)) tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1))
ok_gap = bool(gap) and 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0 ok_gap = bool(gap) and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0
accepted = nvox > 0 and tv_in_ct_frac >= 0.5 and ok_gap accepted = nvox > 0 and tv_in_ct_frac >= 0.5 and ok_gap
return { return {
"skull_ring_gap_mm": gap, "skull_ring_gap_mm": gap,
@ -237,7 +249,9 @@ def _call(pipeline, ct, t1c, tv, args):
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": if pipeline == "niftyreg":
return RP.register_niftyreg(ct, t1c, threads=args.threads) return RP.register_niftyreg(ct, t1c, threads=args.threads,
backend=args.niftyreg_backend,
gpu=args.niftyreg_gpu)
raise KeyError(pipeline) raise KeyError(pipeline)
@ -255,6 +269,11 @@ def _run_case(pipeline, case, args):
"alignment": "ct_incomplete", "accepted": False, "alignment": "ct_incomplete", "accepted": False,
"note": "CT head incomplete (data issue, not a pipeline failure)"}) "note": "CT head incomplete (data issue, not a pipeline failure)"})
return out return out
if args.clip_ct is not None:
# clip BEFORE registration AND scoring so the A/B compares like for like
lo, hi = args.clip_ct
ct = _clip_image(ct, lo, hi)
out["clip_ct"] = [lo, hi]
t0 = time.time() t0 = time.time()
T, meta = _call(pipeline, ct, t1c, tv, args) T, meta = _call(pipeline, ct, t1c, tv, args)
out["reg_wall_s"] = round(time.time() - t0, 1) out["reg_wall_s"] = round(time.time() - t0, 1)
@ -380,8 +399,8 @@ 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)
regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels] regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") 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") or {}).get("median") or 0 for p in labels]
fig, axes = plt.subplots(1, 2, figsize=(10, 4)) fig, axes = plt.subplots(1, 2, figsize=(10, 4))
colors = ["#4878CF", "#6ACC65"] colors = ["#4878CF", "#6ACC65"]
for ax, (title, vals, fmt, ylab) in zip(axes, [ for ax, (title, vals, fmt, ylab) in zip(axes, [
@ -487,6 +506,22 @@ def parse_args():
ap.add_argument("--impact-samples", type=int, default=5000, ap.add_argument("--impact-samples", type=int, default=5000,
help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the " help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
"elastix default)") "elastix default)")
ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"],
default="auto",
help="niftyreg reg_aladin build: docker = official v2.1.1 in the "
"niftyreg:2.1.1-ubuntu image (host glibc is too old for the "
"v2.x binaries), conda = conda-forge 1.5.69 on PATH, "
"auto (default) = docker when the image is present else conda")
ap.add_argument("--niftyreg-gpu", type=int, default=None,
help="niftyreg CUDA device index (docker backend; default: "
"auto-spread over the host GPUs by pid, -1 forces CPU "
"-platf 0; on the GPU the run uses 'docker run --gpus all "
"-platf 1' with linear interpolation, the v2.1.1 CUDA "
"limit)")
ap.add_argument("--clip-ct", nargs=2, type=float, default=None, metavar=("LO", "HI"),
help="clip the CT (HU) to [LO, HI] BEFORE registration and scoring "
"(e.g. --clip-ct 0 80 = brain window, skull/air removed; "
"default: no clipping)")
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")

315
src/registration_pipelines.py
View file

@ -23,15 +23,26 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
Pipelines (listed alphabetically) Pipelines (listed alphabetically)
--------------------------------- ---------------------------------
All pipelines register on the **full volumes** (no head box) and use each Shared rules for every pipeline:
toolkit's **native automatic transform initialization** (no manual
``centroid_offset``): elastix / impact use elastix's - **No mask, no crop / pad**: registration always 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
(no CT present) the first argument is kept as fixed (see
:func:`fixed_moving`).
- **CT brain-window clip for the pipelines that benefit from it**: per the
CT-clipping A/B (see REGISTRATION_BENCHMARK.md), ``niftyreg`` and
``synthmorph`` clip the CT to [0, 80] HU before registration
(``clip_ct`` parameter, default on, MR inputs never clipped); the
Mattes-MI pipelines (antspyx / elastix / simpleitk) register the CT as-is.
Each pipeline uses its toolkit's **native automatic transform initialization**
(no manual ``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, niftyreg default center-of-mass alignment, simpleitk starts from identity, niftyreg
(reg_aladin) initializes from the image centres, and synthmorph needs none (reg_aladin) initializes from the image centres, and synthmorph needs none
(global deep method). All use **fixed = CT, moving = T1c** except IMPACT, (global deep method).
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
@ -42,9 +53,9 @@ TotalSegmentator metric -- see ``register_impact``).
- ``register_impact`` : IMPACT (semantic-similarity metric) inside the - ``register_impact`` : IMPACT (semantic-similarity metric) inside the
``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript
TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF, TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF,
GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix GPU (Jacobian mode). full volumes, fixed=CT/moving=T1c (no mask -- shared
``AutomaticTransformInitialization`` (CenterOfGravity) provides the start. rules above); elastix ``AutomaticTransformInitialization`` (CenterOfGravity)
Needs Docker + the image + a model asset. provides the start. Needs Docker + the image + a model asset.
- ``register_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` -- - ``register_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
symmetric block-matching + trimmed least squares (Modat et al., J. Med. symmetric block-matching + trimmed least squares (Modat et al., J. Med.
Img. 2014), NCC-based block similarity (cross-modality tolerant). Rigid Img. 2014), NCC-based block similarity (cross-modality tolerant). Rigid
@ -123,6 +134,50 @@ def _warp_ct(ct, t1c, T):
return sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear)) return sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
def _is_ct(img):
"""Modality heuristic: CT is in Hounsfield units and carries air at large
negative values; MR intensities are non-negative."""
return float(sitk.GetArrayFromImage(img).min()) < -50.0
def _modality_name(img):
return "CT" if _is_ct(img) else "MR"
def fixed_moving(img_a, img_b):
"""Shared fixed-image rule: the CT is always the fixed image. When both
inputs are MR (or both CT) the given order is kept, i.e. ``img_a`` is the
fixed image. Returns ``(fixed, moving)``.
Note: a CT that has been intensity-clipped to a non-negative window (e.g.
[0, 80] HU) is no longer detectable as CT, so the given order is kept --
callers should pass the CT first (as all pipelines do)."""
if _is_ct(img_b) and not _is_ct(img_a):
return img_b, img_a
return img_a, img_b
def _clip_intensity(img, lo, hi):
"""Clip image intensities to [lo, hi] (geometry preserved)."""
a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
out = sitk.GetImageFromArray(a)
out.CopyInformation(img)
return out
def _maybe_clip_ct(fixed, clip_ct):
"""Apply the pipeline's CT brain-window clip to the fixed image.
The A/B benchmark (``--clip-ct 0 80``) showed clipping the CT to a brain
window helps the block-matching / deep pipelines (niftyreg, synthmorph)
and hurts the Mattes-MI ones (antspyx, elastix, simpleitk), so those two
pipelines clip the CT by default. Returns ``(fixed, clipped: bool)``;
MR fixed images (both-MR pairs) are never clipped."""
if clip_ct and _is_ct(fixed):
return _clip_intensity(fixed, *clip_ct), True
return fixed, False
def fix_orientation(ct, t1c, cand): def fix_orientation(ct, t1c, cand):
"""Return (best_T, inverted) choosing between ``cand`` and ``cand^-1``. """Return (best_T, inverted) choosing between ``cand`` and ``cand^-1``.
@ -203,8 +258,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
except AttributeError: except AttributeError:
pass pass
ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32) ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
# full volumes, fixed = CT, moving = T1c (no head box) # full volumes, no mask / crop / pad; CT is the fixed image (or, if both
fixed, moving = ct, t1c # inputs are MR, the first argument is)
fixed, moving = fixed_moving(ct, t1c)
shape = np.array(ct.GetSize()) shape = np.array(ct.GetSize())
# Identity start (no centroid init); the multi-resolution Mattes-MI pyramid # Identity start (no centroid init); the multi-resolution Mattes-MI pyramid
# recovers the gross inter-session offset. The Euler center is the fixed # recovers the gross inter-session offset. The Euler center is the fixed
@ -236,8 +292,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
"toolkit": f"SimpleITK {sitk.Version_VersionString()}", "toolkit": f"SimpleITK {sitk.Version_VersionString()}",
"metric": "MattesMutualInformation (full volumes)", "metric": "MattesMutualInformation (full volumes)",
"transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, " "transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
f"{iters_full} iters/level, identity start)"), f"{iters_full} iters/level, identity start)"),
"convention": "T1c->CT (benchmark budget)", "auto_inverted": False, "convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
"fixed": _modality_name(fixed), "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm, "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
"skull_ring_gap_mm": gap, "skull_ring_gap_mm": gap,
} }
@ -315,9 +372,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
after a successful registration. after a successful registration.
""" """
itk = _elastix_api() itk = _elastix_api()
# full volumes, fixed = CT, moving = T1c (no head box) # full volumes, no mask / crop / pad; CT is the fixed image
fixed = ct fixed, moving = fixed_moving(ct, t1c)
moving = t1c
params = { params = {
"ImageDimension": 3, "ImageDimension": 3,
"FixedInternalImagePixelType": "float", "FixedInternalImagePixelType": "float",
@ -373,7 +429,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
"metric": "MattesMutualInformation (full volumes)", "metric": "MattesMutualInformation (full volumes)",
"transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)", "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)",
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1), "fixed": _modality_name(fixed), "moving": _modality_name(moving),
"seconds_reg": round(dt, 1),
} }
return T, meta return T, meta
@ -408,8 +465,8 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
aff_iterations = tuple(aff_iterations) if aff_iterations else (500, 500, 200) aff_iterations = tuple(aff_iterations) if aff_iterations else (500, 500, 200)
aff_shrink_factors = tuple(aff_shrink_factors) if aff_shrink_factors else (4, 2, 1) aff_shrink_factors = tuple(aff_shrink_factors) if aff_shrink_factors else (4, 2, 1)
aff_smoothing_sigmas = tuple(aff_smoothing_sigmas) if aff_smoothing_sigmas else (3, 2, 0) aff_smoothing_sigmas = tuple(aff_smoothing_sigmas) if aff_smoothing_sigmas else (3, 2, 0)
# full volumes, fixed = CT, moving = T1c (no head box) # full volumes, no mask / crop / pad; CT is the fixed image
fx, mv = ct, t1c fx, mv = fixed_moving(ct, t1c)
import time import time
t0 = time.time() t0 = time.time()
# No manual centroid init: antspyx's default initial_transform=None computes # No manual centroid init: antspyx's default initial_transform=None computes
@ -431,7 +488,8 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
"pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", "pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
"metric": f"ANTs {aff_metric} (multi-resolution, full volumes)", "metric": f"ANTs {aff_metric} (multi-resolution, full volumes)",
"transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)", "transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)",
"auto_inverted": inverted, "aff_iterations": list(aff_iterations), "auto_inverted": inverted, "fixed": _modality_name(fx), "moving": _modality_name(mv),
"aff_iterations": list(aff_iterations),
"aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1), "aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
} }
return T, meta return T, meta
@ -484,7 +542,8 @@ def _lta_to_sitk(path):
def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None, def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
extent=256, keep_scratch=False): extent=256, keep_scratch=False,
clip_ct=(0.0, 80.0)):
"""SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image """SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image
``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``, ``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``,
pretrained ``synthmorph.rigid.1.h5`` weights). pretrained ``synthmorph.rigid.1.h5`` weights).
@ -492,7 +551,10 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware, SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware,
acquisition-agnostic deep registration: each volume is normalized to an acquisition-agnostic deep registration: each volume is normalized to an
isotropic 1-mm LIA network space and the (pretrained) network estimates the isotropic 1-mm LIA network space and the (pretrained) network estimates the
transform, so no skull-stripping or intensity preprocessing is needed. Being transform, so no skull-stripping or intensity preprocessing is needed --
except the one it benefits from: the CT is clipped to a brain window
(``clip_ct=(0.0, 80.0)`` by default, per the CT-clipping A/B) before
registration; pass ``clip_ct=None`` to register the CT as-is. Being
a global deep method it recovers the gross inter-session offset on its own, a global deep method it recovers the gross inter-session offset on its own,
so **no centroid init is used**. The pipeline so **no centroid init is used**. The pipeline
- runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible - runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible
@ -524,11 +586,15 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
# per-case scratch on a local (non-NFS) tmp fs # per-case scratch on a local (non-NFS) tmp fs
tmp = tempfile.mkdtemp(prefix="synthmorph_reg_") tmp = tempfile.mkdtemp(prefix="synthmorph_reg_")
try: try:
# fixed = CT, moving = T1c (shared convention for all pipelines) # full volumes, no mask / crop / pad; CT is the fixed image
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32), fixed, moving = fixed_moving(ct, t1c)
os.path.join(tmp, "fixed.nii.gz"), True) # CT fixed_name = _modality_name(fixed) # before clipping (a clipped CT
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32), # is no longer intensity-detectable)
os.path.join(tmp, "moving.nii.gz"), True) # T1c fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
os.path.join(tmp, "fixed.nii.gz"), True)
sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
os.path.join(tmp, "moving.nii.gz"), True)
# run as the host uid/gid so the container-written LTA is host-owned # run as the host uid/gid so the container-written LTA is host-owned
# (clean scratch removal) and TF still reaches the GPU. # (clean scratch removal) and TF still reaches the GPU.
@ -570,10 +636,11 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, " f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, "
f"extent {extent})"), f"extent {extent})"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "CT", "moving": "T1c", "fixed": fixed_name, "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "seconds_reg": round(dt, 1),
"docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu, "docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu,
"extent": extent, "extent": extent,
"clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
} }
if keep_scratch: if keep_scratch:
meta["scratch"] = tmp meta["scratch"] = tmp
@ -708,7 +775,8 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
'(FixedInternalImagePixelType "float")', '(FixedInternalImagePixelType "float")',
'(MovingInternalImagePixelType "float")', '(MovingInternalImagePixelType "float")',
'(UseDirectionCosines "true")', '(UseDirectionCosines "true")',
"", ]
L += [
'(Registration "MultiMetricMultiResolutionRegistration")', '(Registration "MultiMetricMultiResolutionRegistration")',
f"(NumberOfResolutions {n})", f"(NumberOfResolutions {n})",
f"(MaximumNumberOfIterations {iters})", f"(MaximumNumberOfIterations {iters})",
@ -748,7 +816,7 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000, 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=CT, moving=T1c) in the Docker image
``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric). ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic* IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
@ -765,12 +833,17 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level ``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level
scale pyramid (no image resampling). scale pyramid (no image resampling).
NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other Like every other pipeline this runs with **fixed=CT, moving=T1c** on the
pipelines). The CT FOV (~300x300x251 mm) is much larger than the T1c crop full volumes: no mask, no crop / pad (shared rules above). The IMPACT
(~220x220x150 mm); with CT as the fixed image, too many CT samples map metric draws random spatial samples from the fixed (CT) buffer and
outside the T1c buffer and the IMPACT metric aborts ("too many samples map evaluates them only where they map inside the moving (T1c) buffer; when
outside moving image buffer"). The returned transform still satisfies the too few samples are in-buffer, elastix aborts ("too many samples map
unified contract (T1c->CT). outside moving image buffer"). On the M6 head-neck CTs (FOV
~300x300x251 mm) vs the tight T1c head crop (~220x220x150 mm) only ~14%
of the samples are in-buffer at the CenterOfGravity auto-init, so the
M6 runs abort -- IMPACT needs a fixed volume whose FOV is comparable to
the moving crop (e.g. a head-only CT). The returned transform satisfies
the unified contract (T1c->CT).
Budget: the default ``iters=2000`` (MaximumNumberOfIterations) / Budget: the default ``iters=2000`` (MaximumNumberOfIterations) /
``samples=5000`` (NumberOfSpatialSamples, the elastix default) runs the ``samples=5000`` (NumberOfSpatialSamples, the elastix default) runs the
full benchmark budget for best results -- roughly 10x the authors' full benchmark budget for best results -- roughly 10x the authors'
@ -792,13 +865,10 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
raise PipelineUnavailable(f"impact: {note}") raise PipelineUnavailable(f"impact: {note}")
model_host, model_rel = _impact_model_host(model) model_host, model_rel = _impact_model_host(model)
# full volumes, fixed = T1c, moving = CT (no head box, no manual init -- # full volumes, no mask / crop / pad; CT is the fixed image (or, if both
# elastix's AutomaticTransformInitialization aligns the image centers). # inputs are MR, the first argument). No manual init -- elastix's
# T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger # AutomaticTransformInitialization aligns the image centers.
# than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples fixed, moving = fixed_moving(ct, t1c)
# map outside the T1c buffer and the IMPACT metric aborts.
fixed = t1c
moving = ct
# GPU: auto-spread across the host GPUs by pid, or forced value # GPU: auto-spread across the host GPUs by pid, or forced value
if gpu is None: if gpu is None:
@ -806,8 +876,8 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
gpu = (os.getpid() % ngpu) if ngpu else -1 gpu = (os.getpid() % ngpu) if ngpu else -1
# per-level feature-extraction scales (recommended 6/3/1.5/1 mm), scaled so # per-level feature-extraction scales (recommended 6/3/1.5/1 mm), scaled so
# the finest level matches the box's (T1c) base spacing # the finest level matches the moving image's base spacing
base = float(np.median(t1c.GetSpacing())) base = float(np.median(moving.GetSpacing()))
n = len(scales) n = len(scales)
voxel_sizes = [base * s for s in scales] voxel_sizes = [base * s for s in scales]
# image pyramid: IMPACT's docs recommend NO image resampling (all-1 schedules) # image pyramid: IMPACT's docs recommend NO image resampling (all-1 schedules)
@ -870,9 +940,9 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}") raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
# the optimized Euler (world frame) is the full combination transform # the optimized Euler (world frame) is the full combination transform
# (AutomaticTransformInitialization is folded in by elastix); it maps # (AutomaticTransformInitialization is folded in by elastix); elastix
# fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT); # maps moving(T1c) -> fixed(CT), which is the contract direction;
# fix_orientation resolves the direction. # fix_orientation verifies it against the skull ring.
T_main = sitk.ReadTransform(itk_tf) T_main = sitk.ReadTransform(itk_tf)
T, inverted = fix_orientation(ct, t1c, T_main) T, inverted = fix_orientation(ct, t1c, T_main)
finally: finally:
@ -886,7 +956,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
"transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; " "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; "
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "T1c", "moving": "CT", "fixed": _modality_name(fixed), "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "seconds_reg": round(dt, 1),
"docker_image": IMPACT_IMAGE, "model": model_rel, "gpu": gpu, "docker_image": IMPACT_IMAGE, "model": model_rel, "gpu": gpu,
"iters": iters, "samples": samples, "iters": iters, "samples": samples,
@ -901,6 +971,9 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
NIFTYREG_DOCKER_IMAGE = "niftyreg:2.1.1-ubuntu"
def _niftyreg_bin(): def _niftyreg_bin():
"""Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable.""" """Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable."""
b = shutil.which("reg_aladin") b = shutil.which("reg_aladin")
@ -913,6 +986,26 @@ def _niftyreg_bin():
"re-run in the activated env.") "re-run in the activated env.")
def _niftyreg_docker_ok():
"""Return (ok, note): is the NiftyReg v2.1.1 Docker image available?
The official KCL-BMEIS v2.x release binaries need GLIBC_2.34, which this
host (Ubuntu 20.04, glibc 2.31) lacks; the ``niftyreg:2.1.1-ubuntu`` image
(built from the release zip, see niftyreg/Dockerfile) runs them instead.
"""
if shutil.which("docker") is None:
return False, "docker CLI not found on PATH"
try:
r = subprocess.run(["docker", "image", "inspect", NIFTYREG_DOCKER_IMAGE],
capture_output=True, text=True, timeout=30)
if r.returncode != 0:
return False, (f"docker image {NIFTYREG_DOCKER_IMAGE} not built "
f"(run: docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/)")
except Exception as e: # noqa
return False, f"docker image check failed ({type(e).__name__}: {e})"
return True, NIFTYREG_DOCKER_IMAGE
def _niftyreg_version(binpath): def _niftyreg_version(binpath):
try: try:
r = subprocess.run([binpath, "--version"], capture_output=True, r = subprocess.run([binpath, "--version"], capture_output=True,
@ -938,7 +1031,8 @@ def _parse_niftyreg_affine(path):
return T return T
def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False): def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False,
backend="auto", gpu=None, clip_ct=(0.0, 80.0)):
"""NiftyReg rigid CT->T1c (reference=CT, floating=T1c) via ``reg_aladin``. """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 reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global
@ -951,32 +1045,96 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
transforms concurrently; ``-maxit 5`` (the tool default) is too coarse for transforms concurrently; ``-maxit 5`` (the tool default) is too coarse for
these repositioned pairs -- the benchmark uses 15. these repositioned pairs -- the benchmark uses 15.
Full volumes (no head box); default initialisation (image centres); The CT is clipped to a brain window (``clip_ct=(0.0, 80.0)`` by default,
``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3`` (cubic) internal per the CT-clipping A/B) before registration; pass ``clip_ct=None`` to
warping. The output affine is documented as ``Affine * Reference = register the CT as-is. MR inputs are never clipped.
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 Full volumes (no head box, no mask / crop / pad); default initialisation
(``conda install -c conda-forge niftyreg``). (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.
``backend`` selects the reg_aladin build: ``"docker"`` runs the official
KCL-BMEIS **v2.1.1** release in the ``niftyreg:2.1.1-ubuntu`` Docker image
(built from NiftyReg-Ubuntu-v2.1.1.zip, see niftyreg/Dockerfile -- the
v2.x binaries need GLIBC_2.34, which this host lacks); ``"conda"`` uses
the conda-forge ``reg_aladin`` on PATH (NiftyReg 1.5.69); ``"auto"``
(default) prefers Docker when the image is present and falls back to
conda.
``gpu`` selects the CUDA device for the Docker backend: ``None`` (default)
auto-spreads over the host GPUs by pid, ``-1`` forces CPU. On the GPU the
run uses ``docker run --gpus all`` + ``-platf 1 -gpuid N`` and
``-interp 1`` (the v2.1.1 CUDA path supports linear interpolation only);
on the CPU it uses ``-platf 0`` with ``-interp 3`` (cubic), as does the
conda backend. The GPU is ~20x faster than CPU on the full volumes.
Preconditions (else :class:`PipelineUnavailable`): the Docker image for
backend docker/auto, or ``reg_aladin`` on PATH (``conda install -c
conda-forge niftyreg``).
""" """
import time import time
binpath = _niftyreg_bin() if backend not in ("auto", "conda", "docker"):
# full volumes, reference = CT, floating = T1c (no head box) raise PipelineUnavailable(
f"niftyreg: unknown backend {backend!r} (use auto|conda|docker)")
dok, dnote = _niftyreg_docker_ok()
use_docker = backend == "docker" or (backend == "auto" and dok)
if use_docker:
if not dok:
raise PipelineUnavailable(f"niftyreg: {dnote}")
binpath = None
if gpu is None:
ngpu = _n_gpus()
gpu = (os.getpid() % ngpu) if ngpu else -1
use_gpu = gpu >= 0
platf = "1" if use_gpu else "0"
# the v2.1.1 CUDA path supports linear interpolation only
interp = "1" if use_gpu else "3"
toolkit = (f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE}, "
f"CUDA)" if use_gpu
else f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE})")
else:
binpath = _niftyreg_bin()
gpu, use_gpu, platf = -1, False, None
interp = "3"
toolkit = _niftyreg_version(binpath)
# full volumes, no mask / crop / pad; reference (fixed) = CT, floating = T1c
fixed, moving = fixed_moving(ct, t1c)
fixed_name = _modality_name(fixed) # before clipping (a clipped CT is no
# longer intensity-detectable)
fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
tmp = tempfile.mkdtemp(prefix="niftyreg_reg_") tmp = tempfile.mkdtemp(prefix="niftyreg_reg_")
dt = 0.0 dt = 0.0
try: try:
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32), sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
os.path.join(tmp, "ref.nii.gz"), True) os.path.join(tmp, "ref.nii.gz"), True)
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32), sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
os.path.join(tmp, "flo.nii.gz"), True) os.path.join(tmp, "flo.nii.gz"), True)
cmd = [binpath, "-ref", "ref.nii.gz", "-flo", "flo.nii.gz", flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
"-rigOnly", "-aff", "aff.txt", "-rigOnly", "-aff", "aff.txt",
"-ln", str(levels), "-maxit", str(maxit), "-ln", str(levels), "-maxit", str(maxit),
"-interp", "3", "-omp", str(threads), "-voff"] "-interp", interp, "-omp", str(threads)]
if platf is not None:
flags += ["-platf", platf]
if use_gpu:
flags += ["-gpuid", str(gpu)]
flags += ["-voff"]
if use_docker:
# reg_aladin is the image ENTRYPOINT; the work dir is /work.
# --gpus all: the v2.1.1 Ubuntu build ships libniftyreg_cuda.so
# (dlopen'd by reg_aladin with -platf 1); the NVIDIA container
# toolkit provides libcuda.so.1, and -gpuid selects the device.
cmd = ["docker", "run", "--rm"]
if use_gpu:
cmd += ["--gpus", "all"]
cmd += ["-v", f"{tmp}:/work", "-w", "/work",
NIFTYREG_DOCKER_IMAGE] + flags
else:
cmd = [binpath] + flags
t0 = time.time() t0 = time.time()
proc = subprocess.run(cmd, capture_output=True, text=True, proc = subprocess.run(cmd, capture_output=True, text=True,
cwd=tmp, timeout=7200) cwd=tmp, timeout=7200)
@ -999,14 +1157,16 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
meta = { meta = {
"pipeline": "niftyreg", "pipeline": "niftyreg",
"toolkit": _niftyreg_version(binpath), "toolkit": toolkit,
"backend": "docker" if use_docker else "conda",
"metric": "NCC block matching + trimmed least squares (symmetric, full volumes)", "metric": "NCC block matching + trimmed least squares (symmetric, full volumes)",
"transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, " "transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, "
f"{maxit} TLS iters/level, cubic interp)"), f"{maxit} TLS iters/level, cubic interp)"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "CT", "moving": "T1c", "fixed": fixed_name, "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "seconds_reg": round(dt, 1),
"levels": levels, "maxit": maxit, "levels": levels, "maxit": maxit, "interp": interp, "gpu": gpu,
"clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
} }
if keep_scratch: if keep_scratch:
meta["scratch"] = tmp meta["scratch"] = tmp
@ -1047,11 +1207,20 @@ def availability():
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": elif name == "niftyreg":
b = shutil.which("reg_aladin") b = shutil.which("reg_aladin")
if b: dok, dnote = _niftyreg_docker_ok()
out[name] = (True, _niftyreg_version(b)) if dok:
ver = f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE})"
if b:
ver += f" | conda fallback: {_niftyreg_version(b)}"
out[name] = (True, ver)
elif b:
out[name] = (True, _niftyreg_version(b)
+ " | docker v2.1.1 image missing (build: "
f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/)")
else: else:
out[name] = (False, "reg_aladin not on PATH " out[name] = (False, "reg_aladin not on PATH and no docker image; "
"(conda install -c conda-forge niftyreg)") "conda install -c conda-forge niftyreg or "
f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/")
elif name == "synthmorph": elif name == "synthmorph":
ok, note = _synthmorph_docker_ok() ok, note = _synthmorph_docker_ok()
if not ok: if not ok: