From 16b9828d15bfc31d45fe14c19e23c0495fb76893 Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Tue, 29 Sep 2026 17:31:24 +0800 Subject: [PATCH] Enhance registration pipelines with CT clipping and NiftyReg Docker support - Added `--clip-ct LO HI` argument to clip CT intensities before registration and scoring, improving A/B comparisons. - Implemented `_clip_image` function to handle intensity clipping of images. - Updated registration functions to utilize CT clipping where beneficial, specifically in `register_niftyreg` and `register_synthmorph`. - Enhanced `fixed_moving` function to ensure CT is always treated as the fixed image. - Introduced Docker support for NiftyReg v2.1.1, including a Dockerfile for building the image. - Updated documentation and help messages to reflect new features and usage instructions. - Adjusted evaluation criteria for registration quality to account for intensity clipping. --- .gitignore | 4 +- REGISTRATION_BENCHMARK.md | 349 ++++++++++++++++++++++----------- niftyreg/Dockerfile | 35 ++++ scripts/02_build_m6_dataset.py | 79 ++++---- scripts/bench_registration.py | 43 +++- src/registration_pipelines.py | 315 ++++++++++++++++++++++------- 6 files changed, 587 insertions(+), 238 deletions(-) create mode 100644 niftyreg/Dockerfile diff --git a/.gitignore b/.gitignore index 4664ead..ba2b644 100644 --- a/.gitignore +++ b/.gitignore @@ -8,4 +8,6 @@ runs_nnu/ runs_monai/ __pycache__/ -elastix/ \ No newline at end of file +elastix/ +niftyreg/NiftyReg-Ubuntu-v2.1.1.zip +niftyreg/NiftyReg-Ubuntu-v2.1.1/ diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index 7f16871..6974684 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/REGISTRATION_BENCHMARK.md @@ -22,36 +22,47 @@ 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 and `impact` aborts on the M6 +pairs ("too many samples map outside moving image buffer", ~14% of samples +in-buffer at the auto-init) — 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** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ⛔ aborts on M6 (no-mask rule; 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 +73,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 +119,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 | 0 | - | - | - (aborts, no mask) | +| **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 | abort | 0.292 | 0.642 | 0.277 | +| m6_23XIMO4M_20230222 | **0.655** | 0.416 | abort | 0.632 | 0.652 | 0.616 | +| m6_24V62QQM_20230223 | **0.704** | 0.701 | abort | 0.260 | 0.663 | 0.256 | +| m6_24V62QQM_20231017 | 0.601 | 0.590 | abort | 0.493 | **0.602** | 0.484 | +| m6_24V62QQM_20240820 | **0.653** | 0.644 | abort | 0.118 | 0.652 | 0.111 | +| m6_26TYQI3C_20221027 | 0.659 | 0.653 | abort | 0.609 | **0.695** | 0.590 | +| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | abort | 0.069 | **0.577** | 0.061 | +| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | abort | 0.233 | 0.634 | 0.231 | +| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | abort | 0.340 | 0.651 | 0.323 | ### Interpretation @@ -163,52 +179,128 @@ 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) does not run on these pairs under the no-mask + rule** — earlier (masked) runs showed it strongest on the hardest cases + (`m6_22IRRKU6` 0.624, `m6_2C4TKI6H` 0.626, `m6_2CDNN7BH` 0.618, median + 0.511) with a stochastic draw per run; it also was by far the slowest + (~14 min/case at the 500/2000 budget, hours/case at 2000/5000). With the CT + head-neck FOV as the fixed volume it aborts ("too many samples map outside + moving image buffer", ~10–20% in-buffer) — see the IMPACT section. - **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 does not run on + M6 under the no-mask rule; when it did — masked, 500/2000 budget — it was + ~833 s/case, GPU feature extraction with backprop per iteration.) - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its per-case numbers carry run-to-run variance (this table is one draw); SimpleITK's MI (and registration) is deterministic at fixed thread count; 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 | - | - (still aborts; clipping does not fix the FOV geometry, 194–388/2000 samples in-buffer) | + +**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 aborts on the M6 pairs under the shared no-mask rule (the CT +# head-neck FOV maps mostly outside the T1c head crop) — 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,14 +341,13 @@ 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.) +(9 images each; the incomplete-CT case is skipped. `impact` produces none on +M6 — it aborts under the no-mask rule, see the IMPACT section.) ## SynthMorph pipeline (Docker) — one-time setup @@ -299,20 +390,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 +426,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 `; 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 @@ -337,11 +443,13 @@ The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.co 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"). +Docker-based (alongside `synthmorph`). It follows the shared rules (full +volumes, no mask, fixed = CT, moving = T1c) — which means it **aborts on the +M6 pairs**: the IMPACT 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 and elastix raises "too many samples map outside moving image +buffer" (see the IMPACT specifics below). Prerequisites (check with `RP.availability()`; a missing piece raises `PipelineUnavailable` with the exact command): @@ -370,9 +478,15 @@ repo for the full list and FOV/layer tables. 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). **On the M6 + pairs this aborts**: the IMPACT 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 2000 samples are in-buffer and + elastix raises "Too many samples map outside moving image buffer: 277 / + 2000". IMPACT therefore needs a fixed volume whose FOV is comparable to + the moving crop (e.g. a head-only CT) to work under the no-mask rule. 2. **Auto init (no `-t0`)** — the parameter map sets `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the @@ -391,13 +505,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 +554,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 @@ -478,20 +598,21 @@ binary directly (`/lib/elastix-install/bin/elastix`): 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). +- **Aborts on M6 under the no-mask rule.** 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 most of the CT buffer maps outside the T1c moving + buffer and the metric aborts ("Too many samples map outside moving image + buffer: 277 / 2000" on M6 at the CenterOfGravity auto-init). The + transform contract (T1c→CT) still holds when a run completes; for M6 + pairs a head-only CT as the fixed volume would be needed. - **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). Under the + current no-mask rule it does not run on M6 at all (see above). - **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 diff --git a/niftyreg/Dockerfile b/niftyreg/Dockerfile new file mode 100644 index 0000000..5716b33 --- /dev/null +++ b/niftyreg/Dockerfile @@ -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"] \ No newline at end of file diff --git a/scripts/02_build_m6_dataset.py b/scripts/02_build_m6_dataset.py index ecd752c..9ab3675 100644 --- a/scripts/02_build_m6_dataset.py +++ b/scripts/02_build_m6_dataset.py @@ -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" diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index aec7571..c84e6e8 100644 --- a/scripts/bench_registration.py +++ b/scripts/bench_registration.py @@ -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") diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index 26f23fa..dcbb2e2 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -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 @@ -708,7 +775,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,7 +816,7 @@ 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* @@ -765,12 +833,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 +865,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 +876,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) @@ -870,9 +940,9 @@ 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: @@ -886,7 +956,7 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000, "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, @@ -901,6 +971,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 +986,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 +1031,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 +1045,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 +1157,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 +1207,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: