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

Compare commits

base: xfr:4ad972b00fec3fa578570a4d410f139d1df272b2
Branches Tags
xfr:main
..
compare: xfr:6b85e950e556941fb6a878585726abbf31e6fb72
Branches Tags
xfr:main

No commits in common. "4ad972b00fec3fa578570a4d410f139d1df272b2" and "6b85e950e556941fb6a878585726abbf31e6fb72" have entirely different histories.
4ad972b00f ... 6b85e950e5

3 changed files with 319 additions and 464 deletions
Show stats Expand all files Collapse all files

327
REGISTRATION_BENCHMARK.md
View file

@ -22,34 +22,21 @@ 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 use destination→source convention this `T` maps **T1c → CT**. All pipelines register
**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c, on the same CT-space head box (MI computed on head tissue, not neck/air) and,
moving = CT** (swapping its roles breaks the TotalSegmentator-based metric, for the single-init tools, from the **same centroid-aligned start**
which aborts with "Error in metric" or converges to a wrong rotation). The (`centroid_offset`), so the comparison is on the optimizer / metric / transform
intensity-metric pipelines register on the same CT-space head box (MI computed handling — not on who finds the ~30–60 mm gross offset first.
on head tissue, not neck/air); `synthmorph` uses the full volumes (it
normalizes each one to an isotropic network space internally). The local
optimizers that cannot recover the gross inter-session offset from identity
(`antspyx`, `impact`) start from a centroid-aligned translation
(`centroid_offset`); the others start from identity, so the comparison is on
each toolkit's native optimizer / metric / global-search capability.
## The five pipelines ## The four pipelines
Listed alphabetically. All use **fixed = CT, moving = T1c** except `impact`, | # | pipeline | toolkit | metric | transform | status in this env |
which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator |---|----------|---------|--------|-----------|--------------------|
metric). A pipeline that is a local optimizer and cannot recover the gross | 1 | `elastix` | ITK-Elastix 0.25.4 (pip `itk` + `itk-elastix`, elastix 5.3.1) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x replacement for `Rigid3DTransform`), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
inter-session offset from identity (`antspyx`, `impact`) starts from a | 2 | `simpleitk` | SimpleITK 2.5.6 | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
centroid-aligned translation; the others start from identity / their own global | 3 | `antspyx` | ANTsPy / ANTs 0.6.3 | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
search. | 4 | `deepreg` | Oxford DeepReg (CycleGAN) | CycleGAN + reg. loss (pretrained) | 3D displacement field | **not installed** (see below) |
| 5 | `impact` | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT) | **IMPACT semantic** (TotalSegmentator `TS/M730.pt` features, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales (no image resampling) | ✅ ran (main env, Docker + GPU) |
| # | pipeline | fixed / moving | centroid init | metric | transform | status in this env |
|---|----------|----------------|---------------|--------|-----------|--------------------|
| 1 | `antspyx` | CT / T1c | **yes** (required) | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
| 2 | `elastix` | CT / T1c | no (identity) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
| 3 | `impact` | **T1c / CT** | **yes** | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) |
| 4 | `simpleitk` | CT / T1c | no (identity) | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
| 5 | `synthmorph` | CT / T1c | no (global deep) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
`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
@ -98,95 +85,76 @@ Production QA, crop-independent:
| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) | | pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|---|---|---|---|---|---|---|---|---| |---|---|---|---|---|---|---|---|---|
| **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 18.9 / 20.1 | 31.7 / 31.3 | 10.9 | 0.94 | **0.436 / 0.393** | | **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.6 | 0.93 | **0.643 / 0.630** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 16.6 / 17.6 | 16.9 / 18.8 | 2.4 | 0.93 | **0.641 / 0.583** | | **antspyx** | ANTsPy 0.6.3 | 9 | **0/10 (0.00)** | 14.9 / 16.5 | 18.4 / 22.7 | 14.5 | 0.93 | **0.598 / 0.529** |
| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 13.1 / 16.8 | 13.2 / 18.6 | **576.0** | 0.93 | **0.469 / 0.482** | | **elastix** | ITK-Elastix 0.25.4 | 9 | **5/10 (0.50)** | 15.3 / 16.1 | 9.9 / 13.8 | 3.2 | 0.93 | **0.562 / 0.549** |
| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.7 / 16.0 | 13.5 / 15.7 | 7.6 | 0.93 | **0.645 / 0.632** | | **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 12.4 / 19.0 | 15.1 / 16.7 | **562.3** | 0.93 | **0.475 / 0.448** |
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 22.9 / 21.7 | 26.0 | 0.93 | **0.241 / 0.308** | | deepreg | — | — | — | — | — | — | — | — |
*(antspyx and impact are stochastic (random metric/patch subsampling), so their *(antspyx and impact are stochastic (random metric/patch subsampling), so their
accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.6 accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.4
min/case on GPU, orders of magnitude above the intensity-metric tools. min/case on GPU, orders of magnitude above the intensity-metric tools.)*
elastix, simpleitk and synthmorph are deterministic at fixed threads.)*
Per-case (mid / std, `OK` = accepted): Per-case (mid / std, `OK` = accepted):
| case | antspyx | elastix | impact | simpleitk | synthmorph | | case | simpleitk | antspyx | elastix | impact |
|---|---|---|---|---|---| |---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 18.9 / 43.8 | 27.0 / 29.3 | 8.3 / 13.2 | 26.2 / 29.6 | 10.4 / 23.7 | | m6_22IRRKU6_20240115 | 25.7 / 29.3 | 24.4 / 27.7 | 26.4 / 30.3 | **8.4 / 12.6** |
| m6_23XIMO4M_20230222 | 22.9 / 24.5 | 22.6 / 23.9 | 23.1 / 25.7 | 22.6 / 24.0 | 20.3 / 23.3 | | m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.4 / 24.8 | 21.1 / 23.3 | 21.4 / 24.1 |
| m6_24V62QQM_20230223 | 18.7 / 47.5 | 11.3 / 16.9 | 9.9 / 7.9 **OK** | 10.2 / 16.2 | 10.4 / 22.9 | | m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.5 / 18.4 | 10.6 / 8.2 **OK** | 10.0 / 8.3 **OK** |
| m6_24V62QQM_20231017 | 15.3 / 37.9 | 16.6 / 11.1 **OK** | 6.6 / 6.9 **OK** | 11.9 / 7.7 **OK** | 11.5 / 8.4 **OK** | | m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 14.9 / 15.5 | 11.0 / 6.1 **OK** | 11.1 / 8.8 **OK** |
| m6_24V62QQM_20240820 | 38.7 / 54.0 | 10.2 / 8.2 **OK** | 37.1 / 45.6 | 10.8 / 9.0 **OK** | 43.4 / 15.8 | | m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.3 / 17.5 | 15.3 / 14.0 | 54.4 / 15.1 |
| m6_26TYQI3C_20221027 | 11.2 / 3.7 **OK** | 11.7 / 3.7 **OK** | 23.1 / 13.0 | 11.1 / 2.6 **OK** | 11.8 / 2.9 **OK** | | m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 17.7 / 37.6 | 11.4 / 3.2 **OK** | 23.9 / 22.6 |
| m6_2BBN5XA4_20250411 | 20.6 / 24.2 | 25.3 / 30.8 | 10.1 / 6.5 **OK** | 19.1 / 8.8 **OK** | 52.4 / 53.0 | | m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** | 10.1 / 5.9 **OK** |
| m6_2C4TKI6H_20240506 | 15.3 / 14.2 | 14.2 / 15.3 | 13.1 / 17.4 | 12.7 / 13.5 | 16.7 / 28.2 | | m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.1 / 14.3 | 12.0 / 6.8 **OK** | 12.4 / 22.3 |
| m6_2CDNN7BH_20240102 | 19.4 / 31.7 | 19.5 / 30.3 | 19.6 / 30.9 | 19.3 / 30.1 | 18.9 / 16.8 | | m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 17.2 / 35.4 | 18.9 / 22.7 | 19.6 / 30.3 |
### Interpretation ### Interpretation
- **Translation is recovered well by all five** (mid gaps mostly 6–20 mm, - **Translation is recovered well by all four** (mid gaps mostly 6–20 mm,
within the `≤ 20 mm` band, except the grossly-drifted cases below). within the `≤ 20 mm` band). `tv-in-CT` is identical (~0.93): the warped GTV
`tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for lands inside the CT head for every pipeline, so the target is in the right
every pipeline, so the target is in the right place. place.
- **The binding failure is rotation** — `std > 12 mm` on the repositioned cases - **The binding failure is rotation** — `std > 12 mm` on the repositioned cases
(the ~30° inter-session tilt produces local minima that a single start (the ~30° inter-session tilt produces local minima that a single centered
(identity or centroid) can lock onto). This is exactly why production start can lock onto). This is exactly why production `02_build_m6_dataset.py`
`02_build_m6_dataset.py` uses a **multi-start tilt grid**. uses a **multi-start tilt grid**.
- **SimpleITK is the most reliable intensity-metric tool in this draw** - **Elastix is the most accurate of the standard single-init configs**
(4/10 accepted, best median mid 12.7), and **Elastix is the fastest** (5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 18.4) and
(2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s
median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c SimpleITK / 14.5 s ANTs): ASGD's adaptive step sizes climb out of small
convention the two settle into different local minima on the hard tilted misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223`
cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017` (std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on
(std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles (std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both
better, and vice-versa on the drift cases. The fixed/moving direction others accept — different optimizers settle in different local minima on the
materially affects where the ASGD/RSGD optimizers converge on the tilted cases. hard tilted cases.
- **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this - **ANTs' Rigid is the weakest of the intensity-metric tools** (0/10 in this
draw; it is stochastic, so the count varies): it frequently locks onto the draw; it is stochastic, so a prior draw gave 1/10): it frequently locks onto
tilted local minimum (median std 31.7, worst of all). It accepts only the tilted local minimum. Note `m6_26TYQI3C`: ANTs nailed the position (mid
`m6_26TYQI3C` (mid 11.2 / std 3.7) and misses `m6_24V62QQM_20230223` 17.7) but missed the orientation (std 37.6) — a textbook tilted local minimum.
(std 47.5) and `m6_24V62QQM_20231017` (std 37.9) that the other tools accept —
textbook tilted local minima.
- **IMPACT (semantic metric) is the standout on the hardest repositioned case** - **IMPACT (semantic metric) is the standout on the hardest repositioned case**
and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every
intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of intensity tool fails, std 28–30 mm) it gets **mid 8.4 / std 12.6** — the best
any pipeline and within ~1 mm of the acceptance threshold — because matching translation of any pipeline and within 0.6 mm of the acceptance threshold —
*anatomical semantics* (TotalSegmentator features) is far more discriminating because matching *anatomical semantics* (TotalSegmentator features) is far
for gross pose than raw cross-modality intensity. But the semantic landscape more discriminating for gross pose than raw cross-modality intensity. But the
is rugged: single-init ASGD lands in bad local minima on several cases, and semantic landscape is rugged: single-init ASGD lands in bad local minima on
two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0). 4 of 9 cases (std 15–30 mm), and one (`m6_24V62QQM_20240820`) even drifts the
Net 3/10 accepted — between ANTs and SimpleITK — at ~9.6 min/case, two orders translation (mid 54 mm). Net 3/10 accepted — between ANTs and SimpleITK — at
of magnitude slower than the intensity tools. ~9.4 min/case, two orders of magnitude slower than the intensity tools.
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.4 s - **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~3.2 s
(warm workers), SimpleITK ~7.6 s, ANTs ~10.9 s (CPU-bound), (warm workers; a worker's first case pays a ~10 s warm-up), SimpleITK ~8.6 s,
**SynthMorph ~26 s** (GPU: a single pretrained-model forward pass per case), ANTs ~14.5 s (CPU-bound), **IMPACT ~562 s** (GPU-bound: every ASGD iteration
**IMPACT ~576 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator runs the TotalSegmentator extractor on the sampled patches with backprop).
extractor on the sampled patches with backprop). - **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**:
- **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**: medians 0.643 (SimpleITK) / 0.598 (ANTs) / 0.562 (Elastix) / 0.475 (IMPACT),
medians 0.645 (SimpleITK) / 0.567 (Elastix) / 0.469 (IMPACT) / 0.436 (ANTs) / and it diverges sharply on the tilted cases — ANTs `m6_26TYQI3C` drops to
0.241 (SynthMorph), ~0.1 nats (rotation ~38 mm off) while SimpleITK holds 0.670. A high MI means
and it diverges sharply on the tilted cases — ANTs `m6_24V62QQM_20230223` the CT skull/tissue intensities fall on the matching MR structures after the
drops to ~0.06 nats (rotation ~47 mm off) while SimpleITK holds 0.703. A high inverse transform. IMPACT's lower median MI is consistent with its less
MI means the CT skull/tissue intensities fall on the matching MR structures reliable single-init convergence.
after the inverse transform. ANTs' and SynthMorph's low median MI is - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
consistent with their tilted-case convergence failures.
- **SynthMorph (deep, pretrained) is near-weakest of the five on these
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.241)
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
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 (RSGD wall time is thread-sensitive, convergence is not); the in-process
@ -199,15 +167,14 @@ Per-case (mid / std, `OK` = accepted):
## Reproducing ## Reproducing
```bash ```bash
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here) # main env — SimpleITK + Elastix + ANTsPy (all pip-installed here); deepreg
# reports unavailable and is recorded as such
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 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 --n-cases 10 \
--n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark --workers 4 --threads 24 --simpleitk-mode benchmark
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT # IMPACT (needs Docker + the pulled image + a model; see the Impact section
# sections below for the one-time setup) # below for the one-time setup)
python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \
--workers 4 --threads 24
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
@ -218,8 +185,7 @@ 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`, `--simpleitk-mode {benchmark,production}`,
`--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`, `--deepreg-ckpt <.ckpt>` (DeepReg), `--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`, `--impact-samples N`,
`--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below). `--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
@ -250,61 +216,35 @@ visuals agree. Generate for the available pipelines 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
python scripts/bench_registration.py \ python scripts/bench_registration.py \
--pipelines simpleitk elastix antspyx --n-cases 10 --qa-images --pipelines simpleitk elastix antspyx --n-cases 10 --qa-images
python scripts/bench_registration.py \
--pipelines synthmorph --n-cases 10 --qa-images
python scripts/bench_registration.py \ python scripts/bench_registration.py \
--pipelines impact --n-cases 10 --qa-images --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.)
## SynthMorph pipeline (Docker) — one-time setup ## Enabling the blocked pipeline
The `synthmorph` pipeline runs the [`freesurfer/synthmorph`](https://hub.docker.com/r/freesurfer/synthmorph) **DeepReg (Oxford)** — the PyPI `deepreg` (0.1.2) is an unrelated stub that pins
image: [SynthMorph](https://synthmorph.io) (Hoffmann et al., IEEE TMI 2022), an `numpy==1.18.5`; do not install it. DeepReg is not a one-shot rigid tool: it is a
anatomy-aware, acquisition-agnostic **deep** registration (TensorFlow CycleGAN image-to-image registration that must be **trained or run from a
`VxmAffineFeatureDetector` with the pretrained `synthmorph.rigid.1.h5` weights). pretrained checkpoint**. To benchmark it:
No skull-stripping or intensity preprocessing is needed: each volume is
normalized to an isotropic 1-mm LIA network space internally. It is the
second Docker-based pipeline and the only deep-learning one.
Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command):
```bash ```bash
# 1. Docker (with NVIDIA container toolkit so --gpus works) git clone https://github.com/OB-CVHU/DeepReg && cd DeepReg && pip install -e .
docker --version python scripts/bench_registration.py --pipelines deepreg --n-cases 10 \
--deepreg-ckpt /path/to/pretrained_3d.ckpt
# 2. Pull the image (~11 GB, one-time)
docker pull freesurfer/synthmorph
``` ```
No extra model assets: the rigid weights ship inside the image The adapter runs `deepreg test` to predict a voxel displacement field, converts
(`/freesurfer/models`). it to a SimpleITK `DisplacementFieldTransform`, and scores it like the others.
### How the adapter drives the image
1. **Rigid matrix model, no init** — `register -m rigid` on GPU when one is
visible (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces
CPU). SynthMorph is a global deep method, so it recovers the gross
inter-session offset on its own and needs **no** initial transform.
2. **Read-back** — the `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in shared world
space; with fixed=CT/moving=T1c the `-T` matrix is the T1c→CT contract
transform) is parsed into a SimpleITK `AffineTransform` and run through the
usual `fix_orientation` gate.
3. **GPU sharing** — the container is run as the host uid (clean scratch
removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
concurrent benchmark workers can share a GPU.
## 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)
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 the only pipeline that
Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps is GPU-based and the only one that is Docker-based.
**fixed = T1c, moving = CT** (the reverse of the others): swapping its roles
breaks the TotalSegmentator metric ("Error in metric" / wrong rotation).
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):
@ -355,77 +295,50 @@ to add the shared centered start:
## Notes / caveats ## Notes / caveats
- The benchmark compares each toolkit in its **native configuration**: the - The benchmark compares **standard single-init** configurations (centered
local optimizers that cannot recover the gross offset from identity start). Production uses `--simpleitk-mode production`'s multi-start tilt grid,
(`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`, which raises SimpleITK's acceptance at ~20× the cost; that is out of scope for
`elastix` and `synthmorph` start from identity / their own global search. the head-to-head but is the config the dataset actually ships with.
Production uses `--simpleitk-mode production`'s multi-start tilt grid, which - `deepreg` could not be executed here, so its row is "not installed", not
raises SimpleITK's acceptance at ~20× the cost; that is out of scope for the measured; its code path is implemented and unit-verified for import/convention.
head-to-head but is the config the dataset actually ships with. - `fix_orientation` (elastix/ants/deepreg/impact) chooses between a transform
- `synthmorph` is deterministic (one fixed pretrained model, no random
subsampling), so its per-case numbers are stable run to run — unlike ANTs /
IMPACT. Its low acceptance (2/10) is a genuine cross-modality / tilt
weakness, not variance.
- `fix_orientation` (elastix/ants/synthmorph/impact) chooses between a transform
and its inverse by skull-ring concentricity — a correctness gate that makes and its inverse by skull-ring concentricity — a correctness gate that makes
the result independent of each toolkit's transform-direction convention. With the result independent of each toolkit's transform-direction convention. In
the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the practice it never fired (all `auto_inverted: False`): elastix/ants/impact
inverse of the contract, so it reports `auto_inverted: True` (the gate returned the T1c→CT direction as expected.
corrects it); antspyx/synthmorph/impact return the contract direction
directly (`auto_inverted: False`).
- **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU):
- **Cross-modality / tilt weakness.** The `rigid` model is trained on brain
*MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly
defeated by the ~30° inter-session tilt (7/9 cases have `std > 15 mm`).
It is a legitimate data point for how far a pretrained MRI deep model
generalizes out-of-domain to CT — not a bug.
- **No centroid init needed.** Being a global deep method it recovers the
gross offset on its own; the `-i` initial-transform path is not used
(which also sidesteps this image's float64/float32 bug on the `-i` code
path).
- **GPU sharing.** The container runs as the host uid with
`TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
- **~26 s/case** (GPU: one pretrained forward pass), between ANTs (~11 s)
and IMPACT (~576 s).
- **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
3D transform is `EulerTransform` (6-DOF, `Rigid3DTransform` is gone from 3D transform is `EulerTransform` (6-DOF, `Rigid3DTransform` is gone from
the component set), the metric key is `Metric` (not `ImageMetric`), the the component set), the metric key is `Metric` (not `ImageMetric`), the
random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x
names and is not compatible with elastix ≤ 4.x parameter files. names and is not compatible with elastix ≤ 4.x parameter files.
- No initial transform is passed (identity start): the ASGD + 3-resolution - The initial transform (the shared `centroid_offset` start) is passed as an
Mattes-MI pyramid recovers the gross inter-session offset on its own, so `InitialTransformParameterObject` (a `TranslationTransform` with
the shared `centroid_offset` start is not required. An **output directory** `NumberOfParameters 3`); the C++ filter only wires it up when an
is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter **output directory is set** (`SetOutputDirectory`), so the adapter uses a
uses a scratch dir under `data/m6reg/elastix/itkelastix/`. scratch dir under `data/m6reg/elastix/itkelastix/`.
- `WriteResultImage` must be the string `"false"` (Python `False`); the - `WriteResultImage` must be the string `"false"` (Python `False`); the
default (absent) makes the filter raise "Could not read result image" default (absent) makes the filter raise "Could not read result image"
after a *successful* registration. after a *successful* registration.
- With fixed=CT/moving=T1c the combination transform is CT→T1c (the inverse - The result transform is recovered by point-mapping
of the contract), so `fix_orientation` reports `auto_inverted: True` and
inverts it. The result is recovered by point-mapping
(`GetCombinationTransform().TransformPoint`) into an affine — independent (`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 - Registration wall time is ~1.8 s/case at full budget (ASGD 500/500/250
500/500/250 over 3 resolutions). over 3 resolutions); each worker's first case pays a ~10 s warm-up, so the
median wall is 3.2 s and the mean is inflated by warm-up.
- **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0, - **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, so the recommended
500-iter / 2000-sample budget is ~9.6 min/case (vs 2–11 s for the other 500-iter / 2000-sample budget is ~9 min/case (vs 3–15 s for the intensity
tools). `--impact-iters` / `--impact-samples` trade accuracy for speed. tools). `--impact-iters` / `--impact-samples` trade accuracy for speed.
- **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT - **Rugged single-init convergence.** The semantic landscape has more local
follow the shared fixed=CT/moving=T1c convention: swapping its roles minima than Mattes-MI: the pipeline nails the hardest repositioned case
breaks the TotalSegmentator metric (elastix aborts with "Error in metric" (`m6_22IRRKU6`, mid 8.4 mm — best of all) but lands in bad minima on some
or converges to a wrong rotation). The returned transform still satisfies others. A multi-start or higher-iteration budget would likely raise its
the contract (T1c→CT). acceptance; that is out of scope for the single-init head-to-head.
- **Rugged single-init convergence.** The semantic landscape has more local
minima than Mattes-MI: the pipeline nails the hardest repositioned case
(`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some
others. A multi-start or higher-iteration budget would likely raise its
acceptance; that is out of scope for the single-init head-to-head.
- **Stochastic.** Jacobian mode resamples patches each iteration, so results - **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

43
scripts/bench_registration.py
View file

@ -1,7 +1,7 @@
"""Benchmark the 5 rigid CT->T1c registration pipelines on N M6 test cases. """Benchmark the 5 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 | deepreg | impact
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
@ -19,22 +19,22 @@ geometric QA used in production (scripts/02_build_m6_dataset.py):
Because some toolkits are not pip-installable into the main Because some toolkits are not pip-installable into the main
`longitudinal-313` env (antspyx is kept in an isolated env, see `longitudinal-313` env (antspyx is kept in an isolated env, see
REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under REGISTRATION_BENCHMARK.md; PyPI deepreg is a broken stub), the benchmark is
whichever conda env has a given toolkit installed, and it records only the env-portable: run it under whichever conda env has a given toolkit installed,
pipelines that are importable there. elastix runs in the main env through the and it records only the pipelines that are importable there. elastix runs in the main env through the ITK-Elastix
ITK-Elastix python API (pip ``itk`` + ``itk-elastix``). ``synthmorph`` and python API (pip ``itk`` + ``itk-elastix``). ``impact`` runs in the
``impact`` run in Docker images (``freesurfer/synthmorph`` and ``vboussot/elastix_impact`` Docker image (needs Docker + the image + a model,
``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see see REGISTRATION_BENCHMARK.md). The `--merge` mode then combines
REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the 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; impact needs Docker + the image;
# deepreg unavailable)
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
# synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md # impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
python scripts/bench_registration.py --pipelines synthmorph impact python scripts/bench_registration.py --pipelines impact
# any pipeline that is not importable in the current env is recorded as # any pipeline that is not importable in the current env is recorded as
# unavailable, so the script is env-portable; --merge combines the # unavailable, so the script is env-portable; --merge combines the
@ -225,10 +225,9 @@ def _call(pipeline, ct, t1c, tv, args):
return RP.register_elastix(ct, t1c, workdir=args.elastix_dir, threads=args.threads) return RP.register_elastix(ct, t1c, workdir=args.elastix_dir, threads=args.threads)
if pipeline == "antspyx": if pipeline == "antspyx":
return RP.register_antspyx(ct, t1c) return RP.register_antspyx(ct, t1c)
if pipeline == "synthmorph": if pipeline == "deepreg":
return RP.register_synthmorph(ct, t1c, model=args.synthmorph_model, return RP.register_deepreg(ct, t1c, workdir=args.deepreg_dir,
gpu=args.synthmorph_gpu, threads=args.threads, checkpoint=args.deepreg_ckpt)
extent=args.synthmorph_extent)
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)
@ -477,15 +476,8 @@ def parse_args():
"production = full multi-start tilt grid (slow)") "production = full multi-start tilt grid (slow)")
ap.add_argument("--out-dir", default=os.path.join(ROOT, "results", "bench_registration")) ap.add_argument("--out-dir", default=os.path.join(ROOT, "results", "bench_registration"))
ap.add_argument("--elastix-dir", default=None, help="elastix scratch dir (default data/m6reg/elastix)") ap.add_argument("--elastix-dir", default=None, help="elastix scratch dir (default data/m6reg/elastix)")
ap.add_argument("--synthmorph-model", choices=["rigid", "affine"], default="rigid", ap.add_argument("--deepreg-dir", default=None, help="deepreg scratch dir")
help="SynthMorph matrix model (default rigid; deformable/joint are " ap.add_argument("--deepreg-ckpt", default=None, help="path to a DeepReg 3D pretrained checkpoint")
"outside the rigid-transform contract)")
ap.add_argument("--synthmorph-gpu", type=int, default=None,
help="GPU device index for SynthMorph (default: auto-spread over "
"the host GPUs by pid; -1 forces CPU)")
ap.add_argument("--synthmorph-extent", type=int, default=256,
help="SynthMorph isotropic registration-space extent in mm "
"(192 or 256; default 256)")
ap.add_argument("--impact-model", default=None, ap.add_argument("--impact-model", default=None,
help="IMPACT TorchScript model, path relative to the HF " help="IMPACT TorchScript model, path relative to the HF "
f"layout (default {RP.IMPACT_DEFAULT_MODEL}) or an absolute host path") f"layout (default {RP.IMPACT_DEFAULT_MODEL}) or an absolute host path")
@ -514,6 +506,7 @@ def main():
merge(args) merge(args)
return return
args.elastix_dir = args.elastix_dir or d("data/m6reg/elastix") args.elastix_dir = args.elastix_dir or d("data/m6reg/elastix")
args.deepreg_dir = args.deepreg_dir or d("data/m6reg/deepreg")
if args.qa_images: if args.qa_images:
args.qa_out = os.path.abspath(args.qa_out or os.path.join(args.out_dir, "qa")) args.qa_out = os.path.abspath(args.qa_out or os.path.join(args.out_dir, "qa"))
os.makedirs(args.qa_out, exist_ok=True) os.makedirs(args.qa_out, exist_ok=True)

413
src/registration_pipelines.py
View file

@ -21,39 +21,23 @@ destination->source convention this transform therefore maps **T1c -> CT**
coordinates. This is exactly the convention the production SimpleITK pipeline coordinates. This is exactly the convention the production SimpleITK pipeline
returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``). returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``).
Pipelines (listed alphabetically) Pipelines
--------------------------------- ---------
All pipelines use the shared convention **fixed = CT, moving = T1c**, except - ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a
IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its CT-space head box, multi-start tilt grid. Reuses the *proven* production
TotalSegmentator metric -- see ``register_impact``). A pipeline that is a ``register_ct_to_t1c`` so the baseline is bit-identical to 02_build_m6_dataset.
local optimizer and cannot recover the gross inter-session offset from - ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` + ``itk-elastix``,
identity (antspyx, impact) starts from a centroid-aligned translation; the elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
others (simpleitk, elastix, synthmorph) do not require it and start from - ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
identity / their own global search. - ``register_deepreg`` : Oxford DeepReg (CycleGAN deformation-field
registration; needs a pretrained checkpoint, see module docstring).
- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). - ``register_impact`` : IMPACT (semantic-similarity metric) inside the
Centroid pre-shift of the moving (T1c) image is *required* (ANTs
Rigid-from-identity collapses on the gross offset).
- ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` +
``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
Identity start (the ASGD + multi-resolution pyramid recovers the offset).
- ``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
GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); centroid init 6-DOF, GPU (Jacobian mode). Needs Docker + the image + a model asset.
kept (local ASGD on the semantic metric benefits from the gross alignment).
Needs Docker + the image + a model asset.
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a
CT-space head box. Reuses the *proven* production ``register_ct_to_t1c`` in
``mode="production"``; the benchmark mode is a single identity-start
multi-resolution pyramid run.
- ``register_synthmorph`` : FreeSurfer SynthMorph (TensorFlow, pretrained
anatomy-aware contrast-invariant model) inside the ``freesurfer/synthmorph``
Docker image. Rigid 6-DOF matrix model, GPU when available. Global deep
method -- no centroid init. Needs Docker + the pulled image.
A toolkit that is not installed (or, for the Docker-based pipelines, missing A toolkit that is not installed (or, for DeepReg, missing its checkpoint, or,
its image) for IMPACT, missing its Docker image / model)
raises :class:`PipelineUnavailable` with the exact install / precondition. The raises :class:`PipelineUnavailable` with the exact install / precondition. The
benchmark catches that and reports the pipeline as *not available* instead of benchmark catches that and reports the pipeline as *not available* instead of
failing the whole run, so the script is env-portable: run it under whichever failing the whole run, so the script is env-portable: run it under whichever
@ -66,10 +50,11 @@ per toolkit. :func:`fix_orientation` resolves this by warping the CT with both
``T`` and ``T^-1`` and keeping whichever gives the smaller skull-ring ``T`` and ``T^-1`` and keeping whichever gives the smaller skull-ring
concentricity gap (a pure correctness gate -- a wrongly-oriented rigid transform concentricity gap (a pure correctness gate -- a wrongly-oriented rigid transform
lands the head on the far side of the box and yields a pathological gap). It is lands the head on the far side of the box and yields a pathological gap). It is
applied to elastix / antspyx / synthmorph / impact, not to the trusted applied to elastix / antspyx / deepreg / impact, not to the trusted production
production SimpleITK path. SimpleITK path.
""" """
import os import os
import json
import shutil import shutil
import subprocess import subprocess
import importlib.util import importlib.util
@ -184,11 +169,10 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
"""SimpleITK rigid CT->T1c, Mattes MI on a CT-space head box. """SimpleITK rigid CT->T1c, Mattes MI on a CT-space 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. centered start, capped at ``iters_full`` iterations at full resolution.
This is the *standard* SimpleITK algorithm and matches the budget the This is the *standard* SimpleITK algorithm and matches the budget the
other toolkits get in the benchmark (one rigid MI registration), so the other toolkits get in the benchmark (one rigid MI registration), so the
pipelines are compared on equal footing. The pyramid recovers the gross pipelines are compared on equal footing.
inter-session offset, so 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
tilt but much slower -- not the default for the benchmark). tilt but much slower -- not the default for the benchmark).
@ -225,11 +209,12 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
ref.SetOrigin(tuple(float(x) for x in lo)) ref.SetOrigin(tuple(float(x) for x in lo))
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
# Identity start: the multi-resolution Mattes-MI pyramid recovers the gross # centered start: align CT head centroid with T1c head centroid (shared init
# inter-session offset on its own, so no centroid init is required here. # with the other single-init toolkits via centroid_offset)
off = centroid_offset(ct, t1c, ct_box=fixed) or (0.0, 0.0, 0.0)
ct_c = np.array([(lo[i] + hi[i]) / 2.0 for i in range(3)]) ct_c = np.array([(lo[i] + hi[i]) / 2.0 for i in range(3)])
init = sitk.Euler3DTransform(); init.SetCenter(tuple(float(x) for x in ct_c)) init = sitk.Euler3DTransform(); init.SetCenter(tuple(float(x) for x in ct_c))
init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) init.SetParameters((0.0, 0.0, 0.0) + off)
import time import time
t0 = time.time() t0 = time.time()
r = sitk.ImageRegistrationMethod() r = sitk.ImageRegistrationMethod()
@ -311,20 +296,21 @@ def _itk_transform_to_sitk_rigid(itk_tf):
def register_elastix(ct, t1c, workdir=None, threads=8): def register_elastix(ct, t1c, workdir=None, threads=8):
"""ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), Mattes MI. """ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI.
Runs in-process through the ITK-Elastix python API (pip ``itk`` + Runs in-process through the ITK-Elastix python API (pip ``itk`` +
``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the ``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the
5.x replacement for ``Rigid3DTransform``) + Mattes MI 5.x replacement for ``Rigid3DTransform``) + Mattes MI
(``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common
CT-space head box, from the identity start (the ASGD + multi-resolution CT-space head box, from the *shared* centroid-aligned start
pyramid recovers the gross inter-session offset, so no centroid init is (``centroid_offset`` as an initial TranslationTransform, same init the
required). Elastix 5.x component names: sampler ``Random`` (not SimpleITK/ANTs benchmark modes use). Elastix 5.x component names: sampler
``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``). An output ``Random`` (not ``RandomSampler``), metric key ``Metric`` (not
directory is required for the C++ filter to run, and ``WriteResultImage`` ``ImageMetric``). An output directory is required so the initial
must be the string ``"false"`` (a Python bool False) or the filter raises transform parameter object is honored (the C++ filter only wires it up
"Could not read result image" after a successful registration. when ``SetOutputDirectory`` is set), and ``WriteResultImage`` must be the
""" string ``"false"`` (a Python bool False) or the filter raises "Could not
read result image" after a successful registration."""
itk = _elastix_api() itk = _elastix_api()
lo, hi = _ct_head_box_mm(ct) lo, hi = _ct_head_box_mm(ct)
if lo is None: if lo is None:
@ -337,8 +323,18 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
ref.SetSpacing(t1c.GetSpacing()) ref.SetSpacing(t1c.GetSpacing())
ref.SetDirection(t1c.GetDirection()) ref.SetDirection(t1c.GetDirection())
ref.SetOrigin(tuple(float(x) for x in lo)) ref.SetOrigin(tuple(float(x) for x in lo))
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = reference fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = moving moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = moving
# shared centered start: elastix transforms map fixed(T1c) -> moving(CT),
# i.e. the inverse direction of centroid_offset, hence the negation.
off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0)
init_tp0 = {
"Transform": "TranslationTransform",
"TransformBase": "ElastixTransformBase",
"TransformParameterReadoutFormat": 1,
"NumberOfParameters": 3,
"TransformParameters": [-float(off[0]), -float(off[1]), -float(off[2])],
}
params = { params = {
"ImageDimension": 3, "ImageDimension": 3,
"FixedInternalImagePixelType": "float", "FixedInternalImagePixelType": "float",
@ -365,6 +361,7 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
reg.SetFixedImage(_sitk_to_itk(fixed, itk)) reg.SetFixedImage(_sitk_to_itk(fixed, itk))
reg.SetMovingImage(_sitk_to_itk(moving, itk)) reg.SetMovingImage(_sitk_to_itk(moving, itk))
reg.SetParameterObject(itk.dict_to_parameter_object(params)) reg.SetParameterObject(itk.dict_to_parameter_object(params))
reg.SetInitialTransformParameterObject(itk.dict_to_parameter_object(init_tp0))
reg.SetOutputDirectory(outdir) reg.SetOutputDirectory(outdir)
reg.SetNumberOfThreads(threads) reg.SetNumberOfThreads(threads)
reg.LogToConsoleOff() reg.LogToConsoleOff()
@ -376,9 +373,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
tail = "\n".join(buf.getvalue().splitlines()[-10:]) tail = "\n".join(buf.getvalue().splitlines()[-10:])
raise RuntimeError(f"elastix registration failed: {e}\n{tail}") from e raise RuntimeError(f"elastix registration failed: {e}\n{tail}") from e
dt = time.time() - t0 dt = time.time() - t0
# The combination transform maps fixed(CT) -> moving(T1c) (elastix # The combination transform maps fixed(T1c) -> moving(CT) == the contract
# convention), i.e. the inverse of the contract (T1c->CT); # (T1c->CT); fix_orientation still verifies it against the skull ring.
# fix_orientation resolves the direction against the skull ring.
cand = _itk_transform_to_sitk_rigid(reg.GetCombinationTransform()) cand = _itk_transform_to_sitk_rigid(reg.GetCombinationTransform())
T, inverted = fix_orientation(ct, t1c, cand) T, inverted = fix_orientation(ct, t1c, cand)
try: try:
@ -389,9 +385,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
meta = { meta = {
"pipeline": "elastix", "toolkit": version, "pipeline": "elastix", "toolkit": version,
"metric": "MattesMutualInformation (head box, CT space)", "metric": "MattesMutualInformation (head box, CT space)",
"transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)", "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, centroid 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": "T1c", "moving": "CT", "seconds_reg": round(dt, 1),
} }
return T, meta return T, meta
@ -403,7 +399,7 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None, def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
aff_shrink_factors=None, aff_smoothing_sigmas=None, aff_shrink_factors=None, aff_smoothing_sigmas=None,
aff_metric="mattes", use_head_box=True, center_init=True): aff_metric="mattes", use_head_box=True, center_init=True):
"""ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c). """ANTs rigid CT->T1c via antspyx (fixed=T1c, moving=CT).
antspyx 0.6.3: ``ants.registration(fixed, moving, type_of_transform='Rigid', antspyx 0.6.3: ``ants.registration(fixed, moving, type_of_transform='Rigid',
aff_metric='mattes', aff_iterations=..., aff_shrink_factors=..., aff_metric='mattes', aff_iterations=..., aff_shrink_factors=...,
@ -411,11 +407,8 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
multi-resolution Mattes-MI rigid registration -- the cross-modality metric multi-resolution Mattes-MI rigid registration -- the cross-modality metric
this CT/MR pair needs. To keep the benchmark tractable and consistent with this CT/MR pair needs. To keep the benchmark tractable and consistent with
the other pipelines (and to avoid neck/air in the metric), both images are the other pipelines (and to avoid neck/air in the metric), both images are
first resampled into the same CT-space head box. ANTs Rigid-from-identity first resampled into the same CT-space head box. Returns the transform that
does NOT recover the gross inter-session offset, so a centroid pre-shift of warps CT into the T1c grid (orientation verified by :func:`fix_orientation`).
the moving (T1c) image is required (see ``center_init``). Returns the
transform that warps CT into the T1c grid (orientation verified by
:func:`fix_orientation`).
""" """
try: try:
import ants # antspyx ships the `ants` module import ants # antspyx ships the `ants` module
@ -437,26 +430,25 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32)) ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection()) ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
ref.SetOrigin(tuple(float(x) for x in lo)) ref.SetOrigin(tuple(float(x) for x in lo))
fx = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (fixed) fx = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (fixed)
mv = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (moving) mv = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (moving)
else: else:
fx, mv = ct, t1c fx, mv = t1c, ct
import time import time
t0 = time.time() t0 = time.time()
# Centered start (REQUIRED for ANTs): antspyx cannot use a # Centered start: antspyx cannot use a TranslationTransform as initial_transform
# TranslationTransform as initial_transform (the ANTs CLI rejects it), and # (the ANTs CLI rejects it), and Rigid-from-identity does not recover the large
# Rigid-from-identity does not recover the large (~100-130 mm) inter-session # (~30-60 mm) inter-session offset. So we pre-shift the CT box so its head
# offset -- without this init ANTs collapses on 8/9 cases. So we pre-shift the # centroid coincides with the T1c head centroid, run Rigid from identity, and
# T1c (moving) box so its head centroid coincides with the CT (fixed) head # compose the shift back into the final transform (undoing it to original CT
# centroid, run Rigid from identity, then compose the shift + inverse back so # coordinates). This gives ANTs the same centered initialization as the other
# the result is expressed in the original T1c/CT coordinates. # single-init toolkits so the pipelines are compared on their optimizer.
off = None off = None
if center_init: if center_init:
off = centroid_offset(ct, t1c, ct_box=fx) off = centroid_offset(ct, t1c, ct_box=mv)
if off is not None: if off is not None:
# shift the T1c content by -off (head c_t1c -> c_ct): resample with +off
tr_off = sitk.TranslationTransform(3) tr_off = sitk.TranslationTransform(3)
tr_off.SetOffset((off[0], off[1], off[2])) tr_off.SetOffset((-off[0], -off[1], -off[2]))
mv = sitk.Resample(mv, mv, tr_off, sitk.sitkLinear, 0.0) mv = sitk.Resample(mv, mv, tr_off, sitk.sitkLinear, 0.0)
af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)) af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32))
am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32))
@ -465,16 +457,16 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
aff_shrink_factors=aff_shrink_factors, aff_shrink_factors=aff_shrink_factors,
aff_smoothing_sigmas=aff_smoothing_sigmas) aff_smoothing_sigmas=aff_smoothing_sigmas)
dt = time.time() - t0 dt = time.time() - t0
# r['fwdtransforms'][0] is the ANTs forward transform: fixed(CT) -> moving # r['fwdtransforms'][0] is a path to the (single) Rigid transform file, in the
# (pre-shifted T1c). We need the contract T1c(orig) -> CT: map T1c(orig) to # shared physical space (but on the pre-shifted CT). Compose the shift back so
# the pre-shifted T1c (translation by -off), then invert the ANTs transform. # the result maps T1c -> original CT; fix_orientation resolves T vs T^-1.
T_ants = sitk.ReadTransform(r["fwdtransforms"][0]) T_ants = sitk.ReadTransform(r["fwdtransforms"][0])
if off is not None: if off is not None:
tr_shift = sitk.TranslationTransform(3) tr_undo = sitk.TranslationTransform(3)
tr_shift.SetOffset((-off[0], -off[1], -off[2])) # T1c orig -> pre tr_undo.SetOffset((-off[0], -off[1], -off[2]))
cand = sitk.CompositeTransform([tr_shift, T_ants.GetInverse()]) cand = sitk.CompositeTransform([T_ants, tr_undo])
else: else:
cand = T_ants.GetInverse() cand = T_ants
T, inverted = fix_orientation(ct, t1c, cand) T, inverted = fix_orientation(ct, t1c, cand)
meta = { meta = {
"pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", "pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
@ -492,143 +484,109 @@ def sys_pyver():
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# 4. SynthMorph (Docker image freesurfer/synthmorph) # 4. DeepReg (Oxford)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
SYNTHMORPH_IMAGE = "freesurfer/synthmorph" def register_deepreg(ct, t1c, workdir=None, checkpoint=None, config=None):
# The image entrypoint is ``mri_synthmorph``; the ``register`` subcommand and """Oxford DeepReg deformation-field registration (CT -> T1c).
# options follow. ``TF_FORCE_GPU_ALLOW_GROWTH`` lets concurrent benchmark
# workers share a GPU instead of each pre-allocating all of it.
DeepReg is not a one-shot rigid tool: it is a CycleGAN-based image-to-image
def _synthmorph_docker_ok(): registration that must be *trained* or run from a *pretrained checkpoint*.
"""Return (ok, note): is the Docker CLI + the synthmorph image available?""" It predicts a voxel displacement field; we convert that to a SimpleITK
if shutil.which("docker") is None: DisplacementFieldTransform to warp CT/GTV onto the T1c grid. Preconditions:
return False, "docker CLI not found on PATH" - ``deepreg`` installed from git (https://github.com/OB-CVHU/DeepReg)
try: -- the PyPI ``deepreg`` (0.1.2) is an unrelated stub that pins
r = subprocess.run(["docker", "image", "inspect", SYNTHMORPH_IMAGE], numpy==1.18.5 and will not work.
capture_output=True, text=True, timeout=30) - a 3D pretrained checkpoint (``checkpoint``) -- typically one trained on
if r.returncode != 0: a CT/MRI template pair; see the repo's pretrained-models note.
return False, (f"docker image {SYNTHMORPH_IMAGE} not pulled "
f"(run: docker pull {SYNTHMORPH_IMAGE})")
except Exception as e: # noqa
return False, f"docker image check failed ({type(e).__name__}: {e})"
return True, SYNTHMORPH_IMAGE
def _lta_to_sitk(path):
"""Parse a FreeSurfer LTA (LINEAR_RAS_TO_RAS) into a SimpleITK AffineTransform.
The 4x4 matrix block lives in the shared RAS world space. Both inputs are
RAS-oriented by ``read_nii_img`` (DICOMOrient), so RAS == SimpleITK's
physical frame and the matrix applies directly to world coordinates."""
lines = [l for l in open(path).read().splitlines()
if l and not l.startswith("#")]
i = lines.index("1 4 4")
M = np.array([[float(x) for x in lines[i + 1 + r].split()] for r in range(4)])
T = sitk.AffineTransform(3)
T.SetMatrix(M[:3, :3].ravel())
T.SetTranslation(M[:3, 3])
return T
def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
extent=256, keep_scratch=False):
"""SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image
``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``,
pretrained ``synthmorph.rigid.1.h5`` weights).
SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware,
acquisition-agnostic deep registration: each volume is normalized to an
isotropic 1-mm LIA network space and the (pretrained) network estimates the
transform, so no skull-stripping or intensity preprocessing is needed. Being
a global deep method it recovers the gross inter-session offset on its own,
so **no centroid init is used**. The pipeline
- runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible
(auto-spread over the host GPUs by pid; ``gpu=-1`` forces CPU),
- reads back the ``-T`` LTA (fixed->moving = the T1c->CT transform the
unified contract requires); :func:`fix_orientation` verifies the
direction against the skull ring.
Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI and
the ``freesurfer/synthmorph`` image. Only the matrix models (``rigid``,
``affine``) are exposed: the deformable/joint models emit displacement
fields, outside the rigid-transform contract.
""" """
import time
ok, note = _synthmorph_docker_ok()
if not ok:
raise PipelineUnavailable(f"synthmorph: {note}")
if model not in ("rigid", "affine"):
raise PipelineUnavailable(
f"synthmorph: only the matrix models (rigid/affine) fit the unified "
f"rigid-transform contract, got {model!r}")
# GPU: auto-spread across the host GPUs by pid, or forced value
if gpu is None:
ngpu = _n_gpus()
gpu = (os.getpid() % ngpu) if ngpu else -1
# per-case scratch on a local (non-NFS) tmp fs
tmp = tempfile.mkdtemp(prefix="synthmorph_reg_")
try: try:
# fixed = CT, moving = T1c (shared convention for all pipelines) import deepreg # noqa
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32), except Exception as e: # noqa
os.path.join(tmp, "fixed.nii.gz"), True) # CT raise PipelineUnavailable(
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32), f"deepreg not importable ({e!r}). Install the Oxford DeepReg from "
os.path.join(tmp, "moving.nii.gz"), True) # T1c "source: `git clone https://github.com/OB-CVHU/DeepReg && "
"cd DeepReg && pip install -e .` (NOT the PyPI `deepreg` stub). "
# run as the host uid/gid so the container-written LTA is host-owned "Also requires a pretrained 3D checkpoint (see --deepreg-ckpt).")
# (clean scratch removal) and TF still reaches the GPU. if not checkpoint or not os.path.exists(checkpoint):
cmd = ["docker", "run", "--rm", raise PipelineUnavailable(
"--user", f"{os.getuid()}:{os.getgid()}", "deepreg: no pretrained checkpoint provided. Pass --deepreg-ckpt "
"-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"] "to a 3D .ckpt produced by `deepreg train` (or a released "
if gpu >= 0: "pretrained model). DeepReg cannot do one-shot registration "
cmd += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"] "without weights.")
cmd += ["-v", f"{tmp}:/data", SYNTHMORPH_IMAGE, workdir = workdir or d("data/m6reg/deepreg")
"register", "-m", model] workdir = os.path.join(workdir, f"pid{os.getpid()}") # per-worker scratch
if gpu >= 0: os.makedirs(workdir, exist_ok=True)
cmd.append("-g") cfg = _deepreg_config(workdir, ct, t1c, checkpoint)
if threads: outdir = os.path.join(workdir, "pred")
cmd += ["-j", str(threads)] subprocess.run(
cmd += ["-e", str(extent), ["python", "-m", "deepreg.main", "test",
"-T", "/data/f2m.lta", "-checkpoint_path", checkpoint, "-config_path", cfg, "-out_dir", outdir],
"/data/moving.nii.gz", "/data/fixed.nii.gz"] check=True, capture_output=True, text=True, cwd=workdir)
t0 = time.time() # predicted deformation field (voxels) on the fixed (T1c) grid
proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200) pred = _find_single_nii(os.path.join(outdir))
dt = time.time() - t0 if pred is None:
if proc.returncode != 0: raise PipelineUnavailable("deepreg: no predicted deformation field found")
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) disp = sitk.GetArrayFromImage(sitk.ReadImage(pred)) # (Z, Y, X, 3)
raise RuntimeError(f"synthmorph: docker failed (rc={proc.returncode}):\n{tail}") T = _displacement_to_sitk(disp, t1c)
if not os.path.exists(os.path.join(tmp, "f2m.lta")):
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
raise RuntimeError(f"synthmorph: no /data/f2m.lta produced:\n{tail}")
T = _lta_to_sitk(os.path.join(tmp, "f2m.lta"))
T, inverted = fix_orientation(ct, t1c, T)
finally:
if not keep_scratch:
shutil.rmtree(tmp, ignore_errors=True)
meta = { meta = {
"pipeline": "synthmorph", "pipeline": "deepreg", "toolkit": "Oxford DeepReg (CycleGAN)",
"toolkit": f"Docker {SYNTHMORPH_IMAGE} (TensorFlow, pretrained)", "metric": "CycleGAN + registration loss (pretrained)",
"metric": "SynthMorph anatomy-aware contrast-invariant (deep, pretrained)", "transform": "3D displacement field", "convention": "T1c->CT (verified)",
"transform": (f"SynthMorph {model} (identity start; " "auto_inverted": False, "checkpoint": checkpoint,
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, "
f"extent {extent})"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "CT", "moving": "T1c",
"seconds_reg": round(dt, 1),
"docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu,
"extent": extent,
} }
if keep_scratch:
meta["scratch"] = tmp
return T, meta return T, meta
def _deepreg_config(workdir, ct, t1c, checkpoint):
"""Minimal DeepReg test config: moving=CT, fixed=T1c, 3D, single channel."""
from deepreg.registry import registry # noqa
cfg = {
"general": {"data_format": "nifti", "device": "cuda:0"},
"data": {"moving": os.path.join(workdir, "moving"),
"fixed": os.path.join(workdir, "fixed"),
"field": "deformation", "moving_fields": "t1c", "fixed_fields": "t1c",
"mask": None, "resample_size": "30 186 30"},
"validation": {},
"training": {"epoch": 1, "batch_size": 1,
"moving_fields": "t1c", "fixed_fields": "t1c",
"resample_size": "30 186 30", "resample_spacing": "1 1 1"},
"testing": {"moving": os.path.join(workdir, "moving"),
"fixed": os.path.join(workdir, "fixed"),
"resample_size": "30 186 30", "resample_spacing": "1 1 1"},
}
p = os.path.join(workdir, "deepreg_test.yaml")
try:
import yaml
with open(p, "w") as f:
yaml.safe_dump(cfg, f)
except Exception:
with open(p, "w") as f:
json.dump(cfg, f)
return p
def _find_single_nii(dirpath):
if not os.path.isdir(dirpath):
return None
ns = [os.path.join(dirpath, x) for x in sorted(os.listdir(dirpath))
if x.endswith(".nii.gz") or x.endswith(".nii")]
return ns[0] if ns else None
def _displacement_to_sitk(disp, ref):
"""Build a sitk.DisplacementFieldTransform from a (Z,Y,X,3) voxel field."""
from scipy.ndimage import map_coordinates
zz, yy, xx = np.mgrid[0:disp.shape[0], 0:disp.shape[1], 0:disp.shape[2]]
# disp[i,j,k] = (dz, dy, dx) to add at grid point (i,j,k) to reach the source
src_z = zz + disp[..., 0]; src_y = yy + disp[..., 1]; src_x = xx + disp[..., 2]
# DisplacementFieldTransform stores, per output voxel, the offset to the input
off = sitk.GetImageFromArray(disp[..., ::-1].astype(np.float32)) # (X,Y,Z,3)
off.CopyInformation(ref)
t = sitk.DisplacementFieldTransform(off)
return t
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# 5. IMPACT (Docker image vboussot/elastix_impact) # 5. IMPACT (Docker image vboussot/elastix_impact)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
@ -660,7 +618,7 @@ def _impact_docker_ok():
_NGPUS_CACHE = None _NGPUS_CACHE = None
def _n_gpus(): def _impact_n_gpus():
global _NGPUS_CACHE global _NGPUS_CACHE
if _NGPUS_CACHE is None: if _NGPUS_CACHE is None:
n = 0 n = 0
@ -794,13 +752,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
"""IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image
``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric). ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
NOTE on the fixed/moving convention: every other pipeline here uses
fixed=CT, moving=T1c. IMPACT is the exception -- swapping its fixed/moving
to CT/T1c breaks the TotalSegmentator-based semantic metric (elastix aborts
with "Error in metric" or converges to a wrong rotation), so IMPACT keeps
the working MRI-as-fixed (T1c) / CT-as-moving convention. The returned
transform still satisfies the unified contract (T1c->CT).
IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic* IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run
on randomly sampled patches of the fixed and moving images and the on randomly sampled patches of the fixed and moving images and the
@ -851,15 +802,14 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference
moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = moving moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = moving
# shared centered start (REQUIRED: the local ASGD optimizer does not recover # shared centered start: elastix maps fixed(T1c)->moving(CT), i.e. the
# the gross offset from identity). Elastix maps fixed(CT)->moving(T1c), the # inverse of centroid_offset, hence the negation.
# same direction as centroid_offset, so no negation.
off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0) off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0)
init = (-float(off[0]), -float(off[1]), -float(off[2])) init = (-float(off[0]), -float(off[1]), -float(off[2]))
# 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:
ngpu = _n_gpus() ngpu = _impact_n_gpus()
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
@ -935,8 +885,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
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 (box/world frame) composed with the fixed -t0 # the optimized Euler (box/world frame) composed with the fixed -t0
# translation: T = T_euler(T_init(p)). Both map fixed(CT)->moving(T1c), # translation: T = T_euler(T_init(p))
# i.e. the inverse of the contract (T1c->CT); fix_orientation resolves it.
T_main = sitk.ReadTransform(itk_tf) T_main = sitk.ReadTransform(itk_tf)
T_init = sitk.TranslationTransform(3) T_init = sitk.TranslationTransform(3)
T_init.SetOffset(init) T_init.SetOffset(init)
@ -964,11 +913,11 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
REGISTERS = { REGISTERS = {
"antspyx": register_antspyx,
"elastix": register_elastix,
"impact": register_impact,
"simpleitk": register_simpleitk, "simpleitk": register_simpleitk,
"synthmorph": register_synthmorph, "elastix": register_elastix,
"antspyx": register_antspyx,
"deepreg": register_deepreg,
"impact": register_impact,
} }
ALL_PIPELINES = list(REGISTERS) ALL_PIPELINES = list(REGISTERS)
@ -994,12 +943,12 @@ 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 == "synthmorph": elif name == "deepreg":
ok, note = _synthmorph_docker_ok() try:
if not ok: import deepreg # noqa
out[name] = (False, note) out[name] = (True, "deepreg importable")
else: except Exception as e: # noqa
out[name] = (True, f"Docker {SYNTHMORPH_IMAGE}") out[name] = (False, "deepreg not importable")
elif name == "impact": elif name == "impact":
ok, note = _impact_docker_ok() ok, note = _impact_docker_ok()
if not ok: if not ok: