From 77653e4018d4ce39680cc56227050aae4e989dde Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Wed, 30 Sep 2026 21:23:18 +0800 Subject: [PATCH] refactor(impact): replace KonfAI pip runtime with Docker elastix_impact image Migrate register_impact from impact-reg-konfai (pip + one-time glibc/LibTorch shim) to the vboussot/elastix_impact Docker image running a rigid 6-DOF EulerTransform with MIND + TotalSegmentator static features and a Mattes-MI secondary metric. Delete scripts/setup_impact_konfai.sh. Vendor ImpactLoss (vboussot/ImpactLoss) into src/impact_loss.py and add impact_loss as a per-case semantic scoring metric in the benchmark runner. Update docs (AGENTS.md, REGISTRATION_BENCHMARK.md) to reflect the all-rigid pipeline set and the new Docker-based IMPACT runtime. --- AGENTS.md | 2 +- REGISTRATION_BENCHMARK.md | 542 +++++++++++++++--------------- scripts/bench_registration.py | 132 +++++--- scripts/setup_impact_konfai.sh | 110 ------- src/impact_loss.py | 152 +++++++++ src/registration_pipelines.py | 580 ++++++++++++++++----------------- 6 files changed, 786 insertions(+), 732 deletions(-) delete mode 100644 scripts/setup_impact_konfai.sh create mode 100644 src/impact_loss.py diff --git a/AGENTS.md b/AGENTS.md index cebb331..117a6e1 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -10,7 +10,7 @@ Activate with: source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 ``` -Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), `impact-reg-konfai` (pip; IMPACT — on this host the elastix binary needs the one-time shim from `scripts/setup_impact_konfai.sh`), plus the Docker image `freesurfer/synthmorph`. +Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), plus the Docker images `freesurfer/synthmorph` and `vboussot/elastix_impact` (IMPACT; rigid MIND + TotalSegmentator config, feature models under `data/impact_models`). The benchmark's ImpactLoss scoring metric uses `src/impact_loss.py` (vendored from vboussot/ImpactLoss; runs on the env's torch). ## Project diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index f8555b8..495aa99 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/REGISTRATION_BENCHMARK.md @@ -1,8 +1,7 @@ # Registration Pipeline Benchmark — CT → T1c (M6) -Benchmarks **six cross-modality registration pipelines** (five rigid; -`impact` is the native deformable B-spline preset of the ImpactReg package) -that align a head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be +Benchmarks **six cross-modality registration pipelines — all rigid** (6-DOF) +— that align a head-and-neck CT to a tightly-cropped T1c MRI so the GTV can be warped onto the T1c (MR) native grid. This is the registration step of `scripts/02_build_m6_dataset.py`. The pairs are hard: different modality (CT Hounsfield vs MR intensity) and the patient is repositioned between the CT and @@ -30,11 +29,10 @@ every pipeline: exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no padding). - **The T1c (MR) is always the fixed image** (moving = CT) for every - pipeline. `impact`'s ImpactReg **MR_CT** convention coincides with this - rule: its output displacement field (defined on the fixed grid) is the - contract transform T1c→CT directly — see the IMPACT section. (A CT clipped - to a non-negative window reads as MR, so the T1c slot decides; both-CT - pairs keep the given order.) + pipeline. (A CT clipped to a non-negative window reads as MR, so the T1c + slot decides; both-CT pairs keep the given order.) +- **Every transform is rigid (6-DOF)** — Euler / rigid-matrix parameterizations + per toolkit. No affine, no deformable. All pipelines also use each toolkit's **native automatic transform initialization** (no manual centroid offset): elastix / impact use @@ -42,13 +40,13 @@ initialization** (no manual centroid offset): elastix / impact use 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. With the ImpactReg MR_CT direction the samples -are drawn from the tight T1c head crop (fixed, ~220×220×150 mm) and map mostly -inside the large CT head-neck FOV (moving, ~300×300×251 mm), so the gate -rarely fires; the pipeline still sets `RequiredRatioOfValidSamples 0.1` (vs -the build default 0.25) as a safety override for the FOV mismatch — see the -IMPACT section. The comparison is on each toolkit's native optimizer / metric -/ global-search capability. +map inside the moving buffer. With fixed=T1c the samples are drawn from the +tight T1c head crop (fixed, ~220×220×150 mm) and map mostly inside the large +CT head-neck FOV (moving, ~300×300×251 mm), so the gate rarely fires; the +pipeline still sets `RequiredRatioOfValidSamples 0.1` (vs the build default +0.25) as a safety override for the FOV mismatch — see the IMPACT section. +The comparison is on each toolkit's native optimizer / metric / global-search +capability. ## The six pipelines @@ -58,16 +56,15 @@ 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 = T1c, moving = CT** (shared rule above; -`impact`'s ImpactReg MR_CT convention coincides with it). `niftyreg` and -`synthmorph` clip the CT (the moving image) to [0, 80] HU before +method). **All use fixed = T1c, moving = CT** (shared rule above). `niftyreg` +and `synthmorph` clip the CT (the moving image) 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` | T1c / CT | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | | 2 | `elastix` | T1c / CT | 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** (MR_CT) | auto CenterOfGravity (override) | **IMPACT semantic** (MIND `R1D2_3D` L1 + TotalSegmentator `M730` deepest-layer Dice, per-level weighted, Static) + Mattes MI (w 0.4) + bending energy | B-spline DVF (preset native), ASGD, 4 scales | ✅ ran (impact-reg-konfai + GPU; `Elastix_IMPACT_Static` preset — see IMPACT section) | +| 3 | `impact` | T1c / CT | auto CenterOfGravity | **IMPACT semantic** (MIND `R1D2_3D` L1 + TotalSegmentator `M852`/`M850` deepest-layer Dice, per-level weighted, Static) + Mattes MI (w 0.4) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (Docker `vboussot/elastix_impact` + GPU — see IMPACT section) | | 4 | `niftyreg` | T1c / CT | 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` | T1c / CT | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | | 6 | `synthmorph` | T1c / CT | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) | @@ -75,16 +72,14 @@ registration (per the CT-clipping A/B below; `clip_ct=None` disables it). Benchmark budget (best-results settings): `elastix` and `simpleitk` run with **MaximumNumberOfIterations = 2000** (per resolution / per pyramid level), with the sampler and metric left at their defaults (**NumberOfSpatialSamples -5000**, **NumberOfHistogramBins 32**). `impact` runs the published ImpactReg -preset **`Elastix_IMPACT_Static`** (the `impact-reg-konfai` package) — an -elastix B-spline (4 scales, FinalGridSpacing 14 mm, per-level -**MaximumNumberOfIterations 400/300/200/200**, **NumberOfSpatialSamples -2000**) driven by the IMPACT metric: a per-level **MIND (L1) + -TotalSegmentator M730 (deepest layer, Dice)** combination, layer weights -0.2/0.8 → 0.7/0.3, Static mode, + Mattes MI weight 0.4 + bending energy — -with CenterOfGravity auto-init and `RequiredRatioOfValidSamples 0.1` applied -as raw-elastix overrides for the M6 pairs. `antspyx`, `niftyreg` and -`synthmorph` keep their native defaults. +5000**, **NumberOfHistogramBins 32**). `impact` runs the MIND + TotalSegmentator +**static** feature combination of the ImpactReg `MR_CT_TS` reference parameter +map (per `VBoussot/ImpactReg`) as a **rigid** `EulerTransform`: 4 levels, +per-level voxel sizes 6/3/2/2 mm (no image-resampling pyramid), per-level +**MaximumNumberOfIterations 400/300/200/200**, **NumberOfSpatialSamples 2000**, ++ Mattes MI weight 0.4, CenterOfGravity auto-init, +`RequiredRatioOfValidSamples 0.1`. `antspyx`, `niftyreg` and `synthmorph` keep +their native defaults. `simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid, identity start — the default, comparable to the other tools) and @@ -108,20 +103,30 @@ m6_2BBN5XA4_20250411 m6_2C4TKI6H_20240506 m6_2CDNN7BH_20240102 ## Metrics (identical for every pipeline) -The headline alignment score is **`mi_ct_mr`** — the registration's *achieved* -Mattes MI, read off as a standalone score: +Alignment is scored two ways, both as the registration's *achieved* metric read +off as a standalone score: the intensity-based **Mattes MI** (`mi_ct_mr`) and +the semantic **ImpactLoss** (`impact_loss`): - **`mi_ct_mr`** — mean Mattes MI (nats) between the **CT inverse-transformed into MR coordinates** (the warped CT on the T1c grid) and the **original MR**, restricted to the MR head and each image percentile-binned independently (SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are better in-sync = better aligned; it collapses toward 0 when the CT is - rotated/translated away from the MR. Tables and plots report this (plus - registration time and GTV-in-CT coverage). -- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT - head tissue (the slab must straddle the lesion). + rotated/translated away from the MR. +- **`impact_loss`** — the IMPACT semantic loss from + [vboussot/ImpactLoss](https://github.com/vboussot/ImpactLoss) (the metric + family the `impact` pipeline registers with): a weighted L1 distance between + the deep feature maps of the **warped CT** and the **original MR** extracted + by the pretrained TotalSegmentator `TS/M730.pt` TorchScript model (two + deepest layers, weights 1/1, no resampling), over the **full co-registered + volumes**. The model standardizes intensity internally (per-image z-score), + so both volumes go in raw and the score is invariant to affine intensity + scaling. Lower = the two volumes carry more in-sync anatomy = better + aligned; it rises when the CT is rotated/translated away from the MR. - **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring). +Tables and plots report `mi_ct_mr` and `impact_loss` (plus registration time). + The production geometric QA is still computed per case and kept in the JSON (it drives dataset acceptance in `02_build_m6_dataset.py`, but its summary numbers were found to correlate poorly with perceived alignment, so they are no @@ -132,111 +137,116 @@ longer reported in the tables/plots): CT) at mid and low levels, plus the angular std. `mid` measures translational concentricity; `std` measures rotational concentricity. - **`tv_nvox`** — warped GTV volume. -- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the - production acceptance rule; aligned ⇔ `|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). +- **`accepted`** — the production acceptance rule of `02_build_m6_dataset.py` + (warped GTV non-empty and concentric with the head; 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, no clip) Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over -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 = T1c / moving = CT (shared rule; `impact`'s MR_CT convention -coincides with it). `niftyreg` and `synthmorph` clip the CT (the moving -image) 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` runs the `impact-reg-konfai` package on the GPU (see the IMPACT -section). +the MR head; higher = better) **and** the **ImpactLoss** semantic score +(warped-CT ~ MR over the full volumes; lower = better), with the +**unclipped** CT, so all rows are comparable on the same score scales. Full +volumes, no mask / crop / pad, fixed = T1c / moving = CT (shared rule). +`niftyreg` and `synthmorph` clip the CT (the moving image) to [0, 80] HU +*before registration* (per the CT-clipping A/B below; `clip_ct=None` +disables it); `antspyx`, `elastix`, `impact` and `simpleitk` register the CT +as-is. `niftyreg` is the v2.1.1 Docker build on the GPU (`-platf 1`); +`impact` runs the `vboussot/elastix_impact` Docker image on the GPU (see the +IMPACT section). All six transforms are rigid (6-DOF). -| pipeline | toolkit | registered | s/reg (med) | tv-in-CT | MI (med/mean, nats) | +| pipeline | toolkit | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) | |---|---|---|---|---|---| -| **antspyx** | ANTsPy 0.63 | 9 | 15.5 | 0.93 | **0.652 / 0.650** | -| **elastix** | ITK-Elastix 0.25.4 | 9 | 9.8 | 0.93 | **0.646 / 0.638** | -| **impact** | impact-reg-konfai (Elastix 5.2 + IMPACT, `Elastix_IMPACT_Static` preset, MIND+M730, deformable) | 9 | 216.4 | 0.93 | **0.578 / 0.558** | -| **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 16.2 | 0.93 | **0.292 / 0.339** | -| **simpleitk** | SimpleITK 2.5.6 | 9 | 9.6 | 0.93 | **0.647 / 0.642** | -| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 27.9 | 0.93 | **0.277 / 0.328** | +| **antspyx** | ANTsPy 0.63 | 9 | 16.7 | 0.157 / 0.155 | **0.652 / 0.650** | +| **elastix** | ITK-Elastix 0.25.4 | 9 | 11.7 | 0.157 / 0.156 | 0.646 / 0.638 | +| **impact** | Docker `vboussot/elastix_impact` (Elastix 5.2 + IMPACT, rigid MIND+TS static config) | 9 | 111.2 | 0.157 / 0.156 | **0.652 / 0.641** | +| **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 15.4 | 0.166 / 0.170 | 0.292 / 0.338 | +| **simpleitk** | SimpleITK 2.5.6 | 9 | 10.8 | **0.156 / 0.156** | 0.647 / 0.642 | +| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 22.2 | 0.167 / 0.170 | 0.277 / 0.328 | Per-case MI (nats, unclipped score; **bold** = best for the case). All rows are re-runs under the **fixed = T1c / moving = CT** convention (the shared -rule); the CT-fixed-era numbers are in the git history. The `impact` column -is the `impact-reg-konfai` deformable `Elastix_IMPACT_Static` preset (the -retired rigid Docker config is archived in +rule); the CT-fixed-era numbers and the retired konfai-preset run are in the +git history / archived JSONs. The `impact` column is the Docker rigid +MIND+TS static config (the CT-fixed-era Docker run is archived in `results/bench_registration/impact_docker_legacy.json`). | 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.648** | 0.644 | 0.542 | 0.292 | 0.647 | 0.277 | -| m6_23XIMO4M_20230222 | **0.654** | 0.652 | 0.543 | 0.632 | **0.654** | 0.616 | -| m6_24V62QQM_20230223 | 0.703 | 0.698 | 0.667 | 0.260 | **0.704** | 0.256 | -| m6_24V62QQM_20231017 | **0.601** | 0.590 | 0.587 | 0.493 | 0.597 | 0.484 | -| m6_24V62QQM_20240820 | **0.652** | 0.646 | 0.612 | 0.118 | 0.638 | 0.111 | -| m6_26TYQI3C_20221027 | **0.718** | 0.650 | 0.717 | 0.609 | 0.670 | 0.590 | -| m6_2BBN5XA4_20250411 | 0.577 | 0.577 | 0.258 | 0.069 | **0.579** | 0.061 | -| m6_2C4TKI6H_20240506 | **0.645** | 0.638 | 0.578 | 0.233 | 0.636 | 0.231 | -| m6_2CDNN7BH_20240102 | **0.655** | 0.650 | 0.523 | 0.340 | 0.652 | 0.323 | +| m6_22IRRKU6_20240115 | **0.648** | 0.644 | 0.645 | 0.292 | 0.647 | 0.277 | +| m6_23XIMO4M_20230222 | **0.655** | 0.652 | 0.654 | 0.632 | 0.654 | 0.616 | +| m6_24V62QQM_20230223 | 0.703 | 0.698 | 0.700 | 0.260 | **0.704** | 0.256 | +| m6_24V62QQM_20231017 | **0.602** | 0.590 | 0.596 | 0.493 | 0.597 | 0.484 | +| m6_24V62QQM_20240820 | **0.652** | 0.646 | 0.652 | 0.118 | 0.638 | 0.111 | +| m6_26TYQI3C_20221027 | **0.718** | 0.650 | 0.658 | 0.609 | 0.670 | 0.590 | +| m6_2BBN5XA4_20250411 | 0.577 | 0.577 | 0.572 | 0.069 | **0.579** | 0.061 | +| m6_2C4TKI6H_20240506 | **0.645** | 0.638 | 0.639 | 0.233 | 0.636 | 0.231 | +| m6_2CDNN7BH_20240102 | **0.655** | 0.650 | 0.655 | 0.340 | 0.652 | 0.323 | + +Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) 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.1507 | 0.1508 | 0.1510 | 0.1656 | **0.1504** | 0.1664 | +| m6_23XIMO4M_20230222 | 0.1565 | **0.1563** | 0.1564 | 0.1570 | 0.1564 | 0.1575 | +| m6_24V62QQM_20230223 | **0.1466** | 0.1470 | 0.1469 | 0.1664 | 0.1467 | 0.1665 | +| m6_24V62QQM_20231017 | **0.1520** | 0.1539 | 0.1534 | 0.1577 | 0.1532 | 0.1576 | +| m6_24V62QQM_20240820 | **0.1566** | 0.1571 | 0.1567 | 0.1840 | 0.1585 | 0.1842 | +| m6_26TYQI3C_20221027 | **0.1554** | 0.1620 | 0.1622 | 0.1623 | 0.1604 | 0.1620 | +| m6_2BBN5XA4_20250411 | 0.1617 | **0.1612** | 0.1615 | 0.1847 | 0.1616 | 0.1856 | +| m6_2C4TKI6H_20240506 | 0.1573 | 0.1579 | 0.1580 | 0.1754 | **0.1547** | 0.1757 | +| m6_2CDNN7BH_20240102 | 0.1591 | **0.1590** | 0.1591 | 0.1734 | 0.1587 | 0.1741 | ### Interpretation -- **ANTsPy leads, with SimpleITK and Elastix within 0.006 of it** (median MI - 0.652 / 0.647 / 0.646): on the well-positioned cases the three cluster at - 0.60–0.70 and all now survive the harder repositioned cases (0.58–0.65 on - `m6_2BBN5XA4` / `m6_22IRRKU6`). -- **Switching to fixed = T1c / moving = CT lifted Elastix sharply** - (0.517 → 0.646 median, and it is now #3 rather than #4). The Mattes-MI - pipelines sample spatial pairs from the fixed buffer and discard the ones - that map outside the moving image; with the small T1c head crop as fixed and - the large CT head+neck FOV as moving, (nearly) all fixed samples are - in-buffer, versus only a minority when the big CT was fixed. Elastix's - repositioned cases jumped (`m6_2BBN5XA4` 0.287→0.577, `m6_2C4TKI6H` - 0.401→0.638, `m6_23XIMO4M` 0.416→0.652). The other pipelines are essentially - unchanged by the swap (antspyx 0.653→0.652, simpleitk 0.651→0.647, niftyreg - 0.292→0.292, synthmorph 0.277→0.277) — they were not sample-coverage - limited. +- **The top four are a statistical tie** (median MI 0.646–0.652, median + ImpactLoss 0.156–0.157): `antspyx`, `impact`, `simpleitk` and `elastix` + all recover the gross offset on every case, including the two ~30° + repositioned ones (`m6_2BBN5XA4` 0.572–0.579, `m6_24V62QQM_20240820` + 0.638–0.652). The two alignment scores — an intensity histogram and a deep + semantic one — agree on the ranking. +- **IMPACT (Docker, rigid, T1c-fixed) is in the top tier** at median MI + 0.652 (0.6523) — best or within 0.006 of the best on 8 of 9 cases + (`m6_24V62QQM_20230223` 0.700, `m6_26TYQI3C` 0.658 where it trails + ANTsPy's 0.718). This is a large step + up from the legacy CT-fixed Docker config (median 0.467, archived in + `impact_docker_legacy.json`): the T1c-fixed direction fixes the + sample-coverage problem the same way it fixed Elastix (the metric draws + samples from the fixed T1c crop and evaluates them where they map into the + large CT FOV, instead of the other way round; previously only ~14% of the + samples were in-buffer at the auto-init). The `m6_2BBN5XA4` ~30° tilt case + that defeated every earlier IMPACT variant (0.223 CT-fixed Docker, 0.258 + konfai preset) now registers at 0.572. Cost: the slowest of the six + (~111 s/case, ~10× the intensity tools). Concentricity is comparable to + the MI tools (4/9 accepted by the production rule). +- **ImpactLoss corroborates the MI leaderboard from a different metric + family.** The top four cluster within 0.001 of each other on the median + (0.1564–0.1571), and on every case the best (lowest) ImpactLoss belongs to + one of them (antspyx 4/9, elastix 3/9, simpleitk 2/9, impact 0/9); + `niftyreg` (0.1664) and `synthmorph` (0.1665) are clearly higher, and the + gap widens on the hard cases (0.184–0.186 for the two tail tools vs + 0.155–0.162 for the top four). A semantic score computed by pretrained + segmentation features (not the registration's own objective) ranking the + pipelines the same way the intensity score does is a strong check that the + leaderboard is not an artifact of the MI scale. - **The grossly repositioned cases still separate the leaders from the tail.** On `m6_2BBN5XA4` / `m6_24V62QQM_20240820` (~30° inter-session tilt) - ANTs/SimpleITK/Elastix hold 0.58–0.65, while niftyreg (0.07–0.12) and - synthmorph (0.06–0.11) — the single-start block-matching and pretrained-MRI + the top four hold 0.57–0.65, while niftyreg (0.069 / 0.118) and synthmorph + (0.061 / 0.111) — the single-start block-matching and pretrained-MRI methods — lock onto local minima. This is exactly why production `02_build_m6_dataset.py` uses a **multi-start tilt grid**. -- **IMPACT now runs the published `Elastix_IMPACT_Static` preset** - (`impact-reg-konfai` package; the old Docker backend is retired): per-level - **MIND + TotalSegmentator M730** combination (L1 + deepest-layer Dice, - weights 0.2/0.8 → 0.7/0.3, Static mode, whole-image features), the preset's - native deformable B-spline (400/300/200/200 iterations, 2000 samples, - + Mattes MI 0.4 + bending energy), registered in the MR_CT direction - (fixed=T1c) with CenterOfGravity auto-init and the 0.1 sample-validity gate - as overrides. It is the best of the IMPACT variants on M6 and now **second - overall at median MI 0.578** (behind SimpleITK/ANTsPy 0.65; ahead of - elastix 0.52), a clear step up from the old rigid Docker config (median - 0.467, archived in `impact_docker_legacy.json`): the deformable refinement - + MR_CT direction lifts the repositioned cases hardest (`m6_2C4TKI6H` - 0.334→0.578, `m6_23XIMO4M` 0.362→0.543, `m6_22IRRKU6` 0.467→0.542, - `m6_2CDNN7BH` 0.332→0.523; `m6_26TYQI3C` 0.717 — now the best tool on that - case) at a small cost on two well-aligned ones (`m6_24V62QQM_20230223` - 0.700→0.667, `m6_24V62QQM_20240820` 0.646→0.612, still accepted). - Concentricity also improved (median mid 12.0 / std 21.9 mm vs 21.1 / 28.0; - 3/9 accepted). It is the slowest of the six (~216 s/case, the deformable - optimization + KonfAI wrapper; the old rigid config was ~78 s). The - `m6_2BBN5XA4` ~30° tilt case is still its weakest (0.258) — the single - CenterOfGravity start remains the limiting factor. Earlier Docker-era - variants: MIND+MRSegmentator median 0.417 (~72 s), single-model MIND - median 0.406 (~37 s), MRSegmentator-only Jacobian configs 0.631 - (500/2000) and 0.591 (2000/5000) at ~5×/~35× the cost (archived: - `results/bench_registration/impact_marseg_*.json`, `impact_2000_5000.json`, - `impact_500_2000_explicit_smoothing.json`, `impact_mind_500_ref.json`). -- **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching - is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class; +- **NiftyReg (reg_aladin) is second-from-bottom**: single-start symmetric + block matching is strong on moderate cases (`m6_23XIMO4M` 0.632; `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). + on the GPU (`-platf 1`) is ~15 s/case (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.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 @@ -244,30 +254,25 @@ retired rigid Docker config is archived in 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, T1c-fixed re-run): - SimpleITK 9.6 s, elastix 9.8 s, **ANTs 15.5 s** (down from 36.6 s — the - smaller fixed image also shaves ANTs' per-level work), **NiftyReg 16.2 s** - (v2.1.1, GPU block matching), SynthMorph 27.9 s (GPU forward pass), - **IMPACT 216.4 s** (`Elastix_IMPACT_Static` preset — Static feature mode, no - backprop, but the deformable B-spline optimization at 1100 total iterations - + the KonfAI wrapper; the retired rigid Docker config was 78.0 s). The - earlier MRSegmentator Jacobian configs were far slower still: ~354 s/case at - 500/2000, 2534 s/case at 2000/5000 (GPU feature extraction with backprop - per iteration). +- **Speed** (median wall, 3–4 cases in parallel, 24 threads, T1c-fixed): + SimpleITK 10.8 s, elastix 11.7 s, NiftyReg 15.4 s (v2.1.1, GPU block + matching), ANTs 16.7 s, SynthMorph 22.2 s (GPU forward pass), **IMPACT + 111.2 s** (Docker, static feature mode, no backprop, 400/300/200/200 + rigid ASGD iterations). - **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 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 in the Jacobian preset - (`Elastix_IMPACT_Jacobian`) is the most stochastic**: it draws random - patches through the network every iteration, so the final pose changes run - to run — treat such numbers as one sample from a distribution. The default - Static preset re-samples coordinates every iteration as well, but the - feature maps are deterministic given the pose. + ~±0.01–0.05 MI on the weak cases). IMPACT in the **Static** config + re-samples spatial coordinates every iteration, but the feature maps are + deterministic given the pose. The earlier Jacobian-mode Docker experiments + (archived in `results/bench_registration/impact_*.json`) resampled patches + through the network each iteration, so their per-case numbers are one + sample from a distribution — e.g. the MRSegmentator-only Jacobian configs + reached 0.631 (500/2000) and 0.591 (2000/5000) at ~5×/~35× the Static + cost, and the MIND+TS/MRSeg static variants 0.417/0.406 (~72/~37 s). ## CT clipping A/B (`--clip-ct 0 80`) @@ -328,15 +333,13 @@ not clip.** ```bash # 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 \ +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 -# impact (pip; GPU; ~120-220 s/case; one-time elastix shim on hosts too old -# for the stock binary — see the IMPACT section). --impact-preset selects -# another published preset (e.g. Elastix_IMPACT_Jacobian, Generic_Rigid); -# --impact-iters / --impact-samples are no-ops for the default preset's -# per-resolution matrix +# impact (Docker vboussot/elastix_impact; GPU; ~110-120 s/case; see the IMPACT +# section). --impact-iters / --impact-samples override the reference map's +# per-level budget (one iters value for all levels) python scripts/bench_registration.py --pipelines impact --n-cases 10 \ --workers 3 --threads 24 @@ -363,10 +366,10 @@ Other flags: `--sample random --seed 0` (reproducible random 10), `--threads T` (also the OpenMP thread count for `reg_aladin`), `--simpleitk-mode {benchmark,production}`, `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu `, -`--synthmorph-extent `, `--impact-preset ` (default -`Elastix_IMPACT_Static`), `--impact-gpu `, `--impact-iters N` -(preset max_iterations override; no-op for matrix presets), -`--impact-samples N` (preset spatial_samples override), +`--synthmorph-extent `, `--impact-gpu `, `--impact-iters N` +(one MaximumNumberOfIterations value for all levels; default the reference +map's 400/300/200/200), `--impact-samples N` (NumberOfSpatialSamples; default +2000), `--out-dir `, `--qa-images` / `--qa-out ` (visual QA, below). ## Visual QA images (checkerboard + color overlay) @@ -395,9 +398,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 niftyreg --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 + --pipelines "synthmorph impact" --n-cases 10 --qa-images ``` (9 images each; the incomplete-CT case is skipped.) @@ -494,96 +497,77 @@ conda install -c conda-forge niftyreg inverts it on the M6 cases (8/9 `auto_inverted: True`) to the contract transform (T1c→CT). -## IMPACT pipeline (impact-reg-konfai) — one-time setup +## IMPACT pipeline (Docker `vboussot/elastix_impact`) — one-time setup The `impact` pipeline runs the -[impact-reg-konfai](https://pypi.org/project/impact-reg-konfai/) package: -[IMPACT](https://arxiv.org/abs/2503.24121) registration presets from the -[`VBoussot/ImpactReg`](https://huggingface.co/VBoussot/ImpactReg) Hugging Face -repo, driven through the [KonfAI](https://github.com/fideus-labs/KonfAI) -runtime — no Docker. The default preset is **`Elastix_IMPACT_Static`**: -Elastix 5.2 with the IMPACT *semantic* similarity metric — a per-level -combination of the **MIND** handcrafted descriptor (`MIND/R1D2_3D.pt`, L1, -subset 32) and **TotalSegmentator** `TS/M730.pt` (deepest layer, Dice, subset -64), layer weights 0.2/0.8 → 0.7/0.3 from coarse to fine, **Static** mode -(whole-image features, no backprop), a 4-scale B-spline (FinalGridSpacing -14 mm), 400/300/200/200 iterations per level, 2000 spatial samples, + Mattes -MI (w 0.4) + bending energy — the published successor of the `MR_CT_TS` -reference parameter map this repo used to hand-build (and replace the -`vboussot/elastix_impact` Docker image). Other presets are selectable with -`--impact-preset` (`Elastix_IMPACT_Jacobian`, `Generic_Rigid`, -`Generic_Rigid_BSpline`, `FireANTs_Anatomix/IMPACT/SyN`, -`ConvexAdam_Coarse/Composite`; the live database is listed by -`impact-reg-konfai register --help`). It is GPU-based (like `synthmorph`), -but pip-installed. +[`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact) +image: Elastix 5.2 with the [IMPACT](https://arxiv.org/abs/2503.24121) +*semantic* similarity metric ([vboussot/ImpactLoss](https://github.com/vboussot/ImpactLoss)) +as a native Elastix metric. The parameter map is the MIND + TotalSegmentator +**static** feature combination of the `MR_CT_TS` reference map in the +[`VBoussot/ImpactReg`](https://huggingface.co/VBoussot/ImpactReg) repo, with a +**rigid 6-DOF `EulerTransform`** in place of the reference's RecursiveBSpline +(its bending-energy penalty is zero for a rigid transform, so it is dropped) — +i.e. the reference metric/budget, rigid-only, on the GPU. Prerequisites (check with `RP.availability()`; a missing piece raises `PipelineUnavailable` with the exact command): ```bash -# in the longitudinal-313 env (Python >= 3.11): -pip install impact-reg-konfai # CLI + KonfAI runtime (pulls h5py etc.) -bash scripts/setup_impact_konfai.sh # one-time, this host only: the elastix shim +# 1. Docker (with NVIDIA container toolkit so --gpus works) +docker --version + +# 2. Pull the image (~38 GB, one-time) +docker pull vboussot/elastix_impact + +# 3. Feature models under the cache (the pipeline prints the exact +# curl commands if any are missing): +# data/impact_models/MIND/R1D2_3D.pt (~15 KB) +# data/impact_models/TS/M852.pt (~200 MB) +# data/impact_models/TS/M850.pt (~200 MB) +# (from HF VBoussot/impact-torchscript-models, or `sh +# Data/Models/download_models.sh` from the ImpactLoss repo) ``` -The CLI then auto-downloads the one-time assets on first `register`: the -preset app (HF `VBoussot/ImpactReg`), the TorchScript feature models (HF -`VBoussot/impact-torchscript-models`: MIND ~15 KB + TS/M730 ~200 MB) and the -**elastix-impact 1.0.0** binary (GitHub `vboussot/ImpactElastix`, ~40 MB) -into `~/.cache` (HF cache + `~/.cache/konfai`). +### How the adapter drives the image -### Why the elastix shim (this host only) - -The stock binary needs **glibc ≥ 2.35** and links **LibTorch 2.8.0**; this -host (Ubuntu 20.04, glibc 2.31, env torch 2.14) satisfies neither, so -`scripts/setup_impact_konfai.sh` builds a self-contained shim at -`~/.cache/konfai/elastix-impact-shim`: the glibc 2.35 loader + libs extracted -from the `niftyreg:2.1.1-ubuntu` image (docker is needed once for that copy -only), the self-contained LibTorch 2.8.0 cu128 from download.pytorch.org -(~3.8 GB zip), and a wrapper `bin/elastix` that runs the stock binary through -the shim loader with the LibTorch 2.8.0 library path (stripping the inherited -`LD_LIBRARY_PATH`, so the env's torch 2.14 cannot leak in). The pipeline -selects the shim via `KONFAI_ELASTIX_DIR` (`_impact_elastix_root`); a host -new enough to run the stock binary directly needs no shim at all. - -### How the adapter drives the CLI - -1. **Full volumes, fixed=T1c (MR), moving=CT** — the ImpactReg **MR_CT** - convention, which coincides with the shared rule (the T1c is the fixed - image): the preset's output is a displacement field on the **fixed** grid - mapping fixed→moving, so the field is the contract transform T1c→CT - **directly** (displacement fields cannot be inverted with `GetInverse`, - which a fixed=CT direction would require). It also keeps - the metric's sample-validity gate healthy: samples drawn from the tight - T1c head crop (~220×220×150 mm) map mostly inside the large CT head-neck - FOV (~300×300×251 mm), unlike fixed=CT where only ~14% were in-buffer at - the auto-init. -2. **M6 robustness overrides** — the preset's map is tuned for well-aligned - pairs, so the adapter forwards raw elastix lines through the package's - `parameter_overrides` knob: `AutomaticTransformInitialization=true` + - `AutomaticTransformInitializationMethod=CenterOfGravity` (the preset's - GeometricalCenter start leaves the intensity centroids ~50-100 mm apart on - the head+neck FOV mismatch) and `RequiredRatioOfValidSamples 0.1` (build - default 0.25; keeps the divergence gate while letting the FOV-mismatched - runs proceed). -3. **`--fields-only`** — the adapter derives the warp itself (the package's - derived "Moved" image would be an extra full-size resample of the same - content). -4. **GPU auto-spread** across the host GPUs by pid (`--impact-gpu -1` forces - CPU; the engine then turns mixed precision off). -5. **Read-back** — the transform comes from `/P000/Transform.h5` (an - ITK-form HDF5 `DisplacementFieldTransform` on the T1c grid; the SimpleITK - build here has no ITK HDF5 transform reader, so it is parsed with h5py) - and direction-checked by `_orient_impact`: the skull-ring gap of the field - vs a sign-flipped (first-order reverse) proxy — the `fix_orientation` - analogue for a field with no `GetInverse`. -6. **Intensity normalization happens inside the metric / model** (no +1. **Full volumes, fixed=T1c (MR), moving=CT** — the shared fixed-image rule: + the metric draws its random spatial samples from the fixed buffer and + evaluates them where they map inside the moving buffer. With the tight T1c + head crop (~220×220×150 mm) as fixed, the samples map mostly inside the + large CT head-neck FOV (~300×300×251 mm), so the + `RequiredRatioOfValidSamples 0.1` gate (build default 0.25) rarely fires. +2. **Parameter map** (the rigid `MR_CT_TS` reference config): + `(Metric "Impact" "AdvancedMattesMutualInformation")` with weights + 1.0 / 0.4; per level (coarse → fine) the **MIND** handcrafted descriptor + (`MIND/R1D2_3D.pt`, L1, subset 32, layer mask `"1"`) combined with + **TotalSegmentator** `TS/M852.pt` (levels 0-1) / `TS/M850.pt` (levels + 2-3, Dice, subset 64, deepest layers), layer weights 0.2/0.8, 0.3/0.7, + 0.6/0.4, 0.7/0.3, at ImpactVoxelSize 6/3/2/2 mm (scaled to the moving + spacing); **Static** mode (whole-image features, no backprop); no + image-resampling pyramid (all-1 schedules — the scale is carried by + ImpactVoxelSize, per the IMPACT docs); `RandomCoordinate` sampler; ASGD + with per-level MaximumNumberOfIterations 400/300/200/200 and + NumberOfSpatialSamples 2000; `EulerTransform` + + `AutomaticTransformInitialization` (CenterOfGravity) + + `AutomaticScalesEstimation`. A copy of the exact map is dumped to + `data/m6reg/impact/paramdump/` per run for inspection. +3. **GPU auto-spread** across the host GPUs by pid (`docker run --gpus all`, + `ImpactGPU` = the worker's device; `--impact-gpu -1` forces CPU and turns + mixed precision off). +4. **Read-back** — the transform comes from + `Out/TransformParameters.0.itk.txt` (an `Euler3DTransform` in the world + physical frame; the full combination transform, with the auto-init folded + in by elastix) and is direction-checked by `fix_orientation` against the + skull ring (with fixed=T1c / moving=CT the gate keeps the file's transform, + `auto_inverted: False`). +5. **Intensity normalization happens inside the metric / models** (no pipeline-side normalization): the MIND descriptor normalizes each 5-voxel neighbourhood to [0, 1] locally (modality-agnostic by construction); the - segmentation models (TS/MRSeg) are standardized to zero mean / unit - variance per the model card. Pre-normalizing the volumes would be an - exact no-op for the z-scored models and would only corrupt the HU scale - used by MI scoring and the QA brain-window images. + segmentation models are standardized to zero mean / unit variance per the + model card. Pre-normalizing the volumes would be an exact no-op for the + z-scored models and would only corrupt the HU scale used by MI scoring and + the QA brain-window images. ## Notes / caveats @@ -598,22 +582,20 @@ new enough to run the stock binary directly needs no shim at all. 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 - genuine cross-modality / tilt weakness, not variance. -- `fix_orientation` (elastix/ants/niftyreg/synthmorph) chooses between a - transform and its inverse by skull-ring concentricity — a correctness gate + IMPACT. Its low MI (median 0.277 with the default CT clip; 0.06–0.11 on the + repositioned cases) is a genuine cross-modality / tilt weakness, not + variance. +- `fix_orientation` (elastix/ants/impact/niftyreg/synthmorph) chooses between + a transform and its inverse by skull-ring concentricity — a correctness gate that makes the result independent of each toolkit's transform-direction convention, so the fixed/moving swap is absorbed by the gate rather than - re-derived per toolkit. Under the shared fixed=T1c / moving=CT rule the raw - transforms come out in the contract direction for elastix (combination - transform = fixed→moving), antspyx (ANTs forward transform) and synthmorph - (the `-T` LTA, model estimates first→second) — completed runs report - `auto_inverted: False` for those; `niftyreg`'s `-aff` file comes out the - other way and the gate inverts it (8/9 `auto_inverted: True`). `impact` - returns a displacement field (no `GetInverse` in ITK), so it uses - `_orient_impact` — the same skull-ring decision against a sign-flipped - (first-order reverse) proxy field; its MR_CT field is in the contract - direction (`auto_inverted: False`). +re-derived per toolkit. Under the shared fixed=T1c / moving=CT rule the raw + transforms come out in the contract direction for elastix (combination + transform = fixed→moving), impact (its `TransformParameters.0.itk.txt`, + the full combination Euler transform) and antspyx (ANTs forward transform) + and synthmorph (the `-T` LTA, model estimates first→second) — completed + runs report `auto_inverted: False` for those; `niftyreg`'s `-aff` file + comes out the other way and the gate inverts it (8/9 `auto_inverted: True`). - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU): - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly @@ -626,8 +608,8 @@ new enough to run the stock binary directly needs no shim at all. path). - **GPU sharing.** The container runs as the host uid with `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU. -- **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK - (~10 s) / ANTs (~37 s) and IMPACT (~216 s, deformable preset). +- **~28 s/case** (GPU: one pretrained forward pass), between SimpleITK + (~11 s) / ANTs (~17 s) and IMPACT (~115 s, rigid docker config). - **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`, symmetric block matching + TLS, CPU/OpenMP): - **Different algorithmic family.** No histogram MI objective: 4-voxel @@ -676,51 +658,39 @@ new enough to run the stock binary directly needs no shim at all. - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution (ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32). -- **IMPACT specifics** (`impact-reg-konfai` 1.8.6, KonfAI runtime, - elastix-impact 5.2, GPU; preset `Elastix_IMPACT_Static`): +- **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2 + + IMPACT metric, GPU; rigid MIND+TS static config of the `MR_CT_TS` reference + map): - **Runs on M6 with the T1c-fixed direction + relaxed sample-validity gate.** Full volumes, no mask, no crop / pad, fixed=T1c/moving=CT (the - shared rule; the ImpactReg MR_CT convention coincides with it): the - preset's output field is on the fixed grid, so that direction is the - contract transform (T1c→CT) directly, and the metric's samples (drawn - from the tight T1c crop) map mostly inside the large CT FOV. The package's generated map still defaults to - `RequiredRatioOfValidSamples 0.25`, so the adapter overrides it to 0.1 - (a genuinely diverging transform, <10% in-buffer, still aborts). In the - old Docker era (fixed=CT) the gate aborted the no-mask M6 pairs until - the same 0.1 relaxation. The transform contract holds; `_orient_impact` - decides the direction as usual (completed runs report - `auto_inverted: False` — the MR_CT field is already the contract). - - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) is - applied as an override (the preset's GeometricalCenter start is too - crude for the head+neck FOV mismatch); no manual translation is used. - - **Deformable, not rigid.** The preset's native transform is a - RecursiveBSpline (FinalGridSpacing 14 mm, 4 scales) regularized by the - bending-energy metric; the old pipeline forced a rigid `EulerTransform` - in place of it. The benchmark contract (warp the CT concentric with the - MR head) is transform-agnostic, so this is a native-config comparison - point like every other row. + shared rule): the metric's samples (drawn from the tight T1c head crop) + map mostly inside the large CT head-neck FOV, so the gate rarely fires. + The map still sets `RequiredRatioOfValidSamples 0.1` (build default + 0.25) as a safety override — a genuinely diverging transform, <10% + in-buffer, still aborts. The combination Euler transform comes out in + the contract direction (completed runs report `auto_inverted: False`). + - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) + aligns the intensity centroids (no manual translation); the reference + map's GeometricalCenter start is too crude for the head+neck FOV + mismatch. + - **Rigid, like every other pipeline.** The reference map's + RecursiveBSpline + bending-energy regularizer is replaced by a 6-DOF + `EulerTransform` (the bending penalty is zero for a rigid transform). - **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 showed the same pattern (0.635 MI on `m6_22IRRKU6` - at 2000/5000 but mid 26.1 / std 28.6 there; accepted on `m6_2BBN5XA4`, - mid 15.1 / std 10.7). The single CenterOfGravity start on the ~30° - tilt pairs remains IMPACT's limiting factor — a multi-start (or - `--tta` flips through the preset app) would be the next step. - - **Stochastic.** The Jacobian preset (`Elastix_IMPACT_Jacobian`) - resamples patches through the network each iteration, so results vary - run to run; the default Static preset re-samples coordinates each - iteration but the feature maps are deterministic given the pose (see - Stochasticity). - - **The elastix shim (this host).** The stock binary needs glibc ≥ 2.35 - + LibTorch 2.8.0; the host (glibc 2.31, torch 2.14) runs it through - `~/.cache/konfai/elastix-impact-shim` (built by - `scripts/setup_impact_konfai.sh`, selected via `KONFAI_ELASTIX_DIR`). - Same glibc story as the niftyreg v2.x binaries. - - **Models / presets are not in the repo.** They auto-download from HF - (`VBoussot/ImpactReg`, `VBoussot/impact-torchscript-models`) into - `~/.cache/huggingface` on first use; the old `data/impact_models/` - cache is obsolete. - - **CPU fallback:** `--impact-gpu -1` (no `--gpu` flag; the engine runs - the feature models on CPU with mixed precision off — much slower). \ No newline at end of file + local minima than Mattes-MI. The T1c-fixed direction fixes most of the + hard cases (the repositioned `m6_22IRRKU6` / `m6_23XIMO4M` now register + as well as the MI tools do), but the single CenterOfGravity start on the + ~30° tilt pair `m6_2BBN5XA4` remains a weakness — a multi-start tilt + grid (as production SimpleITK uses) would be the next step. + - **Stochastic.** The config is **Static** mode: it re-samples spatial + coordinates every iteration, but the feature maps are deterministic + given the pose (see Stochasticity). The older Jacobian-mode configs + (docker-era experiments, archived in `results/bench_registration/impact_*.json`) + resampled patches through the network each iteration, so their + per-case numbers are one sample from a distribution. + - **Models are under the data root, not the repo.** The TorchScript + feature models live in `data/impact_models/` (MIND ~15 KB + TS ~200 MB + each, from HF `VBoussot/impact-torchscript-models`); the pipeline + raises with the exact `curl` commands if one is missing. + - **CPU fallback:** `--impact-gpu -1` (`ImpactGPU -1`, mixed precision + off — much slower). \ No newline at end of file diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index 3a6ab8c..2d6bf6f 100644 --- a/scripts/bench_registration.py +++ b/scripts/bench_registration.py @@ -11,9 +11,10 @@ higher = better in-sync = better aligned. The production geometric QA (scripts/02_build_m6_dataset.py) is still recorded per case in the JSON: - mi_ct_mr : Mattes MI (nats), warped-CT ~ MR over the MR head - (the reported alignment score). - - tv_in_ct_frac : fraction of the (dilated) warped GTV that falls inside - the warped-CT head tissue (the slab must straddle it). + (the reported alignment score, higher = better). + - impact_loss : IMPACT semantic loss (vboussot/ImpactLoss; TotalSeg + M730 deep features, weighted L1), warped-CT ~ MR over + the full co-registered volumes (lower = better). - tv_nvox : warped GTV volume (>0.5 voxels). - skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR head silhouette (mid/low level + angular std). @@ -31,10 +32,11 @@ whichever conda env has a given toolkit installed, and it records only the pipelines that are importable there. elastix runs in the main env through the ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively (``conda install -c conda-forge niftyreg`` provides ``reg_aladin``); -``impact`` runs the ``impact-reg-konfai`` package (pip; GPU; one-time elastix -shim setup on hosts too old for the stock binary — see -REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker image -(``freesurfer/synthmorph``; needs Docker + the pulled image). The `--merge` +``impact`` runs the ``vboussot/elastix_impact`` Docker image (GPU; the rigid +MIND + TotalSegmentator static config, needs Docker + the pulled image + the +model assets — see REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker +image (``freesurfer/synthmorph``; needs Docker + the pulled image). The +`--merge` mode then combines the per-pipeline JSON files (accumulated across envs) into one summary + plot. @@ -99,6 +101,54 @@ def _m6(): return RP._load_m6() +# the ImpactLoss scoring model (the repo's default baseline for MR/CT) +IMPACT_SCORE_MODEL = "TS/M730.pt" +_IMPACT_LOSS_STATE = None + + +def _impact_loss_scorer(): + """Lazily build the ImpactLoss semantic-loss scorer (once per process) on + the GPU auto-picked for this worker. Returns ``(fn, device)`` or None if + torch / the model are unavailable (the score is then omitted).""" + global _IMPACT_LOSS_STATE + if _IMPACT_LOSS_STATE is None: + try: + import warnings + import torch + from src.impact_loss import IMPACTReg + ngpu = torch.cuda.device_count() if torch.cuda.is_available() else 0 + if ngpu: + torch.cuda.set_device(os.getpid() % ngpu) + with warnings.catch_warnings(): + # torch.jit.load is marked deprecated by upstream torch 2.14; + # the IMPACT models are TorchScript and jit.load is required + warnings.simplefilter("ignore", FutureWarning) + fn = IMPACTReg(IMPACT_SCORE_MODEL, shape=[0, 0, 0], in_channels=1, + weights=[1.0, 1.0], model_cache=d("data/impact_models")) + _IMPACT_LOSS_STATE = (fn, torch.device("cuda" if ngpu else "cpu")) + except Exception: # noqa + _IMPACT_LOSS_STATE = False + return _IMPACT_LOSS_STATE or None + + +def _impact_loss_value(warped_ct, mr): + """IMPACT semantic loss (lower = more in-sync) between the co-registered + pair: the CT inverse-transformed into MR coordinates vs the original MR, + over the full volumes (the IMPACT models standardize intensity inside + the metric, so both go in raw).""" + st = _impact_loss_scorer() + if st is None: + return None + fn, dev = st + import torch + a = torch.as_tensor(np.ascontiguousarray(warped_ct, dtype="float32"), + device=dev).unsqueeze(0).unsqueeze(0) + b = torch.as_tensor(np.ascontiguousarray(mr, dtype="float32"), + device=dev).unsqueeze(0).unsqueeze(0) + with torch.no_grad(): + return float(fn(a, b).item()) + + def _clip_image(img, lo, hi): """Clip image intensities to [lo, hi] (geometry preserved).""" a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi) @@ -153,10 +203,12 @@ def evaluate(ct, t1c, tv, T): warped_ct = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear)) warped_tv = sitk.GetArrayFromImage(sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0.0)) gap = m6._skull_ring_gap(warped_ct, ha) - # achieved-metric alignment score: Mattes MI between the CT inverse-transformed - # into MR coordinates (warped_ct) and the original MR (ha), over the MR head. + # achieved-metric alignment scores: Mattes MI between the CT + # inverse-transformed into MR coordinates (warped_ct) and the original MR + # (ha), over the MR head; IMPACT semantic loss over the full volumes. ha_f = ha.astype("float32", copy=False) mi = _mi_masked(warped_ct, ha_f, head_mask_from_image(ha_f)) + il = _impact_loss_value(warped_ct, ha_f) posw = warped_ct[warped_ct > 0] tissue_w = (warped_ct > max(2.0, float(np.percentile(posw, 2)))) if posw.size \ else np.zeros_like(warped_ct, bool) @@ -171,6 +223,7 @@ def evaluate(ct, t1c, tv, T): return { "skull_ring_gap_mm": gap, "mi_ct_mr": round(mi, 4) if mi is not None else None, + "impact_loss": round(il, 4) if il is not None else None, "tv_nvox": nvox, "tv_in_ct_frac": round(tv_in_ct_frac, 3), "alignment": "concentric" if ok_gap else "not_concentric", @@ -248,8 +301,9 @@ def _call(pipeline, ct, t1c, tv, args): gpu=args.synthmorph_gpu, threads=args.threads, extent=args.synthmorph_extent) if pipeline == "impact": - return RP.register_impact(ct, t1c, preset=args.impact_preset, gpu=args.impact_gpu, - iters=args.impact_iters, samples=args.impact_samples) + return RP.register_impact(ct, t1c, gpu=args.impact_gpu, + iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS, + samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES) if pipeline == "niftyreg": return RP.register_niftyreg(ct, t1c, threads=args.threads, backend=args.niftyreg_backend, @@ -268,6 +322,7 @@ def _run_case(pipeline, case, args): out = {"key": key, "subject": case["subject"], "date": case["date"]} if not ct_complete(ct): out.update({"skull_ring_gap_mm": None, "tv_nvox": 0, "tv_in_ct_frac": None, + "impact_loss": None, "alignment": "ct_incomplete", "accepted": False, "note": "CT head incomplete (data issue, not a pipeline failure)"}) return out @@ -313,7 +368,7 @@ def run_pipeline(pipeline, cases, args): r = results[-1] reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT" print(f" {c['key']}: mi={r.get('mi_ct_mr')} " - f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True) + f"impact_loss={r.get('impact_loss')} {reg}", flush=True) except Exception as e: # noqa errors.append({"key": c["key"], "error": repr(e)}) print(f" {c['key']}: ERROR {e!r}", flush=True) @@ -325,7 +380,7 @@ def run_pipeline(pipeline, cases, args): gaps_mid = [r["skull_ring_gap_mm"]["mid_mm"] for r in registered if r.get("skull_ring_gap_mm")] gaps_std = [r["skull_ring_gap_mm"]["std_mm"] for r in registered if r.get("skull_ring_gap_mm")] regt = [r["reg_wall_s"] for r in registered] - tvfr = [r["tv_in_ct_frac"] for r in registered if r.get("tv_in_ct_frac") is not None] + ilv = [r["impact_loss"] for r in registered if r.get("impact_loss") is not None] mifr = [r["mi_ct_mr"] for r in registered if r.get("mi_ct_mr") is not None] first = (registered or results)[0] if (registered or results) else {} row = { @@ -341,7 +396,8 @@ def run_pipeline(pipeline, cases, args): "median": round(float(np.median(gaps_mid)), 1)} if gaps_mid else None, "skull_ring_gap_std_mm": {"mean": round(float(np.mean(gaps_std)), 1), "median": round(float(np.median(gaps_std)), 1)} if gaps_std else None, - "tv_in_ct_frac": {"mean": round(float(np.mean(tvfr)), 3)} if tvfr else None, + "impact_loss": {"mean": round(float(np.mean(ilv)), 4), + "median": round(float(np.median(ilv)), 4)} if ilv else None, "mi_ct_mr": {"mean": round(float(np.mean(mifr)), 4), "median": round(float(np.median(mifr)), 4)} if mifr else None, "cases": results, "errors": errors, @@ -368,24 +424,24 @@ def _fmt(v, spec="{}", default="-"): def print_table(rows): order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)] - print("\n" + "=" * 92) + print("\n" + "=" * 96) hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} " - f"{'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}") + f"{'s/reg':>7s} {'ImpactLoss':>10s} {'MI':>7s}") print(hdr) - print("-" * 92) + print("-" * 96) for p in order: r = next(x for x in rows if x["pipeline"] == p) if not r.get("available"): print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} " - f"{'-':>7s} {'-':>7s} {'-':>7s} [{r['reason']}]") + f"{'-':>7s} {'-':>10s} {'-':>7s} [{r['reason']}]") continue rt = (r.get("reg_wall_s") or {}).get("median") - tv = (r.get("tv_in_ct_frac") or {}).get("mean") + il = (r.get("impact_loss") or {}).get("median") mi = (r.get("mi_ct_mr") or {}).get("median") print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} " f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} " - f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}") - print("=" * 92) + f"{_fmt(il, '{:.3f}'):>10s} {_fmt(mi, '{:.3f}'):>7s}") + print("=" * 96) def write_plot(rows, out_png): @@ -402,13 +458,15 @@ def write_plot(rows, out_png): return labels = list(ok_rows) regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") or 0 for p in labels] + il = [(ok_rows[p].get("impact_loss") 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"] + fig, axes = plt.subplots(1, 3, figsize=(14, 4)) + colors = ["#6ACC65", "#F5A623", "#4878CF"] for ax, (title, vals, fmt, ylab) in zip(axes, [ ("Registration time (s, median)", regt, "{:.0f}", "s"), + ("ImpactLoss: warped-CT ~ MR (semantic, median; lower = better)", il, "{:.3f}", "loss"), ("Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats")]): - bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)]) + bars = ax.bar(range(len(vals)), vals, color=colors[0]) for b, v in zip(bars, vals): ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()), ha="center", va="bottom", fontsize=8) @@ -438,7 +496,9 @@ def merge(args): summary = { "timestamp": datetime.now().isoformat(timespec="seconds"), "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI " - "(nats, warped-CT ~ MR over the MR head, higher = better); reg time " + "(nats, warped-CT ~ MR over the MR head, higher = better) and the " + "IMPACT semantic loss (vboussot/ImpactLoss, TS/M730 features, " + "warped-CT ~ MR over the full volumes, lower = better); reg time " "excludes I/O+scoring; per-case JSON retains the production skull-ring " "gap / acceptance QA; best-of-T vs T^-1 orientation gate for the " "elastix-style toolkits"), @@ -452,7 +512,7 @@ def merge(args): print("summary:", sf) write_plot(rows, os.path.join(out_dir, "bench_registration.png")) # markdown table for the report - md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | tv-in-CT | MI (nats, median) |", + md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |", "|---|---|---|---|---|---|---|---|---|"] for p in [x["pipeline"] for x in rows]: r = next(x for x in rows if x["pipeline"] == p) @@ -463,7 +523,7 @@ def merge(args): f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | " f"{r['n_error']} | " f"{(r.get('reg_wall_s') or {}).get('median')} | " - f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | " + f"{(r.get('impact_loss') or {}).get('median')} | " f"{(r.get('mi_ct_mr') or {}).get('median')} |") mf = os.path.join(out_dir, "summary.md") with open(mf, "w") as f: @@ -496,24 +556,16 @@ def parse_args(): ap.add_argument("--synthmorph-extent", type=int, default=256, help="SynthMorph isotropic registration-space extent in mm " "(192 or 256; default 256)") - ap.add_argument("--impact-preset", default=RP.IMPACT_DEFAULT_PRESET, - help="ImpactReg preset from the published database " - f"(default: {RP.IMPACT_DEFAULT_PRESET}; e.g. " - "Elastix_IMPACT_Jacobian, Generic_Rigid, " - "Generic_Rigid_BSpline, FireANTs_SyN, FireANTs_IMPACT, " - "ConvexAdam_Coarse/Composite; list live with: " - "impact-reg-konfai register --help)") ap.add_argument("--impact-gpu", type=int, default=None, help="GPU device index for IMPACT (default: auto-spread over " "the host GPUs by pid; -1 forces CPU)") ap.add_argument("--impact-iters", type=int, default=None, - help="preset max_iterations override, one value for all " - "resolutions (only effective for presets without a " - "per-resolution model matrix; default: the preset's " - "own budget, 400/300/200/200 for Elastix_IMPACT_Static)") + help="IMPACT MaximumNumberOfIterations, one value for all " + "levels (default: the reference map's per-level budget " + f"{RP.IMPACT_DEFAULT_ITERS})") ap.add_argument("--impact-samples", type=int, default=None, - help="preset spatial_samples override (NumberOfSpatialSamples; " - "default: the preset's own value, 2000)") + help="IMPACT NumberOfSpatialSamples per iteration (default: " + f"{RP.IMPACT_DEFAULT_SAMPLES})") ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"], default="auto", help="niftyreg reg_aladin build: docker = official v2.1.1 in the " diff --git a/scripts/setup_impact_konfai.sh b/scripts/setup_impact_konfai.sh deleted file mode 100644 index 31e2352..0000000 --- a/scripts/setup_impact_konfai.sh +++ /dev/null @@ -1,110 +0,0 @@ -#!/usr/bin/env bash -# One-time host setup for the IMPACT registration pipeline (impact-reg-konfai). -# -# The stock elastix-impact 1.0.0 binary (vboussot/ImpactElastix) needs -# glibc >= 2.35 and links LibTorch 2.8.0. This host (Ubuntu 20.04, glibc 2.31, -# env torch 2.14) satisfies neither, so this script builds a self-contained -# shim at ~/.cache/konfai/elastix-impact-shim: -# * glibc 2.35 loader + libs, extracted from the niftyreg:2.1.1-ubuntu image -# (ubuntu 22.04 base) -- docker is needed ONLY for this one-time copy; -# * LibTorch 2.8.0 cu128, self-contained (its own CUDA runtime), from -# download.pytorch.org; -# * a wrapper bin/elastix that runs the stock binary through the glibc 2.35 -# loader with the LibTorch 2.8.0 library path (ignoring the inherited -# LD_LIBRARY_PATH, so the env's torch 2.14 cannot leak in). -# The pipeline selects the shim via KONFAI_ELASTIX_DIR -# (src/registration_pipelines.py::_impact_elastix_root). -# -# Run from the repo root, inside the conda env: -# source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 -# pip install impact-reg-konfai # the CLI itself -# bash scripts/setup_impact_konfai.sh # this script -# -# First `register` afterwards auto-downloads from HF/GitHub: the preset app -# (VBoussot/ImpactReg), the feature models (VBoussot/impact-torchscript-models) -# and, on hosts that CAN run the stock binary directly, the binary itself. -set -euo pipefail - -KONFAI_CACHE="${KONFAI_CACHE:-$HOME/.cache/konfai}" -SHIM="$KONFAI_CACHE/elastix-impact-shim" -STOCK="$KONFAI_CACHE/elastix-impact" -GLIBC_SRC_DIR="${GLIBC_SRC_DIR:-glibc-2.35}" -GLIBC_IMG="${GLIBC_IMG:-niftyreg:2.1.1-ubuntu}" -LT_VER="2.8.0" -LT_DIR="libtorch-${LT_VER}-cu128" -LT_URL="https://download.pytorch.org/libtorch/cu128/libtorch-shared-with-deps-${LT_VER}%2Bcu128.zip" - -echo "[1/5] stock elastix-impact binary (package auto-install)" -if [ ! -x "$STOCK/bin/elastix" ]; then - python - <<'PY' -from pathlib import Path -from impact_reg_konfai.models.elastix_install import ( - ELASTIX_CACHE, get_elastix_bin, install_elastix_impact) -if not get_elastix_bin(ELASTIX_CACHE).exists(): - # flavor auto-detected (cu128 on this host); downloads + extracts the zip - install_elastix_impact(ELASTIX_CACHE, force_cuda=False, force_cpu=False) -PY -fi -[ -x "$STOCK/bin/elastix" ] || { echo "stock binary missing: $STOCK/bin/elastix"; exit 1; } - -echo "[2/5] glibc 2.35 from $GLIBC_IMG (docker, one-time)" -if [ ! -x "$KONFAI_CACHE/$GLIBC_SRC_DIR/ld-linux-x86-64.so.2" ]; then - docker run --rm --entrypoint sh -v "$KONFAI_CACHE":/dest "$GLIBC_IMG" -c \ - "rm -rf /dest/$GLIBC_SRC_DIR && mkdir -p /dest/$GLIBC_SRC_DIR && \ - cp -a /lib/x86_64-linux-gnu/. /dest/$GLIBC_SRC_DIR/ && \ - cp -a /usr/lib/x86_64-linux-gnu/libstdc++.so.6* /dest/$GLIBC_SRC_DIR/ 2>/dev/null || true" -fi - -echo "[3/5] LibTorch ${LT_VER} cu128 (self-contained)" -if [ ! -f "$SHIM/lib/$LT_DIR/libtorch/lib/libc10.so" ]; then - mkdir -p "$SHIM/lib" - curl -L --fail -o "$KONFAI_CACHE/$LT_DIR.zip" "$LT_URL" - unzip -q "$KONFAI_CACHE/$LT_DIR.zip" -d "$SHIM/lib/$LT_DIR" - rm "$KONFAI_CACHE/$LT_DIR.zip" -fi -# the with-deps zip carries the CUDA libs under hashed names; expose the public -# soname elastix links directly -LT_LIB="$SHIM/lib/$LT_DIR/libtorch/lib" -CUDART="$(ls "$LT_LIB"/libcudart-*.so.12 | head -1)" -ln -sf "$(basename "$CUDART")" "$LT_LIB/libcudart.so.12" - -echo "[4/5] shim layout + wrapper" -mkdir -p "$SHIM/bin" -cp "$STOCK/bin/elastix" "$SHIM/bin/elastix-real" -[ -d "$STOCK/lib" ] && cp -a "$STOCK/lib/." "$SHIM/lib/" -[ -d "$KONFAI_CACHE/$GLIBC_SRC_DIR" ] && cp -a "$KONFAI_CACHE/$GLIBC_SRC_DIR" "$SHIM/lib/$GLIBC_SRC_DIR" -cat > "$SHIM/bin/elastix" <<'WRAP' -#!/bin/sh -# elastix-impact 1.0.0 shim: the host glibc (2.31) is too old for the binary -# (built against glibc 2.35) and the binary expects LibTorch 2.8.0, which the -# conda env's torch (2.14) does not provide. Run it through the bundled -# glibc 2.35 loader with the bundled self-contained LibTorch 2.8.0 (cu128). -ROOT="$(cd "$(dirname "$0")/.." && pwd)" -LT="$ROOT/lib/libtorch-2.8.0-cu128/libtorch/lib" -GLIBC="$ROOT/lib/glibc-2.35" -exec env -u LD_LIBRARY_PATH -u LD_PRELOAD \ - "$GLIBC/ld-linux-x86-64.so.2" \ - --library-path "$LT:$ROOT/lib:$GLIBC" \ - "$ROOT/bin/elastix-real" "$@" -WRAP -chmod +x "$SHIM/bin/elastix" - -echo "[5/5] verify" -KONFAI_ELASTIX_DIR="$SHIM" python - <<'PY' -import os -from pathlib import Path -from impact_reg_konfai.models.elastix_install import get_elastix_bin, loader_env, try_elastix -import subprocess -root = Path(os.environ["KONFAI_ELASTIX_DIR"]) -try_elastix(root) -r = subprocess.run([str(get_elastix_bin(root).resolve()), "--extended-version"], - capture_output=True, text=True, env=loader_env(root)) -assert r.returncode == 0, r.stderr -print(r.stdout.splitlines()[0] if r.stdout else "elastix OK") -print("elastix-impact shim OK:", root) -PY - -echo -echo "Done. The impact pipeline picks the shim up automatically" -echo "(src/registration_pipelines.py::_impact_elastix_root), e.g.:" -echo " python scripts/bench_registration.py --pipelines impact --n-cases 2" \ No newline at end of file diff --git a/src/impact_loss.py b/src/impact_loss.py new file mode 100644 index 0000000..ecb2fd8 --- /dev/null +++ b/src/impact_loss.py @@ -0,0 +1,152 @@ +# Vendored from https://github.com/vboussot/ImpactLoss (IMPACT.py) +# Copyright (c) 2025 Valentin Boussot. Apache-2.0. +# +# IMPACT Reg: A Generic Semantic Loss for Multimodal Image Registration +# (Boussot et al., arXiv:2503.24121). Modification relative to upstream: +# model files are resolved from a local cache dir first (``model_cache`` / +# the ``IMPACT_MODELS_DIR`` env var) before falling back to the Hugging +# Face download, so offline hosts can use the already-cached +# ``data/impact_models`` layout. + +import os + +import torch +from huggingface_hub import hf_hub_download + +IMPACT_HF_REPO = "VBoussot/impact-torchscript-models" + + +def _resolve_model(model_name, model_cache): + """Local cache path if present, else the Hugging Face download.""" + if model_cache: + p = os.path.join(model_cache, model_name) + if os.path.exists(p): + return p + return hf_hub_download(repo_id=IMPACT_HF_REPO, filename=model_name, + repo_type="model", revision=None) # nosec B615 + + +class IMPACTReg(torch.nn.Module): + """ + IMPACT-Reg Loss Module + ------------------ + A similarity metric based on pretrained TorchScript feature extractors. + Computes a weighted L1 distance between feature maps of two images. + + Parameters + ---------- + model_name : str + Name of the TorchScript model (e.g. "TS/M730"). + shape : list[int] + Expected spatial size. Use [0, 0, 0] for no resample. + in_channels : int + Number of input channels (before adapting to the model). + weights : list[float] + Per-layer feature weights. + model_cache : str, optional + Directory the model files are looked up in before downloading. + """ + def __init__( + self, model_name: str, shape: list[int] = [0, 0], in_channels: int = 1, + weights: list[float] = [1], model_cache: str = None, + ) -> None: + super().__init__() + if model_name is None: + return + self.in_channels = in_channels + self.loss = torch.nn.L1Loss() + self.weights = weights + self.nb_layer = len(weights) + + # Load the TorchScript model (local cache first, HF download after) + model_cache = model_cache or os.environ.get("IMPACT_MODELS_DIR", "") + self.model_path = _resolve_model(model_name, model_cache) + + self.model: torch.nn.Module = torch.jit.load(self.model_path, map_location=torch.device("cpu")) # nosec B614 + self.dim = len(shape) + self.shape = shape if all(s > 0 for s in shape) else None + + dummy_input = torch.zeros((1, self.in_channels, *(self.shape if self.shape else [128] * self.dim))).to(0) + try: + out = self.model.to(0)(dummy_input, torch.tensor([self.nb_layer])) + if not isinstance(out, (list, tuple)): + raise TypeError(f"Expected model output to be a list or tuple, but got {type(out)}.") + if self.nb_layer != len(out): + raise ValueError( + f"Loss 'L1': mismatch between the number of weights " + f"({self.nb_layer}) and the number of model outputs " + f"({len(out)}). Each output must have a corresponding weight." + ) + except Exception as e: + msg = ( + f"[Model Sanity Check Failed]\n" + f"Input shape attempted: {dummy_input.shape}\n" + f"Error: {type(e).__name__}: {e}" + ) + raise RuntimeError(msg) from e + self.model = None + + def preprocessing(self, tensor: torch.Tensor) -> torch.Tensor: + """ + Resize and adjust the tensor channels to match model expectations. + + Note: + This method does not perform any intensity normalization. + Input images must be provided in their raw (non-normalized) form, + i.e., in their original intensity space (e.g., HU for CT, native + values for MRI). Intensity normalization / standardization is the + responsibility of the IMPACT models. + """ + if self.shape is not None and not all( + tensor.shape[-i - 1] == size for i, size in enumerate(reversed(self.shape[2:])) + ): + if tensor.dtype == torch.uint8: + mode = "nearest" + elif len(tensor.shape) < 4: + mode = "bilinear" + else: + mode = "trilinear" + tensor = torch.nn.functional.interpolate( + tensor, mode=mode, size=tuple(self.shape), align_corners=False + ).type(torch.float32) + + # Replicate channels if needed (e.g., 1 -> 3) + if tensor.shape[1] != self.in_channels: + tensor = tensor.repeat(tuple([1, self.in_channels] + [1 for _ in range(self.dim)])) + return [tensor, torch.tensor([self.nb_layer]), + torch.tensor([tensor.min(), tensor.max(), tensor.mean(), tensor.std()])] + + def _compute(self, output: torch.Tensor, target: torch.Tensor) -> torch.Tensor: + """Compute the feature-based L1 loss between two tensors.""" + loss = torch.zeros(1, device=output.device, dtype=torch.float32) + + output = self.preprocessing(output) + target = self.preprocessing(target) + + for weight, output_features, target_features in zip( + self.weights, self.model(*output), self.model(*target) + ): + if weight == 0: + continue + loss = loss + weight * self.loss(output_features, target_features) + return loss + + def forward(self, output: torch.Tensor, target: torch.Tensor) -> torch.Tensor: + """ + Main forward method. + If the model is 2D and inputs are 3D, compute slice-wise average loss. + """ + if self.model is None: + self.model = torch.jit.load(self.model_path) # nosec B614 + self.model.to(output.device) + self.model.eval() + + loss = torch.zeros(1, device=output.device, dtype=torch.float32) + # 3D -> 2D model case: iterate over slices + if len(output.shape) == 5 and self.dim == 2: + for i in range(output.shape[2]): + loss = loss + self._compute(output[:, :, i, ...], target[:, :, i, ...]) + loss = loss / output.shape[2] + else: + loss = self._compute(output, target) + return loss \ No newline at end of file diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index edb38ef..438159f 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -31,8 +31,6 @@ Shared rules for every pipeline: - **The T1c (MR) is always the fixed image** (moving = CT). When both inputs are the same class (both CT, or a clipped CT that no longer reads as CT) the T1c-slot (second) argument is kept as fixed; see :func:`fixed_moving`. - Every pipeline -- including ``impact``, whose ImpactReg MR_CT convention - coincides with this rule -- registers fixed=T1c / moving=CT. - **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 (the moving image) to [0, 80] HU before @@ -53,18 +51,13 @@ default center-of-mass alignment, simpleitk starts from identity, niftyreg - ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` + ``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD). ``AutomaticTransformInitialization`` (CenterOfGravity) start. -- ``register_impact`` : IMPACT (semantic-similarity metric) via the - ``impact-reg-konfai`` package (pip; the KonfAI runtime drives the - elastix-impact 5.2 binary -- no Docker). Default preset - ``Elastix_IMPACT_Static``: the published successor of the ImpactReg - ``MR_CT_TS`` reference map (elastix B-spline + per-level MIND + - TotalSegmentator M730, Static mode). Registered fixed=T1c/moving=CT (the - ImpactReg MR_CT convention) so the preset's output field (on the fixed - grid) is directly the T1c->CT contract transform; CenterOfGravity auto - init + the 0.1 sample-validity gate as raw-elastix overrides. GPU. - Full volumes, no mask (shared rules above). First use auto-downloads the - preset app, models and binary into ``~/.cache``; on hosts too old for the - stock binary, ``scripts/setup_impact_konfai.sh`` builds the shim. +- ``register_impact`` : IMPACT (semantic-similarity metric) inside the + ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + TorchScript feature + models as a native Elastix metric). Rigid 6-DOF ``EulerTransform`` driven + by the MIND + TotalSegmentator static feature combination of the ImpactReg + ``MR_CT_TS`` reference map (per-level models / scales / layer weights) + a + Mattes-MI secondary metric (w 0.4). GPU. Full volumes, no mask (shared + rules above). Needs Docker + the image + the model assets under the cache. - ``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 @@ -93,10 +86,8 @@ per toolkit. :func:`fix_orientation` resolves this by warping the CT with both ``T`` and ``T^-1`` and keeping whichever gives the smaller skull-ring concentricity gap (a pure correctness gate -- a wrongly-oriented rigid transform lands the head on the far side of the box and yields a pathological gap). It is -applied to elastix / antspyx / synthmorph, not to the trusted production -SimpleITK path. ``impact`` returns a displacement field (which has no -``GetInverse``), so it uses :func:`_orient_impact` -- the same skull-ring -decision against a sign-flipped (first-order reverse) proxy field. +applied to elastix / antspyx / synthmorph / impact, not to the trusted +production SimpleITK path. """ import os import re @@ -670,91 +661,46 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None, # --------------------------------------------------------------------------- # -# 5. IMPACT (impact-reg-konfai -- KonfAI runtime, no Docker) +# 5. IMPACT (Docker image vboussot/elastix_impact) # --------------------------------------------------------------------------- # -IMPACT_CLI = "impact-reg-konfai" -# default preset from the published ImpactReg database (VBoussot/ImpactReg, -# resolved by the CLI): elastix B-spline driven by the IMPACT semantic metric -# -- per-level MIND (L1) + TotalSegmentator M730 (deepest layer, Dice), Static -# mode (whole-image features, no backprop), 4 scales, 400/300/200/200 -# iterations, 2000 spatial samples, + Mattes MI (w 0.4) + bending energy. -# The published successor of the old MR_CT_TS reference parameter map. -IMPACT_DEFAULT_PRESET = "Elastix_IMPACT_Static" -# M6 robustness overrides -- raw elastix parameter lines the package applies -# on top of the preset's generated map (the ``parameter_overrides`` knob): -# * CenterOfGravity auto-init (instead of the preset's GeometricalCenter): -# aligns the intensity centroids -- the T1c head-crop and the CT head+neck -# FOV have their buffer centres ~50-100 mm apart, and the GeometricalCenter -# start leaves the metric with almost nothing in-buffer; -# * RequiredRatioOfValidSamples 0.1 (build default 0.25): the FOV-mismatched -# pairs keep only a minority of the random spatial samples in-buffer at the -# init; 0.1 keeps the divergence gate while letting the runs proceed (a -# genuinely diverging transform, <10% in-buffer, still aborts). -IMPACT_PARAM_OVERRIDES = ( - "AutomaticTransformInitialization=true", - "AutomaticTransformInitializationMethod=CenterOfGravity", - "RequiredRatioOfValidSamples=0.1", +IMPACT_IMAGE = "vboussot/elastix_impact" +IMPACT_ELASTIX_BIN = "/lib/elastix-install/bin/elastix" +IMPACT_MODEL_CACHE = "data/impact_models" +IMPACT_HF_MODEL_URL = ( + "https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/{rel}" ) -# The stock elastix-impact 1.0.0 binary (vboussot/ImpactElastix) needs -# glibc >= 2.35 and links LibTorch 2.8.0; a host that cannot run it directly -# (this one: glibc 2.31 + env torch 2.14) builds a self-contained shim with -# scripts/setup_impact_konfai.sh, selected via KONFAI_ELASTIX_DIR. -IMPACT_ELASTIX_SHIM = os.path.join( - os.path.expanduser("~"), ".cache", "konfai", "elastix-impact-shim") -IMPACT_ELASTIX_STOCK = os.path.join( - os.path.expanduser("~"), ".cache", "konfai", "elastix-impact") -# orientation-gate threshold: a candidate DVF is kept as-is when its -# skull-ring gap score (std_mm + mid_mm) is below this. A misdirected field -# lands the head on the far side of the box (gap ~100-300 mm), so the margin -# is wide; see _orient_impact. -IMPACT_GAP_ACCEPT = 60.0 +# The handcrafted MIND descriptor (whole-image features, layer mask "1", L1) +# combined per level with a TotalSegmentator feature model (deepest layers, +# Dice) -- the MIND+TS static combination of the ImpactReg ``MR_CT_TS`` +# reference parameter map, coarse -> fine. Per level: (TS model, voxel-size +# scale in mm, MIND layer mask, TS layer mask, MIND weight, TS weight). +IMPACT_MIND_MODEL = "MIND/R1D2_3D.pt" +IMPACT_LEVELS = ( + ("TS/M852.pt", 6.0, "1", "0000001", 0.2, 0.8), + ("TS/M852.pt", 3.0, "1", "0000001", 0.3, 0.7), + ("TS/M850.pt", 2.0, "1", "00000001", 0.6, 0.4), + ("TS/M850.pt", 2.0, "1", "00000001", 0.7, 0.3), +) +IMPACT_DEFAULT_ITERS = (400, 300, 200, 200) +IMPACT_DEFAULT_SAMPLES = 2000 +# secondary Mattes-MI metric weight (the semantic metric gets weight 1.0) +IMPACT_MI_WEIGHT = 0.4 -def _impact_cli_ok(): - """Return (ok, note): is the impact-reg-konfai CLI installed?""" - if shutil.which(IMPACT_CLI) is None: - return False, f"{IMPACT_CLI} not on PATH (pip install impact-reg-konfai)" +def _impact_docker_ok(): + """Return (ok, note): is the Docker CLI + the IMPACT image available?""" + if shutil.which("docker") is None: + return False, "docker CLI not found on PATH" try: - import importlib.metadata - ver = importlib.metadata.version("impact-reg-konfai") - except Exception: # noqa - ver = "?" - return True, f"{IMPACT_CLI} {ver}" - - -_ELASTIX_ROOT_CACHE = None - - -def _impact_elastix_root(): - """The elastix-impact install root that runs on this host, or None. - - Tried in order: a caller-set ``KONFAI_ELASTIX_DIR``, the setup-script shim - (glibc 2.35 loader + LibTorch 2.8.0, for hosts too old for the stock - binary), and the package's own auto-install cache (populated by the first - run when the host can run the stock binary directly). Each candidate is - actually executed (``elastix -h``) so a stale/broken install is skipped.""" - global _ELASTIX_ROOT_CACHE - if _ELASTIX_ROOT_CACHE is not None: - return _ELASTIX_ROOT_CACHE - from pathlib import Path - from impact_reg_konfai.models import elastix_install as _ei - - root = os.environ.get("KONFAI_ELASTIX_DIR", "") - candidates = [c for c in (root, IMPACT_ELASTIX_SHIM, IMPACT_ELASTIX_STOCK) - if c] - found = None - for c in candidates: - if not os.path.exists(os.path.join(c, "bin", "elastix")): - continue - try: - _ei.try_elastix(Path(c)) - found = c - break - except Exception: # noqa - continue - _ELASTIX_ROOT_CACHE = found - return found + r = subprocess.run(["docker", "image", "inspect", IMPACT_IMAGE], + capture_output=True, text=True, timeout=30) + if r.returncode != 0: + return False, (f"docker image {IMPACT_IMAGE} not pulled " + f"(run: docker pull {IMPACT_IMAGE})") + except Exception as e: # noqa + return False, f"docker image check failed ({type(e).__name__}: {e})" + return True, IMPACT_IMAGE _NGPUS_CACHE = None @@ -773,156 +719,196 @@ def _n_gpus(): return _NGPUS_CACHE -def _impact_dvf_from_h5(path): - """An ITK-form HDF5 transform (the ``:itktransform`` form the KonfAI - presets write) read with h5py -- the SimpleITK build in this env has no - ITK HDF5 transform reader. Returns ``(DisplacementFieldTransform, - dvf_image)``; the field is in mm, channel-last ``(z, y, x, 3)``.""" - import h5py - with h5py.File(str(path), "r") as f: - grp = f["TransformGroup/0"] - ttype = grp["TransformType"][()][0] - if isinstance(ttype, bytes): - ttype = ttype.decode() - fixed = grp["TransformFixedParameters"][()] - params = grp["TransformParameters"][()] - if not str(ttype).startswith("DisplacementFieldTransform"): - raise RuntimeError(f"impact: unsupported h5 transform type {ttype!r}") - size = fixed[:3][::-1].astype(int) # stored (x, y, z) -> array (z, y, x) - arr = params.reshape(int(size[0]), int(size[1]), int(size[2]), 3) - img = sitk.GetImageFromArray(arr) - img.SetOrigin(tuple(float(v) for v in fixed[3:6])) - img.SetSpacing(tuple(float(v) for v in fixed[6:9])) - img.SetDirection(tuple(float(v) for v in fixed[9:18])) - # this SimpleITK build invalidates the image wrapper handed to the DVF - # constructor (it comes back as an empty uint8 image), so build the - # transform from a copy and keep ``img`` for the caller - return sitk.DisplacementFieldTransform(sitk.Image(img)), img +def _impact_models_host(): + """Resolve the IMPACT feature models to files under the model cache, or + raise :class:`PipelineUnavailable` with the exact download commands for + the missing ones. Returns the de-duplicated paths relative to the cache + (which the pipeline bind-mounts at ``/Data/Models`` inside the + container).""" + cache = d(IMPACT_MODEL_CACHE) + rels = [IMPACT_MIND_MODEL] + [rel for rel, *_ in IMPACT_LEVELS] + missing = [r for r in dict.fromkeys(rels) + if not os.path.exists(os.path.join(cache, r))] + if missing: + cmds = "\n".join( + f" mkdir -p {os.path.join(cache, os.path.dirname(r))} && " + f"curl -L -o {os.path.join(cache, r)} {IMPACT_HF_MODEL_URL.format(rel=r)}" + for r in missing) + raise PipelineUnavailable( + f"impact: model(s) not found under {cache}.\nDownload once with:\n" + f"{cmds}\n(or `sh Data/Models/download_models.sh` from the " + f"ImpactLoss repo).") + return list(dict.fromkeys(rels)) -def _impact_read_transform(path): - """The preset's single-case transform output as - ``(sitk.Transform, dvf_image_or_None)``. - - KonfAI presets write the registration in whatever form the preset - declares: an ITK-form HDF5 DisplacementFieldTransform (``Transform.h5``) - by default, or an ITK transform file (``.itk.txt`` / ``.tfm``) or a vector - image (``.mha``) for other presets.""" - from pathlib import Path - p = Path(path) - if p.name.endswith(".itk.txt") or p.suffix in (".h5", ".tfm"): - if p.suffix == ".h5": - return _impact_dvf_from_h5(p) - return sitk.ReadTransform(str(p)), None - img = sitk.ReadImage(str(p)) - return (sitk.DisplacementFieldTransform( - sitk.Cast(img, sitk.sitkVectorFloat64)), img) - - -def _orient_impact(ct, t1c, cand, dvf_img): - """``(T, inverted)`` in the contract direction, for an impact candidate. - - ``cand`` is nominally T1c->CT: the preset ran with fixed=T1c / moving=CT - and its output field is on the fixed (T1c) grid. Displacement fields - cannot be inverted with ``GetInverse`` (ITK), so the direction gate does - not compare ``cand`` with ``cand^-1`` (see :func:`fix_orientation`); it - compares the skull-ring gap of ``cand`` against a reverse proxy -- the - sign-flipped field (first-order inverse, exact enough for the decision: - a misdirected DVF lands the head on the far side of the box, gap - ~100+ mm instead of tens of mm) -- and falls back to it only when the - nominal direction is clearly wrong.""" - ha = sitk.GetArrayFromImage(t1c).astype("float32") - s0 = _gap_score(_skull_gap(_warp_ct(ct, t1c, cand), ha)) - if s0 <= IMPACT_GAP_ACCEPT: - return cand, False - if dvf_img is not None: - img = sitk.GetImageFromArray(-sitk.GetArrayViewFromImage(dvf_img)) - img.CopyInformation(dvf_img) - rev = sitk.DisplacementFieldTransform(img) +def _impact_param_map(gpu, iters, samples, voxel_sizes): + """Build the Elastix-5.2 parameter map for a RIGID (Euler) IMPACT + registration: the MIND + TotalSegmentator static feature combination of + the ImpactReg ``MR_CT_TS`` reference map (one model pair per level, + coarse -> fine, see :data:`IMPACT_LEVELS`) + a secondary Mattes-MI metric + (:data:`IMPACT_MI_WEIGHT`). The image pyramid is NOT resampled (all-1 + schedules); the scale comes from the per-level ImpactVoxelSize (as the + IMPACT docs recommend). ``iters`` is a per-level tuple or a single int. + """ + n = len(voxel_sizes) + if isinstance(iters, (int,)): + iters_l = str(iters) else: - try: - rev = cand.GetInverse() - except Exception: # noqa: parametric transform with no inverse - return cand, False - s1 = _gap_score(_skull_gap(_warp_ct(ct, t1c, rev), ha)) - if s1 < s0: - return rev, True - return cand, False + iters_l = " ".join(f"{i:g}" for i in iters) + mixed = "true" if gpu >= 0 else "false" + L = ["// IMPACT semantic-similarity metric (vboussot/ImpactLoss),"] + L.append("// MIND + TotalSegmentator static combination of the") + L.append("// VBoussot/ImpactReg MR_CT_TS reference map (EulerTransform") + L.append("// instead of the reference RecursiveBSpline -- rigid)") + L.append('(Metric "Impact" "AdvancedMattesMutualInformation")') + L.append("(Metric0Weight 1)") + L.append(f"(Metric1Weight {IMPACT_MI_WEIGHT})") + L.append('(UseRelativeWeights "false")') + for i, (ts_rel, _vs, mind_mask, ts_mask, mind_w, ts_w) in enumerate(IMPACT_LEVELS): + vs = voxel_sizes[i] + L += [ + f'(ImpactModelsPath{i} "/Data/Models/{IMPACT_MIND_MODEL}" ' + f'"/Data/Models/{ts_rel}")', + f"(ImpactDimension{i} 3 3)", + f"(ImpactNumberOfChannels{i} 1 1)", + f"(ImpactPatchSize{i} 0 0 0 0 0 0)", + f"(ImpactVoxelSize{i} {vs:g} {vs:g} {vs:g} {vs:g} {vs:g} {vs:g})", + f'(ImpactLayersMask{i} "{mind_mask}" "{ts_mask}")', + f"(ImpactSubsetFeatures{i} 32 64)", + f"(ImpactPCA{i} 0 0)", + f'(ImpactDistance{i} "L1" "Dice")', + f"(ImpactLayersWeight{i} {mind_w} {ts_w})", + ] + L += [ + '(ImpactMode "Static")', + f"(ImpactGPU {gpu})", + f'(ImpactUseMixedPrecision "{mixed}")', + "(ImpactFeaturesMapUpdateInterval -1)", + '(ImpactWriteFeatureMaps "false")', + # safety gate: abort if < ratio of the spatial samples map inside the + # moving buffer. With fixed=T1c / moving=CT the samples (drawn from + # the tight T1c head crop) map mostly inside the large CT head-neck + # FOV, so the gate rarely fires; 0.1 (vs the build default 0.25) + # keeps a divergence gate for the FOV mismatch. + '(CheckNumberOfSamples "true")', + "(RequiredRatioOfValidSamples 0.1)", + "", + "// rigid 6-DOF to be optimized (auto center init, no -t0)", + '(Transform "EulerTransform")', + '(TransformBase "ElastixTransformBase")', + '(TransformParameterReadoutFormat "AllParametersOneLine")', + '(FixedCenteredTransform "true")', + '(HowToCombineTransforms "Compose")', + '(AutomaticTransformInitialization "true")', + # align intensity centroids (heads), not buffer corners -- the T1c + # head crop and the CT head+neck FOV have their buffer centres far + # apart, so GeometricalCenter leaves the heads ~100 mm apart and the + # metric aborts ("too many samples outside buffer") + '(AutomaticTransformInitializationMethod "CenterOfGravity")', + '(AutomaticScalesEstimation "true")', + "", + "(ImageDimension 3)", + '(FixedInternalImagePixelType "float")', + '(MovingInternalImagePixelType "float")', + '(UseDirectionCosines "true")', + "", + '(Registration "MultiMetricMultiResolutionRegistration")', + f"(NumberOfResolutions {n})", + f"(MaximumNumberOfIterations {iters_l})", + f"(NumberOfSpatialSamples {samples})", + "", + '(FixedImagePyramid "FixedGenericImagePyramid")', + '(MovingImagePyramid "MovingGenericImagePyramid")', + f"(FixedImagePyramidRescaleSchedule {' '.join('1' for _ in range(3 * n))})", + f"(MovingImagePyramidRescaleSchedule {' '.join('1' for _ in range(3 * n))})", + "", + '(Optimizer "AdaptiveStochasticGradientDescent")', + "(MaximumNumberOfSamplingAttempts 8)", + '(UseAdaptiveStepSizes "true")', + '(UseMultiThreadingForMetrics "true")', + '(ASGDParameterEstimationMethod "DisplacementDistribution")', + "(SigmoidInitialTime 0.0)", + '(NoiseCompensation "true")', + "(NumberOfSamplesForExactGradient 4096)", + '(AutomaticParameterEstimation "true")', + "", + '(ImageSampler "RandomCoordinate")', + '(NewSamplesEveryIteration "true")', + '(Interpolator "BSplineInterpolator")', + "(BSplineInterpolationOrder 3)", + '(UseRandomSampleRegion "false")', + "", + '(WriteResultImage "false")', + "(DefaultPixelValue -1024)", + '(ResultImageFormat "mha")', + '(ITKTransformOutputFileNameExtension "itk.txt")', + '(WriteITKCompositeTransform "true")', + '(UseBinaryFormatForTransformationParameters "false")', + '(WriteTransformParametersEachIteration "false")', + '(WriteTransformParametersEachResolution "false")', + ] + return "\n".join(L) + "\n" -def register_impact(ct, t1c, preset=None, gpu=None, iters=None, samples=None, +def register_impact(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, + samples=IMPACT_DEFAULT_SAMPLES, threads=12, keep_scratch=False): - """IMPACT (semantic-similarity metric) CT->T1c via the - ``impact-reg-konfai`` package (pip; KonfAI runtime driving the - elastix-impact 5.2 binary; no Docker). + """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) 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. The default preset ``Elastix_IMPACT_Static`` is the published - successor of the old ImpactReg ``MR_CT_TS`` reference parameter map: an - elastix **B-spline** (4 scales, FinalGridSpacing 14 mm, 400/300/200/200 - iterations, 2000 spatial samples, bending-energy regularizer) driven by - the IMPACT metric -- a per-level **MIND** (L1) + **TotalSegmentator M730** - (deepest layer, Dice) combination in **Static** mode (whole-image - features, no backprop), plus a Mattes-MI secondary metric (w 0.4). Other - presets (``Elastix_IMPACT_Jacobian``, ``Generic_Rigid``, - ``Generic_Rigid_BSpline``, ``FireANTs_*``, ``ConvexAdam_*``) are - selectable with ``preset=``; the live database is listed by - ``impact-reg-konfai register --help``. + one: pretrained TorchScript feature models are run on the images and the + registration minimizes the distance between the extracted feature maps. + The metric here is the MIND + TotalSegmentator **static** combination of + the ImpactReg ``MR_CT_TS`` reference map (``MR_CT_TS/ParameterMap_MRI_TS`` + in ``VBoussot/ImpactReg``): per level the MIND handcrafted descriptor + (``MIND/R1D2_3D.pt``, L1) is combined with a TotalSegmentator deepest-layer + feature (``TS/M852.pt`` coarse, ``TS/M850.pt`` fine, Dice), layer weights + 0.2/0.8 -> 0.7/0.3, at voxel sizes 6/3/2/2 mm (scales with the moving + image's spacing), plus a secondary Mattes-MI metric (weight 0.4). The + models standardize intensity inside the metric (zero mean / unit variance + per the model card), so no pipeline-side pre-normalization is applied -- + the volumes go in in their raw intensity space. Intensity standardization + is affine-invariant, so pre-normalizing would be an exact no-op anyway. -Fixed/moving: the pair is registered in the ImpactReg **MR_CT** - convention (fixed=T1c, moving=CT), which coincides with the shared - fixed-image rule (the T1c is the fixed image). The preset's output is a - displacement field on the **fixed** grid mapping fixed->moving, so the - field comes out directly in the contract direction T1c->CT (displacement - fields cannot be inverted with ``GetInverse``, which the fixed=CT - direction would have required). That also keeps the metric's - sample-validity gate healthy: samples drawn from the tight T1c head crop - (fixed) map mostly inside the large CT head+neck FOV (moving), unlike - fixed=CT where only ~14% of the samples are in-buffer at the auto-init. - The transform is the preset's native deformable B-spline (a DVF) -- the - contract (warp the CT concentric with the MR head) is unchanged. + The optimized transform is a rigid 6-DOF ``EulerTransform`` (the 5.x name + for the classic 3D rigid) with ASGD -- no image resampling pyramid, the + scale comes from the per-level ImpactVoxelSize. No manual centroid init: + elastix's ``AutomaticTransformInitialization`` (CenterOfGravity) aligns + the intensity centroids. - M6 robustness overrides (raw elastix lines the package applies on top of - the preset's generated map): CenterOfGravity auto init (the preset's - GeometricalCenter start leaves the intensity centroids ~50-100 mm apart - on the head+neck FOV mismatch) and ``RequiredRatioOfValidSamples 0.1`` - (the build default 0.25 aborts the FOV-mismatched pairs at the init; 0.1 - keeps the divergence gate). + Like every other pipeline this runs on the **full volumes with + fixed=T1c / moving=CT** (the shared fixed-image rule; no mask, no crop / + pad). The metric draws its random spatial samples from the fixed (T1c) + buffer and evaluates them where they map inside the moving (CT) buffer; + the samples of the tight T1c head crop map mostly inside the large CT + head-neck FOV, so the 0.1 sample-validity gate rarely fires. The + returned transform satisfies the unified contract (T1c->CT); + :func:`fix_orientation` verifies the direction against the skull ring. - The package auto-downloads its one-time assets from HF/GitHub on first - use: the preset app (``VBoussot/ImpactReg``), the TorchScript feature - models (``VBoussot/impact-torchscript-models``) and the elastix-impact - binary (``vboussot/ImpactElastix``) into ``~/.cache``. The binary needs - glibc >= 2.35 and LibTorch 2.8.0; on a host that cannot run it directly - (this one: glibc 2.31 + torch 2.14) ``scripts/setup_impact_konfai.sh`` - builds a self-contained shim, selected via ``KONFAI_ELASTIX_DIR``. + Budget: the default ``iters=(400, 300, 200, 200)`` per level / + ``samples=2000`` is the reference map's budget (~75-80 s/case on the GPU). + Pass a single int to ``iters`` to use one budget for every level. A GPU + is used when one is visible (auto-spread over the host GPUs by pid); + pass ``gpu=-1`` to force CPU (mixed precision is then turned off). + The result is read back from ``TransformParameters.0.itk.txt`` (an + ``Euler3DTransform`` in the world physical frame). - Budget: ``iters`` / ``samples`` forward to the preset's - ``max_iterations`` / ``spatial_samples`` knobs (only effective for - presets without a per-resolution model matrix; the default preset's - matrix wins). A GPU is used when one is visible (auto-spread over the - host GPUs by pid); pass ``gpu=-1`` to force CPU (the engine then turns - mixed precision off). The result is read back from - ``/P000/.h5`` (an ITK-form DisplacementFieldTransform on the - T1c grid) and direction-checked by :func:`_orient_impact` (skull-ring - gap vs the reverse-proxy field). + Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI, + the ``vboussot/elastix_impact`` image, and the model assets under the + cache (``data/impact_models``, auto-created; see the exact download + commands in the error). """ import time - ok, note = _impact_cli_ok() + ok, note = _impact_docker_ok() if not ok: raise PipelineUnavailable(f"impact: {note}") - root = _impact_elastix_root() - if root is None: - raise PipelineUnavailable( - "impact: no runnable elastix-impact binary on this host (the stock " - "build needs glibc >= 2.35 + LibTorch 2.8.0). Run " - "scripts/setup_impact_konfai.sh once, or point KONFAI_ELASTIX_DIR " - "at a working install.") - preset = preset or IMPACT_DEFAULT_PRESET + model_rels = _impact_models_host() - # fixed=T1c (MR), moving=CT: see the docstring (contract direction comes - # out of the preset directly; MR_CT convention = the shared rule) + # full volumes, no mask / crop / pad; the T1c is the fixed image (shared + # rule). No manual init -- elastix's AutomaticTransformInitialization + # (CenterOfGravity) aligns the image centers. fixed, moving = fixed_moving(ct, t1c) # GPU: auto-spread across the host GPUs by pid, or forced value @@ -930,78 +916,85 @@ Fixed/moving: the pair is registered in the ImpactReg **MR_CT** ngpu = _n_gpus() gpu = (os.getpid() % ngpu) if ngpu else -1 + # per-level feature-extraction scales (6/3/2/2 mm in the moving image's + # base-spacing units) + base = float(np.median(moving.GetSpacing())) + n = len(IMPACT_LEVELS) + voxel_sizes = [base * scale for _, scale, *_ in IMPACT_LEVELS] + + # per-case scratch on a local (non-NFS) tmp fs: Data/ + Out/ tmp = tempfile.mkdtemp(prefix="impact_reg_") - tail = "" + datadir = os.path.join(tmp, "Data") + outdir = os.path.join(tmp, "Out") + os.makedirs(datadir, exist_ok=True) + os.makedirs(outdir, exist_ok=True) try: - fixed_path = os.path.join(tmp, "Fixed.mha") - moving_path = os.path.join(tmp, "Moving.mha") - outdir = os.path.join(tmp, "Out") - sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), fixed_path, True) - sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), moving_path, True) + sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), + os.path.join(datadir, "Fixed_image.mha"), True) + sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), + os.path.join(datadir, "Moving_image.mha"), True) + param_text = _impact_param_map(gpu, iters, samples, voxel_sizes) + with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: + f.write(param_text) - cmd = [IMPACT_CLI, "register", preset, - "-f", fixed_path, "-m", moving_path, - "-o", outdir, - "--fields-only", # we derive the warp ourselves - "--tmp-dir", os.path.join(tmp, "Work")] + # 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) + param_dump = os.path.join(dump_dir, f"ParameterMap_MIND+TS_{os.getpid()}.txt") + with open(param_dump, "w") as f: + f.write(param_text) + + cmd = ["docker", "run", "--rm"] if gpu >= 0: - cmd += ["--gpu", str(gpu)] - # preset tuning (YAML values, forwarded by the CLI to the preset app): - # the M6 robustness overrides, always; the budget knobs, when set - # (no-op for presets whose per-resolution matrix wins) - ov = "[" + ",".join(f'"{o}"' for o in IMPACT_PARAM_OVERRIDES) + "]" - cmd += ["--set", f"parameter_overrides={ov}"] - if iters is not None: - cmd += ["--set", f"max_iterations={int(iters)}"] - if samples is not None: - cmd += ["--set", f"spatial_samples={int(samples)}"] - - env = os.environ.copy() - env["KONFAI_ELASTIX_DIR"] = root - + cmd += ["--gpus", "all"] + cmd += [ + "-v", f"{datadir}:/Data", + "-v", f"{outdir}:/Out", + "-v", f"{d(IMPACT_MODEL_CACHE)}:/Data/Models", + IMPACT_IMAGE, IMPACT_ELASTIX_BIN, + "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha", + "-p", "/Data/ParameterMap.txt", + "-out", "/Out/", "-threads", str(threads), + ] t0 = time.time() - proc = subprocess.run(cmd, capture_output=True, text=True, - timeout=7200, env=env) + proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200) dt = time.time() - t0 - tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) if proc.returncode != 0: - raise RuntimeError(f"impact: {IMPACT_CLI} failed " - f"(rc={proc.returncode}):\n{tail}") + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"impact: docker elastix failed (rc={proc.returncode}):\n{tail}") - case_dirs = sorted(p for p in os.listdir(outdir) - if os.path.isdir(os.path.join(outdir, p))) - if not case_dirs: - raise RuntimeError(f"impact: no case output under {outdir}:\n{tail}") - case_dir = os.path.join(outdir, case_dirs[0]) - entries = [p for p in os.listdir(case_dir) - if os.path.isfile(os.path.join(case_dir, p)) - and not p.startswith("Moved")] - if len(entries) != 1: - raise RuntimeError( - f"impact: expected exactly one transform output in {case_dir}, " - f"got {entries}:\n{tail}") - T_main, dvf_img = _impact_read_transform(os.path.join(case_dir, entries[0])) - T, inverted = _orient_impact(ct, t1c, T_main, dvf_img) + itk_tf = os.path.join(outdir, "TransformParameters.0.itk.txt") + if not os.path.exists(itk_tf): + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + 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); with + # fixed=T1c / moving=CT the direction it holds is resolved empirically + # by fix_orientation 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) + models = " + ".join(model_rels) meta = { "pipeline": "impact", - "toolkit": (f"impact-reg-konfai (KonfAI; elastix-impact 5.2, " - f"preset {preset})"), - "metric": ("IMPACT semantic (MIND R1D2_3D L1 + TotalSegmentator M730 " - "deepest-layer Dice, per-level weights, Static) + " - "Mattes MI (w 0.4) + bending energy"), - "transform": (f"Deformable (preset-native B-spline DVF on the T1c grid; " - f"CenterOfGravity auto init; " + "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", + "metric": (f"IMPACT semantic ({models}, L1+Dice, Static) + " + f"Mattes MI (w {IMPACT_MI_WEIGHT})"), + "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": _modality_name(fixed), "moving": _modality_name(moving), "seconds_reg": round(dt, 1), - "preset": preset, "gpu": gpu, "iters": iters, "samples": samples, - "parameter_overrides": list(IMPACT_PARAM_OVERRIDES), - "elastix_root": os.path.basename(root), + "docker_image": IMPACT_IMAGE, "model": models, "gpu": gpu, + "iters": list(iters) if not isinstance(iters, int) else iters, + "samples": samples, "param_dump": param_dump, } if keep_scratch: meta["scratch"] = tmp @@ -1272,21 +1265,18 @@ def availability(): else: out[name] = (True, f"Docker {SYNTHMORPH_IMAGE}") elif name == "impact": - ok, note = _impact_cli_ok() + ok, note = _impact_docker_ok() if not ok: out[name] = (False, note) continue - root = _impact_elastix_root() - if root is None: - out[name] = (False, - note + "; no runnable elastix-impact binary " - "(needs glibc >= 2.35 + LibTorch 2.8.0) -- run " - "scripts/setup_impact_konfai.sh once") - else: - out[name] = (True, - note + "; elastix-impact under " - f"{os.path.basename(root)}; preset " - f"{IMPACT_DEFAULT_PRESET}") + try: + rels = _impact_models_host() + except PipelineUnavailable as e: + out[name] = (False, str(e)) + continue + out[name] = (True, f"Docker {IMPACT_IMAGE} + {len(rels)} models " + f"in {d(IMPACT_MODEL_CACHE)} (rigid MIND+TS " + f"static config)") except Exception as e: # noqa out[name] = (False, repr(e)) return out \ No newline at end of file