Compare commits
2 commits
a65189cd33
...
73234db39e
| Author | SHA1 | Date | |
|---|---|---|---|
| 73234db39e | |||
| 16b9828d15 |
6 changed files with 698 additions and 267 deletions
4
.gitignore
vendored
4
.gitignore
vendored
|
|
@ -8,4 +8,6 @@ runs_nnu/
|
|||
runs_monai/
|
||||
__pycache__/
|
||||
|
||||
elastix/
|
||||
elastix/
|
||||
niftyreg/NiftyReg-Ubuntu-v2.1.1.zip
|
||||
niftyreg/NiftyReg-Ubuntu-v2.1.1/
|
||||
|
|
|
|||
|
|
@ -22,36 +22,50 @@ sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0) # GTV onto the MR grid
|
|||
```
|
||||
|
||||
places the CT / GTV concentric with the MR head. In `Resample`'s
|
||||
destination→source convention this `T` maps **T1c → CT**. All pipelines
|
||||
register on the **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.
|
||||
destination→source convention this `T` maps **T1c → CT**. Shared rules for
|
||||
every pipeline:
|
||||
|
||||
- **No mask, no crop / pad**: registration runs on the full input volumes
|
||||
exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
|
||||
padding).
|
||||
- **CT is always the fixed image** (moving = T1c); when both inputs are MR
|
||||
the first argument is kept as fixed.
|
||||
|
||||
All pipelines also 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. One
|
||||
consequence of the no-mask rule: the IMPACT metric only evaluates samples that
|
||||
map inside the moving buffer, so with the CT (head-neck, ~300×300×251 mm) as
|
||||
the fixed image and the tight T1c head crop (~220×220×150 mm) as moving, most
|
||||
of the CT buffer maps outside the T1c buffer (~14% of samples in-buffer at the
|
||||
auto-init). Elastix's sample-validity gate (build default
|
||||
`RequiredRatioOfValidSamples 0.25`) used to abort the M6 pairs on that ("too
|
||||
many samples map outside moving image buffer"); the parameter map now sets
|
||||
`RequiredRatioOfValidSamples 0.1` so the FOV-mismatched runs proceed — see the
|
||||
IMPACT section. The comparison is on each toolkit's native optimizer / metric
|
||||
/ global-search capability.
|
||||
|
||||
## The six pipelines
|
||||
|
||||
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).
|
||||
Listed alphabetically. All register on the **full volumes** (no head box, no
|
||||
mask, no crop / pad) 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** (shared rule above; when both
|
||||
inputs are MR the first argument is fixed). `niftyreg` and `synthmorph`
|
||||
clip the CT to [0, 80] HU before registration (per the CT-clipping A/B
|
||||
below; `clip_ct=None` disables it).
|
||||
|
||||
| # | 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) |
|
||||
| 3 | `impact` | **CT / T1c** | auto CenterOfGravity | **IMPACT semantic** (MRSegmentator `MRSeg/MRSeg.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (Docker + GPU; sample gate 0.1 — see IMPACT section) |
|
||||
| 4 | `niftyreg` | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (Docker `niftyreg:2.1.1-ubuntu`, GPU `-platf 1`) |
|
||||
| 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) |
|
||||
|
||||
|
|
@ -62,10 +76,10 @@ level), with the sampler and metric left at their defaults
|
|||
`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
|
||||
identity start — the default, comparable to the other tools) and
|
||||
`--simpleitk-mode production` (the full multi-start tilt grid from
|
||||
`02_build_m6_dataset.py`, ~20× slower, the most robust but not used for the
|
||||
head-to-head).
|
||||
`02_build_m6_dataset.py`, now also on the full volumes per the no-crop rule,
|
||||
~20× slower, the most robust but not used for the head-to-head).
|
||||
|
||||
## Test cases
|
||||
|
||||
|
|
@ -108,40 +122,45 @@ longer reported in the tables/plots):
|
|||
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`).
|
||||
production acceptance rule; aligned ⇔ `|mid| ≤ 20 mm` and `std ≤ 12 mm` —
|
||||
a small negative `mid` is fine: on full-volume registration the CT skull
|
||||
ring can land slightly inside the MR silhouette, which is equally
|
||||
concentric).
|
||||
|
||||
## Results (10 cases)
|
||||
## Results (10 cases, no clip)
|
||||
|
||||
Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over
|
||||
the MR head; higher = better).
|
||||
the MR head; higher = better) with the **unclipped** CT, so all rows are
|
||||
comparable on the same score scale. Full volumes, no mask / crop / pad,
|
||||
fixed = CT. `niftyreg` and `synthmorph` clip the CT to [0, 80] HU *before
|
||||
registration* (per the CT-clipping A/B below; `clip_ct=None` disables it);
|
||||
`antspyx`, `elastix` and `simpleitk` register the CT as-is. `niftyreg` is the
|
||||
v2.1.1 Docker build on the GPU (`-platf 1`); `impact` cannot run on these
|
||||
pairs under the no-mask rule (see the IMPACT section).
|
||||
|
||||
| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
|
||||
|---|---|---|---|---|---|
|
||||
| **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** |
|
||||
| **antspyx** | ANTsPy 0.6.3 | 9 | 36.6 | 0.93 | **0.653 / 0.644** |
|
||||
| **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.517 / 0.520** |
|
||||
| **impact** | Docker Elastix 5.2 + IMPACT (MRSegmentator) | 8 | 2534 | 0.93 | **0.591 / 0.575** |
|
||||
| **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 20.8 | 0.93 | **0.292 / 0.338** |
|
||||
| **simpleitk** | SimpleITK 2.5.6 | 9 | 10.6 | 0.93 | **0.651 / 0.641** |
|
||||
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 24.5 | 0.93 | **0.277 / 0.328** |
|
||||
|
||||
\* `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):
|
||||
Per-case MI (nats, unclipped score; **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 |
|
||||
| m6_22IRRKU6_20240115 | **0.649** | 0.472 | 0.635 | 0.292 | 0.642 | 0.277 |
|
||||
| m6_23XIMO4M_20230222 | **0.655** | 0.416 | 0.478 | 0.632 | 0.652 | 0.616 |
|
||||
| m6_24V62QQM_20230223 | **0.704** | 0.701 | 0.697 | 0.260 | 0.663 | 0.256 |
|
||||
| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.556 | 0.493 | **0.602** | 0.484 |
|
||||
| m6_24V62QQM_20240820 | **0.653** | 0.644 | 0.531 | 0.118 | 0.652 | 0.111 |
|
||||
| m6_26TYQI3C_20221027 | 0.659 | 0.653 | gate-abort | 0.609 | **0.695** | 0.590 |
|
||||
| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | 0.450 | 0.069 | **0.577** | 0.061 |
|
||||
| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | 0.629 | 0.233 | 0.634 | 0.231 |
|
||||
| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | 0.625 | 0.340 | 0.651 | 0.323 |
|
||||
|
||||
### Interpretation
|
||||
|
||||
|
|
@ -163,52 +182,138 @@ Per-case MI (nats; **bold** = best for the case):
|
|||
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).
|
||||
- **IMPACT (semantic metric) now runs on these pairs** — the no-mask abort was
|
||||
elastix's sample-validity gate (build default `RequiredRatioOfValidSamples`
|
||||
0.25) rejecting the ~14% in-buffer fraction of the CT head-neck FOV vs the
|
||||
tight T1c crop; the parameter map sets it to 0.1 (see the IMPACT section).
|
||||
With MRSegmentator at the 2000/5000 budget it is top-three in MI on 7 of the
|
||||
8 registered cases (median 0.591): 0.635 on the hardest `m6_22IRRKU6`
|
||||
(masked-era 0.624), 0.697 on `m6_24V62QQM_20230223` (near-best), 0.629 on
|
||||
`m6_2C4TKI6H` — and it shares the production-acceptance of the grossly
|
||||
repositioned `m6_2BBN5XA4` (std 10.7 mm) with antspyx/simpleitk, the only
|
||||
non-intensity tool to do so. One run diverges mid-registration
|
||||
(`m6_26TYQI3C`: the in-buffer fraction falls below 10% at iteration 256 and
|
||||
the gate aborts — a stochastic-draw outcome). It remains by far the slowest
|
||||
(~42 min/case at 2000/5000).
|
||||
- **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.
|
||||
is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class;
|
||||
`m6_26TYQI3C` 0.609) but falls into the same tilt local minima
|
||||
(`m6_2BBN5XA4` 0.069, `m6_24V62QQM_20240820` 0.118). The v2.1.1 Docker build
|
||||
on the GPU (`-platf 1`) is ~21 s/case, ~4× slower than elastix (vs ~4 min
|
||||
CPU for the conda 1.5.69 build). It clips the CT to [0,80] HU by default
|
||||
(per the A/B below) — neutral on the unclipped score scale
|
||||
(0.296→0.292 median).
|
||||
- **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.
|
||||
(median 0.277 with its default CT clip; 0.241 unclipped): its `rigid` model
|
||||
is trained on brain *MRI*, so the ~30° tilt defeats it on the repositioned
|
||||
cases (0.06–0.11) even though it matches the intensity tools on the
|
||||
cleanest (`m6_26TYQI3C` 0.590). The [0,80] brain-window clip of the CT is
|
||||
its only genuine A/B win (+0.036 median, 7/9 cases), so it is the default.
|
||||
It needs **no centroid init** (global deep method).
|
||||
- **`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).
|
||||
SimpleITK 10.6 s, **NiftyReg 20.8 s** (v2.1.1, GPU block matching),
|
||||
SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT is by far the
|
||||
slowest — 2534 s/case median at the 2000/5000 benchmark budget, ~276 s/case
|
||||
at the authors' 500/2000; GPU feature extraction with backprop per
|
||||
iteration. No-mask runs are faster per iteration than the masked-era
|
||||
~833 s/case @500/2000 because only the ~14% in-buffer samples are scored.)
|
||||
- **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;
|
||||
the in-process Elastix ASGD is deterministic given its fixed seed;
|
||||
SynthMorph and NiftyReg are deterministic. **IMPACT is the most stochastic**:
|
||||
SynthMorph is deterministic; NiftyReg is deterministic on the CPU
|
||||
(conda 1.5.69) but shows small run-to-run variation on the GPU (v2.1.1,
|
||||
~±0.01–0.05 MI on the weak cases). **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.
|
||||
|
||||
## CT clipping A/B (`--clip-ct 0 80`)
|
||||
|
||||
Question: does clipping the CT to a brain window ([0, 80] HU — skull and air
|
||||
removed) **before registration** improve the result? The harness clips the CT
|
||||
before registration AND before scoring (percentile-binned MI is scale-free,
|
||||
so the A/B compares like for like). Same 10 cases, same budgets.
|
||||
|
||||
Median achieved MI (nats) and wall-clock (9 complete-CT cases):
|
||||
|
||||
| pipeline | no clip | clip [0,80] | Δ (median MI) | clip wins (cases by MI) |
|
||||
|---|---|---|---|---|
|
||||
| antspyx | 0.653 / 36.6 s | 0.466 / 34.7 s | **−0.186** | 0/9 |
|
||||
| elastix | 0.517 / 6.0 s | 0.425 / 5.9 s | **−0.092** | 1/9 |
|
||||
| niftyreg | 0.296 / 23.8 s | 0.320 / 20.2 s | +0.024 | 5/9 |
|
||||
| simpleitk | 0.651 / 10.6 s | 0.436 / 11.1 s | **−0.215** | 0/9 |
|
||||
| synthmorph | 0.241 / 26.0 s | 0.318 / 25.6 s | +0.077 | 5/9 |
|
||||
| impact | abort | abort | - | - (historical: at the time both aborted on the old 0.25 sample-validity gate; the gate is now 0.1, A/B to be re-done) |
|
||||
|
||||
**Verdict: no — clipping the CT to [0, 80] HU does not improve the result; do
|
||||
not clip.**
|
||||
|
||||
- The headline alignment score drops for **all three leading intensity
|
||||
pipelines**: antspyx 0.653→0.466, simpleitk 0.651→0.436, elastix
|
||||
0.517→0.425 (median nats; simpleitk and antspyx lose on **all 9** cases,
|
||||
elastix on 8/9).
|
||||
- The two weaker pipelines show small median-MI gains (niftyreg +0.024,
|
||||
synthmorph +0.077) but only mixed case-by-case (5/9), and their geometric
|
||||
QA does not support a real improvement.
|
||||
- **Geometric concentricity (skull-ring gap) is a wash**: mean mid/std move by
|
||||
~1–5 mm in either direction (niftyreg mid 23.6→18.1 the only clear gain;
|
||||
antspyx std 16.4→18.0 the clear loss). Production acceptance even falls
|
||||
(no-clip 4/4/2/4/2 → clip 1/0/1/0/2).
|
||||
- **Mechanism.** Clipping is a saturating (non-invertible) intensity map: the
|
||||
skull atom (800+ HU → 80) removes the strongest intensity anchor
|
||||
(skull ↔ MR scalp) and compresses the CT rank distribution, so the
|
||||
percentile-binned MI between the warped CT and the MR is inherently lower
|
||||
*for the same pose*. The pose itself does not systematically degrade — the
|
||||
head is rigid, so skull-ring concentricity is its fingerprint, and that is
|
||||
unchanged. The MI drop is therefore mostly a score-scale effect, but it is
|
||||
exactly the score the benchmark reports, and there is no geometric benefit
|
||||
to pay it for.
|
||||
- Scoring note: `_skull_ring_gap` tests the ring with `>= p97` (not `>`) so
|
||||
the orientation gate and the gap QA also work on a clipped CT (where p97
|
||||
sits exactly on the clipped-skull atom).
|
||||
- **Decision / adopted default.** The pipelines were updated: `niftyreg` and
|
||||
`synthmorph` clip the CT to [0, 80] HU before registration
|
||||
(`clip_ct=(0.0, 80.0)` default, `clip_ct=None` disables it; MR inputs are
|
||||
never clipped); `antspyx` / `elastix` / `simpleitk` register the CT as-is.
|
||||
Re-checked on the **unclipped score scale** (the protocol used by the table
|
||||
above): synthmorph keeps a genuine gain (0.241→0.277 median, 7/9 cases),
|
||||
niftyreg is neutral (0.296→0.292; its apparent +0.024 above is mostly the
|
||||
clip-scoring effect). All rows of the main table are scored on the
|
||||
unclipped warped CT, so pipelines remain comparable.
|
||||
|
||||
## Reproducing
|
||||
|
||||
```bash
|
||||
# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (pip/conda-installed here)
|
||||
# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1)
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg \
|
||||
--n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
|
||||
# force the conda 1.5.69 CPU build: --niftyreg-backend conda
|
||||
|
||||
# Docker pipelines (need the pulled image; see the SynthMorph / IMPACT
|
||||
# sections below for the one-time setup)
|
||||
python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \
|
||||
--workers 4 --threads 24
|
||||
# NOTE: impact now runs on the M6 pairs (the parameter map relaxes the
|
||||
# sample-validity gate to RequiredRatioOfValidSamples 0.1) but is by far the
|
||||
# slowest: ~5 min/case at 500/2000, ~42 min/case (2534 s median) at 2000/5000
|
||||
# (--impact-iters / --impact-samples) — see the IMPACT section
|
||||
python scripts/bench_registration.py --pipelines impact --n-cases 10 \
|
||||
--workers 4 --threads 24
|
||||
|
||||
# A/B: does clipping the CT to a brain window help? (--clip-ct clips the CT
|
||||
# for ALL pipelines, before registration and scoring; niftyreg/synthmorph
|
||||
# already clip by default, so this only affects antspyx/elastix/simpleitk)
|
||||
python scripts/bench_registration.py --n-cases 10 --workers 4 --threads 24 \
|
||||
--out-dir results/bench_reg_noclip
|
||||
python scripts/bench_registration.py --n-cases 10 --workers 4 --threads 24 \
|
||||
--clip-ct 0 80 --out-dir results/bench_reg_clip0_80
|
||||
python scripts/bench_registration.py --merge --out-dir results/bench_reg_noclip
|
||||
python scripts/bench_registration.py --merge --out-dir results/bench_reg_clip0_80
|
||||
|
||||
# combine the per-pipeline JSONs accumulated across envs into one table + plot
|
||||
python scripts/bench_registration.py --merge
|
||||
```
|
||||
|
|
@ -249,11 +354,9 @@ available pipelines with:
|
|||
```bash
|
||||
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
|
||||
python scripts/bench_registration.py \
|
||||
--pipelines simpleitk elastix antspyx --n-cases 10 --qa-images
|
||||
--pipelines simpleitk elastix antspyx niftyreg --n-cases 10 --qa-images
|
||||
python scripts/bench_registration.py \
|
||||
--pipelines synthmorph --n-cases 10 --qa-images
|
||||
python scripts/bench_registration.py \
|
||||
--pipelines impact --n-cases 10 --qa-images
|
||||
```
|
||||
|
||||
(9 images each; the incomplete-CT case is skipped.)
|
||||
|
|
@ -299,20 +402,33 @@ No extra model assets: the rigid weights ship inside the image
|
|||
## 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.
|
||||
suite (KCL-BMEIS, the actively-maintained NiftyReg): **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.
|
||||
|
||||
**Default backend: the official v2.1.1 release in Docker.** The v2.x release
|
||||
binaries need GLIBC_2.34, which this host (Ubuntu 20.04, glibc 2.31) lacks, so
|
||||
the v2.1.1 Ubuntu zip is packaged in the `niftyreg:2.1.1-ubuntu` image
|
||||
(nvidia/cuda ubuntu22.04 base for the CUDA plugin). GPU runs use
|
||||
`docker run --gpus all` + `-platf 1 -gpuid N` and are ~20× faster than CPU;
|
||||
the v2.1.1 CUDA path supports **linear interpolation only** (`-interp 1`), so
|
||||
GPU runs use linear and CPU/conda runs use cubic (`-interp 3`) — the poses
|
||||
they recover on these cases agree to ~4 decimal places. The conda-forge
|
||||
build (1.5.69, CPU) remains the fallback when the Docker image is absent
|
||||
(`--niftyreg-backend conda`).
|
||||
|
||||
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)
|
||||
# default: build the v2.1.1 Docker image (one-time)
|
||||
docker build -t niftyreg:2.1.1-ubuntu niftyreg/ # from NiftyReg-Ubuntu-v2.1.1.zip
|
||||
|
||||
# fallback: conda-forge build (provides reg_aladin + reg_f3d, NiftyReg 1.5.69)
|
||||
conda install -c conda-forge niftyreg
|
||||
```
|
||||
|
||||
|
|
@ -322,8 +438,10 @@ conda install -c conda-forge niftyreg
|
|||
(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).
|
||||
coarse), `-interp 1` (GPU) / `-interp 3` (CPU) internal warping, default
|
||||
initialisation (image centres). GPU: `--gpus all`, `-platf 1`,
|
||||
`-gpuid <auto by pid>`; CPU/conda: `-platf 0` (v2.1.1) or no `-platf`
|
||||
(1.5.69).
|
||||
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
|
||||
|
|
@ -336,12 +454,17 @@ conda install -c conda-forge niftyreg
|
|||
The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
|
||||
image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121)
|
||||
*semantic* similarity metric compiled in as a native Elastix metric, plus a
|
||||
TorchScript **TotalSegmentator** feature extractor. It is GPU-based and
|
||||
Docker-based (alongside `synthmorph`). It is also the one pipeline that keeps
|
||||
**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").
|
||||
TorchScript segmentation feature extractor (**MRSegmentator** `MRSeg/MRSeg.pt`
|
||||
by default; TotalSegmentator `TS/M730.pt` and other models from
|
||||
`VBoussot/impact-torchscript-models` are selectable). It is GPU-based and
|
||||
Docker-based (alongside `synthmorph`). It follows the shared rules (full
|
||||
volumes, no mask, fixed = CT, moving = T1c). The metric only evaluates samples
|
||||
that map inside the moving buffer, and the CT head-neck FOV (~300×300×251 mm)
|
||||
is much larger than the T1c head crop (~220×220×150 mm), so ~14% of the samples
|
||||
are in-buffer at the auto-init. The image build's sample-validity gate
|
||||
defaults to `RequiredRatioOfValidSamples 0.25` and used to abort the M6 pairs
|
||||
on that ("too many samples map outside moving image buffer"); the parameter
|
||||
map now sets it to `0.1` (see the IMPACT specifics below).
|
||||
|
||||
Prerequisites (check with `RP.availability()`; a missing piece raises
|
||||
`PipelineUnavailable` with the exact command):
|
||||
|
|
@ -354,25 +477,40 @@ docker --version
|
|||
docker pull vboussot/elastix_impact
|
||||
|
||||
# 3. Fetch the TorchScript model (~123 MB, one-time) into the repo's model
|
||||
# cache (data/impact_models/). M730.pt is the recommended TotalSegmentator.
|
||||
mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/TS
|
||||
curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/TS/M730.pt \
|
||||
https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/TS/M730.pt
|
||||
# cache (data/impact_models/). MRSeg.pt (MRSegmentator) is the default:
|
||||
# it is trained on both MRI and CT, the right prior for this
|
||||
# cross-modality CT->T1c task (all models share the same 2^l+3 FOV
|
||||
# architecture, so the parameter map is unchanged between them).
|
||||
mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg
|
||||
curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg/MRSeg.pt \
|
||||
https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/MRSeg/MRSeg.pt
|
||||
```
|
||||
|
||||
Other models (M258/M291/M850, MIND, SAM2.1, …) are selectable with
|
||||
`--impact-model <HF-relative-path-or-abs>`; see the
|
||||
Other models (TotalSegmentator `TS/M730.pt`, `TS/M258/M291/M850`, MIND, SAM2.1,
|
||||
Anatomix, …) are selectable with `--impact-model <HF-relative-path-or-abs>`;
|
||||
see the
|
||||
[impact-torchscript-models](https://huggingface.co/VBoussot/impact-torchscript-models)
|
||||
repo for the full list and FOV/layer tables.
|
||||
repo for the full list and FOV/layer tables. The model choice does not affect
|
||||
the sample-validity gate (below): the "too many samples outside moving buffer"
|
||||
check is geometric (CT head-neck FOV vs the tight T1c crop) and fires at the
|
||||
first sample evaluation, before feature extraction.
|
||||
|
||||
### How the adapter drives the image
|
||||
|
||||
The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image
|
||||
binary directly (`/lib/elastix-install/bin/elastix`):
|
||||
|
||||
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.
|
||||
1. **Full volumes, fixed=CT, moving=T1c** — both volumes are used as-is (no
|
||||
head box, no mask, no crop / pad, per the shared rules above). The metric
|
||||
only evaluates spatial samples that map inside the moving buffer, and the
|
||||
CT FOV (~300×300×251 mm, head+neck) is much larger than the tight T1c head
|
||||
crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
|
||||
the samples are in-buffer. The image build's default validity gate
|
||||
(`RequiredRatioOfValidSamples 0.25`) aborted on that ("Too many samples map
|
||||
outside moving image buffer: 704 / 5000"); the parameter map therefore sets
|
||||
`(RequiredRatioOfValidSamples 0.1)` so the M6 pairs run on their in-buffer
|
||||
samples while a genuinely diverging transform (<10% in-buffer) still
|
||||
aborts.
|
||||
2. **Auto init (no `-t0`)** — the parameter map sets
|
||||
`AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so
|
||||
elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the
|
||||
|
|
@ -381,9 +519,24 @@ binary directly (`/lib/elastix-install/bin/elastix`):
|
|||
**Jacobian** mode (differentiable, per the recommended config), 4 semantic
|
||||
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
|
||||
ASGD, `Random` sampler, `(RequiredRatioOfValidSamples 0.1)`. A copy of the
|
||||
exact parameter map is dumped to
|
||||
`data/m6reg/impact/paramdump/ParameterMap_<model>_<pid>.txt` before every
|
||||
docker run (also recorded as `param_dump` in the result JSON). GPU
|
||||
auto-spread across the host GPUs by pid (`--impact-gpu -1` forces CPU).
|
||||
4. **Intensity standardization happens inside the metric** (no pipeline-side
|
||||
normalization — the model card's "standardize intensities to zero mean and
|
||||
unit variance" is already applied): at init the metric computes
|
||||
whole-volume `[min, max, mean, sigma]` for each image
|
||||
(`itk::StatisticsImageFilter`) and passes it to the model, whose
|
||||
`normalize` branch rescales every extracted patch with the image's own
|
||||
`(mean, sigma)`. It is computed over the whole buffer (CT air included) and
|
||||
identically for fixed and moving, so both modalities' features sit on the
|
||||
same normalized scale. Pre-standardizing the MHA files in the pipeline
|
||||
would be an exact no-op (z-scoring is affine-invariant: the metric
|
||||
recomputes the stats of the image it receives) and would only corrupt the
|
||||
HU scale used by MI scoring and the QA brain-window images.
|
||||
5. **Read-back** — the result `TransformParameters.0.itk.txt` is an
|
||||
`Euler3DTransform` in the world physical frame (the
|
||||
`AutomaticTransformInitialization` is folded in by elastix); the adapter
|
||||
runs it through the usual `fix_orientation` gate.
|
||||
|
|
@ -391,13 +544,14 @@ binary directly (`/lib/elastix-install/bin/elastix`):
|
|||
## Notes / caveats
|
||||
|
||||
- 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.
|
||||
register on the full volumes (no head box, no mask, no crop / pad) 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 MI (median 0.238; 0.06–0.10 on the repositioned cases) is a
|
||||
|
|
@ -439,11 +593,16 @@ binary directly (`/lib/elastix-install/bin/elastix`):
|
|||
- **`-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.
|
||||
- **GPU is the default now** (v2.1.1, `-platf 1` via
|
||||
`docker run --gpus all`): ~20–25 s/case on the full volumes vs ~4 min
|
||||
CPU (120–440 s, `-omp` thread-bound) for the conda 1.5.69 build.
|
||||
`--niftyreg-gpu -1` forces the CPU platform (`-platf 0`).
|
||||
- **Version note.** The upstream KCL-BMEIS release binaries (v2.x) need
|
||||
glibc ≥ 2.34; this host (Ubuntu 20.04, glibc 2.31) runs them in the
|
||||
`niftyreg:2.1.1-ubuntu` Docker image (nvidia/cuda ubuntu22.04 base;
|
||||
the CUDA plugin needs libcusolver.so.11, provided by the image). The
|
||||
conda-forge 1.5.69 build remains as the no-Docker fallback and gives
|
||||
essentially identical poses (MI matches to ±0.001 on the 10 cases).
|
||||
- **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
|
||||
|
|
@ -470,28 +629,41 @@ binary directly (`/lib/elastix-install/bin/elastix`):
|
|||
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. 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: 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).
|
||||
- **Slow, by design.** Every ASGD iteration runs the segmentation
|
||||
extractor (MRSegmentator by default) 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: measured 2534 s/case (median) on the no-mask M6
|
||||
pairs, ~276 s/case at 500/2000. `--impact-iters` / `--impact-samples`
|
||||
trade accuracy for speed.
|
||||
- **Runs on M6 with a relaxed sample-validity gate.** Like every other
|
||||
pipeline it runs fixed=CT/moving=T1c on the full volumes (no mask, no
|
||||
crop / pad). The CT FOV (~300×300×251 mm) is much larger than the T1c
|
||||
crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of
|
||||
the samples map inside the T1c moving buffer. The metric calls
|
||||
`CheckNumberOfSamples()` after every evaluation and the image build
|
||||
defaults to `RequiredRatioOfValidSamples 0.25`, so it used to abort
|
||||
("Too many samples map outside moving image buffer: 704 / 5000"). The
|
||||
parameter map now sets `(RequiredRatioOfValidSamples 0.1)`: the gate
|
||||
still aborts a genuinely diverging transform (<10% in-buffer) while the
|
||||
FOV-mismatched M6 pairs run on their ~14% in-buffer samples — exactly
|
||||
the native "only in-buffer samples are evaluated" behaviour. In the
|
||||
2000/5000 re-run the gate still fired once, mid-registration: on
|
||||
`m6_26TYQI3C` the transform drifted until only 468/5000 samples (9.4%)
|
||||
stayed in-buffer at iteration 256 and the run aborted. The transform
|
||||
contract (T1c→CT) holds; `fix_orientation` decides the direction as
|
||||
usual (all 8 completed runs report `auto_inverted: True`).
|
||||
- **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: 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.
|
||||
- **History: single-init convergence.** The semantic landscape has more
|
||||
local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer
|
||||
nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0;
|
||||
the pre-mask head-box flow got it to mid 8.3 / std 13.2). The no-mask
|
||||
MRSegmentator runs show the same pattern: strong MI (0.635 on
|
||||
`m6_22IRRKU6` at 2000/5000) but the pose still misses the concentricity
|
||||
acceptance rule there (mid 26.1 / std 28.6); it *is* accepted on
|
||||
`m6_2BBN5XA4` (mid 15.1 / std 10.7).
|
||||
- **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
|
||||
|
|
|
|||
35
niftyreg/Dockerfile
Normal file
35
niftyreg/Dockerfile
Normal file
|
|
@ -0,0 +1,35 @@
|
|||
# NiftyReg v2.1.1 (official KCL-BMEIS Ubuntu release) in Docker.
|
||||
#
|
||||
# The host (Ubuntu 20.04, glibc 2.31) is too old for the v2.x binaries, which
|
||||
# need GLIBC_2.34; the nvidia/cuda ubuntu22.04 base (glibc 2.35) runs them
|
||||
# directly. The CUDA runtime provides libcusolver.so.11, which the
|
||||
# libniftyreg_cuda.so plugin (dlopen'd by reg_aladin with -platf 1) links
|
||||
# against; libcuda.so.1 comes from the host driver via `docker run --gpus all`.
|
||||
# Build:
|
||||
# docker build -t niftyreg:2.1.1-ubuntu niftyreg/
|
||||
# Run (GPU, the pipeline's default):
|
||||
# docker run --rm --gpus all -v /path/to/work:/work -w /work \
|
||||
# niftyreg:2.1.1-ubuntu -ref ref.nii.gz -flo flo.nii.gz \
|
||||
# -rigOnly -aff aff.txt -platf 1
|
||||
# Run (CPU fallback):
|
||||
# docker run --rm -v /path/to/work:/work -w /work \
|
||||
# niftyreg:2.1.1-ubuntu -ref ref.nii.gz -flo flo.nii.gz \
|
||||
# -rigOnly -aff aff.txt -platf 0
|
||||
# reg_aladin is the ENTRYPOINT, so the NiftyReg CLI flags follow directly.
|
||||
# -platf 0|1 = CPU|CUDA, -gpuid selects the device, -omp sets OpenMP threads.
|
||||
FROM nvidia/cuda:12.4.1-runtime-ubuntu22.04
|
||||
|
||||
RUN apt-get update && apt-get install -y --no-install-recommends \
|
||||
libgomp1 libstdc++6 zlib1g \
|
||||
&& rm -rf /var/lib/apt/lists/*
|
||||
|
||||
COPY NiftyReg-Ubuntu-v2.1.1/ /niftyreg/
|
||||
# reg_aladin searches /usr/local/lib for the CUDA plugin
|
||||
RUN chmod +x /niftyreg/* \
|
||||
&& cp /niftyreg/libniftyreg_cuda.so /usr/local/lib/ \
|
||||
&& ldconfig
|
||||
|
||||
ENV PATH=/niftyreg:$PATH \
|
||||
LD_LIBRARY_PATH=/niftyreg:/usr/local/lib
|
||||
|
||||
ENTRYPOINT ["reg_aladin"]
|
||||
|
|
@ -1,8 +1,9 @@
|
|||
"""Build M6-2025 manifests: T1c volumes, labeled when GTV exists.
|
||||
|
||||
GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on a head box
|
||||
located in CT space; multi-start initial tilts because inter-session head tilt
|
||||
reaches ~30 deg and tilted MI local minima can outscore the true alignment).
|
||||
GTV lives in CT space; T1c is registered to CT (rigid, Mattes MI on the full
|
||||
volumes -- no mask, no crop / pad; multi-start initial tilts because
|
||||
inter-session head tilt reaches ~30 deg and tilted MI local minima can
|
||||
outscore the true alignment).
|
||||
The inverse transform then warps CT and GTV (nearest) into T1c native space.
|
||||
Acceptance requires the warped CT skull ring to be concentric with the MR head
|
||||
silhouette (crop-independent check).
|
||||
|
|
@ -105,6 +106,10 @@ def _skull_ring_gap(warped_ct, ha):
|
|||
gap between the MR head boundary and the CT skull ring (p97 of positive
|
||||
CT) at two levels. Returns {"mid_mm", "low_mm", "std_mm"} or None when no
|
||||
skull ring is found. Aligned = small gap with low angular variance.
|
||||
|
||||
The ring test is ``>=`` (not ``>``) so it also works on a CT clipped to a
|
||||
brain window (e.g. [0, 80] HU), where the p97 threshold sits exactly on
|
||||
the clipped-skull atom and a strict comparison would give an empty ring.
|
||||
"""
|
||||
posw = warped_ct[warped_ct > 0]
|
||||
if not posw.size:
|
||||
|
|
@ -119,7 +124,7 @@ def _skull_ring_gap(warped_ct, ha):
|
|||
sl = sil[iz]
|
||||
if sl.sum() < 500:
|
||||
return None, None
|
||||
sk = warped_ct[iz] > bthr
|
||||
sk = warped_ct[iz] >= bthr
|
||||
if not sk.any():
|
||||
return None, None
|
||||
hy, hx = np.nonzero(sl)
|
||||
|
|
@ -171,14 +176,11 @@ def _valid_pyramid(size, shrink, smooth, min_dim=6):
|
|||
|
||||
|
||||
def register_ct_to_t1c(ct, t1c, threads=12):
|
||||
"""Rigid registration with CT as fixed and T1c as moving, Mattes MI on
|
||||
head-box images.
|
||||
"""Rigid registration with CT as fixed and T1c as moving, Mattes MI on the
|
||||
full volumes (no mask, no crop / pad).
|
||||
|
||||
Both inputs share physical space (oriented). The registration grid is a head
|
||||
box located in CT space from the top of the head (the CTs are head-neck
|
||||
exams, so an MR-head crop would show mostly neck/air on the CT); the box
|
||||
excludes most air so MI is computed on tissue. The box is in physical
|
||||
space, so the resulting transform applies directly to the full volumes.
|
||||
Both inputs share physical space (oriented). The transform is computed on
|
||||
the full CT / T1c volumes, so it applies directly to any resampling.
|
||||
Returns (transform, info): transform is the registration transform, mapping
|
||||
T1c (moving) coordinates -> CT (fixed) coordinates; resampling CT/TV onto
|
||||
the T1c grid with it applies the inverse transform to CT-space points,
|
||||
|
|
@ -198,38 +200,9 @@ def register_ct_to_t1c(ct, t1c, threads=12):
|
|||
from src.common import head_mask_from_image
|
||||
ha = sitk.GetArrayFromImage(t1c)
|
||||
hm = head_mask_from_image(ha)
|
||||
# Registration grid = a head box located in CT space. The CTs are head-neck
|
||||
# (often head-neck-chest) exams, so cropping to the MR head box would leave
|
||||
# mostly neck/chest/air on the fixed image and MI locks onto the wrong
|
||||
# structures. Vertex = topmost tissue slice with 10mm of continuity below
|
||||
# (skips small artifacts floating above the head; bone HU is unreliable
|
||||
# here, so tissue > -150 HU is used). x/y from the tissue extent of the
|
||||
# 150mm band below the vertex. Box covers the whole head (x/y: extent +
|
||||
# 15mm margin, z: +10 to -215 from vertex).
|
||||
ca = sitk.GetArrayFromImage(ct)
|
||||
tissue = ca > -150
|
||||
csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
|
||||
vz = _ct_vertex_slice(ct)
|
||||
if vz is not 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)
|
||||
else: # no head found: fall back to the whole CT volume
|
||||
gs = ct.GetSize()
|
||||
box_lo = tuple(float(x) for x in co)
|
||||
box_hi = (co[0] + gs[0] * csp[0], co[1] + gs[1] * csp[1], co[2] + gs[2] * csp[2])
|
||||
# full volumes, no mask / crop / pad: fixed = CT, moving = T1c
|
||||
fixed, moving = ct, t1c
|
||||
tsp = np.array(t1c.GetSpacing())
|
||||
shape = tuple(int(np.ceil((box_hi[d] - box_lo[d]) / tsp[d])) + 1 for d 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 box_lo))
|
||||
fixed = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
|
||||
moving = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)
|
||||
|
||||
stages = []
|
||||
|
||||
|
|
@ -261,8 +234,21 @@ def register_ct_to_t1c(ct, t1c, threads=12):
|
|||
# rank candidates. Run the full pyramid from the centered start; if its
|
||||
# full-resolution geometric concentricity (skull-ring gap) is poor, repeat
|
||||
# from a grid of initial tilts and keep the best by gap.
|
||||
# Euler rotation center = the CT head centroid (the 150 mm band below the
|
||||
# vertex, same head definition as before); off maps it onto the MR head
|
||||
# centroid. Both are init-only quantities; the registration itself runs on
|
||||
# the full volumes.
|
||||
ca = sitk.GetArrayFromImage(ct)
|
||||
csp, co = np.array(ct.GetSpacing()), np.array(ct.GetOrigin())
|
||||
tissue = ca > -150
|
||||
vz = _ct_vertex_slice(ct)
|
||||
if vz is not None:
|
||||
z0 = max(0, vz - int(round(150 / csp[2])))
|
||||
idx = np.array([z0, 0, 0]) + np.argwhere(tissue[z0:vz + 1]).mean(0)
|
||||
else:
|
||||
idx = np.argwhere(tissue).mean(0)
|
||||
ct_c = np.array(co) + idx * csp
|
||||
mr_c = np.array(t1c.GetOrigin()) + np.argwhere(hm).mean(0) * tsp
|
||||
ct_c = np.array([(box_lo[d] + box_hi[d]) / 2.0 for d in range(3)])
|
||||
off = (float(mr_c[0] - ct_c[0]), float(mr_c[1] - ct_c[1]), float(mr_c[2] - ct_c[2]))
|
||||
rad = np.pi / 180.0
|
||||
|
||||
|
|
@ -276,7 +262,9 @@ def register_ct_to_t1c(ct, t1c, threads=12):
|
|||
return tr0, gap
|
||||
|
||||
def is_good(g):
|
||||
return bool(g) and g["std_mm"] <= 12.0 and 0.0 <= g["mid_mm"] <= 20.0
|
||||
# |mid| (not 0 <= mid): full-volume registration can land the CT skull
|
||||
# ring slightly inside the MR silhouette, which is equally concentric
|
||||
return bool(g) and g["std_mm"] <= 12.0 and abs(g["mid_mm"]) <= 20.0
|
||||
|
||||
tr_reg, gap = full_run(0, 0, 0)
|
||||
if not is_good(gap):
|
||||
|
|
@ -325,7 +313,6 @@ def register_ct_to_t1c(ct, t1c, threads=12):
|
|||
"offset": [round(float(x), 3) for x in t0], "scale": scale},
|
||||
"coverage_head_box": round(coverage, 3),
|
||||
"skull_ring_gap_mm": skull_gap,
|
||||
"ct_head_box_mm": [[round(float(x), 1) for x in box_lo], [round(float(x), 1) for x in box_hi]],
|
||||
"seconds_total": round(time.time() - t_start, 1)}
|
||||
print("reg done %.0fs (rigid %.1fs) mi %.3f -> %.3f cov %.2f%s" % (
|
||||
info["seconds_total"], sum(s["seconds"] for s in stages),
|
||||
|
|
@ -373,7 +360,7 @@ def reg_worker(row):
|
|||
tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1))
|
||||
tv_in_ct = tv_in_ct_frac >= 0.5 # slab must straddle the lesion
|
||||
gap = info["skull_ring_gap_mm"]
|
||||
if gap and 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
|
||||
if gap and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0:
|
||||
align = "concentric"
|
||||
else:
|
||||
align = "not_concentric"
|
||||
|
|
|
|||
|
|
@ -20,6 +20,10 @@ higher = better in-sync = better aligned. The production geometric QA
|
|||
- accepted : production rule (nvox>0 and tv_in_ct and concentric).
|
||||
- reg_wall_s : wall-clock registration time (excludes I/O + scoring).
|
||||
|
||||
`--clip-ct LO HI` (e.g. `--clip-ct 0 80`, brain window) clips the CT to
|
||||
[LO, HI] before registration AND before scoring, so a clipped vs unclipped
|
||||
A/B compares the same input to both the toolkit and the MI score.
|
||||
|
||||
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
|
||||
|
|
@ -93,6 +97,14 @@ def _m6():
|
|||
return RP._load_m6()
|
||||
|
||||
|
||||
def _clip_image(img, lo, hi):
|
||||
"""Clip image intensities to [lo, hi] (geometry preserved)."""
|
||||
a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
|
||||
out = sitk.GetImageFromArray(a)
|
||||
out.CopyInformation(img)
|
||||
return out
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# case selection
|
||||
# --------------------------------------------------------------------------- #
|
||||
|
|
@ -152,7 +164,7 @@ def evaluate(ct, t1c, tv, T):
|
|||
if nvox > 0:
|
||||
tvd = ndimage.binary_dilation(warped_tv > 0.5, iterations=3)
|
||||
tv_in_ct_frac = float((tvd & tissue_w).sum() / max(tvd.sum(), 1))
|
||||
ok_gap = bool(gap) and 0.0 <= gap["mid_mm"] <= 20.0 and (gap["std_mm"] or 0) <= 12.0
|
||||
ok_gap = bool(gap) and abs(gap["mid_mm"]) <= 20.0 and (gap["std_mm"] or 0) <= 12.0
|
||||
accepted = nvox > 0 and tv_in_ct_frac >= 0.5 and ok_gap
|
||||
return {
|
||||
"skull_ring_gap_mm": gap,
|
||||
|
|
@ -237,7 +249,9 @@ def _call(pipeline, ct, t1c, tv, args):
|
|||
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)
|
||||
return RP.register_niftyreg(ct, t1c, threads=args.threads,
|
||||
backend=args.niftyreg_backend,
|
||||
gpu=args.niftyreg_gpu)
|
||||
raise KeyError(pipeline)
|
||||
|
||||
|
||||
|
|
@ -255,6 +269,11 @@ def _run_case(pipeline, case, args):
|
|||
"alignment": "ct_incomplete", "accepted": False,
|
||||
"note": "CT head incomplete (data issue, not a pipeline failure)"})
|
||||
return out
|
||||
if args.clip_ct is not None:
|
||||
# clip BEFORE registration AND scoring so the A/B compares like for like
|
||||
lo, hi = args.clip_ct
|
||||
ct = _clip_image(ct, lo, hi)
|
||||
out["clip_ct"] = [lo, hi]
|
||||
t0 = time.time()
|
||||
T, meta = _call(pipeline, ct, t1c, tv, args)
|
||||
out["reg_wall_s"] = round(time.time() - t0, 1)
|
||||
|
|
@ -380,8 +399,8 @@ def write_plot(rows, out_png):
|
|||
print("no available pipelines to plot")
|
||||
return
|
||||
labels = list(ok_rows)
|
||||
regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels]
|
||||
mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels]
|
||||
regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") or 0 for p in labels]
|
||||
mi = [(ok_rows[p].get("mi_ct_mr") or {}).get("median") or 0 for p in labels]
|
||||
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
|
||||
colors = ["#4878CF", "#6ACC65"]
|
||||
for ax, (title, vals, fmt, ylab) in zip(axes, [
|
||||
|
|
@ -487,6 +506,22 @@ def parse_args():
|
|||
ap.add_argument("--impact-samples", type=int, default=5000,
|
||||
help="IMPACT NumberOfSpatialSamples per iteration (default 5000, the "
|
||||
"elastix default)")
|
||||
ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"],
|
||||
default="auto",
|
||||
help="niftyreg reg_aladin build: docker = official v2.1.1 in the "
|
||||
"niftyreg:2.1.1-ubuntu image (host glibc is too old for the "
|
||||
"v2.x binaries), conda = conda-forge 1.5.69 on PATH, "
|
||||
"auto (default) = docker when the image is present else conda")
|
||||
ap.add_argument("--niftyreg-gpu", type=int, default=None,
|
||||
help="niftyreg CUDA device index (docker backend; default: "
|
||||
"auto-spread over the host GPUs by pid, -1 forces CPU "
|
||||
"-platf 0; on the GPU the run uses 'docker run --gpus all "
|
||||
"-platf 1' with linear interpolation, the v2.1.1 CUDA "
|
||||
"limit)")
|
||||
ap.add_argument("--clip-ct", nargs=2, type=float, default=None, metavar=("LO", "HI"),
|
||||
help="clip the CT (HU) to [LO, HI] BEFORE registration and scoring "
|
||||
"(e.g. --clip-ct 0 80 = brain window, skull/air removed; "
|
||||
"default: no clipping)")
|
||||
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")
|
||||
|
|
|
|||
|
|
@ -23,15 +23,26 @@ returns (``tr_reg.GetInverse()`` with ``tr_reg = Execute(fixed=CT, moving=T1c)``
|
|||
|
||||
Pipelines (listed alphabetically)
|
||||
---------------------------------
|
||||
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
|
||||
Shared rules for every pipeline:
|
||||
|
||||
- **No mask, no crop / pad**: registration always runs on the full input
|
||||
volumes exactly as given (no ``FixedMask`` / ``MovingMask``, no head-box
|
||||
resampling, no padding).
|
||||
- **CT is always the fixed image** (moving = T1c). When both inputs are MR
|
||||
(no CT present) the first argument is kept as fixed (see
|
||||
:func:`fixed_moving`).
|
||||
- **CT brain-window clip for the pipelines that benefit from it**: per the
|
||||
CT-clipping A/B (see REGISTRATION_BENCHMARK.md), ``niftyreg`` and
|
||||
``synthmorph`` clip the CT to [0, 80] HU before registration
|
||||
(``clip_ct`` parameter, default on, MR inputs never clipped); the
|
||||
Mattes-MI pipelines (antspyx / elastix / simpleitk) register the CT as-is.
|
||||
|
||||
Each pipeline uses its 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``).
|
||||
(global deep method).
|
||||
|
||||
- ``register_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
|
||||
Default ``initial_transform=None`` computes a translation aligning the image
|
||||
|
|
@ -42,9 +53,9 @@ TotalSegmentator metric -- see ``register_impact``).
|
|||
- ``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); elastix
|
||||
``AutomaticTransformInitialization`` (CenterOfGravity) provides the start.
|
||||
Needs Docker + the image + a model asset.
|
||||
GPU (Jacobian mode). full volumes, fixed=CT/moving=T1c (no mask -- shared
|
||||
rules above); elastix ``AutomaticTransformInitialization`` (CenterOfGravity)
|
||||
provides the start. Needs Docker + the image + a model asset.
|
||||
- ``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
|
||||
|
|
@ -123,6 +134,50 @@ def _warp_ct(ct, t1c, T):
|
|||
return sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
|
||||
|
||||
|
||||
def _is_ct(img):
|
||||
"""Modality heuristic: CT is in Hounsfield units and carries air at large
|
||||
negative values; MR intensities are non-negative."""
|
||||
return float(sitk.GetArrayFromImage(img).min()) < -50.0
|
||||
|
||||
|
||||
def _modality_name(img):
|
||||
return "CT" if _is_ct(img) else "MR"
|
||||
|
||||
|
||||
def fixed_moving(img_a, img_b):
|
||||
"""Shared fixed-image rule: the CT is always the fixed image. When both
|
||||
inputs are MR (or both CT) the given order is kept, i.e. ``img_a`` is the
|
||||
fixed image. Returns ``(fixed, moving)``.
|
||||
|
||||
Note: a CT that has been intensity-clipped to a non-negative window (e.g.
|
||||
[0, 80] HU) is no longer detectable as CT, so the given order is kept --
|
||||
callers should pass the CT first (as all pipelines do)."""
|
||||
if _is_ct(img_b) and not _is_ct(img_a):
|
||||
return img_b, img_a
|
||||
return img_a, img_b
|
||||
|
||||
|
||||
def _clip_intensity(img, lo, hi):
|
||||
"""Clip image intensities to [lo, hi] (geometry preserved)."""
|
||||
a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi)
|
||||
out = sitk.GetImageFromArray(a)
|
||||
out.CopyInformation(img)
|
||||
return out
|
||||
|
||||
|
||||
def _maybe_clip_ct(fixed, clip_ct):
|
||||
"""Apply the pipeline's CT brain-window clip to the fixed image.
|
||||
|
||||
The A/B benchmark (``--clip-ct 0 80``) showed clipping the CT to a brain
|
||||
window helps the block-matching / deep pipelines (niftyreg, synthmorph)
|
||||
and hurts the Mattes-MI ones (antspyx, elastix, simpleitk), so those two
|
||||
pipelines clip the CT by default. Returns ``(fixed, clipped: bool)``;
|
||||
MR fixed images (both-MR pairs) are never clipped."""
|
||||
if clip_ct and _is_ct(fixed):
|
||||
return _clip_intensity(fixed, *clip_ct), True
|
||||
return fixed, False
|
||||
|
||||
|
||||
def fix_orientation(ct, t1c, cand):
|
||||
"""Return (best_T, inverted) choosing between ``cand`` and ``cand^-1``.
|
||||
|
||||
|
|
@ -203,8 +258,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
|
|||
except AttributeError:
|
||||
pass
|
||||
ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32)
|
||||
# full volumes, fixed = CT, moving = T1c (no head box)
|
||||
fixed, moving = ct, t1c
|
||||
# full volumes, no mask / crop / pad; CT is the fixed image (or, if both
|
||||
# inputs are MR, the first argument is)
|
||||
fixed, moving = 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
|
||||
|
|
@ -236,8 +292,9 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000):
|
|||
"toolkit": f"SimpleITK {sitk.Version_VersionString()}",
|
||||
"metric": "MattesMutualInformation (full volumes)",
|
||||
"transform": (f"Rigid (Euler3D, fixed-LR budget, pyramid {sh}, "
|
||||
f"{iters_full} iters/level, identity start)"),
|
||||
f"{iters_full} iters/level, identity start)"),
|
||||
"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,
|
||||
}
|
||||
|
|
@ -315,9 +372,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
|
|||
after a successful registration.
|
||||
"""
|
||||
itk = _elastix_api()
|
||||
# full volumes, fixed = CT, moving = T1c (no head box)
|
||||
fixed = ct
|
||||
moving = t1c
|
||||
# full volumes, no mask / crop / pad; CT is the fixed image
|
||||
fixed, moving = fixed_moving(ct, t1c)
|
||||
params = {
|
||||
"ImageDimension": 3,
|
||||
"FixedInternalImagePixelType": "float",
|
||||
|
|
@ -373,7 +429,8 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
|
|||
"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),
|
||||
"fixed": _modality_name(fixed), "moving": _modality_name(moving),
|
||||
"seconds_reg": round(dt, 1),
|
||||
}
|
||||
return T, meta
|
||||
|
||||
|
|
@ -408,8 +465,8 @@ 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)
|
||||
# full volumes, fixed = CT, moving = T1c (no head box)
|
||||
fx, mv = ct, t1c
|
||||
# full volumes, no mask / crop / pad; CT is the fixed image
|
||||
fx, mv = fixed_moving(ct, t1c)
|
||||
import time
|
||||
t0 = time.time()
|
||||
# No manual centroid init: antspyx's default initial_transform=None computes
|
||||
|
|
@ -431,7 +488,8 @@ def register_antspyx(ct, t1c, transform_type="Rigid", aff_iterations=None,
|
|||
"pipeline": "antspyx", "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}",
|
||||
"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),
|
||||
"auto_inverted": inverted, "fixed": _modality_name(fx), "moving": _modality_name(mv),
|
||||
"aff_iterations": list(aff_iterations),
|
||||
"aff_shrink_factors": list(aff_shrink_factors), "seconds_reg": round(dt, 1),
|
||||
}
|
||||
return T, meta
|
||||
|
|
@ -484,7 +542,8 @@ def _lta_to_sitk(path):
|
|||
|
||||
|
||||
def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
|
||||
extent=256, keep_scratch=False):
|
||||
extent=256, keep_scratch=False,
|
||||
clip_ct=(0.0, 80.0)):
|
||||
"""SynthMorph CT->T1c (fixed=CT, moving=T1c) in the Docker image
|
||||
``freesurfer/synthmorph`` (TensorFlow ``VxmAffineFeatureDetector``,
|
||||
pretrained ``synthmorph.rigid.1.h5`` weights).
|
||||
|
|
@ -492,7 +551,10 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
|
|||
SynthMorph (Hoffmann et al., IEEE TMI 2022) is an anatomy-aware,
|
||||
acquisition-agnostic deep registration: each volume is normalized to an
|
||||
isotropic 1-mm LIA network space and the (pretrained) network estimates the
|
||||
transform, so no skull-stripping or intensity preprocessing is needed. Being
|
||||
transform, so no skull-stripping or intensity preprocessing is needed --
|
||||
except the one it benefits from: the CT is clipped to a brain window
|
||||
(``clip_ct=(0.0, 80.0)`` by default, per the CT-clipping A/B) before
|
||||
registration; pass ``clip_ct=None`` to register the CT as-is. Being
|
||||
a global deep method it recovers the gross inter-session offset on its own,
|
||||
so **no centroid init is used**. The pipeline
|
||||
- runs the ``rigid`` (6-DOF) matrix model, on GPU when one is visible
|
||||
|
|
@ -524,11 +586,15 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
|
|||
# per-case scratch on a local (non-NFS) tmp fs
|
||||
tmp = tempfile.mkdtemp(prefix="synthmorph_reg_")
|
||||
try:
|
||||
# fixed = CT, moving = T1c (shared convention for all pipelines)
|
||||
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "fixed.nii.gz"), True) # CT
|
||||
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "moving.nii.gz"), True) # T1c
|
||||
# full volumes, no mask / crop / pad; CT is the fixed image
|
||||
fixed, moving = fixed_moving(ct, t1c)
|
||||
fixed_name = _modality_name(fixed) # before clipping (a clipped CT
|
||||
# is no longer intensity-detectable)
|
||||
fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
|
||||
sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "fixed.nii.gz"), True)
|
||||
sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "moving.nii.gz"), True)
|
||||
|
||||
# run as the host uid/gid so the container-written LTA is host-owned
|
||||
# (clean scratch removal) and TF still reaches the GPU.
|
||||
|
|
@ -570,10 +636,11 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
|
|||
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, "
|
||||
f"extent {extent})"),
|
||||
"convention": "T1c->CT (verified)", "auto_inverted": inverted,
|
||||
"fixed": "CT", "moving": "T1c",
|
||||
"fixed": fixed_name, "moving": _modality_name(moving),
|
||||
"seconds_reg": round(dt, 1),
|
||||
"docker_image": SYNTHMORPH_IMAGE, "model": model, "gpu": gpu,
|
||||
"extent": extent,
|
||||
"clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
|
||||
}
|
||||
if keep_scratch:
|
||||
meta["scratch"] = tmp
|
||||
|
|
@ -587,7 +654,7 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None,
|
|||
IMPACT_IMAGE = "vboussot/elastix_impact"
|
||||
IMPACT_ELASTIX_BIN = "/lib/elastix-install/bin/elastix"
|
||||
IMPACT_MODEL_CACHE = "data/impact_models"
|
||||
IMPACT_DEFAULT_MODEL = "TS/M730.pt"
|
||||
IMPACT_DEFAULT_MODEL = "MRSeg/MRSeg.pt"
|
||||
IMPACT_HF_MODEL_URL = (
|
||||
"https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/{rel}"
|
||||
)
|
||||
|
|
@ -690,6 +757,14 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
|
|||
f'(ImpactUseMixedPrecision "{mixed}")',
|
||||
"(ImpactFeaturesMapUpdateInterval -1)",
|
||||
'(ImpactWriteFeatureMaps "false")',
|
||||
# safety gate: abort if < ratio of the spatial samples map inside the
|
||||
# moving buffer. The build default is 0.25, but the no-mask M6 pairs
|
||||
# (CT head-neck FOV vs the tight T1c crop) only have ~14% in-buffer at
|
||||
# the CenterOfGravity auto-init, so elastix aborts at the first metric
|
||||
# evaluation. 0.1 keeps the gate (a diverging transform still aborts)
|
||||
# while letting the FOV-mismatched runs proceed -- the metric natively
|
||||
# evaluates only the in-buffer samples.
|
||||
"(RequiredRatioOfValidSamples 0.1)",
|
||||
"",
|
||||
"// rigid 6-DOF to be optimized (auto center init, no -t0)",
|
||||
'(Transform "EulerTransform")',
|
||||
|
|
@ -708,7 +783,8 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
|
|||
'(FixedInternalImagePixelType "float")',
|
||||
'(MovingInternalImagePixelType "float")',
|
||||
'(UseDirectionCosines "true")',
|
||||
"",
|
||||
]
|
||||
L += [
|
||||
'(Registration "MultiMetricMultiResolutionRegistration")',
|
||||
f"(NumberOfResolutions {n})",
|
||||
f"(MaximumNumberOfIterations {iters})",
|
||||
|
|
@ -748,16 +824,25 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes,
|
|||
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
|
||||
"""IMPACT rigid CT->T1c (fixed=CT, moving=T1c) in the Docker image
|
||||
``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
|
||||
|
||||
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
|
||||
one: a TorchScript segmentation model (``MRSeg/MRSeg.pt`` by default,
|
||||
MRSegmentator -- robust multi-modality segmentation on MRI *and* CT, the
|
||||
right prior for this cross-modality CT->T1c task; TotalSegmentator
|
||||
``TS/M730.pt`` and the other models in the HF repo
|
||||
``VBoussot/impact-torchscript-models`` work too via ``model=``) is run on
|
||||
randomly sampled patches of the fixed and moving images and the
|
||||
registration minimizes the distance between the extracted feature maps
|
||||
(L2, per the recommended config). This is a genuinely different metric from
|
||||
the Mattes-MI pipelines above -- it is the cross-modality tool the image was
|
||||
built for.
|
||||
(L2, per the recommended config). Intensity standardization (zero mean /
|
||||
unit variance per the model card) is applied *inside the metric*: at init
|
||||
it computes whole-volume [min, max, mean, sigma] per image and the model's
|
||||
normalize branch rescales every patch with that image's own (mean, sigma),
|
||||
so no pipeline-side pre-normalization is needed (and would be an exact
|
||||
no-op -- z-scoring is affine-invariant). This is a genuinely different
|
||||
metric from the Mattes-MI pipelines above -- it is the cross-modality tool
|
||||
the image was built for.
|
||||
|
||||
The pipeline registers on the full CT / T1c volumes (no head box) with no
|
||||
manual centroid init: elastix's ``AutomaticTransformInitialization``
|
||||
|
|
@ -765,12 +850,17 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
|||
``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).
|
||||
Like every other pipeline this runs with **fixed=CT, moving=T1c** on the
|
||||
full volumes: no mask, no crop / pad (shared rules above). The IMPACT
|
||||
metric draws random spatial samples from the fixed (CT) buffer and
|
||||
evaluates them only where they map inside the moving (T1c) buffer; when
|
||||
too few samples are in-buffer, elastix aborts ("too many samples map
|
||||
outside moving image buffer"). On the M6 head-neck CTs (FOV
|
||||
~300x300x251 mm) vs the tight T1c head crop (~220x220x150 mm) only ~14%
|
||||
of the samples are in-buffer at the CenterOfGravity auto-init, so the
|
||||
M6 runs abort -- IMPACT needs a fixed volume whose FOV is comparable to
|
||||
the moving crop (e.g. a head-only CT). The returned transform 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'
|
||||
|
|
@ -792,13 +882,10 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
|||
raise PipelineUnavailable(f"impact: {note}")
|
||||
model_host, model_rel = _impact_model_host(model)
|
||||
|
||||
# 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
|
||||
# full volumes, no mask / crop / pad; CT is the fixed image (or, if both
|
||||
# inputs are MR, the first argument). No manual init -- elastix's
|
||||
# AutomaticTransformInitialization aligns the image centers.
|
||||
fixed, moving = fixed_moving(ct, t1c)
|
||||
|
||||
# GPU: auto-spread across the host GPUs by pid, or forced value
|
||||
if gpu is None:
|
||||
|
|
@ -806,8 +893,8 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
|||
gpu = (os.getpid() % ngpu) if ngpu else -1
|
||||
|
||||
# per-level feature-extraction scales (recommended 6/3/1.5/1 mm), scaled so
|
||||
# the finest level matches the box's (T1c) base spacing
|
||||
base = float(np.median(t1c.GetSpacing()))
|
||||
# the finest level matches the moving image's base spacing
|
||||
base = float(np.median(moving.GetSpacing()))
|
||||
n = len(scales)
|
||||
voxel_sizes = [base * s for s in scales]
|
||||
# image pyramid: IMPACT's docs recommend NO image resampling (all-1 schedules)
|
||||
|
|
@ -840,10 +927,22 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
|||
os.path.join(datadir, "Fixed_image.mha"), True)
|
||||
sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32),
|
||||
os.path.join(datadir, "Moving_image.mha"), True)
|
||||
with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
|
||||
f.write(_impact_param_map(os.path.join("/Data/Models", model_rel),
|
||||
param_text = _impact_param_map(os.path.join("/Data/Models", model_rel),
|
||||
gpu, iters, samples, voxel_sizes,
|
||||
rescale=schedule, smoothing=smoothe))
|
||||
rescale=schedule, smoothing=smoothe)
|
||||
with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f:
|
||||
f.write(param_text)
|
||||
|
||||
# debug aid: keep a persistent copy of the exact parameter map elastix
|
||||
# is about to run, written BEFORE the docker call so an abort leaves it
|
||||
# behind for inspection (per-pid name: parallel workers do not clobber
|
||||
# each other's dumps)
|
||||
dump_dir = d("data/m6reg/impact/paramdump")
|
||||
os.makedirs(dump_dir, exist_ok=True)
|
||||
family = os.path.basename(os.path.dirname(model_rel)) or "model"
|
||||
param_dump = os.path.join(dump_dir, f"ParameterMap_{family}_{os.getpid()}.txt")
|
||||
with open(param_dump, "w") as f:
|
||||
f.write(param_text)
|
||||
|
||||
cmd = ["docker", "run", "--rm"]
|
||||
if gpu >= 0:
|
||||
|
|
@ -870,26 +969,28 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
|||
raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
|
||||
|
||||
# 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.
|
||||
# (AutomaticTransformInitialization is folded in by elastix); elastix
|
||||
# maps moving(T1c) -> fixed(CT), which is the contract direction;
|
||||
# fix_orientation verifies it against the skull ring.
|
||||
T_main = sitk.ReadTransform(itk_tf)
|
||||
T, inverted = fix_orientation(ct, t1c, T_main)
|
||||
finally:
|
||||
if not keep_scratch:
|
||||
shutil.rmtree(tmp, ignore_errors=True)
|
||||
|
||||
family = os.path.basename(os.path.dirname(model_rel)) \
|
||||
or os.path.splitext(os.path.basename(model_rel))[0]
|
||||
meta = {
|
||||
"pipeline": "impact",
|
||||
"toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
|
||||
"metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)",
|
||||
"metric": f"IMPACT semantic ({family} {os.path.basename(model_rel)}, L2, Jacobian)",
|
||||
"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",
|
||||
"fixed": _modality_name(fixed), "moving": _modality_name(moving),
|
||||
"seconds_reg": round(dt, 1),
|
||||
"docker_image": IMPACT_IMAGE, "model": model_rel, "gpu": gpu,
|
||||
"iters": iters, "samples": samples,
|
||||
"iters": iters, "samples": samples, "param_dump": param_dump,
|
||||
}
|
||||
if keep_scratch:
|
||||
meta["scratch"] = tmp
|
||||
|
|
@ -901,6 +1002,9 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000,
|
|||
# --------------------------------------------------------------------------- #
|
||||
|
||||
|
||||
NIFTYREG_DOCKER_IMAGE = "niftyreg:2.1.1-ubuntu"
|
||||
|
||||
|
||||
def _niftyreg_bin():
|
||||
"""Locate reg_aladin (NiftyReg) on PATH or raise PipelineUnavailable."""
|
||||
b = shutil.which("reg_aladin")
|
||||
|
|
@ -913,6 +1017,26 @@ def _niftyreg_bin():
|
|||
"re-run in the activated env.")
|
||||
|
||||
|
||||
def _niftyreg_docker_ok():
|
||||
"""Return (ok, note): is the NiftyReg v2.1.1 Docker image available?
|
||||
|
||||
The official KCL-BMEIS v2.x release binaries need GLIBC_2.34, which this
|
||||
host (Ubuntu 20.04, glibc 2.31) lacks; the ``niftyreg:2.1.1-ubuntu`` image
|
||||
(built from the release zip, see niftyreg/Dockerfile) runs them instead.
|
||||
"""
|
||||
if shutil.which("docker") is None:
|
||||
return False, "docker CLI not found on PATH"
|
||||
try:
|
||||
r = subprocess.run(["docker", "image", "inspect", NIFTYREG_DOCKER_IMAGE],
|
||||
capture_output=True, text=True, timeout=30)
|
||||
if r.returncode != 0:
|
||||
return False, (f"docker image {NIFTYREG_DOCKER_IMAGE} not built "
|
||||
f"(run: docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/)")
|
||||
except Exception as e: # noqa
|
||||
return False, f"docker image check failed ({type(e).__name__}: {e})"
|
||||
return True, NIFTYREG_DOCKER_IMAGE
|
||||
|
||||
|
||||
def _niftyreg_version(binpath):
|
||||
try:
|
||||
r = subprocess.run([binpath, "--version"], capture_output=True,
|
||||
|
|
@ -938,7 +1062,8 @@ def _parse_niftyreg_affine(path):
|
|||
return T
|
||||
|
||||
|
||||
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)):
|
||||
"""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
|
||||
|
|
@ -951,32 +1076,96 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
|
|||
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.
|
||||
The CT is clipped to a brain window (``clip_ct=(0.0, 80.0)`` by default,
|
||||
per the CT-clipping A/B) before registration; pass ``clip_ct=None`` to
|
||||
register the CT as-is. MR inputs are never clipped.
|
||||
|
||||
Preconditions (else :class:`PipelineUnavailable`): ``reg_aladin`` on PATH
|
||||
(``conda install -c conda-forge niftyreg``).
|
||||
Full volumes (no head box, no mask / crop / pad); 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.
|
||||
|
||||
``backend`` selects the reg_aladin build: ``"docker"`` runs the official
|
||||
KCL-BMEIS **v2.1.1** release in the ``niftyreg:2.1.1-ubuntu`` Docker image
|
||||
(built from NiftyReg-Ubuntu-v2.1.1.zip, see niftyreg/Dockerfile -- the
|
||||
v2.x binaries need GLIBC_2.34, which this host lacks); ``"conda"`` uses
|
||||
the conda-forge ``reg_aladin`` on PATH (NiftyReg 1.5.69); ``"auto"``
|
||||
(default) prefers Docker when the image is present and falls back to
|
||||
conda.
|
||||
|
||||
``gpu`` selects the CUDA device for the Docker backend: ``None`` (default)
|
||||
auto-spreads over the host GPUs by pid, ``-1`` forces CPU. On the GPU the
|
||||
run uses ``docker run --gpus all`` + ``-platf 1 -gpuid N`` and
|
||||
``-interp 1`` (the v2.1.1 CUDA path supports linear interpolation only);
|
||||
on the CPU it uses ``-platf 0`` with ``-interp 3`` (cubic), as does the
|
||||
conda backend. The GPU is ~20x faster than CPU on the full volumes.
|
||||
|
||||
Preconditions (else :class:`PipelineUnavailable`): the Docker image for
|
||||
backend docker/auto, or ``reg_aladin`` on PATH (``conda install -c
|
||||
conda-forge niftyreg``).
|
||||
"""
|
||||
import time
|
||||
|
||||
binpath = _niftyreg_bin()
|
||||
# full volumes, reference = CT, floating = T1c (no head box)
|
||||
if backend not in ("auto", "conda", "docker"):
|
||||
raise PipelineUnavailable(
|
||||
f"niftyreg: unknown backend {backend!r} (use auto|conda|docker)")
|
||||
dok, dnote = _niftyreg_docker_ok()
|
||||
use_docker = backend == "docker" or (backend == "auto" and dok)
|
||||
if use_docker:
|
||||
if not dok:
|
||||
raise PipelineUnavailable(f"niftyreg: {dnote}")
|
||||
binpath = None
|
||||
if gpu is None:
|
||||
ngpu = _n_gpus()
|
||||
gpu = (os.getpid() % ngpu) if ngpu else -1
|
||||
use_gpu = gpu >= 0
|
||||
platf = "1" if use_gpu else "0"
|
||||
# the v2.1.1 CUDA path supports linear interpolation only
|
||||
interp = "1" if use_gpu else "3"
|
||||
toolkit = (f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE}, "
|
||||
f"CUDA)" if use_gpu
|
||||
else f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE})")
|
||||
else:
|
||||
binpath = _niftyreg_bin()
|
||||
gpu, use_gpu, platf = -1, False, None
|
||||
interp = "3"
|
||||
toolkit = _niftyreg_version(binpath)
|
||||
# full volumes, no mask / crop / pad; reference (fixed) = CT, floating = T1c
|
||||
fixed, moving = fixed_moving(ct, t1c)
|
||||
fixed_name = _modality_name(fixed) # before clipping (a clipped CT is no
|
||||
# longer intensity-detectable)
|
||||
fixed, ct_clipped = _maybe_clip_ct(fixed, clip_ct)
|
||||
tmp = tempfile.mkdtemp(prefix="niftyreg_reg_")
|
||||
dt = 0.0
|
||||
try:
|
||||
sitk.WriteImage(sitk.Cast(ct, sitk.sitkFloat32),
|
||||
sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32),
|
||||
os.path.join(tmp, "ref.nii.gz"), True)
|
||||
sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32),
|
||||
sitk.WriteImage(sitk.Cast(moving, 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"]
|
||||
flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz",
|
||||
"-rigOnly", "-aff", "aff.txt",
|
||||
"-ln", str(levels), "-maxit", str(maxit),
|
||||
"-interp", interp, "-omp", str(threads)]
|
||||
if platf is not None:
|
||||
flags += ["-platf", platf]
|
||||
if use_gpu:
|
||||
flags += ["-gpuid", str(gpu)]
|
||||
flags += ["-voff"]
|
||||
if use_docker:
|
||||
# reg_aladin is the image ENTRYPOINT; the work dir is /work.
|
||||
# --gpus all: the v2.1.1 Ubuntu build ships libniftyreg_cuda.so
|
||||
# (dlopen'd by reg_aladin with -platf 1); the NVIDIA container
|
||||
# toolkit provides libcuda.so.1, and -gpuid selects the device.
|
||||
cmd = ["docker", "run", "--rm"]
|
||||
if use_gpu:
|
||||
cmd += ["--gpus", "all"]
|
||||
cmd += ["-v", f"{tmp}:/work", "-w", "/work",
|
||||
NIFTYREG_DOCKER_IMAGE] + flags
|
||||
else:
|
||||
cmd = [binpath] + flags
|
||||
t0 = time.time()
|
||||
proc = subprocess.run(cmd, capture_output=True, text=True,
|
||||
cwd=tmp, timeout=7200)
|
||||
|
|
@ -999,14 +1188,16 @@ def register_niftyreg(ct, t1c, threads=24, levels=3, maxit=15, keep_scratch=Fals
|
|||
|
||||
meta = {
|
||||
"pipeline": "niftyreg",
|
||||
"toolkit": _niftyreg_version(binpath),
|
||||
"toolkit": toolkit,
|
||||
"backend": "docker" if use_docker else "conda",
|
||||
"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",
|
||||
"fixed": fixed_name, "moving": _modality_name(moving),
|
||||
"seconds_reg": round(dt, 1),
|
||||
"levels": levels, "maxit": maxit,
|
||||
"levels": levels, "maxit": maxit, "interp": interp, "gpu": gpu,
|
||||
"clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None,
|
||||
}
|
||||
if keep_scratch:
|
||||
meta["scratch"] = tmp
|
||||
|
|
@ -1047,11 +1238,20 @@ def availability():
|
|||
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))
|
||||
dok, dnote = _niftyreg_docker_ok()
|
||||
if dok:
|
||||
ver = f"NiftyReg 2.1.1 (reg_aladin, Docker {NIFTYREG_DOCKER_IMAGE})"
|
||||
if b:
|
||||
ver += f" | conda fallback: {_niftyreg_version(b)}"
|
||||
out[name] = (True, ver)
|
||||
elif b:
|
||||
out[name] = (True, _niftyreg_version(b)
|
||||
+ " | docker v2.1.1 image missing (build: "
|
||||
f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/)")
|
||||
else:
|
||||
out[name] = (False, "reg_aladin not on PATH "
|
||||
"(conda install -c conda-forge niftyreg)")
|
||||
out[name] = (False, "reg_aladin not on PATH and no docker image; "
|
||||
"conda install -c conda-forge niftyreg or "
|
||||
f"docker build -t {NIFTYREG_DOCKER_IMAGE} niftyreg/")
|
||||
elif name == "synthmorph":
|
||||
ok, note = _synthmorph_docker_ok()
|
||||
if not ok:
|
||||
|
|
|
|||
Loading…
Reference in a new issue