diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index 6974684..565f6a1 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/REGISTRATION_BENCHMARK.md @@ -39,10 +39,13 @@ center-of-mass alignment, simpleitk starts from identity, niftyreg consequence of the no-mask rule: the IMPACT metric only evaluates samples that map inside the moving buffer, so with the CT (head-neck, ~300×300×251 mm) as the fixed image and the tight T1c head crop (~220×220×150 mm) as moving, most -of the CT buffer maps outside the T1c buffer and `impact` aborts on the M6 -pairs ("too many samples map outside moving image buffer", ~14% of samples -in-buffer at the auto-init) — see the IMPACT section. The comparison is on -each toolkit's native optimizer / metric / global-search capability. +of the CT buffer maps outside the T1c buffer (~14% of samples in-buffer at the +auto-init). Elastix's sample-validity gate (build default +`RequiredRatioOfValidSamples 0.25`) used to abort the M6 pairs on that ("too +many samples map outside moving image buffer"); the parameter map now sets +`RequiredRatioOfValidSamples 0.1` so the FOV-mismatched runs proceed — see the +IMPACT section. The comparison is on each toolkit's native optimizer / metric +/ global-search capability. ## The six pipelines @@ -61,7 +64,7 @@ below; `clip_ct=None` disables it). |---|----------|----------------|------|--------|-----------|--------------------| | 1 | `antspyx` | CT / T1c | auto COM | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (main env; stochastic MI subsampling) | | 2 | `elastix` | CT / T1c | auto CenterOfGravity | Mattes MI (32 bins) | `EulerTransform` rigid (5.x), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) | -| 3 | `impact` | **CT / T1c** | auto CenterOfGravity | **IMPACT semantic** (TotalSegmentator `TS/M730.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ⛔ aborts on M6 (no-mask rule; see IMPACT section) | +| 3 | `impact` | **CT / T1c** | auto CenterOfGravity | **IMPACT semantic** (MRSegmentator `MRSeg/MRSeg.pt`, L2, Jacobian/GPU) | `EulerTransform` rigid, ASGD, 4 scales | ✅ ran (Docker + GPU; sample gate 0.1 — see IMPACT section) | | 4 | `niftyreg` | CT / T1c | image centres | NCC block matching + trimmed least squares (symmetric) | `reg_aladin -rigOnly`, 3-level pyramid | ✅ ran (Docker `niftyreg:2.1.1-ubuntu`, GPU `-platf 1`) | | 5 | `simpleitk` | CT / T1c | identity | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran | | 6 | `synthmorph` | CT / T1c | none (global) | **SynthMorph** anatomy-aware contrast-invariant (deep, pretrained) | 6-DOF rigid matrix | ✅ ran (main env, Docker + GPU) | @@ -139,7 +142,7 @@ pairs under the no-mask rule (see the IMPACT section). |---|---|---|---|---|---| | **antspyx** | ANTsPy 0.6.3 | 9 | 36.6 | 0.93 | **0.653 / 0.644** | | **elastix** | ITK-Elastix 0.25.4 | 9 | 6.0 | 0.93 | **0.517 / 0.520** | -| **impact** | Docker Elastix 5.2 + IMPACT | 0 | - | - | - (aborts, no mask) | +| **impact** | Docker Elastix 5.2 + IMPACT (MRSegmentator) | 8 | 2534 | 0.93 | **0.591 / 0.575** | | **niftyreg** | NiftyReg 2.1.1 (reg_aladin, Docker + CUDA) | 9 | 20.8 | 0.93 | **0.292 / 0.338** | | **simpleitk** | SimpleITK 2.5.6 | 9 | 10.6 | 0.93 | **0.651 / 0.641** | | **synthmorph** | Docker SynthMorph (TensorFlow) | 9 | 24.5 | 0.93 | **0.277 / 0.328** | @@ -149,15 +152,15 @@ Per-case MI (nats, unclipped score; **bold** = best for the case): | case | antspyx | elastix | impact | niftyreg | simpleitk | synthmorph | |---|---|---|---|---|---|---| | m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | -| m6_22IRRKU6_20240115 | **0.649** | 0.472 | abort | 0.292 | 0.642 | 0.277 | -| m6_23XIMO4M_20230222 | **0.655** | 0.416 | abort | 0.632 | 0.652 | 0.616 | -| m6_24V62QQM_20230223 | **0.704** | 0.701 | abort | 0.260 | 0.663 | 0.256 | -| m6_24V62QQM_20231017 | 0.601 | 0.590 | abort | 0.493 | **0.602** | 0.484 | -| m6_24V62QQM_20240820 | **0.653** | 0.644 | abort | 0.118 | 0.652 | 0.111 | -| m6_26TYQI3C_20221027 | 0.659 | 0.653 | abort | 0.609 | **0.695** | 0.590 | -| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | abort | 0.069 | **0.577** | 0.061 | -| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | abort | 0.233 | 0.634 | 0.231 | -| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | abort | 0.340 | 0.651 | 0.323 | +| m6_22IRRKU6_20240115 | **0.649** | 0.472 | 0.635 | 0.292 | 0.642 | 0.277 | +| m6_23XIMO4M_20230222 | **0.655** | 0.416 | 0.478 | 0.632 | 0.652 | 0.616 | +| m6_24V62QQM_20230223 | **0.704** | 0.701 | 0.697 | 0.260 | 0.663 | 0.256 | +| m6_24V62QQM_20231017 | 0.601 | 0.590 | 0.556 | 0.493 | **0.602** | 0.484 | +| m6_24V62QQM_20240820 | **0.653** | 0.644 | 0.531 | 0.118 | 0.652 | 0.111 | +| m6_26TYQI3C_20221027 | 0.659 | 0.653 | gate-abort | 0.609 | **0.695** | 0.590 | +| m6_2BBN5XA4_20250411 | **0.577** | 0.287 | 0.450 | 0.069 | **0.577** | 0.061 | +| m6_2C4TKI6H_20240506 | **0.645** | 0.401 | 0.629 | 0.233 | 0.634 | 0.231 | +| m6_2CDNN7BH_20240102 | **0.655** | 0.516 | 0.625 | 0.340 | 0.651 | 0.323 | ### Interpretation @@ -179,13 +182,19 @@ Per-case MI (nats, unclipped score; **bold** = best for the case): single CenterOfGravity start still locks onto local minima on the repositioned ones (`m6_2BBN5XA4` 0.287, `m6_2C4TKI6H` 0.401, `m6_23XIMO4M` 0.416). -- **IMPACT (semantic metric) does not run on these pairs under the no-mask - rule** — earlier (masked) runs showed it strongest on the hardest cases - (`m6_22IRRKU6` 0.624, `m6_2C4TKI6H` 0.626, `m6_2CDNN7BH` 0.618, median - 0.511) with a stochastic draw per run; it also was by far the slowest - (~14 min/case at the 500/2000 budget, hours/case at 2000/5000). With the CT - head-neck FOV as the fixed volume it aborts ("too many samples map outside - moving image buffer", ~10–20% in-buffer) — see the IMPACT section. +- **IMPACT (semantic metric) now runs on these pairs** — the no-mask abort was + elastix's sample-validity gate (build default `RequiredRatioOfValidSamples` + 0.25) rejecting the ~14% in-buffer fraction of the CT head-neck FOV vs the + tight T1c crop; the parameter map sets it to 0.1 (see the IMPACT section). + With MRSegmentator at the 2000/5000 budget it is top-three in MI on 7 of the + 8 registered cases (median 0.591): 0.635 on the hardest `m6_22IRRKU6` + (masked-era 0.624), 0.697 on `m6_24V62QQM_20230223` (near-best), 0.629 on + `m6_2C4TKI6H` — and it shares the production-acceptance of the grossly + repositioned `m6_2BBN5XA4` (std 10.7 mm) with antspyx/simpleitk, the only + non-intensity tool to do so. One run diverges mid-registration + (`m6_26TYQI3C`: the in-buffer fraction falls below 10% at iteration 256 and + the gate aborts — a stochastic-draw outcome). It remains by far the slowest + (~42 min/case at 2000/5000). - **NiftyReg (reg_aladin) is mid-pack**: single-start symmetric block matching is strong on moderate cases (`m6_23XIMO4M` 0.632, best in class; `m6_26TYQI3C` 0.609) but falls into the same tilt local minima @@ -205,9 +214,11 @@ Per-case MI (nats, unclipped score; **bold** = best for the case): CT head everywhere, so coverage is not the discriminator — MI is. - **Speed** (median wall, 4 cases in parallel, 24 threads): elastix 6.0 s, SimpleITK 10.6 s, **NiftyReg 20.8 s** (v2.1.1, GPU block matching), - SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT does not run on - M6 under the no-mask rule; when it did — masked, 500/2000 budget — it was - ~833 s/case, GPU feature extraction with backprop per iteration.) + SynthMorph 24.5 s (GPU forward pass), ANTs 36.6 s. (IMPACT is by far the + slowest — 2534 s/case median at the 2000/5000 benchmark budget, ~276 s/case + at the authors' 500/2000; GPU feature extraction with backprop per + iteration. No-mask runs are faster per iteration than the masked-era + ~833 s/case @500/2000 because only the ~14% in-buffer samples are scored.) - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its per-case numbers carry run-to-run variance (this table is one draw); SimpleITK's MI (and registration) is deterministic at fixed thread count; @@ -235,7 +246,7 @@ Median achieved MI (nats) and wall-clock (9 complete-CT cases): | niftyreg | 0.296 / 23.8 s | 0.320 / 20.2 s | +0.024 | 5/9 | | simpleitk | 0.651 / 10.6 s | 0.436 / 11.1 s | **−0.215** | 0/9 | | synthmorph | 0.241 / 26.0 s | 0.318 / 25.6 s | +0.077 | 5/9 | -| impact | abort | abort | - | - (still aborts; clipping does not fix the FOV geometry, 194–388/2000 samples in-buffer) | +| impact | abort | abort | - | - (historical: at the time both aborted on the old 0.25 sample-validity gate; the gate is now 0.1, A/B to be re-done) | **Verdict: no — clipping the CT to [0, 80] HU does not improve the result; do not clip.** @@ -286,8 +297,10 @@ python scripts/bench_registration.py --pipelines simpleitk elastix antspyx nifty # sections below for the one-time setup) python scripts/bench_registration.py --pipelines synthmorph --n-cases 10 \ --workers 4 --threads 24 -# NOTE: impact aborts on the M6 pairs under the shared no-mask rule (the CT -# head-neck FOV maps mostly outside the T1c head crop) — see the IMPACT section +# NOTE: impact now runs on the M6 pairs (the parameter map relaxes the +# sample-validity gate to RequiredRatioOfValidSamples 0.1) but is by far the +# slowest: ~5 min/case at 500/2000, ~42 min/case (2534 s median) at 2000/5000 +# (--impact-iters / --impact-samples) — see the IMPACT section python scripts/bench_registration.py --pipelines impact --n-cases 10 \ --workers 4 --threads 24 @@ -346,8 +359,7 @@ python scripts/bench_registration.py \ --pipelines synthmorph --n-cases 10 --qa-images ``` -(9 images each; the incomplete-CT case is skipped. `impact` produces none on -M6 — it aborts under the no-mask rule, see the IMPACT section.) +(9 images each; the incomplete-CT case is skipped.) ## SynthMorph pipeline (Docker) — one-time setup @@ -442,14 +454,17 @@ conda install -c conda-forge niftyreg The `impact` pipeline runs the [`vboussot/elastix_impact`](https://hub.docker.com/r/vboussot/elastix_impact) image: **Elastix 5.2.0** with the [IMPACT](https://arxiv.org/abs/2503.24121) *semantic* similarity metric compiled in as a native Elastix metric, plus a -TorchScript **TotalSegmentator** feature extractor. It is GPU-based and +TorchScript segmentation feature extractor (**MRSegmentator** `MRSeg/MRSeg.pt` +by default; TotalSegmentator `TS/M730.pt` and other models from +`VBoussot/impact-torchscript-models` are selectable). It is GPU-based and Docker-based (alongside `synthmorph`). It follows the shared rules (full -volumes, no mask, fixed = CT, moving = T1c) — which means it **aborts on the -M6 pairs**: the IMPACT metric only evaluates samples that map inside the -moving buffer, and the CT head-neck FOV (~300×300×251 mm) is much larger than -the T1c head crop (~220×220×150 mm), so ~14% of the samples are in-buffer at -the auto-init and elastix raises "too many samples map outside moving image -buffer" (see the IMPACT specifics below). +volumes, no mask, fixed = CT, moving = T1c). The metric only evaluates samples +that map inside the moving buffer, and the CT head-neck FOV (~300×300×251 mm) +is much larger than the T1c head crop (~220×220×150 mm), so ~14% of the samples +are in-buffer at the auto-init. The image build's sample-validity gate +defaults to `RequiredRatioOfValidSamples 0.25` and used to abort the M6 pairs +on that ("too many samples map outside moving image buffer"); the parameter +map now sets it to `0.1` (see the IMPACT specifics below). Prerequisites (check with `RP.availability()`; a missing piece raises `PipelineUnavailable` with the exact command): @@ -462,16 +477,23 @@ docker --version docker pull vboussot/elastix_impact # 3. Fetch the TorchScript model (~123 MB, one-time) into the repo's model -# cache (data/impact_models/). M730.pt is the recommended TotalSegmentator. -mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/TS -curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/TS/M730.pt \ - https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/TS/M730.pt +# cache (data/impact_models/). MRSeg.pt (MRSegmentator) is the default: +# it is trained on both MRI and CT, the right prior for this +# cross-modality CT->T1c task (all models share the same 2^l+3 FOV +# architecture, so the parameter map is unchanged between them). +mkdir -p /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg +curl -L -o /mnt/t24/Public/xfr/longitudinal/impact_models/MRSeg/MRSeg.pt \ + https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/MRSeg/MRSeg.pt ``` -Other models (M258/M291/M850, MIND, SAM2.1, …) are selectable with -`--impact-model `; see the +Other models (TotalSegmentator `TS/M730.pt`, `TS/M258/M291/M850`, MIND, SAM2.1, +Anatomix, …) are selectable with `--impact-model `; +see the [impact-torchscript-models](https://huggingface.co/VBoussot/impact-torchscript-models) -repo for the full list and FOV/layer tables. +repo for the full list and FOV/layer tables. The model choice does not affect +the sample-validity gate (below): the "too many samples outside moving buffer" +check is geometric (CT head-neck FOV vs the tight T1c crop) and fires at the +first sample evaluation, before feature extraction. ### How the adapter drives the image @@ -479,14 +501,16 @@ The image's `CMD.sh` is a generic launcher, so the adapter calls the in-image binary directly (`/lib/elastix-install/bin/elastix`): 1. **Full volumes, fixed=CT, moving=T1c** — both volumes are used as-is (no - head box, no mask, no crop / pad, per the shared rules above). **On the M6 - pairs this aborts**: the IMPACT metric only evaluates spatial samples that - map inside the moving buffer, and the CT FOV (~300×300×251 mm, head+neck) - is much larger than the tight T1c head crop (~220×220×150 mm), so at the - CenterOfGravity auto-init only ~14% of the 2000 samples are in-buffer and - elastix raises "Too many samples map outside moving image buffer: 277 / - 2000". IMPACT therefore needs a fixed volume whose FOV is comparable to - the moving crop (e.g. a head-only CT) to work under the no-mask rule. + head box, no mask, no crop / pad, per the shared rules above). The metric + only evaluates spatial samples that map inside the moving buffer, and the + CT FOV (~300×300×251 mm, head+neck) is much larger than the tight T1c head + crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of + the samples are in-buffer. The image build's default validity gate + (`RequiredRatioOfValidSamples 0.25`) aborted on that ("Too many samples map + outside moving image buffer: 704 / 5000"); the parameter map therefore sets + `(RequiredRatioOfValidSamples 0.1)` so the M6 pairs run on their in-buffer + samples while a genuinely diverging transform (<10% in-buffer) still + aborts. 2. **Auto init (no `-t0`)** — the parameter map sets `AutomaticTransformInitialization "true"` with method `CenterOfGravity`, so elastix aligns the intensity centroids; `AutomaticScalesEstimation` sets the @@ -495,9 +519,24 @@ binary directly (`/lib/elastix-install/bin/elastix`): **Jacobian** mode (differentiable, per the recommended config), 4 semantic scales `ImpactVoxelSize` 6/3/1.5/1 mm (scaled to the T1c base spacing) with **no image resampling** (all-1 pyramid schedules, as the IMPACT docs require), - ASGD, `Random` sampler. GPU auto-spread across the host GPUs by pid - (`--impact-gpu -1` forces CPU). -4. **Read-back** — the result `TransformParameters.0.itk.txt` is an + ASGD, `Random` sampler, `(RequiredRatioOfValidSamples 0.1)`. A copy of the + exact parameter map is dumped to + `data/m6reg/impact/paramdump/ParameterMap__.txt` before every + docker run (also recorded as `param_dump` in the result JSON). GPU + auto-spread across the host GPUs by pid (`--impact-gpu -1` forces CPU). +4. **Intensity standardization happens inside the metric** (no pipeline-side + normalization — the model card's "standardize intensities to zero mean and + unit variance" is already applied): at init the metric computes + whole-volume `[min, max, mean, sigma]` for each image + (`itk::StatisticsImageFilter`) and passes it to the model, whose + `normalize` branch rescales every extracted patch with the image's own + `(mean, sigma)`. It is computed over the whole buffer (CT air included) and + identically for fixed and moving, so both modalities' features sit on the + same normalized scale. Pre-standardizing the MHA files in the pipeline + would be an exact no-op (z-scoring is affine-invariant: the metric + recomputes the stats of the image it receives) and would only corrupt the +HU scale used by MI scoring and the QA brain-window images. + 5. **Read-back** — the result `TransformParameters.0.itk.txt` is an `Euler3DTransform` in the world physical frame (the `AutomaticTransformInitialization` is folded in by elastix); the adapter runs it through the usual `fix_orientation` gate. @@ -590,29 +629,41 @@ binary directly (`/lib/elastix-install/bin/elastix`): defaults (`NumberOfSpatialSamples` 5000, `NumberOfHistogramBins` 32). - **IMPACT specifics** (Docker `vboussot/elastix_impact`, Elastix 5.2.0, GPU, Jacobian mode): -- **Slow, by design.** Every ASGD iteration runs the TotalSegmentator - extractor on the sampled patches with backpropagation. The benchmark - budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 = - the elastix default) is ~10× the authors' 500-iter / 2000-sample - Jacobian budget, i.e. on the order of *hours* per case on full volumes — - the results shown for `impact` are from the previous 500/2000 budget and - were not re-run at the new budget. `--impact-iters` / `--impact-samples` - trade accuracy for speed. -- **Aborts on M6 under the no-mask rule.** Like every other pipeline it - runs fixed=CT/moving=T1c on the full volumes (no mask, no crop / pad). - The CT FOV (~300×300×251 mm) is much larger than the T1c crop - (~220×220×150 mm), so most of the CT buffer maps outside the T1c moving - buffer and the metric aborts ("Too many samples map outside moving image - buffer: 277 / 2000" on M6 at the CenterOfGravity auto-init). The - transform contract (T1c→CT) still holds when a run completes; for M6 - pairs a head-only CT as the fixed volume would be needed. +- **Slow, by design.** Every ASGD iteration runs the segmentation + extractor (MRSegmentator by default) on the sampled patches with + backpropagation. The benchmark +budget (MaximumNumberOfIterations 2000, NumberOfSpatialSamples 5000 = + the elastix default) is ~10× the authors' 500-iter / 2000-sample + Jacobian budget: measured 2534 s/case (median) on the no-mask M6 + pairs, ~276 s/case at 500/2000. `--impact-iters` / `--impact-samples` + trade accuracy for speed. +- **Runs on M6 with a relaxed sample-validity gate.** Like every other + pipeline it runs fixed=CT/moving=T1c on the full volumes (no mask, no + crop / pad). The CT FOV (~300×300×251 mm) is much larger than the T1c + crop (~220×220×150 mm), so at the CenterOfGravity auto-init only ~14% of + the samples map inside the T1c moving buffer. The metric calls + `CheckNumberOfSamples()` after every evaluation and the image build + defaults to `RequiredRatioOfValidSamples 0.25`, so it used to abort + ("Too many samples map outside moving image buffer: 704 / 5000"). The + parameter map now sets `(RequiredRatioOfValidSamples 0.1)`: the gate + still aborts a genuinely diverging transform (<10% in-buffer) while the + FOV-mismatched M6 pairs run on their ~14% in-buffer samples — exactly + the native "only in-buffer samples are evaluated" behaviour. In the + 2000/5000 re-run the gate still fired once, mid-registration: on + `m6_26TYQI3C` the transform drifted until only 468/5000 samples (9.4%) + stayed in-buffer at iteration 256 and the run aborted. The transform + contract (T1c→CT) holds; `fix_orientation` decides the direction as + usual (all 8 completed runs report `auto_inverted: True`). - **Auto init.** `AutomaticTransformInitialization` (CenterOfGravity) aligns the intensity centroids; no manual `-t0` translation is used. - **History: single-init convergence.** The semantic landscape has more local minima than Mattes-MI; earlier (masked, fixed=CT) runs no longer nailed the hardest repositioned case (`m6_22IRRKU6` mid 22.7 / std 27.0; - the pre-mask head-box flow got it to mid 8.3 / std 13.2). Under the - current no-mask rule it does not run on M6 at all (see above). + the pre-mask head-box flow got it to mid 8.3 / std 13.2). The no-mask + MRSegmentator runs show the same pattern: strong MI (0.635 on + `m6_22IRRKU6` at 2000/5000) but the pose still misses the concentricity + acceptance rule there (mid 26.1 / std 28.6); it *is* accepted on + `m6_2BBN5XA4` (mid 15.1 / std 10.7). - **Stochastic.** Jacobian mode resamples patches each iteration, so results vary run to run (see Stochasticity). - **Elastix 5.x differences vs the in-process `elastix` pipeline:** the image diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index dcbb2e2..a3a8d1e 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -654,7 +654,7 @@ def register_synthmorph(ct, t1c, model="rigid", gpu=None, threads=None, IMPACT_IMAGE = "vboussot/elastix_impact" IMPACT_ELASTIX_BIN = "/lib/elastix-install/bin/elastix" IMPACT_MODEL_CACHE = "data/impact_models" -IMPACT_DEFAULT_MODEL = "TS/M730.pt" +IMPACT_DEFAULT_MODEL = "MRSeg/MRSeg.pt" IMPACT_HF_MODEL_URL = ( "https://huggingface.co/VBoussot/impact-torchscript-models/resolve/main/{rel}" ) @@ -757,6 +757,14 @@ def _impact_param_map(model_cont_path, gpu, iters, samples, voxel_sizes, f'(ImpactUseMixedPrecision "{mixed}")', "(ImpactFeaturesMapUpdateInterval -1)", '(ImpactWriteFeatureMaps "false")', + # safety gate: abort if < ratio of the spatial samples map inside the + # moving buffer. The build default is 0.25, but the no-mask M6 pairs + # (CT head-neck FOV vs the tight T1c crop) only have ~14% in-buffer at + # the CenterOfGravity auto-init, so elastix aborts at the first metric + # evaluation. 0.1 keeps the gate (a diverging transform still aborts) + # while letting the FOV-mismatched runs proceed -- the metric natively + # evaluates only the in-buffer samples. + "(RequiredRatioOfValidSamples 0.1)", "", "// rigid 6-DOF to be optimized (auto center init, no -t0)", '(Transform "EulerTransform")', @@ -820,12 +828,21 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000, ``vboussot/elastix_impact`` (Elastix 5.2.0 + IMPACT semantic metric). IMPACT (arXiv:2503.24121) replaces the intensity metric with a *semantic* - one: a TorchScript TotalSegmentator model (``TS/M730.pt`` by default) is run - on randomly sampled patches of the fixed and moving images and the + one: a TorchScript segmentation model (``MRSeg/MRSeg.pt`` by default, + MRSegmentator -- robust multi-modality segmentation on MRI *and* CT, the + right prior for this cross-modality CT->T1c task; TotalSegmentator + ``TS/M730.pt`` and the other models in the HF repo + ``VBoussot/impact-torchscript-models`` work too via ``model=``) is run on + randomly sampled patches of the fixed and moving images and the registration minimizes the distance between the extracted feature maps - (L2, per the recommended config). This is a genuinely different metric from - the Mattes-MI pipelines above -- it is the cross-modality tool the image was - built for. + (L2, per the recommended config). Intensity standardization (zero mean / + unit variance per the model card) is applied *inside the metric*: at init + it computes whole-volume [min, max, mean, sigma] per image and the model's + normalize branch rescales every patch with that image's own (mean, sigma), + so no pipeline-side pre-normalization is needed (and would be an exact + no-op -- z-scoring is affine-invariant). This is a genuinely different + metric from the Mattes-MI pipelines above -- it is the cross-modality tool + the image was built for. The pipeline registers on the full CT / T1c volumes (no head box) with no manual centroid init: elastix's ``AutomaticTransformInitialization`` @@ -910,10 +927,22 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000, os.path.join(datadir, "Fixed_image.mha"), True) sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), os.path.join(datadir, "Moving_image.mha"), True) - with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: - f.write(_impact_param_map(os.path.join("/Data/Models", model_rel), + param_text = _impact_param_map(os.path.join("/Data/Models", model_rel), gpu, iters, samples, voxel_sizes, - rescale=schedule, smoothing=smoothe)) + rescale=schedule, smoothing=smoothe) + with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: + f.write(param_text) + + # debug aid: keep a persistent copy of the exact parameter map elastix + # is about to run, written BEFORE the docker call so an abort leaves it + # behind for inspection (per-pid name: parallel workers do not clobber + # each other's dumps) + dump_dir = d("data/m6reg/impact/paramdump") + os.makedirs(dump_dir, exist_ok=True) + family = os.path.basename(os.path.dirname(model_rel)) or "model" + param_dump = os.path.join(dump_dir, f"ParameterMap_{family}_{os.getpid()}.txt") + with open(param_dump, "w") as f: + f.write(param_text) cmd = ["docker", "run", "--rm"] if gpu >= 0: @@ -949,17 +978,19 @@ def register_impact(ct, t1c, model=None, gpu=None, iters=2000, samples=5000, if not keep_scratch: shutil.rmtree(tmp, ignore_errors=True) + family = os.path.basename(os.path.dirname(model_rel)) \ + or os.path.splitext(os.path.basename(model_rel))[0] meta = { "pipeline": "impact", "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", - "metric": f"IMPACT semantic (TotalSegmentator {os.path.basename(model_rel)}, L2, Jacobian)", + "metric": f"IMPACT semantic ({family} {os.path.basename(model_rel)}, L2, Jacobian)", "transform": (f"Rigid (EulerTransform, ASGD, {n} scales, auto center init; " f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), "convention": "T1c->CT (verified)", "auto_inverted": inverted, "fixed": _modality_name(fixed), "moving": _modality_name(moving), "seconds_reg": round(dt, 1), "docker_image": IMPACT_IMAGE, "model": model_rel, "gpu": gpu, - "iters": iters, "samples": samples, + "iters": iters, "samples": samples, "param_dump": param_dump, } if keep_scratch: meta["scratch"] = tmp