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

feat(benchmark): switch to per-toolkit native init and full-volume registration

Browse source 
Replace the shared CT head-box resampling and manual centroid pre-shift with
each toolkit's native automatic transform initialization (elastix/impact use
CenterOfGravity, antspyx uses center-of-mass, simpleitk uses identity). All
pipelines now register on the full volumes instead of a cropped head box.

Drop the _ct_head_box_mm helper and update the benchmark table to document
per-toolkit init methods in a dedicated "init" column.
This commit is contained in:
Furen Xiao 2026-09-29 07:43:27 +08:00
parent 4ad972b00f
commit 8c083abe92
2 changed files with 230 additions and 308 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

263
REGISTRATION_BENCHMARK.md
View file

@ -22,34 +22,35 @@ 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
**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c, register on the **full CT / T1c volumes** (no head box) and use each toolkit's
moving = CT** (swapping its roles breaks the TotalSegmentator-based metric, **native automatic transform initialization** (no manual centroid offset):
which aborts with "Error in metric" or converges to a wrong rotation). The elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
intensity-metric pipelines register on the same CT-space head box (MI computed antspyx its default center-of-mass alignment, simpleitk starts from identity,
on head tissue, not neck/air); `synthmorph` uses the full volumes (it synthmorph needs none. All use **fixed = CT, moving = T1c** except `impact`,
normalizes each one to an isotropic network space internally). The local which keeps **fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much
optimizers that cannot recover the gross inter-session offset from identity larger than the T1c crop (~220×220×150 mm), so with CT as the fixed image too
(`antspyx`, `impact`) start from a centroid-aligned translation many CT samples map outside the T1c buffer and the IMPACT metric aborts
(`centroid_offset`); the others start from identity, so the comparison is on ("too many samples map outside moving image buffer"). The comparison is on
each toolkit's native optimizer / metric / global-search capability. each toolkit's native optimizer / metric / global-search capability.
## The five pipelines ## The five pipelines
Listed alphabetically. All use **fixed = CT, moving = T1c** except `impact`, Listed alphabetically. All register on the **full volumes** (no head box) and
which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator use each toolkit's **native automatic transform initialization** (no manual
metric). A pipeline that is a local optimizer and cannot recover the gross centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
inter-session offset from identity (`antspyx`, `impact`) starts from a (CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
centroid-aligned translation; the others start from identity / their own global starts from identity, and synthmorph needs none (global deep method). All use
search. **fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c, moving
= CT** (swapping breaks its TotalSegmentator metric — see below).
| # | pipeline | fixed / moving | centroid init | metric | transform | status in this env | | # | pipeline | fixed / moving | 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) | | 1 | `antspyx` | CT / T1c | auto COM | 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) | | 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (10 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
| 3 | `impact` | **T1c / CT** | **yes** | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) | | 3 | `impact` | **T1c / CT** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) |
| 4 | `simpleitk` | CT / T1c | no (identity) | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | | 4 | `simpleitk` | CT / T1c | 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) | | 5 | `synthmorph` | CT / T1c | none (global) | **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,80 +99,73 @@ 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** | | **antspyx** | ANTsPy 0.6.3 | 9 | **4/10 (0.40)** | 12.8 / 15.6 | 15.6 / 16.5 | 36.7 | 0.93 | **0.652 / 0.644** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 16.6 / 17.6 | 16.9 / 18.8 | 2.4 | 0.93 | **0.641 / 0.583** | | **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 13.6 / 15.7 | 14.8 / 17.0 | 2.6 | 0.93 | **0.641 / 0.640** |
| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **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 | **4/10 (0.40)** | 11.6 / 11.2 | 16.2 / 17.3 | **833.5** | 0.93 | **0.511 / 0.493** |
| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.7 / 16.0 | 13.5 / 15.7 | 7.6 | 0.93 | **0.645 / 0.632** | | **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 14.4 / 16.6 | 17.9 / 17.5 | 9.7 | 0.93 | **0.651 / 0.641** |
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 22.9 / 21.7 | 26.0 | 0.93 | **0.241 / 0.308** | | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 21.8 / 21.5 | 25.0 | 0.93 | **0.238 / 0.308** |
*(antspyx and impact are stochastic (random metric/patch subsampling), so their *(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 ~14
min/case on GPU, orders of magnitude above the intensity-metric tools. min/case on GPU (full volumes), orders of magnitude above the intensity-metric
elastix, simpleitk and synthmorph are deterministic at fixed threads.)* 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 | antspyx | elastix | impact | simpleitk | synthmorph |
|---|---|---|---|---|---| |---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | m6_22GD5VL2_20230724 | incomplete-CT | 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 | 24.6 / 28.6 | 24.0 / 28.3 | 22.7 / 27.0 | 26.3 / 29.6 | 10.4 / 23.7 |
| 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 | 22.3 / 23.5 | 9.1 / 16.2 | 21.9 / 23.2 | 20.3 / 23.3 |
| 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 | 9.5 / 15.6 | 11.6 / 14.8 | 8.6 / 7.2 **OK** | 14.4 / 23.2 | 11.1 / 21.8 |
| 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.7 / 7.3 **OK** | 12.2 / 7.4 **OK** | 4.6 / 6.6 **OK** | 11.7 / 7.5 **OK** | 11.5 / 8.4 **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 | 10.8 / 9.1 **OK** | 11.4 / 9.0 **OK** | 11.6 / 4.9 **OK** | 11.5 / 9.1 **OK** | 43.4 / 15.8 |
| 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.3 / 4.4 **OK** | 10.6 / 3.1 **OK** | 11.9 / 4.5 **OK** | 11.3 / 4.6 **OK** | 11.8 / 2.9 **OK** |
| m6_2BBN5XA4_20250411 | 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 | 19.1 / 7.7 **OK** | 18.2 / 14.1 | −1.7 / 29.3 | 18.8 / 7.5 **OK** | 52.4 / 53.0 |
| 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.8 / 17.2 | 13.6 / 20.0 | 15.4 / 23.2 | 14.0 / 17.9 | 16.7 / 28.2 |
| 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 | 17.5 / 34.6 | 17.6 / 32.8 | 18.9 / 37.1 | 19.3 / 35.0 | 18.9 / 16.8 |
### Interpretation ### Interpretation
- **Translation is recovered well by all five** (mid gaps mostly 6–20 mm, - **Translation is recovered well by all five** (mid gaps mostly 4–16 mm,
within the `≤ 20 mm` band, except the grossly-drifted cases below). within the `≤ 20 mm` band, except the grossly-drifted cases).
`tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for `tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for
every pipeline, so the target is in the right place. 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, and removing the head box makes it
(the ~30° inter-session tilt produces local minima that a single start harder.** On the repositioned cases `std > 12 mm` — the ~30° inter-session
(identity or centroid) can lock onto). This is exactly why production tilt produces local minima a single auto-start locks onto, and without the
`02_build_m6_dataset.py` uses a **multi-start tilt grid**. head-box crop the full-volume metric sees more neck/air confounders. Notably,
- **SimpleITK is the most reliable intensity-metric tool in this draw** the hardest tilt case `m6_22IRRKU6` is now **rejected by all five** (std
(4/10 accepted, best median mid 12.7), and **Elastix is the fastest** 23–30 mm), which the head-box config got through with IMPACT. This is exactly
(2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best why production `02_build_m6_dataset.py` uses a **multi-start tilt grid**.
median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c - **The easy cases are now accepted by nearly all pipelines.**
convention the two settle into different local minima on the hard tilted `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820` and `m6_26TYQI3C` are accepted
cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017` by ANTs/Elastix/IMPACT/SimpleITK (SynthMorph accepts the two cleanest). The
(std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles native automatic inits (center-of-gravity / center-of-mass) place the heads
better, and vice-versa on the drift cases. The fixed/moving direction close enough that the optimizers converge on the well-positioned cases.
materially affects where the ASGD/RSGD optimizers converge on the tilted cases. - **SimpleITK, ANTs and IMPACT tie at 4/10; Elastix 3/10; SynthMorph 2/10.**
- **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this ANTs is no longer the weakest — on the full volumes its native center-of-mass
draw; it is stochastic, so the count varies): it frequently locks onto the init works well (median std 15.6, up from the head-box+manual-pre-shift
tilted local minimum (median std 31.7, worst of all). It accepts only 31.7). It is stochastic (random MI subsampling), so the count varies by draw.
`m6_26TYQI3C` (mid 11.2 / std 3.7) and misses `m6_24V62QQM_20230223` - **IMPACT (semantic metric) is 4/10 but the slowest (~14 min/case on full
(std 47.5) and `m6_24V62QQM_20231017` (std 37.9) that the other tools accept — volumes) and lost its standout on the hardest tilt** (`m6_22IRRKU6` now
textbook tilted local minima. mid 22.7 / std 27.0 vs mid 8.3 / std 13.2 with the head box). Matching
- **IMPACT (semantic metric) is the standout on the hardest repositioned case** *anatomical semantics* still gives it the best median mid (11.6) and a
and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every competitive MI (0.511), but the full-volume registration is a larger problem
intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of for the local ASGD. It also drifts grossly on `m6_2BBN5XA4` (mid −1.7, the
any pipeline and within ~1 mm of the acceptance threshold — because matching transform flips to the far side).
*anatomical semantics* (TotalSegmentator features) is far more discriminating - **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.6 s
for gross pose than raw cross-modality intensity. But the semantic landscape (warm workers), SimpleITK ~9.7 s, ANTs ~36.7 s (CPU-bound, full volumes),
is rugged: single-init ASGD lands in bad local minima on several cases, and **SynthMorph ~25 s** (GPU: a single pretrained-model forward pass per case),
two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0). **IMPACT ~834 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator
Net 3/10 accepted — between ANTs and SimpleITK — at ~9.6 min/case, two orders
of magnitude slower than the intensity tools.
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.4 s
(warm workers), SimpleITK ~7.6 s, ANTs ~10.9 s (CPU-bound),
**SynthMorph ~26 s** (GPU: a single pretrained-model forward pass per case),
**IMPACT ~576 s** (GPU-bound: every ASGD iteration 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.645 (SimpleITK) / 0.567 (Elastix) / 0.469 (IMPACT) / 0.436 (ANTs) / medians 0.651 (SimpleITK) / 0.641 (Elastix) / 0.511 (IMPACT) / 0.652 (ANTs) /
0.241 (SynthMorph), 0.238 (SynthMorph).
and it diverges sharply on the tilted cases — ANTs `m6_24V62QQM_20230223` The intensity-metric tools cluster at ~0.64 (a genuine alignment signal),
drops to ~0.06 nats (rotation ~47 mm off) while SimpleITK holds 0.703. A high while SynthMorph's 0.238 reflects its tilt failures. A high MI means the CT
MI means the CT skull/tissue intensities fall on the matching MR structures skull/tissue intensities fall on the matching MR structures after the inverse
after the inverse transform. ANTs' and SynthMorph's low median MI is transform.
consistent with their tilted-case convergence failures.
- **SynthMorph (deep, pretrained) is near-weakest of the five on these - **SynthMorph (deep, pretrained) is near-weakest of the five on these
cross-modality CT pairs.** Its rigid model is trained on brain *MRI* cross-modality CT pairs.** Its rigid model is trained on brain *MRI*
(T1/T2, anatomy-aware but intensity-domain MRI), so it recovers gross (T1/T2, anatomy-aware but intensity-domain MRI), so it recovers gross
@ -181,7 +175,7 @@ Per-case (mid / std, `OK` = accepted):
mid 43). It accepts the clean `m6_26TYQI3C` (mid 11.8 / std 2.9, MI 0.60) and 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 `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 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) method that recovers the gross offset on its own). Its low median MI (0.238)
reflects the tilt failures. On same-modality MRI it would likely rank much 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 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 deterministic (one fixed pretrained model, no random subsampling), so unlike
@ -303,8 +297,10 @@ image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121)
*semantic* similarity metric compiled in as a native Elastix metric, plus a *semantic* similarity metric compiled in as a native Elastix metric, plus a
TorchScript **TotalSegmentator** feature extractor. It is GPU-based and TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
**fixed = T1c, moving = CT** (the reverse of the others): swapping its roles **fixed = T1c, moving = CT** (the reverse of the others): the CT FOV
breaks the TotalSegmentator metric ("Error in metric" / wrong rotation). (~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with
CT as the fixed image too many CT samples map outside the T1c buffer and the
IMPACT metric aborts ("too many samples map outside moving image buffer").
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):
@ -330,37 +326,36 @@ repo for the full list and FOV/layer tables.
### How the adapter drives the image ### How the adapter drives the image
The image's `CMD.sh` only runs `elastix -f -m -p -out` (no `-t0`), so the The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`) binary directly (`/lib/elastix-install/bin/elastix`):
to add the shared centered start:
1. **Head box** — resamples fixed=T1c and moving=CT into the same CT-space head 1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box).
box (physical frame) the other pipelines use, so MI/IMPACT is computed on T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too
head tissue, not neck/air. many CT samples map outside the T1c buffer and the metric aborts.
2. **Shared init via `-t0`** — the `centroid_offset` translation is written to a 2. **Auto init (no `-t0`)** — the parameter map sets
`TranslationTransform` parameter file and passed with `-t0`; with `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
`HowToCombineTransforms "Compose"` the final transform is elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
`T = T_euler(T_init(p))`. (Elastix 5.x removed in-file multi-`Transform` optimizer scales.
chains and `UseScalingFactors`, so the initial transform must be a separate
`-t0` file.)
3. **Parameter map** — rigid `EulerTransform` (6-DOF), `Impact` metric in 3. **Parameter map** — rigid `EulerTransform` (6-DOF), `Impact` metric in
**Jacobian** mode (differentiable, per the recommended config), 4 semantic **Jacobian** mode (differentiable, per the recommended config), 4 semantic
scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the box's base spacing) with scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the T1c base spacing) with
**no image resampling** (all-1 pyramid schedules, as the IMPACT docs require), **no image resampling** (all-1 pyramid schedules, as the IMPACT docs require),
ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid
(`--impact-gpu -1` forces CPU). (`--impact-gpu -1` forces CPU).
4. **Read-back** — the result `TransformParameters.0.itk.txt` is an 4. **Read-back** — the result `TransformParameters.0.itk.txt` is an
`Euler3DTransform` in the box/world physical frame; the adapter composes it `Euler3DTransform` in the world physical frame (the
with the `-t0` translation and runs the usual `fix_orientation` gate. `AutomaticTransformInitialization` is folded in by elastix); the adapter
runs it through the usual `fix_orientation` gate.
## Notes / caveats ## Notes / caveats
- The benchmark compares each toolkit in its **native configuration**: the - The benchmark compares each toolkit in its **native configuration**: all
local optimizers that cannot recover the gross offset from identity register on the full volumes (no head box) and use each toolkit's native
(`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`, automatic transform initialization (no manual `centroid_offset`) — elastix /
`elastix` and `synthmorph` start from identity / their own global search. impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's
Production uses `--simpleitk-mode production`'s multi-start tilt grid, which center-of-mass alignment, simpleitk identity, synthmorph none. Production
raises SimpleITK's acceptance at ~20× the cost; that is out of scope for the uses `--simpleitk-mode production`'s multi-start tilt grid, which raises
SimpleITK's acceptance at ~20× the cost; that is out of scope for the
head-to-head but is the config the dataset actually ships with. head-to-head but is the config the dataset actually ships with.
- `synthmorph` is deterministic (one fixed pretrained model, no random - `synthmorph` is deterministic (one fixed pretrained model, no random
subsampling), so its per-case numbers are stable run to run — unlike ANTs / subsampling), so its per-case numbers are stable run to run — unlike ANTs /
@ -385,8 +380,8 @@ to add the shared centered start:
path). path).
- **GPU sharing.** The container runs as the host uid with - **GPU sharing.** The container runs as the host uid with
`TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU. `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
- **~26 s/case** (GPU: one pretrained forward pass), between ANTs (~11 s) - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK
and IMPACT (~576 s). (~10 s) / ANTs (~37 s) and IMPACT (~834 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
@ -394,11 +389,11 @@ 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.
- No initial transform is passed (identity start): the ASGD + 3-resolution - `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the
Mattes-MI pyramid recovers the gross inter-session offset on its own, so start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the
the shared `centroid_offset` start is not required. An **output directory** optimizer scales. An **output directory** is still set
is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter (`SetOutputDirectory`) for the C++ filter to run; 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.
@ -412,28 +407,32 @@ to add the shared centered start:
500/500/250 over 3 resolutions). 500/500/250 over 3 resolutions).
- **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 ~14 min/case on full volumes (vs 2–37 s
tools). `--impact-iters` / `--impact-samples` trade accuracy for speed. for the other tools). `--impact-iters` / `--impact-samples` trade accuracy
for speed.
- **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT - **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT
follow the shared fixed=CT/moving=T1c convention: swapping its roles follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251
breaks the TotalSegmentator metric (elastix aborts with "Error in metric" mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the
or converges to a wrong rotation). The returned transform still satisfies fixed image too many CT samples map outside the T1c buffer and the metric
the contract (T1c→CT). aborts ("too many samples map outside moving image buffer"). The returned
transform still satisfies the contract (T1c→CT).
- **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns
the intensity centroids; no manual `-t0` translation is used.
- **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: on full volumes it no longer nails the hardest
(`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box
others. A multi-start or higher-iteration budget would likely raise its used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration
acceptance; that is out of scope for the single-init head-to-head. budget would likely raise its acceptance; that is out of scope for the
- **Stochastic.** Jacobian mode resamples patches each iteration, so results single-init head-to-head.
vary run to run (see Stochasticity). - **Stochastic.** Jacobian mode resamples patches each iteration, so results
- **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image vary run to run (see Stochasticity).
bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors` - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image
are gone — the shared centered start is a separate `-t0` translation file, bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors`
and the result is read back as an `Euler3DTransform` (`.itk.txt`) and are gone. The result is read back as an `Euler3DTransform` (`.itk.txt`).
composed with `-t0`. The image's own `CMD.sh` doesn't pass `-t0`, so the The image's own `CMD.sh` is a generic launcher, so the adapter invokes
adapter invokes `/lib/elastix-install/bin/elastix` directly. `/lib/elastix-install/bin/elastix` directly.
- **Models are not in the image.** The TorchScript model(s) must be present - **Models are not in the image.** The TorchScript model(s) must be present
under the mounted `/Data/Models` (the repo caches them in under the mounted `/Data/Models` (the repo caches them in
`data/impact_models/`, bind-mounted read-only). `data/impact_models/`, bind-mounted read-only).

275
src/registration_pipelines.py
View file

@ -23,34 +23,35 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
Pipelines (listed alphabetically) Pipelines (listed alphabetically)
--------------------------------- ---------------------------------
All pipelines use the shared convention **fixed = CT, moving = T1c**, except All pipelines register on the **full volumes** (no head box) and use each
IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its toolkit's **native automatic transform initialization** (no manual
TotalSegmentator metric -- see ``register_impact``). A pipeline that is a ``centroid_offset``): elastix / impact use elastix's
local optimizer and cannot recover the gross inter-session offset from ``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its
identity (antspyx, impact) starts from a centroid-aligned translation; the default center-of-mass alignment, simpleitk starts from identity, and
others (simpleitk, elastix, synthmorph) do not require it and start from synthmorph needs none (global deep method). All use **fixed = CT, moving =
identity / their own global search. T1c** except IMPACT, which keeps **fixed = T1c, moving = CT** (swapping breaks
its TotalSegmentator metric -- see ``register_impact``).
- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). - ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
Centroid pre-shift of the moving (T1c) image is *required* (ANTs Default ``initial_transform=None`` computes a translation aligning the image
Rigid-from-identity collapses on the gross offset). centers of mass.
- ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` + - ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` +
``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD). ``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
Identity start (the ASGD + multi-resolution pyramid recovers the offset). ``AutomaticTransformInitialization`` (CenterOfGravity) start.
- ``register_impact`` : IMPACT (semantic-similarity metric) inside the - ``register_impact`` : IMPACT (semantic-similarity metric) inside the
``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript
TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF, TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF,
GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); centroid init GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix
kept (local ASGD on the semantic metric benefits from the gross alignment). ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
Needs Docker + the image + a model asset. Needs Docker + the image + a model asset.
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a - ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on the
CT-space head box. Reuses the *proven* production ``register_ct_to_t1c`` in full volumes. Reuses the *proven* production ``register_ct_to_t1c`` in
``mode="production"``; the benchmark mode is a single identity-start ``mode="production"``; the benchmark mode is a single identity-start
multi-resolution pyramid run. multi-resolution pyramid run.
- ``register_synthmorph`` : FreeSurfer SynthMorph (TensorFlow, pretrained - ``register_synthmorph`` : FreeSurfer SynthMorph (TensorFlow, pretrained
anatomy-aware contrast-invariant model) inside the ``freesurfer/synthmorph`` anatomy-aware contrast-invariant model) inside the ``freesurfer/synthmorph``
Docker image. Rigid 6-DOF matrix model, GPU when available. Global deep Docker image. Rigid 6-DOF matrix model, GPU when available. Global deep
method -- no centroid init. Needs Docker + the pulled image. method -- no 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 the Docker-based pipelines, missing
its image) its image)
@ -157,38 +158,18 @@ def centroid_offset(ct, t1c, ct_box=None):
return (float(t1c_c[0] - c[0]), float(t1c_c[1] - c[1]), float(t1c_c[2] - c[2])) return (float(t1c_c[0] - c[0]), float(t1c_c[1] - c[1]), float(t1c_c[2] - c[2]))
def _ct_head_box_mm(ct):
"""CT-space head box in physical mm: (box_lo, box_hi) tuples, or None."""
m6 = _load_m6()
ca = sitk.GetArrayFromImage(ct)
tissue = ca > -150
csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
vz = m6._ct_vertex_slice(ct)
if vz is None:
return None, None
ztop = co[2] + vz * csp[2]
band = tissue[max(0, vz - int(round(150 / csp[2]))):vz + 1]
cols = np.nonzero(band.any(0))
box_lo = (co[0] + cols[1].min() * csp[0] - 15,
co[1] + cols[0].min() * csp[1] - 15, ztop - 215)
box_hi = (co[0] + cols[1].max() * csp[0] + 15,
co[1] + cols[0].max() * csp[1] + 15, ztop + 10)
return box_lo, box_hi
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# 2. SimpleITK (production baseline -- reused verbatim) # 2. SimpleITK (production baseline -- reused verbatim)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150): def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
"""SimpleITK rigid CT->T1c, Mattes MI on a CT-space head box. """SimpleITK rigid CT->T1c, Mattes MI, full volumes (no head box).
mode="benchmark" (default): a single multi-resolution pyramid run from the mode="benchmark" (default): a single multi-resolution pyramid run from the
identity start, capped at ``iters_full`` iterations at full resolution. identity start, 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. No centroid init is used.
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).
@ -215,20 +196,16 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
except AttributeError: except AttributeError:
pass pass
ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32) ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
lo, hi = _ct_head_box_mm(ct) # full volumes, fixed = CT, moving = T1c (no head box)
if lo is None: fixed, moving = ct, t1c
raise PipelineUnavailable("simpleitk: no CT head box found (corrupt CT?)") shape = np.array(ct.GetSize())
tsp = np.array(t1c.GetSpacing()) # Identity start (no centroid init); the multi-resolution Mattes-MI pyramid
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3)) # recovers the gross inter-session offset. The Euler center is the fixed
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32)) # (CT) geometric center.
ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection()) init = sitk.Euler3DTransform()
ref.SetOrigin(tuple(float(x) for x in lo)) init.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint(
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2,
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) (fixed.GetSize()[2] - 1) // 2))))
# Identity start: the multi-resolution Mattes-MI pyramid recovers the gross
# inter-session offset on its own, so no centroid init is required here.
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.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
import time import time
t0 = time.time() t0 = time.time()
@ -250,9 +227,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
meta = { meta = {
"pipeline": "simpleitk", "mode": "benchmark", "pipeline": "simpleitk", "mode": "benchmark",
"toolkit": f"SimpleITK {sitk.Version_VersionString()}", "toolkit": f"SimpleITK {sitk.Version_VersionString()}",
"metric": "MattesMutualInformation (head box, CT space)", "metric": "MattesMutualInformation (full volumes)",
"transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, " "transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
f"{iters_full} iters/level)"), f"{iters_full} iters/level, identity start)"),
"convention": "T1c->CT (benchmark budget)", "auto_inverted": False, "convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
"seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm, "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
"skull_ring_gap_mm": gap, "skull_ring_gap_mm": gap,
@ -311,34 +288,26 @@ 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=CT, moving=T1c), Mattes MI, full images.
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 full
CT-space head box, from the identity start (the ASGD + multi-resolution CT / T1c volumes (no head box). No manual centroid init: elastix's
pyramid recovers the gross inter-session offset, so no centroid init is ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity
required). Elastix 5.x component names: sampler ``Random`` (not centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales.
``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``). An output Elastix 5.x
directory is required for the C++ filter to run, and ``WriteResultImage`` component names: sampler ``Random`` (not ``RandomSampler``), metric key
must be the string ``"false"`` (a Python bool False) or the filter raises ``Metric`` (not ``ImageMetric``). An output directory is required for the
"Could not read result image" after a successful registration. C++ filter to run, 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) # full volumes, fixed = CT, moving = T1c (no head box)
if lo is None: fixed = ct
raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)") moving = t1c
# Resample both into a common CT-space head box (mirrors the other pipelines)
# so MI is computed on head tissue, not neck/air.
tsp = np.array(t1c.GetSpacing())
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
ref.SetSpacing(t1c.GetSpacing())
ref.SetDirection(t1c.GetDirection())
ref.SetOrigin(tuple(float(x) for x in lo))
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = reference
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = moving
params = { params = {
"ImageDimension": 3, "ImageDimension": 3,
"FixedInternalImagePixelType": "float", "FixedInternalImagePixelType": "float",
@ -348,6 +317,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
"FixedCenteredTransform": True, "FixedCenteredTransform": True,
"TransformParameterReadoutFormat": 1, "TransformParameterReadoutFormat": 1,
"HowToCombineTransforms": "Compose", "HowToCombineTransforms": "Compose",
"AutomaticTransformInitialization": True,
"AutomaticTransformInitializationMethod": "CenterOfGravity",
"AutomaticScalesEstimation": True,
"Optimizer": "AdaptiveStochasticGradientDescent", "Optimizer": "AdaptiveStochasticGradientDescent",
"MaximumNumberOfIterations": [500, 500, 250], "MaximumNumberOfIterations": [500, 500, 250],
"Metric": "AdvancedMattesMutualInformation", "Metric": "AdvancedMattesMutualInformation",
@ -388,8 +360,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
version = "ITK-Elastix (python API)" version = "ITK-Elastix (python API)"
meta = { meta = {
"pipeline": "elastix", "toolkit": version, "pipeline": "elastix", "toolkit": version,
"metric": "MattesMutualInformation (head box, CT space)", "metric": "MattesMutualInformation (full volumes)",
"transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)", "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)",
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1), "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
} }
@ -402,18 +374,16 @@ 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"):
"""ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c). """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c), full volumes.
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=...,
aff_smoothing_sigmas=...)``. ``aff_metric='mattes'`` (the default) gives a aff_smoothing_sigmas=...)``. ``aff_metric='mattes'`` (the default) gives a
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. Registers on the full CT / T1c volumes (no head box).
the other pipelines (and to avoid neck/air in the metric), both images are No manual centroid init: antspyx's default ``initial_transform=None``
first resampled into the same CT-space head box. ANTs Rigid-from-identity computes a translation that aligns the image centers of mass. Returns the
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 transform that warps CT into the T1c grid (orientation verified by
:func:`fix_orientation`). :func:`fix_orientation`).
""" """
@ -428,36 +398,12 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
aff_iterations = tuple(aff_iterations) if aff_iterations else (500, 500, 200) aff_iterations = tuple(aff_iterations) if aff_iterations else (500, 500, 200)
aff_shrink_factors = tuple(aff_shrink_factors) if aff_shrink_factors else (4, 2, 1) aff_shrink_factors = tuple(aff_shrink_factors) if aff_shrink_factors else (4, 2, 1)
aff_smoothing_sigmas = tuple(aff_smoothing_sigmas) if aff_smoothing_sigmas else (3, 2, 0) aff_smoothing_sigmas = tuple(aff_smoothing_sigmas) if aff_smoothing_sigmas else (3, 2, 0)
if use_head_box: # full volumes, fixed = CT, moving = T1c (no head box)
lo, hi = _ct_head_box_mm(ct) fx, mv = ct, t1c
if lo is None:
raise PipelineUnavailable("antspyx: no CT head box found (corrupt CT?)")
tsp = np.array(t1c.GetSpacing())
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
ref.SetOrigin(tuple(float(x) for x in lo))
fx = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (fixed)
mv = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (moving)
else:
fx, mv = ct, t1c
import time import time
t0 = time.time() t0 = time.time()
# Centered start (REQUIRED for ANTs): antspyx cannot use a # No manual centroid init: antspyx's default initial_transform=None computes
# TranslationTransform as initial_transform (the ANTs CLI rejects it), and # a translation aligning the image centers of mass.
# Rigid-from-identity does not recover the large (~100-130 mm) inter-session
# offset -- without this init ANTs collapses on 8/9 cases. So we pre-shift the
# T1c (moving) box so its head centroid coincides with the CT (fixed) head
# centroid, run Rigid from identity, then compose the shift + inverse back so
# the result is expressed in the original T1c/CT coordinates.
off = None
if center_init:
off = centroid_offset(ct, t1c, ct_box=fx)
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.SetOffset((off[0], off[1], off[2]))
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))
r = ants.registration(fixed=af, moving=am, type_of_transform=transform_type, r = ants.registration(fixed=af, moving=am, type_of_transform=transform_type,
@ -466,20 +412,15 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
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 the ANTs forward transform: fixed(CT) -> moving
# (pre-shifted T1c). We need the contract T1c(orig) -> CT: map T1c(orig) to # (T1c). The contract is T1c -> CT, i.e. its inverse; fix_orientation
# the pre-shifted T1c (translation by -off), then invert the ANTs transform. # verifies the direction against the skull ring.
T_ants = sitk.ReadTransform(r["fwdtransforms"][0]) T_ants = sitk.ReadTransform(r["fwdtransforms"][0])
if off is not None: cand = T_ants.GetInverse()
tr_shift = sitk.TranslationTransform(3)
tr_shift.SetOffset((-off[0], -off[1], -off[2])) # T1c orig -> pre
cand = sitk.CompositeTransform([tr_shift, T_ants.GetInverse()])
else:
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__', '?')}",
"metric": f"ANTs {aff_metric} (multi-resolution)", "metric": f"ANTs {aff_metric} (multi-resolution, full volumes)",
"transform": f"ANTs {transform_type} (head box)", "convention": "T1c->CT (verified)", "transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)",
"auto_inverted": inverted, "aff_iterations": list(aff_iterations), "auto_inverted": inverted, "aff_iterations": list(aff_iterations),
"aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1), "aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
} }
@ -740,12 +681,18 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
"(ImpactFeaturesMapUpdateInterval -1)", "(ImpactFeaturesMapUpdateInterval -1)",
'(ImpactWriteFeatureMaps "false")', '(ImpactWriteFeatureMaps "false")',
"", "",
"// rigid 6-DOF to be optimized (the shared centered start comes via -t0)", "// rigid 6-DOF to be optimized (auto center init, no -t0)",
'(Transform "EulerTransform")', '(Transform "EulerTransform")',
'(TransformBase "ElastixTransformBase")', '(TransformBase "ElastixTransformBase")',
'(TransformParameterReadoutFormat "AllParametersOneLine")', '(TransformParameterReadoutFormat "AllParametersOneLine")',
'(FixedCenteredTransform "true")', '(FixedCenteredTransform "true")',
'(HowToCombineTransforms "Compose")', '(HowToCombineTransforms "Compose")',
'(AutomaticTransformInitialization "true")',
# align intensity centroids (heads), not buffer corners -- the CT FOV is
# much larger than the T1c crop, so GeometricalCenter leaves the heads
# ~100 mm apart and the metric aborts ("too many samples outside buffer")
'(AutomaticTransformInitializationMethod "CenterOfGravity")',
'(AutomaticScalesEstimation "true")',
"", "",
"(ImageDimension 3)", "(ImageDimension 3)",
'(FixedInternalImagePixelType "float")', '(FixedInternalImagePixelType "float")',
@ -794,13 +741,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
@ -809,21 +749,26 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
the Mattes-MI pipelines above -- it is the cross-modality tool the image was the Mattes-MI pipelines above -- it is the cross-modality tool the image was
built for. built for.
To stay on the benchmark's equal footing the pipeline: The pipeline registers on the full CT / T1c volumes (no head box) with no
- registers on the same CT-space head box (both images resampled into it), manual centroid init: elastix's ``AutomaticTransformInitialization``
- starts from the *shared* centroid-aligned translation (passed to Elastix (CenterOfGravity) aligns the intensity centroids. It uses a rigid 6-DOF
as a fixed initial transform via ``-t0``; with ``HowToCombineTransforms ``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level
Compose`` the final transform = ``T_main(T_init(p))``), scale pyramid (no image resampling).
- uses a rigid 6-DOF ``EulerTransform`` (the 5.x name for the classic
3D rigid), ASGD, a 4-level scale pyramid (no image resampling). NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other
pipelines). The CT FOV (~300x300x251 mm) is much larger than the T1c crop
(~220x220x150 mm); with CT as the fixed image, too many CT samples map
outside the T1c buffer and the IMPACT metric aborts ("too many samples map
outside moving image buffer"). The returned transform still satisfies the
unified contract (T1c->CT).
Budget: the default ``iters=500`` / ``samples=2000`` is the authors' Budget: the default ``iters=500`` / ``samples=2000`` is the authors'
recommended Jacobian budget. It is *much* slower than the intensity-metric recommended Jacobian budget. It is *much* slower than the intensity-metric
pipelines (minutes per case on GPU vs seconds) because every iteration runs pipelines (minutes per case on GPU vs seconds) because every iteration runs
the TotalSegmentator feature extractor on the sampled patches with the TotalSegmentator feature extractor on the sampled patches with
backpropagation. Lower ``iters``/``samples`` to trade accuracy for speed. backpropagation. Lower ``iters``/``samples`` to trade accuracy for speed.
The result is read back from ``TransformParameters.0.itk.txt`` The result is read back from ``TransformParameters.0.itk.txt``
(an ``Euler3DTransform`` in the box/world physical frame) and composed with (an ``Euler3DTransform`` in the world physical frame) and checked by
the ``-t0`` translation, then checked by :func:`fix_orientation`. :func:`fix_orientation`.
Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, the Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, the
``vboussot/elastix_impact`` image, and the model asset under the cache. A ``vboussot/elastix_impact`` image, and the model asset under the cache. A
@ -837,25 +782,13 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
raise PipelineUnavailable(f"impact: {note}") raise PipelineUnavailable(f"impact: {note}")
model_host, model_rel = _impact_model_host(model) model_host, model_rel = _impact_model_host(model)
lo, hi = _ct_head_box_mm(ct) # full volumes, fixed = T1c, moving = CT (no head box, no manual init --
if lo is None: # elastix's AutomaticTransformInitialization aligns the image centers).
raise PipelineUnavailable("impact: no CT head box found (corrupt CT?)") # T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger
# than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples
# resample both into the shared CT-space head box (mirrors the other pipelines) # map outside the T1c buffer and the IMPACT metric aborts.
tsp = np.array(t1c.GetSpacing()) fixed = t1c
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3)) moving = ct
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
ref.SetSpacing(t1c.GetSpacing())
ref.SetDirection(t1c.GetDirection())
ref.SetOrigin(tuple(float(x) for x in lo))
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
# shared centered start (REQUIRED: the local ASGD optimizer does not recover
# 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)
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:
@ -899,15 +832,8 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
os.path.join(datadir, "Moving_image.mha"), True) os.path.join(datadir, "Moving_image.mha"), True)
with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
f.write(_impact_param_map(os.path.join("/Data/Models", model_rel), f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
gpu, iters, samples, voxel_sizes, gpu, iters, samples, voxel_sizes,
rescale=schedule, smoothing=smoothe)) rescale=schedule, smoothing=smoothe))
with open(os.path.join(datadir, "InitialTransform.txt"), "w") as f:
f.write(
'(Transform "TranslationTransform")\n'
'(TransformBase "ElastixTransformBase")\n'
'(TransformParameterReadoutFormat "AllParametersOneLine")\n'
"(NumberOfParameters 3)\n"
f"(TransformParameters {init[0]:g} {init[1]:g} {init[2]:g})\n")
cmd = ["docker", "run", "--rm"] cmd = ["docker", "run", "--rm"]
if gpu >= 0: if gpu >= 0:
@ -919,7 +845,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
IMPACT_IMAGE, IMPACT_ELASTIX_BIN, IMPACT_IMAGE, IMPACT_ELASTIX_BIN,
"-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha", "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha",
"-p", "/Data/ParameterMap.txt", "-p", "/Data/ParameterMap.txt",
"-t0", "/Data/InitialTransform.txt",
"-out", "/Out/", "-threads", str(threads), "-out", "/Out/", "-threads", str(threads),
] ]
t0 = time.time() t0 = time.time()
@ -934,14 +859,12 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
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 (world frame) is the full combination transform
# translation: T = T_euler(T_init(p)). Both map fixed(CT)->moving(T1c), # (AutomaticTransformInitialization is folded in by elastix); it maps
# i.e. the inverse of the contract (T1c->CT); fix_orientation resolves it. # fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT);
# fix_orientation resolves the direction.
T_main = sitk.ReadTransform(itk_tf) T_main = sitk.ReadTransform(itk_tf)
T_init = sitk.TranslationTransform(3) T, inverted = fix_orientation(ct, t1c, T_main)
T_init.SetOffset(init)
cand = sitk.CompositeTransform([T_init, T_main])
T, inverted = fix_orientation(ct, t1c, cand)
finally: finally:
if not keep_scratch: if not keep_scratch:
shutil.rmtree(tmp, ignore_errors=True) shutil.rmtree(tmp, ignore_errors=True)
@ -950,7 +873,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
"pipeline": "impact", "pipeline": "impact",
"toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
"metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)", "metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)",
"transform": (f"Rigid (EulerTransform, ASGD, {n} scales, centroid init; " "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; "
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": "T1c", "moving": "CT", "fixed": "T1c", "moving": "CT",