From 8c083abe92a320f81b2f9086ae9a82e82bb7e515 Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Tue, 29 Sep 2026 07:43:27 +0800 Subject: [PATCH] feat(benchmark): switch to per-toolkit native init and full-volume registration 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. --- REGISTRATION_BENCHMARK.md | 263 ++++++++++++++++---------------- src/registration_pipelines.py | 275 ++++++++++++---------------------- 2 files changed, 230 insertions(+), 308 deletions(-) diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index 6717c75..4c93ec6 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/REGISTRATION_BENCHMARK.md @@ -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 -destination→source convention this `T` maps **T1c → CT**. All pipelines use -**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c, -moving = CT** (swapping its roles breaks the TotalSegmentator-based metric, -which aborts with "Error in metric" or converges to a wrong rotation). The -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 +destination→source convention this `T` maps **T1c → CT**. All pipelines +register on the **full CT / T1c volumes** (no head box) and use each toolkit's +**native automatic transform initialization** (no manual centroid offset): +elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity), +antspyx its default center-of-mass alignment, simpleitk starts from identity, +synthmorph needs none. All use **fixed = CT, moving = T1c** except `impact`, +which keeps **fixed = T1c, moving = CT**: the CT FOV (~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"). The comparison is on each toolkit's native optimizer / metric / global-search capability. ## The five pipelines -Listed alphabetically. All use **fixed = CT, moving = T1c** except `impact`, -which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator -metric). A pipeline that is a local optimizer and cannot recover the gross -inter-session offset from identity (`antspyx`, `impact`) starts from a -centroid-aligned translation; the others start from identity / their own global -search. +Listed alphabetically. All register on the **full volumes** (no head box) and +use each toolkit's **native automatic transform initialization** (no manual +centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization` +(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk +starts from identity, and synthmorph needs none (global deep method). All use +**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 | -|---|----------|----------------|---------------|--------|-----------|--------------------| -| 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) | +| # | pipeline | fixed / moving | init | metric | transform | status in this env | +|---|----------|----------------|------|--------|-----------|--------------------| +| 1 | `antspyx` | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | +| 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** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) | +| 4 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | +| 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, 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) | |---|---|---|---|---|---|---|---|---| -| **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** | -| **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** | -| **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** | -| **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** | -| **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** | ANTsPy 0.6.3 | 9 | **4/10 (0.40)** | 12.8 / 15.6 | 15.6 / 16.5 | 36.7 | 0.93 | **0.652 / 0.644** | +| **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 13.6 / 15.7 | 14.8 / 17.0 | 2.6 | 0.93 | **0.641 / 0.640** | +| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **4/10 (0.40)** | 11.6 / 11.2 | 16.2 / 17.3 | **833.5** | 0.93 | **0.511 / 0.493** | +| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 14.4 / 16.6 | 17.9 / 17.5 | 9.7 | 0.93 | **0.651 / 0.641** | +| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 21.8 / 21.5 | 25.0 | 0.93 | **0.238 / 0.308** | *(antspyx and impact are stochastic (random metric/patch subsampling), so their -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. -elastix, simpleitk and synthmorph are deterministic at fixed threads.)* +accepted counts are one draw each — see Interpretation. impact's s/reg is ~14 +min/case on GPU (full volumes), orders of magnitude above the intensity-metric +tools. elastix, simpleitk and synthmorph are deterministic at fixed threads.)* Per-case (mid / std, `OK` = accepted): | case | antspyx | elastix | impact | simpleitk | synthmorph | |---|---|---|---|---|---| | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | -| 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.9 / 24.5 | 22.6 / 23.9 | 23.1 / 25.7 | 22.6 / 24.0 | 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_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 | 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 | 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 | 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 | 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 / 31.7 | 19.5 / 30.3 | 19.6 / 30.9 | 19.3 / 30.1 | 18.9 / 16.8 | +| 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.2 / 23.9 | 22.3 / 23.5 | 9.1 / 16.2 | 21.9 / 23.2 | 20.3 / 23.3 | +| 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 | 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 | 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 | 12.3 / 4.4 **OK** | 10.6 / 3.1 **OK** | 11.9 / 4.5 **OK** | 11.3 / 4.6 **OK** | 11.8 / 2.9 **OK** | +| m6_2BBN5XA4_20250411 | 19.1 / 7.7 **OK** | 18.2 / 14.1 | −1.7 / 29.3 | 18.8 / 7.5 **OK** | 52.4 / 53.0 | +| 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 | 17.5 / 34.6 | 17.6 / 32.8 | 18.9 / 37.1 | 19.3 / 35.0 | 18.9 / 16.8 | ### Interpretation -- **Translation is recovered well by all five** (mid gaps mostly 6–20 mm, - within the `≤ 20 mm` band, except the grossly-drifted cases below). +- **Translation is recovered well by all five** (mid gaps mostly 4–16 mm, + 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 every pipeline, so the target is in the right place. -- **The binding failure is rotation** — `std > 12 mm` on the repositioned cases - (the ~30° inter-session tilt produces local minima that a single start - (identity or centroid) can lock onto). This is exactly why production - `02_build_m6_dataset.py` uses a **multi-start tilt grid**. -- **SimpleITK is the most reliable intensity-metric tool in this draw** - (4/10 accepted, best median mid 12.7), and **Elastix is the fastest** - (2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best - median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c - convention the two settle into different local minima on the hard tilted - cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017` - (std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles - better, and vice-versa on the drift cases. The fixed/moving direction - materially affects where the ASGD/RSGD optimizers converge on the tilted cases. -- **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this - draw; it is stochastic, so the count varies): it frequently locks onto the - tilted local minimum (median std 31.7, worst of all). It accepts only - `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** - and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every - intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of - any pipeline and within ~1 mm of the acceptance threshold — because matching - *anatomical semantics* (TotalSegmentator features) is far more discriminating - for gross pose than raw cross-modality intensity. But the semantic landscape - is rugged: single-init ASGD lands in bad local minima on several cases, and - two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0). - 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 +- **The binding failure is rotation, and removing the head box makes it + harder.** On the repositioned cases `std > 12 mm` — the ~30° inter-session + tilt produces local minima a single auto-start locks onto, and without the + head-box crop the full-volume metric sees more neck/air confounders. Notably, + the hardest tilt case `m6_22IRRKU6` is now **rejected by all five** (std + 23–30 mm), which the head-box config got through with IMPACT. This is exactly + why production `02_build_m6_dataset.py` uses a **multi-start tilt grid**. +- **The easy cases are now accepted by nearly all pipelines.** + `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820` and `m6_26TYQI3C` are accepted + by ANTs/Elastix/IMPACT/SimpleITK (SynthMorph accepts the two cleanest). The + native automatic inits (center-of-gravity / center-of-mass) place the heads + close enough that the optimizers converge on the well-positioned cases. +- **SimpleITK, ANTs and IMPACT tie at 4/10; Elastix 3/10; SynthMorph 2/10.** + ANTs is no longer the weakest — on the full volumes its native center-of-mass + init works well (median std 15.6, up from the head-box+manual-pre-shift + 31.7). It is stochastic (random MI subsampling), so the count varies by draw. +- **IMPACT (semantic metric) is 4/10 but the slowest (~14 min/case on full + volumes) and lost its standout on the hardest tilt** (`m6_22IRRKU6` now + mid 22.7 / std 27.0 vs mid 8.3 / std 13.2 with the head box). Matching + *anatomical semantics* still gives it the best median mid (11.6) and a + competitive MI (0.511), but the full-volume registration is a larger problem + for the local ASGD. It also drifts grossly on `m6_2BBN5XA4` (mid −1.7, the + transform flips to the far side). +- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.6 s + (warm workers), SimpleITK ~9.7 s, ANTs ~36.7 s (CPU-bound, full volumes), + **SynthMorph ~25 s** (GPU: a single pretrained-model forward pass per case), + **IMPACT ~834 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**: - medians 0.645 (SimpleITK) / 0.567 (Elastix) / 0.469 (IMPACT) / 0.436 (ANTs) / - 0.241 (SynthMorph), - and it diverges sharply on the tilted cases — ANTs `m6_24V62QQM_20230223` - drops to ~0.06 nats (rotation ~47 mm off) while SimpleITK holds 0.703. A high - MI means the CT skull/tissue intensities fall on the matching MR structures - after the inverse transform. ANTs' and SynthMorph's low median MI is - consistent with their tilted-case convergence failures. + medians 0.651 (SimpleITK) / 0.641 (Elastix) / 0.511 (IMPACT) / 0.652 (ANTs) / + 0.238 (SynthMorph). + The intensity-metric tools cluster at ~0.64 (a genuine alignment signal), + while SynthMorph's 0.238 reflects its tilt failures. A high MI means the CT + skull/tissue intensities fall on the matching MR structures after the inverse + transform. - **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 @@ -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 `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) + method that recovers the gross offset on its own). Its low median MI (0.238) reflects the tilt failures. On same-modality MRI it would likely rank much higher; this CT benchmark is an out-of-distribution stress test for it. It is deterministic (one fixed pretrained model, no random subsampling), so unlike @@ -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 TorchScript **TotalSegmentator** feature extractor. It is GPU-based and 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). +**fixed = T1c, moving = CT** (the reverse of the others): the CT FOV +(~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 `PipelineUnavailable` with the exact command): @@ -330,37 +326,36 @@ repo for the full list and FOV/layer tables. ### How the adapter drives the image -The image's `CMD.sh` only runs `elastix -f -m -p -out` (no `-t0`), so the -adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`) -to add the shared centered start: +The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image +binary directly (`/lib/elastix-install/bin/elastix`): -1. **Head box** — resamples fixed=T1c and moving=CT into the same CT-space head - box (physical frame) the other pipelines use, so MI/IMPACT is computed on - head tissue, not neck/air. -2. **Shared init via `-t0`** — the `centroid_offset` translation is written to a - `TranslationTransform` parameter file and passed with `-t0`; with - `HowToCombineTransforms "Compose"` the final transform is - `T = T_euler(T_init(p))`. (Elastix 5.x removed in-file multi-`Transform` - chains and `UseScalingFactors`, so the initial transform must be a separate - `-t0` file.) +1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box). + T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too + many CT samples map outside the T1c buffer and the metric aborts. +2. **Auto init (no `-t0`)** — the parameter map sets + `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so + elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the + optimizer scales. 3. **Parameter map** — rigid `EulerTransform` (6-DOF), `Impact` metric in **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), ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid (`--impact-gpu -1` forces CPU). 4. **Read-back** — the result `TransformParameters.0.itk.txt` is an - `Euler3DTransform` in the box/world physical frame; the adapter composes it - with the `-t0` translation and runs the usual `fix_orientation` gate. + `Euler3DTransform` in the world physical frame (the + `AutomaticTransformInitialization` is folded in by elastix); the adapter + runs it through the usual `fix_orientation` gate. ## Notes / caveats -- The benchmark compares each toolkit in its **native configuration**: the - local optimizers that cannot recover the gross offset from identity - (`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`, - `elastix` and `synthmorph` start from identity / their own global search. - Production 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 +- The benchmark compares each toolkit in its **native configuration**: all + register on the full volumes (no head box) and use each toolkit's native + automatic transform initialization (no manual `centroid_offset`) — elastix / + impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's + center-of-mass alignment, simpleitk identity, synthmorph none. Production + 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. - `synthmorph` is deterministic (one fixed pretrained model, no random subsampling), so its per-case numbers are stable run to run — unlike ANTs / @@ -385,8 +380,8 @@ to add the shared centered start: 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). + - **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK + (~10 s) / ANTs (~37 s) and IMPACT (~834 s). - **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python API — no CLI binary involved): - 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 random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x names and is not compatible with elastix ≤ 4.x parameter files. - - No initial transform is passed (identity start): the ASGD + 3-resolution - Mattes-MI pyramid recovers the gross inter-session offset on its own, so - the shared `centroid_offset` start is not required. An **output directory** - is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter - uses a scratch dir under `data/m6reg/elastix/itkelastix/`. + - `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the + start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the + optimizer scales. An **output directory** is still set + (`SetOutputDirectory`) for the C++ filter to run; the adapter uses a + scratch dir under `data/m6reg/elastix/itkelastix/`. - `WriteResultImage` must be the string `"false"` (Python `False`); the default (absent) makes the filter raise "Could not read result image" after a *successful* registration. @@ -412,28 +407,32 @@ to add the shared centered start: 500/500/250 over 3 resolutions). - **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0, 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 - 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. + 500-iter / 2000-sample budget is ~14 min/case on full volumes (vs 2–37 s + 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 - 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). + follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251 + mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the + fixed image too many CT samples map outside the T1c buffer and the metric + 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 - minima than Mattes-MI: the pipeline nails the hardest repositioned case - (`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some - others. A multi-start or higher-iteration budget would likely raise its - acceptance; that is out of scope for the single-init head-to-head. - - **Stochastic.** Jacobian mode resamples patches each iteration, so results - vary run to run (see Stochasticity). - - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image - bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors` - are gone — the shared centered start is a separate `-t0` translation file, - and the result is read back as an `Euler3DTransform` (`.itk.txt`) and - composed with `-t0`. The image's own `CMD.sh` doesn't pass `-t0`, so the - adapter invokes `/lib/elastix-install/bin/elastix` directly. + minima than Mattes-MI: on full volumes it no longer nails the hardest + repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box + used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration + budget would likely raise its acceptance; that is out of scope for the + single-init head-to-head. + - **Stochastic.** Jacobian mode resamples patches each iteration, so results + vary run to run (see Stochasticity). + - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image + bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors` + are gone. The result is read back as an `Euler3DTransform` (`.itk.txt`). + The image's own `CMD.sh` is a generic launcher, so the adapter invokes + `/lib/elastix-install/bin/elastix` directly. - **Models are not in the image.** The TorchScript model(s) must be present under the mounted `/Data/Models` (the repo caches them in `data/impact_models/`, bind-mounted read-only). diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index 152297a..243f722 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -23,34 +23,35 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)`` Pipelines (listed alphabetically) --------------------------------- -All pipelines use the shared convention **fixed = CT, moving = T1c**, except -IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its -TotalSegmentator metric -- see ``register_impact``). A pipeline that is a -local optimizer and cannot recover the gross inter-session offset from -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. +All pipelines register on the **full volumes** (no head box) and use each +toolkit's **native automatic transform initialization** (no manual +``centroid_offset``): elastix / impact use elastix's +``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its +default center-of-mass alignment, simpleitk starts from identity, and +synthmorph needs none (global deep method). All use **fixed = CT, moving = +T1c** except IMPACT, which keeps **fixed = T1c, moving = CT** (swapping breaks +its TotalSegmentator metric -- see ``register_impact``). - ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). - Centroid pre-shift of the moving (T1c) image is *required* (ANTs - Rigid-from-identity collapses on the gross offset). + Default ``initial_transform=None`` computes a translation aligning the image + centers of mass. - ``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). + ``AutomaticTransformInitialization`` (CenterOfGravity) start. - ``register_impact`` : IMPACT (semantic-similarity metric) inside the ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF, - 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). + GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix + ``AutomaticTransformInitialization`` (CenterOfGravity) provides the start. 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 +- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on the + full volumes. 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. + method -- no init. Needs Docker + the pulled image. A toolkit that is not installed (or, for the Docker-based pipelines, missing 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])) -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) # --------------------------------------------------------------------------- # 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 identity start, capped at ``iters_full`` iterations at full resolution. This is the *standard* SimpleITK algorithm and matches the budget the other toolkits get in the benchmark (one rigid MI registration), so the - pipelines are compared on equal footing. The pyramid recovers the gross - inter-session offset, so no centroid init is used. + pipelines are compared on equal footing. No centroid init is used. mode="production": the exact multi-start tilt-grid ``register_ct_to_t1c`` from scripts/02_build_m6_dataset.py (robust to the ~30 deg inter-session 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: pass ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32) - lo, hi = _ct_head_box_mm(ct) - if lo is None: - raise PipelineUnavailable("simpleitk: 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)) - fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) - moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) - # Identity start: the multi-resolution Mattes-MI pyramid recovers the gross - # 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)) + # full volumes, fixed = CT, moving = T1c (no head box) + fixed, moving = ct, t1c + shape = np.array(ct.GetSize()) + # Identity start (no centroid init); the multi-resolution Mattes-MI pyramid + # recovers the gross inter-session offset. The Euler center is the fixed + # (CT) geometric center. + init = sitk.Euler3DTransform() + init.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint( + ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2, + (fixed.GetSize()[2] - 1) // 2)))) init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) import time t0 = time.time() @@ -250,9 +227,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150): meta = { "pipeline": "simpleitk", "mode": "benchmark", "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}, " - f"{iters_full} iters/level)"), + f"{iters_full} iters/level, identity start)"), "convention": "T1c->CT (benchmark budget)", "auto_inverted": False, "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm, "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): - """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`` + ``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the 5.x replacement for ``Rigid3DTransform``) + Mattes MI - (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common - CT-space head box, from the identity start (the ASGD + multi-resolution - pyramid recovers the gross inter-session offset, so no centroid init is - required). Elastix 5.x component names: sampler ``Random`` (not - ``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``). An output - directory is required for the 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. + (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the full + CT / T1c volumes (no head box). No manual centroid init: elastix's + ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity + centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales. + Elastix 5.x + component names: sampler ``Random`` (not ``RandomSampler``), metric key + ``Metric`` (not ``ImageMetric``). An output directory is required for the + 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() - lo, hi = _ct_head_box_mm(ct) - if lo is None: - raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)") - # 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 + # full volumes, fixed = CT, moving = T1c (no head box) + fixed = ct + moving = t1c params = { "ImageDimension": 3, "FixedInternalImagePixelType": "float", @@ -348,6 +317,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8): "FixedCenteredTransform": True, "TransformParameterReadoutFormat": 1, "HowToCombineTransforms": "Compose", + "AutomaticTransformInitialization": True, + "AutomaticTransformInitializationMethod": "CenterOfGravity", + "AutomaticScalesEstimation": True, "Optimizer": "AdaptiveStochasticGradientDescent", "MaximumNumberOfIterations": [500, 500, 250], "Metric": "AdvancedMattesMutualInformation", @@ -388,8 +360,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8): version = "ITK-Elastix (python API)" meta = { "pipeline": "elastix", "toolkit": version, - "metric": "MattesMutualInformation (head box, CT space)", - "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)", + "metric": "MattesMutualInformation (full volumes)", + "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)", "convention": "T1c->CT (verified)", "auto_inverted": inverted, "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, aff_shrink_factors=None, aff_smoothing_sigmas=None, - aff_metric="mattes", use_head_box=True, center_init=True): - """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c). + aff_metric="mattes"): + """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c), full volumes. antspyx 0.6.3: ``ants.registration(fixed, moving, type_of_transform='Rigid', aff_metric='mattes', aff_iterations=..., aff_shrink_factors=..., aff_smoothing_sigmas=...)``. ``aff_metric='mattes'`` (the default) gives a multi-resolution Mattes-MI rigid registration -- the cross-modality metric - 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 - first resampled into the same CT-space head box. ANTs Rigid-from-identity - 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 + this CT/MR pair needs. Registers on the full CT / T1c volumes (no head box). + No manual centroid init: antspyx's default ``initial_transform=None`` + computes a translation that aligns the image centers of mass. Returns the transform that warps CT into the T1c grid (orientation verified by :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_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) - if use_head_box: - lo, hi = _ct_head_box_mm(ct) - 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 + # full volumes, fixed = CT, moving = T1c (no head box) + fx, mv = ct, t1c import time t0 = time.time() - # Centered start (REQUIRED for ANTs): antspyx cannot use a - # TranslationTransform as initial_transform (the ANTs CLI rejects it), and - # 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) + # No manual centroid init: antspyx's default initial_transform=None computes + # a translation aligning the image centers of mass. af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)) am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) 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) dt = time.time() - t0 # 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 - # the pre-shifted T1c (translation by -off), then invert the ANTs transform. + # (T1c). The contract is T1c -> CT, i.e. its inverse; fix_orientation + # verifies the direction against the skull ring. T_ants = sitk.ReadTransform(r["fwdtransforms"][0]) - if off is not None: - 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() + cand = T_ants.GetInverse() T, inverted = fix_orientation(ct, t1c, cand) meta = { "pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", - "metric": f"ANTs {aff_metric} (multi-resolution)", - "transform": f"ANTs {transform_type} (head box)", "convention": "T1c->CT (verified)", + "metric": f"ANTs {aff_metric} (multi-resolution, full volumes)", + "transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)", "auto_inverted": inverted, "aff_iterations": list(aff_iterations), "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)", '(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")', '(TransformBase "ElastixTransformBase")', '(TransformParameterReadoutFormat "AllParametersOneLine")', '(FixedCenteredTransform "true")', '(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)", '(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 ``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* one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run 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 built for. - To stay on the benchmark's equal footing the pipeline: - - registers on the same CT-space head box (both images resampled into it), - - starts from the *shared* centroid-aligned translation (passed to Elastix - as a fixed initial transform via ``-t0``; with ``HowToCombineTransforms - Compose`` the final transform = ``T_main(T_init(p))``), - - uses a rigid 6-DOF ``EulerTransform`` (the 5.x name for the classic - 3D rigid), ASGD, a 4-level scale pyramid (no image resampling). + The pipeline registers on the full CT / T1c volumes (no head box) with no + manual centroid init: elastix's ``AutomaticTransformInitialization`` + (CenterOfGravity) aligns the intensity centroids. It 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' recommended Jacobian budget. It is *much* slower than the intensity-metric pipelines (minutes per case on GPU vs seconds) because every iteration runs the TotalSegmentator feature extractor on the sampled patches with backpropagation. Lower ``iters``/``samples`` to trade accuracy for speed. The result is read back from ``TransformParameters.0.itk.txt`` - (an ``Euler3DTransform`` in the box/world physical frame) and composed with - the ``-t0`` translation, then checked by :func:`fix_orientation`. + (an ``Euler3DTransform`` in the world physical frame) and checked by + :func:`fix_orientation`. Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, the ``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}") model_host, model_rel = _impact_model_host(model) - lo, hi = _ct_head_box_mm(ct) - if lo is None: - raise PipelineUnavailable("impact: no CT head box found (corrupt CT?)") - - # resample both into the shared CT-space head box (mirrors the other pipelines) - 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(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])) + # full volumes, fixed = T1c, moving = CT (no head box, no manual init -- + # elastix's AutomaticTransformInitialization aligns the image centers). + # 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 + # map outside the T1c buffer and the IMPACT metric aborts. + fixed = t1c + moving = ct # GPU: auto-spread across the host GPUs by pid, or forced value 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) with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: f.write(_impact_param_map(os.path.join("/Data/Models", model_rel), - gpu, iters, samples, voxel_sizes, - 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") + gpu, iters, samples, voxel_sizes, + rescale=schedule, smoothing=smoothe)) cmd = ["docker", "run", "--rm"] 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, "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha", "-p", "/Data/ParameterMap.txt", - "-t0", "/Data/InitialTransform.txt", "-out", "/Out/", "-threads", str(threads), ] 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:]) raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}") - # the optimized Euler (box/world frame) composed with the fixed -t0 - # 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. + # the optimized Euler (world frame) is the full combination transform + # (AutomaticTransformInitialization is folded in by elastix); it maps + # fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT); + # fix_orientation resolves the direction. T_main = sitk.ReadTransform(itk_tf) - T_init = sitk.TranslationTransform(3) - T_init.SetOffset(init) - cand = sitk.CompositeTransform([T_init, T_main]) - T, inverted = fix_orientation(ct, t1c, cand) + T, inverted = fix_orientation(ct, t1c, T_main) finally: if not keep_scratch: 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", "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", "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'})"), "convention": "T1c->CT (verified)", "auto_inverted": inverted, "fixed": "T1c", "moving": "CT",