refactor(pipelines): remove shared resample helper and clarify init conventions
Drop the _ct_head_box_resample utility that was only needed for centered initialization, since centroid alignment is now scoped to local optimizers (antspyx, impact) rather than all pipelines. Reorder the pipeline table and module docstring alphabetically and add explicit fixed/moving and centroid-init columns so the per-toolkit registration convention is unambiguous.
This commit is contained in:
parent
f80f22dabe
commit
4ad972b00f
2 changed files with 280 additions and 331 deletions
|
|
@ -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).
|
||||
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.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
|
||||
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 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.
|
||||
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`).
|
||||
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.
|
||||
- **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).
|
||||
- **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,30 +394,36 @@ 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
|
||||
- **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
|
||||
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.4 mm — best of all) but lands in bad minima on some
|
||||
(`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some
|
||||
others. A multi-start or higher-iteration budget would likely raise its
|
||||
acceptance; that is out of scope for the single-init head-to-head.
|
||||
- **Stochastic.** Jacobian mode resamples patches each iteration, so results
|
||||
|
|
|
|||
|
|
@ -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).
|
||||
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).
|
||||
- ``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.
|
||||
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)
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue