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

Compare commits

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

2 commits
6b85e950e5 ... 4ad972b00f

Author SHA1 Message Date
Furen Xiao
4ad972b00f refactor(pipelines): remove shared resample helper and clarify init conventions 
Drop the _ct_head_box_resample utility that was only needed for centered
initialization, since centroid alignment is now scoped to local optimizers
(antspyx, impact) rather than all pipelines. Reorder the pipeline table and
module docstring alphabetically and add explicit fixed/moving and centroid-init
columns so the per-toolkit registration convention is unambiguous.
 2026-09-29 04:09:55 +08:00
Furen Xiao
f80f22dabe feat(pipelines): add SynthMorph Docker registration and drop DeepReg 
Replace the unavailable DeepReg benchmark path with the freesurfer/synthmorph
Docker pipeline for rigid CT-to-T1c registration. Add a shared CT head-box
resample helper for centered initialization and update benchmark docs and CLI
for the new SynthMorph options.
 2026-09-28 23:13:44 +08:00
3 changed files with 461 additions and 316 deletions
Show stats Expand all files Collapse all files

307
REGISTRATION_BENCHMARK.md
View file

@ -22,21 +22,34 @@ 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 register destination→source convention this `T` maps **T1c → CT**. All pipelines use
on the same CT-space head box (MI computed on head tissue, not neck/air) and, **fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c,
for the single-init tools, from the **same centroid-aligned start** moving = CT** (swapping its roles breaks the TotalSegmentator-based metric,
(`centroid_offset`), so the comparison is on the optimizer / metric / transform which aborts with "Error in metric" or converges to a wrong rotation). The
handling — not on who finds the ~30–60 mm gross offset first. intensity-metric pipelines register on the same CT-space head box (MI computed
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 four pipelines ## The five pipelines
| # | pipeline | toolkit | metric | transform | status in this env | Listed alphabetically. All use **fixed = CT, moving = T1c** except `impact`,
|---|----------|---------|--------|-----------|--------------------| which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator
| 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) | metric). A pipeline that is a local optimizer and cannot recover the gross
| 2 | `simpleitk` | SimpleITK 2.5.6 | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | inter-session offset from identity (`antspyx`, `impact`) starts from a
| 3 | `antspyx` | ANTsPy / ANTs 0.6.3 | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | centroid-aligned translation; the others start from identity / their own global
| 4 | `deepreg` | Oxford DeepReg (CycleGAN) | CycleGAN + reg. loss (pretrained) | 3D displacement field | **not installed** (see below) | search.
| 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
@ -85,75 +98,94 @@ 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) |
|---|---|---|---|---|---|---|---|---| |---|---|---|---|---|---|---|---|---|
| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.6 | 0.93 | **0.643 / 0.630** | | **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 18.9 / 20.1 | 31.7 / 31.3 | 10.9 | 0.94 | **0.436 / 0.393** |
| **antspyx** | ANTsPy 0.6.3 | 9 | **0/10 (0.00)** | 14.9 / 16.5 | 18.4 / 22.7 | 14.5 | 0.93 | **0.598 / 0.529** | | **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 16.6 / 17.6 | 16.9 / 18.8 | 2.4 | 0.93 | **0.641 / 0.583** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | **5/10 (0.50)** | 15.3 / 16.1 | 9.9 / 13.8 | 3.2 | 0.93 | **0.562 / 0.549** | | **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 13.1 / 16.8 | 13.2 / 18.6 | **576.0** | 0.93 | **0.469 / 0.482** |
| **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** | | **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** |
| deepreg | — | — | — | — | — | — | — | — | | **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** |
*(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.4 accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.6
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 | simpleitk | antspyx | elastix | impact | | case | antspyx | elastix | impact | simpleitk | synthmorph |
|---|---|---|---|---| |---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 24.4 / 27.7 | 26.4 / 30.3 | **8.4 / 12.6** | | m6_22IRRKU6_20240115 | 18.9 / 43.8 | 27.0 / 29.3 | 8.3 / 13.2 | 26.2 / 29.6 | 10.4 / 23.7 |
| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.4 / 24.8 | 21.1 / 23.3 | 21.4 / 24.1 | | m6_23XIMO4M_20230222 | 22.9 / 24.5 | 22.6 / 23.9 | 23.1 / 25.7 | 22.6 / 24.0 | 20.3 / 23.3 |
| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.5 / 18.4 | 10.6 / 8.2 **OK** | 10.0 / 8.3 **OK** | | 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_20231017 | 11.5 / 7.3 **OK** | 14.9 / 15.5 | 11.0 / 6.1 **OK** | 11.1 / 8.8 **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_20240820 | 11.6 / 9.6 **OK** | 14.3 / 17.5 | 15.3 / 14.0 | 54.4 / 15.1 | | 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_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 17.7 / 37.6 | 11.4 / 3.2 **OK** | 23.9 / 22.6 | | 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_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** | 10.1 / 5.9 **OK** | | 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_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.1 / 14.3 | 12.0 / 6.8 **OK** | 12.4 / 22.3 | | m6_2C4TKI6H_20240506 | 15.3 / 14.2 | 14.2 / 15.3 | 13.1 / 17.4 | 12.7 / 13.5 | 16.7 / 28.2 |
| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 17.2 / 35.4 | 18.9 / 22.7 | 19.6 / 30.3 | | m6_2CDNN7BH_20240102 | 19.4 / 31.7 | 19.5 / 30.3 | 19.6 / 30.9 | 19.3 / 30.1 | 18.9 / 16.8 |
### Interpretation ### Interpretation
- **Translation is recovered well by all four** (mid gaps mostly 6–20 mm, - **Translation is recovered well by all five** (mid gaps mostly 6–20 mm,
within the `≤ 20 mm` band). `tv-in-CT` is identical (~0.93): the warped GTV within the `≤ 20 mm` band, except the grossly-drifted cases below).
lands inside the CT head for every pipeline, so the target is in the right `tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for
place. every pipeline, so the target is in the right 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 centered (the ~30° inter-session tilt produces local minima that a single start
start can lock onto). This is exactly why production `02_build_m6_dataset.py` (identity or centroid) can lock onto). This is exactly why production
uses a **multi-start tilt grid**. `02_build_m6_dataset.py` uses a **multi-start tilt grid**.
- **Elastix is the most accurate of the standard single-init configs** - **SimpleITK is the most reliable intensity-metric tool in this draw**
(5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 18.4) and (4/10 accepted, best median mid 12.7), and **Elastix is the fastest**
the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s (2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best
SimpleITK / 14.5 s ANTs): ASGD's adaptive step sizes climb out of small median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c
misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223` convention the two settle into different local minima on the hard tilted
(std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017`
(std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both (std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles
others accept — different optimizers settle in different local minima on the better, and vice-versa on the drift cases. The fixed/moving direction
hard tilted cases. materially affects where the ASGD/RSGD optimizers converge on the tilted cases.
- **ANTs' Rigid is the weakest of the intensity-metric tools** (0/10 in this - **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this
draw; it is stochastic, so a prior draw gave 1/10): it frequently locks onto draw; it is stochastic, so the count varies): it frequently locks onto the
the tilted local minimum. Note `m6_26TYQI3C`: ANTs nailed the position (mid tilted local minimum (median std 31.7, worst of all). It accepts only
17.7) but missed the orientation (std 37.6) — a textbook tilted local minimum. `m6_26TYQI3C` (mid 11.2 / std 3.7) and misses `m6_24V62QQM_20230223`
(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, std 28–30 mm) it gets **mid 8.4 / std 12.6** — the best intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of
translation of any pipeline and within 0.6 mm of the acceptance threshold — any pipeline and within ~1 mm of the acceptance threshold — because matching
because matching *anatomical semantics* (TotalSegmentator features) is far *anatomical semantics* (TotalSegmentator features) is far more discriminating
more discriminating for gross pose than raw cross-modality intensity. But the for gross pose than raw cross-modality intensity. But the semantic landscape
semantic landscape is rugged: single-init ASGD lands in bad local minima on is rugged: single-init ASGD lands in bad local minima on several cases, and
4 of 9 cases (std 15–30 mm), and one (`m6_24V62QQM_20240820`) even drifts the two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0).
translation (mid 54 mm). Net 3/10 accepted — between ANTs and SimpleITK — at Net 3/10 accepted — between ANTs and SimpleITK — at ~9.6 min/case, two orders
~9.4 min/case, two orders of magnitude slower than the intensity tools. of magnitude slower than the intensity tools.
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~3.2 s - **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.4 s
(warm workers; a worker's first case pays a ~10 s warm-up), SimpleITK ~8.6 s, (warm workers), SimpleITK ~7.6 s, ANTs ~10.9 s (CPU-bound),
ANTs ~14.5 s (CPU-bound), **IMPACT ~562 s** (GPU-bound: every ASGD iteration **SynthMorph ~26 s** (GPU: a single pretrained-model forward pass per case),
runs the TotalSegmentator extractor on the sampled patches with backprop). **IMPACT ~576 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator
extractor on the sampled patches with backprop).
- **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**: - **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
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 - **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
@ -167,14 +199,15 @@ Per-case (mid / std, `OK` = accepted):
## Reproducing ## Reproducing
```bash ```bash
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here); deepreg # main env — SimpleITK + Elastix + ANTsPy (all pip-installed here)
# 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 --n-cases 10 \ python scripts/bench_registration.py --pipelines simpleitk elastix antspyx \
--workers 4 --threads 24 --simpleitk-mode benchmark --n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
# IMPACT (needs Docker + the pulled image + a model; see the Impact section # Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
# below for the one-time setup) # sections 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
@ -185,7 +218,8 @@ 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}`,
`--deepreg-ckpt <.ckpt>` (DeepReg), `--impact-model <rel|abs>`, `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`,
`--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).
@ -216,35 +250,61 @@ 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.)
## Enabling the blocked pipeline ## SynthMorph pipeline (Docker) — one-time setup
**DeepReg (Oxford)** — the PyPI `deepreg` (0.1.2) is an unrelated stub that pins The `synthmorph` pipeline runs the [`freesurfer/synthmorph`](https://hub.docker.com/r/freesurfer/synthmorph)
`numpy==1.18.5`; do not install it. DeepReg is not a one-shot rigid tool: it is a image: [SynthMorph](https://synthmorph.io) (Hoffmann et al., IEEE TMI 2022), an
CycleGAN image-to-image registration that must be **trained or run from a anatomy-aware, acquisition-agnostic **deep** registration (TensorFlow
pretrained checkpoint**. To benchmark it: `VxmAffineFeatureDetector` with the pretrained `synthmorph.rigid.1.h5` weights).
No skull-stripping or intensity preprocessing is needed: each volume is
normalized to an isotropic 1-mm LIA network space internally. It is the
second Docker-based pipeline and the only deep-learning one.
Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command):
```bash ```bash
git clone https://github.com/OB-CVHU/DeepReg && cd DeepReg && pip install -e . # 1. Docker (with NVIDIA container toolkit so --gpus works)
python scripts/bench_registration.py --pipelines deepreg --n-cases 10 \ docker --version
--deepreg-ckpt /path/to/pretrained_3d.ckpt
# 2. Pull the image (~11 GB, one-time)
docker pull freesurfer/synthmorph
``` ```
The adapter runs `deepreg test` to predict a voxel displacement field, converts No extra model assets: the rigid weights ship inside the image
it to a SimpleITK `DisplacementFieldTransform`, and scores it like the others. (`/freesurfer/models`).
### How the adapter drives the image
1. **Rigid matrix model, no init** — `register -m rigid` on GPU when one is
visible (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces
CPU). SynthMorph is a global deep method, so it recovers the gross
inter-session offset on its own and needs **no** initial transform.
2. **Read-back** — the `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in shared world
space; with fixed=CT/moving=T1c the `-T` matrix is the T1c→CT contract
transform) is parsed into a SimpleITK `AffineTransform` and run through the
usual `fix_orientation` gate.
3. **GPU sharing** — the container is run as the host uid (clean scratch
removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
concurrent benchmark workers can share a GPU.
## 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 the only pipeline that TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
is GPU-based and the only one that is Docker-based. Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
**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):
@ -295,17 +355,38 @@ to add the shared centered start:
## Notes / caveats ## Notes / caveats
- The benchmark compares **standard single-init** configurations (centered - The benchmark compares each toolkit in its **native configuration**: the
start). Production uses `--simpleitk-mode production`'s multi-start tilt grid, local optimizers that cannot recover the gross offset from identity
which raises SimpleITK's acceptance at ~20× the cost; that is out of scope for (`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`,
the head-to-head but is the config the dataset actually ships with. `elastix` and `synthmorph` start from identity / their own global search.
- `deepreg` could not be executed here, so its row is "not installed", not Production uses `--simpleitk-mode production`'s multi-start tilt grid, which
measured; its code path is implemented and unit-verified for import/convention. raises SimpleITK's acceptance at ~20× the cost; that is out of scope for the
- `fix_orientation` (elastix/ants/deepreg/impact) chooses between a transform head-to-head but is the config the dataset actually ships with.
- `synthmorph` is deterministic (one fixed pretrained model, no random
subsampling), so its per-case numbers are stable run to run — unlike ANTs /
IMPACT. Its low 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. In the result independent of each toolkit's transform-direction convention. With
practice it never fired (all `auto_inverted: False`): elastix/ants/impact the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the
returned the T1c→CT direction as expected. inverse of the contract, so it reports `auto_inverted: True` (the gate
corrects it); antspyx/synthmorph/impact return the contract direction
directly (`auto_inverted: False`).
- **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
@ -313,30 +394,36 @@ to add the shared centered start:
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.
- The initial transform (the shared `centroid_offset` start) is passed as an - No initial transform is passed (identity start): the ASGD + 3-resolution
`InitialTransformParameterObject` (a `TranslationTransform` with Mattes-MI pyramid recovers the gross inter-session offset on its own, so
`NumberOfParameters 3`); the C++ filter only wires it up when an the shared `centroid_offset` start is not required. An **output directory**
**output directory is set** (`SetOutputDirectory`), so the adapter uses a is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter
scratch dir under `data/m6reg/elastix/itkelastix/`. uses a 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.
- The result transform is recovered by point-mapping - With fixed=CT/moving=T1c the combination transform is CT→T1c (the inverse
of the contract), so `fix_orientation` reports `auto_inverted: True` and
inverts it. The result is recovered by point-mapping
(`GetCombinationTransform().TransformPoint`) into an affine — independent (`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 ~1.8 s/case at full budget (ASGD 500/500/250 - Registration wall time is ~2.4 s/case median at full budget (ASGD
over 3 resolutions); each worker's first case pays a ~10 s warm-up, so the 500/500/250 over 3 resolutions).
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 min/case (vs 3–15 s for the intensity 500-iter / 2000-sample budget is ~9.6 min/case (vs 2–11 s for the other
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
follow the shared fixed=CT/moving=T1c convention: swapping its roles
breaks the TotalSegmentator metric (elastix aborts with "Error in metric"
or converges to a wrong rotation). The returned transform still satisfies
the contract (T1c→CT).
- **Rugged single-init convergence.** The semantic landscape has more local - **Rugged single-init convergence.** The semantic landscape has more local
minima than Mattes-MI: the pipeline nails the hardest repositioned case minima than Mattes-MI: the pipeline nails the hardest repositioned case
(`m6_22IRRKU6`, mid 8.4 mm — best of all) but lands in bad minima on some (`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 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. 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

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 | deepreg | impact simpleitk | elastix | antspyx | synthmorph | 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; PyPI deepreg is a broken stub), the benchmark is REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under
env-portable: run it under whichever conda env has a given toolkit installed, whichever conda env has a given toolkit installed, and it records only the
and it records only the pipelines that are importable there. elastix runs in the main env through the ITK-Elastix pipelines that are importable there. elastix runs in the main env through the
python API (pip ``itk`` + ``itk-elastix``). ``impact`` runs in the ITK-Elastix python API (pip ``itk`` + ``itk-elastix``). ``synthmorph`` and
``vboussot/elastix_impact`` Docker image (needs Docker + the image + a model, ``impact`` run in Docker images (``freesurfer/synthmorph`` and
see REGISTRATION_BENCHMARK.md). The `--merge` mode then combines ``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see
the per-pipeline JSON files (accumulated across envs) into one summary + plot. REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the
per-pipeline JSON files (accumulated across envs) into one summary + plot.
Usage (repo root): Usage (repo root):
# main env (SimpleITK + elastix + antspyx; impact needs Docker + the image; # main env (SimpleITK + elastix + antspyx)
# 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
# impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md # synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
python scripts/bench_registration.py --pipelines impact python scripts/bench_registration.py --pipelines synthmorph 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,9 +225,10 @@ 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 == "deepreg": if pipeline == "synthmorph":
return RP.register_deepreg(ct, t1c, workdir=args.deepreg_dir, return RP.register_synthmorph(ct, t1c, model=args.synthmorph_model,
checkpoint=args.deepreg_ckpt) gpu=args.synthmorph_gpu, threads=args.threads,
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)
@ -476,8 +477,15 @@ 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("--deepreg-dir", default=None, help="deepreg scratch dir") ap.add_argument("--synthmorph-model", choices=["rigid", "affine"], default="rigid",
ap.add_argument("--deepreg-ckpt", default=None, help="path to a DeepReg 3D pretrained checkpoint") help="SynthMorph matrix model (default rigid; deformable/joint are "
"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")
@ -506,7 +514,6 @@ 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)

401
src/registration_pipelines.py
View file

@ -21,23 +21,39 @@ 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 Pipelines (listed alphabetically)
--------- ---------------------------------
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a All pipelines use the shared convention **fixed = CT, moving = T1c**, except
CT-space head box, multi-start tilt grid. Reuses the *proven* production IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its
``register_ct_to_t1c`` so the baseline is bit-identical to 02_build_m6_dataset. TotalSegmentator metric -- see ``register_impact``). A pipeline that is a
- ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` + ``itk-elastix``, local optimizer and cannot recover the gross inter-session offset from
elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD). identity (antspyx, impact) starts from a centroid-aligned translation; the
others (simpleitk, elastix, synthmorph) do not require it and start from
identity / their own global search.
- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). - ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
- ``register_deepreg`` : Oxford DeepReg (CycleGAN deformation-field Centroid pre-shift of the moving (T1c) image is *required* (ANTs
registration; needs a pretrained checkpoint, see module docstring). 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 - ``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 TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF,
6-DOF, GPU (Jacobian mode). Needs Docker + the image + a model asset. GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); centroid init
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 DeepReg, missing its checkpoint, or, A toolkit that is not installed (or, for the Docker-based pipelines, missing
for IMPACT, missing its Docker image / model) its image)
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
@ -50,11 +66,10 @@ 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 / deepreg / impact, not to the trusted production applied to elastix / antspyx / synthmorph / impact, not to the trusted
SimpleITK path. production SimpleITK path.
""" """
import os import os
import json
import shutil import shutil
import subprocess import subprocess
import importlib.util import importlib.util
@ -169,10 +184,11 @@ 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
centered start, capped at ``iters_full`` iterations at full resolution. identity 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. pipelines are compared on equal footing. The pyramid recovers the gross
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).
@ -209,12 +225,11 @@ 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)
# centered start: align CT head centroid with T1c head centroid (shared init # Identity start: the multi-resolution Mattes-MI pyramid recovers the gross
# with the other single-init toolkits via centroid_offset) # inter-session offset on its own, so no centroid init is required here.
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) + off) init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
import time import time
t0 = time.time() t0 = time.time()
r = sitk.ImageRegistrationMethod() r = sitk.ImageRegistrationMethod()
@ -296,21 +311,20 @@ 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=T1c, moving=CT), Mattes MI. """ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), 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 *shared* centroid-aligned start CT-space head box, from the identity start (the ASGD + multi-resolution
(``centroid_offset`` as an initial TranslationTransform, same init the pyramid recovers the gross inter-session offset, so no centroid init is
SimpleITK/ANTs benchmark modes use). Elastix 5.x component names: sampler required). Elastix 5.x component names: sampler ``Random`` (not
``Random`` (not ``RandomSampler``), metric key ``Metric`` (not ``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``). An output
``ImageMetric``). An output directory is required so the initial directory is required for the C++ filter to run, and ``WriteResultImage``
transform parameter object is honored (the C++ filter only wires it up must be the string ``"false"`` (a Python bool False) or the filter raises
when ``SetOutputDirectory`` is set), and ``WriteResultImage`` must be the "Could not read result image" after a successful registration.
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:
@ -323,18 +337,8 @@ 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(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = reference
moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = moving moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = 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",
@ -361,7 +365,6 @@ 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()
@ -373,8 +376,9 @@ 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(T1c) -> moving(CT) == the contract # The combination transform maps fixed(CT) -> moving(T1c) (elastix
# (T1c->CT); fix_orientation still verifies it against the skull ring. # convention), i.e. the inverse of the contract (T1c->CT);
# 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:
@ -385,9 +389,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, centroid init)", "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)",
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "T1c", "moving": "CT", "seconds_reg": round(dt, 1), "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
} }
return T, meta return T, meta
@ -399,7 +403,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=T1c, moving=CT). """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c).
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=...,
@ -407,8 +411,11 @@ 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. Returns the transform that first resampled into the same CT-space head box. ANTs Rigid-from-identity
warps CT into the T1c grid (orientation verified by :func:`fix_orientation`). does NOT recover the gross inter-session offset, so a centroid pre-shift of
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
@ -430,25 +437,26 @@ 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(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (fixed) fx = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (fixed)
mv = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (moving) mv = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (moving)
else: else:
fx, mv = t1c, ct fx, mv = ct, t1c
import time import time
t0 = time.time() t0 = time.time()
# Centered start: antspyx cannot use a TranslationTransform as initial_transform # Centered start (REQUIRED for ANTs): antspyx cannot use a
# (the ANTs CLI rejects it), and Rigid-from-identity does not recover the large # TranslationTransform as initial_transform (the ANTs CLI rejects it), and
# (~30-60 mm) inter-session offset. So we pre-shift the CT box so its head # Rigid-from-identity does not recover the large (~100-130 mm) inter-session
# centroid coincides with the T1c head centroid, run Rigid from identity, and # offset -- without this init ANTs collapses on 8/9 cases. So we pre-shift the
# compose the shift back into the final transform (undoing it to original CT # T1c (moving) box so its head centroid coincides with the CT (fixed) head
# coordinates). This gives ANTs the same centered initialization as the other # centroid, run Rigid from identity, then compose the shift + inverse back so
# single-init toolkits so the pipelines are compared on their optimizer. # the result is expressed in the original T1c/CT coordinates.
off = None off = None
if center_init: if center_init:
off = centroid_offset(ct, t1c, ct_box=mv) off = centroid_offset(ct, t1c, ct_box=fx)
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))
@ -457,16 +465,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 a path to the (single) Rigid transform file, in the # r['fwdtransforms'][0] is the ANTs forward transform: fixed(CT) -> moving
# shared physical space (but on the pre-shifted CT). Compose the shift back so # (pre-shifted T1c). We need the contract T1c(orig) -> CT: map T1c(orig) to
# the result maps T1c -> original CT; fix_orientation resolves T vs T^-1. # the pre-shifted T1c (translation by -off), then invert the ANTs transform.
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_undo = sitk.TranslationTransform(3) tr_shift = sitk.TranslationTransform(3)
tr_undo.SetOffset((-off[0], -off[1], -off[2])) tr_shift.SetOffset((-off[0], -off[1], -off[2])) # T1c orig -> pre
cand = sitk.CompositeTransform([T_ants, tr_undo]) cand = sitk.CompositeTransform([tr_shift, T_ants.GetInverse()])
else: else:
cand = T_ants cand = T_ants.GetInverse()
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__', '?')}",
@ -484,107 +492,141 @@ def sys_pyver():
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# 4. DeepReg (Oxford) # 4. SynthMorph (Docker image freesurfer/synthmorph)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
def register_deepreg(ct, t1c, workdir=None, checkpoint=None, config=None): SYNTHMORPH_IMAGE = "freesurfer/synthmorph"
"""Oxford DeepReg deformation-field registration (CT -> T1c). # The image entrypoint is ``mri_synthmorph``; the ``register`` subcommand and
# 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
registration that must be *trained* or run from a *pretrained checkpoint*. def _synthmorph_docker_ok():
It predicts a voxel displacement field; we convert that to a SimpleITK """Return (ok, note): is the Docker CLI + the synthmorph image available?"""
DisplacementFieldTransform to warp CT/GTV onto the T1c grid. Preconditions: if shutil.which("docker") is None:
- ``deepreg`` installed from git (https://github.com/OB-CVHU/DeepReg) return False, "docker CLI not found on PATH"
-- the PyPI ``deepreg`` (0.1.2) is an unrelated stub that pins
numpy==1.18.5 and will not work.
- a 3D pretrained checkpoint (``checkpoint``) -- typically one trained on
a CT/MRI template pair; see the repo's pretrained-models note.
"""
try: try:
import deepreg # noqa r = subprocess.run(["docker", "image", "inspect", SYNTHMORPH_IMAGE],
capture_output=True, text=True, timeout=30)
if r.returncode != 0:
return False, (f"docker image {SYNTHMORPH_IMAGE} not pulled "
f"(run: docker pull {SYNTHMORPH_IMAGE})")
except Exception as e: # noqa except Exception as e: # noqa
raise PipelineUnavailable( return False, f"docker image check failed ({type(e).__name__}: {e})"
f"deepreg not importable ({e!r}). Install the Oxford DeepReg from " return True, SYNTHMORPH_IMAGE
"source: `git clone https://github.com/OB-CVHU/DeepReg && "
"cd DeepReg && pip install -e .` (NOT the PyPI `deepreg` stub). "
"Also requires a pretrained 3D checkpoint (see --deepreg-ckpt).")
if not checkpoint or not os.path.exists(checkpoint):
raise PipelineUnavailable(
"deepreg: no pretrained checkpoint provided. Pass --deepreg-ckpt "
"to a 3D .ckpt produced by `deepreg train` (or a released "
"pretrained model). DeepReg cannot do one-shot registration "
"without weights.")
workdir = workdir or d("data/m6reg/deepreg")
workdir = os.path.join(workdir, f"pid{os.getpid()}") # per-worker scratch
os.makedirs(workdir, exist_ok=True)
cfg = _deepreg_config(workdir, ct, t1c, checkpoint)
outdir = os.path.join(workdir, "pred")
subprocess.run(
["python", "-m", "deepreg.main", "test",
"-checkpoint_path", checkpoint, "-config_path", cfg, "-out_dir", outdir],
check=True, capture_output=True, text=True, cwd=workdir)
# predicted deformation field (voxels) on the fixed (T1c) grid
pred = _find_single_nii(os.path.join(outdir))
if pred is None:
raise PipelineUnavailable("deepreg: no predicted deformation field found")
disp = sitk.GetArrayFromImage(sitk.ReadImage(pred)) # (Z, Y, X, 3)
T = _displacement_to_sitk(disp, t1c)
meta = {
"pipeline": "deepreg", "toolkit": "Oxford DeepReg (CycleGAN)",
"metric": "CycleGAN + registration loss (pretrained)",
"transform": "3D displacement field", "convention": "T1c->CT (verified)",
"auto_inverted": False, "checkpoint": checkpoint,
}
return T, meta
def _deepreg_config(workdir, ct, t1c, checkpoint): def _lta_to_sitk(path):
"""Minimal DeepReg test config: moving=CT, fixed=T1c, 3D, single channel.""" """Parse a FreeSurfer LTA (LINEAR_RAS_TO_RAS) into a SimpleITK AffineTransform.
from deepreg.registry import registry # noqa
cfg = { The 4x4 matrix block lives in the shared RAS world space. Both inputs are
"general": {"data_format": "nifti", "device": "cuda:0"}, RAS-oriented by ``read_nii_img`` (DICOMOrient), so RAS == SimpleITK's
"data": {"moving": os.path.join(workdir, "moving"), physical frame and the matrix applies directly to world coordinates."""
"fixed": os.path.join(workdir, "fixed"), lines = [l for l in open(path).read().splitlines()
"field": "deformation", "moving_fields": "t1c", "fixed_fields": "t1c", if l and not l.startswith("#")]
"mask": None, "resample_size": "30 186 30"}, i = lines.index("1 4 4")
"validation": {}, M = np.array([[float(x) for x in lines[i + 1 + r].split()] for r in range(4)])
"training": {"epoch": 1, "batch_size": 1, T = sitk.AffineTransform(3)
"moving_fields": "t1c", "fixed_fields": "t1c", T.SetMatrix(M[:3, :3].ravel())
"resample_size": "30 186 30", "resample_spacing": "1 1 1"}, T.SetTranslation(M[:3, 3])
"testing": {"moving": os.path.join(workdir, "moving"), return T
"fixed": os.path.join(workdir, "fixed"),
"resample_size": "30 186 30", "resample_spacing": "1 1 1"},
} def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
p = os.path.join(workdir, "deepreg_test.yaml") 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:
import yaml # fixed = CT, moving = T1c (shared convention for all pipelines)
with open(p, "w") as f: sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
yaml.safe_dump(cfg, f) os.path.join(tmp, "fixed.nii.gz"), True) # CT
except Exception: sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
with open(p, "w") as f: os.path.join(tmp, "moving.nii.gz"), True) # T1c
json.dump(cfg, f)
return p
# run as the host uid/gid so the container-written LTA is host-owned
# (clean scratch removal) and TF still reaches the GPU.
cmd = ["docker", "run", "--rm",
"--user", f"{os.getuid()}:{os.getgid()}",
"-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"]
if gpu >= 0:
cmd += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"]
cmd += ["-v", f"{tmp}:/data", SYNTHMORPH_IMAGE,
"register", "-m", model]
if gpu >= 0:
cmd.append("-g")
if threads:
cmd += ["-j", str(threads)]
cmd += ["-e", str(extent),
"-T", "/data/f2m.lta",
"/data/moving.nii.gz", "/data/fixed.nii.gz"]
t0 = time.time()
proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200)
dt = time.time() - t0
if proc.returncode != 0:
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
raise RuntimeError(f"synthmorph: docker failed (rc={proc.returncode}):\n{tail}")
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}")
def _find_single_nii(dirpath): T = _lta_to_sitk(os.path.join(tmp, "f2m.lta"))
if not os.path.isdir(dirpath): T, inverted = fix_orientation(ct, t1c, T)
return None finally:
ns = [os.path.join(dirpath, x) for x in sorted(os.listdir(dirpath)) if not keep_scratch:
if x.endswith(".nii.gz") or x.endswith(".nii")] shutil.rmtree(tmp, ignore_errors=True)
return ns[0] if ns else None
meta = {
def _displacement_to_sitk(disp, ref): "pipeline": "synthmorph",
"""Build a sitk.DisplacementFieldTransform from a (Z,Y,X,3) voxel field.""" "toolkit": f"Docker {SYNTHMORPH_IMAGE} (TensorFlow, pretrained)",
from scipy.ndimage import map_coordinates "metric": "SynthMorph anatomy-aware contrast-invariant (deep, pretrained)",
zz, yy, xx = np.mgrid[0:disp.shape[0], 0:disp.shape[1], 0:disp.shape[2]] "transform": (f"SynthMorph {model} (identity start; "
# disp[i,j,k] = (dz, dy, dx) to add at grid point (i,j,k) to reach the source f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, "
src_z = zz + disp[..., 0]; src_y = yy + disp[..., 1]; src_x = xx + disp[..., 2] f"extent {extent})"),
# DisplacementFieldTransform stores, per output voxel, the offset to the input "convention": "T1c->CT (verified)", "auto_inverted": inverted,
off = sitk.GetImageFromArray(disp[..., ::-1].astype(np.float32)) # (X,Y,Z,3) "fixed": "CT", "moving": "T1c",
off.CopyInformation(ref) "seconds_reg": round(dt, 1),
t = sitk.DisplacementFieldTransform(off) "docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu,
return t "extent": extent,
}
if keep_scratch:
meta["scratch"] = tmp
return T, meta
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
@ -618,7 +660,7 @@ def _impact_docker_ok():
_NGPUS_CACHE = None _NGPUS_CACHE = None
def _impact_n_gpus(): def _n_gpus():
global _NGPUS_CACHE global _NGPUS_CACHE
if _NGPUS_CACHE is None: if _NGPUS_CACHE is None:
n = 0 n = 0
@ -752,6 +794,13 @@ 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
@ -802,14 +851,15 @@ 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: elastix maps fixed(T1c)->moving(CT), i.e. the # shared centered start (REQUIRED: the local ASGD optimizer does not recover
# inverse of centroid_offset, hence the negation. # the gross offset from identity). Elastix maps fixed(CT)->moving(T1c), the
# 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 = _impact_n_gpus() ngpu = _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
@ -885,7 +935,8 @@ 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)) # translation: T = T_euler(T_init(p)). Both map fixed(CT)->moving(T1c),
# 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)
@ -913,11 +964,11 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
REGISTERS = { REGISTERS = {
"simpleitk": register_simpleitk,
"elastix": register_elastix,
"antspyx": register_antspyx, "antspyx": register_antspyx,
"deepreg": register_deepreg, "elastix": register_elastix,
"impact": register_impact, "impact": register_impact,
"simpleitk": register_simpleitk,
"synthmorph": register_synthmorph,
} }
ALL_PIPELINES = list(REGISTERS) ALL_PIPELINES = list(REGISTERS)
@ -943,12 +994,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 == "deepreg": elif name == "synthmorph":
try: ok, note = _synthmorph_docker_ok()
import deepreg # noqa if not ok:
out[name] = (True, "deepreg importable") out[name] = (False, note)
except Exception as e: # noqa else:
out[name] = (False, "deepreg not importable") out[name] = (True, f"Docker {SYNTHMORPH_IMAGE}")
elif name == "impact": elif name == "impact":
ok, note = _impact_docker_ok() ok, note = _impact_docker_ok()
if not ok: if not ok: