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

refactor(impact): replace KonfAI pip runtime with Docker elastix_impact image

Browse source 
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.
This commit is contained in:
Furen Xiao 2026-09-30 21:23:18 +08:00
parent 5e525c07c8
commit 77653e4018
6 changed files with 786 additions and 732 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
AGENTS.md
View file

@ -10,7 +10,7 @@ Activate with:
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 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 ## Project

538
REGISTRATION_BENCHMARK.md
View file

@ -1,8 +1,7 @@
# Registration Pipeline Benchmark — CT → T1c (M6) # Registration Pipeline Benchmark — CT → T1c (M6)
Benchmarks **six cross-modality registration pipelines** (five rigid; Benchmarks **six cross-modality registration pipelines — all rigid** (6-DOF)
`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
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 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 `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 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 exactly as given (no `FixedMask` / `MovingMask`, no head-box resampling, no
padding). padding).
- **The T1c (MR) is always the fixed image** (moving = CT) for every - **The T1c (MR) is always the fixed image** (moving = CT) for every
pipeline. `impact`'s ImpactReg **MR_CT** convention coincides with this pipeline. (A CT clipped to a non-negative window reads as MR, so the T1c
rule: its output displacement field (defined on the fixed grid) is the slot decides; both-CT pairs keep the given order.)
contract transform T1c→CT directly — see the IMPACT section. (A CT clipped - **Every transform is rigid (6-DOF)** — Euler / rigid-matrix parameterizations
to a non-negative window reads as MR, so the T1c slot decides; both-CT per toolkit. No affine, no deformable.
pairs keep the given order.)
All pipelines also use each toolkit's **native automatic transform All pipelines also use each toolkit's **native automatic transform
initialization** (no manual centroid offset): elastix / impact use 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 center-of-mass alignment, simpleitk starts from identity, niftyreg
(reg_aladin) initializes from the image centres, synthmorph needs none. One (reg_aladin) initializes from the image centres, synthmorph needs none. One
consequence of the no-mask rule: the IMPACT metric only evaluates samples that 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 map inside the moving buffer. With fixed=T1c the samples are drawn from the
are drawn from the tight T1c head crop (fixed, ~220×220×150 mm) and map mostly tight T1c head crop (fixed, ~220×220×150 mm) and map mostly inside the large
inside the large CT head-neck FOV (moving, ~300×300×251 mm), so the gate CT head-neck FOV (moving, ~300×300×251 mm), so the gate rarely fires; the
rarely fires; the pipeline still sets `RequiredRatioOfValidSamples 0.1` (vs pipeline still sets `RequiredRatioOfValidSamples 0.1` (vs the build default
the build default 0.25) as a safety override for the FOV mismatch — see the 0.25) as a safety override for the FOV mismatch — see the IMPACT section.
IMPACT section. The comparison is on each toolkit's native optimizer / metric The comparison is on each toolkit's native optimizer / metric / global-search
/ global-search capability. capability.
## The six pipelines ## The six pipelines
@ -58,16 +56,15 @@ initialization** (no manual centroid offset): elastix / impact use elastix's
`AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default `AutomaticTransformInitialization` (CenterOfGravity), antspyx uses its default
center-of-mass alignment, simpleitk starts from identity, niftyreg center-of-mass alignment, simpleitk starts from identity, niftyreg
initializes from the image centres, and synthmorph needs none (global deep initializes from the image centres, and synthmorph needs none (global deep
method). **All use fixed = T1c, moving = CT** (shared rule above; method). **All use fixed = T1c, moving = CT** (shared rule above). `niftyreg`
`impact`'s ImpactReg MR_CT convention coincides with it). `niftyreg` and and `synthmorph` clip the CT (the moving image) to [0, 80] HU before
`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). registration (per the CT-clipping A/B below; `clip_ct=None` disables it).
| # | pipeline | fixed / moving | init | metric | transform | status in this env | | # | 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) | | 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) | | 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`) | | 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 | | 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) | | 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 Benchmark budget (best-results settings): `elastix` and `simpleitk` run with
**MaximumNumberOfIterations = 2000** (per resolution / per pyramid level), **MaximumNumberOfIterations = 2000** (per resolution / per pyramid level),
with the sampler and metric left at their defaults (**NumberOfSpatialSamples with the sampler and metric left at their defaults (**NumberOfSpatialSamples
5000**, **NumberOfHistogramBins 32**). `impact` runs the published ImpactReg 5000**, **NumberOfHistogramBins 32**). `impact` runs the MIND + TotalSegmentator
preset **`Elastix_IMPACT_Static`** (the `impact-reg-konfai` package) — an **static** feature combination of the ImpactReg `MR_CT_TS` reference parameter
elastix B-spline (4 scales, FinalGridSpacing 14 mm, per-level map (per `VBoussot/ImpactReg`) as a **rigid** `EulerTransform`: 4 levels,
**MaximumNumberOfIterations 400/300/200/200**, **NumberOfSpatialSamples per-level voxel sizes 6/3/2/2 mm (no image-resampling pyramid), per-level
2000**) driven by the IMPACT metric: a per-level **MIND (L1) + **MaximumNumberOfIterations 400/300/200/200**, **NumberOfSpatialSamples 2000**,
TotalSegmentator M730 (deepest layer, Dice)** combination, layer weights + Mattes MI weight 0.4, CenterOfGravity auto-init,
0.2/0.8 → 0.7/0.3, Static mode, + Mattes MI weight 0.4 + bending energy — `RequiredRatioOfValidSamples 0.1`. `antspyx`, `niftyreg` and `synthmorph` keep
with CenterOfGravity auto-init and `RequiredRatioOfValidSamples 0.1` applied their native defaults.
as raw-elastix overrides for the M6 pairs. `antspyx`, `niftyreg` and
`synthmorph` keep their native defaults.
`simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid, `simpleitk` has two budgets: `--simpleitk-mode benchmark` (single pyramid,
identity start — the default, comparable to the other tools) and 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) ## Metrics (identical for every pipeline)
The headline alignment score is **`mi_ct_mr`** — the registration's *achieved* Alignment is scored two ways, both as the registration's *achieved* metric read
Mattes MI, read off as a standalone score: 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 - **`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**, 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 restricted to the MR head and each image percentile-binned independently
(SimpleITK Mattes-MI style). Higher = the warped-CT and MR intensities are (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 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 rotated/translated away from the MR.
registration time and GTV-in-CT coverage). - **`impact_loss`** — the IMPACT semantic loss from
- **`tv_in_ct_frac`** — fraction of the dilated warped GTV inside the warped-CT [vboussot/ImpactLoss](https://github.com/vboussot/ImpactLoss) (the metric
head tissue (the slab must straddle the lesion). 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). - **`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 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 (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 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 CT) at mid and low levels, plus the angular std. `mid` measures
translational concentricity; `std` measures rotational concentricity. translational concentricity; `std` measures rotational concentricity.
- **`tv_nvox`** — warped GTV volume. - **`tv_nvox`** — warped GTV volume.
- **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the - **`accepted`** — the production acceptance rule of `02_build_m6_dataset.py`
production acceptance rule; aligned ⇔ `|mid| ≤ 20 mm` and `std ≤ 12 mm` — (warped GTV non-empty and concentric with the head; aligned ⇔ `|mid| ≤ 20 mm`
a small negative `mid` is fine: on full-volume registration the CT skull and `std ≤ 12 mm` — a small negative `mid` is fine: on full-volume
ring can land slightly inside the MR silhouette, which is equally registration the CT skull ring can land slightly inside the MR silhouette,
concentric). which is equally concentric).
## Results (10 cases, no clip) ## Results (10 cases, no clip)
Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ MR over 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 the MR head; higher = better) **and** the **ImpactLoss** semantic score
comparable on the same score scale. Full volumes, no mask / crop / pad, (warped-CT ~ MR over the full volumes; lower = better), with the
fixed = T1c / moving = CT (shared rule; `impact`'s MR_CT convention **unclipped** CT, so all rows are comparable on the same score scales. Full
coincides with it). `niftyreg` and `synthmorph` clip the CT (the moving volumes, no mask / crop / pad, fixed = T1c / moving = CT (shared rule).
image) to [0, 80] HU *before registration* (per the CT-clipping A/B below; `niftyreg` and `synthmorph` clip the CT (the moving image) to [0, 80] HU
`clip_ct=None` disables it); `antspyx`, `elastix` and `simpleitk` register *before registration* (per the CT-clipping A/B below; `clip_ct=None`
the CT as-is. `niftyreg` is the v2.1.1 Docker build on the GPU (`-platf 1`); disables it); `antspyx`, `elastix`, `impact` and `simpleitk` register the CT
`impact` runs the `impact-reg-konfai` package on the GPU (see the IMPACT as-is. `niftyreg` is the v2.1.1 Docker build on the GPU (`-platf 1`);
section). `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** | | **antspyx** | ANTsPy 0.63 | 9 | 16.7 | 0.157 / 0.155 | **0.652 / 0.650** |
| **elastix** | ITK-Elastix 0.25.4 | 9 | 9.8 | 0.93 | **0.646 / 0.638** | | **elastix** | ITK-Elastix 0.25.4 | 9 | 11.7 | 0.157 / 0.156 | 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** | | **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 | 16.2 | 0.93 | **0.292 / 0.339** | | **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 | 9.6 | 0.93 | **0.647 / 0.642** | | **simpleitk** | SimpleITK 2.5.6 | 9 | 10.8 | **0.156 / 0.156** | 0.647 / 0.642 |
| **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 27.9 | 0.93 | **0.277 / 0.328** | | **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 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 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 rule); the CT-fixed-era numbers and the retired konfai-preset run are in the
is the `impact-reg-konfai` deformable `Elastix_IMPACT_Static` preset (the git history / archived JSONs. The `impact` column is the Docker rigid
retired rigid Docker config is archived in MIND+TS static config (the CT-fixed-era Docker run is archived in
`results/bench_registration/impact_docker_legacy.json`). `results/bench_registration/impact_docker_legacy.json`).
| case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph | | case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph |
|---|---|---|---|---|---|---| |---|---|---|---|---|---|---|
| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT |
| m6_22IRRKU6_20240115 | **0.648** | 0.644 | 0.542 | 0.292 | 0.647 | 0.277 | | m6_22IRRKU6_20240115 | **0.648** | 0.644 | 0.645 | 0.292 | 0.647 | 0.277 |
| m6_23XIMO4M_20230222 | **0.654** | 0.652 | 0.543 | 0.632 | **0.654** | 0.616 | | m6_23XIMO4M_20230222 | **0.655** | 0.652 | 0.654 | 0.632 | 0.654 | 0.616 |
| m6_24V62QQM_20230223 | 0.703 | 0.698 | 0.667 | 0.260 | **0.704** | 0.256 | | m6_24V62QQM_20230223 | 0.703 | 0.698 | 0.700 | 0.260 | **0.704** | 0.256 |
| m6_24V62QQM_20231017 | **0.601** | 0.590 | 0.587 | 0.493 | 0.597 | 0.484 | | m6_24V62QQM_20231017 | **0.602** | 0.590 | 0.596 | 0.493 | 0.597 | 0.484 |
| m6_24V62QQM_20240820 | **0.652** | 0.646 | 0.612 | 0.118 | 0.638 | 0.111 | | m6_24V62QQM_20240820 | **0.652** | 0.646 | 0.652 | 0.118 | 0.638 | 0.111 |
| m6_26TYQI3C_20221027 | **0.718** | 0.650 | 0.717 | 0.609 | 0.670 | 0.590 | | m6_26TYQI3C_20221027 | **0.718** | 0.650 | 0.658 | 0.609 | 0.670 | 0.590 |
| m6_2BBN5XA4_20250411 | 0.577 | 0.577 | 0.258 | 0.069 | **0.579** | 0.061 | | m6_2BBN5XA4_20250411 | 0.577 | 0.577 | 0.572 | 0.069 | **0.579** | 0.061 |
| m6_2C4TKI6H_20240506 | **0.645** | 0.638 | 0.578 | 0.233 | 0.636 | 0.231 | | m6_2C4TKI6H_20240506 | **0.645** | 0.638 | 0.639 | 0.233 | 0.636 | 0.231 |
| m6_2CDNN7BH_20240102 | **0.655** | 0.650 | 0.523 | 0.340 | 0.652 | 0.323 | | 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 ### Interpretation
- **ANTsPy leads, with SimpleITK and Elastix within 0.006 of it** (median MI - **The top four are a statistical tie** (median MI 0.646–0.652, median
0.652 / 0.647 / 0.646): on the well-positioned cases the three cluster at ImpactLoss 0.156–0.157): `antspyx`, `impact`, `simpleitk` and `elastix`
0.60–0.70 and all now survive the harder repositioned cases (0.58–0.65 on all recover the gross offset on every case, including the two ~30°
`m6_2BBN5XA4` / `m6_22IRRKU6`). repositioned ones (`m6_2BBN5XA4` 0.572–0.579, `m6_24V62QQM_20240820`
- **Switching to fixed = T1c / moving = CT lifted Elastix sharply** 0.638–0.652). The two alignment scores — an intensity histogram and a deep
(0.517 → 0.646 median, and it is now #3 rather than #4). The Mattes-MI semantic one — agree on the ranking.
pipelines sample spatial pairs from the fixed buffer and discard the ones - **IMPACT (Docker, rigid, T1c-fixed) is in the top tier** at median MI
that map outside the moving image; with the small T1c head crop as fixed and 0.652 (0.6523) — best or within 0.006 of the best on 8 of 9 cases
the large CT head+neck FOV as moving, (nearly) all fixed samples are (`m6_24V62QQM_20230223` 0.700, `m6_26TYQI3C` 0.658 where it trails
in-buffer, versus only a minority when the big CT was fixed. Elastix's ANTsPy's 0.718). This is a large step
repositioned cases jumped (`m6_2BBN5XA4` 0.287→0.577, `m6_2C4TKI6H` up from the legacy CT-fixed Docker config (median 0.467, archived in
0.401→0.638, `m6_23XIMO4M` 0.416→0.652). The other pipelines are essentially `impact_docker_legacy.json`): the T1c-fixed direction fixes the
unchanged by the swap (antspyx 0.653→0.652, simpleitk 0.651→0.647, niftyreg sample-coverage problem the same way it fixed Elastix (the metric draws
0.292→0.292, synthmorph 0.277→0.277) — they were not sample-coverage samples from the fixed T1c crop and evaluates them where they map into the
limited. 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 - **The grossly repositioned cases still separate the leaders from the
tail.** On `m6_2BBN5XA4` / `m6_24V62QQM_20240820` (~30° inter-session tilt) 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 the top four hold 0.57–0.65, while niftyreg (0.069 / 0.118) and synthmorph
synthmorph (0.06–0.11) — the single-start block-matching and pretrained-MRI (0.061 / 0.111) — the single-start block-matching and pretrained-MRI
methods — lock onto local minima. This is exactly why production methods — lock onto local minima. This is exactly why production
`02_build_m6_dataset.py` uses a **multi-start tilt grid**. `02_build_m6_dataset.py` uses a **multi-start tilt grid**.
- **IMPACT now runs the published `Elastix_IMPACT_Static` preset** - **NiftyReg (reg_aladin) is second-from-bottom**: single-start symmetric
(`impact-reg-konfai` package; the old Docker backend is retired): per-level block matching is strong on moderate cases (`m6_23XIMO4M` 0.632;
**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;
`m6_26TYQI3C` 0.609) but falls into the same tilt local minima `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 (`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 on the GPU (`-platf 1`) is ~15 s/case (vs ~4 min CPU for the conda 1.5.69
CPU for the conda 1.5.69 build). It clips the CT to [0,80] HU by default build). It clips the CT to [0,80] HU by default (per the A/B below) —
(per the A/B below) — neutral on the unclipped score scale neutral on the unclipped score scale (0.296→0.292 median).
(0.296→0.292 median).
- **SynthMorph is weakest of the six on these cross-modality CT pairs** - **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 (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 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 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. 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). It needs **no centroid init** (global deep method).
- **`tv-in-CT` is ~0.93 for every pipeline**: the warped GTV lands inside the - **Speed** (median wall, 3–4 cases in parallel, 24 threads, T1c-fixed):
CT head everywhere, so coverage is not the discriminator — MI is. SimpleITK 10.8 s, elastix 11.7 s, NiftyReg 15.4 s (v2.1.1, GPU block
- **Speed** (median wall, 4 cases in parallel, 24 threads, T1c-fixed re-run): matching), ANTs 16.7 s, SynthMorph 22.2 s (GPU forward pass), **IMPACT
SimpleITK 9.6 s, elastix 9.8 s, **ANTs 15.5 s** (down from 36.6 s — the 111.2 s** (Docker, static feature mode, no backprop, 400/300/200/200
smaller fixed image also shaves ANTs' per-level work), **NiftyReg 16.2 s** rigid ASGD iterations).
(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).
- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
per-case numbers carry run-to-run variance (this table is one draw); per-case numbers carry run-to-run variance (this table is one draw);
SimpleITK's MI (and registration) is deterministic at fixed thread count; SimpleITK's MI (and registration) is deterministic at fixed thread count;
the in-process Elastix ASGD is deterministic given its fixed seed; the in-process Elastix ASGD is deterministic given its fixed seed;
SynthMorph is deterministic; NiftyReg is deterministic on the CPU 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, (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 ~±0.01–0.05 MI on the weak cases). IMPACT in the **Static** config
(`Elastix_IMPACT_Jacobian`) is the most stochastic**: it draws random re-samples spatial coordinates every iteration, but the feature maps are
patches through the network every iteration, so the final pose changes run deterministic given the pose. The earlier Jacobian-mode Docker experiments
to run — treat such numbers as one sample from a distribution. The default (archived in `results/bench_registration/impact_*.json`) resampled patches
Static preset re-samples coordinates every iteration as well, but the through the network each iteration, so their per-case numbers are one
feature maps are deterministic given the pose. 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`) ## CT clipping A/B (`--clip-ct 0 80`)
@ -328,15 +333,13 @@ not clip.**
```bash ```bash
# main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1) # main env — SimpleITK + Elastix + ANTsPy + NiftyReg (niftyreg via Docker v2.1.1)
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 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 --n-cases 10 --workers 4 --threads 24 --simpleitk-mode benchmark
# force the conda 1.5.69 CPU build: --niftyreg-backend conda # 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 # impact (Docker vboussot/elastix_impact; GPU; ~110-120 s/case; see the IMPACT
# for the stock binary — see the IMPACT section). --impact-preset selects # section). --impact-iters / --impact-samples override the reference map's
# another published preset (e.g. Elastix_IMPACT_Jacobian, Generic_Rigid); # per-level budget (one iters value for all levels)
# --impact-iters / --impact-samples are no-ops for the default preset's
# per-resolution matrix
python scripts/bench_registration.py --pipelines impact --n-cases 10 \ python scripts/bench_registration.py --pipelines impact --n-cases 10 \
--workers 3 --threads 24 --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`), `--threads T` (also the OpenMP thread count for `reg_aladin`),
`--simpleitk-mode {benchmark,production}`, `--simpleitk-mode {benchmark,production}`,
`--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`, `--synthmorph-model {rigid,affine}`, `--synthmorph-gpu <int|-1>`,
`--synthmorph-extent <int>`, `--impact-preset <name>` (default `--synthmorph-extent <int>`, `--impact-gpu <int|-1>`, `--impact-iters N`
`Elastix_IMPACT_Static`), `--impact-gpu <int|-1>`, `--impact-iters N` (one MaximumNumberOfIterations value for all levels; default the reference
(preset max_iterations override; no-op for matrix presets), map's 400/300/200/200), `--impact-samples N` (NumberOfSpatialSamples; default
`--impact-samples N` (preset spatial_samples override), 2000),
`--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below). `--out-dir <dir>`, `--qa-images` / `--qa-out <dir>` (visual QA, below).
## Visual QA images (checkerboard + color overlay) ## Visual QA images (checkerboard + color overlay)
@ -395,9 +398,9 @@ available pipelines with:
```bash ```bash
source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313
python scripts/bench_registration.py \ 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 \ 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.) (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 inverts it on the M6 cases (8/9 `auto_inverted: True`) to the contract
transform (T1c→CT). 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 The `impact` pipeline runs the
[impact-reg-konfai](https://pypi.org/project/impact-reg-konfai/) package: [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact)
[IMPACT](https://arxiv.org/abs/2503.24121) registration presets from the image: Elastix 5.2 with the [IMPACT](https://arxiv.org/abs/2503.24121)
[`VBoussot/ImpactReg`](https://huggingface.co/VBoussot/ImpactReg) Hugging Face *semantic* similarity metric ([vboussot/ImpactLoss](https://github.com/vboussot/ImpactLoss))
repo, driven through the [KonfAI](https://github.com/fideus-labs/KonfAI) as a native Elastix metric. The parameter map is the MIND + TotalSegmentator
runtime — no Docker. The default preset is **`Elastix_IMPACT_Static`**: **static** feature combination of the `MR_CT_TS` reference map in the
Elastix 5.2 with the IMPACT *semantic* similarity metric — a per-level [`VBoussot/ImpactReg`](https://huggingface.co/VBoussot/ImpactReg) repo, with a
combination of the **MIND** handcrafted descriptor (`MIND/R1D2_3D.pt`, L1, **rigid 6-DOF `EulerTransform`** in place of the reference's RecursiveBSpline
subset 32) and **TotalSegmentator** `TS/M730.pt` (deepest layer, Dice, subset (its bending-energy penalty is zero for a rigid transform, so it is dropped) —
64), layer weights 0.2/0.8 → 0.7/0.3 from coarse to fine, **Static** mode i.e. the reference metric/budget, rigid-only, on the GPU.
(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.
Prerequisites (check with `RP.availability()`; a missing piece raises Prerequisites (check with `RP.availability()`; a missing piece raises
`PipelineUnavailable` with the exact command): `PipelineUnavailable` with the exact command):
```bash ```bash
# in the longitudinal-313 env (Python >= 3.11): # 1. Docker (with NVIDIA container toolkit so --gpus works)
pip install impact-reg-konfai # CLI + KonfAI runtime (pulls h5py etc.) docker --version
bash scripts/setup_impact_konfai.sh # one-time, this host only: the elastix shim
# 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 ### How the adapter drives the image
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`).
### Why the elastix shim (this host only) 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
The stock binary needs **glibc ≥ 2.35** and links **LibTorch 2.8.0**; this evaluates them where they map inside the moving buffer. With the tight T1c
host (Ubuntu 20.04, glibc 2.31, env torch 2.14) satisfies neither, so head crop (~220×220×150 mm) as fixed, the samples map mostly inside the
`scripts/setup_impact_konfai.sh` builds a self-contained shim at large CT head-neck FOV (~300×300×251 mm), so the
`~/.cache/konfai/elastix-impact-shim`: the glibc 2.35 loader + libs extracted `RequiredRatioOfValidSamples 0.1` gate (build default 0.25) rarely fires.
from the `niftyreg:2.1.1-ubuntu` image (docker is needed once for that copy 2. **Parameter map** (the rigid `MR_CT_TS` reference config):
only), the self-contained LibTorch 2.8.0 cu128 from download.pytorch.org `(Metric "Impact" "AdvancedMattesMutualInformation")` with weights
(~3.8 GB zip), and a wrapper `bin/elastix` that runs the stock binary through 1.0 / 0.4; per level (coarse → fine) the **MIND** handcrafted descriptor
the shim loader with the LibTorch 2.8.0 library path (stripping the inherited (`MIND/R1D2_3D.pt`, L1, subset 32, layer mask `"1"`) combined with
`LD_LIBRARY_PATH`, so the env's torch 2.14 cannot leak in). The pipeline **TotalSegmentator** `TS/M852.pt` (levels 0-1) / `TS/M850.pt` (levels
selects the shim via `KONFAI_ELASTIX_DIR` (`_impact_elastix_root`); a host 2-3, Dice, subset 64, deepest layers), layer weights 0.2/0.8, 0.3/0.7,
new enough to run the stock binary directly needs no shim at all. 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
### How the adapter drives the CLI image-resampling pyramid (all-1 schedules — the scale is carried by
ImpactVoxelSize, per the IMPACT docs); `RandomCoordinate` sampler; ASGD
1. **Full volumes, fixed=T1c (MR), moving=CT** — the ImpactReg **MR_CT** with per-level MaximumNumberOfIterations 400/300/200/200 and
convention, which coincides with the shared rule (the T1c is the fixed NumberOfSpatialSamples 2000; `EulerTransform` +
image): the preset's output is a displacement field on the **fixed** grid `AutomaticTransformInitialization` (CenterOfGravity) +
mapping fixed→moving, so the field is the contract transform T1c→CT `AutomaticScalesEstimation`. A copy of the exact map is dumped to
**directly** (displacement fields cannot be inverted with `GetInverse`, `data/m6reg/impact/paramdump/` per run for inspection.
which a fixed=CT direction would require). It also keeps 3. **GPU auto-spread** across the host GPUs by pid (`docker run --gpus all`,
the metric's sample-validity gate healthy: samples drawn from the tight `ImpactGPU` = the worker's device; `--impact-gpu -1` forces CPU and turns
T1c head crop (~220×220×150 mm) map mostly inside the large CT head-neck mixed precision off).
FOV (~300×300×251 mm), unlike fixed=CT where only ~14% were in-buffer at 4. **Read-back** — the transform comes from
the auto-init. `Out/TransformParameters.0.itk.txt` (an `Euler3DTransform` in the world
2. **M6 robustness overrides** — the preset's map is tuned for well-aligned physical frame; the full combination transform, with the auto-init folded
pairs, so the adapter forwards raw elastix lines through the package's in by elastix) and is direction-checked by `fix_orientation` against the
`parameter_overrides` knob: `AutomaticTransformInitialization=true` + skull ring (with fixed=T1c / moving=CT the gate keeps the file's transform,
`AutomaticTransformInitializationMethod=CenterOfGravity` (the preset's `auto_inverted: False`).
GeometricalCenter start leaves the intensity centroids ~50-100 mm apart on 5. **Intensity normalization happens inside the metric / models** (no
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 `<out>/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
pipeline-side normalization): the MIND descriptor normalizes each 5-voxel pipeline-side normalization): the MIND descriptor normalizes each 5-voxel
neighbourhood to [0, 1] locally (modality-agnostic by construction); the neighbourhood to [0, 1] locally (modality-agnostic by construction); the
segmentation models (TS/MRSeg) are standardized to zero mean / unit segmentation models are standardized to zero mean / unit variance per the
variance per the model card. Pre-normalizing the volumes would be an model card. Pre-normalizing the volumes would be an exact no-op for the
exact no-op for the z-scored models and would only corrupt the HU scale z-scored models and would only corrupt the HU scale used by MI scoring and
used by MI scoring and the QA brain-window images. the QA brain-window images.
## Notes / caveats ## Notes / caveats
@ -598,22 +582,20 @@ new enough to run the stock binary directly needs no shim at all.
dataset actually ships with. dataset actually ships with.
- `synthmorph` is deterministic (one fixed pretrained model, no random - `synthmorph` is deterministic (one fixed pretrained model, no random
subsampling), so its per-case numbers are stable run to run — unlike ANTs / 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 IMPACT. Its low MI (median 0.277 with the default CT clip; 0.06–0.11 on the
genuine cross-modality / tilt weakness, not variance. repositioned cases) is a genuine cross-modality / tilt weakness, not
- `fix_orientation` (elastix/ants/niftyreg/synthmorph) chooses between a variance.
transform and its inverse by skull-ring concentricity — a correctness gate - `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 that makes the result independent of each toolkit's transform-direction
convention, so the fixed/moving swap is absorbed by the gate rather than 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 re-derived per toolkit. Under the shared fixed=T1c / moving=CT rule the raw
transforms come out in the contract direction for elastix (combination transforms come out in the contract direction for elastix (combination
transform = fixed→moving), antspyx (ANTs forward transform) and synthmorph transform = fixed→moving), impact (its `TransformParameters.0.itk.txt`,
(the `-T` LTA, model estimates first→second) — completed runs report the full combination Euler transform) and antspyx (ANTs forward transform)
`auto_inverted: False` for those; `niftyreg`'s `-aff` file comes out the and synthmorph (the `-T` LTA, model estimates first→second) — completed
other way and the gate inverts it (8/9 `auto_inverted: True`). `impact` runs report `auto_inverted: False` for those; `niftyreg`'s `-aff` file
returns a displacement field (no `GetInverse` in ITK), so it uses comes out the other way and the gate inverts it (8/9 `auto_inverted: True`).
`_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`).
- **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU): - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU):
- **Cross-modality / tilt weakness.** The `rigid` model is trained on brain - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain
*MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly *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). path).
- **GPU sharing.** The container runs as the host uid with - **GPU sharing.** The container runs as the host uid with
`TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU. `TF_FORCE_GPU_ALLOW_GROWTH` so concurrent benchmark workers share a GPU.
- **~25 s/case** (GPU: one pretrained forward pass), between SimpleITK - **~28 s/case** (GPU: one pretrained forward pass), between SimpleITK
(~10 s) / ANTs (~37 s) and IMPACT (~216 s, deformable preset). (~11 s) / ANTs (~17 s) and IMPACT (~115 s, rigid docker config).
- **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`, - **NiftyReg specifics** (conda-forge `niftyreg` 1.5.69, `reg_aladin`,
symmetric block matching + TLS, CPU/OpenMP): symmetric block matching + TLS, CPU/OpenMP):
- **Different algorithmic family.** No histogram MI objective: 4-voxel - **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 - Benchmark budget: `MaximumNumberOfIterations` 2000 per resolution
(ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their (ASGD 2000/2000/2000 over 3 resolutions) with the sampler/metric at their
defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32). defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32).
- **IMPACT specifics** (`impact-reg-konfai` 1.8.6, KonfAI runtime, - **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2 +
elastix-impact 5.2, GPU; preset `Elastix_IMPACT_Static`): 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 - **Runs on M6 with the T1c-fixed direction + relaxed sample-validity
gate.** Full volumes, no mask, no crop / pad, fixed=T1c/moving=CT (the gate.** Full volumes, no mask, no crop / pad, fixed=T1c/moving=CT (the
shared rule; the ImpactReg MR_CT convention coincides with it): the shared rule): the metric's samples (drawn from the tight T1c head crop)
preset's output field is on the fixed grid, so that direction is the map mostly inside the large CT head-neck FOV, so the gate rarely fires.
contract transform (T1c→CT) directly, and the metric's samples (drawn The map still sets `RequiredRatioOfValidSamples 0.1` (build default
from the tight T1c crop) map mostly inside the large CT FOV. The package's generated map still defaults to 0.25) as a safety override — a genuinely diverging transform, <10%
`RequiredRatioOfValidSamples 0.25`, so the adapter overrides it to 0.1 in-buffer, still aborts. The combination Euler transform comes out in
(a genuinely diverging transform, <10% in-buffer, still aborts). In the the contract direction (completed runs report `auto_inverted: False`).
old Docker era (fixed=CT) the gate aborted the no-mask M6 pairs until - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity)
the same 0.1 relaxation. The transform contract holds; `_orient_impact` aligns the intensity centroids (no manual translation); the reference
decides the direction as usual (completed runs report map's GeometricalCenter start is too crude for the head+neck FOV
`auto_inverted: False` — the MR_CT field is already the contract). mismatch.
- **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) is - **Rigid, like every other pipeline.** The reference map's
applied as an override (the preset's GeometricalCenter start is too RecursiveBSpline + bending-energy regularizer is replaced by a 6-DOF
crude for the head+neck FOV mismatch); no manual translation is used. `EulerTransform` (the bending penalty is zero for a rigid transform).
- **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.
- **History: single-init convergence.** The semantic landscape has more - **History: single-init convergence.** The semantic landscape has more
local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer local minima than Mattes-MI. The T1c-fixed direction fixes most of the
nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0; hard cases (the repositioned `m6_22IRRKU6` / `m6_23XIMO4M` now register
the pre-mask head-box flow got it to mid 8.3 / std 13.2). The no-mask as well as the MI tools do), but the single CenterOfGravity start on the
MRSegmentator runs showed the same pattern (0.635 MI on `m6_22IRRKU6` ~30° tilt pair `m6_2BBN5XA4` remains a weakness — a multi-start tilt
at 2000/5000 but mid 26.1 / std 28.6 there; accepted on `m6_2BBN5XA4`, grid (as production SimpleITK uses) would be the next step.
mid 15.1 / std 10.7). The single CenterOfGravity start on the ~30° - **Stochastic.** The config is **Static** mode: it re-samples spatial
tilt pairs remains IMPACT's limiting factor — a multi-start (or coordinates every iteration, but the feature maps are deterministic
`--tta` flips through the preset app) would be the next step. given the pose (see Stochasticity). The older Jacobian-mode configs
- **Stochastic.** The Jacobian preset (`Elastix_IMPACT_Jacobian`) (docker-era experiments, archived in `results/bench_registration/impact_*.json`)
resamples patches through the network each iteration, so results vary resampled patches through the network each iteration, so their
run to run; the default Static preset re-samples coordinates each per-case numbers are one sample from a distribution.
iteration but the feature maps are deterministic given the pose (see - **Models are under the data root, not the repo.** The TorchScript
Stochasticity). feature models live in `data/impact_models/` (MIND ~15 KB + TS ~200 MB
- **The elastix shim (this host).** The stock binary needs glibc ≥ 2.35 each, from HF `VBoussot/impact-torchscript-models`); the pipeline
+ LibTorch 2.8.0; the host (glibc 2.31, torch 2.14) runs it through raises with the exact `curl` commands if one is missing.
`~/.cache/konfai/elastix-impact-shim` (built by - **CPU fallback:** `--impact-gpu -1` (`ImpactGPU -1`, mixed precision
`scripts/setup_impact_konfai.sh`, selected via `KONFAI_ELASTIX_DIR`). off — much slower).
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).

132
scripts/bench_registration.py
View file

@ -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: (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 - mi_ct_mr : Mattes MI (nats), warped-CT ~ MR over the MR head
(the reported alignment score). (the reported alignment score, higher = better).
- tv_in_ct_frac : fraction of the (dilated) warped GTV that falls inside - impact_loss : IMPACT semantic loss (vboussot/ImpactLoss; TotalSeg
the warped-CT head tissue (the slab must straddle it). M730 deep features, weighted L1), warped-CT ~ MR over
the full co-registered volumes (lower = better).
- tv_nvox : warped GTV volume (>0.5 voxels). - tv_nvox : warped GTV volume (>0.5 voxels).
- skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR - skull_ring_gap_mm : concentricity of the warped-CT skull ring vs the MR
head silhouette (mid/low level + angular std). 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 pipelines that are importable there. elastix runs in the main env through the
ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively
(``conda install -c conda-forge niftyreg`` provides ``reg_aladin``); (``conda install -c conda-forge niftyreg`` provides ``reg_aladin``);
``impact`` runs the ``impact-reg-konfai`` package (pip; GPU; one-time elastix ``impact`` runs the ``vboussot/elastix_impact`` Docker image (GPU; the rigid
shim setup on hosts too old for the stock binary — see MIND + TotalSegmentator static config, needs Docker + the pulled image + the
REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker image model assets — see REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker
(``freesurfer/synthmorph``; needs Docker + the pulled image). The `--merge` image (``freesurfer/synthmorph``; needs Docker + the pulled image). The
`--merge`
mode then combines the per-pipeline JSON files (accumulated across envs) into mode then combines the per-pipeline JSON files (accumulated across envs) into
one summary + plot. one summary + plot.
@ -99,6 +101,54 @@ def _m6():
return RP._load_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): def _clip_image(img, lo, hi):
"""Clip image intensities to [lo, hi] (geometry preserved).""" """Clip image intensities to [lo, hi] (geometry preserved)."""
a = np.clip(sitk.GetArrayFromImage(img).astype("float32"), lo, hi) 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_ct = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
warped_tv = sitk.GetArrayFromImage(sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0.0)) warped_tv = sitk.GetArrayFromImage(sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0.0))
gap = m6._skull_ring_gap(warped_ct, ha) gap = m6._skull_ring_gap(warped_ct, ha)
# achieved-metric alignment score: Mattes MI between the CT inverse-transformed # achieved-metric alignment scores: Mattes MI between the CT
# into MR coordinates (warped_ct) and the original MR (ha), over the MR head. # 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) ha_f = ha.astype("float32", copy=False)
mi = _mi_masked(warped_ct, ha_f, head_mask_from_image(ha_f)) 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] posw = warped_ct[warped_ct > 0]
tissue_w = (warped_ct > max(2.0, float(np.percentile(posw, 2)))) if posw.size \ tissue_w = (warped_ct > max(2.0, float(np.percentile(posw, 2)))) if posw.size \
else np.zeros_like(warped_ct, bool) else np.zeros_like(warped_ct, bool)
@ -171,6 +223,7 @@ def evaluate(ct, t1c, tv, T):
return { return {
"skull_ring_gap_mm": gap, "skull_ring_gap_mm": gap,
"mi_ct_mr": round(mi, 4) if mi is not None else None, "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_nvox": nvox,
"tv_in_ct_frac": round(tv_in_ct_frac, 3), "tv_in_ct_frac": round(tv_in_ct_frac, 3),
"alignment": "concentric" if ok_gap else "not_concentric", "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, gpu=args.synthmorph_gpu, threads=args.threads,
extent=args.synthmorph_extent) extent=args.synthmorph_extent)
if pipeline == "impact": if pipeline == "impact":
return RP.register_impact(ct, t1c, preset=args.impact_preset, gpu=args.impact_gpu, return RP.register_impact(ct, t1c, gpu=args.impact_gpu,
iters=args.impact_iters, samples=args.impact_samples) iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS,
samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES)
if pipeline == "niftyreg": if pipeline == "niftyreg":
return RP.register_niftyreg(ct, t1c, threads=args.threads, return RP.register_niftyreg(ct, t1c, threads=args.threads,
backend=args.niftyreg_backend, backend=args.niftyreg_backend,
@ -268,6 +322,7 @@ def _run_case(pipeline, case, args):
out = {"key": key, "subject": case["subject"], "date": case["date"]} out = {"key": key, "subject": case["subject"], "date": case["date"]}
if not ct_complete(ct): if not ct_complete(ct):
out.update({"skull_ring_gap_mm": None, "tv_nvox": 0, "tv_in_ct_frac": None, out.update({"skull_ring_gap_mm": None, "tv_nvox": 0, "tv_in_ct_frac": None,
"impact_loss": None,
"alignment": "ct_incomplete", "accepted": False, "alignment": "ct_incomplete", "accepted": False,
"note": "CT head incomplete (data issue, not a pipeline failure)"}) "note": "CT head incomplete (data issue, not a pipeline failure)"})
return out return out
@ -313,7 +368,7 @@ def run_pipeline(pipeline, cases, args):
r = results[-1] r = results[-1]
reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT" 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')} " 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 except Exception as e: # noqa
errors.append({"key": c["key"], "error": repr(e)}) errors.append({"key": c["key"], "error": repr(e)})
print(f" {c['key']}: ERROR {e!r}", flush=True) 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_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")] 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] 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] 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 {} first = (registered or results)[0] if (registered or results) else {}
row = { row = {
@ -341,7 +396,8 @@ def run_pipeline(pipeline, cases, args):
"median": round(float(np.median(gaps_mid)), 1)} if gaps_mid else None, "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), "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, "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), "mi_ct_mr": {"mean": round(float(np.mean(mifr)), 4),
"median": round(float(np.median(mifr)), 4)} if mifr else None, "median": round(float(np.median(mifr)), 4)} if mifr else None,
"cases": results, "errors": errors, "cases": results, "errors": errors,
@ -368,24 +424,24 @@ def _fmt(v, spec="{}", default="-"):
def print_table(rows): def print_table(rows):
order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in 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} " 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(hdr)
print("-" * 92) print("-" * 96)
for p in order: for p in order:
r = next(x for x in rows if x["pipeline"] == p) r = next(x for x in rows if x["pipeline"] == p)
if not r.get("available"): if not r.get("available"):
print(f"{p:10s} {('no'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} " 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 continue
rt = (r.get("reg_wall_s") or {}).get("median") 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") 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} " 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"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {_fmt(rt, '{:.1f}'):>7s} "
f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}") f"{_fmt(il, '{:.3f}'):>10s} {_fmt(mi, '{:.3f}'):>7s}")
print("=" * 92) print("=" * 96)
def write_plot(rows, out_png): def write_plot(rows, out_png):
@ -402,13 +458,15 @@ def write_plot(rows, out_png):
return return
labels = list(ok_rows) labels = list(ok_rows)
regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") or 0 for p in labels] regt = [(ok_rows[p].get("reg_wall_s") or {}).get("median") or 0 for p in labels]
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] 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)) fig, axes = plt.subplots(1, 3, figsize=(14, 4))
colors = ["#4878CF", "#6ACC65"] colors = ["#6ACC65", "#F5A623", "#4878CF"]
for ax, (title, vals, fmt, ylab) in zip(axes, [ for ax, (title, vals, fmt, ylab) in zip(axes, [
("Registration time (s, median)", regt, "{:.0f}", "s"), ("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")]): ("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): for b, v in zip(bars, vals):
ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()), ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()),
ha="center", va="bottom", fontsize=8) ha="center", va="bottom", fontsize=8)
@ -438,7 +496,9 @@ def merge(args):
summary = { summary = {
"timestamp": datetime.now().isoformat(timespec="seconds"), "timestamp": datetime.now().isoformat(timespec="seconds"),
"protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI " "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 " "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 " "gap / acceptance QA; best-of-T vs T^-1 orientation gate for the "
"elastix-style toolkits"), "elastix-style toolkits"),
@ -452,7 +512,7 @@ def merge(args):
print("summary:", sf) print("summary:", sf)
write_plot(rows, os.path.join(out_dir, "bench_registration.png")) write_plot(rows, os.path.join(out_dir, "bench_registration.png"))
# markdown table for the report # 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]: for p in [x["pipeline"] for x in rows]:
r = next(x for x in rows if x["pipeline"] == p) 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"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | "
f"{r['n_error']} | " f"{r['n_error']} | "
f"{(r.get('reg_wall_s') or {}).get('median')} | " 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')} |") f"{(r.get('mi_ct_mr') or {}).get('median')} |")
mf = os.path.join(out_dir, "summary.md") mf = os.path.join(out_dir, "summary.md")
with open(mf, "w") as f: with open(mf, "w") as f:
@ -496,24 +556,16 @@ def parse_args():
ap.add_argument("--synthmorph-extent", type=int, default=256, ap.add_argument("--synthmorph-extent", type=int, default=256,
help="SynthMorph isotropic registration-space extent in mm " help="SynthMorph isotropic registration-space extent in mm "
"(192 or 256; default 256)") "(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, ap.add_argument("--impact-gpu", type=int, default=None,
help="GPU device index for IMPACT (default: auto-spread over " help="GPU device index for IMPACT (default: auto-spread over "
"the host GPUs by pid; -1 forces CPU)") "the host GPUs by pid; -1 forces CPU)")
ap.add_argument("--impact-iters", type=int, default=None, ap.add_argument("--impact-iters", type=int, default=None,
help="preset max_iterations override, one value for all " help="IMPACT MaximumNumberOfIterations, one value for all "
"resolutions (only effective for presets without a " "levels (default: the reference map's per-level budget "
"per-resolution model matrix; default: the preset's " f"{RP.IMPACT_DEFAULT_ITERS})")
"own budget, 400/300/200/200 for Elastix_IMPACT_Static)")
ap.add_argument("--impact-samples", type=int, default=None, ap.add_argument("--impact-samples", type=int, default=None,
help="preset spatial_samples override (NumberOfSpatialSamples; " help="IMPACT NumberOfSpatialSamples per iteration (default: "
"default: the preset's own value, 2000)") f"{RP.IMPACT_DEFAULT_SAMPLES})")
ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"], ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"],
default="auto", default="auto",
help="niftyreg reg_aladin build: docker = official v2.1.1 in the " help="niftyreg reg_aladin build: docker = official v2.1.1 in the "

110
scripts/setup_impact_konfai.sh
View file

@ -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"

152
src/impact_loss.py Normal file
View file

@ -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

580
src/registration_pipelines.py
View file

@ -31,8 +31,6 @@ Shared rules for every pipeline:
- **The T1c (MR) is always the fixed image** (moving = CT). When both inputs - **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) 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`. 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 brain-window clip for the pipelines that benefit from it**: per the
CT-clipping A/B (see REGISTRATION_BENCHMARK.md), ``niftyreg`` and CT-clipping A/B (see REGISTRATION_BENCHMARK.md), ``niftyreg`` and
``synthmorph`` clip the CT (the moving image) to [0, 80] HU before ``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`` + - ``register_elastix`` : ITK-Elastix Python API (pip ``itk`` +
``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD). ``itk-elastix``, elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
``AutomaticTransformInitialization`` (CenterOfGravity) start. ``AutomaticTransformInitialization`` (CenterOfGravity) start.
- ``register_impact`` : IMPACT (semantic-similarity metric) via the - ``register_impact`` : IMPACT (semantic-similarity metric) inside the
``impact-reg-konfai`` package (pip; the KonfAI runtime drives the ``vboussot/elastix_impact`` Docker image (Elastix 5.2 + TorchScript feature
elastix-impact 5.2 binary -- no Docker). Default preset models as a native Elastix metric). Rigid 6-DOF ``EulerTransform`` driven
``Elastix_IMPACT_Static``: the published successor of the ImpactReg by the MIND + TotalSegmentator static feature combination of the ImpactReg
``MR_CT_TS`` reference map (elastix B-spline + per-level MIND + ``MR_CT_TS`` reference map (per-level models / scales / layer weights) + a
TotalSegmentator M730, Static mode). Registered fixed=T1c/moving=CT (the Mattes-MI secondary metric (w 0.4). GPU. Full volumes, no mask (shared
ImpactReg MR_CT convention) so the preset's output field (on the fixed rules above). Needs Docker + the image + the model assets under the cache.
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_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` -- - ``register_niftyreg`` : NiftyReg (KCL-BMEIS suite) ``reg_aladin`` --
symmetric block-matching + trimmed least squares (Modat et al., J. Med. symmetric block-matching + trimmed least squares (Modat et al., J. Med.
Img. 2014), NCC-based block similarity (cross-modality tolerant). Rigid 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 ``T`` and ``T^-1`` and keeping whichever gives the smaller skull-ring
concentricity gap (a pure correctness gate -- a wrongly-oriented rigid transform 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 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 applied to elastix / antspyx / synthmorph / impact, not to the trusted
SimpleITK path. ``impact`` returns a displacement field (which has no production SimpleITK path.
``GetInverse``), so it uses :func:`_orient_impact` -- the same skull-ring
decision against a sign-flipped (first-order reverse) proxy field.
""" """
import os import os
import re 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" IMPACT_IMAGE = "vboussot/elastix_impact"
# default preset from the published ImpactReg database (VBoussot/ImpactReg, IMPACT_ELASTIX_BIN = "/lib/elastix-install/bin/elastix"
# resolved by the CLI): elastix B-spline driven by the IMPACT semantic metric IMPACT_MODEL_CACHE = "data/impact_models"
# -- per-level MIND (L1) + TotalSegmentator M730 (deepest layer, Dice), Static IMPACT_HF_MODEL_URL = (
# mode (whole-image features, no backprop), 4 scales, 400/300/200/200 "https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/{rel}"
# 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",
) )
# The stock elastix-impact 1.0.0 binary (vboussot/ImpactElastix) needs # The handcrafted MIND descriptor (whole-image features, layer mask "1", L1)
# glibc >= 2.35 and links LibTorch 2.8.0; a host that cannot run it directly # combined per level with a TotalSegmentator feature model (deepest layers,
# (this one: glibc 2.31 + env torch 2.14) builds a self-contained shim with # Dice) -- the MIND+TS static combination of the ImpactReg ``MR_CT_TS``
# scripts/setup_impact_konfai.sh, selected via KONFAI_ELASTIX_DIR. # reference parameter map, coarse -> fine. Per level: (TS model, voxel-size
IMPACT_ELASTIX_SHIM = os.path.join( # scale in mm, MIND layer mask, TS layer mask, MIND weight, TS weight).
os.path.expanduser("~"), ".cache", "konfai", "elastix-impact-shim") IMPACT_MIND_MODEL = "MIND/R1D2_3D.pt"
IMPACT_ELASTIX_STOCK = os.path.join( IMPACT_LEVELS = (
os.path.expanduser("~"), ".cache", "konfai", "elastix-impact") ("TS/M852.pt", 6.0, "1", "0000001", 0.2, 0.8),
# orientation-gate threshold: a candidate DVF is kept as-is when its ("TS/M852.pt", 3.0, "1", "0000001", 0.3, 0.7),
# skull-ring gap score (std_mm + mid_mm) is below this. A misdirected field ("TS/M850.pt", 2.0, "1", "00000001", 0.6, 0.4),
# lands the head on the far side of the box (gap ~100-300 mm), so the margin ("TS/M850.pt", 2.0, "1", "00000001", 0.7, 0.3),
# is wide; see _orient_impact. )
IMPACT_GAP_ACCEPT = 60.0 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(): def _impact_docker_ok():
"""Return (ok, note): is the impact-reg-konfai CLI installed?""" """Return (ok, note): is the Docker CLI + the IMPACT image available?"""
if shutil.which(IMPACT_CLI) is None: if shutil.which("docker") is None:
return False, f"{IMPACT_CLI} not on PATH (pip install impact-reg-konfai)" return False, "docker CLI not found on PATH"
try: try:
import importlib.metadata r = subprocess.run(["docker", "image", "inspect", IMPACT_IMAGE],
ver = importlib.metadata.version("impact-reg-konfai") capture_output=True, text=True, timeout=30)
except Exception: # noqa if r.returncode != 0:
ver = "?" return False, (f"docker image {IMPACT_IMAGE} not pulled "
return True, f"{IMPACT_CLI} {ver}" f"(run: docker pull {IMPACT_IMAGE})")
except Exception as e: # noqa
return False, f"docker image check failed ({type(e).__name__}: {e})"
_ELASTIX_ROOT_CACHE = None return True, IMPACT_IMAGE
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
_NGPUS_CACHE = None _NGPUS_CACHE = None
@ -773,156 +719,196 @@ def _n_gpus():
return _NGPUS_CACHE return _NGPUS_CACHE
def _impact_dvf_from_h5(path): def _impact_models_host():
"""An ITK-form HDF5 transform (the ``:itktransform`` form the KonfAI """Resolve the IMPACT feature models to files under the model cache, or
presets write) read with h5py -- the SimpleITK build in this env has no raise :class:`PipelineUnavailable` with the exact download commands for
ITK HDF5 transform reader. Returns ``(DisplacementFieldTransform, the missing ones. Returns the de-duplicated paths relative to the cache
dvf_image)``; the field is in mm, channel-last ``(z, y, x, 3)``.""" (which the pipeline bind-mounts at ``/Data/Models`` inside the
import h5py container)."""
with h5py.File(str(path), "r") as f: cache = d(IMPACT_MODEL_CACHE)
grp = f["TransformGroup/0"] rels = [IMPACT_MIND_MODEL] + [rel for rel, *_ in IMPACT_LEVELS]
ttype = grp["TransformType"][()][0] missing = [r for r in dict.fromkeys(rels)
if isinstance(ttype, bytes): if not os.path.exists(os.path.join(cache, r))]
ttype = ttype.decode() if missing:
fixed = grp["TransformFixedParameters"][()] cmds = "\n".join(
params = grp["TransformParameters"][()] f" mkdir -p {os.path.join(cache, os.path.dirname(r))} && "
if not str(ttype).startswith("DisplacementFieldTransform"): f"curl -L -o {os.path.join(cache, r)} {IMPACT_HF_MODEL_URL.format(rel=r)}"
raise RuntimeError(f"impact: unsupported h5 transform type {ttype!r}") for r in missing)
size = fixed[:3][::-1].astype(int) # stored (x, y, z) -> array (z, y, x) raise PipelineUnavailable(
arr = params.reshape(int(size[0]), int(size[1]), int(size[2]), 3) f"impact: model(s) not found under {cache}.\nDownload once with:\n"
img = sitk.GetImageFromArray(arr) f"{cmds}\n(or `sh Data/Models/download_models.sh` from the "
img.SetOrigin(tuple(float(v) for v in fixed[3:6])) f"ImpactLoss repo).")
img.SetSpacing(tuple(float(v) for v in fixed[6:9])) return list(dict.fromkeys(rels))
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_read_transform(path): def _impact_param_map(gpu, iters, samples, voxel_sizes):
"""The preset's single-case transform output as """Build the Elastix-5.2 parameter map for a RIGID (Euler) IMPACT
``(sitk.Transform, dvf_image_or_None)``. registration: the MIND + TotalSegmentator static feature combination of
the ImpactReg ``MR_CT_TS`` reference map (one model pair per level,
KonfAI presets write the registration in whatever form the preset coarse -> fine, see :data:`IMPACT_LEVELS`) + a secondary Mattes-MI metric
declares: an ITK-form HDF5 DisplacementFieldTransform (``Transform.h5``) (:data:`IMPACT_MI_WEIGHT`). The image pyramid is NOT resampled (all-1
by default, or an ITK transform file (``.itk.txt`` / ``.tfm``) or a vector schedules); the scale comes from the per-level ImpactVoxelSize (as the
image (``.mha``) for other presets.""" IMPACT docs recommend). ``iters`` is a per-level tuple or a single int.
from pathlib import Path """
p = Path(path) n = len(voxel_sizes)
if p.name.endswith(".itk.txt") or p.suffix in (".h5", ".tfm"): if isinstance(iters, (int,)):
if p.suffix == ".h5": iters_l = str(iters)
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)
else: else:
try: iters_l = " ".join(f"{i:g}" for i in iters)
rev = cand.GetInverse() mixed = "true" if gpu >= 0 else "false"
except Exception: # noqa: parametric transform with no inverse L = ["// IMPACT semantic-similarity metric (vboussot/ImpactLoss),"]
return cand, False L.append("// MIND + TotalSegmentator static combination of the")
s1 = _gap_score(_skull_gap(_warp_ct(ct, t1c, rev), ha)) L.append("// VBoussot/ImpactReg MR_CT_TS reference map (EulerTransform")
if s1 < s0: L.append("// instead of the reference RecursiveBSpline -- rigid)")
return rev, True L.append('(Metric "Impact" "AdvancedMattesMutualInformation")')
return cand, False 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): keep_scratch=False):
"""IMPACT (semantic-similarity metric) CT->T1c via the """IMPACT rigid CT->T1c (fixed=T1c, moving=CT) in the Docker image
``impact-reg-konfai`` package (pip; KonfAI runtime driving the ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric).
elastix-impact 5.2 binary; no Docker).
IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic* IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic*
one. The default preset ``Elastix_IMPACT_Static`` is the published one: pretrained TorchScript feature models are run on the images and the
successor of the old ImpactReg ``MR_CT_TS`` reference parameter map: an registration minimizes the distance between the extracted feature maps.
elastix **B-spline** (4 scales, FinalGridSpacing 14 mm, 400/300/200/200 The metric here is the MIND + TotalSegmentator **static** combination of
iterations, 2000 spatial samples, bending-energy regularizer) driven by the ImpactReg ``MR_CT_TS`` reference map (``MR_CT_TS/ParameterMap_MRI_TS``
the IMPACT metric -- a per-level **MIND** (L1) + **TotalSegmentator M730** in ``VBoussot/ImpactReg``): per level the MIND handcrafted descriptor
(deepest layer, Dice) combination in **Static** mode (whole-image (``MIND/R1D2_3D.pt``, L1) is combined with a TotalSegmentator deepest-layer
features, no backprop), plus a Mattes-MI secondary metric (w 0.4). Other feature (``TS/M852.pt`` coarse, ``TS/M850.pt`` fine, Dice), layer weights
presets (``Elastix_IMPACT_Jacobian``, ``Generic_Rigid``, 0.2/0.8 -> 0.7/0.3, at voxel sizes 6/3/2/2 mm (scales with the moving
``Generic_Rigid_BSpline``, ``FireANTs_*``, ``ConvexAdam_*``) are image's spacing), plus a secondary Mattes-MI metric (weight 0.4). The
selectable with ``preset=``; the live database is listed by models standardize intensity inside the metric (zero mean / unit variance
``impact-reg-konfai register --help``. 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** The optimized transform is a rigid 6-DOF ``EulerTransform`` (the 5.x name
convention (fixed=T1c, moving=CT), which coincides with the shared for the classic 3D rigid) with ASGD -- no image resampling pyramid, the
fixed-image rule (the T1c is the fixed image). The preset's output is a scale comes from the per-level ImpactVoxelSize. No manual centroid init:
displacement field on the **fixed** grid mapping fixed->moving, so the elastix's ``AutomaticTransformInitialization`` (CenterOfGravity) aligns
field comes out directly in the contract direction T1c->CT (displacement the intensity centroids.
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.
M6 robustness overrides (raw elastix lines the package applies on top of Like every other pipeline this runs on the **full volumes with
the preset's generated map): CenterOfGravity auto init (the preset's fixed=T1c / moving=CT** (the shared fixed-image rule; no mask, no crop /
GeometricalCenter start leaves the intensity centroids ~50-100 mm apart pad). The metric draws its random spatial samples from the fixed (T1c)
on the head+neck FOV mismatch) and ``RequiredRatioOfValidSamples 0.1`` buffer and evaluates them where they map inside the moving (CT) buffer;
(the build default 0.25 aborts the FOV-mismatched pairs at the init; 0.1 the samples of the tight T1c head crop map mostly inside the large CT
keeps the divergence gate). 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 Budget: the default ``iters=(400, 300, 200, 200)`` per level /
use: the preset app (``VBoussot/ImpactReg``), the TorchScript feature ``samples=2000`` is the reference map's budget (~75-80 s/case on the GPU).
models (``VBoussot/impact-torchscript-models``) and the elastix-impact Pass a single int to ``iters`` to use one budget for every level. A GPU
binary (``vboussot/ImpactElastix``) into ``~/.cache``. The binary needs is used when one is visible (auto-spread over the host GPUs by pid);
glibc >= 2.35 and LibTorch 2.8.0; on a host that cannot run it directly pass ``gpu=-1`` to force CPU (mixed precision is then turned off).
(this one: glibc 2.31 + torch 2.14) ``scripts/setup_impact_konfai.sh`` The result is read back from ``TransformParameters.0.itk.txt`` (an
builds a self-contained shim, selected via ``KONFAI_ELASTIX_DIR``. ``Euler3DTransform`` in the world physical frame).
Budget: ``iters`` / ``samples`` forward to the preset's Preconditions (else :class:`PipelineUnavailable`): the ``docker`` CLI,
``max_iterations`` / ``spatial_samples`` knobs (only effective for the ``vboussot/elastix_impact`` image, and the model assets under the
presets without a per-resolution model matrix; the default preset's cache (``data/impact_models``, auto-created; see the exact download
matrix wins). A GPU is used when one is visible (auto-spread over the commands in the error).
host GPUs by pid); pass ``gpu=-1`` to force CPU (the engine then turns
mixed precision off). The result is read back from
``<out>/P000/<group>.h5`` (an ITK-form DisplacementFieldTransform on the
T1c grid) and direction-checked by :func:`_orient_impact` (skull-ring
gap vs the reverse-proxy field).
""" """
import time import time
ok, note = _impact_cli_ok() ok, note = _impact_docker_ok()
if not ok: if not ok:
raise PipelineUnavailable(f"impact: {note}") raise PipelineUnavailable(f"impact: {note}")
root = _impact_elastix_root() model_rels = _impact_models_host()
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
# fixed=T1c (MR), moving=CT: see the docstring (contract direction comes # full volumes, no mask / crop / pad; the T1c is the fixed image (shared
# out of the preset directly; MR_CT convention = the shared rule) # rule). No manual init -- elastix's AutomaticTransformInitialization
# (CenterOfGravity) aligns the image centers.
fixed, moving = fixed_moving(ct, t1c) fixed, moving = fixed_moving(ct, t1c)
# GPU: auto-spread across the host GPUs by pid, or forced value # 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() ngpu = _n_gpus()
gpu = (os.getpid() % ngpu) if ngpu else -1 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_") tmp = tempfile.mkdtemp(prefix="impact_reg_")
tail = "" datadir = os.path.join(tmp, "Data")
try:
fixed_path = os.path.join(tmp, "Fixed.mha")
moving_path = os.path.join(tmp, "Moving.mha")
outdir = os.path.join(tmp, "Out") outdir = os.path.join(tmp, "Out")
sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), fixed_path, True) os.makedirs(datadir, exist_ok=True)
sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), moving_path, True) os.makedirs(outdir, exist_ok=True)
try:
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, # debug aid: keep a persistent copy of the exact parameter map elastix
"-f", fixed_path, "-m", moving_path, # is about to run, written BEFORE the docker call so an abort leaves it
"-o", outdir, # behind for inspection (per-pid name: parallel workers do not clobber
"--fields-only", # we derive the warp ourselves # each other's dumps)
"--tmp-dir", os.path.join(tmp, "Work")] 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: if gpu >= 0:
cmd += ["--gpu", str(gpu)] cmd += ["--gpus", "all"]
# preset tuning (YAML values, forwarded by the CLI to the preset app): cmd += [
# the M6 robustness overrides, always; the budget knobs, when set "-v", f"{datadir}:/Data",
# (no-op for presets whose per-resolution matrix wins) "-v", f"{outdir}:/Out",
ov = "[" + ",".join(f'"{o}"' for o in IMPACT_PARAM_OVERRIDES) + "]" "-v", f"{d(IMPACT_MODEL_CACHE)}:/Data/Models",
cmd += ["--set", f"parameter_overrides={ov}"] IMPACT_IMAGE, IMPACT_ELASTIX_BIN,
if iters is not None: "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha",
cmd += ["--set", f"max_iterations={int(iters)}"] "-p", "/Data/ParameterMap.txt",
if samples is not None: "-out", "/Out/", "-threads", str(threads),
cmd += ["--set", f"spatial_samples={int(samples)}"] ]
env = os.environ.copy()
env["KONFAI_ELASTIX_DIR"] = root
t0 = time.time() t0 = time.time()
proc = subprocess.run(cmd, capture_output=True, text=True, proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200)
timeout=7200, env=env)
dt = time.time() - t0 dt = time.time() - t0
tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
if proc.returncode != 0: if proc.returncode != 0:
raise RuntimeError(f"impact: {IMPACT_CLI} failed " tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
f"(rc={proc.returncode}):\n{tail}") raise RuntimeError(f"impact: docker elastix failed (rc={proc.returncode}):\n{tail}")
case_dirs = sorted(p for p in os.listdir(outdir) itk_tf = os.path.join(outdir, "TransformParameters.0.itk.txt")
if os.path.isdir(os.path.join(outdir, p))) if not os.path.exists(itk_tf):
if not case_dirs: tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:])
raise RuntimeError(f"impact: no case output under {outdir}:\n{tail}") raise RuntimeError(f"impact: no TransformParameters.0.itk.txt produced:\n{tail}")
case_dir = os.path.join(outdir, case_dirs[0])
entries = [p for p in os.listdir(case_dir) # the optimized Euler (world frame) is the full combination transform
if os.path.isfile(os.path.join(case_dir, p)) # (AutomaticTransformInitialization is folded in by elastix); with
and not p.startswith("Moved")] # fixed=T1c / moving=CT the direction it holds is resolved empirically
if len(entries) != 1: # by fix_orientation against the skull ring.
raise RuntimeError( T_main = sitk.ReadTransform(itk_tf)
f"impact: expected exactly one transform output in {case_dir}, " T, inverted = fix_orientation(ct, t1c, T_main)
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)
finally: finally:
if not keep_scratch: if not keep_scratch:
shutil.rmtree(tmp, ignore_errors=True) shutil.rmtree(tmp, ignore_errors=True)
models = " + ".join(model_rels)
meta = { meta = {
"pipeline": "impact", "pipeline": "impact",
"toolkit": (f"impact-reg-konfai (KonfAI; elastix-impact 5.2, " "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)",
f"preset {preset})"), "metric": (f"IMPACT semantic ({models}, L1+Dice, Static) + "
"metric": ("IMPACT semantic (MIND R1D2_3D L1 + TotalSegmentator M730 " f"Mattes MI (w {IMPACT_MI_WEIGHT})"),
"deepest-layer Dice, per-level weights, Static) + " "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; "
"Mattes MI (w 0.4) + bending energy"),
"transform": (f"Deformable (preset-native B-spline DVF on the T1c grid; "
f"CenterOfGravity auto init; "
f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"),
"convention": "T1c->CT (verified)", "auto_inverted": inverted, "convention": "T1c->CT (verified)", "auto_inverted": inverted,
"fixed": _modality_name(fixed), "moving": _modality_name(moving), "fixed": _modality_name(fixed), "moving": _modality_name(moving),
"seconds_reg": round(dt, 1), "seconds_reg": round(dt, 1),
"preset": preset, "gpu": gpu, "iters": iters, "samples": samples, "docker_image": IMPACT_IMAGE, "model": models, "gpu": gpu,
"parameter_overrides": list(IMPACT_PARAM_OVERRIDES), "iters": list(iters) if not isinstance(iters, int) else iters,
"elastix_root": os.path.basename(root), "samples": samples, "param_dump": param_dump,
} }
if keep_scratch: if keep_scratch:
meta["scratch"] = tmp meta["scratch"] = tmp
@ -1272,21 +1265,18 @@ def availability():
else: else:
out[name] = (True, f"Docker {SYNTHMORPH_IMAGE}") out[name] = (True, f"Docker {SYNTHMORPH_IMAGE}")
elif name == "impact": elif name == "impact":
ok, note = _impact_cli_ok() ok, note = _impact_docker_ok()
if not ok: if not ok:
out[name] = (False, note) out[name] = (False, note)
continue continue
root = _impact_elastix_root() try:
if root is None: rels = _impact_models_host()
out[name] = (False, except PipelineUnavailable as e:
note + "; no runnable elastix-impact binary " out[name] = (False, str(e))
"(needs glibc >= 2.35 + LibTorch 2.8.0) -- run " continue
"scripts/setup_impact_konfai.sh once") out[name] = (True, f"Docker {IMPACT_IMAGE} + {len(rels)} models "
else: f"in {d(IMPACT_MODEL_CACHE)} (rigid MIND+TS "
out[name] = (True, f"static config)")
note + "; elastix-impact under "
f"{os.path.basename(root)}; preset "
f"{IMPACT_DEFAULT_PRESET}")
except Exception as e: # noqa except Exception as e: # noqa
out[name] = (False, repr(e)) out[name] = (False, repr(e))
return out return out