xfr/longitudinal
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add visual QA script for TV label warping and registration methods

Browse source 
- Implemented `qa_label_interpolation.py` to evaluate the registration of T1c and CT images using various methods (rigid, affine, and SyN).
- The script generates warped TV images, visualizations, and summary metrics for each registration method.
- Outputs include PNG figures for visual comparison, a summary JSON file with metrics, and cached registration transforms.
- Added functionality for volume preservation checks and detailed logging of processing steps.
This commit is contained in:
Furen Xiao 2026-10-04 08:32:56 +08:00
parent bec0e9834f
commit 89bb1e6331
5 changed files with 1071 additions and 91 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

231
REGISTRATION_BENCHMARK_AFFINE.md Normal file
View file

@ -0,0 +1,231 @@
# Registration Pipeline Benchmark — CT → T1c (M6) — Affine
Companion to [REGISTRATION_BENCHMARK_RIGID.md](REGISTRATION_BENCHMARK_RIGID.md)
and [REGISTRATION_BENCHMARK_DEFORMABLE.md](REGISTRATION_BENCHMARK_DEFORMABLE.md):
the **affine** (12-DOF) variant of the same benchmark. The rigid benchmark
asked "which toolkit aligns the head best with 6 DOF?"; this one asks "how
much does the full 12-DOF affine buy on top of rigid — and do the toolkits'
affine workflows still rank the same way?" — on the same 10 M6 CT→T1c test
cases, the same unified contract, and the same metrics.
Run date: 2026-10-03, 13:50–14:40 local (single run, 6 workers, 24 threads,
2 workers per GPU; ~50 min wall, dominated by `impact`'s ~6 min phase).
10 M6 test cases, 3× NVIDIA A100 80GB PCIe, 96-core host, conda env
`longitudinal-313`.
## Unified contract
Identical to the rigid benchmark: the returned transform `T` must satisfy
`sitk.Resample(ct, t1c, T, ...)` = CT on the T1c grid, i.e. **T maps T1c →
CT coordinates**. fixed = T1c, moving = CT, full volumes (no head box, no
mask, no crop, no pad), each toolkit's native automatic transform
initialization, CT clipped to [0, 80] HU only for the pipelines where the
rigid A/B showed clipping helps (niftyreg, synthmorph).
One change vs the rigid benchmark: every pipeline estimates a **full affine
(12-DOF)** transform instead of a rigid (6-DOF) one — each toolkit's native
affine workflow (or, for the antspyx rows, a stock linear-only preset).
## The nine pipelines (affine configurations)
| pipeline | affine configuration | metric |
|---|---|---|
| **antspyx** | antspyx stock `Affine` preset (the toolkit's native affine): Affine[0.25], 2100×1200×1200×10 iters, shrink 6×4×2×1, smoothing 3×2×1×0, center-of-mass auto init | Mattes MI (32 bins, 20 % regular sampling) |
| **ants_trsaa** | antspyx stock `TRSAA` preset: Translation[1] → Rigid[1] → Similarity[1] (each 2000×2000×0 iters, 4×2×1, 2×1×0) → Affine[1] × 2 (each 40×20×0, 4×2×1, 2×1×0), COM init — the full TRSA cascade | Mattes MI (32 bins), all stages |
| **ants_reproa** | `antsRegistrationSyNRepro[a]`: Rigid[0.1] + Affine[0.1] (each 1000×500×250×100, 8×4×2×1, 3×2×1×0vox), COM init — the linear-only half of the Repro workflow | GC (1/1/Regular/0.25), both stages |
| **ants_quickreproa** | `antsRegistrationSyNQuickRepro[a]`: Rigid[0.1] + Affine[0.1] (each 1000×500×250×0 — no full-res level —, 8×4×2×1, 3×2×1×0vox), COM init | GC (1/1/Regular/0.25), both stages |
| **elastix** | AffineTransform (12-DOF, the rigid benchmark's map with Euler→Affine), ASGD 2000 iters × 3 resolutions, CenterOfGravity auto init | Mattes MI (32 bins), 5000 samples/level |
| **impact** | the ImpactReg `MR_CT_TS` reference map with **AffineTransform** in place of EulerTransform (same MIND + TotalSegmentator static features, 4 scales 400/300/200/200 iters, 2000 samples, same CoG-init / 0.1 sample-gate adaptations) | IMPACT semantic (L1 + Dice) + Mattes MI |
| **niftyreg** | reg_aladin **full affine** (no `-rigOnly`), 3-level pyramid, 15 TLS iters/level | NCC block matching (init) + TLS |
| **simpleitk** | two-stage (the toolkit's canonical affine workflow): Euler3D stage (2000 iters/level, pyramid [4,2,1], identity start) → AffineTransform stage (same budget) initialized to reproduce the rigid stage exactly — a single-stage affine from identity collapses its scale on these pairs | Mattes MI (32 bins), both stages |
| **synthmorph** | `affine` matrix model (the toolkit's native 12-DOF deep model) | anatomy-aware contrast-invariant network (pretrained) |
The `ants_trsaa` / `ants_reproa` / `ants_quickreproa` rows are the linear-only
(`[a]`) antspyx presets — the same stage configurations the
`antsRegistrationSyNRepro*.sh` / `antsRegistrationSyNQuickRepro*.sh`
workflows use for their Rigid + Affine head.
## Test cases
Same 10 as the rigid / deformable benchmarks (first 10 M6 labelled cases of
`m6.discover()`, all with a GTV label):
`m6_22GD5VL2_20230724`, `m6_22IRRKU6_20240115`, `m6_23XIMO4M_20230222`,
`m6_24V62QQM_20230223`, `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820`,
`m6_26TYQI3C_20221027`, `m6_2BBN5XA4_20250411`, `m6_2C4TKI6H_20240506`,
`m6_2CDNN7BH_20240102`.
(`m6_22GD5VL2_20230724` has an incomplete CT head — a known data issue,
not a pipeline failure; it is scored as ct_incomplete in every pipeline.)
## Metrics (identical for every pipeline)
- **Mattes MI** (nats) between the warped CT and the T1c over the MR head
(largest connected component of the T1c mask) — achieved alignment, higher
= better.
- **IMPACT loss** (vboussot/ImpactLoss; pretrained TotalSegmentator
`TS/M730` deep features — weighted L1 over the two deepest layers, full
volumes; the same scorer as the other benchmarks) between the warped CT
and the T1c — semantic-anatomy similarity, lower = better.
- **Registration wall time** (seconds, excludes I/O + scoring).
Per-case JSONs additionally retain the production skull-ring QA (gap +
concentricity, and the acceptance verdict) for reference.
### Result
Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ T1c
over the MR head; higher = better) **and** the **ImpactLoss** semantic score
(warped-CT ~ T1c over the full volumes; lower = better), on the same
unclipped scoring volumes as the other benchmarks, so all rows are
comparable on the same scales and against the rigid / deformable numbers.
`niftyreg` and `synthmorph` clip the CT (the moving image) to [0, 80] HU
*before registration* (per the rigid benchmark's clipping A/B); the others
register the CT as-is. All nine transforms are affine (12-DOF).
All nine pipelines registered 9/9 cases (the 10th is the known
incomplete-CT data issue, unscored in every pipeline). **Bold** = best for
the column (MI higher is better, ImpactLoss lower is better).
| pipeline | toolkit | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) |
|---|---|---|---|---|---|
| **antspyx** | ANTsPy 0.6.3 (`Affine` preset) | 9 | 12.6 | 0.156 / 0.155 | **0.660 / 0.660** |
| **ants_trsaa** | ANTsPy 0.6.3 (`TRSAA` preset) | 9 | 30.4 | 0.156 / 0.155 | 0.658 / 0.659 |
| **ants_reproa** | ANTsPy 0.6.3 (`Repro[a]`) | 9 | 63.1 | 0.160 / 0.159 | 0.525 / 0.514 |
| **ants_quickreproa** | ANTsPy 0.6.3 (`QuickRepro[a]`) | 9 | 28.4 | 0.161 / 0.159 | 0.523 / 0.514 |
| **elastix** | ITK-Elastix 0.25.4 | 9 | 12.3 | 0.156 / 0.156 | 0.651 / 0.645 |
| **impact** | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT, affine MIND+TS reference map) | 9 | 215.5 | 0.159 / 0.158 | 0.623 / 0.567 |
| **niftyreg** | NiftyReg 2.1.1 (reg_aladin full affine, Docker + CUDA) | 9 | 19.8 | 0.167 / 0.170 | 0.278 / 0.330 |
| **simpleitk** | SimpleITK 2.5.6 | 9 | 23.3 | 0.156 / 0.156 | 0.654 / 0.641 |
| **synthmorph** | Docker `freesurfer/synthmorph` (TensorFlow, pretrained, `affine` model) | 9 | 27.1 | 0.167 / 0.171 | 0.282 / 0.320 |
Per-case MI (nats, unclipped score; **bold** = best for the case):
| case | antspyx | ants_trsaa | ants_reproa | ants_quickreproa | elastix | impact | niftyreg | simpleitk | synthmorph |
|---|---|---|---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | **0.655** | 0.647 | 0.547 | 0.546 | 0.638 | 0.391 | 0.278 | 0.654 | 0.282 |
| m6_23XIMO4M_20230222 | **0.661** | 0.660 | 0.519 | 0.521 | 0.656 | 0.657 | 0.611 | 0.660 | 0.520 |
| m6_24V62QQM_20230223 | **0.710** | 0.708 | 0.568 | 0.569 | 0.684 | 0.708 | 0.254 | **0.710** | 0.252 |
| m6_24V62QQM_20231017 | 0.628 | **0.635** | 0.425 | 0.423 | 0.622 | 0.623 | 0.473 | 0.619 | 0.479 |
| m6_24V62QQM_20240820 | 0.670 | **0.673** | 0.459 | 0.464 | 0.670 | 0.659 | 0.121 | 0.545 | 0.122 |
| m6_26TYQI3C_20221027 | **0.719** | 0.717 | 0.525 | 0.523 | 0.675 | 0.676 | 0.593 | 0.699 | 0.570 |
| m6_2BBN5XA4_20250411 | **0.590** | 0.588 | 0.505 | 0.507 | 0.568 | 0.360 | 0.069 | **0.590** | 0.071 |
| m6_2C4TKI6H_20240506 | **0.646** | 0.645 | 0.546 | 0.542 | 0.636 | 0.451 | 0.241 | 0.635 | 0.256 |
| m6_2CDNN7BH_20240102 | **0.660** | 0.658 | 0.535 | 0.534 | 0.651 | 0.575 | 0.327 | 0.656 | 0.328 |
Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest)
for the case):
| case | antspyx | ants_trsaa | ants_reproa | ants_quickreproa | elastix | impact | niftyreg | simpleitk | synthmorph |
|---|---|---|---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 0.1514 | 0.1516 | 0.1527 | 0.1528 | **0.1509** | 0.1607 | 0.1667 | 0.1513 | 0.1666 |
| m6_23XIMO4M_20230222 | 0.1565 | 0.1564 | 0.1598 | 0.1599 | 0.1565 | 0.1564 | 0.1585 | **0.1563** | 0.1616 |
| m6_24V62QQM_20230223 | **0.1466** | 0.1468 | 0.1527 | 0.1529 | 0.1470 | 0.1468 | 0.1670 | **0.1466** | 0.1672 |
| m6_24V62QQM_20231017 | 0.1519 | **0.1517** | 0.1599 | 0.1599 | 0.1527 | 0.1533 | 0.1587 | 0.1520 | 0.1581 |
| m6_24V62QQM_20240820 | 0.1560 | **0.1559** | 0.1618 | 0.1617 | 0.1564 | 0.1570 | 0.1840 | 0.1599 | 0.1843 |
| m6_26TYQI3C_20221027 | **0.1563** | 0.1570 | 0.1604 | 0.1615 | 0.1615 | 0.1657 | 0.1621 | 0.1585 | 0.1624 |
| m6_2BBN5XA4_20250411 | 0.1617 | 0.1617 | 0.1622 | 0.1622 | **0.1593** | 0.1662 | 0.1845 | 0.1617 | 0.1849 |
| m6_2C4TKI6H_20240506 | **0.1559** | 0.1561 | 0.1612 | 0.1613 | 0.1561 | 0.1594 | 0.1759 | **0.1559** | 0.1766 |
| m6_2CDNN7BH_20240102 | 0.1591 | 0.1590 | 0.1609 | 0.1617 | **0.1587** | **0.1587** | 0.1740 | **0.1587** | 0.1741 |
### Interpretation
- **The four intensity tools form a statistical tie, and affine barely
moves them.** `antspyx` (0.660), `ants_trsaa` (0.658), `simpleitk`
(0.654) and `elastix` (0.651) all sit within 0.009 nats of each other on
the median (ImpactLoss medians 0.156), and per-case MI bests split
antspyx 5/9, trsaa 2/9, simpleitk 2/9 (ties). Against the rigid
benchmark's per-pipeline medians, affine buys only
`antspyx` 0.652 → 0.660 (+0.008), `simpleitk` 0.647 → 0.654 (+0.007),
`elastix` 0.646 → 0.651 (+0.005) — the M6 head is essentially rigid
between sessions, so the extra 6 DOF mostly absorb inter-session
anisotropic scaling (≤ a few %), not misalignment.
- **`ants_trsaa` is the full TRSA cascade for free quality — and it is not
free.** Translation → Rigid → Similarity → Affine → Affine (5 linear
stages) reaches 0.658, i.e. the same as the single-stage `Affine` preset
(0.660) at 2.4× the wall time (30.4 s vs 12.6 s): the cascaded 6-DOF and
9-DOF stages just re-solve what the one 12-DOF stage finds.
- **`ants_quickreproa` ≈ `ants_reproa` quality at 45 % of the time**
(0.523 vs 0.525 MI; 28.4 s vs 63.1 s): dropping the full-resolution rigid
level costs nothing. But both are clearly the weakest of the *top-tier
eligible* antspyx presets — the GC-based Rigid+Affine head (the linear
half of the Repro T1w-T1w workflow) locks a cross-modality-weak pose, the
same story as the Repro family in the deformable benchmark (0.530–0.549
MI there).
- **`impact`'s affine run degrades vs its own rigid.** Median 0.623
(mean 0.567) vs 0.652 for rigid: the semantic 12-DOF landscape has more
local minima, and the single CenterOfGravity start falls into two of them
(`m6_22IRRKU6` 0.391, `m6_2BBN5XA4` 0.360 — both 0.645 / 0.572 in the
rigid run). It also roughly doubles in wall time (215.5 s vs 111.2 s)
for the extra 6 DOF.
- **`niftyreg` full affine does not fix its gross-pose failures** (0.278 vs
0.292 rigid; the two repositioned cases still at 0.069 / 0.121): reg_
aladin's block-matching + TLS is a local optimizer either way, and the
output affine is still flipped by the orientation gate on 8/9 cases.
- **`synthmorph`'s `affine` model confirms the rigid-benchmark
investigation** (0.282 vs 0.277; within ±0.03 per case): the linear part
of its output is a near-pure rotation — the extra 6 DOF are essentially
unused, and the pretrained matrix head's gross-pose estimation weakness
stands for 12 DOF as for 6.
- **The deformable stage is still worth it.** The best affine median here
is 0.660 (`antspyx`); the best deformable in the companion benchmark is
0.753 (`ants_rigid_synonly`; 0.723 for the stock `SyN` preset) — the
deformable stage adds +0.06 to +0.09 MI on top of the best affine, which
is fine-tissue alignment, not gross pose.
- **Production skull-ring QA (reference only, kept in the per-case JSONs):**
acceptance 4/9 `antspyx` / `ants_reproa` / `ants_quickreproa` /
`elastix` / `simpleitk`, 3/9 `ants_trsaa`, 2/9 `impact`, 1/9 `niftyreg` /
`synthmorph`. Concentricity: `ants_reproa` tightest mid-gap (5.6 mm),
`niftyreg` / `impact` / `synthmorph` worst on the rotational term
(23.5–25.3 mm std) — the same weak correlation with perceived alignment
noted in the other benchmarks, so it is reported for reference only.
## Transform conventions (affine-specific)
- **antspyx linear-only presets write one merged affine per consecutive
linear run.** The `TRSAA` cascade (5 linear stages) comes back as a single
*binary* `0GenericAffine.mat` that SimpleITK cannot parse; the pipeline
recovers it by point-probing the transform (exact for affine maps,
verified to NCC 1.0 against ANTs' own `warpedmovout` on the pilot cases).
`Repro[a]` / `QuickRepro[a]` (Rigid + Affine) merge the same way.
- **`niftyreg`'s `-aff` file** is documented as `Affine * Reference =
Floating`; the `fix_orientation` gate resolves the direction empirically
and inverts it on 8/9 M6 cases (as in the rigid benchmark).
- **`elastix` / `impact`** return the full combination transform (auto-init
folded in) in the contract direction (`auto_inverted: False`); the impact
affine read-back is `TransformParameters.0.itk.txt` (parsed by
SimpleITK, with a manual 12-parameter fallback).
- **`simpleitk`** returns the two-stage composite (rigid-init affine o
optimized affine) condensed to a single AffineTransform by point mapping.
- **`synthmorph`**'s `-T` LTA (LINEAR_RAS_TO_RAS 4×4) is the contract
transform directly, direction-checked by the usual gate.
## Reproducing
```bash
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
cd /mnt/b4/xfr/git26/longitudinal
# all nine affine pipelines (results go to results/bench_registration_affine/;
# 2026-10-03 run, 6 workers, ~50 min wall dominated by impact's ~6 min phase)
python scripts/bench_registration.py --affine --pipelines "antspyx ants_trsaa ants_reproa ants_quickreproa elastix impact niftyreg simpleitk synthmorph" --workers 6 --threads 24 --qa-images
# single pipeline / case subset
python scripts/bench_registration.py --affine --pipelines antspyx simpleitk --n-cases 3
# combine the per-pipeline JSONs into summary.json / summary.md / plot
python scripts/bench_registration.py --affine --merge
```
## Visual QA images
`--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to
`results/bench_registration_affine/qa/<pipeline>__<case>.png` (9 per
pipeline; the incomplete-CT case is skipped) — same format as the other
benchmarks: a 3×2 grid of the mid-axial / mid-coronal / mid-sagittal slices
through the head center, **left column** checkerboard blend of the warped CT
vs the T1c (clean registration = seamless mosaic), **right column** additive
color overlay (CT red, T1c green; agreement = yellow, misregistration =
red/green fringing). CT is brain-windowed 0–80 HU, T1c 2–98 % head window;
the title carries the pipeline, case key and the achieved Mattes MI.

100
REGISTRATION_BENCHMARK_DEFORMABLE.md
View file

@ -245,10 +245,14 @@ for the case):
## ANTs variants: direct SyNRA / SyN and the Repro / QuickRepro workflows ## ANTs variants: direct SyNRA / SyN and the Repro / QuickRepro workflows
The current `antspyx` row above is antspyx's stock **`SyN` preset** (Affine The `ants_syn` row below is antspyx's stock **`SyN` preset** (Affine
2100×1200×1200×0 + SyN[0.2,3,0] 40×20×0, Mattes MI 32, center-of-mass init). 2100×1200×1200×0 + SyN[0.2,3,0] 40×20×0, Mattes MI 32, center-of-mass
As a follow-up (2026-10-03, 6 workers as the main run) seven further init) — the same configuration as the `antspyx` row of the main table
ANTs configurations were benchmarked on the same 10 cases and the same above, re-run in this follow-up batch and reproduced case-for-case within
0.002 nats (results in `ants_syn.json`), listed under its own name to
avoid confusing it with the `ants_synra` / `ants_rigid_synonly` variants.
As a follow-up (2026-10-03, 6 workers as the main run) six further ANTs
configurations were benchmarked on the same 10 cases and the same
contract — the question: is there a better stock ANTs deformable workflow contract — the question: is there a better stock ANTs deformable workflow
for CT→T1c than the `SyN` preset? for CT→T1c than the `SyN` preset?
@ -256,7 +260,6 @@ for CT→T1c than the `SyN` preset?
|---|---|---| |---|---|---|
| **ants_rigid_synonly** | legacy two-call (`Rigid` + `SyNOnly`) | Rigid[0.25] (Mattes 32, 500×500×200, 4×2×1, 3×2×0) → SyNOnly SyN[0.2,3,0] (Mattes 32, 100×100×0, 4×2×1, 2×1×0) initialized from the Rigid stage; `Composite.h5` read-back | | **ants_rigid_synonly** | legacy two-call (`Rigid` + `SyNOnly`) | Rigid[0.25] (Mattes 32, 500×500×200, 4×2×1, 3×2×0) → SyNOnly SyN[0.2,3,0] (Mattes 32, 100×100×0, 4×2×1, 2×1×0) initialized from the Rigid stage; `Composite.h5` read-back |
| **ants_synra** | `SyNRA` | Rigid[0.25] + Affine[0.25] (each Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1) | | **ants_synra** | `SyNRA` | Rigid[0.25] + Affine[0.25] (each Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1) |
| **ants_syn** | `SyN` | Affine[0.25] (Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1); COM init — the stock preset, identical configuration to the existing `antspyx` row |
| **ants_synreprosr** | `antsRegistrationSyNRepro[sr]` | Rigid[0.1] (GC 1/1/Regular/0.25, 1000×500×250×100, 8×4×2×1, 3×2×1×0vox) + SyN[0.2,3,0] (CC 1/2, 50×20, 2×1, 1×0vox) | | **ants_synreprosr** | `antsRegistrationSyNRepro[sr]` | Rigid[0.1] (GC 1/1/Regular/0.25, 1000×500×250×100, 8×4×2×1, 3×2×1×0vox) + SyN[0.2,3,0] (CC 1/2, 50×20, 2×1, 1×0vox) |
| **ants_synreprobr** | `antsRegistrationSyNRepro[br]` | the `Repro` Rigid[0.1] stage + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×20, 2×1, 1×0vox) | | **ants_synreprobr** | `antsRegistrationSyNRepro[br]` | the `Repro` Rigid[0.1] stage + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×20, 2×1, 1×0vox) |
| **ants_synquickreprosr** | `antsRegistrationSyNQuickRepro[sr]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + SyN[0.2,3,0] (CC 1/2, 50×0, 2×1) | | **ants_synquickreprosr** | `antsRegistrationSyNQuickRepro[sr]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + SyN[0.2,3,0] (CC 1/2, 50×0, 2×1) |
@ -280,16 +283,14 @@ linear stages (`SyNRA`), it writes the merged linear stage as a *binary*
by point-probing the transform (exact for affine maps, verified to by point-probing the transform (exact for affine maps, verified to
2×10⁻⁴ mm). 2×10⁻⁴ mm).
All six variants registered 9/9 (the 10th is the known incomplete-CT data All seven rows registered 9/9 (the 10th is the known incomplete-CT data
issue). **Bold** = best for the column; `antspyx` is the existing row, issue). **Bold** = best for the column.
re-listed for reference.
| pipeline | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) | | pipeline | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) |
|---|---|---|---|---| |---|---|---|---|---|
| **ants_rigid_synonly** | 9 | 208.4 | 0.156 / 0.155 | **0.753 / 0.763** | | **ants_rigid_synonly** | 9 | 208.4 | 0.156 / 0.155 | **0.753 / 0.763** |
| **antspyx** (existing) | 9 | 31.7 | 0.158 / 0.156 | 0.723 / 0.729 |
| **ants_synra** | 9 | 43.6 | 0.158 / 0.156 | 0.722 / 0.729 |
| **ants_syn** | 9 | 33.8 | 0.158 / 0.156 | 0.722 / 0.729 | | **ants_syn** | 9 | 33.8 | 0.158 / 0.156 | 0.722 / 0.729 |
| **ants_synra** | 9 | 43.6 | 0.158 / 0.156 | 0.722 / 0.729 |
| **ants_synreprosr** | 9 | 294.7 | **0.151 / 0.153** | 0.530 / 0.529 | | **ants_synreprosr** | 9 | 294.7 | **0.151 / 0.153** | 0.530 / 0.529 |
| **ants_synreprobr** | 9 | 370.3 | 0.156 / 0.157 | 0.541 / 0.537 | | **ants_synreprobr** | 9 | 370.3 | 0.156 / 0.157 | 0.541 / 0.537 |
| **ants_synquickreprosr** | 9 | 83.7 | 0.154 / 0.156 | 0.542 / 0.534 | | **ants_synquickreprosr** | 9 | 83.7 | 0.154 / 0.156 | 0.542 / 0.534 |
@ -297,61 +298,56 @@ re-listed for reference.
Per-case MI (nats, unclipped score; **bold** = best for the case): Per-case MI (nats, unclipped score; **bold** = best for the case):
| case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | | case | ants_syn | ants_rigid_synonly | ants_synra | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr |
|---|---|---|---|---|---|---|---|---| |---|---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 0.679 | **0.741** | 0.681 | 0.681 | 0.623 | 0.646 | 0.630 | 0.627 | | m6_22IRRKU6_20240115 | 0.681 | **0.741** | 0.681 | 0.623 | 0.646 | 0.630 | 0.627 |
| m6_23XIMO4M_20230222 | 0.723 | **0.760** | 0.722 | 0.722 | 0.530 | 0.544 | 0.529 | 0.530 | | m6_23XIMO4M_20230222 | 0.722 | **0.760** | 0.722 | 0.530 | 0.544 | 0.529 | 0.530 |
| m6_24V62QQM_20230223 | 0.789 | **0.824** | 0.789 | 0.787 | 0.520 | 0.528 | 0.542 | 0.549 | | m6_24V62QQM_20230223 | 0.787 | **0.824** | 0.789 | 0.520 | 0.528 | 0.542 | 0.549 |
| m6_24V62QQM_20231017 | **0.720** | 0.716 | **0.720** | **0.720** | 0.473 | 0.475 | 0.469 | 0.462 | | m6_24V62QQM_20231017 | **0.720** | 0.716 | **0.720** | 0.473 | 0.475 | 0.469 | 0.462 |
| m6_24V62QQM_20240820 | 0.738 | **0.753** | 0.738 | 0.738 | 0.524 | 0.525 | 0.517 | 0.524 | | m6_24V62QQM_20240820 | 0.738 | **0.753** | 0.738 | 0.524 | 0.525 | 0.517 | 0.524 |
| m6_26TYQI3C_20221027 | 0.799 | **0.857** | 0.798 | 0.799 | 0.557 | 0.541 | 0.599 | 0.632 | | m6_26TYQI3C_20221027 | 0.799 | **0.857** | 0.798 | 0.557 | 0.541 | 0.599 | 0.632 |
| m6_2BBN5XA4_20250411 | 0.665 | **0.696** | 0.665 | 0.666 | 0.230 | 0.303 | 0.235 | 0.297 | | m6_2BBN5XA4_20250411 | 0.666 | **0.696** | 0.665 | 0.230 | 0.303 | 0.235 | 0.297 |
| m6_2C4TKI6H_20240506 | 0.719 | **0.751** | 0.719 | 0.719 | 0.672 | 0.649 | 0.654 | 0.621 | | m6_2C4TKI6H_20240506 | 0.719 | **0.751** | 0.719 | 0.672 | 0.649 | 0.654 | 0.621 |
| m6_2CDNN7BH_20240102 | 0.726 | **0.765** | 0.726 | 0.725 | 0.629 | 0.622 | 0.631 | 0.599 | | m6_2CDNN7BH_20240102 | 0.725 | **0.765** | 0.726 | 0.629 | 0.622 | 0.631 | 0.599 |
Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) for Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) for
the case): the case):
| case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | | case | ants_syn | ants_rigid_synonly | ants_synra | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr |
|---|---|---|---|---|---|---|---|---| |---|---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | 0.1541 | 0.1530 | 0.1542 | 0.1540 | **0.1442** | 0.1482 | 0.1474 | 0.1493 | | m6_22IRRKU6_20240115 | 0.1540 | 0.1530 | 0.1542 | **0.1442** | 0.1482 | 0.1474 | 0.1493 |
| m6_23XIMO4M_20230222 | 0.1581 | 0.1565 | 0.1581 | 0.1581 | **0.1505** | 0.1549 | 0.1535 | 0.1585 | | m6_23XIMO4M_20230222 | 0.1581 | 0.1565 | 0.1581 | **0.1505** | 0.1549 | 0.1535 | 0.1585 |
| m6_24V62QQM_20230223 | 0.1473 | **0.1464** | 0.1472 | 0.1472 | 0.1515 | 0.1553 | 0.1535 | 0.1560 | | m6_24V62QQM_20230223 | 0.1472 | **0.1464** | 0.1472 | 0.1515 | 0.1553 | 0.1535 | 0.1560 |
| m6_24V62QQM_20231017 | 0.1518 | **0.1512** | 0.1516 | 0.1519 | 0.1559 | 0.1575 | 0.1577 | 0.1581 | | m6_24V62QQM_20231017 | 0.1519 | **0.1512** | 0.1516 | 0.1559 | 0.1575 | 0.1577 | 0.1581 |
| m6_24V62QQM_20240820 | **0.1545** | **0.1545** | 0.1547 | 0.1546 | 0.1553 | 0.1578 | 0.1586 | 0.1586 | | m6_24V62QQM_20240820 | 0.1546 | **0.1545** | 0.1547 | 0.1553 | 0.1578 | 0.1586 | 0.1586 |
| m6_26TYQI3C_20221027 | 0.1580 | 0.1562 | 0.1578 | 0.1579 | **0.1561** | 0.1574 | 0.1602 | 0.1610 | | m6_26TYQI3C_20221027 | 0.1579 | 0.1562 | 0.1578 | **0.1561** | 0.1574 | 0.1602 | 0.1610 |
| m6_2BBN5XA4_20250411 | 0.1624 | **0.1611** | 0.1624 | 0.1624 | 0.1680 | 0.1695 | 0.1690 | 0.1707 | | m6_2BBN5XA4_20250411 | 0.1624 | **0.1611** | 0.1624 | 0.1680 | 0.1695 | 0.1690 | 0.1707 |
| m6_2C4TKI6H_20240506 | 0.1594 | 0.1582 | 0.1594 | 0.1593 | **0.1490** | 0.1526 | 0.1517 | 0.1545 | | m6_2C4TKI6H_20240506 | 0.1593 | 0.1582 | 0.1594 | **0.1490** | 0.1526 | 0.1517 | 0.1545 |
| m6_2CDNN7BH_20240102 | 0.1590 | 0.1582 | 0.1589 | 0.1589 | **0.1493** | 0.1560 | 0.1538 | 0.1580 | | m6_2CDNN7BH_20240102 | 0.1589 | 0.1582 | 0.1589 | **0.1493** | 0.1560 | 0.1538 | 0.1580 |
### Interpretation (ANTs variants) ### Interpretation (ANTs variants)
- **`ants_rigid_synonly` is the best ANTs configuration of the eight** — - **`ants_rigid_synonly` is the best ANTs configuration of the seven** —
the only one that beats the current `antspyx` row: median MI **0.753** the only one that beats the current `ants_syn` row: median MI **0.753**
(per-case best on 8 of 9; its 0.857 on `m6_26TYQI3C` is the highest MI (per-case best on 8 of 9; its 0.857 on `m6_26TYQI3C` is the highest MI
of the whole benchmark, ahead of the previous 0.801) at ImpactLoss of the whole benchmark, ahead of the previous 0.801) at ImpactLoss
0.156 (second-best median after `synthmorph`'s 0.156). The gain over 0.156 (second-best median after `synthmorph`'s 0.156). The gain over
the `SyN` preset comes from the *deeper* SyN stage: 100 iterations at the `SyN` preset comes from the *deeper* SyN stage: 100 iterations at
every level **including full resolution** (100×100×0 vs 40×20×0) every level **including full resolution** (100×100×0 vs 40×20×0)
refining the Rigid stage's pose. The cost is 208.4 s/case (~6.6× the refining the Rigid stage's pose. The cost is 208.4 s/case (~6.2× the
preset) — the full-resolution SyN iterations dominate. `ants_syn` time) — the full-resolution SyN iterations dominate.
- **`ants_synra` is indistinguishable from the current `antspyx` row** on - **`ants_synra` is indistinguishable from the `ants_syn` row** on
quality (median MI 0.722 vs 0.723, ImpactLoss 0.158 vs 0.158; per-case quality (median MI 0.722 vs 0.722, ImpactLoss 0.158 vs 0.158; per-case
deltas within 0.002 nats, ties on 4 of 9) at **+38 % wall time** deltas within 0.002 nats, identical on 6 of 9) at **+29 % wall time**
(43.6 s vs 31.7 s): the explicit Rigid stage ahead of the Affine stage (43.6 s vs 33.8 s): the explicit Rigid stage ahead of the Affine stage
buys nothing — on these pairs the 12-DOF stage from the center-of-mass buys nothing — on these pairs the 12-DOF stage from the center-of-mass
init already finds the gross pose, so the extra 6-DOF stage just init already finds the gross pose, so the extra 6-DOF stage just
re-solves it. re-solves it.
- **`ants_syn` reproduces the existing `antspyx` row** (median MI 0.722
vs 0.723, ImpactLoss 0.158, s/reg 33.8 vs 31.7; per-case deltas within
0.002 nats) — as expected, since it is the *same stock `SyN` preset*
(Affine + SyN[0.2,3,0], COM init) re-run in this follow-up batch: it
serves as a consistency check that the shared read-back / scoring path
is configuration-stable.
- **The Repro / QuickRepro family is semantically competitive but weaker on - **The Repro / QuickRepro family is semantically competitive but weaker on
intensity:** median MI 0.530–0.549 (vs 0.723 for the `SyN` preset), intensity:** median MI 0.530–0.549 (vs 0.722 for the `SyN` preset, i.e.
`ants_syn`),
driven by the GC-based Rigid stage and the short 2-level SyN/BSplineSyN driven by the GC-based Rigid stage and the short 2-level SyN/BSplineSyN
refinement (50×20 or 50×0 iters, 2×1 pyramid, CC 2-bin metric) — a refinement (50×20 or 50×0 iters, 2×1 pyramid, CC 2-bin metric) — a
workflow tuned for intra-modality T1w-T1w. `ants_synreprosr` still takes workflow tuned for intra-modality T1w-T1w. `ants_synreprosr` still takes
@ -363,10 +359,10 @@ the case):
slowest ANTs configuration (370 s/case) at the same (low) MI as the slowest ANTs configuration (370 s/case) at the same (low) MI as the
84–96 s `QuickRepro` variants (which drop the 4th rigid level and the 84–96 s `QuickRepro` variants (which drop the 4th rigid level and the
2nd SyN level, changing the MI by < 0.02). Only `ants_rigid_synonly` 2nd SyN level, changing the MI by < 0.02). Only `ants_rigid_synonly`
challenges the `SyN` preset's quality (better MI at 6.6× the time); challenges `ants_syn`'s quality (better MI at ≈6× the time); every
every other variant is worse for more or for the same. other variant is worse for more or for the same.
- **Production skull-ring QA (reference only):** acceptance 4/9 - **Production skull-ring QA (reference only):** acceptance 4/9
`ants_syn` / `ants_synra` / `ants_rigid_synonly` (= baseline), 2/9 `ants_syn` / `ants_synra` / `ants_rigid_synonly`, 2/9
`ants_synreprobr`, 1/9 the rest; rotational dispersion worst in the `ants_synreprobr`, 1/9 the rest; rotational dispersion worst in the
`Repro` / `QuickRepro` family (24–29 mm std) vs 14–15 mm for `Repro` / `QuickRepro` family (24–29 mm std) vs 14–15 mm for
`ants_syn` / `ants_synra` / `ants_rigid_synonly` — consistent with `ants_syn` / `ants_synra` / `ants_rigid_synonly` — consistent with

55
scripts/bench_registration.py
View file

@ -293,6 +293,28 @@ def make_qa_png(warped_ct, mr, out_png, title, block=24, ct_brain_window=(0.0, 8
def _call(pipeline, ct, t1c, tv, args): def _call(pipeline, ct, t1c, tv, args):
if args.affine:
if pipeline in ("antspyx", "ants_trsaa", "ants_reproa", "ants_quickreproa"):
return RP.register_ants_affine_variant(ct, t1c, preset_name=pipeline)
if pipeline == "elastix":
return RP.register_elastix_affine(ct, t1c, workdir=args.elastix_dir,
threads=args.threads)
if pipeline == "impact":
return RP.register_impact_affine(ct, t1c, gpu=args.impact_gpu,
iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS,
samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES,
threads=args.threads)
if pipeline == "niftyreg":
return RP.register_niftyreg_affine(ct, t1c, threads=args.threads,
backend=args.niftyreg_backend,
gpu=args.niftyreg_gpu)
if pipeline == "simpleitk":
return RP.register_simpleitk_affine(ct, t1c, threads=args.threads)
if pipeline == "synthmorph":
return RP.register_synthmorph_affine(ct, t1c, gpu=args.synthmorph_gpu,
threads=args.threads,
extent=args.synthmorph_extent)
raise KeyError(pipeline)
if args.deformable: if args.deformable:
if pipeline == "simpleitk": if pipeline == "simpleitk":
return RP.register_simpleitk_deformable(ct, t1c, threads=args.threads) return RP.register_simpleitk_deformable(ct, t1c, threads=args.threads)
@ -526,9 +548,24 @@ def merge(args):
raise SystemExit(f"no pipeline results found in {out_dir}") raise SystemExit(f"no pipeline results found in {out_dir}")
print(f"[merge] {len(rows)} pipeline result files in {out_dir}") print(f"[merge] {len(rows)} pipeline result files in {out_dir}")
print_table(rows) print_table(rows)
deformable = (os.path.basename(os.path.normpath(args.out_dir)) out_base = os.path.basename(os.path.normpath(args.out_dir))
== "bench_registration_deformable") affine = out_base == "bench_registration_affine"
if deformable: deformable = out_base == "bench_registration_deformable"
if affine:
protocol = ("affine CT->T1c; cross-modality; each toolkit estimates a "
"full affine (12-DOF) transform with its native workflow "
"(antspyx: stock Affine preset + TRSAA / Repro[a] / "
"QuickRepro[a] linear-only presets; elastix AffineTransform; "
"impact MIND+TS reference map with AffineTransform; "
"niftyreg reg_aladin full affine; simpleitk AffineTransform; "
"synthmorph affine matrix model); alignment scored by "
"achieved Mattes MI (nats, warped-CT ~ MR over the MR head, "
"higher = better) and the IMPACT semantic loss "
"(vboussot/ImpactLoss, TS/M730 features, warped-CT ~ MR over "
"the full volumes, lower = better); reg time excludes "
"I/O+scoring; per-case JSON retains the production "
"skull-ring gap / acceptance QA")
elif deformable:
protocol = ("deformable CT->T1c; cross-modality; each pipeline runs its " protocol = ("deformable CT->T1c; cross-modality; each pipeline runs its "
"toolkit's native deformable workflow (linear init + deformable " "toolkit's native deformable workflow (linear init + deformable "
"stage, or a single deformable stage where that is native); the " "stage, or a single deformable stage where that is native); the "
@ -559,7 +596,8 @@ def merge(args):
json.dump(summary, f, indent=1) json.dump(summary, f, indent=1)
print("summary:", sf) print("summary:", sf)
write_plot(rows, os.path.join(out_dir, "bench_registration.png"), write_plot(rows, os.path.join(out_dir, "bench_registration.png"),
label="deformable" if deformable else "rigid") label="deformable" if deformable
else ("affine" if affine else "rigid"))
# markdown table for the report # markdown table for the report
md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |", md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |",
"|---|---|---|---|---|---|---|---|---|"] "|---|---|---|---|---|---|---|---|---|"]
@ -582,6 +620,12 @@ def merge(args):
def parse_args(): def parse_args():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--affine", action="store_true",
help="run the affine (12-DOF) variants of the pipelines "
"(each toolkit's affine workflow + the antspyx "
"Affine / TRSAA / Repro[a] / QuickRepro[a] presets; "
"results go to results/bench_registration_affine "
"by default)")
ap.add_argument("--deformable", action="store_true", ap.add_argument("--deformable", action="store_true",
help="run the deformable variants of the pipelines " help="run the deformable variants of the pipelines "
"(each toolkit's native deformable workflow; results go " "(each toolkit's native deformable workflow; results go "
@ -650,6 +694,9 @@ def main():
if args.deformable: if args.deformable:
if args.out_dir == os.path.join(ROOT, "results", "bench_registration"): if args.out_dir == os.path.join(ROOT, "results", "bench_registration"):
args.out_dir = os.path.join(ROOT, "results", "bench_registration_deformable") args.out_dir = os.path.join(ROOT, "results", "bench_registration_deformable")
if args.affine:
if args.out_dir == os.path.join(ROOT, "results", "bench_registration"):
args.out_dir = os.path.join(ROOT, "results", "bench_registration_affine")
if args.merge: if args.merge:
merge(args) merge(args)
return return

436
scripts/qa_label_interpolation.py Normal file
View file

@ -0,0 +1,436 @@
"""Visual QA of TV label warping: registration method x label interpolation method.
For one M6 case (default: 2BBN5XA4 / 20250411):
- fixed = T1c (MR, picked by scripts/02_build_m6_dataset.pick_t1c),
moving = CT (RT/ct_image.nii.gz); the TV (RT/TV-*.nii.gz) lives in the CT grid.
- registration via ANTsPy presets (src/registration_pipelines.py), each returning
a transform T that maps T1c -> CT, so sitk.Resample(tv, t1c, T, ...) places
the TV in MR (T1c) coordinates:
ants_rigid : Rigid[0.25] (Mattes, 500x500x200, 4x2x1)
ants_affine : Affine[0.25] (Mattes, 2100x1200x1200x10, 6x4x2x1)
ants_syn : SyN (Affine[0.25] + SyN[0.2,3,0] stages)
- the binary TV is then warped to the T1c grid by nearest-neighbor (nn) and
by SignedMaurerDistanceMap (mm, inside<0) -> {linear, B-spline order-3}
interp -> threshold (smd <= thr, default 0.5)
- one PNG per registration: 3 rows (axial / coronal / sagittal) x 4 columns
(leftmost = bare MR, then nn + each SMD method at its threshold); all
columns show the SAME slice
(the nn warp's largest-TV cross-section per view), ZOOMED to the TV
bounding box + margin (the NN box, shared by all columns); only the
1-voxel TV contour is drawn (yellow, no fill), plus the volume +
% of the original CT-grid TV volume in the column header.
Outputs (default results/label_interp_qa/<subject>_<date>/):
<reg>_T.h5 registration transform (cached; re-runs skip it)
<reg>__<warp_key>_tv.nii.gz warped TV on the T1c grid (warp_key includes thr)
<subject>_<date>_label_qa_<reg>.png one 3x5 figure per registration
summary.json per-combination metrics (nvox, mm3, dice, shifts)
+ the SMD thresholding volume-preservation check
Usage (repo root):
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
python scripts/qa_label_interpolation.py
python scripts/qa_label_interpolation.py --subject 2BBN5XA4 --date 20250411 \
--regs ants_rigid ants_affine ants_syn
"""
import argparse
import json
import os
import sys
import time
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
sys.path.insert(0, ROOT)
import numpy as np
import SimpleITK as sitk
from scipy import ndimage
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from src.common import d, read_nii_img
from src import registration_pipelines as RP
M6 = RP._load_m6()
REGS = ["ants_rigid", "ants_affine", "ants_syn"]
SMD_THRESHOLDS = [0.5] # SMD mask = (resampled smd <= thr); 0.0 / 1.0 were
# tested earlier (36% under / 40% over volume) and dropped from the figure
METHOD_LABELS = {"nn": "NN", "smd_linear": "SMD + linear", "smd_bspline": "SMD + B-spline(3)"}
VIEWS = ["axial", "coronal", "sagittal"]
SMD_DEFAULT = 1e5 # out-of-grid fill for the SMD volume: far OUTSIDE (inside < 0)
# figure / warp columns: NN, then each SMD method at each threshold
COLS = ([("nn", None)]
+ [(m, t) for m in ("smd_linear", "smd_bspline") for t in SMD_THRESHOLDS])
ZOOM_MARGIN = 12 # voxels of context around the TV bbox in the displayed slices
def warp_key(method, thr):
return method if thr is None else f"{method}_thr{thr:g}"
def col_label(method, thr):
if thr is None:
return METHOD_LABELS[method]
return f"{METHOD_LABELS[method]} (smd<={thr:g})"
def reg_call(reg, ct, t1c):
"""ANTsPy registration, fixed = T1c / moving = CT (the shared RP rule)."""
if reg == "ants_rigid":
return RP.register_antspyx(ct, t1c)
if reg == "ants_affine":
return RP.register_ants_affine_variant(ct, t1c, "antspyx")
if reg == "ants_syn":
return RP.register_ants_variant_deformable(ct, t1c, "ants_syn")
raise KeyError(f"unknown registration {reg!r} (one of {REGS})")
def find_case(subject, date):
labeled, _unlabeled = M6.discover()
for r in labeled:
if r.get("subject") == subject and r.get("date") == date:
return r
raise SystemExit(f"case {subject}/{date} not found among m6 labeled rows")
def get_transform(reg, ct, t1c, case_dir, use_cache):
"""Register (or load the cached transform), return (T, meta)."""
tf = os.path.join(case_dir, reg + "_T.h5")
if use_cache and os.path.exists(tf):
print(f"[{reg}] loading cached transform {tf}", flush=True)
return sitk.ReadTransform(tf), {"cached": True}
t0 = time.time()
T, meta = reg_call(reg, ct, t1c)
dt = time.time() - t0
sitk.WriteTransform(T, tf)
print(f"[{reg}] registered in {dt:.0f}s ({meta.get('toolkit')}) "
f"inverted={meta.get('auto_inverted')} gap={meta.get('skull_ring_gap_mm')}",
flush=True)
return T, {**meta, "seconds_wall": round(dt, 1), "cached": False}
def make_tv_bin(tv):
"""Binary 0/1 TV image on the TV's own (CT) grid."""
a = (sitk.GetArrayFromImage(tv) > 0.5).astype(np.uint8)
img = sitk.GetImageFromArray(a)
img.CopyInformation(tv)
return img
def warp_tv(tv_bin, t1c, T, method, thr=0.0):
"""Warp the binary TV onto the t1c grid with T (T1c -> CT). Returns a
(z, y, x) uint8 0/1 mask (the interpolated label in MR coordinates).
For the SMD methods the threshold is mask = (smd <= thr)."""
rs = sitk.ResampleImageFilter()
rs.SetReferenceImage(t1c)
rs.SetTransform(T)
if method == "nn":
rs.SetInterpolator(sitk.sitkNearestNeighbor)
rs.SetDefaultPixelValue(0)
rs.SetOutputPixelType(sitk.sitkUInt8)
return (sitk.GetArrayFromImage(rs.Execute(tv_bin)) > 0).astype(np.uint8)
# signed distance to the TV boundary in mm: negative inside, positive
# outside. Note the SimpleITK 2.5 signature is
# (image1, insideIsPositive=False, squaredDistance=True, useImageSpacing=False,
# backgroundValue=0.0) -- keywords are mandatory: the defaults would give
# a SQUARED distance map, which interpolates wrongly. On this build the
# boundary (outer shell) voxels of the foreground get exactly 0, so the
# mask is recovered with <= 0, not < 0.
smd = sitk.SignedMaurerDistanceMap(
tv_bin, insideIsPositive=False, squaredDistance=False, useImageSpacing=True)
smd = sitk.Cast(smd, sitk.sitkFloat32)
if method == "smd_linear":
rs.SetInterpolator(sitk.sitkLinear)
else:
rs.SetInterpolator(sitk.sitkBSpline3) # B-spline order 3
rs.SetDefaultPixelValue(SMD_DEFAULT)
rs.SetOutputPixelType(sitk.sitkFloat32)
return (sitk.GetArrayFromImage(rs.Execute(smd)) <= thr).astype(np.uint8)
def best_slices(mask, shape):
"""(idx, area) of the largest TV cross-section per orthogonal view."""
if not mask.any():
return {"axial": (shape[0] // 2, 0),
"coronal": (shape[1] // 2, 0),
"sagittal": (shape[2] // 2, 0)}
az = mask.sum(axis=(1, 2)) # per z (axial)
ay = mask.sum(axis=(0, 2)) # per y (coronal)
ax = mask.sum(axis=(0, 1)) # per x (sagittal)
return {"axial": (int(np.argmax(az)), int(az.max())),
"coronal": (int(np.argmax(ay)), int(ay.max())),
"sagittal": (int(np.argmax(ax)), int(ax.max()))}
def tv_crop_box(mask, shape, margin):
"""(lo, hi) 3-D index box around the mask (per axis of (z, y, x)) with a
margin; centered box when the mask is empty."""
shape = np.array(shape)
if mask.any():
idx = np.argwhere(mask)
lo, hi = idx.min(0), idx.max(0)
else:
lo = hi = shape // 2
lo = np.clip(lo - margin, 0, shape - 1)
hi = np.clip(hi + margin, 0, shape - 1)
return lo, hi
def view_crop(arr, kind, idx, lo_b, hi_b):
"""Zoomed slice display around the TV box (bench_registration.make_qa_png
convention: coronal / sagittal flipped so superior is up)."""
if kind == "axial":
s = arr[idx, lo_b[1]:hi_b[1] + 1, lo_b[2]:hi_b[2] + 1]
elif kind == "coronal":
s = arr[lo_b[0]:hi_b[0] + 1, idx, lo_b[2]:hi_b[2] + 1][::-1, :]
else:
s = arr[lo_b[0]:hi_b[0] + 1, lo_b[1]:hi_b[1] + 1, idx][::-1, :]
return np.ascontiguousarray(s)
def dice3(a, b):
a, b = a > 0, b > 0
if not a.any() and not b.any():
return 1.0
return float(2 * (a & b).sum() / max(a.sum() + b.sum(), 1))
def centroid_mm(mask, spacing):
c = np.array(ndimage.center_of_mass(mask)) * np.array(spacing)
return c
def check_volume(tv_bin, t1c, source_vol_mm3):
"""Does thresholding the SMD preserve the TV volume? Resampled with the
IDENTITY transform (CT grid -> MR grid, no registration), which isolates
the interpolation + thresholding step from registration error. One row
per figure column (NN + each SMD method x threshold)."""
T = sitk.Transform(3, sitk.sitkIdentity)
sp = [float(s) for s in t1c.GetSpacing()]
vv = sp[0] * sp[1] * sp[2]
out = {}
for method, thr in COLS:
mask = warp_tv(tv_bin, t1c, T, method, thr or 0.0)
out[warp_key(method, thr)] = round(float(mask.sum() * vv), 1)
print(f"volume preservation (identity resample CT grid -> MR grid): "
f"source TV = {source_vol_mm3:.0f} mm3", flush=True)
for tag, v in out.items():
print(f" {tag:20s}: {v:8.1f} mm3 ({100.0 * v / source_vol_mm3:5.1f}%)",
flush=True)
return out
def make_png(mr, combos, out_png, title, spacing, window, source_vol):
"""One registration: 3 rows (views) x (len(combos)+1) columns; the
leftmost column is the bare MR (no overlay). All columns show the SAME
slices (the NN warp's largest-area slice per view); volumes go in the
column headers (row-independent): the original CT TV volume under the
bare-MR column, "% of it" under each warp column."""
lo, hi = window
ncol, nrow = len(combos) + 1, len(VIEWS)
cw, ch, gap, m = 3.6, 3.4, 0.15, 0.62
top = 1.5
fig_w = ncol * cw + (ncol - 1) * gap + 2 * m
fig_h = nrow * ch + (nrow - 1) * gap + m + top
fig = plt.figure(figsize=(fig_w, fig_h))
fig.patch.set_facecolor("k")
aspects = {"axial": spacing[1] / spacing[0],
"coronal": spacing[2] / spacing[0],
"sagittal": spacing[2] / spacing[1]}
cmap_edge = ListedColormap(["#ffee00"])
def w(a):
return np.clip((a - lo) / max(hi - lo, 1e-9), 0.0, 1.0)
for r, kind in enumerate(VIEWS):
for c in range(ncol):
x0 = (m + c * (cw + gap)) / fig_w
y0 = (m + (nrow - 1 - r) * (ch + gap)) / fig_h
axi = fig.add_axes([x0, y0, cw / fig_w, ch / fig_h])
# every column (incl. the bare-MR one) shows the same slice/box
ref = combos[0]["views"][kind]
axi.imshow(w(ref["mr"]), cmap="gray", aspect=aspects[kind])
if c > 0:
v = combos[c - 1]["views"][kind]
m2 = (v["mask"] > 0)
if m2.any():
cont = m2 & ~ndimage.binary_erosion(m2)
axi.imshow(np.ma.masked_where(~cont, cont), cmap=cmap_edge,
aspect=aspects[kind])
else:
axi.text(0.5, 0.5, "empty warp", transform=axi.transAxes,
ha="center", va="center", color="r", fontsize=8)
axi.text(0.02, 0.02, f"{kind}={ref['idx']}", transform=axi.transAxes,
color="w", fontsize=7,
bbox=dict(facecolor="k", alpha=0.55, pad=1.0))
axi.set_xticks([]); axi.set_yticks([])
for s_ in axi.spines.values():
s_.set_edgecolor("0.4")
# column headers: bare MR + original CT TV volume, then method + % volume
cx0 = (m + cw / 2) / fig_w
fig.text(cx0, 1 - 0.74 / fig_h, "MR (no overlay)", color="#ffffff",
fontsize=9.5, ha="center")
fig.text(cx0, 1 - 1.02 / fig_h, f"CT TV: {source_vol:.0f} mm3",
color="#b6f0b6", fontsize=9, ha="center")
for c, comb in enumerate(combos):
cx = (m + (c + 1) * (cw + gap) + cw / 2) / fig_w
fig.text(cx, 1 - 0.74 / fig_h, comb["ilabel"], color="#9fd0ff",
fontsize=9.5, ha="center")
fig.text(cx, 1 - 1.02 / fig_h,
f"{comb['vol_mm3']:.0f} mm3 ({comb['pct_of_source']:.0f}%)",
color="#b6f0b6", fontsize=9, ha="center")
# view labels (left)
for r, kind in enumerate(VIEWS):
y0 = (m + (nrow - 1 - r) * (ch + gap)) / fig_h
fig.text(0.012, y0 + ch / 2 / fig_h, kind.capitalize(), color="w", fontsize=10,
ha="center", va="center", rotation=90)
fig.suptitle(title, color="w", fontsize=11, y=1 - 0.16 / fig_h)
fig.savefig(out_png, facecolor="k", dpi=100)
plt.close(fig)
def main():
ap = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--subject", default="2BBN5XA4")
ap.add_argument("--date", default="20250411")
ap.add_argument("--regs", nargs="*", default=REGS,
help="subset of: " + " ".join(REGS))
ap.add_argument("--tv", default=None,
help="explicit TV file (default: first TV-*.nii.gz in RT)")
ap.add_argument("--out", default=None,
help="output dir (default results/label_interp_qa/<subject>_<date>)")
ap.add_argument("--no-cache", action="store_true",
help="re-register and re-warp even if cached outputs exist")
args = ap.parse_args()
row = find_case(args.subject, args.date)
case_dir = args.out or os.path.join(ROOT, "results", "label_interp_qa",
f"{args.subject}_{args.date}")
os.makedirs(case_dir, exist_ok=True)
ct_path = row["ct"]
t1c_path = row["img"]
tv_path = args.tv or row["label"]
print(f"case {row['key']}: ct={os.path.basename(ct_path)} "
f"t1c={os.path.basename(t1c_path)} tv={os.path.basename(tv_path)}", flush=True)
ct = read_nii_img(ct_path)
t1c = read_nii_img(t1c_path)
tv = read_nii_img(tv_path)
tv_bin = make_tv_bin(tv)
mr = sitk.GetArrayFromImage(t1c).astype("float32")
spacing = [float(s) for s in t1c.GetSpacing()] # (x, y, z)
print(f"t1c {t1c.GetSize()} @ {spacing}mm, tv nvox={int((sitk.GetArrayFromImage(tv_bin) > 0).sum())} "
f"({int((sitk.GetArrayFromImage(tv_bin) > 0).sum()) * ct.GetSpacing()[0] * ct.GetSpacing()[1] * ct.GetSpacing()[2]:.0f} mm3 on the CT grid)",
flush=True)
# MR display window from the head
pos = mr > 0
thr = np.percentile(mr[pos], 2) if pos.any() else 0.0
head = mr > max(thr, 1e-6)
lo, hi = (np.percentile(mr[head], (2, 98)) if head.any() else (0.0, 1.0))
# volume-preservation check for the SMD thresholding (no registration)
sp_ct = [float(s) for s in tv_bin.GetSpacing()]
source_vol = float((sitk.GetArrayFromImage(tv_bin) > 0).sum()) * sp_ct[0] * sp_ct[1] * sp_ct[2]
vol_check = check_volume(tv_bin, t1c, source_vol)
combos = []
reg_meta = {}
for reg in args.regs:
T, meta = get_transform(reg, ct, t1c, case_dir, use_cache=not args.no_cache)
reg_meta[reg] = meta
reg_combos = []
# warp every column (NN + each SMD method x threshold)
for method, thr in COLS:
key = warp_key(method, thr)
p = os.path.join(case_dir, f"{reg}__{key}_tv.nii.gz")
if not args.no_cache and os.path.exists(p):
print(f"[{reg}/{key}] loading cached warp {p}", flush=True)
mask = (sitk.GetArrayFromImage(read_nii_img(p)) > 0).astype(np.uint8)
dt = 0.0
else:
t0 = time.time()
mask = warp_tv(tv_bin, t1c, T, method, thr or 0.0)
dt = time.time() - t0
mimg = sitk.GetImageFromArray(mask)
mimg.CopyInformation(t1c)
sitk.WriteImage(mimg, p, True)
nvox = int(mask.sum())
reg_combos.append({"reg": reg, "method": method, "thr": thr,
"key": key, "ilabel": col_label(method, thr),
"mask": mask, "nvox": nvox,
"vol_mm3": round(nvox * spacing[0] * spacing[1] * spacing[2], 1),
"seconds_warp": round(dt, 1)})
print(f"[{reg}/{key}] nvox={nvox} ({reg_combos[-1]['vol_mm3']} mm3) ({dt:.1f}s)",
flush=True)
# all columns show the SAME slices (the NN warp's largest-area slice
# per view) and the SAME zoom box (NN mask bbox + margin)
nn = next(c for c in reg_combos if c["method"] == "nn")
sl = best_slices(nn["mask"], mr.shape)
lo_b, hi_b = tv_crop_box(nn["mask"], mr.shape, ZOOM_MARGIN)
for c in reg_combos:
c["slices"] = [sl[k][0] for k in VIEWS]
c["views"] = {kind: {"mr": view_crop(mr, kind, sl[kind][0], lo_b, hi_b),
"mask": view_crop(c["mask"], kind, sl[kind][0], lo_b, hi_b),
"idx": sl[kind][0], "area": sl[kind][1]} for kind in VIEWS}
c["pct_of_source"] = 100.0 * c["vol_mm3"] / source_vol
c["pct_vs_nn"] = (100.0 * c["vol_mm3"] / nn["vol_mm3"]
if nn["vol_mm3"] else 0.0)
if c["method"] == "nn":
c["dice_vs_nn"] = 1.0
c["centroid_shift_mm_vs_nn"] = 0.0
else:
c["dice_vs_nn"] = round(dice3(nn["mask"], c["mask"]), 4)
cs = centroid_mm(c["mask"], spacing) - centroid_mm(nn["mask"], spacing)
c["centroid_shift_mm_vs_nn"] = [round(float(x), 2) for x in cs]
combos.extend(reg_combos)
# warped-volume table per registration (% of the original CT-grid TV volume)
parts = " | ".join(f"{c['key']}={c['vol_mm3']:.0f} mm3 ({c['pct_of_source']:.0f}%)"
for c in reg_combos)
print(f"warped volumes {reg} (CT TV={source_vol:.0f} mm3, "
f"slices ax/cor/sag={nn['slices']}): {parts}", flush=True)
pngs = {}
for reg in args.regs:
title = (f"{args.subject} {args.date} -- {reg}: TV warping QA "
f"(same slices = NN largest-area per view)\n"
f"fixed=T1c ({os.path.basename(t1c_path)[:40]}...) moving=CT | "
f"yellow = TV contour | MR window {lo:.0f}..{hi:.0f}")
png = os.path.join(case_dir, f"{args.subject}_{args.date}_label_qa_{reg}.png")
make_png(mr, [c for c in combos if c["reg"] == reg], png, title, spacing,
(float(lo), float(hi)), source_vol)
pngs[reg] = png
print(f"png: {png}", flush=True)
summary = {
"case": row["key"], "subject": args.subject, "date": args.date,
"inputs": {"ct": ct_path, "t1c": t1c_path, "tv": tv_path},
"t1c_geometry": {"size": list(t1c.GetSize()), "spacing": spacing,
"origin": [float(o) for o in t1c.GetOrigin()]},
"regs": {k: {kk: vv for kk, vv in v.items() if kk != "stages"}
for k, v in reg_meta.items()},
"volume_check": {"source_tv_mm3": round(source_vol, 1),
"identity_resample_mm3": vol_check},
"pngs": pngs,
"combos": {f"{c['reg']}/{c['key']}": {
"nvox": c["nvox"], "vol_mm3": c["vol_mm3"],
"pct_of_source_tv": round(c["pct_of_source"], 1),
"slices_ax_cor_sag": c["slices"],
"dice_vs_nn": c["dice_vs_nn"],
"centroid_shift_mm_vs_nn": c["centroid_shift_mm_vs_nn"],
"seconds_warp": c["seconds_warp"]} for c in combos},
}
sf = os.path.join(case_dir, "summary.json")
with open(sf, "w") as f:
json.dump(summary, f, indent=1)
print(f"summary: {sf}", flush=True)
print("done", flush=True)
if __name__ == "__main__":
main()

340
src/registration_pipelines.py
View file

@ -359,12 +359,15 @@ def _itk_transform_to_sitk_rigid(itk_tf):
return t return t
def register_elastix(ct, t1c, workdir=None, threads=8): def register_elastix(ct, t1c, workdir=None, threads=8,
"""ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI, full images. transform="EulerTransform"):
"""ITK-Elastix rigid/affine CT->T1c (fixed=T1c, moving=CT), Mattes MI,
full images.
Runs in-process through the ITK-Elastix python API (pip ``itk`` + Runs in-process through the ITK-Elastix python API (pip ``itk`` +
``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the ``itk-elastix``, elastix 5.x). ``transform`` = ``EulerTransform`` (rigid
5.x replacement for ``Rigid3DTransform``) + Mattes MI 6-DOF, the 5.x replacement for ``Rigid3DTransform``) or
``AffineTransform`` (12-DOF, the affine benchmark) + Mattes MI
(``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the full (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the full
CT / T1c volumes (no head box). No manual centroid init: elastix's CT / T1c volumes (no head box). No manual centroid init: elastix's
``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity ``AutomaticTransformInitialization`` (CenterOfGravity) aligns the intensity
@ -386,7 +389,7 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
"ImageDimension": 3, "ImageDimension": 3,
"FixedInternalImagePixelType": "float", "FixedInternalImagePixelType": "float",
"MovingInternalImagePixelType": "float", "MovingInternalImagePixelType": "float",
"Transform": "EulerTransform", "Transform": transform,
"TransformBase": "ElastixTransformBase", "TransformBase": "ElastixTransformBase",
"FixedCenteredTransform": True, "FixedCenteredTransform": True,
"TransformParameterReadoutFormat": 1, "TransformParameterReadoutFormat": 1,
@ -402,7 +405,9 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
"ImageSamplerPatchSize": [25, 25, 25], "ImageSamplerPatchSize": [25, 25, 25],
"WriteResultImage": False, "WriteResultImage": False,
} }
outdir = os.path.join(workdir or d("data/m6reg/elastix"), "itkelastix") outdir = os.path.join(workdir or d("data/m6reg/elastix"),
"itkelastix" if transform == "EulerTransform"
else "itkelastix_aff")
os.makedirs(outdir, exist_ok=True) os.makedirs(outdir, exist_ok=True)
import time, io, contextlib import time, io, contextlib
t0 = time.time() t0 = time.time()
@ -435,7 +440,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
meta = { meta = {
"pipeline": "elastix", "toolkit": version, "pipeline": "elastix", "toolkit": version,
"metric": "MattesMutualInformation (full volumes)", "metric": "MattesMutualInformation (full volumes)",
"transform": "Rigid (EulerTransform, ASGD, 3 resolutions, auto center init)", "transform": (f"{'Rigid' if transform == 'EulerTransform' else 'Affine'} "
f"({transform}, ASGD, 3 resolutions, auto center init)"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": _modality_name(fixed), "moving": _modality_name(moving), "fixed": _modality_name(fixed), "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "seconds_reg": round(dt, 1),
@ -443,6 +449,13 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
return T, meta return T, meta
def register_elastix_affine(ct, t1c, workdir=None, threads=8):
"""Affine (12-DOF) variant of :func:`register_elastix` — same budget,
``AffineTransform`` instead of ``EulerTransform``."""
return register_elastix(ct, t1c, workdir=workdir, threads=threads,
transform="AffineTransform")
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# 3. ANTsPy # 3. ANTsPy
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
@ -738,14 +751,40 @@ def _impact_models_host():
return list(dict.fromkeys(rels)) return list(dict.fromkeys(rels))
def _impact_param_map(gpu, iters, samples, voxel_sizes): def _itk_affine_txt_to_sitk(path):
"""Build the Elastix-5.2 parameter map for a RIGID (Euler) IMPACT """Parse an ITK ``AffineTransform`` .txt file into a SimpleITK
registration: the MIND + TotalSegmentator static feature combination of AffineTransform (fallback for when SimpleITK's own reader does not accept
the ImpactReg ``MR_CT_TS`` reference map (one model pair per level, the file). ``Parameters`` = 9 matrix values (row-major) + 3 translation,
coarse -> fine, see :data:`IMPACT_LEVELS`) + a secondary Mattes-MI metric plus a ``Center`` line; the mapping is p -> M (p - c) + c + t."""
(:data:`IMPACT_MI_WEIGHT`). The image pyramid is NOT resampled (all-1 text = open(path).read()
schedules); the scale comes from the per-level ImpactVoxelSize (as the if "AffineTransform" not in text:
IMPACT docs recommend). ``iters`` is a per-level tuple or a single int. raise RuntimeError(f"{path} is not an AffineTransform file")
body = text.split("Parameters:", 1)[1].split("FixedParameters:", 1)[0]
params = [float(x) for x in body.split()]
if len(params) != 12:
raise RuntimeError(f"expected 12 affine parameters in {path}, "
f"got {len(params)}")
m = re.search(r"^Center:\s*([\d\.\-e+ ]+)", text, re.M)
center = ([float(x) for x in m.group(1).split()[:3]] if m
else [0.0, 0.0, 0.0])
t = sitk.AffineTransform(3)
t.SetCenter(tuple(center))
t.SetMatrix(params[0:9])
t.SetTranslation(params[9:12])
return t
def _impact_param_map(gpu, iters, samples, voxel_sizes,
transform="EulerTransform"):
"""Build the Elastix-5.2 parameter map for a RIGID (Euler) or AFFINE
IMPACT registration: the MIND + TotalSegmentator static feature
combination of the ImpactReg ``MR_CT_TS`` reference map (one model pair
per level, coarse -> fine, see :data:`IMPACT_LEVELS`) + a secondary
Mattes-MI metric (:data:`IMPACT_MI_WEIGHT`). The image pyramid is NOT
resampled (all-1 schedules); the scale comes from the per-level
ImpactVoxelSize (as the IMPACT docs recommend). ``iters`` is a per-level
tuple or a single int; ``transform`` = ``EulerTransform`` (6-DOF) or
``AffineTransform`` (12-DOF, the affine benchmark).
""" """
n = len(voxel_sizes) n = len(voxel_sizes)
if isinstance(iters, (int,)): if isinstance(iters, (int,)):
@ -790,8 +829,9 @@ def _impact_param_map(gpu, iters, samples, voxel_sizes):
'(CheckNumberOfSamples "true")', '(CheckNumberOfSamples "true")',
"(RequiredRatioOfValidSamples 0.1)", "(RequiredRatioOfValidSamples 0.1)",
"", "",
"// rigid 6-DOF to be optimized (auto center init, no -t0)", f"// {'rigid 6-DOF' if transform == 'EulerTransform' else 'affine 12-DOF'} "
'(Transform "EulerTransform")', "to be optimized (auto center init, no -t0)",
f'(Transform "{transform}")',
'(TransformBase "ElastixTransformBase")', '(TransformBase "ElastixTransformBase")',
'(TransformParameterReadoutFormat "AllParametersOneLine")', '(TransformParameterReadoutFormat "AllParametersOneLine")',
'(FixedCenteredTransform "true")', '(FixedCenteredTransform "true")',
@ -849,8 +889,8 @@ def _impact_param_map(gpu, iters, samples, voxel_sizes):
def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS,
samples=IMPACT_DEFAULT_SAMPLES, threads=12, samples=IMPACT_DEFAULT_SAMPLES, threads=12,
keep_scratch=False): keep_scratch=False, transform="EulerTransform"):
"""IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image """IMPACT rigid/affine CT->T1c (fixed=T1c, moving=CT) in the Docker image
``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric). ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic* IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
@ -930,7 +970,7 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS,
os.path.join(datadir, "Fixed_image.mha"), True) os.path.join(datadir, "Fixed_image.mha"), True)
sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
os.path.join(datadir, "Moving_image.mha"), True) os.path.join(datadir, "Moving_image.mha"), True)
param_text = _impact_param_map(gpu, iters, samples, voxel_sizes) param_text = _impact_param_map(gpu, iters, samples, voxel_sizes, transform)
with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
f.write(param_text) f.write(param_text)
@ -968,11 +1008,14 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS,
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}") raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
# the optimized Euler (world frame) is the full combination transform # the optimized transform (world frame) is the full combination
# (AutomaticTransformInitialization is folded in by elastix); with # (AutomaticTransformInitialization is folded in by elastix); with
# fixed=T1c / moving=CT the direction it holds is resolved empirically # fixed=T1c / moving=CT the direction it holds is resolved
# by fix_orientation against the skull ring. # empirically by fix_orientation against the skull ring.
T_main = sitk.ReadTransform(itk_tf) try:
T_main = sitk.ReadTransform(itk_tf)
except Exception: # noqa: SimpleITK reader lacks this ITK format
T_main = _itk_affine_txt_to_sitk(itk_tf)
T, inverted = fix_orientation(ct, t1c, T_main) T, inverted = fix_orientation(ct, t1c, T_main)
finally: finally:
if not keep_scratch: if not keep_scratch:
@ -984,7 +1027,8 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS,
"toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
"metric": (f"IMPACT semantic ({models}, L1+Dice, Static) + " "metric": (f"IMPACT semantic ({models}, L1+Dice, Static) + "
f"Mattes MI (w {IMPACT_MI_WEIGHT})"), f"Mattes MI (w {IMPACT_MI_WEIGHT})"),
"transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; " "transform": (f"{'Rigid' if transform == 'EulerTransform' else 'Affine'} "
f"({transform}, ASGD, {n} scales, auto center init; "
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": _modality_name(fixed), "moving": _modality_name(moving), "fixed": _modality_name(fixed), "moving": _modality_name(moving),
@ -998,6 +1042,17 @@ def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS,
return T, meta return T, meta
def register_impact_affine(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS,
samples=IMPACT_DEFAULT_SAMPLES, threads=12,
keep_scratch=False):
"""Affine (12-DOF) variant of :func:`register_impact` — the MIND+TS
reference map with ``AffineTransform`` instead of ``EulerTransform``,
same budgets and M6 adaptations."""
return register_impact(ct, t1c, gpu=gpu, iters=iters, samples=samples,
threads=threads, keep_scratch=keep_scratch,
transform="AffineTransform")
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# 6. NiftyReg (KCL-BMEIS suite, reg_aladin) # 6. NiftyReg (KCL-BMEIS suite, reg_aladin)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
@ -1064,7 +1119,8 @@ def _parse_niftyreg_affine(path):
def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False, def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=False,
backend="auto", gpu=None, clip_ct=(0.0, 80.0)): backend="auto", gpu=None, clip_ct=(0.0, 80.0),
rigid_only=True):
"""NiftyReg rigid CT->T1c (reference=T1c, floating=CT) via ``reg_aladin``. """NiftyReg rigid CT->T1c (reference=T1c, floating=CT) via ``reg_aladin``.
reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global reg_aladin (Modat et al., J. Med. Img. 1(2) 024003, 2014) is a global
@ -1147,10 +1203,12 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
os.path.join(tmp, "ref.nii.gz"), True) os.path.join(tmp, "ref.nii.gz"), True)
sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
os.path.join(tmp, "flo.nii.gz"), True) os.path.join(tmp, "flo.nii.gz"), True)
flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz", flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz"]
"-rigOnly", "-aff", "aff.txt", if rigid_only:
"-ln", str(levels), "-maxit", str(maxit), flags.append("-rigOnly")
"-interp", interp, "-omp", str(threads)] flags += ["-aff", "aff.txt",
"-ln", str(levels), "-maxit", str(maxit),
"-interp", interp, "-omp", str(threads)]
if platf is not None: if platf is not None:
flags += ["-platf", platf] flags += ["-platf", platf]
if use_gpu: if use_gpu:
@ -1194,8 +1252,9 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
"toolkit": toolkit, "toolkit": toolkit,
"backend": "docker" if use_docker else "conda", "backend": "docker" if use_docker else "conda",
"metric": "NCC block matching + trimmed least squares (symmetric, full volumes)", "metric": "NCC block matching + trimmed least squares (symmetric, full volumes)",
"transform": (f"Rigid (reg_aladin -rigOnly, {levels}-level pyramid, " "transform": (f"{'Rigid' if rigid_only else 'Affine'} (reg_aladin"
f"{maxit} TLS iters/level, cubic interp)"), f"{' -rigOnly' if rigid_only else ''}, {levels}-level "
f"pyramid, {maxit} TLS iters/level)"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": fixed_name, "moving": moving_name, "fixed": fixed_name, "moving": moving_name,
"seconds_reg": round(dt, 1), "seconds_reg": round(dt, 1),
@ -1207,6 +1266,17 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
return T, meta return T, meta
def register_niftyreg_affine(ct, t1c, threads=24, levels=3, maxit=15,
keep_scratch=False, backend="auto", gpu=None,
clip_ct=(0.0, 80.0)):
"""Affine (12-DOF) variant of :func:`register_niftyreg` — ``reg_aladin``
full affine (no ``-rigOnly``), same pyramid / TLS budget."""
return register_niftyreg(ct, t1c, threads=threads, levels=levels,
maxit=maxit, keep_scratch=keep_scratch,
backend=backend, gpu=gpu, clip_ct=clip_ct,
rigid_only=False)
REGISTERS = { REGISTERS = {
"antspyx": register_antspyx, "antspyx": register_antspyx,
"elastix": register_elastix, "elastix": register_elastix,
@ -1218,6 +1288,203 @@ REGISTERS = {
ALL_PIPELINES = list(REGISTERS) ALL_PIPELINES = list(REGISTERS)
# --------------------------------------------------------------------------- #
# Affine variants (see REGISTRATION_BENCHMARK_AFFINE.md)
# --------------------------------------------------------------------------- #
#
# Same unified contract as the rigid pipelines, but every toolkit estimates a
# full affine (12-DOF) transform. The antspyx rows run stock linear-only
# presets; the other toolkits keep their rigid-benchmark budgets with an
# affine transform in place of the rigid one. All return (T, meta) with T an
# affine (or composite of affines) mapping T1c -> CT, direction resolved by
# fix_orientation like the rigid pipelines.
ANT_AFFINE_PRESETS = {
# pipeline name -> (antspyx type_of_transform, description)
"antspyx": (
"Affine",
"direct Affine preset: Affine[0.25] (Mattes 32, 2100x1200x1200x10, "
"6x4x2x1, 3x2x1x0); COM init"),
"ants_trsaa": (
"TRSAA",
"TRSA cascade: Translation[1] -> Rigid[1] -> Similarity[1] (each "
"Mattes 32, 2000x2000x0, 4x2x1, 2x1x0) -> Affine[1] x2 (Mattes 32, "
"40x20x0, 4x2x1, 2x1x0); COM init"),
"ants_reproa": (
"antsRegistrationSyNRepro[a]",
"Repro (a): Rigid[0.1] + Affine[0.1] (each GC 1/1/Regular/0.25, "
"1000x500x250x100, 8x4x2x1, 3x2x1x0vox); COM init"),
"ants_quickreproa": (
"antsRegistrationSyNQuickRepro[a]",
"QuickRepro (a): Rigid[0.1] + Affine[0.1] (each GC 1/1/Regular/0.25, "
"1000x500x250x0, 8x4x2x1, 3x2x1x0vox); COM init"),
}
def register_ants_affine_variant(ct, t1c, preset_name):
"""Run one stock antspyx linear-only preset (:data:`ANT_AFFINE_PRESETS`)
as an affine pipeline under the unified contract (fixed = T1c, moving =
CT, full volumes, no clip). All stages are linear, so the stage files
ANTs writes (in numeric order = the order it applies them) are composed
into the contract affine transform; binary merged-linear .mat files
(Ants combines consecutive linear stages) are recovered by point-probing
(:func:`_ants_stage_to_sitk`, exact for affine maps). OMP threads are
left at the system default, as for the other antspyx rows."""
if preset_name not in ANT_AFFINE_PRESETS:
raise KeyError(f"unknown antspyx affine preset {preset_name!r} "
f"(one of {list(ANT_AFFINE_PRESETS)})")
tx_type, desc = ANT_AFFINE_PRESETS[preset_name]
try:
import ants
except Exception as e: # noqa
raise PipelineUnavailable(f"antspyx: {e!r}")
import glob
import time
fx, mv = fixed_moving(ct, t1c)
af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32))
am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32))
tmp = tempfile.mkdtemp(prefix=f"ants_aff_{preset_name}_")
dt = 0.0
stage_names = []
try:
t0 = time.time()
r = ants.registration(fixed=af, moving=am, type_of_transform=tx_type,
outprefix=tmp + "/")
dt = time.time() - t0
stage_files = [f for f in glob.glob(os.path.join(tmp, "*"))
if f.endswith(".mat")]
if not stage_files:
raise RuntimeError(f"{preset_name}: no .mat stage files in {tmp}")
stage_files.sort(key=lambda p: int(re.match(
r"(\d+)", os.path.basename(p)).group(1)))
stage_names = [os.path.basename(p) for p in stage_files]
stages = [_ants_stage_to_sitk(p, fx) for p in stage_files]
T = stages[0] if len(stages) == 1 else sitk.CompositeTransform(stages)
T, inverted = fix_orientation(ct, t1c, T)
finally:
shutil.rmtree(tmp, ignore_errors=True)
meta = {
"pipeline": preset_name,
"toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
"metric": "ANTs stock preset metrics (Mattes MI / GC per stage)",
"transform": f"Affine (antspyx preset '{tx_type}'): {desc}",
"convention": "T1c->CT (verified)", "auto_inverted": bool(inverted),
"fixed": _modality_name(fx), "moving": _modality_name(mv),
"seconds_reg": round(dt, 1),
"stage_files": stage_names,
}
return T, meta
def register_simpleitk_affine(ct, t1c, threads=12, iters_full=2000):
"""SimpleITK affine CT->T1c, Mattes MI, full volumes: the rigid
benchmark's Euler3D stage (``iters_full`` iters/level on the [4,2,1]
pyramid, identity start) followed by an ``AffineTransform`` stage with
the same budget, initialised to reproduce the rigid stage exactly so the
affine refines the pose instead of re-estimating it (a single-stage affine
from the identity collapses its scale on these pairs). The returned
transform is the composite condensed to a single affine by point mapping
(T1c -> CT)."""
m6 = _load_m6()
_sitk_set_threads(threads)
ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
fixed, moving = fixed_moving(ct, t1c)
shape = np.array(ct.GetSize())
sh = [s for s in [4, 2, 1] if min(shape) / s >= 6] or [1]
sm = [2, 1, 0][:len(sh)]
center = tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint(
((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2,
(fixed.GetSize()[2] - 1) // 2)))
def _method():
rm = sitk.ImageRegistrationMethod()
rm.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32)
rm.SetOptimizerAsRegularStepGradientDescent(
learningRate=10.0, minStep=0.2, numberOfIterations=iters_full,
relaxationFactor=0.5, gradientMagnitudeTolerance=1e-6)
rm.SetOptimizerScalesFromPhysicalShift()
rm.SetShrinkFactorsPerLevel(sh); rm.SetSmoothingSigmasPerLevel(sm)
rm.SetInterpolator(sitk.sitkLinear)
return rm
def _condense_affine(tr):
"""Condense an affine (or composite of affines) to a single
AffineTransform by point mapping (exact for affine maps)."""
o = np.array(tr.TransformPoint((0.0, 0.0, 0.0)), dtype=float)
M = np.empty((3, 3), dtype=float)
for i in range(3):
p = [0.0, 0.0, 0.0]
p[i] = 1.0
M[:, i] = np.array(tr.TransformPoint(tuple(p)), dtype=float) - o
a = sitk.AffineTransform(3)
a.SetMatrix(M.ravel())
a.SetTranslation(o)
return a
import time
t0 = time.time()
init1 = sitk.Euler3DTransform()
init1.SetCenter(center)
init1.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
r1 = _method()
r1.SetInitialTransform(init1, inPlace=False)
tr1 = _condense_affine(r1.Execute(fixed, moving))
init2 = sitk.AffineTransform(3)
init2.SetCenter(tr1.GetCenter())
init2.SetMatrix(tr1.GetMatrix())
init2.SetTranslation(tr1.GetTranslation())
r2 = _method()
r2.SetInitialTransform(init2, inPlace=False)
tr = _condense_affine(r2.Execute(fixed, moving))
dt = time.time() - t0
ha = sitk.GetArrayFromImage(t1c)
gap = m6._skull_ring_gap(_warp_ct(ct, t1c, tr), ha)
meta = {
"pipeline": "simpleitk", "mode": "benchmark",
"toolkit": f"SimpleITK {sitk.Version_VersionString()}",
"metric": "MattesMutualInformation (both stages, full volumes)",
"transform": (f"Affine (Euler3D rigid init + AffineTransform, "
f"{iters_full} iters/level on pyramid {sh} each, "
f"identity start on the rigid, affine initialised to "
f"reproduce the rigid stage)"),
"convention": "T1c->CT (benchmark budget)", "auto_inverted": False,
"fixed": _modality_name(fixed), "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm,
"skull_ring_gap_mm": gap,
}
return tr, meta
def _ants_affine_register(name):
def _reg(ct, t1c):
return register_ants_affine_variant(ct, t1c, preset_name=name)
return _reg
def register_synthmorph_affine(ct, t1c, gpu=None, threads=None, extent=256,
keep_scratch=False, clip_ct=(0.0, 80.0)):
"""Affine (12-DOF) variant of :func:`register_synthmorph` — the
``-m affine`` matrix model."""
return register_synthmorph(ct, t1c, model="affine", gpu=gpu, threads=threads,
extent=extent, keep_scratch=keep_scratch,
clip_ct=clip_ct)
REGISTERS_AFFINE = {
"antspyx": _ants_affine_register("antspyx"),
"ants_trsaa": _ants_affine_register("ants_trsaa"),
"ants_reproa": _ants_affine_register("ants_reproa"),
"ants_quickreproa": _ants_affine_register("ants_quickreproa"),
"elastix": register_elastix_affine,
"impact": register_impact_affine,
"niftyreg": register_niftyreg_affine,
"simpleitk": register_simpleitk_affine,
"synthmorph": register_synthmorph_affine,
}
ALL_PIPELINES_AFFINE = list(REGISTERS_AFFINE)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
# Deformable variants (see REGISTRATION_BENCHMARK_DEFORMABLE.md) # Deformable variants (see REGISTRATION_BENCHMARK_DEFORMABLE.md)
# --------------------------------------------------------------------------- # # --------------------------------------------------------------------------- #
@ -2311,8 +2578,10 @@ REGISTERS_DEFORMABLE.update({n: _ants_variant_register(n) for n in ANT_VARIANTS}
REGISTERS_DEFORMABLE["ants_rigid_synonly"] = register_ants_rigid_synonly_deformable REGISTERS_DEFORMABLE["ants_rigid_synonly"] = register_ants_rigid_synonly_deformable
ALL_PIPELINES_DEFORMABLE = list(REGISTERS_DEFORMABLE) ALL_PIPELINES_DEFORMABLE = list(REGISTERS_DEFORMABLE)
# every pipeline name the harness may dispatch / merge on (rigid + variants) # every pipeline name the harness may dispatch / merge on (rigid + variants)
ALL_PIPELINE_NAMES = (ALL_PIPELINES + list(ANT_VARIANTS) ALL_PIPELINE_NAMES = (
+ ["ants_rigid_synonly"]) ALL_PIPELINES + list(ANT_VARIANTS) + ["ants_rigid_synonly"]
+ [n for n in ALL_PIPELINES_AFFINE if n not in ALL_PIPELINES]
)
def availability(): def availability():
@ -2374,6 +2643,7 @@ def availability():
except Exception as e: # noqa except Exception as e: # noqa
out[name] = (False, repr(e)) out[name] = (False, repr(e))
# the ANTs variants all run through antspyx (same availability as antspyx) # the ANTs variants all run through antspyx (same availability as antspyx)
for name in list(ANT_VARIANTS) + ["ants_rigid_synonly"]: for name in (list(ANT_VARIANTS) + ["ants_rigid_synonly"]
+ [n for n in ALL_PIPELINES_AFFINE if n not in ALL_PIPELINES]):
out[name] = out.get("antspyx", (False, "antspyx not checked")) out[name] = out.get("antspyx", (False, "antspyx not checked"))
return out return out