Compare commits
2 commits
4ad972b00f
...
a65189cd33
| Author | SHA1 | Date | |
|---|---|---|---|
| a65189cd33 | |||
| 8c083abe92 |
4 changed files with 564 additions and 452 deletions
|
|
@ -10,7 +10,7 @@ Activate with:
|
|||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
```
|
||||
|
||||
Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK.
|
||||
Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), plus Docker images `freesurfer/synthmorph` and `vboussot/elastix_impact`.
|
||||
|
||||
## Project
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
# Registration Pipeline Benchmark — CT → T1c (M6)
|
||||
|
||||
Benchmarks **five rigid, cross-modality registration pipelines** that align a
|
||||
Benchmarks **six rigid, cross-modality registration pipelines** that align a
|
||||
head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be warped onto the
|
||||
T1c (MR) native grid. This is the registration step of
|
||||
`scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT
|
||||
|
|
@ -22,34 +22,44 @@ sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0) # GTV onto the MR grid
|
|||
```
|
||||
|
||||
places the CT / GTV concentric with the MR head. In `Resample`'s
|
||||
destination→source convention this `T` maps **T1c → CT**. All pipelines use
|
||||
**fixed = CT, moving = T1c** except `impact`, which keeps **fixed = T1c,
|
||||
moving = CT** (swapping its roles breaks the TotalSegmentator-based metric,
|
||||
which aborts with "Error in metric" or converges to a wrong rotation). The
|
||||
intensity-metric pipelines register on the same CT-space head box (MI computed
|
||||
on head tissue, not neck/air); `synthmorph` uses the full volumes (it
|
||||
normalizes each one to an isotropic network space internally). The local
|
||||
optimizers that cannot recover the gross inter-session offset from identity
|
||||
(`antspyx`, `impact`) start from a centroid-aligned translation
|
||||
(`centroid_offset`); the others start from identity, so the comparison is on
|
||||
each toolkit's native optimizer / metric / global-search capability.
|
||||
destination→source convention this `T` maps **T1c → CT**. All pipelines
|
||||
register on the **full CT / T1c volumes** (no head box) and use each toolkit's
|
||||
**native automatic transform initialization** (no manual centroid offset):
|
||||
elastix / impact use `AutomaticTransformInitialization` (CenterOfGravity),
|
||||
antspyx its default center-of-mass alignment, simpleitk starts from identity,
|
||||
niftyreg (reg_aladin) initializes from the image centres, synthmorph needs
|
||||
none. All use **fixed = CT, moving = T1c** except `impact`, which keeps
|
||||
**fixed = T1c, moving = CT**: the CT FOV (~300×300×251 mm) is much larger than
|
||||
the T1c crop (~220×220×150 mm), so with CT as the fixed image too many CT
|
||||
samples map outside the T1c buffer and the IMPACT metric aborts ("too many
|
||||
samples map outside moving image buffer"). The comparison is on each
|
||||
toolkit's native optimizer / metric / global-search capability.
|
||||
|
||||
## The five pipelines
|
||||
## The six pipelines
|
||||
|
||||
Listed alphabetically. All use **fixed = CT, moving = T1c** except `impact`,
|
||||
which keeps **fixed = T1c, moving = CT** (swapping breaks its TotalSegmentator
|
||||
metric). A pipeline that is a local optimizer and cannot recover the gross
|
||||
inter-session offset from identity (`antspyx`, `impact`) starts from a
|
||||
centroid-aligned translation; the others start from identity / their own global
|
||||
search.
|
||||
Listed alphabetically. All register on the **full volumes** (no head box) and
|
||||
use each toolkit's **native automatic transform initialization** (no manual
|
||||
centroid offset): elastix / impact use elastix's `AutomaticTransformInitialization`
|
||||
(CenterOfGravity), antspyx uses its default center-of-mass alignment, simpleitk
|
||||
starts from identity, niftyreg initializes from the image centres, and
|
||||
synthmorph needs none (global deep method). All use **fixed = CT, moving =
|
||||
T1c** except `impact`, which keeps **fixed = T1c, moving = CT** (swapping
|
||||
breaks its TotalSegmentator metric — see below).
|
||||
|
||||
| # | pipeline | fixed / moving | centroid init | metric | transform | status in this env |
|
||||
|---|----------|----------------|---------------|--------|-----------|--------------------|
|
||||
| 1 | `antspyx` | CT / T1c | **yes** (required) | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
|
||||
| 2 | `elastix` | CT / T1c | no (identity) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
|
||||
| 3 | `impact` | **T1c / CT** | **yes** | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) |
|
||||
| 4 | `simpleitk` | CT / T1c | no (identity) | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
|
||||
| 5 | `synthmorph` | CT / T1c | no (global deep) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
|
||||
| # | pipeline | fixed / moving | init | metric | transform | status in this env |
|
||||
|---|----------|----------------|------|--------|-----------|--------------------|
|
||||
| 1 | `antspyx` | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) |
|
||||
| 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (32 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
|
||||
| 3 | `impact` | **T1c / CT** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (main env, Docker + GPU) |
|
||||
| 4 | `niftyreg` | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (main env, conda-forge `reg_aladin` 1.5.69) |
|
||||
| 5 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
|
||||
| 6 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) |
|
||||
|
||||
Benchmark budget (best-results settings): `elastix`, `impact` and `simpleitk`
|
||||
run with **MaximumNumberOfIterations = 2000** (per resolution / per pyramid
|
||||
level), with the sampler and metric left at their defaults
|
||||
(**NumberOfSpatialSamples 5000**, **NumberOfHistogramBins 32**). `antspyx`,
|
||||
`niftyreg` and `synthmorph` keep their native defaults.
|
||||
|
||||
`simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid,
|
||||
centered start — the default, comparable to the other tools) and
|
||||
|
|
@ -73,135 +83,123 @@ m6_2BBN5XA4_20250411 m6_2C4TKI6H_20240506 m6_2CDNN7BH_20240102
|
|||
|
||||
## Metrics (identical for every pipeline)
|
||||
|
||||
Production QA, crop-independent:
|
||||
The headline alignment score is **`mi_ct_mr`** — the registration's *achieved*
|
||||
Mattes MI, read off as a standalone score:
|
||||
|
||||
- **`skull_ring_gap_mm`** — march rays from the MR head center; the gap between
|
||||
the MR head silhouette boundary and the warped-CT skull ring (p97 of positive
|
||||
CT) at mid and low levels, plus the angular std. **Aligned ⇔
|
||||
`0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`.** `mid` measures translational
|
||||
concentricity; `std` measures rotational concentricity.
|
||||
- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT
|
||||
head tissue (the slab must straddle the lesion). Acceptance needs `≥ 0.5`.
|
||||
- **`tv_nvox`** — warped GTV volume.
|
||||
- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the
|
||||
production acceptance rule).
|
||||
- **`mi_ct_mr`** — mean Mattes MI (nats) between the **CT inverse-transformed
|
||||
into MR coordinates** (the warped CT on the T1c grid) and the **original MR**,
|
||||
restricted to the MR head and each image percentile-binned independently
|
||||
(SimpleITK Mattes-MI style). This is the registration's *achieved* metric read
|
||||
off as a standalone alignment score: higher = the warped-CT and MR intensities
|
||||
are better in-sync = better aligned; it collapses toward 0 when the CT is
|
||||
rotated/translated away from the MR.
|
||||
(SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are
|
||||
better in-sync = better aligned; it collapses toward 0 when the CT is
|
||||
rotated/translated away from the MR. Tables and plots report this (plus
|
||||
registration time and GTV-in-CT coverage).
|
||||
- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT
|
||||
head tissue (the slab must straddle the lesion).
|
||||
- **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring).
|
||||
|
||||
The production geometric QA is still computed per case and kept in the JSON
|
||||
(it drives dataset acceptance in `02_build_m6_dataset.py`, but its summary
|
||||
numbers were found to correlate poorly with perceived alignment, so they are no
|
||||
longer reported in the tables/plots):
|
||||
|
||||
- **`skull_ring_gap_mm`** — march rays from the MR head center; the gap between
|
||||
the MR head silhouette boundary and the warped-CT skull ring (p97 of positive
|
||||
CT) at mid and low levels, plus the angular std. `mid` measures
|
||||
translational concentricity; `std` measures rotational concentricity.
|
||||
- **`tv_nvox`** — warped GTV volume.
|
||||
- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the
|
||||
production acceptance rule; aligned ⇔ `0 ≤ mid ≤ 20 mm` and `std ≤ 12 mm`).
|
||||
|
||||
## Results (10 cases)
|
||||
|
||||
| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|
||||
|---|---|---|---|---|---|---|---|---|
|
||||
| **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 18.9 / 20.1 | 31.7 / 31.3 | 10.9 | 0.94 | **0.436 / 0.393** |
|
||||
| **elastix** | ITK-Elastix 0.25.4 | 9 | **3/10 (0.30)** | 16.6 / 17.6 | 16.9 / 18.8 | 2.4 | 0.93 | **0.641 / 0.583** |
|
||||
| **impact** | Docker Elastix 5.2 + IMPACT | 9 | **3/10 (0.30)** | 13.1 / 16.8 | 13.2 / 18.6 | **576.0** | 0.93 | **0.469 / 0.482** |
|
||||
| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.7 / 16.0 | 13.5 / 15.7 | 7.6 | 0.93 | **0.645 / 0.632** |
|
||||
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | **2/10 (0.20)** | 16.7 / 21.8 | 22.9 / 21.7 | 26.0 | 0.93 | **0.241 / 0.308** |
|
||||
Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
|
||||
the MR head; higher = better).
|
||||
|
||||
*(antspyx and impact are stochastic (random metric/patch subsampling), so their
|
||||
accepted counts are one draw each — see Interpretation. impact's s/reg is ~9.6
|
||||
min/case on GPU, orders of magnitude above the intensity-metric tools.
|
||||
elastix, simpleitk and synthmorph are deterministic at fixed threads.)*
|
||||
|
||||
Per-case (mid / std, `OK` = accepted):
|
||||
|
||||
| case | antspyx | elastix | impact | simpleitk | synthmorph |
|
||||
| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|
||||
|---|---|---|---|---|---|
|
||||
| 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 |
|
||||
| **antspyx** | ANTsPy 0.6.3 | 9 | 35.0 | 0.93 | **0.652 / 0.644** |
|
||||
| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.516 / 0.520** |
|
||||
| **impact**\* | Docker Elastix 5.2 + IMPACT | 9 | 833.5 | 0.93 | **0.511 / 0.493** |
|
||||
| **niftyreg** | NiftyReg 1.5.69 (reg_aladin) | 9 | 240.1 | 0.93 | **0.296 / 0.339** |
|
||||
| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.4 | 0.93 | **0.651 / 0.641** |
|
||||
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 25.1 | 0.93 | **0.238 / 0.308** |
|
||||
|
||||
\* `impact` at its **previous** 500-iter / 2000-sample budget; the new default
|
||||
(MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000) is ~10× slower
|
||||
(hours/case) and was not re-run.
|
||||
|
||||
Per-case MI (nats; **bold** = best for the case):
|
||||
|
||||
| case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
|
||||
|---|---|---|---|---|---|---|
|
||||
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
|
||||
| m6_22IRRKU6_20240115 | 0.649 | 0.472 | **0.624** | 0.296 | 0.642 | 0.235 |
|
||||
| m6_23XIMO4M_20230222 | 0.655 | 0.416 | 0.474 | **0.634** | 0.652 | 0.510 |
|
||||
| m6_24V62QQM_20230223 | 0.703 | **0.701** | 0.554 | 0.260 | 0.663 | 0.238 |
|
||||
| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.419 | 0.493 | **0.602** | 0.516 |
|
||||
| m6_24V62QQM_20240820 | 0.652 | 0.643 | 0.458 | 0.118 | **0.652** | 0.096 |
|
||||
| m6_26TYQI3C_20221027 | 0.657 | 0.653 | 0.511 | 0.608 | **0.695** | 0.604 |
|
||||
| m6_2BBN5XA4_20250411 | 0.576 | 0.287 | 0.155 | 0.069 | **0.577** | 0.060 |
|
||||
| m6_2C4TKI6H_20240506 | 0.645 | 0.401 | **0.626** | 0.231 | 0.634 | 0.209 |
|
||||
| m6_2CDNN7BH_20240102 | 0.655 | 0.516 | **0.618** | 0.344 | 0.651 | 0.301 |
|
||||
|
||||
### Interpretation
|
||||
|
||||
- **Translation is recovered well by all five** (mid gaps mostly 6–20 mm,
|
||||
within the `≤ 20 mm` band, except the grossly-drifted cases below).
|
||||
`tv-in-CT` is identical (~0.93): the warped GTV lands inside the CT head for
|
||||
every pipeline, so the target is in the right place.
|
||||
- **The binding failure is rotation** — `std > 12 mm` on the repositioned cases
|
||||
(the ~30° inter-session tilt produces local minima that a single start
|
||||
(identity or centroid) can lock onto). This is exactly why production
|
||||
`02_build_m6_dataset.py` uses a **multi-start tilt grid**.
|
||||
- **SimpleITK is the most reliable intensity-metric tool in this draw**
|
||||
(4/10 accepted, best median mid 12.7), and **Elastix is the fastest**
|
||||
(2.4 s/case median wall vs 7.6 s SimpleITK / 10.9 s ANTs) with the best
|
||||
median MI (0.641) but 3/10 accepted. Under the fixed=CT/moving=T1c
|
||||
convention the two settle into different local minima on the hard tilted
|
||||
cases: Elastix accepts `m6_24V62QQM_20240820` (std 8.2) and `m6_24V62QQM_20231017`
|
||||
(std 11.1) but misses `m6_24V62QQM_20230223` (std 16.9) that SimpleITK handles
|
||||
better, and vice-versa on the drift cases. The fixed/moving direction
|
||||
materially affects where the ASGD/RSGD optimizers converge on the tilted cases.
|
||||
- **ANTs' Rigid is the weakest of the intensity-metric tools** (1/10 in this
|
||||
draw; it is stochastic, so the count varies): it frequently locks onto the
|
||||
tilted local minimum (median std 31.7, worst of all). It accepts only
|
||||
`m6_26TYQI3C` (mid 11.2 / std 3.7) and misses `m6_24V62QQM_20230223`
|
||||
(std 47.5) and `m6_24V62QQM_20231017` (std 37.9) that the other tools accept —
|
||||
textbook tilted local minima.
|
||||
- **IMPACT (semantic metric) is the standout on the hardest repositioned case**
|
||||
and the least reliable overall. On `m6_22IRRKU6` (the ~30° tilt every
|
||||
intensity tool fails) it gets **mid 8.3 / std 13.2** — the best translation of
|
||||
any pipeline and within ~1 mm of the acceptance threshold — because matching
|
||||
*anatomical semantics* (TotalSegmentator features) is far more discriminating
|
||||
for gross pose than raw cross-modality intensity. But the semantic landscape
|
||||
is rugged: single-init ASGD lands in bad local minima on several cases, and
|
||||
two drift grossly (`m6_24V62QQM_20240820` mid 37, `m6_26TYQI3C` std 13.0).
|
||||
Net 3/10 accepted — between ANTs and SimpleITK — at ~9.6 min/case, two orders
|
||||
of magnitude slower than the intensity tools.
|
||||
- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~2.4 s
|
||||
(warm workers), SimpleITK ~7.6 s, ANTs ~10.9 s (CPU-bound),
|
||||
**SynthMorph ~26 s** (GPU: a single pretrained-model forward pass per case),
|
||||
**IMPACT ~576 s** (GPU-bound: every ASGD iteration runs the TotalSegmentator
|
||||
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.
|
||||
- **SimpleITK and ANTsPy tie for best** (median MI 0.651 / 0.652): on the
|
||||
well-positioned cases they cluster at 0.60–0.70 and both survive the harder
|
||||
repositioned cases (0.58–0.65 on `m6_2BBN5XA4` / `m6_22IRRKU6`). SimpleITK's
|
||||
per-case MI is **unchanged** from the previous 150-iteration budget — RSGD
|
||||
converges well before the 2000-iteration cap, so the cap is now headroom
|
||||
rather than a constraint.
|
||||
- **The grossly repositioned cases define the ranking.** On
|
||||
`m6_2BBN5XA4` and `m6_24V62QQM_20240820` (~30° inter-session tilt) the tools
|
||||
split: ANTs/SimpleITK hold 0.58–0.65, elastix drops to 0.29–0.64, niftyreg
|
||||
0.07–0.12, synthmorph 0.06–0.10 — single-start optimizers locking onto local
|
||||
minima. This is exactly why production `02_build_m6_dataset.py` uses a
|
||||
**multi-start tilt grid**.
|
||||
- **Elastix is second at the new budget** (2000 iters per resolution, 32 bins,
|
||||
5000 samples): median 0.516, below the 0.641 recorded at the old 500/500/250
|
||||
+ 10-bin setting. The long ASGD run nails the clean cases (0.64–0.70) but its
|
||||
single CenterOfGravity start still locks onto local minima on the
|
||||
repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M`
|
||||
0.416).
|
||||
- **IMPACT (semantic metric) is the strongest on the hardest cases** — top MI
|
||||
on `m6_22IRRKU6` (0.624), `m6_2C4TKI6H` (0.626) and `m6_2CDNN7BH` (0.618),
|
||||
where anatomical-segmentation similarity out-resolves intensity MI — but its
|
||||
median (0.511) is dragged by stochastic draws (`m6_2BBN5XA4` 0.155). It is
|
||||
also by far the slowest (~14 min/case at the shown budget; the new 2000/5000
|
||||
budget would be on the order of hours/case).
|
||||
- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching
|
||||
is strong on moderate cases (`m6_23XIMO4M` 0.634, best in class;
|
||||
`m6_26TYQI3C` 0.608) but falls into the same tilt local minima
|
||||
(`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). CPU-only and ~40× slower
|
||||
than elastix at full volumes.
|
||||
- **SynthMorph is weakest of the six on these cross-modality CT pairs**
|
||||
(median 0.238): its `rigid` model is trained on brain *MRI*, so the ~30°
|
||||
tilt defeats it on the repositioned cases (0.06–0.10) even though it matches
|
||||
the intensity tools on the cleanest (`m6_26TYQI3C` 0.604). It needs **no
|
||||
centroid init** (global deep method) and is deterministic.
|
||||
- **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the
|
||||
CT head everywhere, so coverage is not the discriminator — MI is.
|
||||
- **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s,
|
||||
SimpleITK 9.4 s, SynthMorph 25.1 s (GPU forward pass), ANTs 35.0 s,
|
||||
**NiftyReg 240.1 s** (CPU block matching), **IMPACT 833.5 s** (previous
|
||||
budget; GPU feature extraction with backprop per iteration).
|
||||
- **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
|
||||
Elastix ASGD is deterministic given its fixed seed. **IMPACT is the most
|
||||
stochastic**: its Jacobian mode draws random patches every iteration, so the
|
||||
final pose — and even *which* cases get accepted — changes run to run (two
|
||||
independent 10-case draws both accepted 3/10 but on different cases). Treat
|
||||
its per-case numbers as one sample from a distribution.
|
||||
SimpleITK's MI (and registration) is deterministic at fixed thread count;
|
||||
the in-process Elastix ASGD is deterministic given its fixed seed;
|
||||
SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**:
|
||||
its Jacobian mode draws random patches every iteration, so the final pose
|
||||
changes run to run. Treat its per-case numbers as one sample from a
|
||||
distribution.
|
||||
|
||||
## Reproducing
|
||||
|
||||
```bash
|
||||
# main env — SimpleITK + Elastix + ANTsPy (all pip-installed here)
|
||||
# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (pip/conda-installed here)
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx \
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg \
|
||||
--n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
|
||||
|
||||
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
|
||||
|
|
@ -217,10 +215,12 @@ 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}`,
|
||||
`--threads T` (also the OpenMP thread count for `reg_aladin`),
|
||||
`--simpleitk-mode {benchmark,production}`,
|
||||
`--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`,
|
||||
`--impact-gpu <int|-1>`, `--impact-iters N` (default 2000),
|
||||
`--impact-samples N` (default 5000),
|
||||
`--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
|
||||
|
||||
## Visual QA images (checkerboard + color overlay)
|
||||
|
|
@ -242,9 +242,9 @@ parenchyma is visible and comparable to the MR (the skull, ~800+ HU, clips to
|
|||
white); the MR uses a 2–98 % head window. The coronal and sagittal views are
|
||||
oriented superior-up to match the axial view.
|
||||
|
||||
The title carries the pipeline, case key, Mattes MI, and the skull-ring gap
|
||||
(mid/std), so a glance tells you both the achieved-metric score and whether the
|
||||
visuals agree. Generate for the available pipelines with:
|
||||
The title carries the pipeline, case key and the achieved Mattes MI, so a
|
||||
glance tells you whether the visuals agree with the score. Generate for the
|
||||
available pipelines with:
|
||||
|
||||
```bash
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
|
|
@ -296,6 +296,41 @@ No extra model assets: the rigid weights ship inside the image
|
|||
removal, TF still reaches the GPU) with `TF_FORCE_GPU_ALLOW_GROWTH` set so
|
||||
concurrent benchmark workers can share a GPU.
|
||||
|
||||
## NiftyReg pipeline — one-time setup
|
||||
|
||||
The `niftyreg` pipeline runs `reg_aladin` from the [NiftyReg](https://github.com/KCL-BMEIS/niftyreg)
|
||||
suite (KCL-BMEIS, the actively-maintained NiftyReg — the release zips for
|
||||
v2.x need a newer glibc, so this env uses the conda-forge build):
|
||||
**symmetric block matching + trimmed least squares** (Modat et al.,
|
||||
J. Med. Img. 1(2) 024003, 2014). Instead of a histogram MI objective it
|
||||
matches 4-voxel image blocks by **normalised cross-correlation** at each level
|
||||
of a 3-level Gaussian pyramid and solves a robust (trimmed) least-squares for
|
||||
the rigid transform — a genuinely different algorithmic family, and
|
||||
cross-modality tolerant because NCC is invariant to intensity scaling.
|
||||
|
||||
Prerequisites (check with `RP.availability()`; a missing piece raises
|
||||
`PipelineUnavailable` with the exact command):
|
||||
|
||||
```bash
|
||||
# one-time, in the benchmark env (provides reg_aladin + reg_f3d, NiftyReg 1.5.69)
|
||||
conda install -c conda-forge niftyreg
|
||||
```
|
||||
|
||||
### How the adapter drives the tool
|
||||
|
||||
1. **Rigid 6-DOF, full volumes** — `reg_aladin -ref CT -flo T1c -rigOnly`
|
||||
(symmetric variant on by default, which matters: the non-symmetric mode
|
||||
drifts grossly on these repositioned pairs). `-ln 3` pyramid,
|
||||
`-maxit 15` TLS iterations per level (the tool's default of 5 is too
|
||||
coarse), `-interp 3` cubic internal warping, default initialisation
|
||||
(image centres).
|
||||
2. **Read-back** — the output affine (`-aff`) is a 4×4 row-major matrix,
|
||||
documented as `Affine * Reference = Floating` (CT→T1c); the adapter parses
|
||||
it into a SimpleITK `AffineTransform` and runs it through the usual
|
||||
`fix_orientation` gate (skull-ring concentricity). On the M6 cases the gate
|
||||
keeps the file's affine as-is (`auto_inverted: False`), i.e. it is already
|
||||
the contract transform (T1c→CT).
|
||||
|
||||
## IMPACT pipeline (Docker) — one-time setup
|
||||
|
||||
The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
|
||||
|
|
@ -303,8 +338,10 @@ image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121)
|
|||
*semantic* similarity metric compiled in as a native Elastix metric, plus a
|
||||
TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
|
||||
Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
|
||||
**fixed = T1c, moving = CT** (the reverse of the others): swapping its roles
|
||||
breaks the TotalSegmentator metric ("Error in metric" / wrong rotation).
|
||||
**fixed = T1c, moving = CT** (the reverse of the others): the CT FOV
|
||||
(~300×300×251 mm) is much larger than the T1c crop (~220×220×150 mm), so with
|
||||
CT as the fixed image too many CT samples map outside the T1c buffer and the
|
||||
IMPACT metric aborts ("too many samples map outside moving image buffer").
|
||||
|
||||
Prerequisites (check with `RP.availability()`; a missing piece raises
|
||||
`PipelineUnavailable` with the exact command):
|
||||
|
|
@ -330,49 +367,50 @@ repo for the full list and FOV/layer tables.
|
|||
|
||||
### How the adapter drives the image
|
||||
|
||||
The image's `CMD.sh` only runs `elastix -f -m -p -out` (no `-t0`), so the
|
||||
adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`)
|
||||
to add the shared centered start:
|
||||
The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
|
||||
binary directly (`/lib/elastix-install/bin/elastix`):
|
||||
|
||||
1. **Head box** — resamples fixed=T1c and moving=CT into the same CT-space head
|
||||
box (physical frame) the other pipelines use, so MI/IMPACT is computed on
|
||||
head tissue, not neck/air.
|
||||
2. **Shared init via `-t0`** — the `centroid_offset` translation is written to a
|
||||
`TranslationTransform` parameter file and passed with `-t0`; with
|
||||
`HowToCombineTransforms "Compose"` the final transform is
|
||||
`T = T_euler(T_init(p))`. (Elastix 5.x removed in-file multi-`Transform`
|
||||
chains and `UseScalingFactors`, so the initial transform must be a separate
|
||||
`-t0` file.)
|
||||
1. **Full volumes** — fixed=T1c and moving=CT are used as-is (no head box).
|
||||
T1c must be the fixed image: the CT FOV is much larger, so with CT fixed too
|
||||
many CT samples map outside the T1c buffer and the metric aborts.
|
||||
2. **Auto init (no `-t0`)** — the parameter map sets
|
||||
`AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
|
||||
elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
|
||||
optimizer scales.
|
||||
3. **Parameter map** — rigid `EulerTransform` (6-DOF), `Impact` metric in
|
||||
**Jacobian** mode (differentiable, per the recommended config), 4 semantic
|
||||
scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the box's base spacing) with
|
||||
scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the T1c base spacing) with
|
||||
**no image resampling** (all-1 pyramid schedules, as the IMPACT docs require),
|
||||
ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid
|
||||
(`--impact-gpu -1` forces CPU).
|
||||
4. **Read-back** — the result `TransformParameters.0.itk.txt` is an
|
||||
`Euler3DTransform` in the box/world physical frame; the adapter composes it
|
||||
with the `-t0` translation and runs the usual `fix_orientation` gate.
|
||||
`Euler3DTransform` in the world physical frame (the
|
||||
`AutomaticTransformInitialization` is folded in by elastix); the adapter
|
||||
runs it through the usual `fix_orientation` gate.
|
||||
|
||||
## Notes / caveats
|
||||
|
||||
- The benchmark compares each toolkit in its **native configuration**: the
|
||||
local optimizers that cannot recover the gross offset from identity
|
||||
(`antspyx`, `impact`) get a centroid-aligned start, while `simpleitk`,
|
||||
`elastix` and `synthmorph` start from identity / their own global search.
|
||||
Production uses `--simpleitk-mode production`'s multi-start tilt grid, which
|
||||
raises SimpleITK's acceptance at ~20× the cost; that is out of scope for the
|
||||
- The benchmark compares each toolkit in its **native configuration**: all
|
||||
register on the full volumes (no head box) and use each toolkit's native
|
||||
automatic transform initialization (no manual `centroid_offset`) — elastix /
|
||||
impact `AutomaticTransformInitialization` (CenterOfGravity), antspyx's
|
||||
center-of-mass alignment, simpleitk identity, synthmorph none. Production
|
||||
uses `--simpleitk-mode production`'s multi-start tilt grid, which raises
|
||||
SimpleITK's acceptance at ~20× the cost; that is out of scope for the
|
||||
head-to-head but is the config the dataset actually ships with.
|
||||
- `synthmorph` is deterministic (one fixed pretrained model, no random
|
||||
subsampling), so its per-case numbers are stable run to run — unlike ANTs /
|
||||
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. 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`).
|
||||
IMPACT. Its low MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
|
||||
genuine cross-modality / tilt weakness, not variance.
|
||||
- `fix_orientation` (elastix/ants/niftyreg/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. 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`). `niftyreg`'s
|
||||
`-aff` file is documented as CT→T1c but the gate keeps it as-is on these
|
||||
cases (`auto_inverted: False`) — the empirical gate, not the docs, decides.
|
||||
- **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
|
||||
|
|
@ -385,8 +423,27 @@ to add the shared centered start:
|
|||
path).
|
||||
- **GPU sharing.** The container runs as the host uid with
|
||||
`TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
|
||||
- **~26 s/case** (GPU: one pretrained forward pass), between ANTs (~11 s)
|
||||
and IMPACT (~576 s).
|
||||
- **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK
|
||||
(~10 s) / ANTs (~37 s) and IMPACT (~834 s).
|
||||
- **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`,
|
||||
symmetric block matching + TLS, CPU/OpenMP):
|
||||
- **Different algorithmic family.** No histogram MI objective: 4-voxel
|
||||
blocks are matched by normalised cross-correlation and the rigid
|
||||
transform is solved by a trimmed least squares, iterated per pyramid
|
||||
level (`-ln 3 -maxit 15`). NCC's scale invariance is what makes it
|
||||
cross-modality tolerant.
|
||||
- **Symmetric mode matters.** The default symmetric variant (optimising
|
||||
forward and backward transforms concurrently) is the robust one on these
|
||||
repositioned pairs; `-noSym` drifts grossly (gap std 45 mm on a clean
|
||||
case every other tool passes).
|
||||
- **`-maxit 5` (the tool default) is too coarse** — on the clean case
|
||||
`m6_26TYQI3C` it reaches MI 0.41 vs 0.61 at 15 TLS iterations/level;
|
||||
15 is the benchmark setting.
|
||||
- **CPU-only, thread-bound** (`-omp`); ~4 min/case median (120–440 s range)
|
||||
at full volumes with 24 threads.
|
||||
- **Version note.** The upstream KCL-BMEIS release zips (v2.x) need
|
||||
glibc ≥ 2.33; this host (Ubuntu 20.04, glibc 2.31) uses the conda-forge
|
||||
build (1.5.69) instead.
|
||||
- **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python
|
||||
API — no CLI binary involved):
|
||||
- Elastix 5.x renamed components vs the classic parameter files: the rigid
|
||||
|
|
@ -394,11 +451,11 @@ to add the shared centered start:
|
|||
the component set), the metric key is `Metric` (not `ImageMetric`), the
|
||||
random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x
|
||||
names and is not compatible with elastix ≤ 4.x parameter files.
|
||||
- No initial transform is passed (identity start): the ASGD + 3-resolution
|
||||
Mattes-MI pyramid recovers the gross inter-session offset on its own, so
|
||||
the shared `centroid_offset` start is not required. An **output directory**
|
||||
is still set (`SetOutputDirectory`) for the C++ filter to run; the adapter
|
||||
uses a scratch dir under `data/m6reg/elastix/itkelastix/`.
|
||||
- `AutomaticTransformInitialization "true"` (CenterOfGravity) provides the
|
||||
start (no manual `centroid_offset`); `AutomaticScalesEstimation` sets the
|
||||
optimizer scales. An **output directory** is still set
|
||||
(`SetOutputDirectory`) for the C++ filter to run; the adapter uses a
|
||||
scratch dir under `data/m6reg/elastix/itkelastix/`.
|
||||
- `WriteResultImage` must be the string `"false"` (Python `False`); the
|
||||
default (absent) makes the filter raise "Could not read result image"
|
||||
after a *successful* registration.
|
||||
|
|
@ -408,32 +465,40 @@ to add the shared centered start:
|
|||
(`GetCombinationTransform().TransformPoint`) into an affine — independent
|
||||
of elastix's Euler-angle parameterization — and wrapped in a SimpleITK
|
||||
`AffineTransform(3)`.
|
||||
- Registration wall time is ~2.4 s/case median at full budget (ASGD
|
||||
500/500/250 over 3 resolutions).
|
||||
- Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
|
||||
(ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their
|
||||
defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32).
|
||||
- **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0,
|
||||
GPU, Jacobian mode):
|
||||
- **Slow, by design.** Every ASGD iteration runs the TotalSegmentator
|
||||
extractor on the sampled patches with backpropagation, so the recommended
|
||||
500-iter / 2000-sample budget is ~9.6 min/case (vs 2–11 s for the other
|
||||
tools). `--impact-iters` / `--impact-samples` trade accuracy for speed.
|
||||
extractor on the sampled patches with backpropagation. The benchmark
|
||||
budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 =
|
||||
the elastix default) is ~10× the authors' 500-iter / 2000-sample
|
||||
Jacobian budget, i.e. on the order of *hours* per case on full volumes —
|
||||
the results shown for `impact` are from the previous 500/2000 budget and
|
||||
were not re-run at the new budget. `--impact-iters` / `--impact-samples`
|
||||
trade accuracy for speed.
|
||||
- **Keeps fixed=T1c/moving=CT.** IMPACT is the one pipeline that does NOT
|
||||
follow the shared fixed=CT/moving=T1c convention: swapping its roles
|
||||
breaks the TotalSegmentator metric (elastix aborts with "Error in metric"
|
||||
or converges to a wrong rotation). The returned transform still satisfies
|
||||
the contract (T1c→CT).
|
||||
follow the shared fixed=CT/moving=T1c convention: the CT FOV (~300×300×251
|
||||
mm) is much larger than the T1c crop (~220×220×150 mm), so with CT as the
|
||||
fixed image too many CT samples map outside the T1c buffer and the metric
|
||||
aborts ("too many samples map outside moving image buffer"). The returned
|
||||
transform still satisfies the contract (T1c→CT).
|
||||
- **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns
|
||||
the intensity centroids; no manual `-t0` translation is used.
|
||||
- **Rugged single-init convergence.** The semantic landscape has more local
|
||||
minima than Mattes-MI: the pipeline nails the hardest repositioned case
|
||||
(`m6_22IRRKU6`, mid 8.3 mm — best of all) but lands in bad minima on some
|
||||
others. A multi-start or higher-iteration budget would likely raise its
|
||||
acceptance; that is out of scope for the single-init head-to-head.
|
||||
minima than Mattes-MI: on full volumes it no longer nails the hardest
|
||||
repositioned case (`m6_22IRRKU6`, now mid 22.7 / std 27.0 — the head box
|
||||
used to get it to mid 8.3 / std 13.2). A multi-start or higher-iteration
|
||||
budget would likely raise its acceptance; that is out of scope for the
|
||||
single-init head-to-head.
|
||||
- **Stochastic.** Jacobian mode resamples patches each iteration, so results
|
||||
vary run to run (see Stochasticity).
|
||||
- **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image
|
||||
bundles 5.2.0, so the in-file multi-`Transform` chain and `UseScalingFactors`
|
||||
are gone — the shared centered start is a separate `-t0` translation file,
|
||||
and the result is read back as an `Euler3DTransform` (`.itk.txt`) and
|
||||
composed with `-t0`. The image's own `CMD.sh` doesn't pass `-t0`, so the
|
||||
adapter invokes `/lib/elastix-install/bin/elastix` directly.
|
||||
are gone. The result is read back as an `Euler3DTransform` (`.itk.txt`).
|
||||
The image's own `CMD.sh` is a generic launcher, so the adapter invokes
|
||||
`/lib/elastix-install/bin/elastix` directly.
|
||||
- **Models are not in the image.** The TorchScript model(s) must be present
|
||||
under the mounted `/Data/Models` (the repo caches them in
|
||||
`data/impact_models/`, bind-mounted read-only).
|
||||
|
|
|
|||
|
|
@ -1,37 +1,41 @@
|
|||
"""Benchmark the 5 rigid CT->T1c registration pipelines on N M6 test cases.
|
||||
"""Benchmark the 6 rigid CT->T1c registration pipelines on N M6 test cases.
|
||||
|
||||
Pipelines (see src/registration_pipelines.py for the shared transform contract):
|
||||
simpleitk | elastix | antspyx | synthmorph | impact
|
||||
simpleitk | elastix | antspyx | synthmorph | impact | niftyreg
|
||||
|
||||
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
|
||||
the T1c native grid. Every pipeline is scored with the SAME crop-independent
|
||||
geometric QA used in production (scripts/02_build_m6_dataset.py):
|
||||
the T1c native grid. The headline alignment score is the achieved Mattes MI
|
||||
(mi_ct_mr, nats) between the warped-CT and the original MR over the MR head --
|
||||
higher = better in-sync = better aligned. The production geometric QA
|
||||
(scripts/02_build_m6_dataset.py) is still recorded per case in the JSON:
|
||||
|
||||
- skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR
|
||||
head silhouette (mid/low level + angular std).
|
||||
Aligned = 0 <= mid <= 20 mm and std <= 12 mm.
|
||||
- mi_ct_mr : Mattes MI (nats), warped-CT ~ MR over the MR head
|
||||
(the reported alignment score).
|
||||
- tv_in_ct_frac : fraction of the (dilated) warped GTV that falls inside
|
||||
the warped-CT head tissue (the slab must straddle it).
|
||||
- tv_nvox : warped GTV volume (>0.5 voxels).
|
||||
- skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR
|
||||
head silhouette (mid/low level + angular std).
|
||||
- accepted : production rule (nvox>0 and tv_in_ct and concentric).
|
||||
- reg_wall_s : wall-clock registration time (excludes I/O + scoring).
|
||||
|
||||
Because some toolkits are not pip-installable into the main
|
||||
`longitudinal-313` env (antspyx is kept in an isolated env, see
|
||||
Because some toolkits are not pip-installable into every
|
||||
`longitudinal-313` env (antspyx needs a cp310-cp313 wheel, see
|
||||
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
|
||||
ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively
|
||||
(``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
|
||||
``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)
|
||||
# main env (SimpleITK + elastix + antspyx + niftyreg)
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg
|
||||
|
||||
# synthmorph / impact (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md
|
||||
python scripts/bench_registration.py --pipelines synthmorph impact
|
||||
|
|
@ -232,6 +236,8 @@ def _call(pipeline, ct, t1c, tv, args):
|
|||
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)
|
||||
if pipeline == "niftyreg":
|
||||
return RP.register_niftyreg(ct, t1c, threads=args.threads)
|
||||
raise KeyError(pipeline)
|
||||
|
||||
|
||||
|
|
@ -255,10 +261,7 @@ def _run_case(pipeline, case, args):
|
|||
out.update(meta)
|
||||
out.update(evaluate(ct, t1c, tv, T))
|
||||
if getattr(args, "qa_images", False):
|
||||
g = out.get("skull_ring_gap_mm")
|
||||
title = (f"{pipeline} | {key} | MI={out.get('mi_ct_mr')} | "
|
||||
f"gap mid {g['mid_mm']}mm std {g['std_mm']}mm" if g
|
||||
else f"{pipeline} | {key}")
|
||||
title = f"{pipeline} | {key} | MI={out.get('mi_ct_mr')}"
|
||||
qp = os.path.join(args.qa_out, f"{pipeline}__{key}.png")
|
||||
try:
|
||||
wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
|
||||
|
|
@ -287,12 +290,8 @@ def run_pipeline(pipeline, cases, args):
|
|||
try:
|
||||
results.append(fu.result())
|
||||
r = results[-1]
|
||||
g = r.get("skull_ring_gap_mm")
|
||||
reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT"
|
||||
print(f" {c['key']}: accepted={r['accepted']} "
|
||||
f"gap={g['mid_mm'] if g else 'None'}mm "
|
||||
f"std={g['std_mm'] if g else 'None'}mm "
|
||||
f"mi={r.get('mi_ct_mr')} "
|
||||
print(f" {c['key']}: mi={r.get('mi_ct_mr')} "
|
||||
f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True)
|
||||
except Exception as e: # noqa
|
||||
errors.append({"key": c["key"], "error": repr(e)})
|
||||
|
|
@ -348,27 +347,24 @@ def _fmt(v, spec="{}", default="-"):
|
|||
|
||||
def print_table(rows):
|
||||
order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)]
|
||||
print("\n" + "=" * 100)
|
||||
print("\n" + "=" * 92)
|
||||
hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} "
|
||||
f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
|
||||
f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
|
||||
print(hdr)
|
||||
print("-" * 100)
|
||||
print("-" * 92)
|
||||
for p in order:
|
||||
r = next(x for x in rows if x["pipeline"] == p)
|
||||
if not r.get("available"):
|
||||
print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} "
|
||||
f"{'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]")
|
||||
f"{'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]")
|
||||
continue
|
||||
gm = (r.get("skull_ring_gap_mid_mm") or {}).get("median")
|
||||
gs = (r.get("skull_ring_gap_std_mm") or {}).get("median")
|
||||
rt = (r.get("reg_wall_s") or {}).get("median")
|
||||
tv = (r.get("tv_in_ct_frac") or {}).get("mean")
|
||||
mi = (r.get("mi_ct_mr") or {}).get("median")
|
||||
print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} "
|
||||
f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {r['accepted_rate']:>7.2f} "
|
||||
f"{_fmt(gm, '{:.1f}'):>7s} {_fmt(gs, '{:.1f}'):>7s} {_fmt(rt, '{:.1f}'):>7s} "
|
||||
f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} "
|
||||
f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
|
||||
print("=" * 100)
|
||||
print("=" * 92)
|
||||
|
||||
|
||||
def write_plot(rows, out_png):
|
||||
|
|
@ -384,19 +380,13 @@ def write_plot(rows, out_png):
|
|||
print("no available pipelines to plot")
|
||||
return
|
||||
labels = list(ok_rows)
|
||||
gap_mid = [ok_rows[p].get("skull_ring_gap_mid_mm", {}).get("median") or 0 for p in labels]
|
||||
gap_std = [ok_rows[p].get("skull_ring_gap_std_mm", {}).get("median") or 0 for p in labels]
|
||||
regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels]
|
||||
acc = [ok_rows[p].get("accepted_rate") or 0 for p in labels]
|
||||
mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels]
|
||||
fig, axes = plt.subplots(1, 5, figsize=(20, 4))
|
||||
colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD", "#E08E45"]
|
||||
for ax, (title, vals, fmt, ylab, inv) in zip(axes, [
|
||||
("Skull-ring gap mid (mm, median)", gap_mid, "{:.1f}", "mm", True),
|
||||
("Skull-ring gap std (mm, median)", gap_std, "{:.1f}", "mm", True),
|
||||
("Registration time (s, median)", regt, "{:.0f}", "s", True),
|
||||
("Acceptance rate", acc, "{:.2f}", "frac", False),
|
||||
("Mean Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats", False)]):
|
||||
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
|
||||
colors = ["#4878CF", "#6ACC65"]
|
||||
for ax, (title, vals, fmt, ylab) in zip(axes, [
|
||||
("Registration time (s, median)", regt, "{:.0f}", "s"),
|
||||
("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]):
|
||||
bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)])
|
||||
for b, v in zip(bars, vals):
|
||||
ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()),
|
||||
|
|
@ -426,10 +416,11 @@ def merge(args):
|
|||
print_table(rows)
|
||||
summary = {
|
||||
"timestamp": datetime.now().isoformat(timespec="seconds"),
|
||||
"protocol": ("rigid CT->T1c; cross-modality; scored by production skull-ring "
|
||||
"concentricity (aligned: 0<=mid<=20mm, std<=12mm) + GTV-in-CT fraction "
|
||||
"(>=0.5); reg time excludes I/O+scoring; per-pipeline best-of-T vs T^-1 "
|
||||
"orientation gate for elastix/antspyx"),
|
||||
"protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI "
|
||||
"(nats, warped-CT ~ MR over the MR head, higher = better); reg time "
|
||||
"excludes I/O+scoring; per-case JSON retains the production skull-ring "
|
||||
"gap / acceptance QA; best-of-T vs T^-1 orientation gate for the "
|
||||
"elastix-style toolkits"),
|
||||
"n_cases": next((r["n_cases"] for r in rows if r.get("available")), None),
|
||||
"case_keys": sorted({c["key"] for r in rows if r.get("available") for c in r.get("cases", [])}),
|
||||
"results": rows,
|
||||
|
|
@ -440,18 +431,16 @@ def merge(args):
|
|||
print("summary:", sf)
|
||||
write_plot(rows, os.path.join(out_dir, "bench_registration.png"))
|
||||
# markdown table for the report
|
||||
md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT | MI (nats) |",
|
||||
"|---|---|---|---|---|---|---|---|---|---|---|---|"]
|
||||
md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | tv-in-CT | MI (nats, median) |",
|
||||
"|---|---|---|---|---|---|---|---|---|"]
|
||||
for p in [x["pipeline"] for x in rows]:
|
||||
r = next(x for x in rows if x["pipeline"] == p)
|
||||
if not r.get("available"):
|
||||
md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | - | - | - | ({r['reason']}) |")
|
||||
md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | ({r['reason']}) |")
|
||||
continue
|
||||
md.append(
|
||||
f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | "
|
||||
f"{r['n_error']} | {r['accepted_rate']:.2f} | "
|
||||
f"{(r.get('skull_ring_gap_mid_mm') or {}).get('median')} | "
|
||||
f"{(r.get('skull_ring_gap_std_mm') or {}).get('median')} | "
|
||||
f"{r['n_error']} | "
|
||||
f"{(r.get('reg_wall_s') or {}).get('median')} | "
|
||||
f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | "
|
||||
f"{(r.get('mi_ct_mr') or {}).get('median')} |")
|
||||
|
|
@ -492,12 +481,12 @@ def parse_args():
|
|||
ap.add_argument("--impact-gpu", type=int, default=None,
|
||||
help="GPU device index for IMPACT (default: auto-spread over "
|
||||
"the host GPUs by pid; -1 forces CPU)")
|
||||
ap.add_argument("--impact-iters", type=int, default=500,
|
||||
help="IMPACT per-resolution max iterations (default 500, the authors' "
|
||||
"recommended Jacobian budget; lower to trade accuracy for speed)")
|
||||
ap.add_argument("--impact-samples", type=int, default=2000,
|
||||
help="IMPACT random spatial samples per iteration (default 2000, "
|
||||
"the recommended budget)")
|
||||
ap.add_argument("--impact-iters", type=int, default=2000,
|
||||
help="IMPACT MaximumNumberOfIterations (default 2000, the benchmark "
|
||||
"budget; ~10x the authors' 500-iter Jacobian budget, much slower)")
|
||||
ap.add_argument("--impact-samples", type=int, default=5000,
|
||||
help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
|
||||
"elastix default)")
|
||||
ap.add_argument("--qa-images", action="store_true",
|
||||
help="also write a checkreg-style PNG (checkerboard blend + color "
|
||||
"overlay of the inverse-warped CT vs the original MR) per case")
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
"""Five rigid CT -> T1c registration pipelines behind one unified contract.
|
||||
"""Six rigid CT -> T1c registration pipelines behind one unified contract.
|
||||
|
||||
Problem (shared with ``scripts/02_build_m6_dataset.py``)
|
||||
--------------------------------------------------------
|
||||
|
|
@ -23,34 +23,41 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
|
|||
|
||||
Pipelines (listed alphabetically)
|
||||
---------------------------------
|
||||
All pipelines use the shared convention **fixed = CT, moving = T1c**, except
|
||||
IMPACT which keeps **fixed = T1c, moving = CT** (swapping breaks its
|
||||
TotalSegmentator metric -- see ``register_impact``). A pipeline that is a
|
||||
local optimizer and cannot recover the gross inter-session offset from
|
||||
identity (antspyx, impact) starts from a centroid-aligned translation; the
|
||||
others (simpleitk, elastix, synthmorph) do not require it and start from
|
||||
identity / their own global search.
|
||||
All pipelines register on the **full volumes** (no head box) and use each
|
||||
toolkit's **native automatic transform initialization** (no manual
|
||||
``centroid_offset``): elastix / impact use elastix's
|
||||
``AutomaticTransformInitialization`` (CenterOfGravity), antspyx uses its
|
||||
default center-of-mass alignment, simpleitk starts from identity, niftyreg
|
||||
(reg_aladin) initializes from the image centres, and synthmorph needs none
|
||||
(global deep method). All use **fixed = CT, moving = T1c** except IMPACT,
|
||||
which keeps **fixed = T1c, moving = CT** (swapping breaks its
|
||||
TotalSegmentator metric -- see ``register_impact``).
|
||||
|
||||
- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
|
||||
Centroid pre-shift of the moving (T1c) image is *required* (ANTs
|
||||
Rigid-from-identity collapses on the gross offset).
|
||||
Default ``initial_transform=None`` computes a translation aligning the image
|
||||
centers of mass.
|
||||
- ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` +
|
||||
``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
|
||||
Identity start (the ASGD + multi-resolution pyramid recovers the offset).
|
||||
``AutomaticTransformInitialization`` (CenterOfGravity) start.
|
||||
- ``register_impact`` : IMPACT (semantic-similarity metric) inside the
|
||||
``vboussot/elastix_impact`` Docker image (Elastix 5.2 + a TorchScript
|
||||
TotalSegmentator feature extractor as a native Elastix metric). Rigid 6-DOF,
|
||||
GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); centroid init
|
||||
kept (local ASGD on the semantic metric benefits from the gross alignment).
|
||||
GPU (Jacobian mode). Kept as fixed=T1c/moving=CT (see above); elastix
|
||||
``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
|
||||
Needs Docker + the image + a model asset.
|
||||
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a
|
||||
CT-space head box. Reuses the *proven* production ``register_ct_to_t1c`` in
|
||||
- ``register_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
|
||||
symmetric block-matching + trimmed least squares (Modat et al., J. Med.
|
||||
Img. 2014), NCC-based block similarity (cross-modality tolerant). Rigid
|
||||
6-DOF (``-rigOnly``), 3-level pyramid, 15 TLS iterations/level. Needs
|
||||
``conda install -c conda-forge niftyreg`` (provides ``reg_aladin``).
|
||||
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on the
|
||||
full volumes. Reuses the *proven* production ``register_ct_to_t1c`` in
|
||||
``mode="production"``; the benchmark mode is a single identity-start
|
||||
multi-resolution pyramid run.
|
||||
- ``register_synthmorph`` : FreeSurfer SynthMorph (TensorFlow, pretrained
|
||||
anatomy-aware contrast-invariant model) inside the ``freesurfer/synthmorph``
|
||||
Docker image. Rigid 6-DOF matrix model, GPU when available. Global deep
|
||||
method -- no centroid init. Needs Docker + the pulled image.
|
||||
method -- no init. Needs Docker + the pulled image.
|
||||
|
||||
A toolkit that is not installed (or, for the Docker-based pipelines, missing
|
||||
its image)
|
||||
|
|
@ -70,6 +77,7 @@ applied to elastix / antspyx / synthmorph / impact, not to the trusted
|
|||
production SimpleITK path.
|
||||
"""
|
||||
import os
|
||||
import re
|
||||
import shutil
|
||||
import subprocess
|
||||
import importlib.util
|
||||
|
|
@ -157,38 +165,18 @@ def centroid_offset(ct, t1c, ct_box=None):
|
|||
return (float(t1c_c[0] - c[0]), float(t1c_c[1] - c[1]), float(t1c_c[2] - c[2]))
|
||||
|
||||
|
||||
def _ct_head_box_mm(ct):
|
||||
"""CT-space head box in physical mm: (box_lo, box_hi) tuples, or None."""
|
||||
m6 = _load_m6()
|
||||
ca = sitk.GetArrayFromImage(ct)
|
||||
tissue = ca > -150
|
||||
csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
|
||||
vz = m6._ct_vertex_slice(ct)
|
||||
if vz is None:
|
||||
return None, None
|
||||
ztop = co[2] + vz * csp[2]
|
||||
band = tissue[max(0, vz - int(round(150 / csp[2]))):vz + 1]
|
||||
cols = np.nonzero(band.any(0))
|
||||
box_lo = (co[0] + cols[1].min() * csp[0] - 15,
|
||||
co[1] + cols[0].min() * csp[1] - 15, ztop - 215)
|
||||
box_hi = (co[0] + cols[1].max() * csp[0] + 15,
|
||||
co[1] + cols[0].max() * csp[1] + 15, ztop + 10)
|
||||
return box_lo, box_hi
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# 2. SimpleITK (production baseline -- reused verbatim)
|
||||
# --------------------------------------------------------------------------- #
|
||||
|
||||
def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
|
||||
"""SimpleITK rigid CT->T1c, Mattes MI on a CT-space head box.
|
||||
def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
|
||||
"""SimpleITK rigid CT->T1c, Mattes MI, full volumes (no head box).
|
||||
|
||||
mode="benchmark" (default): a single multi-resolution pyramid run from the
|
||||
identity start, capped at ``iters_full`` iterations at full resolution.
|
||||
This is the *standard* SimpleITK algorithm and matches the budget the
|
||||
other toolkits get in the benchmark (one rigid MI registration), so the
|
||||
pipelines are compared on equal footing. The pyramid recovers the gross
|
||||
inter-session offset, so no centroid init is used.
|
||||
identity start, with ``iters_full`` (=2000) optimizer iterations at every
|
||||
pyramid level. This is the *standard* SimpleITK algorithm and matches the
|
||||
2000-iteration budget the other toolkits get in the benchmark, so the
|
||||
pipelines are compared on equal footing. No centroid init is used.
|
||||
mode="production": the exact multi-start tilt-grid ``register_ct_to_t1c``
|
||||
from scripts/02_build_m6_dataset.py (robust to the ~30 deg inter-session
|
||||
tilt but much slower -- not the default for the benchmark).
|
||||
|
|
@ -215,20 +203,16 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
|
|||
except AttributeError:
|
||||
pass
|
||||
ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
|
||||
lo, hi = _ct_head_box_mm(ct)
|
||||
if lo is None:
|
||||
raise PipelineUnavailable("simpleitk: no CT head box found (corrupt CT?)")
|
||||
tsp = np.array(t1c.GetSpacing())
|
||||
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
|
||||
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
|
||||
ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
|
||||
ref.SetOrigin(tuple(float(x) for x in lo))
|
||||
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
|
||||
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
|
||||
# Identity start: the multi-resolution Mattes-MI pyramid recovers the gross
|
||||
# inter-session offset on its own, so no centroid init is required here.
|
||||
ct_c = np.array([(lo[i] + hi[i]) / 2.0 for i in range(3)])
|
||||
init = sitk.Euler3DTransform(); init.SetCenter(tuple(float(x) for x in ct_c))
|
||||
# full volumes, fixed = CT, moving = T1c (no head box)
|
||||
fixed, moving = ct, t1c
|
||||
shape = np.array(ct.GetSize())
|
||||
# Identity start (no centroid init); the multi-resolution Mattes-MI pyramid
|
||||
# recovers the gross inter-session offset. The Euler center is the fixed
|
||||
# (CT) geometric center.
|
||||
init = sitk.Euler3DTransform()
|
||||
init.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint(
|
||||
((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2,
|
||||
(fixed.GetSize()[2] - 1) // 2))))
|
||||
init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
|
||||
import time
|
||||
t0 = time.time()
|
||||
|
|
@ -250,9 +234,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
|
|||
meta = {
|
||||
"pipeline": "simpleitk", "mode": "benchmark",
|
||||
"toolkit": f"SimpleITK {sitk.Version_VersionString()}",
|
||||
"metric": "MattesMutualInformation (head box, CT space)",
|
||||
"metric": "MattesMutualInformation (full volumes)",
|
||||
"transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
|
||||
f"{iters_full} iters/level)"),
|
||||
f"{iters_full} iters/level, identity start)"),
|
||||
"convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
|
||||
"seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
|
||||
"skull_ring_gap_mm": gap,
|
||||
|
|
@ -311,34 +295,29 @@ def _itk_transform_to_sitk_rigid(itk_tf):
|
|||
|
||||
|
||||
def register_elastix(ct, t1c, workdir=None, threads=8):
|
||||
"""ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), Mattes MI.
|
||||
"""ITK-Elastix rigid CT->T1c (fixed=CT, moving=T1c), Mattes MI, full images.
|
||||
|
||||
Runs in-process through the ITK-Elastix python API (pip ``itk`` +
|
||||
``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the
|
||||
5.x replacement for ``Rigid3DTransform``) + Mattes MI
|
||||
(``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common
|
||||
CT-space head box, from the identity start (the ASGD + multi-resolution
|
||||
pyramid recovers the gross inter-session offset, so no centroid init is
|
||||
required). Elastix 5.x component names: sampler ``Random`` (not
|
||||
``RandomSampler``), metric key ``Metric`` (not ``ImageMetric``). An output
|
||||
directory is required for the C++ filter to run, and ``WriteResultImage``
|
||||
must be the string ``"false"`` (a Python bool False) or the filter raises
|
||||
"Could not read result image" after a successful registration.
|
||||
(``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the full
|
||||
CT / T1c volumes (no head box). No manual centroid init: elastix's
|
||||
``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity
|
||||
centroids, and ``AutomaticScalesEstimation`` sets the optimizer scales.
|
||||
Budget: ``MaximumNumberOfIterations`` 2000 per resolution; the sampler /
|
||||
metric run at their defaults (``NumberOfSpatialSamples`` 5000,
|
||||
``NumberOfHistogramBins`` 32).
|
||||
Elastix 5.x
|
||||
component names: sampler ``Random`` (not ``RandomSampler``), metric key
|
||||
``Metric`` (not ``ImageMetric``). An output directory is required for the
|
||||
C++ filter to run, and ``WriteResultImage`` must be the string ``"false"``
|
||||
(a Python bool False) or the filter raises "Could not read result image"
|
||||
after a successful registration.
|
||||
"""
|
||||
itk = _elastix_api()
|
||||
lo, hi = _ct_head_box_mm(ct)
|
||||
if lo is None:
|
||||
raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)")
|
||||
# Resample both into a common CT-space head box (mirrors the other pipelines)
|
||||
# so MI is computed on head tissue, not neck/air.
|
||||
tsp = np.array(t1c.GetSpacing())
|
||||
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
|
||||
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
|
||||
ref.SetSpacing(t1c.GetSpacing())
|
||||
ref.SetDirection(t1c.GetDirection())
|
||||
ref.SetOrigin(tuple(float(x) for x in lo))
|
||||
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = reference
|
||||
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = moving
|
||||
# full volumes, fixed = CT, moving = T1c (no head box)
|
||||
fixed = ct
|
||||
moving = t1c
|
||||
params = {
|
||||
"ImageDimension": 3,
|
||||
"FixedInternalImagePixelType": "float",
|
||||
|
|
@ -348,10 +327,13 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
|
|||
"FixedCenteredTransform": True,
|
||||
"TransformParameterReadoutFormat": 1,
|
||||
"HowToCombineTransforms": "Compose",
|
||||
"AutomaticTransformInitialization": True,
|
||||
"AutomaticTransformInitializationMethod": "CenterOfGravity",
|
||||
"AutomaticScalesEstimation": True,
|
||||
"Optimizer": "AdaptiveStochasticGradientDescent",
|
||||
"MaximumNumberOfIterations": [500, 500, 250],
|
||||
"MaximumNumberOfIterations": [2000, 2000, 2000],
|
||||
"Metric": "AdvancedMattesMutualInformation",
|
||||
"NumberOfHistogramBins": 10,
|
||||
"NumberOfHistogramBins": 32,
|
||||
"ImageSampler": "Random",
|
||||
"ImageSamplerPatchSize": [25, 25, 25],
|
||||
"WriteResultImage": False,
|
||||
|
|
@ -388,8 +370,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
|
|||
version = "ITK-Elastix (python API)"
|
||||
meta = {
|
||||
"pipeline": "elastix", "toolkit": version,
|
||||
"metric": "MattesMutualInformation (head box, CT space)",
|
||||
"transform": "Rigid (EulerTransform, ASGD, 3 resolutions, identity start)",
|
||||
"metric": "MattesMutualInformation (full volumes)",
|
||||
"transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)",
|
||||
"convention": "T1c->CT (verified)", "auto_inverted": inverted,
|
||||
"fixed": "CT", "moving": "T1c", "seconds_reg": round(dt, 1),
|
||||
}
|
||||
|
|
@ -402,18 +384,16 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
|
|||
|
||||
def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
|
||||
aff_shrink_factors=None, aff_smoothing_sigmas=None,
|
||||
aff_metric="mattes", use_head_box=True, center_init=True):
|
||||
"""ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c).
|
||||
aff_metric="mattes"):
|
||||
"""ANTs rigid CT->T1c via antspyx (fixed=CT, moving=T1c), full volumes.
|
||||
|
||||
antspyx 0.6.3: ``ants.registration(fixed, moving, type_of_transform='Rigid',
|
||||
aff_metric='mattes', aff_iterations=..., aff_shrink_factors=...,
|
||||
aff_smoothing_sigmas=...)``. ``aff_metric='mattes'`` (the default) gives a
|
||||
multi-resolution Mattes-MI rigid registration -- the cross-modality metric
|
||||
this CT/MR pair needs. To keep the benchmark tractable and consistent with
|
||||
the other pipelines (and to avoid neck/air in the metric), both images are
|
||||
first resampled into the same CT-space head box. ANTs Rigid-from-identity
|
||||
does NOT recover the gross inter-session offset, so a centroid pre-shift of
|
||||
the moving (T1c) image is required (see ``center_init``). Returns the
|
||||
this CT/MR pair needs. Registers on the full CT / T1c volumes (no head box).
|
||||
No manual centroid init: antspyx's default ``initial_transform=None``
|
||||
computes a translation that aligns the image centers of mass. Returns the
|
||||
transform that warps CT into the T1c grid (orientation verified by
|
||||
:func:`fix_orientation`).
|
||||
"""
|
||||
|
|
@ -428,36 +408,12 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
|
|||
aff_iterations = tuple(aff_iterations) if aff_iterations else (500, 500, 200)
|
||||
aff_shrink_factors = tuple(aff_shrink_factors) if aff_shrink_factors else (4, 2, 1)
|
||||
aff_smoothing_sigmas = tuple(aff_smoothing_sigmas) if aff_smoothing_sigmas else (3, 2, 0)
|
||||
if use_head_box:
|
||||
lo, hi = _ct_head_box_mm(ct)
|
||||
if lo is None:
|
||||
raise PipelineUnavailable("antspyx: no CT head box found (corrupt CT?)")
|
||||
tsp = np.array(t1c.GetSpacing())
|
||||
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
|
||||
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
|
||||
ref.SetSpacing(t1c.GetSpacing()); ref.SetDirection(t1c.GetDirection())
|
||||
ref.SetOrigin(tuple(float(x) for x in lo))
|
||||
fx = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT box (fixed)
|
||||
mv = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c box (moving)
|
||||
else:
|
||||
# full volumes, fixed = CT, moving = T1c (no head box)
|
||||
fx, mv = ct, t1c
|
||||
import time
|
||||
t0 = time.time()
|
||||
# Centered start (REQUIRED for ANTs): antspyx cannot use a
|
||||
# TranslationTransform as initial_transform (the ANTs CLI rejects it), and
|
||||
# Rigid-from-identity does not recover the large (~100-130 mm) inter-session
|
||||
# offset -- without this init ANTs collapses on 8/9 cases. So we pre-shift the
|
||||
# T1c (moving) box so its head centroid coincides with the CT (fixed) head
|
||||
# centroid, run Rigid from identity, then compose the shift + inverse back so
|
||||
# the result is expressed in the original T1c/CT coordinates.
|
||||
off = None
|
||||
if center_init:
|
||||
off = centroid_offset(ct, t1c, ct_box=fx)
|
||||
if off is not None:
|
||||
# shift the T1c content by -off (head c_t1c -> c_ct): resample with +off
|
||||
tr_off = sitk.TranslationTransform(3)
|
||||
tr_off.SetOffset((off[0], off[1], off[2]))
|
||||
mv = sitk.Resample(mv, mv, tr_off, sitk.sitkLinear, 0.0)
|
||||
# No manual centroid init: antspyx's default initial_transform=None computes
|
||||
# a translation aligning the image centers of mass.
|
||||
af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32))
|
||||
am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32))
|
||||
r = ants.registration(fixed=af, moving=am, type_of_transform=transform_type,
|
||||
|
|
@ -466,20 +422,15 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
|
|||
aff_smoothing_sigmas=aff_smoothing_sigmas)
|
||||
dt = time.time() - t0
|
||||
# r['fwdtransforms'][0] is the ANTs forward transform: fixed(CT) -> moving
|
||||
# (pre-shifted T1c). We need the contract T1c(orig) -> CT: map T1c(orig) to
|
||||
# the pre-shifted T1c (translation by -off), then invert the ANTs transform.
|
||||
# (T1c). The contract is T1c -> CT, i.e. its inverse; fix_orientation
|
||||
# verifies the direction against the skull ring.
|
||||
T_ants = sitk.ReadTransform(r["fwdtransforms"][0])
|
||||
if off is not None:
|
||||
tr_shift = sitk.TranslationTransform(3)
|
||||
tr_shift.SetOffset((-off[0], -off[1], -off[2])) # T1c orig -> pre
|
||||
cand = sitk.CompositeTransform([tr_shift, T_ants.GetInverse()])
|
||||
else:
|
||||
cand = T_ants.GetInverse()
|
||||
T, inverted = fix_orientation(ct, t1c, cand)
|
||||
meta = {
|
||||
"pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
|
||||
"metric": f"ANTs {aff_metric} (multi-resolution)",
|
||||
"transform": f"ANTs {transform_type} (head box)", "convention": "T1c->CT (verified)",
|
||||
"metric": f"ANTs {aff_metric} (multi-resolution, full volumes)",
|
||||
"transform": f"ANTs {transform_type} (auto COM init)", "convention": "T1c->CT (verified)",
|
||||
"auto_inverted": inverted, "aff_iterations": list(aff_iterations),
|
||||
"aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
|
||||
}
|
||||
|
|
@ -740,12 +691,18 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
|
|||
"(ImpactFeaturesMapUpdateInterval -1)",
|
||||
'(ImpactWriteFeatureMaps "false")',
|
||||
"",
|
||||
"// rigid 6-DOF to be optimized (the shared centered start comes via -t0)",
|
||||
"// rigid 6-DOF to be optimized (auto center init, no -t0)",
|
||||
'(Transform "EulerTransform")',
|
||||
'(TransformBase "ElastixTransformBase")',
|
||||
'(TransformParameterReadoutFormat "AllParametersOneLine")',
|
||||
'(FixedCenteredTransform "true")',
|
||||
'(HowToCombineTransforms "Compose")',
|
||||
'(AutomaticTransformInitialization "true")',
|
||||
# align intensity centroids (heads), not buffer corners -- the CT FOV is
|
||||
# much larger than the T1c crop, so GeometricalCenter leaves the heads
|
||||
# ~100 mm apart and the metric aborts ("too many samples outside buffer")
|
||||
'(AutomaticTransformInitializationMethod "CenterOfGravity")',
|
||||
'(AutomaticScalesEstimation "true")',
|
||||
"",
|
||||
"(ImageDimension 3)",
|
||||
'(FixedInternalImagePixelType "float")',
|
||||
|
|
@ -788,19 +745,12 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
|
|||
return "\n".join(L) + "\n"
|
||||
|
||||
|
||||
def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
||||
def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
||||
threads=12, keep_scratch=False, scales=(6.0, 3.0, 1.5, 1.0),
|
||||
downsample=False):
|
||||
"""IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image
|
||||
``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
|
||||
|
||||
NOTE on the fixed/moving convention: every other pipeline here uses
|
||||
fixed=CT, moving=T1c. IMPACT is the exception -- swapping its fixed/moving
|
||||
to CT/T1c breaks the TotalSegmentator-based semantic metric (elastix aborts
|
||||
with "Error in metric" or converges to a wrong rotation), so IMPACT keeps
|
||||
the working MRI-as-fixed (T1c) / CT-as-moving convention. The returned
|
||||
transform still satisfies the unified contract (T1c->CT).
|
||||
|
||||
IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
|
||||
one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run
|
||||
on randomly sampled patches of the fixed and moving images and the
|
||||
|
|
@ -809,21 +759,26 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
the Mattes-MI pipelines above -- it is the cross-modality tool the image was
|
||||
built for.
|
||||
|
||||
To stay on the benchmark's equal footing the pipeline:
|
||||
- registers on the same CT-space head box (both images resampled into it),
|
||||
- starts from the *shared* centroid-aligned translation (passed to Elastix
|
||||
as a fixed initial transform via ``-t0``; with ``HowToCombineTransforms
|
||||
Compose`` the final transform = ``T_main(T_init(p))``),
|
||||
- uses a rigid 6-DOF ``EulerTransform`` (the 5.x name for the classic
|
||||
3D rigid), ASGD, a 4-level scale pyramid (no image resampling).
|
||||
Budget: the default ``iters=500`` / ``samples=2000`` is the authors'
|
||||
recommended Jacobian budget. It is *much* slower than the intensity-metric
|
||||
pipelines (minutes per case on GPU vs seconds) because every iteration runs
|
||||
the TotalSegmentator feature extractor on the sampled patches with
|
||||
backpropagation. Lower ``iters``/``samples`` to trade accuracy for speed.
|
||||
The pipeline registers on the full CT / T1c volumes (no head box) with no
|
||||
manual centroid init: elastix's ``AutomaticTransformInitialization``
|
||||
(CenterOfGravity) aligns the intensity centroids. It uses a rigid 6-DOF
|
||||
``EulerTransform`` (the 5.x name for the classic 3D rigid), ASGD, a 4-level
|
||||
scale pyramid (no image resampling).
|
||||
|
||||
NOTE: IMPACT keeps **fixed=T1c, moving=CT** (the reverse of the other
|
||||
pipelines). The CT FOV (~300x300x251 mm) is much larger than the T1c crop
|
||||
(~220x220x150 mm); with CT as the fixed image, too many CT samples map
|
||||
outside the T1c buffer and the IMPACT metric aborts ("too many samples map
|
||||
outside moving image buffer"). The returned transform still satisfies the
|
||||
unified contract (T1c->CT).
|
||||
Budget: the default ``iters=2000`` (MaximumNumberOfIterations) /
|
||||
``samples=5000`` (NumberOfSpatialSamples, the elastix default) runs the
|
||||
full benchmark budget for best results -- roughly 10x the authors'
|
||||
500/2000 Jacobian budget, i.e. hours per case on full volumes. Lower
|
||||
``iters``/``samples`` to trade accuracy for speed.
|
||||
The result is read back from ``TransformParameters.0.itk.txt``
|
||||
(an ``Euler3DTransform`` in the box/world physical frame) and composed with
|
||||
the ``-t0`` translation, then checked by :func:`fix_orientation`.
|
||||
(an ``Euler3DTransform`` in the world physical frame) and checked by
|
||||
:func:`fix_orientation`.
|
||||
|
||||
Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, the
|
||||
``vboussot/elastix_impact`` image, and the model asset under the cache. A
|
||||
|
|
@ -837,25 +792,13 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
raise PipelineUnavailable(f"impact: {note}")
|
||||
model_host, model_rel = _impact_model_host(model)
|
||||
|
||||
lo, hi = _ct_head_box_mm(ct)
|
||||
if lo is None:
|
||||
raise PipelineUnavailable("impact: no CT head box found (corrupt CT?)")
|
||||
|
||||
# resample both into the shared CT-space head box (mirrors the other pipelines)
|
||||
tsp = np.array(t1c.GetSpacing())
|
||||
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
|
||||
ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
|
||||
ref.SetSpacing(t1c.GetSpacing())
|
||||
ref.SetDirection(t1c.GetDirection())
|
||||
ref.SetOrigin(tuple(float(x) for x in lo))
|
||||
fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference
|
||||
moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # CT = moving
|
||||
|
||||
# shared centered start (REQUIRED: the local ASGD optimizer does not recover
|
||||
# the gross offset from identity). Elastix maps fixed(CT)->moving(T1c), the
|
||||
# same direction as centroid_offset, so no negation.
|
||||
off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0)
|
||||
init = (-float(off[0]), -float(off[1]), -float(off[2]))
|
||||
# full volumes, fixed = T1c, moving = CT (no head box, no manual init --
|
||||
# elastix's AutomaticTransformInitialization aligns the image centers).
|
||||
# T1c must be the fixed image: the CT FOV (~300x300x251 mm) is much larger
|
||||
# than the T1c crop (~220x220x150 mm), so with CT fixed, too many CT samples
|
||||
# map outside the T1c buffer and the IMPACT metric aborts.
|
||||
fixed = t1c
|
||||
moving = ct
|
||||
|
||||
# GPU: auto-spread across the host GPUs by pid, or forced value
|
||||
if gpu is None:
|
||||
|
|
@ -901,13 +844,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
|
||||
gpu, iters, samples, voxel_sizes,
|
||||
rescale=schedule, smoothing=smoothe))
|
||||
with open(os.path.join(datadir, "InitialTransform.txt"), "w") as f:
|
||||
f.write(
|
||||
'(Transform "TranslationTransform")\n'
|
||||
'(TransformBase "ElastixTransformBase")\n'
|
||||
'(TransformParameterReadoutFormat "AllParametersOneLine")\n'
|
||||
"(NumberOfParameters 3)\n"
|
||||
f"(TransformParameters {init[0]:g} {init[1]:g} {init[2]:g})\n")
|
||||
|
||||
cmd = ["docker", "run", "--rm"]
|
||||
if gpu >= 0:
|
||||
|
|
@ -919,7 +855,6 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
IMPACT_IMAGE, IMPACT_ELASTIX_BIN,
|
||||
"-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha",
|
||||
"-p", "/Data/ParameterMap.txt",
|
||||
"-t0", "/Data/InitialTransform.txt",
|
||||
"-out", "/Out/", "-threads", str(threads),
|
||||
]
|
||||
t0 = time.time()
|
||||
|
|
@ -934,14 +869,12 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
|
||||
raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
|
||||
|
||||
# the optimized Euler (box/world frame) composed with the fixed -t0
|
||||
# translation: T = T_euler(T_init(p)). Both map fixed(CT)->moving(T1c),
|
||||
# i.e. the inverse of the contract (T1c->CT); fix_orientation resolves it.
|
||||
# the optimized Euler (world frame) is the full combination transform
|
||||
# (AutomaticTransformInitialization is folded in by elastix); it maps
|
||||
# fixed(CT) -> moving(T1c), the inverse of the contract (T1c->CT);
|
||||
# fix_orientation resolves the direction.
|
||||
T_main = sitk.ReadTransform(itk_tf)
|
||||
T_init = sitk.TranslationTransform(3)
|
||||
T_init.SetOffset(init)
|
||||
cand = sitk.CompositeTransform([T_init, T_main])
|
||||
T, inverted = fix_orientation(ct, t1c, cand)
|
||||
T, inverted = fix_orientation(ct, t1c, T_main)
|
||||
finally:
|
||||
if not keep_scratch:
|
||||
shutil.rmtree(tmp, ignore_errors=True)
|
||||
|
|
@ -950,7 +883,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
"pipeline": "impact",
|
||||
"toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
|
||||
"metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)",
|
||||
"transform": (f"Rigid (EulerTransform, ASGD, {n} scales, centroid init; "
|
||||
"transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; "
|
||||
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
|
||||
"convention": "T1c->CT (verified)", "auto_inverted": inverted,
|
||||
"fixed": "T1c", "moving": "CT",
|
||||
|
|
@ -963,10 +896,128 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=500, samples=2000,
|
|||
return T, meta
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# 6. NiftyReg (KCL-BMEIS suite, reg_aladin)
|
||||
# --------------------------------------------------------------------------- #
|
||||
|
||||
|
||||
def _niftyreg_bin():
|
||||
"""Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable."""
|
||||
b = shutil.which("reg_aladin")
|
||||
if b:
|
||||
return b
|
||||
raise PipelineUnavailable(
|
||||
"niftyreg: reg_aladin not found on PATH. Install NiftyReg with "
|
||||
"`conda install -c conda-forge niftyreg` (NiftyReg 1.5.69, the "
|
||||
"KCL-BMEIS suite: reg_aladin rigid/affine + reg_f3d non-linear) and "
|
||||
"re-run in the activated env.")
|
||||
|
||||
|
||||
def _niftyreg_version(binpath):
|
||||
try:
|
||||
r = subprocess.run([binpath, "--version"], capture_output=True,
|
||||
text=True, timeout=60)
|
||||
m = re.search(r"(\d+\.\d+\.\d+)", r.stdout + r.stderr)
|
||||
if m:
|
||||
return f"NiftyReg {m.group(1)} (reg_aladin)"
|
||||
except Exception: # noqa
|
||||
pass
|
||||
return "NiftyReg (reg_aladin)"
|
||||
|
||||
|
||||
def _parse_niftyreg_affine(path):
|
||||
"""Parse reg_aladin's output affine (4x4, row-major) into a SimpleITK
|
||||
AffineTransform. The file is *documented* as
|
||||
``Affine * Reference = Floating``; the direction that actually satisfies
|
||||
the contract is verified empirically by :func:`fix_orientation` (on the
|
||||
M6 cases the gate keeps the file's affine as-is)."""
|
||||
M = np.array([float(x) for x in open(path).read().split()]).reshape(4, 4)
|
||||
T = sitk.AffineTransform(3)
|
||||
T.SetMatrix(M[:3, :3].ravel())
|
||||
T.SetTranslation(M[:3, 3])
|
||||
return T
|
||||
|
||||
|
||||
def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False):
|
||||
"""NiftyReg rigid CT->T1c (reference=CT, floating=T1c) via ``reg_aladin``.
|
||||
|
||||
reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global
|
||||
block-matching registration: at each level of a ``levels``-level Gaussian
|
||||
pyramid it matches 4-voxel blocks by normalised cross-correlation
|
||||
(NCC -- a cross-modality-tolerant block similarity, unlike the histogram
|
||||
MI objectives of the other pipelines), then solves a trimmed least squares
|
||||
(TLS) problem for the rigid transform, iterating ``maxit`` times per
|
||||
level. The symmetric variant (default) optimises forward and backward
|
||||
transforms concurrently; ``-maxit 5`` (the tool default) is too coarse for
|
||||
these repositioned pairs -- the benchmark uses 15.
|
||||
|
||||
Full volumes (no head box); default initialisation (image centres);
|
||||
``-rigOnly`` for the 6-DOF rigid contract; ``-interp 3`` (cubic) internal
|
||||
warping. The output affine is documented as ``Affine * Reference =
|
||||
Floating``; :func:`fix_orientation` resolves the direction that actually
|
||||
satisfies the contract against the skull ring (on the M6 cases it keeps
|
||||
the file's affine as-is, ``auto_inverted: False``) and returns the T1c ->
|
||||
CT transform.
|
||||
|
||||
Preconditions (else :class:`PipelineUnavailable`): ``reg_aladin`` on PATH
|
||||
(``conda install -c conda-forge niftyreg``).
|
||||
"""
|
||||
import time
|
||||
|
||||
binpath = _niftyreg_bin()
|
||||
# full volumes, reference = CT, floating = T1c (no head box)
|
||||
tmp = tempfile.mkdtemp(prefix="niftyreg_reg_")
|
||||
dt = 0.0
|
||||
try:
|
||||
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "ref.nii.gz"), True)
|
||||
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "flo.nii.gz"), True)
|
||||
cmd = [binpath, "-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
|
||||
"-rigOnly", "-aff", "aff.txt",
|
||||
"-ln", str(levels), "-maxit", str(maxit),
|
||||
"-interp", "3", "-omp", str(threads), "-voff"]
|
||||
t0 = time.time()
|
||||
proc = subprocess.run(cmd, capture_output=True, text=True,
|
||||
cwd=tmp, timeout=7200)
|
||||
dt = time.time() - t0
|
||||
if proc.returncode != 0:
|
||||
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
|
||||
raise RuntimeError(f"niftyreg: reg_aladin failed (rc={proc.returncode}):\n{tail}")
|
||||
aff = os.path.join(tmp, "aff.txt")
|
||||
if not os.path.exists(aff):
|
||||
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
|
||||
raise RuntimeError(f"niftyreg: no affine output produced:\n{tail}")
|
||||
T = _parse_niftyreg_affine(aff)
|
||||
# aff.txt is documented as CT -> T1c (Affine*Reference=Floating);
|
||||
# the contract is T1c -> CT. fix_orientation resolves the direction
|
||||
# empirically (on these cases it keeps the file's affine as-is).
|
||||
T, inverted = fix_orientation(ct, t1c, T)
|
||||
finally:
|
||||
if not keep_scratch:
|
||||
shutil.rmtree(tmp, ignore_errors=True)
|
||||
|
||||
meta = {
|
||||
"pipeline": "niftyreg",
|
||||
"toolkit": _niftyreg_version(binpath),
|
||||
"metric": "NCC block matching + trimmed least squares (symmetric, full volumes)",
|
||||
"transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, "
|
||||
f"{maxit} TLS iters/level, cubic interp)"),
|
||||
"convention": "T1c->CT (verified)", "auto_inverted": inverted,
|
||||
"fixed": "CT", "moving": "T1c",
|
||||
"seconds_reg": round(dt, 1),
|
||||
"levels": levels, "maxit": maxit,
|
||||
}
|
||||
if keep_scratch:
|
||||
meta["scratch"] = tmp
|
||||
return T, meta
|
||||
|
||||
|
||||
REGISTERS = {
|
||||
"antspyx": register_antspyx,
|
||||
"elastix": register_elastix,
|
||||
"impact": register_impact,
|
||||
"niftyreg": register_niftyreg,
|
||||
"simpleitk": register_simpleitk,
|
||||
"synthmorph": register_synthmorph,
|
||||
}
|
||||
|
|
@ -994,6 +1045,13 @@ 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 == "niftyreg":
|
||||
b = shutil.which("reg_aladin")
|
||||
if b:
|
||||
out[name] = (True, _niftyreg_version(b))
|
||||
else:
|
||||
out[name] = (False, "reg_aladin not on PATH "
|
||||
"(conda install -c conda-forge niftyreg)")
|
||||
elif name == "synthmorph":
|
||||
ok, note = _synthmorph_docker_ok()
|
||||
if not ok:
|
||||
|
|
|
|||
Loading…
Reference in a new issue