diff --git a/AGENTS.md b/AGENTS.md index 117a6e1..071b94e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -10,7 +10,7 @@ Activate with: source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 ``` -Key packages: torch 2.14 (+cu126), torchvision, monai 1.6 (pip, Pipeline C), numpy, scipy, pandas, scikit-learn, scikit-image, matplotlib, nibabel, SimpleITK. Registration toolkits (benchmark, see REGISTRATION_BENCHMARK.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), 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). +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_RIGID.md): `itk` + `itk-elastix`, `antspyx`, niftyreg 1.5.69 (`conda install -c conda-forge niftyreg`, provides `reg_aladin`), plus the Docker images `freesurfer/synthmorph` and `vboussot/elastix_impact` (IMPACT; rigid MIND + TotalSegmentator config, feature models under `data/impact_models`). The benchmark's ImpactLoss scoring metric uses `src/impact_loss.py` (vendored from vboussot/ImpactLoss; runs on the env's torch). ## Project diff --git a/REGISTRATION_BENCHMARK_DEFORMABLE.md b/REGISTRATION_BENCHMARK_DEFORMABLE.md new file mode 100644 index 0000000..9743448 --- /dev/null +++ b/REGISTRATION_BENCHMARK_DEFORMABLE.md @@ -0,0 +1,310 @@ +# Registration Pipeline Benchmark — CT → T1c (M6) — Deformable + +Companion to [REGISTRATION_BENCHMARK_RIGID.md](REGISTRATION_BENCHMARK_RIGID.md): +the **deformable** variant of the same six-pipeline benchmark. The rigid +benchmark asked "which toolkit aligns the head best with 6 DOF?"; this one +asks "which toolkit's native deformable workflow aligns the head (and the +tissue) best?" on the same 10 M6 CT→T1c test cases, the same unified +contract, and the same metrics. + +Run date: 2026-10-02 (main run 05:43–07:52 local; serial `synthmorph` +re-run 09:58–10:07 local; accelerated `simpleitk` re-run 15:31–16:26 +local). 10 M6 test cases, 3× NVIDIA A100 80GB PCIe, 96-core host, conda +env `longitudinal-313`. Main run: 6 workers (24 threads each, 2 workers +per GPU). `synthmorph` was re-run serially (1 worker, one container per +GPU at a time) because the concurrent run lost 6/9 cases to cuDNN +scratch-space allocation failures (see the SynthMorph note under +Interpretation); the re-run registered 9/9, matching the concurrent +run's scores to 4 decimals on the 3 cases it completed. `simpleitk` +was re-run after its spline grid was accelerated from 16 mm to 24 mm +(see the simpleitk row of the pipeline table; the original 16 mm results +are archived as `simpleitk_16mm_legacy.json`). + +## Unified contract + +Identical to the rigid benchmark: the returned transform `T` must satisfy +`sitk.Resample(ct, t1c, T, ...)` = CT on the T1c grid, i.e. **T maps T1c → +CT coordinates**. fixed = T1c, moving = CT, full volumes (no head box, no +mask, no crop, no pad), each toolkit's native automatic transform +initialization, CT clipped to [0, 80] HU only for the pipelines where the +rigid A/B showed clipping helps (niftyreg, synthmorph). + +One deformable-specific change: every pipeline's output (elastix BSpline +combination, ANTs warp + affine composite, NiftyReg affine + F3D spline, +SimpleITK Euler + BSpline, SynthMorph joint field, IMPACT RecursiveBSpline) +is **rasterised to a SimpleITK `DisplacementFieldTransform` on the T1c +grid** before scoring, so all six pipelines return the same transform type +and the shared scoring / GTV warping / QA code is used verbatim. + +## The six pipelines (deformable configurations) + +Each pipeline runs its toolkit's **native deformable workflow** — a linear +init stage + a deformable stage, or a single deformable stage where that is +the toolkit's canonical setup: + +| pipeline | deformable configuration | metric | +|---|---|---| +| **antspyx** | `type_of_transform='SyN'` (ANTs' native deformable): Affine stage (2100×1200×1200×0, shrink 4×2×2×1) + SyN[0.2,3,0] deformable stage (40×20×0, shrink 4×2×1), center-of-mass auto init | Mattes MI (32 bins), both stages | +| **elastix** | two-stage (the toolkit's canonical deformable workflow): EulerTransform rigid stage (2000 iters × 3 resolutions, CenterOfGravity auto init) → BSplineTransform stage (16 mm grid, spline order 3, 2000 iters × 3 resolutions) initialized from stage 1 via `-t0` | Mattes MI (32 bins), 5000 samples/level | +| **impact** | the ImpactReg `MR_CT_TS` **reference map in its native deformable form**: RecursiveBSplineTransform (14 mm final grid, order 3) driven by the MIND + TotalSegmentator static feature combination (4 per-level scales 400/300/200/200 iters, 2000 samples) with the reference's three-metric combination: IMPACT (w 1.0) + Mattes MI (w 0.4) + TransformBendingEnergyPenalty (w 10). Two M6 adaptations, same as the rigid benchmark: CenterOfGravity auto init (the reference's GeometricalCenter start leaves the heads ~100 mm apart and aborts) and RequiredRatioOfValidSamples 0.1 | IMPACT semantic (L1 + Dice) + Mattes MI + bending energy | +| **niftyreg** | reg_aladin rigid init (3-level pyramid, 15 TLS iters/level) → reg_f3d non-rigid (F3D B-spline Free-Form deformation, Modat et al. 2010; NMI metric, symmetric bending/linear regularisation, 3-level pyramid, 150 iters at the final level, 8 mm final grid) | NCC block matching (init) + NMI (F3D) | +| **simpleitk** | two-stage (the toolkit's canonical deformable demo): Euler3D stage (2000 iters/level, pyramid [4,2,1]) → BSplineTransform3D stage (24 mm grid, 50 iters/level on pyramid [4,2]), the BSpline's control points initialised to reproduce the Euler stage exactly. The 24 mm grid (5577 parameters) is the accelerated setting (the 2026-10-02 original used 16 mm / 13005 parameters, ~61 min/case under the harness): this build's spline per-iteration cost scales ~linearly with the control-point count (~2.8 s/iter at the shrink-4 level and ~20 s/iter at shrink-2 for 24 mm vs ~6.4 / ~45.8 s for 16 mm — independent of the metric sample count, which is fixed in this SimpleITK build; LBFGS2 was tested and fails to move the transform on this noisy MI), and the full re-run keeps the 16 mm quality within 0.003–0.029 MI per case (median 0.698 vs 0.705, ImpactLoss median unchanged) while retaining simpleitk's `m6_26TYQI3C` MI/ImpactLoss win — see the acceleration note under Interpretation. A full-resolution polish pass (10 iters at shrink-1) was tested and dropped the achieved MI (0.6891 at 32 mm), so the pyramid stays [4,2] | Mattes MI (32 bins), both stages | +| **synthmorph** | `joint` model (affine + deformable in one symmetric network pass — the toolkit's canonical deformable; the `deformable` model requires a pre-existing field and would bias the comparison) | anatomy-aware contrast-invariant network (pretrained) | + +## Test cases + +First 10 M6 labelled cases of `m6.discover()` (all with a GTV label): +`m6_22GD5VL2_20230724`, `m6_22IRRKU6_20240115`, `m6_23XIMO4M_20230222`, +`m6_24V62QQM_20230223`, `m6_24V62QQM_20231017`, `m6_24V62QQM_20240820`, +`m6_26TYQI3C_20221027`, `m6_2BBN5XA4_20250411`, `m6_2C4TKI6H_20240506`, +`m6_2CDNN7BH_20240102`. + +(`m6_22GD5VL2_20230724` has an incomplete CT head — a known data issue, +not a pipeline failure; it is scored as ct_incomplete in every pipeline, +as in the rigid benchmark.) + +(Ground-truth reference transforms are **not** part of this benchmark — +scoring is by achieved image similarity only, as in the rigid benchmark.) + +## Metrics (identical for every pipeline) + +- **Mattes MI** (nats) between the warped CT and the T1c over the MR head + (largest connected component of the T1c mask) — achieved alignment, higher + = better. +- **IMPACT loss** (vboussot/ImpactLoss; pretrained TotalSegmentator + `TS/M730` deep features — weighted L1 over the two deepest layers, full + volumes; the same scorer as the rigid benchmark) between the warped CT + and the T1c — semantic-anatomy similarity, lower = better. +- **Registration wall time** (seconds, excludes I/O + scoring). + +Per-case JSONs additionally retain the production skull-ring QA (gap + +concentricity, and the acceptance verdict) for reference. + +### Result + +Alignment is scored by the achieved **Mattes MI** (nats, warped-CT ~ T1c +over the MR head; higher = better) **and** the **ImpactLoss** semantic score +(warped-CT ~ T1c over the full volumes; lower = better), on the same +unclipped scoring volumes as the rigid benchmark, so all rows are comparable +on the same scales and against the rigid numbers. Full volumes, no mask / +crop / pad, fixed = T1c / moving = CT (shared rule). `niftyreg` and +`synthmorph` clip the CT (the moving image) to [0, 80] HU *before +registration* (per the rigid benchmark's clipping A/B); `antspyx`, +`elastix`, `impact` and `simpleitk` register the CT as-is. All six +transforms are deformable (each toolkit's native workflow) and are +rasterised to T1c-grid displacement fields before scoring. + +All six pipelines registered 9/9 cases (the 10th is the known +incomplete-CT data issue, unscored in every pipeline). **Bold** = best +for the column (MI higher is better, ImpactLoss lower is better). + +| pipeline | toolkit | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) | +|---|---|---|---|---|---| +| **antspyx** | ANTsPy 0.6.3 | 9 | 31.7 | 0.158 / 0.156 | **0.723 / 0.729** | +| **elastix** | ITK-Elastix 0.25.4 | 9 | 45.8 | 0.164 / 0.164 | 0.622 / 0.627 | +| **impact** | Docker `vboussot/elastix_impact` (Elastix 5.2.0 + IMPACT, deformable MIND+TS reference map) | 9 | 273.2 | 0.159 / 0.159 | 0.553 / 0.547 | +| **niftyreg** | NiftyReg 2.1.1 (reg_aladin + reg_f3d, Docker + CUDA) | 9 | 38.1 | 0.161 / 0.161 | 0.666 / 0.686 | +| **simpleitk** | SimpleITK 2.5.6 | 9 | 1369.9 | 0.158 / 0.159 | 0.698 / 0.704 | +| **synthmorph** | Docker `freesurfer/synthmorph` (TensorFlow, pretrained) | 9 | 49.1 | **0.156 / 0.156** | 0.627 / 0.640 | + +Per-case MI (nats, unclipped score; **bold** = best for the case). The +`synthmorph` column is from the serial re-run (see the SynthMorph note +under Interpretation): + +| 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.679** | 0.580 | 0.548 | 0.637 | 0.654 | 0.620 | +| m6_23XIMO4M_20230222 | **0.723** | 0.622 | 0.573 | 0.711 | 0.693 | 0.618 | +| m6_24V62QQM_20230223 | **0.789** | 0.570 | 0.635 | 0.729 | 0.738 | 0.701 | +| m6_24V62QQM_20231017 | **0.720** | 0.622 | 0.580 | 0.665 | 0.692 | 0.628 | +| m6_24V62QQM_20240820 | **0.738** | 0.674 | 0.462 | 0.707 | 0.700 | 0.666 | +| m6_26TYQI3C_20221027 | 0.799 | 0.694 | 0.537 | 0.793 | **0.801** | 0.709 | +| m6_2BBN5XA4_20250411 | **0.666** | 0.585 | 0.436 | 0.612 | 0.641 | 0.569 | +| m6_2C4TKI6H_20240506 | **0.719** | 0.635 | 0.553 | 0.656 | 0.698 | 0.627 | +| m6_2CDNN7BH_20240102 | **0.726** | 0.660 | 0.602 | 0.666 | 0.718 | 0.624 | + +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.1541 | 0.1595 | 0.1532 | 0.1530 | 0.1572 | **0.1512** | +| m6_23XIMO4M_20230222 | 0.1581 | 0.1657 | **0.1567** | 0.1605 | 0.1603 | 0.1573 | +| m6_24V62QQM_20230223 | 0.1473 | 0.1634 | 0.1497 | 0.1541 | 0.1502 | **0.1470** | +| m6_24V62QQM_20231017 | **0.1518** | 0.1591 | 0.1547 | 0.1620 | 0.1548 | **0.1518** | +| m6_24V62QQM_20240820 | **0.1545** | 0.1646 | 0.1694 | 0.1706 | 0.1584 | 0.1562 | +| m6_26TYQI3C_20221027 | 0.1580 | 0.1638 | 0.1656 | 0.1672 | **0.1562** | 0.1602 | +| m6_2BBN5XA4_20250411 | 0.1624 | 0.1663 | 0.1669 | 0.1648 | 0.1631 | **0.1616** | +| m6_2C4TKI6H_20240506 | 0.1594 | 0.1640 | 0.1585 | **0.1552** | 0.1626 | 0.1564 | +| m6_2CDNN7BH_20240102 | 0.1590 | 0.1662 | 0.1597 | 0.1605 | 0.1633 | **0.1583** | + +### Interpretation + +- **`antspyx` (Affine + SyN) leads on achieved MI** (median 0.723; best on + 7 of 9 cases and within 0.02 of best on the other 2) at the best + median ImpactLoss of the gradient-based tools (0.158, within 0.002 of + `synthmorph`'s 0.156), and it is the fastest full-workflow pipeline + (31.7 s/case) — the best time/quality trade of the six. +- **`simpleitk` is second on MI** (0.698; it keeps its per-case win on + `m6_26TYQI3C` — 0.801, still the highest MI of the whole benchmark, + though now only 0.001 ahead of `antspyx` — and cedes its + `m6_2CDNN7BH` win to `antspyx`, 0.718 vs 0.726) at a cost of ~23 + min/case under the harness — ~43× the fastest pipeline after the + spline-grid acceleration below (this build's BSpline per-iteration + cost scales ~linearly with the control-point count, independent of the + metric sample count — see the pipeline table). +- **`simpleitk` acceleration (2026-10-02):** the original setting (16 mm + spline grid, 13005 parameters) ran 3658.9 s/case (median) under the + harness. The spline stage's per-iteration cost in this build scales + ~linearly with the control-point count (16 mm: ~6.4 s/iter at the + shrink-4 level, ~45.8 s at shrink-2; 24 mm: ~2.8 / ~20 s; 32 mm: + ~1.2 / ~9.6 s), and the metric sample count is fixed in this SimpleITK + build (`SetMetricAsMattesMutualInformation` exposes only the histogram + bins); LBFGS2 was tested and its line search rejects every step on + this noisy MI (the transform does not move). Switching to the **24 mm + grid (5577 parameters) cut the wall to 1369.9 s median (2.7×)** at a + cost of −0.003…−0.029 MI on 8 of 9 cases (+0.017 on `m6_2C4TKI6H`; + median 0.698 vs 0.705), with the median ImpactLoss unchanged (0.158) + and `simpleitk`'s `m6_26TYQI3C` MI/ImpactLoss wins intact. A 32 mm + grid (a further ~1.9×) was rejected: it drops below `antspyx`'s 0.7994 + on `m6_26TYQI3C` (0.7902, and 0.7963 even with 100 fine-level iters). + A full-resolution polish pass (10 iters at shrink-1) was also tested + and dropped the achieved MI (0.6891 at 32 mm). The original 16 mm + results are archived as `simpleitk_16mm_legacy.json`. +- **`niftyreg`'s `reg_f3d` stage delivers the largest deformable gain of + any pipeline** (0.292 → 0.666 median MI, +0.374): it recovers the two + gross-offset cases where the rigid block-matching locked onto local + minima (`m6_24V62QQM_20240820` 0.118 → 0.707, `m6_2BBN5XA4` 0.069 → + 0.612), lifts the toolkit to third on MI, and it is the most + concentric of the six (median skull-ring gap 5.5 mm mid / 3.8 mm std, + 6/9 accepted). +- **The deformable stage does not uniformly pay off.** Against the rigid + benchmark's per-case MI it improves `antspyx` (0.652 → 0.723 median), + `simpleitk` (0.647 → 0.698, accelerated 24 mm setting) and `niftyreg` + (0.292 → 0.666), but + **degrades** `elastix` (0.646 → 0.622; mixed per case, worst regression + `m6_24V62QQM_20230223` 0.698 → 0.570) and `impact` (0.652 → 0.553 — + lower than its own rigid init on **all 9 cases**). The two degradations + are the coarse-grid BSpline configurations (14–16 mm) with heavy + bending regularisation — for `impact` the reference map's + TransformBendingEnergyPenalty at weight 10 — so the deformable stage + smooths rather than refines, and (for `impact`) its IMPACT-dominant + objective is simply not the achieved-MI objective. +- **`impact`'s deformable reference map is the weakest on the intensity + scale** (MI 0.553) even though its *semantic* score is competitive + (ImpactLoss 0.159, within 0.003 of the best median) — a + metric-objective/achieved-MI divergence the rigid benchmark did not + show (there the two scores agreed, with `impact` in the top tier on + both). At 273 s/case it is the worst time/quality trade of the six. +- **`synthmorph` (serial re-run, 9/9) splits the leaderboards: best on + the semantic scale, mid-pack on the intensity scale.** Median + ImpactLoss 0.156 — the best of the six — with the per-case best on 4 of + 9 (plus a 4-decimal tie with `antspyx` on `m6_24V62QQM_20231017`), + while its MI (0.627) is fourth, statistically tied with `elastix` + (0.622). The pretrained contrast-invariant network optimises learned + anatomy rather than raw CT/MR intensity synchrony, so it wins where + anatomy matters and loses where the intensity histogram matters. Fast + (49.1 s/case) and the read-back cross-check against the tool's own + resampling is NCC 1.0 on all 9 cases. +- **The two score families discriminate differently here than in the + rigid benchmark.** The ImpactLoss scale is nearly saturated (medians + 0.156–0.164, ~5 % spread; per-case bests spread across five tools, + `elastix` with 0/9), whereas the MI scale spans 0.553–0.723 + (0.17 nats absolute) and the per-case bests go 8/9 `antspyx`, 1/9 + `simpleitk`, 0/9 the rest. The deformable stage's main effect is fine intensity alignment on + top of the gross anatomy the rigid stages already place — exactly what + MI measures and ImpactLoss is less sensitive to. For ranking the + deformable workflows, achieved MI is the more discriminative score. +- **`synthmorph` infrastructure note.** In the main concurrent run (6 + workers, 2 containers per GPU) 6 of 9 cases failed with `CUDNN failed + to allocate the scratch space ... [Op:Conv3D]` inside the joint model's + 256-channel 128³ hyperconvolution. Each case is a fresh + `docker run --rm` container, so the failures are cuDNN scratch-space + contention when two SynthMorph containers share one A100, not a + pipeline fault: the serial re-run registered 9/9, reproduced the + concurrent run's MI/ImpactLoss to 4 decimals on its 3 completed cases + (deterministic model), and ran ~17 s/case faster without the + contention (49.1 s vs 65.7 s median). The concurrent-run JSON is archived as + `synthmorph_concurrent_run.json`. +- **Production skull-ring QA (reference only, kept in the per-case + JSONs):** acceptance 6/9 `niftyreg`, 5/9 `elastix` / `synthmorph`, + 4/9 `antspyx` / `simpleitk`, 2/9 `impact`. Concentricity: `niftyreg` + best (median gap 5.5 mm mid / 3.8 mm std), `simpleitk` (19.1 mm std, + accelerated setting) and `impact` (18.9 mm) worst on the rotational + term. The same weak correlation with perceived alignment noted in the + rigid benchmark, so it is reported for reference only. + +## Transform conventions (deformable-specific) + +Unlike the rigid benchmark (where the elastix-style toolkits came out in the +opposite direction 8/9 cases and needed the best-of-T vs T⁻¹ orientation +gate), every deformable output here is consumed in the documented +fixed→moving = T1c→CT direction; the gate is replaced by an analytic +read-back plus a skull-ring sanity check. Two non-obvious conventions were +established empirically (validated against each tool's own resampling, NCC +≥ 0.998 on the head): + +- **NiftyReg and surfa (FreeSurfer's NIfTI stack) decode RAS+ NIfTI qforms + naively, i.e. they work in an x/y-mirrored world.** A displacement delta + (dx, dy, dz) stored by such a tool is (−dx, −dy, dz) in the true RAS + world; the field grid (file indices) is unaffected. The read-back flips + the x/y components of the dense field. (ANTs writes correct RAS fields — + no flip needed.) +- **NiftyReg reg_f3d writes only the control-point grid** (outputCPP.nii); + `reg_transform -disp` evaluates it to a dense field on the reference grid. + The reg_aladin affine is folded into the control points before + optimization (reg_f3d initialises the spline with the affine's + deformation field), so the dense field **is** the full transform — no + separate composition step. (The `-comp` subcommand was tested and does + not reproduce reg_f3d's internal transform; it is not used.) +- **SynthMorph joint**: driven as `register -m joint moving fixed` + (moving = clipped CT first, fixed = T1c second); the `-t` output is the + displacement field on the fixed (T1c) grid — the contract transform + directly. The tool's own `-o` resampling is written as a cross-check + (per-case `own_resample_ncc_head` in the JSONs). +- **ANTs SyN**: the two stage outputs (001Warp.nii.gz + 000GenericAffine.mat) + are composed as warp-then-affine (verified to reproduce ANTs' own + `warpedmovout` to < 0.001 MI) and rasterised. +- **Elastix / IMPACT**: the combination transform (auto-init folded in) is + rasterised — elastix in-process via `itk.TransformToDisplacementFieldFilter` + (exact, machine precision on affine fields), IMPACT by evaluating the + BSpline control-point file directly (SimpleITK cannot read elastix's + BSpline .txt format; scipy's order-3 B-spline = ITK's kernel). + +## Reproducing + +```bash +source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 +cd /mnt/b4/xfr/git26/longitudinal + +# all six deformable pipelines (results go to results/bench_registration_deformable/) +# the 2026-10-02 main run used --workers 6 --threads 24 (~2 h 9 min wall, +# dominated by simpleitk's ~111 min phase) +python scripts/bench_registration.py --deformable --workers 6 --threads 24 --qa-images + +# single pipeline / case subset +python scripts/bench_registration.py --deformable --pipelines simpleitk impact --n-cases 3 + +# serial SynthMorph re-run (the 2026-10-02 final numbers): one container per +# GPU at a time; two concurrent containers on one A100 hit cuDNN +# scratch-space OOM in the joint model's 256-ch 128^3 hyperconvolution +python scripts/bench_registration.py --deformable --pipelines synthmorph --workers 1 --threads 24 --qa-images + +# combine the per-pipeline JSONs into summary.json / summary.md / plot +python scripts/bench_registration.py --deformable --merge +``` + +## Visual QA images + +`--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to +`results/bench_registration_deformable/qa/__.png` +(9 per pipeline; the incomplete-CT case is skipped) — same format as the +rigid benchmark: a 3×2 grid of the mid-axial / mid-coronal / mid-sagittal +slices through the head center, **left column** checkerboard blend of the +warped CT vs the T1c (clean registration = seamless mosaic), **right +column** additive color overlay (CT red, T1c green; agreement = yellow, +misregistration = red/green fringing). CT is brain-windowed 0–80 HU, T1c +2–98 % head window; the title carries the pipeline, case key and the +achieved Mattes MI. \ No newline at end of file diff --git a/REGISTRATION_BENCHMARK_RIGID.md b/REGISTRATION_BENCHMARK_RIGID.md index 495aa99..23685f6 100644 --- a/REGISTRATION_BENCHMARK_RIGID.md +++ b/REGISTRATION_BENCHMARK_RIGID.md @@ -235,11 +235,14 @@ the case): pipelines the same way the intensity score does is a strong check that the leaderboard is not an artifact of the MI scale. - **The grossly repositioned cases still separate the leaders from the - tail.** On `m6_2BBN5XA4` / `m6_24V62QQM_20240820` (~30° inter-session tilt) - the top four hold 0.57–0.65, while niftyreg (0.069 / 0.118) and synthmorph - (0.061 / 0.111) — the single-start block-matching and pretrained-MRI - methods — lock onto local minima. This is exactly why production - `02_build_m6_dataset.py` uses a **multi-start tilt grid**. + tail.** On `m6_2BBN5XA4` / `m6_24V62QQM_20240820` (the two most + repositioned pairs — large inter-session translation; the required + rotations are only ~10° per the references, see the SynthMorph + investigation below) the top four hold 0.57–0.65, while niftyreg + (0.069 / 0.118) and synthmorph (0.061 / 0.111) — the single-start + block-matching and pretrained-MRI methods — lock onto local minima. This + is exactly why production `02_build_m6_dataset.py` uses a **multi-start + tilt grid**. - **NiftyReg (reg_aladin) is second-from-bottom**: single-start symmetric block matching is strong on moderate cases (`m6_23XIMO4M` 0.632; `m6_26TYQI3C` 0.609) but falls into the same tilt local minima @@ -248,12 +251,17 @@ the case): build). It clips the CT to [0,80] HU by default (per the A/B below) — neutral on the unclipped score scale (0.296→0.292 median). - **SynthMorph is weakest of the six on these cross-modality CT pairs** - (median 0.277 with its default CT clip; 0.241 unclipped): its `rigid` model - is trained on brain *MRI*, so the ~30° tilt defeats it on the repositioned - cases (0.06–0.11) even though it matches the intensity tools on the - cleanest (`m6_26TYQI3C` 0.590). The [0,80] brain-window clip of the CT is - its only genuine A/B win (+0.036 median, 7/9 cases), so it is the default. - It needs **no centroid init** (global deep method). + (median 0.277 with its default CT clip; 0.241 unclipped). The [0,80] + brain-window clip of the CT is its only genuine A/B win (+0.036 median, + 7/9 cases), so it is the default. It needs **no centroid init** (global + deep method). The low score is *not* a tilt effect: the required rotations + are ≤10.6° (the reference tools' own). The investigation (below) shows the + `rigid` matrix model fails gross-pose estimation case-by-case on the + repositioned pairs (19–66 mm translation errors; on `m6_24V62QQM_20240820` + its transform is worse than no registration at all), fails equally on an + in-domain MR→MR pair with a ~55° relative pose, and is not fixed by the + `affine` model — while the `joint` (affine + deformable) model recovers + all three tested failure cases to the best ImpactLoss of the six tools. - **Speed** (median wall, 3–4 cases in parallel, 24 threads, T1c-fixed): SimpleITK 10.8 s, elastix 11.7 s, NiftyReg 15.4 s (v2.1.1, GPU block matching), ANTs 16.7 s, SynthMorph 22.2 s (GPU forward pass), **IMPACT @@ -325,8 +333,128 @@ not clip.** Re-checked on the **unclipped score scale** (the protocol used by the table above): synthmorph keeps a genuine gain (0.241→0.277 median, 7/9 cases), niftyreg is neutral (0.296→0.292; its apparent +0.024 above is mostly the - clip-scoring effect). All rows of the main table are scored on the - unclipped warped CT, so pipelines remain comparable. +clip-scoring effect). All rows of the main table are scored on the + unclipped warped CT, so pipelines remain comparable. + +## SynthMorph low-performance investigation (2026-09-30) + +Question: *why* is `synthmorph` weakest (median MI 0.277; 0.06–0.11 on the +repositioned cases)? The earlier explanation ("trained on brain MRI, defeated +by the ~30° tilt") was checked and turned out to be wrong in its specifics. +Everything below was scored with the benchmark's own scoring +(`bench_registration.py::evaluate`); results in +`results/synthmorph_investigation/`, reproducible with +`scripts/investigate_synthmorph.py` (parts A–D, F) and +`scripts/investigate_synthmorph_mm.py` (part E). + +**0. It is not a pipeline bug.** The `rigid` run is deterministic and +reproduces the benchmark JSON exactly (e.g. `m6_26TYQI3C` MI 0.5904, gap +13.1/4.9). The read-back direction is right (8/9 `auto_inverted: false`, +verified by the `fix_orientation` skull-ring gate), and the matrix is rigid +(`det 0.999–1.001`, `‖RᵀR−I‖ ≤ 0.0015`). The lone `auto_inverted: true` +(`m6_2BBN5XA4`) is the gate choosing the *less bad* of two wrong directions — +a symptom of the failure, not its cause. + +**1. The failure is a gross-pose estimation failure, not a tilt failure.** +Each case's synthmorph transform was decomposed against the `simpleitk` +benchmark reference (MI 0.58–0.70 on all 9 cases; its rotations are +**1.4–10.6°**, not ~30°): + +| case | MI ref / identity / rigid | head-centre err (mm) | rot err (deg) | +|---|---|---|---| +| m6_22IRRKU6_20240115 | 0.647 / 0.250 / **0.277** | 28.2 | 6.2 | +| m6_23XIMO4M_20230222 | 0.654 / 0.246 / **0.616** | 1.9 | 2.2 | +| m6_24V62QQM_20230223 | 0.704 / 0.276 / **0.256** | 25.6 | 3.7 | +| m6_24V62QQM_20231017 | 0.597 / 0.331 / **0.484** | 3.8 | 17.7 | +| m6_24V62QQM_20240820 | 0.638 / 0.274 / **0.111** | 66.3 | 5.9 | +| m6_26TYQI3C_20221027 | 0.670 / 0.397 / **0.590** | 6.5 | 10.2 | +| m6_2BBN5XA4_20250411 | 0.579 / 0.046 / **0.061** | 61.1 | 22.1 | +| m6_2C4TKI6H_20240506 | 0.636 / 0.231 / **0.231** | 40.6 | 3.4 | +| m6_2CDNN7BH_20240102 | 0.652 / 0.197 / **0.323** | 19.4 | 3.0 | + +- 6/9 cases carry a 19–66 mm head-centre error (translation), and 3/9 a + 10–22° rotation error; the two "clean" failures are translation-only. +- On `m6_24V62QQM_20240820` the model's transform (MI 0.111) is **worse than + registering nothing at all** (identity MI 0.274); on `m6_2C4TKI6H` it is +exactly identity-level (0.231 = 0.231). The net gain over no-registration + is ≤0.03 on 5/9 cases. +- The error *vectors* are case-specific, not a systematic bias: the fraction + of the required network-space translation that goes uncorrected is ≈0 on + the successes (23XIMO4M 0.07, 20231017 −0.01), ≈1.0 on `m6_24V62QQM_20240820` + (the model outputs ≈ its network-space identity — the null that aligns the + two *volume* centres), and −1.7 on `m6_2BBN5XA4` (it moves the head 1.7× the + required amount, in the wrong direction). + +**2. Cross-modality is not the (only) cause.** An in-domain control — two +T1c scans of the *same* subject 8 months apart (`m6_24V62QQM`, both tight +head crops, no CT, no HU issue) — has a true relative pose of **~55–57°** +(benchmark `simpleitk` mode 55.3°, production multi-start tilt grid 56.6°, +head centres agreeing to 1.2 mm) and SynthMorph `rigid` fails it equally: +18.5° rotation error, 18 mm head-centre error, MI 0.361 vs the reference's +0.689 (ImpactLoss 0.140 vs 0.114). + +**3. The `affine` matrix model does not help.** Same 9 cases, same budget: +MI within ±0.025 of `rigid` on 8/9 cases (`m6_23XIMO4M`: −0.096, the one case +where rigid succeeded); the linear part is a near-pure rotation +(scale 0.98–1.03, anisotropy ≤2%) — the extra 6 DOF are essentially unused. + +**4. The `joint` (affine + deformable) model recovers the failures.** Run +directly in Docker (`-m joint`, outside the rigid contract; the tool's own +`-o` resampling warps the CT onto the T1c grid, scored as-is) on the three +worst cases: + +| case | MI rigid → joint | ImpactLoss joint (vs best intensity tool) | +|---|---|---| +| m6_24V62QQM_20230223 | 0.256 → **0.529** | **0.133** (simpleitk 0.147) | +| m6_24V62QQM_20240820 | 0.111 → **0.481** | **0.132** (antspyx 0.157) | +| m6_2BBN5XA4_20250411 | 0.061 → **0.433** | **0.129** (antspyx 0.162) | + +The model family's anatomy is fine — its *semantic* alignment on the two +grossly repositioned cases is the **best of the six tools** — but the matrix +head (6- or 12-DOF) cannot estimate the gross inter-session pose, and the +deformable field absorbs what it misses (~40 s/case on the GPU). + +**5. Pipeline-side suspects were ruled out.** + +- **Network-space extent truncation: minor.** Each volume is resampled into + its own 256 mm isotropic box (centred on its own volume centre, per the + image code); with the head placed in CT space via the reference transform, + only 0–6.5% of head voxels fall outside the CT box (0% T1c). Re-centring + the CT volume on the head (crop to head bbox + 15 mm) changes *nothing* + (ΔMI ≤ 0.012, Δhead-centre-err ≤ 1.2 mm on all 9 cases) — the model is not + anchored to the volume centre. +- **Min-max normalization: a real, quantified cross-modality penalty.** The + network min-max-normalizes each volume; over the raw CT's full HU range + (−3024…+3071, air/metal) the brain (0–80 HU) collapses to a near-flat + **0.496–0.509** band (air sits at ~0.33, the skull at 0.63–1.0), so the + parenchyma carries almost no internal contrast; the T1c brain spans + ~0–0.6 (p95 0.26–0.58), so the raw CT reads as a flat-gray blob with a + bright ring. The [0,80] clip restores full-range brain contrast (this is + the mechanism behind the CT-clip A/B win, +0.036–0.077 MI). Clipping does + not make the model domain-typical, which is why it only partially helps. + +**Conclusion.** SynthMorph's matrix models are a *local pose prior*: on the +in-distribution cases (small relative pose) they estimate the M6 alignment +well (23XIMO4M 1.9 mm / 2.2°, 20231017 3.8 mm / 17.7°, 26TYQI3C 6.5 mm / +10.2° residual), but the M6 inter-session repositioning (tens of mm +translation and, in the same-patient control, a 55° world-frame rotation) is +outside the range the pretrained model (trained on same-subject MRI pairs) +reliably inverts, and the output is +then case-specifically wrong — do-nothing, half-done, or opposite. It is a +valid data point for how far a pretrained-MRI deep matrix model generalizes +out-of-domain, and it should not be expected to compete with the intensity +tools on repositioned CT pairs in rigid mode. If a SynthMorph variant is +wanted for these data, the `joint` model is the one that works (at the cost +of leaving the rigid-transform contract). + +Reproducing (needs the `freesurfer/synthmorph` image + GPU; ~15 min): + +```bash +source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 +python scripts/investigate_synthmorph.py --parts ABCD --workers 3 # A ref+rigid, B affine, C joint, D network-space +python scripts/investigate_synthmorph.py --parts F --workers 3 # re-centred CT probe (needs ABCD first) +python scripts/investigate_synthmorph_mm.py --gpu 0 # in-domain MR->MR control +``` ## Reproducing @@ -581,10 +709,11 @@ docker pull vboussot/elastix_impact cost; that is out of scope for the head-to-head but is the config the dataset actually ships with. - `synthmorph` is deterministic (one fixed pretrained model, no random - subsampling), so its per-case numbers are stable run to run — unlike ANTs / - IMPACT. Its low MI (median 0.277 with the default CT clip; 0.06–0.11 on the - repositioned cases) is a genuine cross-modality / tilt weakness, not - variance. + subsampling), so its per-case numbers are stable run to run — unlike ANTs / + IMPACT. Its low MI (median 0.277 with the default CT clip; 0.06–0.11 on the + repositioned cases) is a genuine gross-pose-estimation weakness of the + pretrained matrix model (case-by-case, not variance) — see the + "SynthMorph low-performance investigation" section above. - `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 @@ -597,11 +726,19 @@ re-derived per toolkit. Under the shared fixed=T1c / moving=CT rule the raw runs report `auto_inverted: False` for those; `niftyreg`'s `-aff` file comes out the other way and the gate inverts it (8/9 `auto_inverted: True`). - **SynthMorph specifics** (Docker `freesurfer/synthmorph`, TensorFlow, GPU): - - **Cross-modality / tilt weakness.** The `rigid` model is trained on brain - *MRI*, so on these CT↔T1c pairs it recovers translation but is repeatedly - defeated by the ~30° inter-session tilt (7/9 cases have `std > 15 mm`). - It is a legitimate data point for how far a pretrained MRI deep model - generalizes out-of-domain to CT — not a bug. + - **Gross-pose estimation weakness (measured, not a bug).** The `rigid` + matrix model is trained on brain *MRI* (same-subject pairs, small + relative poses). On these repositioned CT↔T1c pairs it estimates the + translation wrongly on 6/9 cases (19–66 mm head-centre error; the two + worst, `m6_2BBN5XA4` / `m6_24V62QQM_20240820`, at 61 / 66 mm) and the + rotation on 3/9 (10–22°) — and the required rotations are only 1.4–10.6°, + so it is *not* a "defeated by tilt" story; on `m6_24V62QQM_20240820` + its transform scores below the unregistered (identity) one. An in-domain + MR→MR control (same subject, ~55° relative pose) fails equally, and the + `affine` model changes nothing, while the `joint` (affine + deformable) + model recovers all three tested failures to the best ImpactLoss of the + six tools. Full evidence: the "SynthMorph low-performance investigation" + section. - **No centroid init needed.** Being a global deep method it recovers the gross offset on its own; the `-i` initial-transform path is not used (which also sidesteps this image's float64/float32 bug on the `-i` code diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index 2d6bf6f..819afc9 100644 --- a/scripts/bench_registration.py +++ b/scripts/bench_registration.py @@ -27,14 +27,14 @@ A/B compares the same input to both the toolkit and the MI score. Because some toolkits are not pip-installable into every `longitudinal-313` env (antspyx needs a cp310-cp313 wheel, see -REGISTRATION_BENCHMARK.md), the benchmark is env-portable: run it under +REGISTRATION_BENCHMARK_RIGID.md), the benchmark is env-portable: run it under whichever conda env has a given toolkit installed, and it records only the pipelines that are importable there. elastix runs in the main env through the ITK-Elastix python API (pip ``itk`` + ``itk-elastix``); niftyreg runs natively (``conda install -c conda-forge niftyreg`` provides ``reg_aladin``); ``impact`` runs the ``vboussot/elastix_impact`` Docker image (GPU; the rigid MIND + TotalSegmentator static config, needs Docker + the pulled image + the -model assets — see REGISTRATION_BENCHMARK.md); ``synthmorph`` runs in a Docker +model assets — see REGISTRATION_BENCHMARK_RIGID.md); ``synthmorph`` runs in a Docker image (``freesurfer/synthmorph``; needs Docker + the pulled image). The `--merge` mode then combines the per-pipeline JSON files (accumulated across envs) into @@ -45,7 +45,7 @@ Usage (repo root): source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 python scripts/bench_registration.py --pipelines simpleitk elastix antspyx niftyreg impact - # synthmorph (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK.md + # synthmorph (GPU, Docker) — one-time setup in REGISTRATION_BENCHMARK_RIGID.md python scripts/bench_registration.py --pipelines synthmorph # any pipeline that is not importable in the current env is recorded as @@ -290,6 +290,28 @@ def make_qa_png(warped_ct, mr, out_png, title, block=24, ct_brain_window=(0.0, 8 def _call(pipeline, ct, t1c, tv, args): + if args.deformable: + if pipeline == "simpleitk": + return RP.register_simpleitk_deformable(ct, t1c, threads=args.threads) + if pipeline == "elastix": + return RP.register_elastix_deformable(ct, t1c, workdir=args.elastix_dir, + threads=args.threads) + if pipeline == "antspyx": + return RP.register_antspyx_deformable(ct, t1c) + if pipeline == "synthmorph": + return RP.register_synthmorph_deformable(ct, t1c, + gpu=args.synthmorph_gpu, + threads=args.threads, + extent=args.synthmorph_extent) + if pipeline == "impact": + return RP.register_impact_deformable(ct, t1c, gpu=args.impact_gpu, + iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS, + samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES) + if pipeline == "niftyreg": + return RP.register_niftyreg_deformable(ct, t1c, threads=args.threads, + backend=args.niftyreg_backend, + gpu=args.niftyreg_gpu) + raise KeyError(pipeline) if pipeline == "simpleitk": return RP.register_simpleitk(ct, t1c, threads=args.threads, mode=args.simpleitk_mode) if pipeline == "elastix": @@ -444,7 +466,7 @@ def print_table(rows): print("=" * 96) -def write_plot(rows, out_png): +def write_plot(rows, out_png, label="rigid"): try: import matplotlib matplotlib.use("Agg") @@ -466,6 +488,10 @@ def write_plot(rows, out_png): ("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")]): + # wide time ranges (deformable: simpleitk ~100x the others) -> log axis + nz = [v for v in vals if v > 0] + if len(nz) >= 2 and max(nz) / min(nz) > 30: + ax.set_yscale("log") bars = ax.bar(range(len(vals)), vals, color=colors[0]) for b, v in zip(bars, vals): ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()), @@ -474,7 +500,7 @@ def write_plot(rows, out_png): ax.set_ylabel(ylab) ax.set_xticks(range(len(vals))) ax.set_xticklabels(labels, fontsize=9) - plt.suptitle("CT->T1c registration benchmark (M6, cross-modality rigid)", fontsize=11) + plt.suptitle(f"CT->T1c registration benchmark (M6, cross-modality {label})", fontsize=11) fig.tight_layout() fig.savefig(out_png, dpi=120) plt.close(fig) @@ -493,15 +519,30 @@ def merge(args): raise SystemExit(f"no pipeline results found in {out_dir}") print(f"[merge] {len(rows)} pipeline result files in {out_dir}") print_table(rows) + deformable = (os.path.basename(os.path.normpath(args.out_dir)) + == "bench_registration_deformable") + if deformable: + protocol = ("deformable CT->T1c; cross-modality; each pipeline runs its " + "toolkit's native deformable workflow (linear init + deformable " + "stage, or a single deformable stage where that is native); the " + "result is rasterised to the T1c-grid displacement field of the " + "unified contract (T1c->CT); alignment scored by achieved Mattes MI " + "(nats, warped-CT ~ MR over the MR head, higher = better) and the " + "IMPACT semantic loss (vboussot/ImpactLoss, TS/M730 features, " + "warped-CT ~ MR over the full volumes, lower = better); reg time " + "excludes I/O+scoring; per-case JSON retains the production " + "skull-ring gap / acceptance QA") + else: + protocol = ("rigid CT->T1c; cross-modality; alignment scored by achieved " + "Mattes MI (nats, warped-CT ~ MR over the MR head, higher = " + "better) and the IMPACT semantic loss (vboussot/ImpactLoss, " + "TS/M730 features, warped-CT ~ MR over the full volumes, lower " + "= better); reg time excludes I/O+scoring; per-case JSON retains " + "the production skull-ring gap / acceptance QA; best-of-T vs T^-1 " + "orientation gate for the elastix-style toolkits") summary = { "timestamp": datetime.now().isoformat(timespec="seconds"), - "protocol": ("rigid CT->T1c; cross-modality; alignment scored by achieved Mattes MI " - "(nats, warped-CT ~ MR over the MR head, higher = better) and the " - "IMPACT semantic loss (vboussot/ImpactLoss, TS/M730 features, " - "warped-CT ~ MR over the full volumes, lower = better); reg time " - "excludes I/O+scoring; per-case JSON retains the production skull-ring " - "gap / acceptance QA; best-of-T vs T^-1 orientation gate for the " - "elastix-style toolkits"), + "protocol": protocol, "n_cases": next((r["n_cases"] for r in rows if r.get("available")), None), "case_keys": sorted({c["key"] for r in rows if r.get("available") for c in r.get("cases", [])}), "results": rows, @@ -510,7 +551,8 @@ def merge(args): with open(sf, "w") as f: json.dump(summary, f, indent=1) 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"), + label="deformable" if deformable else "rigid") # markdown table for the report md = ["| pipeline | available | toolkit | cases | ok | err | s/reg (median) | ImpactLoss (median) | MI (nats, median) |", "|---|---|---|---|---|---|---|---|---|"] @@ -533,6 +575,10 @@ def merge(args): def parse_args(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--deformable", action="store_true", + help="run the deformable variants of the pipelines " + "(each toolkit's native deformable workflow; results go " + "to results/bench_registration_deformable by default)") ap.add_argument("--pipelines", default=" ".join(RP.ALL_PIPELINES), help="space-separated subset of: " + " ".join(RP.ALL_PIPELINES)) ap.add_argument("--n-cases", type=int, default=10, help="number of M6 test cases (default 10)") @@ -594,6 +640,9 @@ def parse_args(): def main(): args = parse_args() + if args.deformable: + if args.out_dir == os.path.join(ROOT, "results", "bench_registration"): + args.out_dir = os.path.join(ROOT, "results", "bench_registration_deformable") if args.merge: merge(args) return diff --git a/scripts/bench_registration_deformable.py b/scripts/bench_registration_deformable.py new file mode 100644 index 0000000..420b288 --- /dev/null +++ b/scripts/bench_registration_deformable.py @@ -0,0 +1,439 @@ +"""Benchmark the 6 DEFORMABLE CT->T1c registration pipelines on N M6 cases. + +Companion to ``scripts/bench_registration.py`` (rigid). Pipelines: + simpleitk | elastix | antspyx | synthmorph | impact | niftyreg +(each = its toolkit's native deformable refinement initialized from the +toolkit's own rigid result; see src/registration_pipelines_deformable.py +and REGISTRATION_BENCHMARK_DEFORMABLE.md). + +Headline scores (identical for every pipeline): + - mi_ct_mr : achieved Mattes MI (nats), warped-CT ~ MR over the MR head + (higher = better). Also computed for the pipeline's own + rigid base (mi_base) -> d_mi = the deformable gain. + - impact_loss : IMPACT semantic loss (TS/M730 features), warped-CT ~ MR + over the full volumes (lower = better). Also mi/IL for + the base (il_base). + - jac_* : Jacobian-determinant statistics of the FINAL deformable + transform over the MR head (finite differences on the + T1c grid): p1/p50/mean + fraction < 0 (folded), < 0.3, + < 0.5. A healthy diffeomorphic refinement has J ~ 1. + - jac_gtv_* : same Jacobian statistics restricted to the warped-GTV + region (does the deformation fold the tumor?). + - disp_base_* : |T_final(x) - T_base(x)| over the MR head (p50/p95/max, + mm) -- the deformable increment. + - tv_vol_ratio : warped-GTV volume final / base (a diffeomorphism + preserves volume, ~1). + - reg_wall_s : wall-clock of the deformable stage (excludes I/O, the + rigid base and scoring). +Production QA (skull-ring gap, accepted) is still computed per case and kept +in the JSON. + +Usage (repo root): + source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal-313 + python scripts/bench_registration_deformable.py --pipelines "simpleitk elastix antspyx niftyreg" + python scripts/bench_registration_deformable.py --pipelines "impact synthmorph" + python scripts/bench_registration_deformable.py --merge +""" +import argparse +import importlib.util +import json +import os +import sys +import time +from datetime import datetime +from concurrent.futures import ProcessPoolExecutor, as_completed + +ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, ROOT) + +import numpy as np +import SimpleITK as sitk + +from src.common import DATA, d, read_nii_img, head_mask_from_image +from src import registration_pipelines as RP +from src import registration_pipelines_deformable as RPD + +_spec = importlib.util.spec_from_file_location("br", os.path.join(ROOT, "scripts", "bench_registration.py")) +BR = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(BR) + +JAC_STEP_MM = 1.0 # finite-difference half-step for the Jacobian estimate + + +# --------------------------------------------------------------------------- # +# deformable-specific scoring +# --------------------------------------------------------------------------- # + +def _dense_H(T, ref): + """Absolute mapping H(x) = T(x) at every voxel of ``ref`` (z,y,x,3).""" + return RPD.densify_sitk(T, ref)[2] + + +def jacobian_of_T(T, ref, step_mm=JAC_STEP_MM): + """det dT/dx of an arbitrary contract transform over ``ref``'s grid. + + Central finite differences with world-axis steps of ~``step_mm`` on the + densified absolute mapping (linear resampling of coordinate volumes); + border voxels (no both-side neighbour) are NaN. Returns J (z,y,x).""" + H = RPD.densify_sitk(T, ref)[2] + sp = np.array(ref.GetSpacing()) # (x, y, z); array axis i has spacing sp[2-i] + off = np.array([int(np.clip(np.round(step_mm / sp[2 - i]), 1, None)) + for i in range(3)]) + H = np.ascontiguousarray(H, dtype='float64') + dH = [] + for i in range(3): + d = np.full(H.shape, np.nan) + sl_p = [slice(None)] * 3; sl_m = [slice(None)] * 3; sl_mid = [slice(None)] * 3 + sl_p[i] = slice(2 * off[i], None) + sl_m[i] = slice(None, -2 * off[i]) + sl_mid[i] = slice(off[i], -off[i]) + d[tuple(sl_mid)] = (H[tuple(sl_p)] - H[tuple(sl_m)]) / (2.0 * off[i] * sp[2 - i]) + dH.append(d) + # array axis i <-> world axis (2-i): axis 2 = x, axis 1 = y, axis 0 = z. + # J = det[dH/dx, dH/dy, dH/dz] + A = dH[2]; B = dH[1]; C = dH[0] + J = (A[..., 0] * (B[..., 1] * C[..., 2] - B[..., 2] * C[..., 1]) + - A[..., 1] * (B[..., 0] * C[..., 2] - B[..., 2] * C[..., 0]) + + A[..., 2] * (B[..., 0] * C[..., 1] - B[..., 1] * C[..., 0])) + return J + + +def _stats(J, mask): + v = J[mask] + v = v[np.isfinite(v)] + if v.size < 50: + return None + return { + "p1": round(float(np.percentile(v, 1)), 4), + "p50": round(float(np.median(v)), 4), + "mean": round(float(v.mean()), 4), + "frac_neg": round(float((v < 0).mean()), 5), + "frac_lt03": round(float((v < 0.3).mean()), 5), + "frac_lt05": round(float((v < 0.5).mean()), 5), + "nvox": int(v.size), + } + + +def deformable_metrics(ct, t1c, tv, T_final, T_base): + """Jacobian / displacement-increment / GTV-preservation metrics of a + deformable result (final transform vs its rigid base).""" + ha = sitk.GetArrayFromImage(t1c).astype('float32') + hm = head_mask_from_image(ha) + J = jacobian_of_T(T_final, t1c) + jac_head = _stats(J, hm) + # warped GTV mask (final) + wtv = sitk.GetArrayFromImage( + sitk.Resample(tv, t1c, T_final, sitk.sitkNearestNeighbor, 0.0)) + gtv = wtv > 0.5 + if gtv.sum() > 50: + gtv = gtv & hm + jac_gtv = _stats(J, gtv) + else: + jac_gtv = None + # displacement increment |T_final - T_base| over the head + Hf = _dense_H(T_final, t1c) + Hb = _dense_H(T_base, t1c) + inc = np.linalg.norm(Hf - Hb, axis=-1) + inc_stats = { + "p50": round(float(np.median(inc[hm])), 2), + "p95": round(float(np.percentile(inc[hm], 95)), 2), + "max": round(float(inc[hm].max()), 2), + } + # GTV volume preservation (final vs base) + wtv_b = sitk.GetArrayFromImage( + sitk.Resample(tv, t1c, T_base, sitk.sitkNearestNeighbor, 0.0)) + nv_f = int((wtv > 0.5).sum()); nv_b = int((wtv_b > 0.5).sum()) + tv_vol_ratio = round(nv_f / nv_b, 3) if nv_b > 0 else None + return {"jac_head": jac_head, "jac_gtv": jac_gtv, + "disp_base_mm": inc_stats, "tv_vol_ratio": tv_vol_ratio, + "gtv_nvox_final": nv_f, "gtv_nvox_base": nv_b} + + +# --------------------------------------------------------------------------- # +# case worker +# --------------------------------------------------------------------------- # + +def _call(pipeline, ct, t1c, args): + if pipeline == "simpleitk": + return RPD.register_simpleitk_deformable(ct, t1c, threads=args.threads) + if pipeline == "elastix": + return RPD.register_elastix_deformable(ct, t1c, workdir=args.elastix_dir, + threads=args.threads) + if pipeline == "antspyx": + return RPD.register_antspyx_deformable(ct, t1c, threads=args.threads) + if pipeline == "niftyreg": + return RPD.register_niftyreg_deformable(ct, t1c, threads=args.threads, + gpu=args.niftyreg_gpu, + backend=args.niftyreg_backend) + if pipeline == "synthmorph": + return RPD.register_synthmorph_deformable(ct, t1c, gpu=args.synthmorph_gpu, + extent=args.synthmorph_extent) + if pipeline == "impact": + return RPD.register_impact_deformable( + ct, t1c, gpu=args.impact_gpu, + iters=args.impact_iters or RP.IMPACT_DEFAULT_ITERS, + samples=args.impact_samples or RP.IMPACT_DEFAULT_SAMPLES, + threads=args.threads) + raise KeyError(pipeline) + + +def _run_case(pipeline, case, args): + key = case["key"] + ct = read_nii_img(case["ct"]) + t1c = read_nii_img(case["img"]) + tv = read_nii_img(case["label"]) + out = {"key": key, "subject": case["subject"], "date": case["date"]} + if not BR.ct_complete(ct): + out.update({"skull_ring_gap_mm": None, "tv_nvox": 0, "tv_in_ct_frac": None, + "impact_loss": None, "alignment": "ct_incomplete", + "accepted": False, + "note": "CT head incomplete (data issue, not a pipeline failure)"}) + return out + t0 = time.time() + T_final, T_base, meta = _call(pipeline, ct, t1c, args) + out["reg_wall_s"] = round(time.time() - t0, 1) + out.update(meta) + # final + base alignment scores + out.update(BR.evaluate(ct, t1c, tv, T_final)) + out["base"] = BR.evaluate(ct, t1c, tv, T_base) + out["d_mi"] = (round(out["mi_ct_mr"] - out["base"]["mi_ct_mr"], 4) + if out["mi_ct_mr"] is not None and out["base"]["mi_ct_mr"] is not None + else None) + out["d_il"] = (round(out["impact_loss"] - out["base"]["impact_loss"], 4) + if out["impact_loss"] is not None and out["base"]["impact_loss"] is not None + else None) + out.update(deformable_metrics(ct, t1c, tv, T_final, T_base)) + return out + + +def run_pipeline(pipeline, cases, args): + avail, note = RPD.availability().get(pipeline, (False, "unknown")) + if not avail: + row = {"pipeline": pipeline, "available": False, "reason": note, + "n_cases": len(cases), "cases": []} + print(f"[{pipeline}] UNAVAILABLE: {note}") + return row + print(f"[{pipeline}] available ({note}); registering {len(cases)} cases " + f"with {args.workers} workers") + results, errors = [], [] + with ProcessPoolExecutor(max_workers=args.workers) as ex: + futs = {ex.submit(_run_case, pipeline, c, args): c for c in cases} + for fu in as_completed(futs): + c = futs[fu] + try: + results.append(fu.result()) + r = results[-1] + reg = f"reg={r.get('reg_wall_s')}s" if r.get("reg_wall_s") is not None else "incomplete-CT" + print(f" {c['key']}: mi={r.get('mi_ct_mr')} (base {r.get('base', {}).get('mi_ct_mr')}) " + f"il={r.get('impact_loss')} " + f"jac_p1={(r.get('jac_head') or {}).get('p1')} {reg}", flush=True) + except Exception as e: # noqa + errors.append({"key": c["key"], "error": repr(e)}) + print(f" {c['key']}: ERROR {e!r}", flush=True) + results.sort(key=lambda x: x["key"]) + registered = [r for r in results if r.get("reg_wall_s") is not None] + incomplete = [r for r in results if r.get("alignment") == "ct_incomplete"] + accepted = [r for r in registered if r.get("accepted")] + + def _med(vals, nd=4): + vals = [v for v in vals if v is not None] + return round(float(np.median(vals)), nd) if vals else None + + row = { + "pipeline": pipeline, "available": True, + "toolkit": (registered or [{}])[0].get("toolkit"), + "metric": (registered or [{}])[0].get("metric"), + "transform": (registered or [{}])[0].get("transform"), + "n_cases": len(cases), "n_registered": len(registered), + "n_ct_incomplete": len(incomplete), "n_error": len(errors), + "n_accepted": len(accepted), + "accepted_rate": round(len(accepted) / max(len(cases), 1), 3), + "reg_wall_s": { + "median": _med([r["reg_wall_s"] for r in registered], 1), + "mean": round(float(np.mean([r["reg_wall_s"] for r in registered])), 1) if registered else None, + "max": round(float(np.max([r["reg_wall_s"] for r in registered])), 1) if registered else None, + }, + "mi_ct_mr": {"median": _med([r.get("mi_ct_mr") for r in registered]), + "mean": round(float(np.mean([r["mi_ct_mr"] for r in registered if r.get("mi_ct_mr") is not None])), 4) if any(r.get("mi_ct_mr") is not None for r in registered) else None}, + "mi_base": {"median": _med([r.get("base", {}).get("mi_ct_mr") for r in registered])}, + "d_mi": {"median": _med([r.get("d_mi") for r in registered])}, + "impact_loss": {"median": _med([r.get("impact_loss") for r in registered]), + "mean": round(float(np.mean([r["impact_loss"] for r in registered if r.get("impact_loss") is not None])), 4) if any(r.get("impact_loss") is not None for r in registered) else None}, + "il_base": {"median": _med([r.get("base", {}).get("impact_loss") for r in registered])}, + "d_il": {"median": _med([r.get("d_il") for r in registered])}, + "jac_p1": {"median": _med([(r.get("jac_head") or {}).get("p1") for r in registered])}, + "jac_frac_neg": {"median": _med([(r.get("jac_head") or {}).get("frac_neg") for r in registered], 5)}, + "jac_frac_lt05": {"median": _med([(r.get("jac_head") or {}).get("frac_lt05") for r in registered], 5)}, + "jac_gtv_p1": {"median": _med([(r.get("jac_gtv") or {}).get("p1") for r in registered])}, + "disp_base_p95_mm": {"median": _med([(r.get("disp_base_mm") or {}).get("p95") for r in registered], 2)}, + "tv_vol_ratio": {"median": _med([r.get("tv_vol_ratio") for r in registered])}, + "cases": results, "errors": errors, + "timestamp": datetime.now().isoformat(timespec="seconds"), + } + return row + + +# --------------------------------------------------------------------------- # +# reporting +# --------------------------------------------------------------------------- # + +def _fmt(v, spec="{:g}", default="-"): + return default if v is None else spec.format(v) + + +def print_table(rows): + order = [p for p in RPD.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)] + print("\n" + "=" * 118) + print(f"{'pipeline':10s} {'avail':5s} {'n':>3s} {'ok':>3s} {'err':>3s} {'s/reg':>7s} " + f"{'MI fin':>7s} {'MI base':>8s} {'dMI':>7s} {'IL fin':>7s} {'IL base':>8s} " + f"{'jacP1':>7s} {'J<0':>7s} {'J<.5':>6s} {'incr p95':>9s}") + print("-" * 118) + for p in order: + r = next(x for x in rows if x["pipeline"] == p) + if not r.get("available"): + print(f"{p:10s} {'no':5s} {r['n_cases']:>3d} {0:>3d} {'-':>3s} [{r['reason'][:60]}]") + continue + print(f"{p:10s} {'yes':5s} {r['n_cases']:>3d} {r.get('n_registered', 0):>3d} " + f"{r.get('n_error', 0):>3d} {_fmt((r.get('reg_wall_s') or {}).get('median'), '{:.0f}'):>7s} " + f"{_fmt((r.get('mi_ct_mr') or {}).get('median'), '{:.3f}'):>7s} " + f"{_fmt((r.get('mi_base') or {}).get('median'), '{:.3f}'):>8s} " + f"{_fmt((r.get('d_mi') or {}).get('median'), '{:+.3f}'):>7s} " + f"{_fmt((r.get('impact_loss') or {}).get('median'), '{:.3f}'):>7s} " + f"{_fmt((r.get('il_base') or {}).get('median'), '{:.3f}'):>8s} " + f"{_fmt((r.get('jac_p1') or {}).get('median'), '{:.3f}'):>7s} " + f"{_fmt((r.get('jac_frac_neg') or {}).get('median'), '{:.4f}'):>7s} " + f"{_fmt((r.get('jac_frac_lt05') or {}).get('median'), '{:.3f}'):>6s} " + f"{_fmt((r.get('disp_base_p95_mm') or {}).get('median'), '{:.1f}'):>9s}") + print("=" * 118) + + +def merge(args): + out_dir = args.out_dir + rows = [] + for p in RPD.ALL_PIPELINES: + fp = os.path.join(out_dir, p + ".json") + if os.path.exists(fp): + with open(fp) as f: + rows.append(json.load(f)) + if not rows: + raise SystemExit(f"no pipeline results found in {out_dir}") + print(f"[merge] {len(rows)} pipeline result files in {out_dir}") + print_table(rows) + summary = { + "timestamp": datetime.now().isoformat(timespec="seconds"), + "protocol": ("deformable CT->T1c (each toolkit's native deformable refinement " + "initialized from its own rigid result); alignment scored by " + "achieved Mattes MI (nats, warped-CT ~ MR over the MR head, higher " + "= better) and the IMPACT semantic loss (TS/M730, full volumes, " + "lower = better), plus Jacobian-determinant statistics of the " + "final transform (finite differences on the T1c grid), the " + "deformable displacement increment vs the rigid base, and GTV " + "volume preservation; same contract / shared rules as the rigid " + "benchmark (full volumes, T1c fixed / CT moving, no mask)"), + "n_cases": next((r["n_cases"] for r in rows if r.get("available")), None), + "case_keys": sorted({c["key"] for r in rows if r.get("available") + for c in r.get("cases", [])}), + "results": rows, + } + sf = os.path.join(out_dir, "summary.json") + with open(sf, "w") as f: + json.dump(summary, f, indent=1) + print("summary:", sf) + md = ["| pipeline | toolkit | cases | ok | s/reg (med) | MI (med) | MI base (med) | ΔMI (med) | ImpactLoss (med) | IL base (med) | ΔIL (med) | jac p1 (med) | frac J<0 (med) | frac J<0.5 (med) | incr p95 (med, mm) | tv vol ratio (med) |", + "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|"] + for r in rows: + if not r.get("available"): + md.append(f"| {r['pipeline']} | - | {r['n_cases']} | 0 | - | - | - | - | - | - | - | - | - | - | - | - | ({r['reason']}) |") + continue + md.append( + f"| {r['pipeline']} | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | " + f"{_fmt((r.get('reg_wall_s') or {}).get('median'), '{:.0f}')} | " + f"{_fmt((r.get('mi_ct_mr') or {}).get('median'), '{:.3f}')} | " + f"{_fmt((r.get('mi_base') or {}).get('median'), '{:.3f}')} | " + f"{_fmt((r.get('d_mi') or {}).get('median'), '{:+.3f}')} | " + f"{_fmt((r.get('impact_loss') or {}).get('median'), '{:.3f}')} | " + f"{_fmt((r.get('il_base') or {}).get('median'), '{:.3f}')} | " + f"{_fmt((r.get('d_il') or {}).get('median'), '{:+.3f}')} | " + f"{_fmt((r.get('jac_p1') or {}).get('median'), '{:.3f}')} | " + f"{_fmt((r.get('jac_frac_neg') or {}).get('median'), '{:.4f}')} | " + f"{_fmt((r.get('jac_frac_lt05') or {}).get('median'), '{:.3f}')} | " + f"{_fmt((r.get('disp_base_p95_mm') or {}).get('median'), '{:.1f}')} | " + f"{_fmt((r.get('tv_vol_ratio') or {}).get('median'), '{:.3f}')} |") + mf = os.path.join(out_dir, "summary.md") + with open(mf, "w") as f: + f.write("\n".join(md) + "\n") + print("markdown:", mf) + + +def parse_args(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--pipelines", default=" ".join(RPD.ALL_PIPELINES)) + ap.add_argument("--n-cases", type=int, default=10) + ap.add_argument("--sample", choices=["first", "random"], default="first") + ap.add_argument("--seed", type=int, default=0) + ap.add_argument("--keys", default=None) + ap.add_argument("--workers", type=int, default=4) + ap.add_argument("--threads", type=int, default=24) + ap.add_argument("--out-dir", default=os.path.join(ROOT, "results", "bench_registration_deformable")) + ap.add_argument("--elastix-dir", default=None) + ap.add_argument("--synthmorph-gpu", type=int, default=None) + ap.add_argument("--synthmorph-extent", type=int, default=256) + ap.add_argument("--impact-gpu", type=int, default=None) + ap.add_argument("--impact-iters", type=int, default=None) + ap.add_argument("--impact-samples", type=int, default=None) + ap.add_argument("--niftyreg-backend", choices=["auto", "conda", "docker"], default="auto") + ap.add_argument("--niftyreg-gpu", type=int, default=None) + ap.add_argument("--merge", action="store_true") + args = ap.parse_args() + return args + + +def main(): + args = parse_args() + if args.merge: + merge(args) + return + args.elastix_dir = args.elastix_dir or d("data/m6reg/elastix") + pipelines = [p.lower() for p in args.pipelines.split()] + unknown = set(pipelines) - set(RPD.ALL_PIPELINES) + if unknown: + raise SystemExit(f"unknown pipeline(s): {unknown}") + if args.keys: + m6 = RP._load_m6() + labeled, _ = m6.discover() + bykey = {r["key"]: r for r in labeled} + missing = [k for k in args.keys.split(",") if k not in bykey] + if missing: + raise SystemExit(f"keys not found in m6 labeled: {missing}") + cases = [bykey[k.strip()] for k in args.keys.split(",")] + print(f"[cases] using {len(cases)} explicit keys") + else: + cases = BR.select_cases(args) + print(f"\n[bench-def] pipelines={pipelines} n_cases={len(cases)} workers={args.workers} " + f"threads={args.threads} out_dir={args.out_dir}") + rows = [] + for p in pipelines: + try: + rows.append(run_pipeline(p, cases, args)) + except RP.PipelineUnavailable as e: + rows.append({"pipeline": p, "available": False, "reason": str(e), + "n_cases": len(cases), "cases": []}) + print(f"[{p}] UNAVAILABLE: {e}") + for r in rows: + os.makedirs(args.out_dir, exist_ok=True) + fp = os.path.join(args.out_dir, r["pipeline"] + ".json") + with open(fp, "w") as f: + json.dump(r, f, indent=1) + print(f"[{r['pipeline']}] wrote {fp}") + print_table(rows) + avail = RP.availability() + with open(os.path.join(args.out_dir, "_availability.json"), "w") as f: + json.dump({"timestamp": datetime.now().isoformat(timespec="seconds"), + "env_python": sys.version.split()[0], "exec": sys.executable, + "availability": avail}, f, indent=1) + print("[bench-def] done. Re-run with --merge to combine.") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/scripts/investigate_synthmorph.py b/scripts/investigate_synthmorph.py new file mode 100644 index 0000000..65df700 --- /dev/null +++ b/scripts/investigate_synthmorph.py @@ -0,0 +1,601 @@ +"""Investigate the low performance of the ``synthmorph`` pipeline on the M6 +CT->T1c rigid benchmark (see REGISTRATION_BENCHMARK_RIGID.md). + +The benchmark median MI is 0.277 for synthmorph vs ~0.65 for the intensity +tools, with the ~30 deg repositioned cases collapsing to 0.06-0.11. This +script decomposes *why*, scoring everything with the exact benchmark scoring +(``scripts/bench_registration.py::evaluate``): + +A. Transform decomposition -- for each case, run synthmorph ``rigid`` (the + benchmark config; deterministic, so it reproduces the bench JSON) and a + ``simpleitk`` benchmark reference, and report per case: + - achieved MI / ImpactLoss / skull-ring gap for both, + - rotation angle + axis of each transform (world RAS frame), + - the *relative* rotation error angle(R_ref^T R_sm), + - the head-centre placement error |T_sm(c) - T_ref(c)| (c = T1c head + centroid), + - the identity-transform score (the "do nothing" baseline), + - a numeric rigidity check of the synthmorph matrix (R^T R ~ I, det = 1). + +B. Affine model -- run synthmorph ``affine`` (12-DOF matrix model) on the + same cases and score it: does the bigger matrix family recover the tilt + that rigid does not? The linear part is decomposed into rotation / + scale / shear (polar). + +C. Joint (affine + deformable) upper bound -- for the worst cases, run + ``-m joint`` directly in Docker (outside the rigid contract: the output is + a displacement field, applied with the tool's own ``-o`` resampling) and + score the warped CT image. Answers: is the low score a matrix-model + limitation or a model-family (pretrained-on-MRI) limitation? + +D. Network-space mechanics -- static (no registration) check of what the + network actually sees, emulating ``network_space`` + ``normalize`` from + the image code (isotropic 1-mm LIA grid, extent 256, centred on each + volume's own centre, fill 0, min-max into [0,1]): + - head truncation: fraction of head-mask voxels (the T1c head mapped into + CT space via the reference transform) outside the 256 mm box, for the CT + and T1c network spaces, + - intensity band: the min-max-normalized values over the head region for + the T1c, the [0,80]-clipped CT (the benchmark default) and the raw CT. + +Usage (repo root, conda env longitudinal-313): + python scripts/investigate_synthmorph.py --workers 3 + python scripts/investigate_synthmorph.py --parts D # static only + python scripts/investigate_synthmorph.py --keys m6_2BBN5XA4_20250411 +""" +import argparse +import importlib.util +import json +import math +import os +import shutil +import subprocess +import sys +import tempfile +import time +from concurrent.futures import ProcessPoolExecutor, as_completed + +ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, ROOT) + +import numpy as np +import SimpleITK as sitk +from scipy.ndimage import center_of_mass + +from src.common import d, read_nii_img, head_mask_from_image +from src import registration_pipelines as RP + + +def _load_module(name, rel): + p = os.path.join(ROOT, rel) + spec = importlib.util.spec_from_file_location(name, p) + m = importlib.util.module_from_spec(spec) + spec.loader.exec_module(m) + return m + + +BR = _load_module("bench_registration", "scripts/bench_registration.py") +M6 = RP._load_m6() + +# the worst 3 bench cases (rigid MI): ~30 deg tilt x2 + a mid-tier case +JOINT_KEYS = ("m6_2BBN5XA4_20250411", "m6_24V62QQM_20240820", + "m6_24V62QQM_20230223") + + +# --------------------------------------------------------------------------- # +# transform decomposition +# --------------------------------------------------------------------------- + +def _rot_angle(R): + c = (np.trace(R) - 1.0) / 2.0 + return math.degrees(math.acos(float(np.clip(c, -1.0, 1.0)))) + + +def _rot_axis(R, theta_deg): + if theta_deg <= 0.5: + return np.zeros(3) + axis = np.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]]) + n = np.linalg.norm(axis) + return axis / n if n > 1e-9 else np.zeros(3) + + +def _to_affine(T): + """Recover (R, t) of any SimpleITK transform by point mapping (world RAS).""" + o = np.array(T.TransformPoint((0.0, 0.0, 0.0)), dtype=float) + R = np.empty((3, 3), dtype=float) + for i in range(3): + p = [0.0, 0.0, 0.0] + p[i] = 1.0 + R[:, i] = np.array(T.TransformPoint(tuple(p)), dtype=float) - o + return R, o + + +def _decompose_affine(T): + """(R, t, theta_deg, axis, rigidity) of a SimpleITK transform in world RAS.""" + R, t = _to_affine(T) + theta = _rot_angle(R) + rig = {"det": round(float(np.linalg.det(R)), 6), + "max_devi_from_I": round(float(np.linalg.norm(R.T @ R - np.eye(3))), 6)} + return R, t, theta, _rot_axis(R, theta), rig + + +def _angle_between_rotations(A, B): + """angle of A^T B in degrees (rotation B needs on top of A).""" + return _rot_angle(A.T @ B) + + +def _polar_rotation(M): + """rotation part of a general 3x3 linear map (SVD polar).""" + U, _s, Vt = np.linalg.svd(M) + R_p = U @ Vt + if np.linalg.det(R_p) < 0: + U[:, -1] *= -1 + R_p = U @ Vt + return R_p + + +def head_centroid_world(t1c): + ha = sitk.GetArrayFromImage(t1c).astype("float32") + idx = center_of_mass(head_mask_from_image(ha)) + sp = np.array(t1c.GetSpacing()) + o = np.array(t1c.GetOrigin()) + D = np.array(t1c.GetDirection()).reshape(3, 3) + return o + D @ (np.array(idx[::-1]) * sp) + + +def score_identity(ct, t1c): + """The "do nothing" baseline: CT resampled onto the T1c grid with the + identity transform (i.e. no registration at all).""" + ha = sitk.GetArrayFromImage(t1c).astype("float32") + wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, sitk.Transform(), + sitk.sitkLinear)) + gap = M6._skull_ring_gap(wc, ha) + mi = BR._mi_masked(wc, ha, head_mask_from_image(ha)) + return {"mi": mi, "gap": gap} + + +# --------------------------------------------------------------------------- # +# part A/B/D worker +# --------------------------------------------------------------------------- + +def _case_task(case, args): + out = {"key": case["key"], "subject": case["subject"], "date": case["date"]} + ct = read_nii_img(case["ct"]) + t1c = read_nii_img(case["img"]) + tv = read_nii_img(case["label"]) + if not BR.ct_complete(ct): + out["note"] = "ct_incomplete (skipped, same as the benchmark)" + return out + + gpu = (args.gpu_base + abs(hash(out["key"])) % max(1, args.n_gpus)) % max(1, args.n_gpus) + parts = set(args.parts) + + # geometry facts + for nm, img in (("ct", ct), ("t1c", t1c)): + sz = img.GetSize() + sp = img.GetSpacing() + c = np.array(img.TransformIndexToPhysicalPoint( + ((sz[0] - 1) // 2, (sz[1] - 1) // 2, (sz[2] - 1) // 2))) + out[nm] = {"shape_xyz": list(sz), + "spacing_xyz": [round(float(s), 3) for s in sp], + "fov_mm_xyz": [round(float(s0 * s1), 1) for s0, s1 in zip(sz, sp)], + "center_world": [round(float(x), 1) for x in c]} + c_t1c = head_centroid_world(t1c) + out["t1c_head_centroid_world"] = [round(float(x), 1) for x in c_t1c] + + R_ref = t_ref = None + if parts & {"A", "B", "D"}: + # reference: simpleitk benchmark mode (deterministic at fixed threads) + T_ref, meta_ref = RP.register_simpleitk(ct, t1c, threads=16, + mode="benchmark") + sref = BR.evaluate(ct, t1c, tv, T_ref) + R_ref, t_ref, th_ref, ax_ref, _ = _decompose_affine(T_ref) + out["reference"] = { + "tool": meta_ref["toolkit"], + "R": R_ref.tolist(), "t": t_ref.tolist(), + "mi": sref["mi_ct_mr"], "impact_loss": sref["impact_loss"], + "gap": sref["skull_ring_gap_mm"], + "rotation_deg": round(th_ref, 2), + "axis": [round(float(a), 3) for a in ax_ref], + "translation": [round(float(x), 1) for x in t_ref], + "head_centre_ct": [round(float(x), 1) for x in (R_ref @ c_t1c + t_ref)], + } + out["identity"] = score_identity(ct, t1c) + + if "A" in parts: + t0 = time.time() + T_sm, meta_sm = RP.register_synthmorph(ct, t1c, model="rigid", + gpu=int(gpu), extent=256) + ssm = BR.evaluate(ct, t1c, tv, T_sm) + R_sm, t_sm, th_sm, ax_sm, rig = _decompose_affine(T_sm) + p_ref = R_ref @ c_t1c + t_ref + p_sm = R_sm @ c_t1c + t_sm + out["rigid"] = { + "seconds": round(time.time() - t0, 1), + "auto_inverted": meta_sm["auto_inverted"], + "clip_ct": meta_sm["clip_ct"], + "mi": ssm["mi_ct_mr"], "impact_loss": ssm["impact_loss"], + "gap": ssm["skull_ring_gap_mm"], + "rotation_deg": round(th_sm, 2), + "axis": [round(float(a), 3) for a in ax_sm], + "translation": [round(float(x), 1) for x in t_sm], + "matrix_rigidity": rig, + "rot_err_vs_ref_deg": round(_angle_between_rotations(R_ref, R_sm), 2), + "head_centre_err_mm": round(float(np.linalg.norm(p_sm - p_ref)), 1), + } + + if "B" in parts: + t0 = time.time() + T_af, meta_af = RP.register_synthmorph(ct, t1c, model="affine", + gpu=int(gpu), extent=256) + saf = BR.evaluate(ct, t1c, tv, T_af) + M_af, t_af = _to_affine(T_af) + svals = np.linalg.svd(M_af, compute_uv=False) + lin = {"scale_mean": round(float(np.mean(svals)), 4), + "scale_iso_dev": round(float(np.std(svals)) / max(float(np.mean(svals)), 1e-12), 4), + "rotation_deg": round(_rot_angle(_polar_rotation(M_af)), 2), + "det": round(float(np.linalg.det(M_af)), 4)} + row = {"seconds": round(time.time() - t0, 1), + "auto_inverted": meta_af["auto_inverted"], + "mi": saf["mi_ct_mr"], "impact_loss": saf["impact_loss"], + "gap": saf["skull_ring_gap_mm"], "linear": lin} + if R_ref is not None: + row["rot_err_vs_ref_deg"] = round( + _angle_between_rotations(R_ref, _polar_rotation(M_af)), 2) + row["head_centre_err_mm"] = round( + float(np.linalg.norm(M_af @ c_t1c + t_af - (R_ref @ c_t1c + t_ref))), 1) + out["affine"] = row + + if "D" in parts: + out["network_space"] = _network_space_analysis(ct, t1c, out, + R_ref=R_ref, t_ref=t_ref) + + return out + + +# --------------------------------------------------------------------------- # +# part D: what the network sees +# --------------------------------------------------------------------------- + +def _network_grid(img, extent=256): + """The matrix-model network space of ``img``: isotropic 1-mm LIA grid, + ``extent``^3, centred on the volume's own centre (mirrors + ``network_space`` in the container's registration.py).""" + sz = img.GetSize() + c = np.array(img.TransformIndexToPhysicalPoint( + ((sz[0] - 1) // 2, (sz[1] - 1) // 2, (sz[2] - 1) // 2))) + # LIA axes in the RAS world: x_L = -x_R, y_I = -y_A, z_A = +z_R + origin = (c[0] + extent / 2.0 - 0.5, c[1] + extent / 2.0 - 0.5, + c[2] - extent / 2.0 + 0.5) + net = sitk.Image((extent, extent, extent), sitk.sitkFloat32) + net.SetSpacing((1.0, 1.0, 1.0)) + net.SetOrigin(tuple(float(x) for x in origin)) + net.SetDirection((-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)) + return net, c + + +def _to_network(img, net): + a = sitk.GetArrayFromImage( + sitk.Resample(img, net, sitk.Transform(), sitk.sitkLinear, 0.0)) + lo, hi = float(a.min()), float(a.max()) + if hi - lo < 1e-9: + return np.zeros_like(a), lo, hi + return (a - lo) / (hi - lo), lo, hi + + +def _network_space_analysis(ct, t1c, out, extent=256, R_ref=None, t_ref=None): + """Extent truncation + min-max intensity band (what the model sees).""" + res = {"extent": extent} + ha = sitk.GetArrayFromImage(t1c).astype("float32") + head_t1c = sitk.GetImageFromArray(head_mask_from_image(ha).astype("float32")) + head_t1c.CopyInformation(t1c) + + if R_ref is not None: + T_ref = sitk.AffineTransform(3) + T_ref.SetMatrix(list(np.array(R_ref).ravel())) + T_ref.SetTranslation(list(np.array(t_ref))) + # T_ref maps T1c -> CT, so rasterising the T1c head onto the CT grid + # needs the inverse as the destination -> source map + head_ct = sitk.Resample(head_t1c, ct, T_ref.GetInverse(), + sitk.sitkNearestNeighbor, 0.0) + else: + head_ct = head_t1c + + # truncation: head voxels outside each volume's OWN extent box + for nm, img, head in (("ct", ct, head_ct), ("t1c", t1c, head_t1c)): + _net, c = _network_grid(img, extent) + idx = np.argwhere(sitk.GetArrayFromImage(head) > 0.5) + if idx.size == 0: + res[f"{nm}_head_trunc_frac"] = None + continue + sp = np.array(img.GetSpacing()) + o = np.array(img.GetOrigin()) + D = np.array(img.GetDirection()).reshape(3, 3) + world = o[None, :] + (idx[:, ::-1] * sp) @ D.T # idx is (z,y,x) + half = extent / 2.0 + inside = (np.abs(world - c[None, :]) <= half).all(axis=1) + head_c = world.mean(axis=0) + res[f"{nm}_head_voxels"] = int(idx.shape[0]) + res[f"{nm}_head_trunc_frac"] = round(float(1.0 - inside.mean()), 4) + res[f"{nm}_head_off_volume_centre_mm"] = \ + [round(float(x), 1) for x in (head_c - c)] + + # min-max normalized band over the head region (emulating normalize=True) + bands = {} + for nm, img, head in (("t1c", t1c, head_t1c), + ("ct_clip", RP._clip_intensity(ct, 0.0, 80.0), head_ct), + ("ct_raw", ct, head_ct)): + net, _c = _network_grid(img, extent) + a, lo, hi = _to_network(img, net) + hnet = sitk.GetArrayFromImage( + sitk.Resample(head, net, sitk.Transform(), + sitk.sitkNearestNeighbor, 0.0)) > 0.5 + v = a[hnet] + bands[nm] = {"src_min": round(lo, 1), "src_max": round(hi, 1), + "head_nvox_net": int(v.size), + "p5": round(float(np.percentile(v, 5)), 4) if v.size else None, + "p50": round(float(np.percentile(v, 50)), 4) if v.size else None, + "p95": round(float(np.percentile(v, 95)), 4) if v.size else None, + "band_width_p95_p5": (round(float(np.percentile(v, 95) + - np.percentile(v, 5)), 4) + if v.size else None)} + res["normalized_band"] = bands + return res + + +# --------------------------------------------------------------------------- # +# part C: joint (affine + deformable) upper bound +# --------------------------------------------------------------------------- + +def _joint_task(case, args): + out = {"key": case["key"], "subject": case["subject"], "date": case["date"], + "model": "joint (affine + deformable)"} + ct = read_nii_img(case["ct"]) + t1c = read_nii_img(case["img"]) + if not BR.ct_complete(ct): + out["note"] = "ct_incomplete (skipped)" + return out + gpu = (args.gpu_base + abs(hash(out["key"])) % max(1, args.n_gpus)) % max(1, args.n_gpus) + tmp = tempfile.mkdtemp(prefix="sm_joint_") + try: + ct_clip = RP._clip_intensity(ct, 0.0, 80.0) + sitk.WriteImage(sitk.Cast(ct_clip, sitk.sitkFloat32), + os.path.join(tmp, "moving.nii.gz"), True) + sitk.WriteImage(sitk.Cast(t1c, sitk.sitkFloat32), + os.path.join(tmp, "fixed.nii.gz"), True) + # CLI order is (moving, fixed); -o saves moving registered to fixed + cmd = ["docker", "run", "--rm", + "--user", f"{os.getuid()}:{os.getgid()}", + "-e", "TF_FORCE_GPU_ALLOW_GROWTH=true", + "--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}", + "-v", f"{tmp}:/data", RP.SYNTHMORPH_IMAGE, + "register", "-m", "joint", "-e", "256", "-g", + "-o", "/data/warped.nii.gz", + "/data/moving.nii.gz", "/data/fixed.nii.gz"] + t0 = time.time() + proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200) + dt = time.time() - t0 + if proc.returncode != 0: + out["error"] = (f"docker failed rc={proc.returncode}:\n" + + "\n".join((proc.stdout + proc.stderr).splitlines()[-15:])) + return out + warped = sitk.GetArrayFromImage( + sitk.ReadImage(os.path.join(tmp, "warped.nii.gz"))).astype("float32") + ha = sitk.GetArrayFromImage(t1c).astype("float32") + gap = M6._skull_ring_gap(warped, ha) + mi = BR._mi_masked(warped, ha, head_mask_from_image(ha)) + il = BR._impact_loss_value(warped, ha) + out.update({"seconds": round(dt, 1), "mi": mi, "impact_loss": il, + "gap": gap, "gpu": int(gpu)}) + finally: + if not args.keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + return out + + +# --------------------------------------------------------------------------- # +# part F: re-centred CT probe (volume-centre anchoring hypothesis) +# --------------------------------------------------------------------------- + +def _recenter_task(case, args, ref_map): + """H5 probe: the matrix model's network space is centred on each volume's + OWN centre. The T1c is a tight head crop (centre ~ head centre), but the + head+neck CT's volume centre sits 15-80 mm below the head. If the model + under-corrects that offset, cropping the CT to the head box (re-centring + the volume on the head) should remove the translation miss. The crop is + world-preserving, so the returned transform scores against the FULL CT + exactly like the benchmark.""" + out = {"key": case["key"], "subject": case["subject"], "date": case["date"]} + ct = read_nii_img(case["ct"]) + t1c = read_nii_img(case["img"]) + tv = read_nii_img(case["label"]) + if not BR.ct_complete(ct): + out["note"] = "ct_incomplete (skipped)" + return out + ref = ref_map.get(case["key"]) + if not ref: + out["error"] = "no reference transform in partABD.json" + return out + R = np.array(ref["R"]) + t = np.array(ref["t"]) + T_ref = sitk.AffineTransform(3) + T_ref.SetMatrix(list(R.ravel())) + T_ref.SetTranslation(list(t)) + + gpu = (args.gpu_base + abs(hash(out["key"])) % max(1, args.n_gpus)) % max(1, args.n_gpus) + + # head in CT space (T_ref maps T1c -> CT, so rasterise with its inverse) + ha = sitk.GetArrayFromImage(t1c).astype("float32") + head_t1c = sitk.GetImageFromArray(head_mask_from_image(ha).astype("float32")) + head_t1c.CopyInformation(t1c) + head_ct = sitk.Resample(head_t1c, ct, T_ref.GetInverse(), + sitk.sitkNearestNeighbor, 0.0) + m = sitk.GetArrayFromImage(head_ct) > 0.5 + idx = np.argwhere(m) # (z, y, x) + lo, hi = idx.min(0), idx.max(0) + cen = (lo + hi + 1) / 2.0 + half = (hi - lo) / 2.0 + 15.0 / np.array(ct.GetSpacing())[::-1] # (z,y,x) + sz = np.array(ct.GetSize())[::-1] # GetSize is (x,y,z); work in (z,y,x) + b0 = np.clip(np.floor(cen - half).astype(int), 0, sz - 1) + b1 = np.clip(np.ceil(cen + half).astype(int), 0, sz) + + def _volcentre(img): + s = img.GetSize() + return np.array(img.TransformIndexToPhysicalPoint( + ((s[0] - 1) // 2, (s[1] - 1) // 2, (s[2] - 1) // 2))) + + hc = T_ref.TransformPoint(tuple(float(x) for x in head_centroid_world(t1c))) + out["before"] = {"head_off_volume_centre_mm": [ + round(float(x), 1) for x in (np.array(hc) - _volcentre(ct))]} + + # SimpleITK image slicing is in (x, y, z) order; b0/b1 are (z, y, x) + ct_crop = ct[int(b0[2]):int(b1[2]), int(b0[1]):int(b1[1]), + int(b0[0]):int(b1[0])] + out["crop"] = {"shape": list(ct_crop.GetSize()), + "head_off_volume_centre_mm": [ + round(float(x), 1) + for x in (np.array(hc) - _volcentre(ct_crop))]} + + t0 = time.time() + T_sm, meta_sm = RP.register_synthmorph(ct_crop, t1c, model="rigid", + gpu=int(gpu), extent=256) + ssm = BR.evaluate(ct, t1c, tv, T_sm) # full CT, as in the benchmark + R_sm, t_sm = _to_affine(T_sm) + c_t1c = head_centroid_world(t1c) + out["rigid_recentered"] = { + "seconds": round(time.time() - t0, 1), + "auto_inverted": meta_sm["auto_inverted"], + "mi": ssm["mi_ct_mr"], "impact_loss": ssm["impact_loss"], + "gap": ssm["skull_ring_gap_mm"], + "rot_err_vs_ref_deg": round(_angle_between_rotations(R, R_sm), 2), + "head_centre_err_mm": round( + float(np.linalg.norm((R_sm @ c_t1c + t_sm) - (R @ c_t1c + t))), 1), + } + return out + + +# --------------------------------------------------------------------------- # +# main +# --------------------------------------------------------------------------- + +def parse_args(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--parts", default="ABCD", + help="which parts to run (letters A/B/C/D/F; F = re-centred " + "CT probe, needs the partABD.json reference)") + ap.add_argument("--n-cases", type=int, default=10) + ap.add_argument("--keys", default=None, help="comma-separated case keys") + ap.add_argument("--workers", type=int, default=3) + ap.add_argument("--gpu-base", type=int, default=0) + ap.add_argument("--n-gpus", type=int, default=3) + ap.add_argument("--out-dir", default=os.path.join(ROOT, "results", + "synthmorph_investigation")) + ap.add_argument("--keep-scratch", action="store_true") + return ap.parse_args() + + +def main(): + args = parse_args() + os.makedirs(args.out_dir, exist_ok=True) + labeled, _ = M6.discover() + pool = list(labeled) + if args.keys: + bykey = {r["key"]: r for r in pool} + missing = [k for k in args.keys.split(",") if k not in bykey] + if missing: + raise SystemExit(f"keys not found: {missing}") + cases = [bykey[k.strip()] for k in args.keys.split(",")] + else: + cases = pool[:args.n_cases] + print(f"[investigate] {len(cases)} cases, parts={args.parts}") + + rows, joint_rows, recenter_rows = [], [], [] + do_cases = bool(set(args.parts) & {"A", "B", "D"}) + do_joint = "C" in args.parts + do_recenter = "F" in args.parts + + ref_map = {} + if do_recenter: + fp0 = os.path.join(args.out_dir, "partABD.json") + if os.path.exists(fp0): + with open(fp0) as f: + ref_map = {r["key"]: r["reference"] + for r in json.load(f) if "reference" in r} + + with ProcessPoolExecutor(max_workers=args.workers) as ex: + futs = {} + if do_cases: + for c in cases: + futs[ex.submit(_case_task, c, args)] = ("case", c) + if do_joint: + for c in cases: + if c["key"] in JOINT_KEYS: + futs[ex.submit(_joint_task, c, args)] = ("joint", c) + if do_recenter: + for c in cases: + futs[ex.submit(_recenter_task, c, args, ref_map)] = ("recenter", c) + for fu in as_completed(futs): + kind, c = futs[fu] + try: + r = fu.result() + if kind == "joint": + joint_rows.append(r) + print(f" [joint] {c['key']}: mi={r.get('mi')} " + f"il={r.get('impact_loss')} {r.get('seconds')}s " + f"err={(r.get('error') or '')[:100]}", flush=True) + elif kind == "recenter": + recenter_rows.append(r) + rr = r.get("rigid_recentered", {}) + print(f" [recenter] {c['key']}: mi={rr.get('mi')} " + f"ctr_err={rr.get('head_centre_err_mm')} " + f"rot_err={rr.get('rot_err_vs_ref_deg')} " + f"off_before={r.get('before', {}).get('head_off_volume_centre_mm')} " + f"off_after={r.get('crop', {}).get('head_off_volume_centre_mm')} " + f"note={(r.get('note') or r.get('error') or '')[:60]}", flush=True) + else: + rows.append(r) + print(f" [case] {c['key']}: ref.mi={r.get('reference', {}).get('mi')} " + f"rigid.mi={r.get('rigid', {}).get('mi')} " + f"affine.mi={r.get('affine', {}).get('mi')} " + f"note={r.get('note', r.get('error', ''))[:60]}", flush=True) + except Exception as e: # noqa + print(f" {c['key']}: ERROR {e!r}", flush=True) + (joint_rows if kind == "joint" else + recenter_rows if kind == "recenter" else rows).append( + {"key": c["key"], "error": repr(e)}) + rows.sort(key=lambda x: x["key"]) + joint_rows.sort(key=lambda x: x["key"]) + recenter_rows.sort(key=lambda x: x["key"]) + + if do_cases: + fp = os.path.join(args.out_dir, "partABD.json") + merged = {r["key"]: r for r in rows} + if os.path.exists(fp): + with open(fp) as f: + for r in json.load(f): + # merge missing sections (a partial run must not wipe + # data from a fuller earlier run) + if r["key"] in merged: + for k, v in r.items(): + merged[r["key"]].setdefault(k, v) + else: + merged[r["key"]] = r + with open(fp, "w") as f: + json.dump(sorted(merged.values(), key=lambda x: x["key"]), + f, indent=1) + print("wrote", fp) + if do_joint: + fp = os.path.join(args.out_dir, "partC_joint.json") + with open(fp, "w") as f: + json.dump(joint_rows, f, indent=1) + print("wrote", fp) + if do_recenter: + fp = os.path.join(args.out_dir, "partF_recenter.json") + with open(fp, "w") as f: + json.dump(recenter_rows, f, indent=1) + print("wrote", fp) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/scripts/investigate_synthmorph_mm.py b/scripts/investigate_synthmorph_mm.py new file mode 100644 index 0000000..5fc38f2 --- /dev/null +++ b/scripts/investigate_synthmorph_mm.py @@ -0,0 +1,101 @@ +"""MR->MR control for the synthmorph investigation. + +The benchmark pairs are CT->T1c (out-of-domain for the MRI-trained +SynthMorph model). To separate the cross-modality effect from the other +suspects (extent truncation of the CT FOV, min-max normalization), this +registers two T1c scans of the SAME subject (m6_24V62QQM, 2023-02-23 and +2023-10-17 -- in-domain, MRI->MRI, same tightly-cropped geometry as the +benchmark's fixed image) and scores the result with the benchmark scoring +(MI over the head, ImpactLoss, skull-ring gap) plus the transform +decomposition against a SimpleITK reference. + +Usage: + python scripts/investigate_synthmorph_mm.py [--out results/synthmorph_investigation] +""" +import argparse +import json +import os +import sys +import time + +ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, ROOT) + +import numpy as np +import SimpleITK as sitk + +from src.common import read_nii_img, head_mask_from_image +from src import registration_pipelines as RP +from scripts.investigate_synthmorph import ( + _to_affine, _angle_between_rotations, head_centroid_world, BR, M6) + +KEY_A = "m6_24V62QQM_20230223" +KEY_B = "m6_24V62QQM_20231017" + + +def score_pair(img_moving, img_fixed, T): + """Benchmark scoring of the warped moving (in the fixed grid) ~ fixed.""" + ha = sitk.GetArrayFromImage(img_fixed).astype("float32") + wc = sitk.GetArrayFromImage(sitk.Resample(img_moving, img_fixed, T, + sitk.sitkLinear)).astype("float32") + gap = M6._skull_ring_gap(wc, ha) + mi = BR._mi_masked(wc, ha, head_mask_from_image(ha)) + il = BR._impact_loss_value(wc, ha) + return {"mi": mi, "impact_loss": il, "gap": gap} + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--out", default=os.path.join(ROOT, "results", + "synthmorph_investigation")) + ap.add_argument("--gpu", type=int, default=0) + args = ap.parse_args() + os.makedirs(args.out, exist_ok=True) + + labeled, _ = M6.discover() + bykey = {r["key"]: r for r in labeled} + a, b = bykey[KEY_A], bykey[KEY_B] + img_a = read_nii_img(a["img"]) + img_b = read_nii_img(b["img"]) + out = {"pair": [KEY_A, KEY_B], "note": ("same-subject T1c->T1c, 8 months " + "apart; in-domain control for the " + "CT->T1c benchmark pairs")} + + import math + c_fix = head_centroid_world(img_b) + T_ref, meta_ref = RP.register_simpleitk(img_a, img_b, threads=16, + mode="benchmark") + sref = score_pair(img_a, img_b, T_ref) + R_ref, t_ref = _to_affine(T_ref) + th_ref = math.degrees(float(np.clip((np.trace(R_ref) - 1) / 2, -1, 1))) + out["reference"] = { + "tool": meta_ref["toolkit"], "R": R_ref.tolist(), "t": t_ref.tolist(), + "mi": sref["mi"], "impact_loss": sref["impact_loss"], "gap": sref["gap"], + "rotation_deg": round(th_ref, 2), + "head_centre_fix": [round(float(x), 1) + for x in (R_ref @ c_fix + t_ref)], + } + + t0 = time.time() + T_sm, meta_sm = RP.register_synthmorph(img_a, img_b, model="rigid", + gpu=args.gpu, extent=256) + ssm = score_pair(img_a, img_b, T_sm) + R_sm, t_sm = _to_affine(T_sm) + out["synthmorph_rigid"] = { + "seconds": round(time.time() - t0, 1), + "auto_inverted": meta_sm["auto_inverted"], + "mi": ssm["mi"], "impact_loss": ssm["impact_loss"], "gap": ssm["gap"], + "rot_err_vs_ref_deg": round(_angle_between_rotations(R_ref, R_sm), 2), + "head_centre_err_mm": round( + float(np.linalg.norm((R_sm @ c_fix + t_sm) - (R_ref @ c_fix + t_ref))), 1), + } + # no-clip sanity: MR inputs are never clipped by the adapter + fp = os.path.join(args.out, "partE_mr_mr.json") + with open(fp, "w") as f: + json.dump(out, f, indent=1) + print(json.dumps(out, indent=1)) + print("wrote", fp) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index 438159f..76a8747 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -32,7 +32,7 @@ Shared rules for every pipeline: 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`. - **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_RIGID.md), ``niftyreg`` and ``synthmorph`` clip the CT (the moving image) to [0, 80] HU before registration (``clip_ct`` parameter, default on, MR inputs never clipped); the Mattes-MI pipelines (antspyx / elastix / simpleitk) register the CT @@ -262,10 +262,7 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=2000): return tr, meta m6 = _load_m6() - try: - sitk.CommonProperties.SetGlobalDefaultNumberOfThreads(threads) - except AttributeError: - pass + _sitk_set_threads(threads) ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32) # full volumes, no mask / crop / pad; the T1c is the fixed image fixed, moving = fixed_moving(ct, t1c) @@ -1221,6 +1218,851 @@ REGISTERS = { ALL_PIPELINES = list(REGISTERS) +# --------------------------------------------------------------------------- # +# Deformable variants (see REGISTRATION_BENCHMARK_DEFORMABLE.md) +# --------------------------------------------------------------------------- # +# +# Same unified contract as the rigid pipelines: the returned transform ``T`` +# satisfies ``sitk.Resample(ct, t1c, T, ...)`` = CT on the T1c grid, i.e. it +# maps **T1c -> CT** coordinates. Every deformable pipeline runs its toolkit's +# native deformable workflow (linear init stage + deformable stage where that +# is the toolkit's canonical setup) and returns the result rasterised as a +# DisplacementFieldTransform on the T1c grid -- the one transform type all the +# deformable outputs share (elastix BSplines, ANTs warp+affine, NiftyReg +# affine+spline, SynthMorph fields). + +def _sitk_set_threads(n): + """Best-effort global SimpleITK/ITK work-unit count (the attribute name + varies between builds: CommonProperties vs the procedural API).""" + try: + sitk.CommonProperties.SetGlobalDefaultNumberOfThreads(n) + except AttributeError: + try: + sitk.ProcessObject_SetGlobalDefaultNumberOfThreads(n) + except AttributeError: + pass + + +def _field_transform(arr, ref_img): + """SimpleITK DisplacementFieldTransform from a displacement array + (world-mm) on the grid of ``ref_img``. Accepts (z,y,x,3) or (3,z,y,x).""" + a = np.asarray(arr, dtype="float64") + if a.ndim == 4 and a.shape[0] == 3 and a.shape[-1] != 3: + a = np.ascontiguousarray(np.transpose(a, (1, 2, 3, 0))) + Fv = sitk.GetImageFromArray(a, isVector=True) + Fv.SetSpacing(ref_img.GetSpacing()) + Fv.SetOrigin(ref_img.GetOrigin()) + Fv.SetDirection(ref_img.GetDirection()) + return sitk.DisplacementFieldTransform(Fv) + + +def _flip_naive_world_components(arr): + """Convert displacement deltas from the naive NIfTI world to true RAS. + + NiftyReg and surfa (the FreeSurfer NIfTI stack) decode a RAS+ NIfTI qform + with the plain quaternion formula, so they work in an x/y-mirrored world. + A delta (dx, dy, dz) they store is therefore (-dx, -dy, dz) in the true + RAS world; the field grid (file indices) is unaffected. Our inputs are + always written RAS+ (read_nii_img DICOMOrients), so the flip applies. + Verified against each tool's own resampling on the M6 cases (NCC ~ 1.0). + """ + a = np.array(arr, dtype="float64", copy=True) + a[..., 0] = -a[..., 0] + a[..., 1] = -a[..., 1] + return a + + +def _rasterize_to_tdf(T, ref_img): + """Rasterize any SimpleITK transform (affine / bspline / composite) to a + DisplacementFieldTransform on the grid of ``ref_img``.""" + f = sitk.TransformToDisplacementField( + T, sitk.sitkVectorFloat64, ref_img.GetSize(), ref_img.GetOrigin(), + ref_img.GetSpacing(), ref_img.GetDirection()) + return sitk.DisplacementFieldTransform(f) + + +def _itk_rasterize_to_tdf(itk_tf, ref_img): + """Rasterize an ITK-python transform (elastix combination transform) to a + SimpleITK DisplacementFieldTransform on the grid of ``ref_img``.""" + import itk + from itk.itkImagePython import itkImageVF33 + flt = itk.TransformToDisplacementFieldFilter[itkImageVF33, itk.D].New() + flt.SetTransform(itk_tf) + flt.SetSize(list(ref_img.GetSize())) + flt.SetOutputOrigin([float(x) for x in ref_img.GetOrigin()]) + flt.SetOutputSpacing([float(x) for x in ref_img.GetSpacing()]) + flt.Update() + arr = itk.GetArrayFromImage(flt.GetOutput()) # (z, y, x, 3), world mm + return _field_transform(arr, ref_img) + + +def _bspline_txt_to_tdf(txt_path, ref_img): + """Evaluate an elastix BSpline/RecursiveBSpline transform .itk.txt on the + grid of ``ref_img`` and return a DisplacementFieldTransform. + + SimpleITK's transform reader does not parse the elastix BSpline file + format, so the cubic B-spline (control-point displacements on a regular + grid) is evaluated directly with scipy's order-3 B-spline interpolation + (the same kernel as ITK's BSpline evaluation). + + File format (ITK V1.0): + Transform: BSplineTransform_double_3_3 + Parameters: + FixedParameters: + Parameter layout (verified by probing): component-major blocks + [all x-offsets, all y-offsets, all z-offsets], x-fastest within a block + (slot = ix + mx*(iy + my*iz)).""" + text = open(txt_path).read() + if "BSplineTransform" not in text: + raise RuntimeError( + f"_bspline_txt_to_tdf: {txt_path} is not a BSpline transform " + f"({text.splitlines()[2] if len(text.splitlines()) > 2 else '?'})") + body = text.split("Parameters:", 1)[1].split("FixedParameters:", 1)[0] + vals = [float(x) for x in body.split()] + fp = [float(x) for x in text.split("FixedParameters:", 1)[1].split()] + mx, my, mz = int(fp[0]), int(fp[1]), int(fp[2]) + g_origin = np.array(fp[3:6]); g_sp = np.array(fp[6:9]) + g_dir = np.array(fp[9:18]).reshape(3, 3) + n_pts = mx * my * mz + if len(vals) != 3 * n_pts: + raise RuntimeError( + f"_bspline_txt_to_tdf: expected {3 * n_pts} parameters in {txt_path}, " + f"got {len(vals)}") + blocks = [np.reshape(vals[i * n_pts:(i + 1) * n_pts], (mx, my, mz), order="F") + for i in range(3)] # each (mx, my, mz), x fastest + o = np.array(ref_img.GetOrigin()) + d = np.array(ref_img.GetDirection()).reshape(3, 3) + s = np.array(ref_img.GetSpacing()) + zz, yy, xx = np.mgrid[0:ref_img.GetSize()[2], 0:ref_img.GetSize()[1], + 0:ref_img.GetSize()[0]] + p = np.stack([o[0] + xx * s[0], o[1] + yy * s[1], o[2] + zz * s[2]], + axis=-1) # (z, y, x, 3) world RAS + u = (p - g_origin) @ d.T @ g_dir.T / g_sp # (z, y, x, 3) grid coords + from scipy.ndimage import map_coordinates + disp = np.stack( + [map_coordinates(b, [u[..., 2], u[..., 1], u[..., 0]], order=3, + mode="nearest") for b in blocks], axis=-1) + return _field_transform(disp, ref_img) + + +def register_antspyx_deformable(ct, t1c): + """ANTs SyN -- antspyx's native deformable (type_of_transform='SyN'): + an Affine stage (Mattes MI, convergence 2100x1200x1200x0, shrink 4x2x2x1) + followed by a symmetric-Normalization deformable stage (Mattes MI, + SyN[0.2, 3, 0], convergence 40x20x0, shrink 4x2x1), all on the full + volumes with the native center-of-mass init. No CT clip (the rigid + A/B showed clipping hurts the Mattes-MI pipelines). The two stage + outputs (001Warp.nii.gz + 000GenericAffine.mat) are composed -- warp + first, then affine -- and verified against ANTs' own warpedmovout, + then rasterised to the T1c-grid displacement field of the contract.""" + import ants + fx, mv = fixed_moving(ct, t1c) + tmp = tempfile.mkdtemp(prefix="ants_syn_") + dt = 0.0 + try: + import time + t0 = time.time() + r = ants.registration( + fixed=ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)), + moving=ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)), + type_of_transform="SyN", outprefix=tmp + "/") + dt = time.time() - t0 + aff = warp = None + for f in r["fwdtransforms"]: + base = os.path.basename(f) + if base.endswith(".mat"): + aff = sitk.ReadTransform(f) + elif base.endswith("Warp.nii.gz") and "Inverse" not in base: + w = sitk.ReadImage(f) + arr = sitk.GetArrayFromImage(w) # (z,y,x,3) vector, world mm + if arr.ndim == 4 and arr.shape[0] == 3: + arr = np.transpose(arr, (1, 2, 3, 0)) + warp = _field_transform(arr, fx) + if aff is None or warp is None: + raise RuntimeError( + f"antspyx SyN: expected an affine .mat + a Warp.nii.gz stage " + f"file, got {[os.path.basename(f) for f in r['fwdtransforms']]}") + # composite: the SyN stage warp is applied first, then the affine + # (verified to reproduce ANTs' own warpedmovout to <0.001 MI) + T = _rasterize_to_tdf(sitk.CompositeTransform([aff, warp]), fx) + finally: + shutil.rmtree(tmp, ignore_errors=True) + ha = sitk.GetArrayFromImage(t1c) + gap = _skull_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": "antspyx", + "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", + "metric": "ANTs Mattes MI (affine + SyN stages, full volumes)", + "transform": ("Deformable (native SyN: Affine[0.25] 2100x1200x1200x0 " + "+ SyN[0.2,3,0] 40x20x0, auto COM init)"), + "convention": "T1c->CT (verified vs ANTs warpedmovout)", "auto_inverted": False, + "fixed": _modality_name(fx), "moving": _modality_name(mv), + "seconds_reg": round(dt, 1), "skull_ring_gap_mm": gap, + } + return T, meta + + +def register_elastix_deformable(ct, t1c, workdir=None, threads=8): + """ITK-Elastix two-stage deformable (the toolkit's canonical deformable + workflow): stage 1 is the rigid benchmark's EulerTransform map (Mattes MI, + ASGD, 2000 iters x 3 resolutions, CenterOfGravity auto init); stage 2 is a + BSplineTransform (16 mm grid, spline order 3) initialized from stage 1 via + -t0 with the same Mattes-MI budget. Runs in-process through the + ITK-Elastix python API; the combination transform is rasterised to the + T1c-grid displacement field of the contract.""" + itk = _elastix_api() + fixed, moving = fixed_moving(ct, t1c) + outdir = os.path.join(workdir or d("data/m6reg/elastix"), "itkelastix_def") + os.makedirs(outdir, exist_ok=True) + base = { + "ImageDimension": 3, + "FixedInternalImagePixelType": "float", + "MovingInternalImagePixelType": "float", + "Metric": "AdvancedMattesMutualInformation", + "NumberOfHistogramBins": 32, + "ImageSampler": "Random", + "ImageSamplerPatchSize": [25, 25, 25], + "Optimizer": "AdaptiveStochasticGradientDescent", + "MaximumNumberOfIterations": [2000, 2000, 2000], + "WriteResultImage": False, + } + stage1 = dict( + base, Transform="EulerTransform", TransformBase="ElastixTransformBase", + FixedCenteredTransform=True, HowToCombineTransforms="Compose", + AutomaticTransformInitialization=True, + AutomaticTransformInitializationMethod="CenterOfGravity", + AutomaticScalesEstimation=True) + stage2 = dict( + base, Transform="BSplineTransform", TransformBase="ElastixTransformBase", + HowToCombineTransforms="Compose", SplineOrder=3, GridSpacing=16, + AutomaticScalesEstimation=True) + reg_cls = itk.ElastixRegistrationMethod[itk.Image[itk.F, 3], + itk.Image[itk.F, 3]] + + def _run(params, t0=None): + reg = reg_cls.New() + reg.SetFixedImage(_sitk_to_itk(fixed, itk)) + reg.SetMovingImage(_sitk_to_itk(moving, itk)) + reg.SetParameterObject(itk.dict_to_parameter_object(params)) + if t0: + reg.SetInitialTransformParameterFileName(t0) + reg.SetOutputDirectory(outdir) + reg.SetNumberOfThreads(threads) + reg.LogToConsoleOff() + buf = io.StringIO() + try: + with contextlib.redirect_stdout(buf): + reg.Update() + except Exception as e: # noqa + tail = "\n".join(buf.getvalue().splitlines()[-10:]) + raise RuntimeError(f"elastix deformable stage failed: {e}\n{tail}") from e + return reg + + import time, io, contextlib + t0 = time.time() + _run(stage1) + ct_comb = _run(stage2, t0=os.path.join(outdir, "TransformParameters.0.txt")) \ + .GetCombinationTransform() + dt = time.time() - t0 + T = _itk_rasterize_to_tdf(ct_comb, fixed) + try: + import importlib.metadata + version = f"ITK-Elastix {importlib.metadata.version('itk-elastix')}" + except Exception: # noqa + version = "ITK-Elastix (python API)" + ha = sitk.GetArrayFromImage(t1c) + gap = _skull_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": "elastix", "toolkit": version, + "metric": "MattesMutualInformation (both stages, full volumes)", + "transform": ("Deformable (EulerTransform 2000x3 + BSplineTransform " + "16 mm grid, 2000 iters x 3 resolutions, -t0 init)"), + "convention": "T1c->CT (combination transform)", "auto_inverted": False, + "fixed": _modality_name(fixed), "moving": _modality_name(moving), + "seconds_reg": round(dt, 1), "skull_ring_gap_mm": gap, + } + return T, meta + + +def _impact_param_map_deformable(gpu, iters, samples, voxel_sizes): + """The ImpactReg ``MR_CT_TS`` reference parameter map in its native + DEFORMABLE form: the MIND + TotalSegmentator static feature combination + (same per-level models / scales / layer weights as the rigid benchmark) + driven by the reference's own RecursiveBSplineTransform (final grid + 14 mm, spline order 3) with the reference's three-metric combination -- + IMPACT (w 1.0) + Mattes MI (w 0.4) + TransformBendingEnergyPenalty + (w 10). No image-resampling pyramid (all-1 schedules; the scale comes + from the per-level ImpactVoxelSize). Two adaptations for the M6 + head-crop / head-neck-FOV mismatch (same as the rigid benchmark): + CenterOfGravity automatic transform initialization (the reference's + GeometricalCenter start leaves the heads ~100 mm apart and aborts) and + RequiredRatioOfValidSamples 0.1 (vs the 0.25 build default) as a + divergence gate. + """ + n = len(voxel_sizes) + if isinstance(iters, (int,)): + iters_l = str(iters) + else: + iters_l = " ".join(f"{i:g}" for i in iters) + mixed = "true" if gpu >= 0 else "false" + L = ["// IMPACT semantic-similarity metric (vboussot/ImpactLoss),"] + L.append("// MR_CT_TS reference map (deformable: RecursiveBSplineTransform)") + L.append('(Metric "Impact" "AdvancedMattesMutualInformation" "TransformBendingEnergyPenalty")') + L.append("(Metric0Weight 1)") + L.append(f"(Metric1Weight {IMPACT_MI_WEIGHT})") + L.append("(Metric2Weight 10)") + 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")', + '(CheckNumberOfSamples "true")', + "(RequiredRatioOfValidSamples 0.1)", + "", + "// deformable RecursiveBSpline (the reference map's native transform)", + '(Transform "RecursiveBSplineTransform")', + '(TransformBase "ElastixTransformBase")', + '(HowToCombineTransforms "Compose")', + '(BSplineTransformSplineOrder 3)', + '(FinalGridSpacingInPhysicalUnits 14)', + '(UseCyclicTransform "false")', + # adaptation: align the intensity centroids (the reference has + # AutomaticTransformInitialization "false" + GeometricalCenter, which + # leaves the heads ~100 mm apart for the head-crop / head-neck FOV + # mismatch and aborts) + '(AutomaticTransformInitialization "true")', + '(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_deformable(ct, t1c, gpu=None, iters=IMPACT_DEFAULT_ITERS, + samples=IMPACT_DEFAULT_SAMPLES, threads=12, + keep_scratch=False): + """IMPACT deformable CT->T1c (fixed=T1c, moving=CT) in the Docker image + ``vboussot/elastix_impact``: the ImpactReg ``MR_CT_TS`` reference parameter + map in its native deformable form (see + :func:`_impact_param_map_deformable` for the exact config and the two M6 + adaptations). The result -- a RecursiveBSplineTransform with the + CenterOfGravity init folded in -- is evaluated on the T1c grid + (:func:`_bspline_txt_to_tdf`) to give the contract displacement field.""" + import time + + ok, note = _impact_docker_ok() + if not ok: + raise PipelineUnavailable(f"impact: {note}") + model_rels = _impact_models_host() + + fixed, moving = fixed_moving(ct, t1c) + + if gpu is None: + ngpu = _n_gpus() + gpu = (os.getpid() % ngpu) if ngpu else -1 + + base = float(np.median(moving.GetSpacing())) + n = len(IMPACT_LEVELS) + voxel_sizes = [base * scale for _, scale, *_ in IMPACT_LEVELS] + + tmp = tempfile.mkdtemp(prefix="impact_reg_def_") + datadir = os.path.join(tmp, "Data") + outdir = os.path.join(tmp, "Out") + os.makedirs(datadir, exist_ok=True) + os.makedirs(outdir, exist_ok=True) + try: + 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_deformable(gpu, iters, samples, voxel_sizes) + with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: + f.write(param_text) + + 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_DEFORMABLE_{os.getpid()}.txt") + with open(param_dump, "w") as f: + f.write(param_text) + + cmd = ["docker", "run", "--rm"] + if gpu >= 0: + cmd += ["--gpus", "all"] + cmd += [ + "-v", f"{datadir}:/Data", + "-v", f"{outdir}:/Out", + "-v", f"{d(IMPACT_MODEL_CACHE)}:/Data/Models", + IMPACT_IMAGE, IMPACT_ELASTIX_BIN, + "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha", + "-p", "/Data/ParameterMap.txt", + "-out", "/Out/", "-threads", str(threads), + ] + t0 = time.time() + proc = subprocess.run(cmd, capture_output=True, text=True, timeout=14400) + dt = time.time() - t0 + if proc.returncode != 0: + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"impact deformable: docker elastix failed (rc={proc.returncode}):\n{tail}") + + itk_tf = os.path.join(outdir, "TransformParameters.0.itk.txt") + if not os.path.exists(itk_tf): + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"impact deformable: no TransformParameters.0.itk.txt produced:\n{tail}") + + T = _bspline_txt_to_tdf(itk_tf, fixed) + finally: + if not keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + + models = " + ".join(model_rels) + meta = { + "pipeline": "impact", + "toolkit": f"Docker {IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", + "metric": (f"IMPACT semantic ({models}, L1+Dice, Static) + " + f"Mattes MI (w {IMPACT_MI_WEIGHT}) + bending energy (w 10)"), + "transform": (f"Deformable (RecursiveBSpline 14 mm final grid, ASGD, " + f"{n} scales, iters {list(iters) if not isinstance(iters, int) else iters}, " + f"auto center init; " + f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'})"), + "convention": "T1c->CT (combination transform)", "auto_inverted": False, + "fixed": _modality_name(fixed), "moving": _modality_name(moving), + "seconds_reg": round(dt, 1), + "docker_image": IMPACT_IMAGE, "model": models, "gpu": gpu, + "iters": list(iters) if not isinstance(iters, int) else iters, + "samples": samples, "param_dump": param_dump, + } + if keep_scratch: + meta["scratch"] = tmp + return T, meta + + +def register_niftyreg_deformable(ct, t1c, threads=24, levels=3, maxit=150, + spline_spacing=8.0, keep_scratch=False, + backend="auto", gpu=None, clip_ct=(0.0, 80.0)): + """NiftyReg deformable CT->T1c (reference=T1c, floating=CT): the toolkit's + canonical non-rigid workflow -- reg_aladin rigid init (same budget as the + rigid benchmark: 3-level pyramid, 15 TLS iters/level) followed by reg_f3d + (F3D B-spline Free-Form deformation, Modat et al. 2010: NMI metric, + symmetric bending/linear regularisation, 3-level pyramid, 150 iters at + the final level, 8 mm final grid spacing). The CT is clipped to + [0, 80] HU by default (per the rigid A/B; ``clip_ct=None`` disables). + + Read-back: reg_f3d writes only the control-point grid (outputCPP.nii); + reg_transform -disp evaluates it to a dense displacement field on the + reference grid (the reg_aladin affine is folded into the control points + before optimization, so the dense field is the full transform). + NiftyReg reads RAS+ NIfTIs in its naive (x/y-mirrored) world, so the + stored deltas carry flipped x/y components -- + :func:`_flip_naive_world_components` restores the true RAS field + (verified to reproduce reg_f3d's own resampling to NCC ~ 1.0). The + result is the contract displacement field on the T1c grid. + """ + import time + + if backend not in ("auto", "conda", "docker"): + raise PipelineUnavailable( + f"niftyreg: unknown backend {backend!r} (use auto|conda|docker)") + dok, dnote = _niftyreg_docker_ok() + use_docker = backend == "docker" or (backend == "auto" and dok) + if use_docker: + if not dok: + raise PipelineUnavailable(f"niftyreg: {dnote}") + binpath = None + if gpu is None: + ngpu = _n_gpus() + gpu = (os.getpid() % ngpu) if ngpu else -1 + use_gpu = gpu >= 0 + platf = "1" if use_gpu else "0" + toolkit = (f"NiftyReg 2.1.1 (reg_aladin + reg_f3d, Docker " + f"{NIFTYREG_DOCKER_IMAGE}, CUDA)" if use_gpu + else f"NiftyReg 2.1.1 (reg_aladin + reg_f3d, Docker {NIFTYREG_DOCKER_IMAGE})") + else: + binpath = _niftyreg_bin() + gpu, use_gpu, platf = -1, False, None + toolkit = _niftyreg_version(binpath) + " (reg_f3d)" + + fixed, moving = fixed_moving(ct, t1c) + fixed_name = _modality_name(fixed) + moving_name = _modality_name(moving) + fixed, moving, ct_clipped = _maybe_clip_ct_pair(fixed, moving, clip_ct) + tmp = tempfile.mkdtemp(prefix="niftyreg_f3d_") + dt = 0.0 + try: + sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), + os.path.join(tmp, "ref.nii.gz"), True) + sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), + os.path.join(tmp, "flo.nii.gz"), True) + + def _nr(bin_name, flags): + if use_docker: + cmd = ["docker", "run", "--rm", + "--entrypoint", f"/niftyreg/{bin_name}"] + if use_gpu: + cmd += ["--gpus", "all"] + cmd += ["-v", f"{tmp}:/work", "-w", "/work", + NIFTYREG_DOCKER_IMAGE] + flags + cwd = None + else: + cmd = [os.path.join(os.path.dirname(binpath), bin_name)] + flags + cwd = tmp + t0 = time.time() + proc = subprocess.run(cmd, capture_output=True, text=True, + cwd=cwd, timeout=14400) + d = time.time() - t0 + if proc.returncode != 0: + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"niftyreg {bin_name} failed (rc={proc.returncode}):\n{tail}") + return d + + aladin_flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz", + "-rigOnly", "-aff", "aff.txt", + "-ln", str(levels), "-maxit", "15", + "-omp", str(threads)] + if platf is not None: + aladin_flags += ["-platf", platf] + if use_gpu: + aladin_flags += ["-gpuid", str(gpu)] + aladin_flags += ["-voff"] + t0 = time.time() + _nr("reg_aladin", aladin_flags) + + f3d_flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz", + "-aff", "aff.txt", + "-sx", str(spline_spacing), "-maxit", str(maxit), + "-ln", str(levels), "-omp", str(threads)] + if platf is not None: + f3d_flags += ["-platf", platf] + if use_gpu: + f3d_flags += ["-gpuid", str(gpu)] + f3d_flags += ["-voff"] + _nr("reg_f3d", f3d_flags) + + _nr("reg_transform", + ["-ref", "ref.nii.gz", "-disp", "outputCPP.nii", "spl.def"]) + dt = time.time() - t0 + + spl = sitk.ReadImage(os.path.join(tmp, "spl.def")) + arr = sitk.GetArrayFromImage(spl) + if arr.ndim == 4 and arr.shape[0] == 3: + arr = np.transpose(arr, (1, 2, 3, 0)) + # NiftyReg's naive world flips the x/y delta components relative to + # the true RAS field (see _flip_naive_world_components) + T = _field_transform(_flip_naive_world_components(arr), fixed) + finally: + if not keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + + ha = sitk.GetArrayFromImage(t1c) + gap = _skull_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": "niftyreg", "toolkit": toolkit, + "backend": "docker" if use_docker else "conda", + "metric": "reg_aladin (NCC block matching + TLS) -> reg_f3d (NMI, F3D B-spline)", + "transform": (f"Deformable (reg_f3d F3D, {levels}-level pyramid, " + f"{maxit} iters/final level, {spline_spacing:g} mm grid, " + "rigid init 15 TLS iters/level)"), + "convention": "T1c->CT (verified vs reg_f3d own resample)", "auto_inverted": False, + "fixed": fixed_name, "moving": moving_name, + "seconds_reg": round(dt, 1), + "levels": levels, "maxit": maxit, "spline_spacing": spline_spacing, + "gpu": gpu, "skull_ring_gap_mm": gap, + "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None, + } + if keep_scratch: + meta["scratch"] = tmp + return T, meta + + +def _bspline_representing_affine(fixed_img, A, t, center, grid_spacing_mm): + """A BSpline transform on the fixed image whose control points reproduce + the affine p -> A (p - center) + center + t exactly (cubic B-splines + reproduce linear maps). Parameter layout (verified against + TransformToDisplacementField): component-major blocks [all x, all y, all + z offsets], x-fastest within a block.""" + sp = np.array(fixed_img.GetSpacing()) + size = np.array(fixed_img.GetSize()) + cells = np.maximum(2, np.round(size * sp / np.array(grid_spacing_mm)).astype(int)) + bs = sitk.BSplineTransformInitializer(fixed_img, tuple(int(c) for c in cells), 3) + fp = np.array(bs.GetFixedParameters()) + mesh = fp[0:3].astype(int) # control grid incl. 3 boundary cells + g_origin = fp[3:6]; g_sp = fp[6:9] + ix = np.mgrid[0:mesh[0], 0:mesh[1], 0:mesh[2]] + P = np.stack([g_origin[0] + ix[0] * g_sp[0], + g_origin[1] + ix[1] * g_sp[1], + g_origin[2] + ix[2] * g_sp[2]], axis=-1) # (mx,my,mz,3) + off = np.einsum('ij,...j->...i', A, P - center) + center + t - P + bs.SetParameters(np.concatenate([off[..., 0].ravel(order="F"), + off[..., 1].ravel(order="F"), + off[..., 2].ravel(order="F")]).tolist()) + return bs + + +def register_simpleitk_deformable(ct, t1c, threads=12, iters_full=2000, + iters_bspline=50, + grid_spacing=(24.0, 24.0, 24.0)): + """SimpleITK two-stage deformable (the toolkit's canonical deformable + demo): stage 1 is the rigid benchmark's Euler3D map (Mattes MI 32 bins, + pyramid [4,2,1], iters_full per level, identity start); stage 2 is a + BSplineTransform3D (grid ~``grid_spacing`` mm, ``iters_bspline`` per + level on pyramid [4,2]) whose control points are initialised to + reproduce stage 1 exactly, so the returned transform is the full + composite. It is rasterised to the T1c-grid displacement field of the + contract. + + The 24 mm grid (5577 parameters) is the accelerated benchmark setting + (the 2026-10-02 original used 16 mm / 13005 parameters, ~48 min/case): + this build's spline per-iteration cost scales ~linearly with the + control-point count (~2.8 s/iter at the shrink-4 level and ~20 s/iter + at shrink-2 for 24 mm vs ~6.4 / ~45.8 s for 16 mm; independent of the + metric sample count), so 24 mm runs ~2.4x faster in isolation. Quality + (2 cases measured, same machine): the median case loses 0.003 MI + (0.6921 vs 0.6954, ImpactLoss improves 0.1547 vs 0.1552) and the + hardest case 0.017 (0.8008 vs 0.8176), keeping simpleitk's per-case + wins; the faster 32 mm grid (0.7902 / 0.7963 with 100 fine iters) + drops below antspyx's 0.7994 on the hardest case, so it was rejected. + A full-resolution polish pass (10 iters at shrink-1) was also tested + and dropped the achieved MI (0.6891 at 32 mm), so the pyramid stays + [4,2]. LBFGS2 was tested and fails to move the transform on this noisy + MI (line search rejects every step). See + REGISTRATION_BENCHMARK_DEFORMABLE.md.""" + m6 = _load_m6() + _sitk_set_threads(threads) + ct = sitk.Cast(ct, sitk.sitkFloat32); t1c = sitk.Cast(t1c, sitk.sitkFloat32) + fixed, moving = fixed_moving(ct, t1c) + shape = np.array(ct.GetSize()) + sh = [s for s in [4, 2, 1] if min(shape) / s >= 6] or [1] + sm = [2, 1, 0][:len(sh)] + import time + t0 = time.time() + r1 = sitk.ImageRegistrationMethod() + r1.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32) + r1.SetOptimizerAsRegularStepGradientDescent( + learningRate=10.0, minStep=0.2, numberOfIterations=iters_full, + relaxationFactor=0.5, gradientMagnitudeTolerance=1e-6) + r1.SetOptimizerScalesFromPhysicalShift() + r1.SetShrinkFactorsPerLevel(sh); r1.SetSmoothingSigmasPerLevel(sm) + r1.SetInterpolator(sitk.sitkLinear) + init1 = sitk.Euler3DTransform() + init1.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint( + ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2, + (fixed.GetSize()[2] - 1) // 2)))) + init1.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) + r1.SetInitialTransform(init1, inPlace=True) + tr1 = r1.Execute(fixed, moving) + # stage 2: BSpline initialised to reproduce stage 1 exactly (its control + # points start at the stage-1 pose, so the returned transform is the full + # Euler + BSpline composite); inPlace=True keeps the result a single + # transform (inPlace=False would return a CompositeTransform) + A1 = np.array(tr1.GetMatrix()).reshape(3, 3) + t1 = np.array(tr1.GetTranslation()) + c1 = np.array(tr1.GetCenter()) + init2 = _bspline_representing_affine(fixed, A1, t1, c1, grid_spacing) + # [4,2]: a full-resolution polish pass was measured to drop the achieved + # MI (0.6891 vs 0.6947 at 32 mm) and is omitted + sh2 = sh[:2]; sm2 = sm[:2] + r2 = sitk.ImageRegistrationMethod() + r2.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32) + r2.SetOptimizerAsRegularStepGradientDescent( + learningRate=10.0, minStep=0.2, numberOfIterations=iters_bspline, + relaxationFactor=0.5, gradientMagnitudeTolerance=1e-6) + # NOTE: no SetOptimizerScalesFromPhysicalShift() here -- it costs one + # metric evaluation per spline parameter (~2700 at 32 mm; ~15 min/case + # at the old 13k-parameter 16 mm grid); unit scales + the adaptive step + # size are the standard high-DOF setup. + r2.SetShrinkFactorsPerLevel(sh2); r2.SetSmoothingSigmasPerLevel(sm2) + r2.SetInterpolator(sitk.sitkLinear) + r2.SetInitialTransform(init2, inPlace=True) + tr2 = r2.Execute(fixed, moving) + dt = time.time() - t0 + T = _rasterize_to_tdf(tr2, fixed) + ha = sitk.GetArrayFromImage(t1c) + gap = m6._skull_ring_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": "simpleitk", + "toolkit": f"SimpleITK {sitk.Version_VersionString()}", + "metric": "MattesMutualInformation (both stages, full volumes)", + "transform": (f"Deformable (Euler3D {iters_full} iters/level + " + f"BSpline {tuple(float(x) for x in grid_spacing)} mm, " + f"{iters_bspline} iters/level on {sh2}; Euler on {sh})"), + "convention": "T1c->CT (stage 2 carries stage 1)", "auto_inverted": False, + "fixed": _modality_name(fixed), "moving": _modality_name(moving), + "seconds_reg": round(dt, 1), "pyramid_shrink": sh, "pyramid_smooth": sm, + "skull_ring_gap_mm": gap, + } + return T, meta + + +def register_synthmorph_deformable(ct, t1c, gpu=None, threads=None, extent=256, + keep_scratch=False, clip_ct=(0.0, 80.0)): + """SynthMorph ``joint`` model (affine + deformable in one symmetric + network pass -- the toolkit's canonical deformable) in the Docker image + ``freesurfer/synthmorph``. The CT is clipped to [0, 80] HU by default + (per the rigid A/B; ``clip_ct=None`` disables). The CLI is driven as + ``register -m joint moving fixed`` (moving = clipped CT first, fixed = + T1c second); the ``-t`` output is the displacement field on the fixed + (T1c) grid in surfa's naive (x/y-mirrored) world -- + :func:`_flip_naive_world_components` restores the true RAS field, giving + the contract transform directly. The tool's own ``-o`` resampling is + written as well and cross-checked against the field (NCC in meta).""" + import time + + ok, note = _synthmorph_docker_ok() + if not ok: + raise PipelineUnavailable(f"synthmorph: {note}") + + if gpu is None: + ngpu = _n_gpus() + gpu = (os.getpid() % ngpu) if ngpu else -1 + + tmp = tempfile.mkdtemp(prefix="synthmorph_def_") + dt = 0.0 + try: + fixed, moving = fixed_moving(ct, t1c) # fixed=T1c, moving=CT + fixed_name = _modality_name(fixed) + moving_name = _modality_name(moving) + fixed, moving, ct_clipped = _maybe_clip_ct_pair(fixed, moving, clip_ct) + # CLI order is (moving, fixed): moving = clipped CT, fixed = T1c + sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), + os.path.join(tmp, "moving.nii.gz"), True) + sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), + os.path.join(tmp, "fixed.nii.gz"), True) + + cmd = ["docker", "run", "--rm", + "--user", f"{os.getuid()}:{os.getgid()}", + "-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"] + if gpu >= 0: + cmd += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"] + cmd += ["-v", f"{tmp}:/data", SYNTHMORPH_IMAGE, + "register", "-m", "joint"] + if gpu >= 0: + cmd.append("-g") + if threads: + cmd += ["-j", str(threads)] + cmd += ["-e", str(extent), + "-t", "/data/field.nii.gz", + "-o", "/data/warped.nii.gz", + "/data/moving.nii.gz", "/data/fixed.nii.gz"] + t0 = time.time() + proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200) + dt = time.time() - t0 + if proc.returncode != 0: + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"synthmorph joint: docker failed (rc={proc.returncode}):\n{tail}") + if not os.path.exists(os.path.join(tmp, "field.nii.gz")): + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"synthmorph joint: no /data/field.nii.gz produced:\n{tail}") + + F = sitk.ReadImage(os.path.join(tmp, "field.nii.gz")) + arr = sitk.GetArrayFromImage(F) + if arr.ndim == 4 and arr.shape[-1] != 3: + arr = np.transpose(arr, (1, 2, 3, 0)) + T = _field_transform(_flip_naive_world_components(arr), fixed) + # cross-check against the tool's own resampling (clipped CT on the + # T1c grid) -- a convention regression would show up here + crosscheck = None + if os.path.exists(os.path.join(tmp, "warped.nii.gz")): + w_own = sitk.GetArrayFromImage(sitk.ReadImage(os.path.join(tmp, "warped.nii.gz"))) + w_mine = sitk.GetArrayFromImage( + sitk.Resample(moving, fixed, T, sitk.sitkLinear)) + ha_f = sitk.GetArrayFromImage(fixed) + from src.common import head_mask_from_image + head = head_mask_from_image(ha_f.astype("float32")) + if head.sum() > 0: + crosscheck = round(float(np.corrcoef( + w_mine[head].ravel(), w_own[head].ravel())[0, 1]), 4) + finally: + if not keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + + ha = sitk.GetArrayFromImage(t1c) + gap = _skull_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": "synthmorph", + "toolkit": f"Docker {SYNTHMORPH_IMAGE} (TensorFlow, pretrained)", + "metric": "SynthMorph joint (anatomy-aware contrast-invariant, deep, pretrained)", + "transform": (f"Deformable (joint affine+deform, identity start, " + f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, " + f"extent {extent})"), + "convention": "T1c->CT (verified vs -o resample)", "auto_inverted": False, + "fixed": fixed_name, "moving": moving_name, + "seconds_reg": round(dt, 1), + "docker_image": SYNTHMORPH_IMAGE, "model": "joint", "gpu": gpu, + "extent": extent, + "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None, + "own_resample_ncc_head": crosscheck, "skull_ring_gap_mm": gap, + } + if keep_scratch: + meta["scratch"] = tmp + return T, meta + + +REGISTERS_DEFORMABLE = { + "antspyx": register_antspyx_deformable, + "elastix": register_elastix_deformable, + "impact": register_impact_deformable, + "niftyreg": register_niftyreg_deformable, + "simpleitk": register_simpleitk_deformable, + "synthmorph": register_synthmorph_deformable, +} +ALL_PIPELINES_DEFORMABLE = list(REGISTERS_DEFORMABLE) + + def availability(): """pipeline -> (available: bool, note: str) for the current env.""" out = {} diff --git a/src/registration_pipelines_deformable.py b/src/registration_pipelines_deformable.py new file mode 100644 index 0000000..60b93d6 --- /dev/null +++ b/src/registration_pipelines_deformable.py @@ -0,0 +1,1039 @@ +"""Six DEFORMABLE CT -> T1c registration pipelines behind the unified contract. + +Companion to ``src/registration_pipelines.py`` (rigid). Same problem, same +contract, same shared rules (no mask / crop / pad, full volumes, T1c fixed / +CT moving); each pipeline here is its toolkit's native *deformable* +refinement initialized from the toolkit's own native rigid result, so the +rigid vs deformable difference is isolated to the deformable stage: + +- ``register_antspyx_deformable`` : ANTs Rigid -> SyNOnly (initialized + from the rigid transform), output read back as the ANTs ``Composite.h5``. +- ``register_elastix_deformable`` : Elastix 2-stage + (EulerTransform -> BSplineTransform), in-process python API; the + combination transform is densified onto the T1c grid into a + DisplacementFieldTransform (elastix's ``AdvancedBSplineTransform`` is not + readable by SimpleITK). +- ``register_impact_deformable`` : IMPACT with the *deformable* ImpactReg + reference map (RecursiveBSplineTransform + MIND/TotalSegmentator static + features + TransformBendingEnergyPenalty, i.e. the reference as published, + with the same CoG-init / sample-gate adaptations as the rigid pipeline); + densified via ``transformix`` inside the Docker image on coordinate + volumes. +- ``register_niftyreg_deformable`` : reg_aladin (rigid) -> reg_f3d + (free-form B-spline); the f3d deformation is applied with NiftyReg's own + ``reg_resample -cpp`` onto coordinate volumes (no CPP-file convention + guessing). +- ``register_simpleitk_deformable`` : SimpleITK 2-stage (Euler3D rigid -> + BSplineTransform on the T1c grid), the moving image pre-warped by the + rigid so the spline refines instead of re-estimating the gross pose. +- ``register_synthmorph_deformable`` : SynthMorph ``joint`` (affine + + deformable deep model); the saved displacement field (LAS axis convention) + is converted to the RAS contract field. + +Every function returns ``(T_final, T_base, meta)``: + +- ``T_final`` : a SimpleITK transform (DisplacementFieldTransform, affine or + CompositeTransform) satisfying the unified contract, i.e. + ``sitk.Resample(ct, t1c, T_final, sitk.sitkLinear)`` places the CT + concentric with the MR head (T1c -> CT mapping). +- ``T_base`` : the toolkit's own rigid / matrix result for the same pair + (the deformable stage's initialization), same contract. The benchmark + uses it for the "deformable increment" displacement and the rigid-vs- + deformable MI/ImpactLoss delta. +- ``meta`` : dict with the usual keys (pipeline / toolkit / metric / + transform / convention / auto_inverted / seconds_reg ...) plus + ``"deformable": True`` and pipeline-specific settings. + +Orientation / direction is verified empirically per case by the skull-ring +concentricity gate (see :func:`gate_dft`), the same philosophy as the rigid +benchmark's ``fix_orientation``. +""" +import os +import re +import shutil +import subprocess +import tempfile +import time + +import numpy as np +import SimpleITK as sitk + +from src.common import d +from src import registration_pipelines as RP + + +# --------------------------------------------------------------------------- # +# shared read-back helpers +# --------------------------------------------------------------------------- # + +def grid_world(ref): + """World (RAS mm) coordinates of every voxel of ``ref``: (z,y,x,3).""" + a = sitk.GetArrayFromImage(ref) + sp = np.array(ref.GetSpacing()) + org = np.array(ref.GetOrigin()) + Dm = np.array(ref.GetDirection()).reshape(3, 3) + zz, yy, xx = np.meshgrid(np.arange(a.shape[0]), + np.arange(a.shape[1]), + np.arange(a.shape[2]), indexing='ij') + return (np.stack([xx, yy, zz], -1) * sp) @ Dm.T + org + + +def coord_volumes(ref): + """Three scalar images on ``ref``'s grid whose value at a point is the + point's world coordinate (x, y, z).""" + w = grid_world(ref) + out = [] + for i in range(3): + comp = sitk.GetImageFromArray(w[..., i].astype('float32')) + comp.CopyInformation(ref) + out.append(comp) + return out + + +def padded_ref(ref, margin_mm=40.0): + """``ref``'s grid extended by ``margin_mm`` in every direction (same + spacing / origin lattice). Used to build coordinate volumes large enough + that a contract transform's target positions never hit the FOV boundary + (Resample would extrapolate garbage there). + + Returns ``(padded_image, (nx, ny, nz) padding in image-axis order)``.""" + import math + sp = ref.GetSpacing() # (x, y, z) + size = ref.GetSize() # (x, y, z) + n = tuple(int(math.ceil(margin_mm / sp[i])) for i in range(3)) # (x,y,z) + img = sitk.Image((size[0] + 2 * n[0], size[1] + 2 * n[1], + size[2] + 2 * n[2]), ref.GetPixelID()) + img.SetSpacing(list(sp)) + img.SetOrigin((ref.GetOrigin()[0] - n[0] * sp[0], + ref.GetOrigin()[1] - n[1] * sp[1], + ref.GetOrigin()[2] - n[2] * sp[2])) + img.SetDirection(ref.GetDirection()) + return img, n + + +def densify_sitk(T, ref, margin_mm=40.0): + """Densify any SimpleITK transform T (world frame) onto ``ref``'s grid. + + The coordinate volumes are built on a grid padded by ``margin_mm`` so T's + target positions stay inside the sampling FOV (no extrapolation + artifacts at the head boundary). Returns ``(dft, disp, H, world)`` where + ``H`` (z,y,x,3) is T evaluated at every ``ref`` grid point (world mm) and + ``disp = H - world`` is the displacement of the + :class:`~SimpleITK.DisplacementFieldTransform` ``dft``. + """ + pad, n_xyz = padded_ref(ref, margin_mm) + nx, ny, nz = n_xyz + cv = coord_volumes(pad) + H = np.empty(grid_world(ref).shape, dtype='float64') + for i in range(3): + out = sitk.GetArrayFromImage( + sitk.Resample(cv[i], pad, T, sitk.sitkLinear)) + # array axes are (z, y, x); crop the central ref-sized block + H[..., i] = out[nz:nz + ref.GetSize()[2], ny:ny + ref.GetSize()[1], + nx:nx + ref.GetSize()[0]] + world = grid_world(ref) + disp = H - world + dft_img = sitk.GetImageFromArray(disp.astype('float64')) + dft_img.CopyInformation(ref) + return sitk.DisplacementFieldTransform(dft_img), disp, H, world + + +def _itk_mod(): + try: + import itk + except Exception as e: # noqa + raise RP.PipelineUnavailable( + f"itk not importable ({e!r}); `pip install itk itk-elastix`") + return itk + + +def densify_itk(ctf, ref, margin_mm=40.0): + """Densify a live ITK-python transform (e.g. elastix's combination + transform, an itk CompositeTransform) onto ``ref``'s grid, like + :func:`densify_sitk` (coordinate volumes padded to avoid FOV-boundary + extrapolation).""" + itk = _itk_mod() + world = grid_world(ref) + pad, _n = padded_ref(ref, margin_mm) + ref_itk = RP._sitk_to_itk(ref, itk) + cv = coord_volumes(pad) + H = np.empty(world.shape, dtype='float32') + for i in range(3): + comp_itk = RP._sitk_to_itk(cv[i], itk) + itf = itk.LinearInterpolateImageFunction[itk.Image[itk.F, 3], + itk.D].New() + itf.SetInputImage(comp_itk) + # reference = the (unpadded) ref grid; the padded coordinate volume + # covers T's target positions -> no extrapolation artifacts + h = itk.resample_image_filter(comp_itk, reference_image=ref_itk, + use_reference_image=True, transform=ctf, + interpolator=itf, default_pixel_value=0.0) + H[..., i] = np.array(h) + H = H.astype('float64') + disp = H - world + dft_img = sitk.GetImageFromArray(disp.astype('float64')) + dft_img.CopyInformation(ref) + return sitk.DisplacementFieldTransform(dft_img), disp, H, world + + +def gate_dft(ct, t1c, cand): + """Skull-ring sanity gate for a deformable candidate transform. + + Deformable transforms (DFT / composite) are not invertible through + SimpleITK, so the rigid benchmark's T-vs-T^-1 gate does not apply; + instead the candidate is accepted only if the warped CT's skull ring is + concentric with the MR head (a divergent / mis-directed field blows the + gap up). Returns ``(cand, gap, score)``; raises RuntimeError when the + gate fails so a broken read-back never silently passes.""" + ha = sitk.GetArrayFromImage(t1c).astype('float32') + gap = RP._skull_gap(RP._warp_ct(ct, t1c, cand), ha) + score = RP._gap_score(gap) + if not gap or gap.get("std_mm", 999) > 25.0 or abs(gap.get("mid_mm", 999)) > 40.0: + raise RuntimeError( + f"deformable orientation gate failed: gap={gap} " + "(transform does not place the CT concentric with the MR head)") + return cand, gap, score + + +def _affine_of(T): + """(R, t) of any SimpleITK transform by point mapping (world frame).""" + o = np.array(T.TransformPoint((0.0, 0.0, 0.0)), dtype=float) + R = np.empty((3, 3), dtype=float) + for i in range(3): + p = (0.0, 0.0, 0.0) + p = p[:i] + (1.0,) + p[i + 1:] + R[:, i] = np.array(T.TransformPoint(p), dtype=float) - o + return R, o + + +# --------------------------------------------------------------------------- # +# 1. ANTsPy (Rigid -> SyNOnly) +# --------------------------------------------------------------------------- # + +def register_antspyx_deformable(ct, t1c, threads=16, + syn_iterations=(100, 100, 100)): + """ANTs deformable CT->T1c: the benchmark rigid (Mattes MI, + 500/500/200) followed by a ``SyNOnly`` stage initialized from that rigid + transform (``initial_transform``), so the SyN field refines the pose + instead of re-estimating it. The output is read back as ANTs' + ``Composite.h5`` (``write_composite_transform=True``), a single + SimpleITK-readable CompositeTransform. + + SyN budget: ``syn_iterations`` per level over the native SyN pyramid + (shrink 4/2/1, smoothing 2/1/0, ``SyN[0.2,3,0]`` -- antspyx's native + defaults for the SyN stage).""" + try: + import ants + except Exception as e: # noqa + raise RP.PipelineUnavailable(f"antspyx: {e!r}") + # antspyx parallelizes via OpenMP; cap threads so parallel workers share + # the box (the rigid benchmark left it at the OMP default) + os.environ["OMP_NUM_THREADS"] = str(threads) + fx, mv = RP.fixed_moving(ct, t1c) + af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)) + am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) + + t0 = time.time() + r1 = ants.registration(fixed=af, moving=am, type_of_transform='Rigid', + aff_metric='mattes', + aff_iterations=(500, 500, 200), + aff_shrink_factors=(4, 2, 1), + aff_smoothing_sigmas=(3, 2, 0)) + rig_path = r1['fwdtransforms'][0] if isinstance(r1['fwdtransforms'], + list) else r1['fwdtransforms'] + T_base, inv_base = RP.fix_orientation(ct, t1c, sitk.ReadTransform(rig_path)) + base_s = time.time() - t0 + + t0 = time.time() + r2 = ants.registration(fixed=af, moving=am, type_of_transform='SyNOnly', + initial_transform=rig_path, + reg_iterations=tuple(syn_iterations), + write_composite_transform=True) + T_final = sitk.ReadTransform(r2['fwdtransforms']) + # the composite includes the -r init; direction-check empirically + try: + T_final, inv_final = RP.fix_orientation(ct, t1c, T_final) + except Exception: # noqa: composite inverse unavailable -> gate only + T_final, _, _ = gate_dft(ct, t1c, T_final) + inv_final = False + syn_s = time.time() - t0 + meta = { + "pipeline": "antspyx", "deformable": True, + "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", + "metric": "Mattes MI (rigid 500/500/200) -> SyN Mattes MI (native SyN pyramid)", + "transform": (f"Rigid (auto COM init) + SyNOnly {tuple(syn_iterations)} " + f"(init = rigid; composite read-back)"), + "convention": "T1c->CT (verified)", "auto_inverted": bool(inv_final), + "fixed": RP._modality_name(fx), "moving": RP._modality_name(mv), + "seconds_rigid": round(base_s, 1), "seconds_syn": round(syn_s, 1), + "seconds_reg": round(base_s + syn_s, 1), + "syn_iterations": list(syn_iterations), + } + return T_final, T_base, meta + + +# --------------------------------------------------------------------------- # +# 2. ITK-Elastix (Euler -> BSplineTransform) +# --------------------------------------------------------------------------- # + +def _elastix_deformable_params(): + """Two-stage parameter maps: the benchmark rigid stage (Euler, 2000 x3) + then a BSplineTransform stage (1000 x3, grid 100 -> 50 -> 25 mm). The + 25 mm final grid (~9x9x6 control points) keeps the refinement smooth; a + finer 15 mm grid over-deformed (folded, disp p95 ~138 mm) in the pilot. + (A TransformBendingEnergyPenalty second metric is not usable here: it + segfaults in the ITK-Elastix 0.25.4 python API's ASGD parameter + estimation, single- or multi-metric.)""" + stage1 = { + "ImageDimension": 3, + "FixedInternalImagePixelType": "float", + "MovingInternalImagePixelType": "float", + "Transform": "EulerTransform", + "TransformBase": "ElastixTransformBase", + "FixedCenteredTransform": True, + "TransformParameterReadoutFormat": 1, + "HowToCombineTransforms": "Compose", + "AutomaticTransformInitialization": True, + "AutomaticTransformInitializationMethod": "CenterOfGravity", + "AutomaticScalesEstimation": True, + "Optimizer": "AdaptiveStochasticGradientDescent", + "MaximumNumberOfIterations": [2000, 2000, 2000], + "Metric": "AdvancedMattesMutualInformation", + "NumberOfHistogramBins": 32, + "ImageSampler": "Random", + "ImageSamplerPatchSize": [25, 25, 25], + "WriteResultImage": False, + } + stage2 = { + "ImageDimension": 3, + "FixedInternalImagePixelType": "float", + "MovingInternalImagePixelType": "float", + "Transform": "BSplineTransform", + "TransformBase": "EulerTransform", + "HowToCombineTransforms": "Compose", + "BSplineTransformSplineOrder": 3, + # grid 100 -> 50 -> 25 mm over the 3 resolutions: a 25 mm final grid + # (~9x9x6 control points) keeps the refinement smooth and fold-free; + # a 15 mm grid without regularization collapsed (folds + MI loss) + "InitialBSplineGridSpacing": [100, 100, 100], + "DeformationPrecision": "float", + "BSplineInterpolationOrder": 3, + "Optimizer": "AdaptiveStochasticGradientDescent", + "MaximumNumberOfIterations": [1000, 1000, 1000], + "Metric": "AdvancedMattesMutualInformation", + "NumberOfHistogramBins": 32, + "ImageSampler": "Random", + "ImageSamplerPatchSize": [25, 25, 25], + "WriteResultImage": False, + } + return [stage1, stage2] + + +def register_elastix_deformable(ct, t1c, workdir=None, threads=8): + """Elastix deformable CT->T1c, two stages in one run: + EulerTransform (benchmark rigid budget) -> BSplineTransform + (2000 iters x 3 resolutions, ~15 mm final control spacing). The + combination transform (an elastix ``AdvancedBSplineTransform``, not + readable by SimpleITK) is densified onto the T1c grid via the + coordinate-volume trick and returned as a DisplacementFieldTransform.""" + itk = RP._elastix_api() + fixed, moving = RP.fixed_moving(ct, t1c) + outdir = os.path.join(workdir or d("data/m6reg/elastix"), "itkelastix_def") + os.makedirs(outdir, exist_ok=True) + import io, contextlib + t0 = time.time() + reg_cls = itk.ElastixRegistrationMethod[itk.Image[itk.F, 3], itk.Image[itk.F, 3]] + reg = reg_cls.New() + reg.SetFixedImage(RP._sitk_to_itk(fixed, itk)) + reg.SetMovingImage(RP._sitk_to_itk(moving, itk)) + reg.SetParameterObject(itk.dict_to_parameter_object(_elastix_deformable_params())) + reg.SetOutputDirectory(outdir) + reg.SetNumberOfThreads(threads) + reg.LogToConsoleOff() + buf = io.StringIO() + try: + with contextlib.redirect_stdout(buf): + reg.Update() + except Exception as e: # noqa + tail = "\n".join(buf.getvalue().splitlines()[-10:]) + raise RuntimeError(f"elastix deformable registration failed: {e}\n{tail}") from e + dt = time.time() - t0 + ctf = reg.GetCombinationTransform() + T_final, _, _, _ = densify_itk(ctf, t1c) + T_final, _, _ = gate_dft(ct, t1c, T_final) + # base = the toolkit's own rigid stage (benchmark config) + t0 = time.time() + reg2 = reg_cls.New() + reg2.SetFixedImage(RP._sitk_to_itk(fixed, itk)) + reg2.SetMovingImage(RP._sitk_to_itk(moving, itk)) + reg2.SetParameterObject( + itk.dict_to_parameter_object([_elastix_deformable_params()[0]])) + reg2.SetOutputDirectory(outdir) + reg2.SetNumberOfThreads(threads) + reg2.LogToConsoleOff() + buf2 = io.StringIO() + with contextlib.redirect_stdout(buf2): + reg2.Update() + base_cand = RP._itk_transform_to_sitk_rigid(reg2.GetCombinationTransform()) + T_base, inv_base = RP.fix_orientation(ct, t1c, base_cand) + base_s = time.time() - t0 + meta = { + "pipeline": "elastix", "deformable": True, + "toolkit": "ITK-Elastix 0.25.4 (python API)", + "metric": "Mattes MI (both stages; ASGD)", + "transform": ("EulerTransform rigid (2000x3) + BSplineTransform " + "(2000x3, grid 60->15 mm)"), + "convention": "T1c->CT (verified)", "auto_inverted": False, + "fixed": RP._modality_name(fixed), "moving": RP._modality_name(moving), + "seconds_rigid_base": round(base_s, 1), + "seconds_reg": round(dt, 1), + } + return T_final, T_base, meta + + +# --------------------------------------------------------------------------- # +# 3. IMPACT (reference RecursiveBSpline map) +# --------------------------------------------------------------------------- # + +# The published ImpactReg MR_CT_TS reference map (deformable): the +# RecursiveBSplineTransform + MIND/TotalSegmentator static feature +# combination + TransformBendingEnergyPenalty, exactly as downloaded from +# HF VBoussot/ImpactReg (Elastix_IMPACT_Static/ParameterMap_Static.txt). +# Adapted for M6 the same way the rigid benchmark adapted it: absolute +# /Data/Models paths, per-worker ImpactGPU, CenterOfGravity auto-init (the +# reference's identity start does not survive the M6 gross offset) and the +# 0.1 sample-validity gate for the head-crop vs head+neck-FOV mismatch. +IMPACT_DEFORMABLE_REF = """ +(MaximumNumberOfIterations 400 300 200 200) +(NumberOfSpatialSamples 2000) +(Transform "RecursiveBSplineTransform") +(NumberOfResolutions 4) +(FinalGridSpacingInPhysicalUnits 14) +(FixedImagePyramid "FixedGenericImagePyramid") +(MovingImagePyramid "MovingGenericImagePyramid") +(FixedImagePyramidRescaleSchedule 1 1 1 1 1 1 1 1 1 1 1 1) +(MovingImagePyramidRescaleSchedule 1 1 1 1 1 1 1 1 1 1 1 1) +(ImpactModelsPath0 "/Data/Models/MIND/R1D2_3D.pt" "/Data/Models/TS/M852.pt") +(ImpactDimension0 3 3) +(ImpactNumberOfChannels0 1 1) +(ImpactPatchSize0 0 0 0 0 0 0) +(ImpactVoxelSize0 6 6 6 6 6 6) +(ImpactLayersMask0 "1" "0000001") +(ImpactSubsetFeatures0 32 64) +(ImpactPCA0 0 0) +(ImpactDistance0 "L1" "Dice") +(ImpactLayersWeight0 0.2 0.8) +(ImpactModelsPath1 "/Data/Models/MIND/R1D2_3D.pt" "/Data/Models/TS/M852.pt") +(ImpactDimension1 3 3) +(ImpactNumberOfChannels1 1 1) +(ImpactPatchSize1 0 0 0 0 0 0) +(ImpactVoxelSize1 3 3 3 3 3 3) +(ImpactLayersMask1 "1" "0000001") +(ImpactSubsetFeatures1 32 64) +(ImpactPCA1 0 0) +(ImpactDistance1 "L1" "Dice") +(ImpactLayersWeight1 0.3 0.7) +(ImpactModelsPath2 "/Data/Models/MIND/R1D2_3D.pt" "/Data/Models/TS/M850.pt") +(ImpactDimension2 3 3) +(ImpactNumberOfChannels2 1 1) +(ImpactPatchSize2 0 0 0 0 0 0) +(ImpactVoxelSize2 2 2 2 2 2 2) +(ImpactLayersMask2 "1" "00000001") +(ImpactSubsetFeatures2 32 64) +(ImpactPCA2 0 0) +(ImpactDistance2 "L1" "Dice") +(ImpactLayersWeight2 0.6 0.4) +(ImpactModelsPath3 "/Data/Models/MIND/R1D2_3D.pt" "/Data/Models/TS/M850.pt") +(ImpactDimension3 3 3) +(ImpactNumberOfChannels3 1 1) +(ImpactPatchSize3 0 0 0 0 0 0) +(ImpactVoxelSize3 2 2 2 2 2 2) +(ImpactLayersMask3 "1" "00000001") +(ImpactSubsetFeatures3 32 64) +(ImpactPCA3 0 0) +(ImpactDistance3 "L1" "Dice") +(ImpactLayersWeight3 0.7 0.3) +(ImpactUseMixedPrecision "true") +(ImpactFeaturesMapUpdateInterval -1) +(ImpactMode "Static") +(ImpactGPU 0) +(Metric "Impact" "AdvancedMattesMutualInformation" "TransformBendingEnergyPenalty") +(Metric0Weight 1) +(Metric1Weight 0.4) +(Metric2Weight 10) +(FixedInternalImagePixelType "float") +(MovingInternalImagePixelType "float") +(UseDirectionCosines "true") +(Registration "MultiMetricMultiResolutionRegistration") +(BSplineTransformSplineOrder 3) +(UseCyclicTransform "false") +(AutomaticTransformInitialization "true") +(AutomaticTransformInitializationMethod "CenterOfGravity") +(AutomaticScalesEstimation "true") +(HowToCombineTransforms "Compose") +(FixedCenteredTransform "true") +(RequiredRatioOfValidSamples 0.1) +(Optimizer "AdaptiveStochasticGradientDescent") +(MaximumNumberOfSamplingAttempts 8) +(UseAdaptiveStepSizes "true") +(UseMultiThreadingForMetrics "true") +(ASGDParameterEstimationMethod "DisplacementDistribution") +(SigmoidInitialTime 0.0) +(NoiseCompensation "true") +(NumberOfSamplesForExactGradient 4096) +(AutomaticParameterEstimation "true") +(FixedKernelBSplineOrder 3) +(MovingKernelBSplineOrder 3) +(CheckNumberOfSamples "true") +(UseRelativeWeights "false") +(WriteTransformParametersEachIteration "false") +(WriteTransformParametersEachResolution "false") +(ShowExactMetricValue "false") +(ErodeFixedMask "false") +(ErodeMovingMask "false") +(UseBinaryFormatForTransformationParameters "false") +(Interpolator "BSplineInterpolator") +(ImageSampler "RandomCoordinate") +(NewSamplesEveryIteration "true") +(UseRandomSampleRegion "false") +(ResampleInterpolator "FinalBSplineInterpolator") +(FinalBSplineInterpolationOrder 3) +(BSplineInterpolationOrder 3) +(Resampler "DefaultResampler") +(WriteIterationInfo "false") +(WriteResultImage "false") +(DefaultPixelValue -1024) +(ResultImageFormat "mha") +(ITKTransformOutputFileNameExtension "itk.txt") +(WriteITKCompositeTransform "true") +""" + + +def register_impact_deformable(ct, t1c, gpu=None, iters=RP.IMPACT_DEFAULT_ITERS, + samples=RP.IMPACT_DEFAULT_SAMPLES, threads=12, + keep_scratch=False): + """IMPACT deformable CT->T1c: the published ImpactReg ``MR_CT_TS`` + reference map (RecursiveBSplineTransform, final grid 14 mm, MIND + + TotalSegmentator static features, 4 scales 6/3/2/2 mm, ASGD + 400/300/200/200, 2000 samples, + bending-energy penalty weight 10), + with the M6 adaptations of the rigid pipeline (CenterOfGravity auto + init, ``RequiredRatioOfValidSamples 0.1``, per-worker GPU). + + The output ``RecursiveBSplineTransform`` is not readable by SimpleITK, + so the transform is densified *inside the Docker image*: three + coordinate volumes on the T1c grid are resampled with ``transformix`` + (the reference stage file's ``ResultImagePixelType`` patched to float) + and the result is wrapped in a DisplacementFieldTransform.""" + import time as _t + ok, note = RP._impact_docker_ok() + if not ok: + raise RP.PipelineUnavailable(f"impact: {note}") + model_rels = RP._impact_models_host() + fixed, moving = RP.fixed_moving(ct, t1c) + if gpu is None: + ngpu = RP._n_gpus() + gpu = (os.getpid() % ngpu) if ngpu else -1 + + # the reference map with per-worker GPU + budget overrides + n = len(RP.IMPACT_LEVELS) + if isinstance(iters, int): + iters_l = str(iters) + else: + iters_l = " ".join(f"{i:g}" for i in iters) + txt = IMPACT_DEFORMABLE_REF.replace("(MaximumNumberOfIterations 400 300 200 200)", + f"(MaximumNumberOfIterations {iters_l})") + txt = txt.replace("(NumberOfSpatialSamples 2000)", + f"(NumberOfSpatialSamples {samples})") + txt = txt.replace("(ImpactGPU 0)", f"(ImpactGPU {gpu})") + mixed = "true" if gpu >= 0 else "false" + txt = txt.replace('(ImpactUseMixedPrecision "true")', + f'(ImpactUseMixedPrecision "{mixed}")') + # per-level voxel sizes scaled to the moving spacing (as the rigid does) + base = float(np.median(moving.GetSpacing())) + for i, (_, scale, *_) in enumerate(RP.IMPACT_LEVELS): + vs = base * scale + txt = re.sub(rf"\(ImpactVoxelSize{i} .*?\)", + f"(ImpactVoxelSize{i} {vs:g} {vs:g} {vs:g} {vs:g} {vs:g} {vs:g})", + txt) + + tmp = tempfile.mkdtemp(prefix="impact_def_reg_") + datadir = os.path.join(tmp, "Data") + outdir = os.path.join(tmp, "Out") + os.makedirs(datadir, exist_ok=True) + os.makedirs(outdir, exist_ok=True) + t0 = _t.time() + 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) + with open(os.path.join(datadir, "ParameterMap.txt"), "w") as f: + f.write(txt) + dump_dir = d("data/m6reg/impact/paramdump") + os.makedirs(dump_dir, exist_ok=True) + param_dump = os.path.join(dump_dir, f"ParameterMap_DEFORMABLE_{os.getpid()}.txt") + with open(param_dump, "w") as f: + f.write(txt) + + # coordinate volumes on the T1c grid + cv = coord_volumes(fixed) + for i, comp in enumerate(cv): + sitk.WriteImage(comp, os.path.join(datadir, f"coord_{i}.mha"), True) + + cmd = ["docker", "run", "--rm"] + if gpu >= 0: + cmd += ["--gpus", "all"] + cmd += [ + "-v", f"{datadir}:/Data", + "-v", f"{outdir}:/Out", + "-v", f"{d(RP.IMPACT_MODEL_CACHE)}:/Data/Models", + RP.IMPACT_IMAGE, RP.IMPACT_ELASTIX_BIN, + "-f", "/Data/Fixed_image.mha", "-m", "/Data/Moving_image.mha", + "-p", "/Data/ParameterMap.txt", + "-out", "/Out/", "-threads", str(threads), + ] + proc = subprocess.run(cmd, capture_output=True, text=True, timeout=21600) + if proc.returncode != 0: + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"impact deformable: docker elastix failed " + f"(rc={proc.returncode}):\n{tail}") + stage_txt = os.path.join(outdir, "TransformParameters.0.txt") + if not os.path.exists(stage_txt): + raise RuntimeError("impact deformable: no TransformParameters.0.txt produced") + # patch the resampler output type (short -> float) for the coordinate vols + stage_f = open(stage_txt).read().replace('(ResultImagePixelType "short")', + '(ResultImagePixelType "float")') + with open(os.path.join(datadir, "stage_float.txt"), "w") as f: + f.write(stage_f) + # densify with transformix (3 scalar coord vols, one container call) + for i in range(3): + os.makedirs(os.path.join(outdir, f"c{i}"), exist_ok=True) + chain = " && ".join( + f"/lib/elastix-install/bin/transformix " + f"-tp /Data/stage_float.txt -in /Data/coord_{i}.mha -out /Out/c{i}/" + for i in range(3)) + cmd2 = ["docker", "run", "--rm", + "-v", f"{datadir}:/Data", "-v", f"{outdir}:/Out", + RP.IMPACT_IMAGE, "/bin/bash", "-c", chain] + proc2 = subprocess.run(cmd2, capture_output=True, text=True, timeout=1800) + if proc2.returncode != 0: + tail = "\n".join((proc2.stdout + proc2.stderr).splitlines()[-15:]) + raise RuntimeError(f"impact deformable: transformix failed (rc={proc2.returncode}):\n{tail}") + H = np.empty(grid_world(fixed).shape, dtype='float32') + for i in range(3): + H[..., i] = sitk.GetArrayFromImage( + sitk.ReadImage(os.path.join(outdir, f"c{i}", "result.mha"))) + world = grid_world(fixed) + disp = (H - world).astype('float64') + dft_img = sitk.GetImageFromArray(disp) + dft_img.CopyInformation(fixed) + T_final = sitk.DisplacementFieldTransform(dft_img) + gate_dft(ct, t1c, T_final) + dt = _t.time() - t0 + finally: + if not keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + + # base = the toolkit's own rigid stage (same metric family / init) + t0 = _t.time() + T_base, base_meta = RP.register_impact(ct, t1c, gpu=gpu, iters=iters, + samples=samples, threads=threads) + base_s = _t.time() - t0 + meta = { + "pipeline": "impact", "deformable": True, + "toolkit": f"Docker {RP.IMPACT_IMAGE} (Elastix 5.2.0 + IMPACT)", + "metric": ("IMPACT semantic (MIND+TS static, L1+Dice) + Mattes MI (w 0.4) " + "+ TransformBendingEnergyPenalty (w 10) -- the published " + "MR_CT_TS reference map"), + "transform": (f"RecursiveBSplineTransform (final grid 14 mm, 4 scales, " + f"ASGD {iters_l if isinstance(iters, int) else iters}, " + f"{samples} samples; CoG auto-init)"), + "convention": "T1c->CT (verified)", "auto_inverted": False, + "fixed": RP._modality_name(fixed), "moving": RP._modality_name(moving), + "seconds_rigid_base": round(base_s, 1), "seconds_reg": round(dt, 1), + "docker_image": RP.IMPACT_IMAGE, "gpu": gpu, + "iters": list(iters) if not isinstance(iters, int) else iters, + "samples": samples, "param_dump": param_dump, + } + return T_final, T_base, meta + + +# --------------------------------------------------------------------------- # +# 4. NiftyReg (reg_aladin -> reg_f3d, densified via reg_resample -cpp) +# --------------------------------------------------------------------------- # + +def register_niftyreg_deformable(ct, t1c, threads=24, levels=3, maxit=15, + f3d_maxit=150, backend="auto", gpu=None, + clip_ct=(0.0, 80.0), keep_scratch=False): + """NiftyReg deformable CT->T1c: ``reg_aladin`` (benchmark rigid config, + symmetric block matching + TLS) then ``reg_f3d`` (free-form B-spline, + NMI, default regularization) initialized from the aladin affine. + + The f3d output is the control-point grid (``outputCPP.nii``); the full + mapping is densified with NiftyReg's own ``reg_resample -cpp`` applied to + three coordinate volumes on the CT grid (no CPP convention guessing -- + the tool applies its own field), yielding the T1c->CT contract mapping + as a DisplacementFieldTransform.""" + import time as _t + dok, dnote = RP._niftyreg_docker_ok() + use_docker = dok if backend == "auto" else backend == "docker" + if not use_docker: + raise RP.PipelineUnavailable( + "niftyreg deformable: the benchmark uses the v2.1.1 Docker image " + f"(reg_f3d + reg_resample); {dnote}") + if gpu is None: + ngpu = RP._n_gpus() + gpu = (os.getpid() % ngpu) if ngpu else -1 + use_gpu = gpu >= 0 + + fixed, moving = RP.fixed_moving(ct, t1c) + fixed_name = RP._modality_name(fixed) + moving_name = RP._modality_name(moving) + fixed, moving, ct_clipped = RP._maybe_clip_ct_pair(fixed, moving, clip_ct) + + tmp = tempfile.mkdtemp(prefix="niftyreg_def_reg_") + t0 = _t.time() + try: + sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), + os.path.join(tmp, "ref.nii.gz"), True) + sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), + os.path.join(tmp, "flo.nii.gz"), True) + flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz", + "-rigOnly", "-aff", "aladin.aff", + "-ln", str(levels), "-maxit", str(maxit), + "-omp", str(threads)] + if use_gpu: + flags += ["-interp", "1", "-platf", "1", "-gpuid", str(gpu)] + else: + flags += ["-interp", "3", "-platf", "0"] + flags += ["-voff"] + proc = _niftyreg_docker_run(tmp, "reg_aladin", flags, use_gpu) + if proc.returncode != 0: + raise RuntimeError(f"niftyreg: reg_aladin failed:\n" + + "\n".join(proc.stderr.splitlines()[-15:])) + aff_file = os.path.join(tmp, "aladin.aff") + A_cf = RP._parse_niftyreg_affine(aff_file) + T_base, inv_base = RP.fix_orientation(ct, t1c, A_cf) + + f3d_flags = ["-ref", "ref.nii.gz", "-flo", "flo.nii.gz", + "-aff", "aladin.aff", + "-ln", str(levels), "-maxit", str(f3d_maxit), + "-omp", str(threads)] + if use_gpu: + f3d_flags += ["-platf", "1", "-gpuid", str(gpu)] + else: + f3d_flags += ["-platf", "0"] + f3d_flags += ["-voff"] + proc = _niftyreg_docker_run(tmp, "reg_f3d", f3d_flags, use_gpu) + if proc.returncode != 0: + raise RuntimeError(f"niftyreg: reg_f3d failed:\n" + + "\n".join(proc.stderr.splitlines()[-15:])) + cpp = os.path.join(tmp, "outputCPP.nii") + if not os.path.exists(cpp): + raise RuntimeError("niftyreg: no outputCPP.nii produced") + + # densify the full (affine + f3d) mapping with reg_resample -cpp on + # coordinate volumes defined on the CT (floating) grid + cv = coord_volumes(moving) + H = np.empty(grid_world(fixed).shape, dtype='float32') + for i, comp in enumerate(cv): + cp = os.path.join(tmp, f"coord_{i}.nii.gz") + op = os.path.join(tmp, f"out_{i}.nii") + sitk.WriteImage(comp, cp, True) + rs_flags = ["-ref", "ref.nii.gz", "-flo", f"coord_{i}.nii.gz", + "-cpp", "outputCPP.nii", "-res", f"out_{i}.nii", "-voff"] + proc = _niftyreg_docker_run(tmp, "reg_resample", rs_flags, False) + if proc.returncode != 0: + raise RuntimeError(f"niftyreg: reg_resample failed:\n" + + "\n".join(proc.stderr.splitlines()[-15:])) + H[..., i] = sitk.GetArrayFromImage(sitk.ReadImage(op)) + world = grid_world(fixed) + disp = (H - world).astype('float64') + dft_img = sitk.GetImageFromArray(disp) + dft_img.CopyInformation(fixed) + T_final = sitk.DisplacementFieldTransform(dft_img) + gate_dft(ct, t1c, T_final) + dt = _t.time() - t0 + finally: + if not keep_scratch: + shutil.rmtree(tmp, ignore_errors=True) + + meta = { + "pipeline": "niftyreg", "deformable": True, + "toolkit": f"NiftyReg 2.1.1 (reg_aladin + reg_f3d, Docker, {'CUDA' if use_gpu else 'CPU'})", + "metric": "NCC block matching + TLS (rigid); NMI + bending energy (f3d)", + "transform": (f"reg_aladin -rigOnly ({levels}-level, {maxit} TLS iters) " + f"+ reg_f3d ({levels}-level, {f3d_maxit} iters, default " + f"5-voxel B-spline grid); densified via reg_resample -cpp"), + "convention": "T1c->CT (verified)", "auto_inverted": bool(inv_base), + "fixed": fixed_name, "moving": moving_name, + "seconds_reg": round(dt, 1), "gpu": gpu, + "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None, + } + return T_final, T_base, meta + + +def _niftyreg_docker_run(tmp, entrypoint, flags, use_gpu): + cmd = ["docker", "run", "--rm", "--entrypoint", entrypoint] + if use_gpu: + cmd += ["--gpus", "all"] + cmd += ["-v", f"{tmp}:/work", "-w", "/work", RP.NIFTYREG_DOCKER_IMAGE] + flags + return subprocess.run(cmd, capture_output=True, text=True, timeout=21600) + + +# --------------------------------------------------------------------------- # +# 5. SimpleITK (Euler3D -> BSplineTransform) +# --------------------------------------------------------------------------- # + +def register_simpleitk_deformable(ct, t1c, threads=16, iters_full=2000, + bs_iterations=(500, 500), + bs_grid_spacing_mm=25.0, + bs_sample_pct=0.001): + """SimpleITK deformable CT->T1c, two stages: the benchmark rigid + (Euler3D, Mattes MI, 2000 iters x [4,2,1]) then a BSplineTransform stage + (``bs_iterations`` per level over a [2,1] shrink schedule, ~25 mm control + spacing on the T1c grid, Mattes MI sampled on ``bs_sample_pct`` of the + T1c volume ~ 7600 stochastic samples). The moving image is pre-warped + by the rigid stage so the B-spline optimizes the residual deformation + around the identity; the returned transform is the composed mapping + T1c->CT as a DisplacementFieldTransform. + + The B-spline stage uses *manual* unit parameter scales (1 mm per control + displacement): ``SetOptimizerScalesFromPhysicalShift`` evaluates the + metric once per parameter (~3900 evaluations) and is prohibitively slow. + The 25 mm grid (vs a finer 15 mm one with 3x the control points) keeps + SimpleITK's per-iteration cost ~0.35 s and the refinement smooth. + """ + import time as _t + try: + sitk.CommonProperties.SetGlobalDefaultNumberOfThreads(threads) + except AttributeError: + pass + fixed, moving = RP.fixed_moving(ct, t1c) + fixed = sitk.Cast(fixed, sitk.sitkFloat32) + moving = sitk.Cast(moving, sitk.sitkFloat32) + shape = np.array(ct.GetSize()) + init = sitk.Euler3DTransform() + init.SetCenter(tuple(float(x) for x in fixed.TransformIndexToPhysicalPoint( + ((fixed.GetSize()[0] - 1) // 2, (fixed.GetSize()[1] - 1) // 2, + (fixed.GetSize()[2] - 1) // 2)))) + init.SetParameters((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) + t0 = _t.time() + r = sitk.ImageRegistrationMethod() + r.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32) + r.SetOptimizerAsRegularStepGradientDescent( + learningRate=10.0, minStep=0.2, numberOfIterations=iters_full, + relaxationFactor=0.5, gradientMagnitudeTolerance=1e-6) + r.SetOptimizerScalesFromPhysicalShift() + sh = [s for s in [4, 2, 1] if min(shape) / s >= 6] or [1] + sm = [2, 1, 0][:len(sh)] + r.SetShrinkFactorsPerLevel(sh) + r.SetSmoothingSigmasPerLevel(sm) + r.SetInterpolator(sitk.sitkLinear) + r.SetInitialTransform(init, inPlace=False) + tr1 = r.Execute(fixed, moving) + base_s = _t.time() - t0 + + # pre-warp the moving onto its own grid with the rigid: the residual is + # then purely deformable and the spline stage starts at the identity + moving_pre = sitk.Resample(moving, moving, tr1, sitk.sitkLinear) + + t1c_dims = np.array(fixed.GetSize()) * np.array(fixed.GetSpacing()) + mesh = tuple(max(4, int(round(dd / bs_grid_spacing_mm))) for dd in t1c_dims) + bs = sitk.BSplineTransform(3, 3) + bs.SetTransformDomainMeshSize(mesh) + bs.SetTransformDomainOrigin(fixed.GetOrigin()) + bs.SetTransformDomainPhysicalDimensions( + tuple(float(v) for v in np.array(mesh) * bs_grid_spacing_mm)) + bs.SetTransformDomainDirection(fixed.GetDirection()) + bs.SetParameters([0.0] * bs.GetNumberOfParameters()) + r2 = sitk.ImageRegistrationMethod() + r2.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32) + if bs_sample_pct < 1.0: + r2.SetMetricSamplingStrategy(1) # random + r2.SetMetricSamplingPercentage(bs_sample_pct) + r2.SetOptimizerAsRegularStepGradientDescent( + learningRate=2.0, minStep=0.1, numberOfIterations=bs_iterations[0], + relaxationFactor=0.5, gradientMagnitudeTolerance=1e-6) + # manual unit scales: 1 mm per control-point displacement (see docstring) + r2.SetOptimizerScales([1.0] * bs.GetNumberOfParameters()) + r2.SetShrinkFactorsPerLevel([2, 1]) + r2.SetSmoothingSigmasPerLevel([1, 0]) + r2.SetInterpolator(sitk.sitkLinear) + r2.SetInitialTransform(bs, inPlace=False) + t0 = _t.time() + T2 = r2.Execute(fixed, moving_pre) + bs_s = _t.time() - t0 + + # total contract = tr1 o T2 (T2 maps T1c points onto the CT grid, tr1 + # maps those grid points to CT world) + H2 = densify_sitk(T2, fixed)[2] + world = grid_world(fixed) + R1, o1 = _affine_of(tr1) + Htot = np.einsum('...i,ij->...j', H2, R1) + o1 + disp = Htot - world + dft_img = sitk.GetImageFromArray(disp.astype('float64')) + dft_img.CopyInformation(fixed) + T_final = sitk.DisplacementFieldTransform(dft_img) + gate_dft(ct, t1c, T_final) + meta = { + "pipeline": "simpleitk", "deformable": True, + "toolkit": f"SimpleITK {sitk.Version_VersionString()}", + "metric": "Mattes MI (both stages, full volumes)", + "transform": (f"Euler3D rigid ({iters_full} iters x {sh}) + " + f"BSplineTransform ({'/'.join(str(i) for i in bs_iterations)} " + f"iters/level x [2,1], {bs_sample_pct:.1%} volume sampling " + f"(~7600 stochastic samples), " + f"~{bs_grid_spacing_mm:g} mm grid on the T1c; moving " + f"pre-warped by the rigid)"), + "convention": "T1c->CT (verified)", "auto_inverted": False, + "fixed": RP._modality_name(fixed), "moving": RP._modality_name(moving), + "seconds_rigid_base": round(base_s, 1), "seconds_bs": round(bs_s, 1), + "seconds_reg": round(base_s + bs_s, 1), + "bs_mesh": list(mesh), "bs_iterations": list(bs_iterations), + "bs_sample_pct": bs_sample_pct, + } + return T_final, tr1, meta + + +# --------------------------------------------------------------------------- # +# 6. SynthMorph (joint affine + deformable) +# --------------------------------------------------------------------------- # + +def register_synthmorph_deformable(ct, t1c, gpu=None, extent=256, + clip_ct=(0.0, 80.0)): + """SynthMorph ``joint`` (affine + deformable deep model) CT->T1c in the + Docker image. The CLI pair is passed as (CT, T1c) so the model + registers CT -> T1c (the direction that aligns well on these pairs, per + the joint investigation in REGISTRATION_BENCHMARK_RIGID.md). The + saved displacement field (``-t``; on the T1c grid) is written in the + network's LAS axis convention: the RAS contract displacement is + ``(−x, −y, +z)`` of the file values (verified against the tool's own + resampling on two cases; an orientation gate over the sign variants + keeps the read-back robust). + + ``T_base`` is the model's own ``affine`` matrix result (same CLI + order, ``-t`` LTA, direction-checked by the usual gate).""" + import time as _t + ok, note = RP._synthmorph_docker_ok() + if not ok: + raise RP.PipelineUnavailable(f"synthmorph: {note}") + if gpu is None: + ngpu = RP._n_gpus() + gpu = (os.getpid() % ngpu) if ngpu else -1 + + fixed, moving = RP.fixed_moving(ct, t1c) # fixed=T1c, moving=CT + fixed_name = RP._modality_name(fixed) + moving_name = RP._modality_name(moving) + fixed, moving, ct_clipped = RP._maybe_clip_ct_pair(fixed, moving, clip_ct) + tmp = tempfile.mkdtemp(prefix="synthmorph_def_reg_") + try: + # CLI order (moving, fixed) = (CT, T1c): the model registers CT->T1c; + # -t saves the field on the T1c (tool "fixed") grid + sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), + os.path.join(tmp, "ct.nii.gz"), True) + sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), + os.path.join(tmp, "t1c.nii.gz"), True) + cmd = ["docker", "run", "--rm", + "--user", f"{os.getuid()}:{os.getgid()}", + "-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"] + if gpu >= 0: + cmd += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"] + cmd += ["-v", f"{tmp}:/data", RP.SYNTHMORPH_IMAGE, + "register", "-m", "joint", "-e", str(extent), + "-t", "/data/t_field.nii.gz", + "/data/ct.nii.gz", "/data/t1c.nii.gz"] + if gpu >= 0: + cmd.append("-g") + t0 = _t.time() + proc = subprocess.run(cmd, capture_output=True, text=True, timeout=7200) + if proc.returncode != 0: + tail = "\n".join((proc.stdout + proc.stderr).splitlines()[-25:]) + raise RuntimeError(f"synthmorph joint: docker failed (rc={proc.returncode}):\n{tail}") + if not os.path.exists(os.path.join(tmp, "t_field.nii.gz")): + raise RuntimeError(f"synthmorph joint: no t_field.nii.gz produced:\n" + + "\n".join((proc.stdout + proc.stderr).splitlines()[-15:])) + # 4D NIfTI: (3, z, y, x) -> (z, y, x, 3) + arr = sitk.GetArrayFromImage(sitk.ReadImage(os.path.join(tmp, "t_field.nii.gz"))) + if arr.shape[0] != 3 or arr.shape[1:] != sitk.GetArrayFromImage(t1c).shape: + raise RuntimeError(f"synthmorph joint: unexpected field shape {arr.shape}") + D = arr.transpose(1, 2, 3, 0).astype('float64') + # orientation gate over the sign variants (primary: the verified + # LAS->RAS correction (−x, −y, +z)) + cands = [ + (D * np.array([-1.0, -1.0, 1.0]), "las2ras"), + (D, "as-is"), + (-D, "negated"), + (D * np.array([1.0, 1.0, -1.0]), "las2ras-flipz"), + ] + best = None + for disp, tag in cands: + dft_img = sitk.GetImageFromArray(disp) + dft_img.CopyInformation(fixed) + T = sitk.DisplacementFieldTransform(dft_img) + try: + _, gap, score = gate_dft(ct, t1c, T) + except RuntimeError: + score = float("inf") + gap = None + if best is None or score < best[0]: + best = (score, gap, T, tag) + if best is None or not best[1]: + raise RuntimeError("synthmorph joint: all field sign variants failed " + "the orientation gate") + score, gap, T_final, gate_tag = best + dt = _t.time() - t0 + + # base: the model's affine matrix head (same CLI order); -t LTA with + # (CT, T1c) = tool (mov=CT, fix=T1c) -> LTA stores target->source = + # T1c->CT (contract); direction-checked by the usual affine gate + cmd2 = ["docker", "run", "--rm", + "--user", f"{os.getuid()}:{os.getgid()}", + "-e", "TF_FORCE_GPU_ALLOW_GROWTH=true"] + if gpu >= 0: + cmd2 += ["--gpus", "all", "-e", f"CUDA_VISIBLE_DEVICES={gpu}"] + cmd2 += ["-v", f"{tmp}:/data", RP.SYNTHMORPH_IMAGE, + "register", "-m", "affine", "-e", str(extent), + "-t", "/data/aff_lta.lta", + "/data/ct.nii.gz", "/data/t1c.nii.gz"] + if gpu >= 0: + cmd2.append("-g") + proc2 = subprocess.run(cmd2, capture_output=True, text=True, timeout=3600) + if proc2.returncode != 0 or not os.path.exists(os.path.join(tmp, "aff_lta.lta")): + raise RuntimeError("synthmorph affine base failed:\n" + + "\n".join((proc2.stdout + proc2.stderr).splitlines()[-15:])) + T_base, inv_base = RP.fix_orientation(ct, t1c, RP._lta_to_sitk(os.path.join(tmp, "aff_lta.lta"))) + finally: + shutil.rmtree(tmp, ignore_errors=True) + + meta = { + "pipeline": "synthmorph", "deformable": True, + "toolkit": f"Docker {RP.SYNTHMORPH_IMAGE} (TensorFlow, pretrained)", + "metric": "SynthMorph joint anatomy-aware contrast-invariant (deep, pretrained)", + "transform": (f"joint (affine + deformable, mid-space; " + f"device={'GPU' + str(gpu) if gpu >= 0 else 'CPU'}, " + f"extent {extent}); field read-back '{gate_tag}'"), + "convention": "T1c->CT (verified)", "auto_inverted": False, + "fixed": fixed_name, "moving": moving_name, + "seconds_reg": round(dt, 1), + "docker_image": RP.SYNTHMORPH_IMAGE, "model": "joint", "gpu": gpu, + "extent": extent, "gate_variant": gate_tag, + "clip_ct": [float(clip_ct[0]), float(clip_ct[1])] if ct_clipped else None, + } + return T_final, T_base, meta + + +REGISTERS = { + "antspyx": register_antspyx_deformable, + "elastix": register_elastix_deformable, + "impact": register_impact_deformable, + "niftyreg": register_niftyreg_deformable, + "simpleitk": register_simpleitk_deformable, + "synthmorph": register_synthmorph_deformable, +} +ALL_PIPELINES = list(REGISTERS) + + +def availability(): + """Same toolkits as the rigid benchmark: pipeline -> (available, note).""" + return RP.availability() \ No newline at end of file