diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index a2dd5b5..6717c75 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/REGISTRATION_BENCHMARK.md @@ -22,23 +22,34 @@ sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0) # GTV onto the MR grid ``` places the CT / GTV concentric with the MR head. In `Resample`'s -destination→source convention this `T` maps **T1c → CT**. 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). All single-init tools start from the -**same centroid-aligned start** (`centroid_offset`), so the comparison is on -the optimizer / metric / transform handling — not on who finds the ~30–60 mm -gross offset first. +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 +each toolkit's native optimizer / metric / global-search capability. ## The five pipelines -| # | pipeline | toolkit | metric | transform | status in this env | -|---|----------|---------|--------|-----------|--------------------| -| 1 | `elastix` | ITK-Elastix 0.25.4 (pip `itk` + `itk-elastix`, elastix 5.3.1) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x replacement for `Rigid3DTransform`), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) | -| 2 | `simpleitk` | SimpleITK 2.5.6 | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | -| 3 | `antspyx` | ANTsPy / ANTs 0.6.3 | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | -| 4 | `synthmorph` | Docker `freesurfer/synthmorph` (TensorFlow, pretrained) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix (centroid init) | ✅ ran (main env, Docker + GPU) | -| 5 | `impact` | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT) | **IMPACT semantic** (TotalSegmentator `TS/M730.pt` features, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales (no image resampling) | ✅ ran (main env, Docker + GPU) | +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. + +| # | pipeline | fixed / moving | centroid init | metric | transform | status in this env | +|---|----------|----------------|---------------|--------|-----------|--------------------| +| 1 | `antspyx` | CT / T1c | **yes** (required) | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | +| 2 | `elastix` | CT / T1c | no (identity) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) | +| 3 | `impact` | **T1c / CT** | **yes** | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) | +| 4 | `simpleitk` | CT / T1c | no (identity) | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | +| 5 | `synthmorph` | CT / T1c | no (global deep) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) | `simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid, centered start — the default, comparable to the other tools) and @@ -87,32 +98,31 @@ 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) | |---|---|---|---|---|---|---|---|---| -| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.6 | 0.93 | **0.643 / 0.630** | -| **antspyx** | ANTsPy 0.6.3 | 9 | **0/10 (0.00)** | 14.9 / 16.5 | 18.4 / 22.7 | 14.5 | 0.93 | **0.598 / 0.529** | -| **elastix** | ITK-Elastix 0.25.4 | 9 | **5/10 (0.50)** | 15.3 / 16.1 | 9.9 / 13.8 | 3.2 | 0.93 | **0.562 / 0.549** | -| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **1/10 (0.10)** | 19.8 / 25.1 | 23.0 / 24.1 | 28.2 | 0.93 | **0.215 / 0.309** | -| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 12.4 / 19.0 | 15.1 / 16.7 | **562.3** | 0.93 | **0.475 / 0.448** | +| **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 and impact are stochastic (random metric/patch subsampling), so their -accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.4 +accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.6 min/case on GPU, orders of magnitude above the intensity-metric tools. -synthmorph is deterministic at fixed threads; its low acceptance is a real -cross-modality/tilt weakness, not run-to-run variance.)* +elastix, simpleitk and synthmorph are deterministic at fixed threads.)* Per-case (mid / std, `OK` = accepted): -| case | simpleitk | antspyx | elastix | synthmorph | impact | +| case | antspyx | elastix | impact | simpleitk | synthmorph | |---|---|---|---|---|---| | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | -| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 24.4 / 27.7 | 26.4 / 30.3 | 20.0 / 25.2 | **8.4 / 12.6** | -| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.4 / 24.8 | 21.1 / 23.3 | 19.8 / 23.2 | 21.4 / 24.1 | -| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.5 / 18.4 | 10.6 / 8.2 **OK** | 12.0 / 23.0 | 10.0 / 8.3 **OK** | -| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 14.9 / 15.5 | 11.0 / 6.1 **OK** | 13.0 / 12.6 | 11.1 / 8.8 **OK** | -| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.3 / 17.5 | 15.3 / 14.0 | 53.8 / 20.2 | 54.4 / 15.1 | -| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 17.7 / 37.6 | 11.4 / 3.2 **OK** | 11.3 / 5.0 **OK** | 23.9 / 22.6 | -| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** | 49.3 / 58.9 | 10.1 / 5.9 **OK** | -| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.1 / 14.3 | 12.0 / 6.8 **OK** | 18.7 / 19.3 | 12.4 / 22.3 | -| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 17.2 / 35.4 | 18.9 / 22.7 | 27.6 / 29.1 | 19.6 / 30.3 | +| 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 | ### Interpretation @@ -121,57 +131,61 @@ Per-case (mid / std, `OK` = accepted): `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 centered - start can lock onto). This is exactly why production `02_build_m6_dataset.py` - uses a **multi-start tilt grid**. -- **Elastix is the most accurate of the standard single-init configs** - (5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 18.4) and - the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s - SimpleITK / 14.5 s ANTs): ASGD's adaptive step sizes climb out of small - misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223` - (std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on - (std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both - others accept — different optimizers settle in different local minima on the - hard tilted cases. -- **ANTs' Rigid is the weakest of the intensity-metric tools** (0/10 in this - draw; it is stochastic, so a prior draw gave 1/10): it frequently locks onto - the tilted local minimum. Note `m6_26TYQI3C`: ANTs nailed the position (mid - 17.7) but missed the orientation (std 37.6) — a textbook tilted local minimum. + (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, std 28–30 mm) it gets **mid 8.4 / std 12.6** — the best - translation of any pipeline and within 0.6 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 - 4 of 9 cases (std 15–30 mm), and one (`m6_24V62QQM_20240820`) even drifts the - translation (mid 54 mm). Net 3/10 accepted — between ANTs and SimpleITK — at - ~9.4 min/case, two orders of magnitude slower than the intensity tools. -- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~3.2 s - (warm workers; a worker's first case pays a ~10 s warm-up), SimpleITK ~8.6 s, - ANTs ~14.5 s (CPU-bound), **SynthMorph ~28 s** (GPU: a single pretrained-model - forward pass per case), **IMPACT ~562 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.643 (SimpleITK) / 0.598 (ANTs) / 0.562 (Elastix) / 0.475 (IMPACT) / - 0.215 (SynthMorph), - and it diverges sharply on the tilted cases — ANTs `m6_26TYQI3C` drops to - ~0.1 nats (rotation ~38 mm off) while SimpleITK holds 0.670. A high MI means - the CT skull/tissue intensities fall on the matching MR structures after the - inverse transform. IMPACT's lower median MI is consistent with its less -reliable single-init convergence. - - **SynthMorph (deep, pretrained) is the weakest of the five on these - cross-modality CT pairs.** Its rigid model is trained on brain *MRI* - (T1/T2, anatomy-aware but intensity-domain MRI), so it recovers gross - translation well (mid 12–20 mm on most cases) but is repeatedly defeated by - the ~30° inter-session tilt: 6 of 9 cases have `std > 20 mm` (rotation error) - and two drift grossly (`m6_2BBN5XA4` mid 49 / std 59, `m6_24V62QQM_20240820` - mid 54). It accepts only the clean `m6_26TYQI3C` (mid 11.3 / std 5.0, MI - 0.64 — where it is competitive with the intensity tools). Its low median MI - (0.215) reflects the tilt failures. On same-modality MRI it would likely - rank much higher; this CT benchmark is an out-of-distribution stress test for - it. It is deterministic (one fixed pretrained model, no random subsampling), - so unlike ANTs/IMPACT the per-case numbers are stable run to run. + 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 + 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. +- **SynthMorph (deep, pretrained) is near-weakest of the five on these + cross-modality CT pairs.** Its rigid model is trained on brain *MRI* + (T1/T2, anatomy-aware but intensity-domain MRI), so it recovers gross + translation well (mid 10–20 mm on most cases) but is repeatedly defeated by + the ~30° inter-session tilt: 7 of 9 cases have `std > 15 mm` (rotation error) + and two drift grossly (`m6_2BBN5XA4` mid 52 / std 53, `m6_24V62QQM_20240820` + mid 43). It accepts the clean `m6_26TYQI3C` (mid 11.8 / std 2.9, MI 0.60) and + `m6_24V62QQM_20231017` (mid 11.5 / std 8.4) — where it is competitive with the + intensity tools. Notably it needs **no centroid init** (it is a global deep + method that recovers the gross offset on its own). Its low median MI (0.241) + reflects the tilt failures. On same-modality MRI it would likely rank much + higher; this CT benchmark is an out-of-distribution stress test for it. It is + deterministic (one fixed pretrained model, no random subsampling), so unlike + ANTs/IMPACT the per-case numbers are stable run to run. - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its per-case numbers carry run-to-run variance (this table is one draw); SimpleITK's MI (and registration) is deterministic at fixed thread count @@ -270,29 +284,27 @@ No extra model assets: the rigid weights ship inside the image ### How the adapter drives the image -1. **Shared init** — the `centroid_offset` (CT→T1c head translation) is built - *inside* the container as an initial LTA (`-i`), because surfa requires the - init transform's embedded src/dst geometries to match the input images - exactly. The container also writes the LTA as the host uid so the scratch - removes cleanly. -2. **Rigid matrix model** — `register -m rigid -g` on GPU when one is visible - (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces CPU). -3. **Read-back** — the `-t` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in shared world - space) is parsed into a SimpleITK `AffineTransform` and run through the +1. **Rigid matrix model, no init** — `register -m rigid` on GPU when one is + visible (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces + CPU). SynthMorph is a global deep method, so it recovers the gross + inter-session offset on its own and needs **no** initial transform. +2. **Read-back** — the `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in shared world + space; with fixed=CT/moving=T1c the `-T` matrix is the T1c→CT contract + transform) is parsed into a SimpleITK `AffineTransform` and run through the usual `fix_orientation` gate. -4. **Image-build workarounds** — this image's `-i` code path has a - float64/float32 bug (the dense 256³ mesh is built as `double` and the final - matrix product mismatches dtypes), so the adapter injects a small wrapper - that casts the transform operands to float32, and sets - `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent workers can share a GPU. +3. **GPU sharing** — the container is run as the host uid (clean scratch + removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so + concurrent benchmark workers can share a GPU. ## IMPACT pipeline (Docker) — one-time setup The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact) image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121) *semantic* similarity metric compiled in as a native Elastix metric, plus a -TorchScript **TotalSegmentator** feature extractor. It is the only pipeline that -is GPU-based and the only one that is Docker-based. +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). Prerequisites (check with `RP.availability()`; a missing piece raises `PipelineUnavailable` with the exact command): @@ -343,38 +355,38 @@ to add the shared centered start: ## Notes / caveats -- The benchmark compares **standard single-init** configurations (centered - start). 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. +- 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 + head-to-head but is the config the dataset actually ships with. - `synthmorph` is deterministic (one fixed pretrained model, no random subsampling), so its per-case numbers are stable run to run — unlike ANTs / - IMPACT. Its low acceptance (1/10) is a genuine cross-modality / tilt + IMPACT. Its low acceptance (2/10) is a genuine cross-modality / tilt weakness, not variance. - `fix_orientation` (elastix/ants/synthmorph/impact) chooses between a transform and its inverse by skull-ring concentricity — a correctness gate that makes - the result independent of each toolkit's transform-direction convention. It - rarely fired: elastix/ants/impact returned the T1c→CT direction as expected - (`auto_inverted: False`) and synthmorph did so on 8/9 cases (its one -`True` was the grossly-drifted `m6_2BBN5XA4`, where the gate flipped a bad - pose). + the result independent of each toolkit's transform-direction convention. With + the fixed=CT/moving=T1c swap, `elastix`'s combination transform is now the + inverse of the contract, so it reports `auto_inverted: True` (the gate + corrects it); antspyx/synthmorph/impact return the contract direction + directly (`auto_inverted: False`). - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU): - - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain - *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly - defeated by the ~30° inter-session tilt (6/9 cases have `std > 20 mm`). - It is a legitimate data point for how far a pretrained MRI deep model - generalizes out-of-domain to CT — not a bug. - - **Shared init is built in-container.** surfa rejects an `-i` transform - whose src/dst geometries don't match the inputs (tol 1e-3), so the - adapter has the container write the centroid-translation LTA from the - mounted images themselves. - - **Image-build dtype/GPU workarounds.** This image's `-i` path builds the - dense 256³ mesh as float64 (2× memory → GPU OOM under contention) and - mismatches dtypes in the final matrix product; the injected wrapper casts - transform operands to float32, and `TF_FORCE_GPU_ALLOW_GROWTH` lets - concurrent workers share a GPU. - - **~28 s/case** (GPU: one pretrained forward pass), between ANTs (~15 s) - and IMPACT (~562 s). + - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain + *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly + defeated by the ~30° inter-session tilt (7/9 cases have `std > 15 mm`). + It is a legitimate data point for how far a pretrained MRI deep model + generalizes out-of-domain to CT — not a bug. + - **No centroid init needed.** Being a global deep method it recovers the + gross offset on its own; the `-i` initial-transform path is not used + (which also sidesteps this image's float64/float32 bug on the `-i` code + path). + - **GPU sharing.** The container runs as the host uid with + `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU. + - **~26 s/case** (GPU: one pretrained forward pass), between ANTs (~11 s) + and IMPACT (~576 s). - **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python API — no CLI binary involved): - Elastix 5.x renamed components vs the classic parameter files: the rigid @@ -382,32 +394,38 @@ 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. - - The initial transform (the shared `centroid_offset` start) is passed as an - `InitialTransformParameterObject` (a `TranslationTransform` with - `NumberOfParameters 3`); the C++ filter only wires it up when an - **output directory is set** (`SetOutputDirectory`), so the adapter uses a - scratch dir under `data/m6reg/elastix/itkelastix/`. + - 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/`. - `WriteResultImage` must be the string `"false"` (Python `False`); the default (absent) makes the filter raise "Could not read result image" after a *successful* registration. - - The result transform is recovered by point-mapping + - With fixed=CT/moving=T1c the combination transform is CT→T1c (the inverse + of the contract), so `fix_orientation` reports `auto_inverted: True` and + inverts it. The result is recovered by point-mapping (`GetCombinationTransform().TransformPoint`) into an affine — independent of elastix's Euler-angle parameterization — and wrapped in a SimpleITK `AffineTransform(3)`. - - Registration wall time is ~1.8 s/case at full budget (ASGD 500/500/250 - over 3 resolutions); each worker's first case pays a ~10 s warm-up, so the - median wall is 3.2 s and the mean is inflated by warm-up. + - Registration wall time is ~2.4 s/case median at full budget (ASGD + 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 - extractor on the sampled patches with backpropagation, so the recommended - 500-iter / 2000-sample budget is ~9 min/case (vs 3–15 s for the intensity - tools). `--impact-iters` / `--impact-samples` trade accuracy for speed. - - **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.4 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. +- **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. + - **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT + follow the shared fixed=CT/moving=T1c convention: swapping its roles + breaks the TotalSegmentator metric (elastix aborts with "Error in metric" + or converges to a wrong rotation). The returned transform still satisfies + the contract (T1c→CT). + - **Rugged single-init convergence.** The semantic landscape has more local + 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 diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index e6163c5..152297a 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -21,22 +21,36 @@ destination->source convention this transform therefore maps **T1c -> CT** coordinates. This is exactly the convention the production SimpleITK pipeline returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``). -Pipelines ---------- -- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a - CT-space head box, multi-start tilt grid. Reuses the *proven* production - ``register_ct_to_t1c`` so the baseline is bit-identical to 02_build_m6_dataset. -- ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` + ``itk-elastix``, - elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD). -- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid). -- ``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. Needs Docker - + the pulled image. -- ``register_impact`` : IMPACT (semantic-similarity metric) inside the +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. + +- ``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). +- ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` + + ``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD). + Identity start (the ASGD + multi-resolution pyramid recovers the offset). +- ``register_impact`` : IMPACT (semantic-similarity metric) inside the ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript - TotalSegmentator feature extractor as a native Elastix metric). Rigid - 6-DOF, GPU (Jacobian mode). Needs Docker + the image + a model asset. + 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). + Needs Docker + the image + a model asset. +- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a + CT-space head box. Reuses the *proven* production ``register_ct_to_t1c`` in + ``mode="production"``; the benchmark mode is a single identity-start + multi-resolution pyramid run. +- ``register_synthmorph`` : FreeSurfer SynthMorph (TensorFlow, pretrained + anatomy-aware contrast-invariant model) inside the ``freesurfer/synthmorph`` + Docker image. Rigid 6-DOF matrix model, GPU when available. Global deep + method -- no centroid init. Needs Docker + the pulled image. A toolkit that is not installed (or, for the Docker-based pipelines, missing its image) @@ -162,22 +176,6 @@ def _ct_head_box_mm(ct): return box_lo, box_hi -def _ct_head_box_resample(ct, t1c): - """Resample the CT into the shared CT-space head box (mirrors the other - pipelines), or None when no head box is found. Used only to get the CT - *head* tissue centroid for the shared centered start.""" - lo, hi = _ct_head_box_mm(ct) - if lo is None: - return None - 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)) - return sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) - - # --------------------------------------------------------------------------- # # 2. SimpleITK (production baseline -- reused verbatim) # --------------------------------------------------------------------------- # @@ -186,10 +184,11 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150): """SimpleITK rigid CT->T1c, Mattes MI on a CT-space head box. mode="benchmark" (default): a single multi-resolution pyramid run from the - centered start, capped at ``iters_full`` iterations at full resolution. + identity start, capped at ``iters_full`` iterations at full resolution. This is the *standard* SimpleITK algorithm and matches the budget the other toolkits get in the benchmark (one rigid MI registration), so the - pipelines are compared on equal footing. + pipelines are compared on equal footing. The pyramid recovers the gross + inter-session offset, so no centroid init is used. mode="production": the exact multi-start tilt-grid ``register_ct_to_t1c`` from scripts/02_build_m6_dataset.py (robust to the ~30 deg inter-session tilt but much slower -- not the default for the benchmark). @@ -226,12 +225,11 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150): 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) - # centered start: align CT head centroid with T1c head centroid (shared init - # with the other single-init toolkits via centroid_offset) - off = centroid_offset(ct, t1c, ct_box=fixed) or (0.0, 0.0, 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)) - init.SetParameters((0.0, 0.0, 0.0) + off) + init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) import time t0 = time.time() r = sitk.ImageRegistrationMethod() @@ -313,21 +311,20 @@ def _itk_transform_to_sitk_rigid(itk_tf): def register_elastix(ct, t1c, workdir=None, threads=8): - """ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI. + """ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), Mattes MI. Runs in-process through the ITK-Elastix python API (pip ``itk`` + ``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 *shared* centroid-aligned start - (``centroid_offset`` as an initial TranslationTransform, same init the - SimpleITK/ANTs benchmark modes use). Elastix 5.x component names: sampler - ``Random`` (not ``RandomSampler``), metric key ``Metric`` (not - ``ImageMetric``). An output directory is required so the initial - transform parameter object is honored (the C++ filter only wires it up - when ``SetOutputDirectory`` is set), and ``WriteResultImage`` must be the - string ``"false"`` (a Python bool False) or the filter raises "Could not - read result image" after a successful registration.""" + 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. + """ itk = _elastix_api() lo, hi = _ct_head_box_mm(ct) if lo is None: @@ -340,18 +337,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8): 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: elastix transforms map fixed(T1c) -> moving(CT), - # i.e. the inverse direction of centroid_offset, hence the negation. - off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0) - init_tp0 = { - "Transform": "TranslationTransform", - "TransformBase": "ElastixTransformBase", - "TransformParameterReadoutFormat": 1, - "NumberOfParameters": 3, - "TransformParameters": [-float(off[0]), -float(off[1]), -float(off[2])], - } + 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 = { "ImageDimension": 3, "FixedInternalImagePixelType": "float", @@ -378,7 +365,6 @@ def register_elastix(ct, t1c, workdir=None, threads=8): reg.SetFixedImage(_sitk_to_itk(fixed, itk)) reg.SetMovingImage(_sitk_to_itk(moving, itk)) reg.SetParameterObject(itk.dict_to_parameter_object(params)) - reg.SetInitialTransformParameterObject(itk.dict_to_parameter_object(init_tp0)) reg.SetOutputDirectory(outdir) reg.SetNumberOfThreads(threads) reg.LogToConsoleOff() @@ -390,8 +376,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8): tail = "\n".join(buf.getvalue().splitlines()[-10:]) raise RuntimeError(f"elastix registration failed: {e}\n{tail}") from e dt = time.time() - t0 - # The combination transform maps fixed(T1c) -> moving(CT) == the contract - # (T1c->CT); fix_orientation still verifies it against the skull ring. + # The combination transform maps fixed(CT) -> moving(T1c) (elastix + # convention), i.e. the inverse of the contract (T1c->CT); + # fix_orientation resolves the direction against the skull ring. cand = _itk_transform_to_sitk_rigid(reg.GetCombinationTransform()) T, inverted = fix_orientation(ct, t1c, cand) try: @@ -402,9 +389,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8): meta = { "pipeline": "elastix", "toolkit": version, "metric": "MattesMutualInformation (head box, CT space)", - "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, centroid init)", + "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)", "convention": "T1c->CT (verified)", "auto_inverted": inverted, - "fixed": "T1c", "moving": "CT", "seconds_reg": round(dt, 1), + "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1), } return T, meta @@ -416,7 +403,7 @@ 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=T1c, moving=CT). + """ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c). antspyx 0.6.3: ``ants.registration(fixed, moving, type_of_transform='Rigid', aff_metric='mattes', aff_iterations=..., aff_shrink_factors=..., @@ -424,8 +411,11 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None, 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. Returns the transform that - warps CT into the T1c grid (orientation verified by :func:`fix_orientation`). + 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 + transform that warps CT into the T1c grid (orientation verified by + :func:`fix_orientation`). """ try: import ants # antspyx ships the `ants` module @@ -447,25 +437,26 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None, 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(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (fixed) - mv = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (moving) + 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 = t1c, ct + fx, mv = ct, t1c import time t0 = time.time() - # Centered start: antspyx cannot use a TranslationTransform as initial_transform -# (the ANTs CLI rejects it), and Rigid-from-identity does not recover the large -# (~30-60 mm) inter-session offset. So we pre-shift the CT box so its head -# centroid coincides with the T1c head centroid, run Rigid from identity, and -# compose the shift back into the final transform (undoing it to original CT -# coordinates). This gives ANTs the same centered initialization as the other -# single-init toolkits so the pipelines are compared on their optimizer. + # 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=mv) + 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])) + 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)) am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) @@ -474,16 +465,16 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None, aff_shrink_factors=aff_shrink_factors, aff_smoothing_sigmas=aff_smoothing_sigmas) dt = time.time() - t0 - # r['fwdtransforms'][0] is a path to the (single) Rigid transform file, in the - # shared physical space (but on the pre-shifted CT). Compose the shift back so - # the result maps T1c -> original CT; fix_orientation resolves T vs T^-1. + # 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. T_ants = sitk.ReadTransform(r["fwdtransforms"][0]) if off is not None: - tr_undo = sitk.TranslationTransform(3) - tr_undo.SetOffset((-off[0], -off[1], -off[2])) - cand = sitk.CompositeTransform([T_ants, tr_undo]) + 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 + cand = T_ants.GetInverse() T, inverted = fix_orientation(ct, t1c, cand) meta = { "pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", @@ -505,63 +496,9 @@ def sys_pyver(): # --------------------------------------------------------------------------- # SYNTHMORPH_IMAGE = "freesurfer/synthmorph" - -# The image's `-i` (initial-matrix) code path has a float64/float32 bug in this -# build: it computes `net_to_mov = net_to_mov @ tf.linalg.inv(init)` where `init` -# is float64, so `net_to_mov` becomes a TF *float64* tensor. Downstream, (a) -# voxelmorph's compose() only casts non-tensor inputs to float32, so the final -# matrix product fails ("expected a double tensor but is a float tensor"), and -# (b) affine_to_dense_shift() builds the 256^3 dense mesh as float64 (2x memory), -# which OOMs a GPU under contention. This wrapper monkey-patches tf.linalg.inv -# to return float32 (the root cause) and compose() to cast every operand to -# float32 (a safety net), then drives the mri_synthmorph CLI unchanged. It is -# written into the per-case scratch dir and run inside the container. -_SYNTHMORPH_WRAPPER = '''\ -import os -os.environ["TF_FORCE_GPU_ALLOW_GROWTH"] = "true" # let concurrent workers share a GPU -import sys -sys.path.insert(0, "/freesurfer") -import tensorflow as tf -import voxelmorph as vxm - -def _c32(t): - if not tf.is_tensor(t): - return tf.constant(t, dtype=tf.float32) - if t.dtype != tf.float32: - return tf.cast(t, tf.float32) - return t - -_orig_inv = tf.linalg.inv -def inv(x, *a, **k): - return _orig_inv(_c32(x), *a, **k) -tf.linalg.inv = inv - -_orig_compose = vxm.utils.compose -def compose(transforms, *a, **k): - return _orig_compose(tuple(_c32(t) for t in transforms), *a, **k) -vxm.utils.compose = compose - -sys.argv = ["mri_synthmorph"] + sys.argv[1:] -import runpy -runpy.run_path("/freesurfer/mri_synthmorph", run_name="__main__") -''' - -# surfa refuses an init transform whose embedded src/dst geometries do not -# match the input images (tol 1e-3), and reproducing that geometry on the host -# is fragile. So the container builds the -i initial LTA itself from the -# mounted images: a pure CT->T1c world translation (the shared centroid start) -# with source=moving, target=fixed. The offset is passed via SYNTH_INIT_T. -_SYNTHMORPH_BUILD_INIT = '''\ -import os -import numpy as np -import surfa as sf -mov = sf.load_volume("/data/moving.nii.gz") -fix = sf.load_volume("/data/fixed.nii.gz") -t = [float(x) for x in os.environ["SYNTH_INIT_T"].split()] -M = np.eye(4) -M[:3, 3] = t -sf.Affine(M, source=mov, target=fix, space="world").save("/data/init.lta") -''' +# The image entrypoint is ``mri_synthmorph``; the ``register`` subcommand and +# options follow. ``TF_FORCE_GPU_ALLOW_GROWTH`` lets concurrent benchmark +# workers share a GPU instead of each pre-allocating all of it. def _synthmorph_docker_ok(): @@ -597,24 +534,21 @@ def _lta_to_sitk(path): def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None, extent=256, keep_scratch=False): - """SynthMorph CT->T1c (moving=CT, fixed=T1c) in the Docker image + """SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image ``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``, pretrained ``synthmorph.rigid.1.h5`` weights). SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware, acquisition-agnostic deep registration: each volume is normalized to an isotropic 1-mm LIA network space and the (pretrained) network estimates the - transform, so no skull-stripping or intensity preprocessing is needed. To - stay on the benchmark's equal footing the pipeline - - starts from the *shared* centroid-aligned translation (``centroid_offset``, - built inside the container as an initial LTA -- see - ``_SYNTHMORPH_BUILD_INIT`` -- because surfa requires the init transform's - geometries to match the inputs exactly), + transform, so no skull-stripping or intensity preprocessing is needed. Being + a global deep method it recovers the gross inter-session offset on its own, + so **no centroid init is used**. The pipeline - runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible (auto-spread over the host GPUs by pid; ``gpu=-1`` forces CPU), - - reads back the ``-t`` LTA. In RAS world space that matrix is the - T1c->CT transform the unified contract requires (its ``-T`` sibling is - the inverse); :func:`fix_orientation` still verifies the direction. + - reads back the ``-T`` LTA (fixed->moving = the T1c->CT transform the + unified contract requires); :func:`fix_orientation` verifies the + direction against the skull ring. Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI and the ``freesurfer/synthmorph`` image. Only the matrix models (``rigid``, @@ -636,54 +570,42 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None, ngpu = _n_gpus() gpu = (os.getpid() % ngpu) if ngpu else -1 - # shared centered start: the CT->T1c head-centroid translation. As with the - # other single-init pipelines, resample the CT to the head box so the tissue - # centroid is the CT *head* (not head+neck). - box = _ct_head_box_resample(ct, t1c) - off = centroid_offset(ct, t1c, ct_box=box) or (0.0, 0.0, 0.0) - # per-case scratch on a local (non-NFS) tmp fs tmp = tempfile.mkdtemp(prefix="synthmorph_reg_") try: + # fixed = CT, moving = T1c (shared convention for all pipelines) sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32), - os.path.join(tmp, "moving.nii.gz"), True) # CT + os.path.join(tmp, "fixed.nii.gz"), True) # CT sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32), - os.path.join(tmp, "fixed.nii.gz"), True) # T1c - with open(os.path.join(tmp, "build_init.py"), "w") as f: - f.write(_SYNTHMORPH_BUILD_INIT) - with open(os.path.join(tmp, "run_synthmorph.py"), "w") as f: - f.write(_SYNTHMORPH_WRAPPER) + os.path.join(tmp, "moving.nii.gz"), True) # T1c # run as the host uid/gid so the container-written LTA is host-owned # (clean scratch removal) and TF still reaches the GPU. cmd = ["docker", "run", "--rm", - "--user", f"{os.getuid()}:{os.getgid()}"] + "--user", f"{os.getuid()}:{os.getgid()}", + "-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"] if gpu >= 0: cmd += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"] - gpu_flag = " -g" if gpu >= 0 else "" - thread_flag = f" -j {threads}" if threads else "" - cmd += [ - "-e", f"SYNTH_INIT_T={off[0]:.6f} {off[1]:.6f} {off[2]:.6f}", - "-v", f"{tmp}:/data", - "--entrypoint", "sh", - SYNTHMORPH_IMAGE, - "-c", ("python3 /data/build_init.py && " - "python3 /data/run_synthmorph.py register " - f"-m {model}{gpu_flag}{thread_flag} -e {extent} " - "-i /data/init.lta -t /data/t.lta " - "/data/moving.nii.gz /data/fixed.nii.gz"), - ] + cmd += ["-v", f"{tmp}:/data", SYNTHMORPH_IMAGE, + "register", "-m", model] + if gpu >= 0: + cmd.append("-g") + if threads: + cmd += ["-j", str(threads)] + cmd += ["-e", str(extent), + "-T", "/data/f2m.lta", + "/data/moving.nii.gz", "/data/fixed.nii.gz"] t0 = time.time() proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200) dt = time.time() - t0 if proc.returncode != 0: tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) raise RuntimeError(f"synthmorph: docker failed (rc={proc.returncode}):\n{tail}") - if not os.path.exists(os.path.join(tmp, "t.lta")): + if not os.path.exists(os.path.join(tmp, "f2m.lta")): tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) - raise RuntimeError(f"synthmorph: no /data/t.lta produced:\n{tail}") + raise RuntimeError(f"synthmorph: no /data/f2m.lta produced:\n{tail}") - T = _lta_to_sitk(os.path.join(tmp, "t.lta")) + T = _lta_to_sitk(os.path.join(tmp, "f2m.lta")) T, inverted = fix_orientation(ct, t1c, T) finally: if not keep_scratch: @@ -693,11 +615,11 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None, "pipeline": "synthmorph", "toolkit": f"Docker {SYNTHMORPH_IMAGE} (TensorFlow, pretrained)", "metric": "SynthMorph anatomy-aware contrast-invariant (deep, pretrained)", - "transform": (f"SynthMorph {model} (centroid init; " + "transform": (f"SynthMorph {model} (identity start; " f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, " f"extent {extent})"), "convention": "T1c->CT (verified)", "auto_inverted": inverted, - "fixed": "T1c", "moving": "CT", + "fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1), "docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu, "extent": extent, @@ -872,6 +794,13 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image ``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 @@ -922,8 +851,9 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = moving - # shared centered start: elastix maps fixed(T1c)->moving(CT), i.e. the - # inverse of centroid_offset, hence the negation. + # 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])) @@ -1005,7 +935,8 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}") # the optimized Euler (box/world frame) composed with the fixed -t0 - # translation: T = T_euler(T_init(p)) + # translation: T = T_euler(T_init(p)). Both map fixed(CT)->moving(T1c), + # i.e. the inverse of the contract (T1c->CT); fix_orientation resolves it. T_main = sitk.ReadTransform(itk_tf) T_init = sitk.TranslationTransform(3) T_init.SetOffset(init) @@ -1033,11 +964,11 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000, REGISTERS = { - "simpleitk": register_simpleitk, - "elastix": register_elastix, "antspyx": register_antspyx, - "synthmorph": register_synthmorph, + "elastix": register_elastix, "impact": register_impact, + "simpleitk": register_simpleitk, + "synthmorph": register_synthmorph, } ALL_PIPELINES = list(REGISTERS)