Compare commits
2 commits
6b85e950e5
...
4ad972b00f
| Author | SHA1 | Date | |
|---|---|---|---|
| 4ad972b00f | |||
| f80f22dabe |
3 changed files with 461 additions and 316 deletions
|
|
@ -22,21 +22,34 @@ sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0) # GTV onto the MR grid
|
|||
```
|
||||
|
||||
places the CT / GTV concentric with the MR head. In `Resample`'s
|
||||
destination→source convention this `T` maps **T1c → CT**. All pipelines register
|
||||
on the same CT-space head box (MI computed on head tissue, not neck/air) and,
|
||||
for the single-init tools, 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 four pipelines
|
||||
## 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 | `deepreg` | Oxford DeepReg (CycleGAN) | CycleGAN + reg. loss (pretrained) | 3D displacement field | **not installed** (see below) |
|
||||
| 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
|
||||
|
|
@ -85,76 +98,95 @@ 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** |
|
||||
| **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** |
|
||||
| deepreg | — | — | — | — | — | — | — | — |
|
||||
| **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
|
||||
min/case on GPU, orders of magnitude above the intensity-metric tools.)*
|
||||
accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.6
|
||||
min/case on GPU, orders of magnitude above the intensity-metric tools.
|
||||
elastix, simpleitk and synthmorph are deterministic at fixed threads.)*
|
||||
|
||||
Per-case (mid / std, `OK` = accepted):
|
||||
|
||||
| case | simpleitk | antspyx | elastix | impact |
|
||||
|---|---|---|---|---|
|
||||
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
|
||||
| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 24.4 / 27.7 | 26.4 / 30.3 | **8.4 / 12.6** |
|
||||
| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.4 / 24.8 | 21.1 / 23.3 | 21.4 / 24.1 |
|
||||
| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.5 / 18.4 | 10.6 / 8.2 **OK** | 10.0 / 8.3 **OK** |
|
||||
| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 14.9 / 15.5 | 11.0 / 6.1 **OK** | 11.1 / 8.8 **OK** |
|
||||
| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.3 / 17.5 | 15.3 / 14.0 | 54.4 / 15.1 |
|
||||
| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 17.7 / 37.6 | 11.4 / 3.2 **OK** | 23.9 / 22.6 |
|
||||
| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** | 10.1 / 5.9 **OK** |
|
||||
| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.1 / 14.3 | 12.0 / 6.8 **OK** | 12.4 / 22.3 |
|
||||
| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 17.2 / 35.4 | 18.9 / 22.7 | 19.6 / 30.3 |
|
||||
| case | antspyx | elastix | impact | simpleitk | synthmorph |
|
||||
|---|---|---|---|---|---|
|
||||
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
|
||||
| m6_22IRRKU6_20240115 | 18.9 / 43.8 | 27.0 / 29.3 | 8.3 / 13.2 | 26.2 / 29.6 | 10.4 / 23.7 |
|
||||
| m6_23XIMO4M_20230222 | 22.9 / 24.5 | 22.6 / 23.9 | 23.1 / 25.7 | 22.6 / 24.0 | 20.3 / 23.3 |
|
||||
| m6_24V62QQM_20230223 | 18.7 / 47.5 | 11.3 / 16.9 | 9.9 / 7.9 **OK** | 10.2 / 16.2 | 10.4 / 22.9 |
|
||||
| m6_24V62QQM_20231017 | 15.3 / 37.9 | 16.6 / 11.1 **OK** | 6.6 / 6.9 **OK** | 11.9 / 7.7 **OK** | 11.5 / 8.4 **OK** |
|
||||
| m6_24V62QQM_20240820 | 38.7 / 54.0 | 10.2 / 8.2 **OK** | 37.1 / 45.6 | 10.8 / 9.0 **OK** | 43.4 / 15.8 |
|
||||
| m6_26TYQI3C_20221027 | 11.2 / 3.7 **OK** | 11.7 / 3.7 **OK** | 23.1 / 13.0 | 11.1 / 2.6 **OK** | 11.8 / 2.9 **OK** |
|
||||
| m6_2BBN5XA4_20250411 | 20.6 / 24.2 | 25.3 / 30.8 | 10.1 / 6.5 **OK** | 19.1 / 8.8 **OK** | 52.4 / 53.0 |
|
||||
| m6_2C4TKI6H_20240506 | 15.3 / 14.2 | 14.2 / 15.3 | 13.1 / 17.4 | 12.7 / 13.5 | 16.7 / 28.2 |
|
||||
| m6_2CDNN7BH_20240102 | 19.4 / 31.7 | 19.5 / 30.3 | 19.6 / 30.9 | 19.3 / 30.1 | 18.9 / 16.8 |
|
||||
|
||||
### Interpretation
|
||||
|
||||
- **Translation is recovered well by all four** (mid gaps mostly 6–20 mm,
|
||||
within the `≤ 20 mm` band). `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.
|
||||
- **Translation is recovered well by all five** (mid gaps mostly 6–20 mm,
|
||||
within the `≤ 20 mm` band, except the grossly-drifted cases below).
|
||||
`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), **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),
|
||||
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.
|
||||
- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
|
||||
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
|
||||
(RSGD wall time is thread-sensitive, convergence is not); the in-process
|
||||
|
|
@ -167,14 +199,15 @@ Per-case (mid / std, `OK` = accepted):
|
|||
## Reproducing
|
||||
|
||||
```bash
|
||||
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here); deepreg
|
||||
# reports unavailable and is recorded as such
|
||||
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here)
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py --n-cases 10 \
|
||||
--workers 4 --threads 24 --simpleitk-mode benchmark
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx \
|
||||
--n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
|
||||
|
||||
# IMPACT (needs Docker + the pulled image + a model; see the Impact section
|
||||
# below for the one-time setup)
|
||||
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
|
||||
# sections below for the one-time setup)
|
||||
python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \
|
||||
--workers 4 --threads 24
|
||||
python scripts/bench_registration.py --pipelines impact --n-cases 10 \
|
||||
--workers 4 --threads 24
|
||||
|
||||
|
|
@ -185,7 +218,8 @@ python scripts/bench_registration.py --merge
|
|||
Other flags: `--sample random --seed 0` (reproducible random 10),
|
||||
`--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`,
|
||||
`--threads T`, `--simpleitk-mode {benchmark,production}`,
|
||||
`--deepreg-ckpt <.ckpt>` (DeepReg), `--impact-model <rel|abs>`,
|
||||
`--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`,
|
||||
`--synthmorph-extent <int>`, `--impact-model <rel|abs>`,
|
||||
`--impact-gpu <int|-1>`, `--impact-iters N`, `--impact-samples N`,
|
||||
`--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
|
||||
|
||||
|
|
@ -216,35 +250,61 @@ visuals agree. Generate for the available pipelines with:
|
|||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py \
|
||||
--pipelines simpleitk elastix antspyx --n-cases 10 --qa-images
|
||||
python scripts/bench_registration.py \
|
||||
--pipelines synthmorph --n-cases 10 --qa-images
|
||||
python scripts/bench_registration.py \
|
||||
--pipelines impact --n-cases 10 --qa-images
|
||||
```
|
||||
|
||||
(9 images each; the incomplete-CT case is skipped.)
|
||||
|
||||
## Enabling the blocked pipeline
|
||||
## SynthMorph pipeline (Docker) — one-time setup
|
||||
|
||||
**DeepReg (Oxford)** — the PyPI `deepreg` (0.1.2) is an unrelated stub that pins
|
||||
`numpy==1.18.5`; do not install it. DeepReg is not a one-shot rigid tool: it is a
|
||||
CycleGAN image-to-image registration that must be **trained or run from a
|
||||
pretrained checkpoint**. To benchmark it:
|
||||
The `synthmorph` pipeline runs the [`freesurfer/synthmorph`](https://hub.docker.com/r/freesurfer/synthmorph)
|
||||
image: [SynthMorph](https://synthmorph.io) (Hoffmann et al., IEEE TMI 2022), an
|
||||
anatomy-aware, acquisition-agnostic **deep** registration (TensorFlow
|
||||
`VxmAffineFeatureDetector` with the pretrained `synthmorph.rigid.1.h5` weights).
|
||||
No skull-stripping or intensity preprocessing is needed: each volume is
|
||||
normalized to an isotropic 1-mm LIA network space internally. It is the
|
||||
second Docker-based pipeline and the only deep-learning one.
|
||||
|
||||
Prerequisites (check with `RP.availability()`; a missing piece raises
|
||||
`PipelineUnavailable` with the exact command):
|
||||
|
||||
```bash
|
||||
git clone https://github.com/OB-CVHU/DeepReg && cd DeepReg && pip install -e .
|
||||
python scripts/bench_registration.py --pipelines deepreg --n-cases 10 \
|
||||
--deepreg-ckpt /path/to/pretrained_3d.ckpt
|
||||
# 1. Docker (with NVIDIA container toolkit so --gpus works)
|
||||
docker --version
|
||||
|
||||
# 2. Pull the image (~11 GB, one-time)
|
||||
docker pull freesurfer/synthmorph
|
||||
```
|
||||
|
||||
The adapter runs `deepreg test` to predict a voxel displacement field, converts
|
||||
it to a SimpleITK `DisplacementFieldTransform`, and scores it like the others.
|
||||
No extra model assets: the rigid weights ship inside the image
|
||||
(`/freesurfer/models`).
|
||||
|
||||
### How the adapter drives the image
|
||||
|
||||
1. **Rigid matrix model, no init** — `register -m rigid` on GPU when one is
|
||||
visible (auto-spread over the host GPUs by pid; `--synthmorph-gpu -1` forces
|
||||
CPU). SynthMorph is a global deep method, so it recovers the gross
|
||||
inter-session offset on its own and needs **no** initial transform.
|
||||
2. **Read-back** — the `-T` LTA (a `LINEAR_RAS_TO_RAS` 4×4 in shared world
|
||||
space; with fixed=CT/moving=T1c the `-T` matrix is the T1c→CT contract
|
||||
transform) is parsed into a SimpleITK `AffineTransform` and run through the
|
||||
usual `fix_orientation` gate.
|
||||
3. **GPU sharing** — the container is run as the host uid (clean scratch
|
||||
removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
|
||||
concurrent benchmark workers can share a GPU.
|
||||
|
||||
## IMPACT pipeline (Docker) — one-time setup
|
||||
|
||||
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):
|
||||
|
|
@ -295,17 +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.
|
||||
- `deepreg` could not be executed here, so its row is "not installed", not
|
||||
measured; its code path is implemented and unit-verified for import/convention.
|
||||
- `fix_orientation` (elastix/ants/deepreg/impact) chooses between a transform
|
||||
- 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 (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. In
|
||||
practice it never fired (all `auto_inverted: False`): elastix/ants/impact
|
||||
returned the T1c→CT direction as expected.
|
||||
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 (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
|
||||
|
|
@ -313,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
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
"""Benchmark the 5 rigid CT->T1c registration pipelines on N M6 test cases.
|
||||
|
||||
Pipelines (see src/registration_pipelines.py for the shared transform contract):
|
||||
simpleitk | elastix | antspyx | deepreg | impact
|
||||
simpleitk | elastix | antspyx | synthmorph | impact
|
||||
|
||||
Each pipeline registers a head-and-neck CT to a tightly-cropped T1c MRI
|
||||
(cross-modality, inter-session tilt up to ~30 deg) and the GTV is warped onto
|
||||
|
|
@ -19,22 +19,22 @@ geometric QA used in production (scripts/02_build_m6_dataset.py):
|
|||
|
||||
Because some toolkits are not pip-installable into the main
|
||||
`longitudinal-313` env (antspyx is kept in an isolated env, see
|
||||
REGISTRATION_BENCHMARK.md; PyPI deepreg is a broken stub), the benchmark is
|
||||
env-portable: run it under whichever conda env has a given toolkit installed,
|
||||
and it records only the pipelines that are importable there. elastix runs in the main env through the ITK-Elastix
|
||||
python API (pip ``itk`` + ``itk-elastix``). ``impact`` runs in the
|
||||
``vboussot/elastix_impact`` Docker image (needs Docker + the image + a model,
|
||||
see REGISTRATION_BENCHMARK.md). The `--merge` mode then combines
|
||||
the per-pipeline JSON files (accumulated across envs) into one summary + plot.
|
||||
REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under
|
||||
whichever conda env has a given toolkit installed, and it records only the
|
||||
pipelines that are importable there. elastix runs in the main env through the
|
||||
ITK-Elastix python API (pip ``itk`` + ``itk-elastix``). ``synthmorph`` and
|
||||
``impact`` run in Docker images (``freesurfer/synthmorph`` and
|
||||
``vboussot/elastix_impact`` respectively; needs Docker + the pulled image, see
|
||||
REGISTRATION_BENCHMARK.md). The `--merge` mode then combines the
|
||||
per-pipeline JSON files (accumulated across envs) into one summary + plot.
|
||||
|
||||
Usage (repo root):
|
||||
# main env (SimpleITK + elastix + antspyx; impact needs Docker + the image;
|
||||
# deepreg unavailable)
|
||||
# main env (SimpleITK + elastix + antspyx)
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx
|
||||
|
||||
# impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
|
||||
python scripts/bench_registration.py --pipelines impact
|
||||
# synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
|
||||
python scripts/bench_registration.py --pipelines synthmorph impact
|
||||
|
||||
# any pipeline that is not importable in the current env is recorded as
|
||||
# unavailable, so the script is env-portable; --merge combines the
|
||||
|
|
@ -225,9 +225,10 @@ def _call(pipeline, ct, t1c, tv, args):
|
|||
return RP.register_elastix(ct, t1c, workdir=args.elastix_dir, threads=args.threads)
|
||||
if pipeline == "antspyx":
|
||||
return RP.register_antspyx(ct, t1c)
|
||||
if pipeline == "deepreg":
|
||||
return RP.register_deepreg(ct, t1c, workdir=args.deepreg_dir,
|
||||
checkpoint=args.deepreg_ckpt)
|
||||
if pipeline == "synthmorph":
|
||||
return RP.register_synthmorph(ct, t1c, model=args.synthmorph_model,
|
||||
gpu=args.synthmorph_gpu, threads=args.threads,
|
||||
extent=args.synthmorph_extent)
|
||||
if pipeline == "impact":
|
||||
return RP.register_impact(ct, t1c, model=args.impact_model, gpu=args.impact_gpu,
|
||||
iters=args.impact_iters, samples=args.impact_samples)
|
||||
|
|
@ -476,8 +477,15 @@ def parse_args():
|
|||
"production = full multi-start tilt grid (slow)")
|
||||
ap.add_argument("--out-dir", default=os.path.join(ROOT, "results", "bench_registration"))
|
||||
ap.add_argument("--elastix-dir", default=None, help="elastix scratch dir (default data/m6reg/elastix)")
|
||||
ap.add_argument("--deepreg-dir", default=None, help="deepreg scratch dir")
|
||||
ap.add_argument("--deepreg-ckpt", default=None, help="path to a DeepReg 3D pretrained checkpoint")
|
||||
ap.add_argument("--synthmorph-model", choices=["rigid", "affine"], default="rigid",
|
||||
help="SynthMorph matrix model (default rigid; deformable/joint are "
|
||||
"outside the rigid-transform contract)")
|
||||
ap.add_argument("--synthmorph-gpu", type=int, default=None,
|
||||
help="GPU device index for SynthMorph (default: auto-spread over "
|
||||
"the host GPUs by pid; -1 forces CPU)")
|
||||
ap.add_argument("--synthmorph-extent", type=int, default=256,
|
||||
help="SynthMorph isotropic registration-space extent in mm "
|
||||
"(192 or 256; default 256)")
|
||||
ap.add_argument("--impact-model", default=None,
|
||||
help="IMPACT TorchScript model, path relative to the HF "
|
||||
f"layout (default {RP.IMPACT_DEFAULT_MODEL}) or an absolute host path")
|
||||
|
|
@ -506,7 +514,6 @@ def main():
|
|||
merge(args)
|
||||
return
|
||||
args.elastix_dir = args.elastix_dir or d("data/m6reg/elastix")
|
||||
args.deepreg_dir = args.deepreg_dir or d("data/m6reg/deepreg")
|
||||
if args.qa_images:
|
||||
args.qa_out = os.path.abspath(args.qa_out or os.path.join(args.out_dir, "qa"))
|
||||
os.makedirs(args.qa_out, exist_ok=True)
|
||||
|
|
|
|||
|
|
@ -21,23 +21,39 @@ 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_deepreg`` : Oxford DeepReg (CycleGAN deformation-field
|
||||
registration; needs a pretrained checkpoint, see module docstring).
|
||||
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 DeepReg, missing its checkpoint, or,
|
||||
for IMPACT, missing its Docker image / model)
|
||||
A toolkit that is not installed (or, for the Docker-based pipelines, missing
|
||||
its image)
|
||||
raises :class:`PipelineUnavailable` with the exact install / precondition. The
|
||||
benchmark catches that and reports the pipeline as *not available* instead of
|
||||
failing the whole run, so the script is env-portable: run it under whichever
|
||||
|
|
@ -50,11 +66,10 @@ per toolkit. :func:`fix_orientation` resolves this by warping the CT with both
|
|||
``T`` and ``T^-1`` and keeping whichever gives the smaller skull-ring
|
||||
concentricity gap (a pure correctness gate -- a wrongly-oriented rigid transform
|
||||
lands the head on the far side of the box and yields a pathological gap). It is
|
||||
applied to elastix / antspyx / deepreg / impact, not to the trusted production
|
||||
SimpleITK path.
|
||||
applied to elastix / antspyx / synthmorph / impact, not to the trusted
|
||||
production SimpleITK path.
|
||||
"""
|
||||
import os
|
||||
import json
|
||||
import shutil
|
||||
import subprocess
|
||||
import importlib.util
|
||||
|
|
@ -169,10 +184,11 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
|
|||
"""SimpleITK rigid CT->T1c, Mattes MI on a CT-space head box.
|
||||
|
||||
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).
|
||||
|
|
@ -209,12 +225,11 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
|
|||
ref.SetOrigin(tuple(float(x) for x in lo))
|
||||
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()
|
||||
|
|
@ -296,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:
|
||||
|
|
@ -323,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",
|
||||
|
|
@ -361,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()
|
||||
|
|
@ -373,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:
|
||||
|
|
@ -385,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
|
||||
|
||||
|
|
@ -399,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=...,
|
||||
|
|
@ -407,8 +411,11 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
|
|||
multi-resolution Mattes-MI rigid registration -- the cross-modality metric
|
||||
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
|
||||
|
|
@ -430,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))
|
||||
|
|
@ -457,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__', '?')}",
|
||||
|
|
@ -484,107 +492,141 @@ def sys_pyver():
|
|||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# 4. DeepReg (Oxford)
|
||||
# 4. SynthMorph (Docker image freesurfer/synthmorph)
|
||||
# --------------------------------------------------------------------------- #
|
||||
|
||||
def register_deepreg(ct, t1c, workdir=None, checkpoint=None, config=None):
|
||||
"""Oxford DeepReg deformation-field registration (CT -> T1c).
|
||||
SYNTHMORPH_IMAGE = "freesurfer/synthmorph"
|
||||
# The image entrypoint is ``mri_synthmorph``; the ``register`` subcommand and
|
||||
# options follow. ``TF_FORCE_GPU_ALLOW_GROWTH`` lets concurrent benchmark
|
||||
# workers share a GPU instead of each pre-allocating all of it.
|
||||
|
||||
DeepReg is not a one-shot rigid tool: it is a CycleGAN-based image-to-image
|
||||
registration that must be *trained* or run from a *pretrained checkpoint*.
|
||||
It predicts a voxel displacement field; we convert that to a SimpleITK
|
||||
DisplacementFieldTransform to warp CT/GTV onto the T1c grid. Preconditions:
|
||||
- ``deepreg`` installed from git (https://github.com/OB-CVHU/DeepReg)
|
||||
-- the PyPI ``deepreg`` (0.1.2) is an unrelated stub that pins
|
||||
numpy==1.18.5 and will not work.
|
||||
- a 3D pretrained checkpoint (``checkpoint``) -- typically one trained on
|
||||
a CT/MRI template pair; see the repo's pretrained-models note.
|
||||
"""
|
||||
|
||||
def _synthmorph_docker_ok():
|
||||
"""Return (ok, note): is the Docker CLI + the synthmorph image available?"""
|
||||
if shutil.which("docker") is None:
|
||||
return False, "docker CLI not found on PATH"
|
||||
try:
|
||||
import deepreg # noqa
|
||||
r = subprocess.run(["docker", "image", "inspect", SYNTHMORPH_IMAGE],
|
||||
capture_output=True, text=True, timeout=30)
|
||||
if r.returncode != 0:
|
||||
return False, (f"docker image {SYNTHMORPH_IMAGE} not pulled "
|
||||
f"(run: docker pull {SYNTHMORPH_IMAGE})")
|
||||
except Exception as e: # noqa
|
||||
raise PipelineUnavailable(
|
||||
f"deepreg not importable ({e!r}). Install the Oxford DeepReg from "
|
||||
"source: `git clone https://github.com/OB-CVHU/DeepReg && "
|
||||
"cd DeepReg && pip install -e .` (NOT the PyPI `deepreg` stub). "
|
||||
"Also requires a pretrained 3D checkpoint (see --deepreg-ckpt).")
|
||||
if not checkpoint or not os.path.exists(checkpoint):
|
||||
raise PipelineUnavailable(
|
||||
"deepreg: no pretrained checkpoint provided. Pass --deepreg-ckpt "
|
||||
"to a 3D .ckpt produced by `deepreg train` (or a released "
|
||||
"pretrained model). DeepReg cannot do one-shot registration "
|
||||
"without weights.")
|
||||
workdir = workdir or d("data/m6reg/deepreg")
|
||||
workdir = os.path.join(workdir, f"pid{os.getpid()}") # per-worker scratch
|
||||
os.makedirs(workdir, exist_ok=True)
|
||||
cfg = _deepreg_config(workdir, ct, t1c, checkpoint)
|
||||
outdir = os.path.join(workdir, "pred")
|
||||
subprocess.run(
|
||||
["python", "-m", "deepreg.main", "test",
|
||||
"-checkpoint_path", checkpoint, "-config_path", cfg, "-out_dir", outdir],
|
||||
check=True, capture_output=True, text=True, cwd=workdir)
|
||||
# predicted deformation field (voxels) on the fixed (T1c) grid
|
||||
pred = _find_single_nii(os.path.join(outdir))
|
||||
if pred is None:
|
||||
raise PipelineUnavailable("deepreg: no predicted deformation field found")
|
||||
disp = sitk.GetArrayFromImage(sitk.ReadImage(pred)) # (Z, Y, X, 3)
|
||||
T = _displacement_to_sitk(disp, t1c)
|
||||
meta = {
|
||||
"pipeline": "deepreg", "toolkit": "Oxford DeepReg (CycleGAN)",
|
||||
"metric": "CycleGAN + registration loss (pretrained)",
|
||||
"transform": "3D displacement field", "convention": "T1c->CT (verified)",
|
||||
"auto_inverted": False, "checkpoint": checkpoint,
|
||||
}
|
||||
return T, meta
|
||||
return False, f"docker image check failed ({type(e).__name__}: {e})"
|
||||
return True, SYNTHMORPH_IMAGE
|
||||
|
||||
|
||||
def _deepreg_config(workdir, ct, t1c, checkpoint):
|
||||
"""Minimal DeepReg test config: moving=CT, fixed=T1c, 3D, single channel."""
|
||||
from deepreg.registry import registry # noqa
|
||||
cfg = {
|
||||
"general": {"data_format": "nifti", "device": "cuda:0"},
|
||||
"data": {"moving": os.path.join(workdir, "moving"),
|
||||
"fixed": os.path.join(workdir, "fixed"),
|
||||
"field": "deformation", "moving_fields": "t1c", "fixed_fields": "t1c",
|
||||
"mask": None, "resample_size": "30 186 30"},
|
||||
"validation": {},
|
||||
"training": {"epoch": 1, "batch_size": 1,
|
||||
"moving_fields": "t1c", "fixed_fields": "t1c",
|
||||
"resample_size": "30 186 30", "resample_spacing": "1 1 1"},
|
||||
"testing": {"moving": os.path.join(workdir, "moving"),
|
||||
"fixed": os.path.join(workdir, "fixed"),
|
||||
"resample_size": "30 186 30", "resample_spacing": "1 1 1"},
|
||||
}
|
||||
p = os.path.join(workdir, "deepreg_test.yaml")
|
||||
def _lta_to_sitk(path):
|
||||
"""Parse a FreeSurfer LTA (LINEAR_RAS_TO_RAS) into a SimpleITK AffineTransform.
|
||||
|
||||
The 4x4 matrix block lives in the shared RAS world space. Both inputs are
|
||||
RAS-oriented by ``read_nii_img`` (DICOMOrient), so RAS == SimpleITK's
|
||||
physical frame and the matrix applies directly to world coordinates."""
|
||||
lines = [l for l in open(path).read().splitlines()
|
||||
if l and not l.startswith("#")]
|
||||
i = lines.index("1 4 4")
|
||||
M = np.array([[float(x) for x in lines[i + 1 + r].split()] for r in range(4)])
|
||||
T = sitk.AffineTransform(3)
|
||||
T.SetMatrix(M[:3, :3].ravel())
|
||||
T.SetTranslation(M[:3, 3])
|
||||
return T
|
||||
|
||||
|
||||
def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
|
||||
extent=256, keep_scratch=False):
|
||||
"""SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image
|
||||
``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``,
|
||||
pretrained ``synthmorph.rigid.1.h5`` weights).
|
||||
|
||||
SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware,
|
||||
acquisition-agnostic deep registration: each volume is normalized to an
|
||||
isotropic 1-mm LIA network space and the (pretrained) network estimates the
|
||||
transform, so no skull-stripping or intensity preprocessing is needed. Being
|
||||
a global deep method it recovers the gross inter-session offset on its own,
|
||||
so **no centroid init is used**. The pipeline
|
||||
- runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible
|
||||
(auto-spread over the host GPUs by pid; ``gpu=-1`` forces CPU),
|
||||
- reads back the ``-T`` LTA (fixed->moving = the T1c->CT transform the
|
||||
unified contract requires); :func:`fix_orientation` verifies the
|
||||
direction against the skull ring.
|
||||
|
||||
Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI and
|
||||
the ``freesurfer/synthmorph`` image. Only the matrix models (``rigid``,
|
||||
``affine``) are exposed: the deformable/joint models emit displacement
|
||||
fields, outside the rigid-transform contract.
|
||||
"""
|
||||
import time
|
||||
|
||||
ok, note = _synthmorph_docker_ok()
|
||||
if not ok:
|
||||
raise PipelineUnavailable(f"synthmorph: {note}")
|
||||
if model not in ("rigid", "affine"):
|
||||
raise PipelineUnavailable(
|
||||
f"synthmorph: only the matrix models (rigid/affine) fit the unified "
|
||||
f"rigid-transform contract, got {model!r}")
|
||||
|
||||
# GPU: auto-spread across the host GPUs by pid, or forced value
|
||||
if gpu is None:
|
||||
ngpu = _n_gpus()
|
||||
gpu = (os.getpid() % ngpu) if ngpu else -1
|
||||
|
||||
# per-case scratch on a local (non-NFS) tmp fs
|
||||
tmp = tempfile.mkdtemp(prefix="synthmorph_reg_")
|
||||
try:
|
||||
import yaml
|
||||
with open(p, "w") as f:
|
||||
yaml.safe_dump(cfg, f)
|
||||
except Exception:
|
||||
with open(p, "w") as f:
|
||||
json.dump(cfg, f)
|
||||
return p
|
||||
# fixed = CT, moving = T1c (shared convention for all pipelines)
|
||||
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "fixed.nii.gz"), True) # CT
|
||||
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
|
||||
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()}",
|
||||
"-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"]
|
||||
if gpu >= 0:
|
||||
cmd += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"]
|
||||
cmd += ["-v", f"{tmp}:/data", SYNTHMORPH_IMAGE,
|
||||
"register", "-m", model]
|
||||
if gpu >= 0:
|
||||
cmd.append("-g")
|
||||
if threads:
|
||||
cmd += ["-j", str(threads)]
|
||||
cmd += ["-e", str(extent),
|
||||
"-T", "/data/f2m.lta",
|
||||
"/data/moving.nii.gz", "/data/fixed.nii.gz"]
|
||||
t0 = time.time()
|
||||
proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200)
|
||||
dt = time.time() - t0
|
||||
if proc.returncode != 0:
|
||||
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
|
||||
raise RuntimeError(f"synthmorph: docker failed (rc={proc.returncode}):\n{tail}")
|
||||
if not os.path.exists(os.path.join(tmp, "f2m.lta")):
|
||||
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
|
||||
raise RuntimeError(f"synthmorph: no /data/f2m.lta produced:\n{tail}")
|
||||
|
||||
def _find_single_nii(dirpath):
|
||||
if not os.path.isdir(dirpath):
|
||||
return None
|
||||
ns = [os.path.join(dirpath, x) for x in sorted(os.listdir(dirpath))
|
||||
if x.endswith(".nii.gz") or x.endswith(".nii")]
|
||||
return ns[0] if ns else None
|
||||
T = _lta_to_sitk(os.path.join(tmp, "f2m.lta"))
|
||||
T, inverted = fix_orientation(ct, t1c, T)
|
||||
finally:
|
||||
if not keep_scratch:
|
||||
shutil.rmtree(tmp, ignore_errors=True)
|
||||
|
||||
|
||||
def _displacement_to_sitk(disp, ref):
|
||||
"""Build a sitk.DisplacementFieldTransform from a (Z,Y,X,3) voxel field."""
|
||||
from scipy.ndimage import map_coordinates
|
||||
zz, yy, xx = np.mgrid[0:disp.shape[0], 0:disp.shape[1], 0:disp.shape[2]]
|
||||
# disp[i,j,k] = (dz, dy, dx) to add at grid point (i,j,k) to reach the source
|
||||
src_z = zz + disp[..., 0]; src_y = yy + disp[..., 1]; src_x = xx + disp[..., 2]
|
||||
# DisplacementFieldTransform stores, per output voxel, the offset to the input
|
||||
off = sitk.GetImageFromArray(disp[..., ::-1].astype(np.float32)) # (X,Y,Z,3)
|
||||
off.CopyInformation(ref)
|
||||
t = sitk.DisplacementFieldTransform(off)
|
||||
return t
|
||||
meta = {
|
||||
"pipeline": "synthmorph",
|
||||
"toolkit": f"Docker {SYNTHMORPH_IMAGE} (TensorFlow, pretrained)",
|
||||
"metric": "SynthMorph anatomy-aware contrast-invariant (deep, pretrained)",
|
||||
"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": "CT", "moving": "T1c",
|
||||
"seconds_reg": round(dt, 1),
|
||||
"docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu,
|
||||
"extent": extent,
|
||||
}
|
||||
if keep_scratch:
|
||||
meta["scratch"] = tmp
|
||||
return T, meta
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
|
|
@ -618,7 +660,7 @@ def _impact_docker_ok():
|
|||
_NGPUS_CACHE = None
|
||||
|
||||
|
||||
def _impact_n_gpus():
|
||||
def _n_gpus():
|
||||
global _NGPUS_CACHE
|
||||
if _NGPUS_CACHE is None:
|
||||
n = 0
|
||||
|
|
@ -752,6 +794,13 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
"""IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image
|
||||
``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
|
||||
|
|
@ -802,14 +851,15 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference
|
||||
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]))
|
||||
|
||||
# GPU: auto-spread across the host GPUs by pid, or forced value
|
||||
if gpu is None:
|
||||
ngpu = _impact_n_gpus()
|
||||
ngpu = _n_gpus()
|
||||
gpu = (os.getpid() % ngpu) if ngpu else -1
|
||||
|
||||
# per-level feature-extraction scales (recommended 6/3/1.5/1 mm), scaled so
|
||||
|
|
@ -885,7 +935,8 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
|
||||
|
||||
# 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)
|
||||
|
|
@ -913,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,
|
||||
"deepreg": register_deepreg,
|
||||
"elastix": register_elastix,
|
||||
"impact": register_impact,
|
||||
"simpleitk": register_simpleitk,
|
||||
"synthmorph": register_synthmorph,
|
||||
}
|
||||
ALL_PIPELINES = list(REGISTERS)
|
||||
|
||||
|
|
@ -943,12 +994,12 @@ def availability():
|
|||
out[name] = (True, f"ants importable v{getattr(ants, '__version__', '?')}")
|
||||
except Exception as e: # noqa
|
||||
out[name] = (False, f"ants(antspyx) not importable ({type(e).__name__})")
|
||||
elif name == "deepreg":
|
||||
try:
|
||||
import deepreg # noqa
|
||||
out[name] = (True, "deepreg importable")
|
||||
except Exception as e: # noqa
|
||||
out[name] = (False, "deepreg not importable")
|
||||
elif name == "synthmorph":
|
||||
ok, note = _synthmorph_docker_ok()
|
||||
if not ok:
|
||||
out[name] = (False, note)
|
||||
else:
|
||||
out[name] = (True, f"Docker {SYNTHMORPH_IMAGE}")
|
||||
elif name == "impact":
|
||||
ok, note = _impact_docker_ok()
|
||||
if not ok:
|
||||
|
|
|
|||
Loading…
Reference in a new issue