From bec0e9834f9da643f7a69016e45e314ae2945247 Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Sat, 3 Oct 2026 13:26:12 +0800 Subject: [PATCH] feat(benchmark): add ANTs deformable variant pipelines Extend registration_pipelines and bench_registration with stock antspyx ANTs presets (SyNRA, SyN, Repro sr/br, QuickRepro sr/br) plus the legacy Rigid -> SyNOnly workflow, and update the deformable benchmark documentation with follow-up ANTs-variant runs and results. --- REGISTRATION_BENCHMARK_DEFORMABLE.md | 160 ++++++++++++++++- scripts/bench_registration.py | 13 +- src/registration_pipelines.py | 255 +++++++++++++++++++++++++++ 3 files changed, 418 insertions(+), 10 deletions(-) diff --git a/REGISTRATION_BENCHMARK_DEFORMABLE.md b/REGISTRATION_BENCHMARK_DEFORMABLE.md index 9743448..843973d 100644 --- a/REGISTRATION_BENCHMARK_DEFORMABLE.md +++ b/REGISTRATION_BENCHMARK_DEFORMABLE.md @@ -9,7 +9,12 @@ 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 +local). Follow-up ANTs-variant runs 2026-10-03 (08:44–09:12 local: five +stock antspyx presets — direct `SyNRA` / `SyN`, `Repro[sr|br]`, +`QuickRepro[sr|br]`; 11:00–11:26 local: the legacy two-call +`Rigid → SyNOnly` workflow; 12:44–12:49 local: re-run of the stock `SyN` +preset as `ants_syn`, replacing the initial `SyNOnly` row) on the same +10 cases (see the *ANTs variants* section). 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 @@ -143,8 +148,9 @@ for the case): ### 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 + 8 of 9 cases and within 0.002 of best on the other — + `m6_26TYQI3C`, 0.799 vs `simpleitk`'s 0.801) 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 @@ -213,10 +219,11 @@ for the case): 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. +`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 @@ -236,6 +243,135 @@ for the case): term. The same weak correlation with perceived alignment noted in the rigid benchmark, so it is reported for reference only. +## ANTs variants: direct SyNRA / SyN and the Repro / QuickRepro workflows + +The current `antspyx` row above is antspyx's stock **`SyN` preset** (Affine +2100×1200×1200×0 + SyN[0.2,3,0] 40×20×0, Mattes MI 32, center-of-mass init). +As a follow-up (2026-10-03, 6 workers as the main run) seven further +ANTs configurations were benchmarked on the same 10 cases and the same +contract — the question: is there a better stock ANTs deformable workflow +for CT→T1c than the `SyN` preset? + +| pipeline | antspyx preset | stages (all COM init unless noted) | +|---|---|---| +| **ants_rigid_synonly** | legacy two-call (`Rigid` + `SyNOnly`) | Rigid[0.25] (Mattes 32, 500×500×200, 4×2×1, 3×2×0) → SyNOnly SyN[0.2,3,0] (Mattes 32, 100×100×0, 4×2×1, 2×1×0) initialized from the Rigid stage; `Composite.h5` read-back | +| **ants_synra** | `SyNRA` | Rigid[0.25] + Affine[0.25] (each Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1) | +| **ants_syn** | `SyN` | Affine[0.25] (Mattes 32, 2100×1200×1200×0, 4×2×2×1, 3×2×1×0) + SyN[0.2,3,0] (Mattes 32, 40×20×0, 4×2×1); COM init — the stock preset, identical configuration to the existing `antspyx` row | +| **ants_synreprosr** | `antsRegistrationSyNRepro[sr]` | Rigid[0.1] (GC 1/1/Regular/0.25, 1000×500×250×100, 8×4×2×1, 3×2×1×0vox) + SyN[0.2,3,0] (CC 1/2, 50×20, 2×1, 1×0vox) | +| **ants_synreprobr** | `antsRegistrationSyNRepro[br]` | the `Repro` Rigid[0.1] stage + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×20, 2×1, 1×0vox) | +| **ants_synquickreprosr** | `antsRegistrationSyNQuickRepro[sr]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + SyN[0.2,3,0] (CC 1/2, 50×0, 2×1) | +| **ants_synquickreprobr** | `antsRegistrationSyNQuickRepro[br]` | Rigid[0.1] (GC, 1000×500×250×0, 8×4×2×1) + BSplineSyN[0.2,26,0,3] (CC 1/2, 50×0, 2×1) | + +`sr`/`br` select the deformable stage (SyN vs BSplineSyN, antspyx preset +naming — not a brain mask); these are the same stage configurations the +`antsRegistrationSyNRepro*.sh` / `antsRegistrationSyNQuickRepro*.sh` +workflows use. `ants_rigid_synonly` is the two-call workflow that an early +version of this benchmark used: a Rigid stage, then a +*deeper* SyNOnly stage (100 iterations at every level incl. full +resolution, vs the `SyN` preset's 40×20×0) initialized from the Rigid +stage's transform. Full volumes, no CT clip, each toolkit preset's own +metrics (GC/CC/Mattes as above). As for the other rows, the stage files +are rasterised to the T1c-grid displacement field and scored with the +shared code; the composite read-back was cross-checked per case against +ANTs' own resampling (**NCC 1.0 on all 63 scored cases**, +`own_warp_ncc_head` in the JSONs). One read-back subtlety: when ANTs runs Rigid + Affine as separate +linear stages (`SyNRA`), it writes the merged linear stage as a *binary* +`0GenericAffine.mat` that SimpleITK cannot parse; the pipeline recovers it +by point-probing the transform (exact for affine maps, verified to +2×10⁻⁴ mm). + +All six variants registered 9/9 (the 10th is the known incomplete-CT data +issue). **Bold** = best for the column; `antspyx` is the existing row, +re-listed for reference. + +| pipeline | registered | s/reg (med) | ImpactLoss (med/mean) | MI (med/mean, nats) | +|---|---|---|---|---| +| **ants_rigid_synonly** | 9 | 208.4 | 0.156 / 0.155 | **0.753 / 0.763** | +| **antspyx** (existing) | 9 | 31.7 | 0.158 / 0.156 | 0.723 / 0.729 | +| **ants_synra** | 9 | 43.6 | 0.158 / 0.156 | 0.722 / 0.729 | +| **ants_syn** | 9 | 33.8 | 0.158 / 0.156 | 0.722 / 0.729 | +| **ants_synreprosr** | 9 | 294.7 | **0.151 / 0.153** | 0.530 / 0.529 | +| **ants_synreprobr** | 9 | 370.3 | 0.156 / 0.157 | 0.541 / 0.537 | +| **ants_synquickreprosr** | 9 | 83.7 | 0.154 / 0.156 | 0.542 / 0.534 | +| **ants_synquickreprobr** | 9 | 96.0 | 0.158 / 0.158 | 0.549 / 0.538 | + +Per-case MI (nats, unclipped score; **bold** = best for the case): + +| case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | +|---|---|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | 0.679 | **0.741** | 0.681 | 0.681 | 0.623 | 0.646 | 0.630 | 0.627 | +| m6_23XIMO4M_20230222 | 0.723 | **0.760** | 0.722 | 0.722 | 0.530 | 0.544 | 0.529 | 0.530 | +| m6_24V62QQM_20230223 | 0.789 | **0.824** | 0.789 | 0.787 | 0.520 | 0.528 | 0.542 | 0.549 | +| m6_24V62QQM_20231017 | **0.720** | 0.716 | **0.720** | **0.720** | 0.473 | 0.475 | 0.469 | 0.462 | +| m6_24V62QQM_20240820 | 0.738 | **0.753** | 0.738 | 0.738 | 0.524 | 0.525 | 0.517 | 0.524 | +| m6_26TYQI3C_20221027 | 0.799 | **0.857** | 0.798 | 0.799 | 0.557 | 0.541 | 0.599 | 0.632 | +| m6_2BBN5XA4_20250411 | 0.665 | **0.696** | 0.665 | 0.666 | 0.230 | 0.303 | 0.235 | 0.297 | +| m6_2C4TKI6H_20240506 | 0.719 | **0.751** | 0.719 | 0.719 | 0.672 | 0.649 | 0.654 | 0.621 | +| m6_2CDNN7BH_20240102 | 0.726 | **0.765** | 0.726 | 0.725 | 0.629 | 0.622 | 0.631 | 0.599 | + +Per-case ImpactLoss (semantic, unclipped score; **bold** = best (lowest) for +the case): + +| case | antspyx | ants_rigid_synonly | ants_synra | ants_syn | ants_synreprosr | ants_synreprobr | ants_synquickreprosr | ants_synquickreprobr | +|---|---|---|---|---|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | incomplete-CT | +| m6_22IRRKU6_20240115 | 0.1541 | 0.1530 | 0.1542 | 0.1540 | **0.1442** | 0.1482 | 0.1474 | 0.1493 | +| m6_23XIMO4M_20230222 | 0.1581 | 0.1565 | 0.1581 | 0.1581 | **0.1505** | 0.1549 | 0.1535 | 0.1585 | +| m6_24V62QQM_20230223 | 0.1473 | **0.1464** | 0.1472 | 0.1472 | 0.1515 | 0.1553 | 0.1535 | 0.1560 | +| m6_24V62QQM_20231017 | 0.1518 | **0.1512** | 0.1516 | 0.1519 | 0.1559 | 0.1575 | 0.1577 | 0.1581 | +| m6_24V62QQM_20240820 | **0.1545** | **0.1545** | 0.1547 | 0.1546 | 0.1553 | 0.1578 | 0.1586 | 0.1586 | +| m6_26TYQI3C_20221027 | 0.1580 | 0.1562 | 0.1578 | 0.1579 | **0.1561** | 0.1574 | 0.1602 | 0.1610 | +| m6_2BBN5XA4_20250411 | 0.1624 | **0.1611** | 0.1624 | 0.1624 | 0.1680 | 0.1695 | 0.1690 | 0.1707 | +| m6_2C4TKI6H_20240506 | 0.1594 | 0.1582 | 0.1594 | 0.1593 | **0.1490** | 0.1526 | 0.1517 | 0.1545 | +| m6_2CDNN7BH_20240102 | 0.1590 | 0.1582 | 0.1589 | 0.1589 | **0.1493** | 0.1560 | 0.1538 | 0.1580 | + +### Interpretation (ANTs variants) + +- **`ants_rigid_synonly` is the best ANTs configuration of the eight** — + the only one that beats the current `antspyx` row: median MI **0.753** + (per-case best on 8 of 9; its 0.857 on `m6_26TYQI3C` is the highest MI + of the whole benchmark, ahead of the previous 0.801) at ImpactLoss + 0.156 (second-best median after `synthmorph`'s 0.156). The gain over + the `SyN` preset comes from the *deeper* SyN stage: 100 iterations at + every level **including full resolution** (100×100×0 vs 40×20×0) + refining the Rigid stage's pose. The cost is 208.4 s/case (~6.6× the + preset) — the full-resolution SyN iterations dominate. +- **`ants_synra` is indistinguishable from the current `antspyx` row** on + quality (median MI 0.722 vs 0.723, ImpactLoss 0.158 vs 0.158; per-case + deltas within 0.002 nats, ties on 4 of 9) at **+38 % wall time** + (43.6 s vs 31.7 s): the explicit Rigid stage ahead of the Affine stage + buys nothing — on these pairs the 12-DOF stage from the center-of-mass + init already finds the gross pose, so the extra 6-DOF stage just + re-solves it. +- **`ants_syn` reproduces the existing `antspyx` row** (median MI 0.722 + vs 0.723, ImpactLoss 0.158, s/reg 33.8 vs 31.7; per-case deltas within + 0.002 nats) — as expected, since it is the *same stock `SyN` preset* + (Affine + SyN[0.2,3,0], COM init) re-run in this follow-up batch: it + serves as a consistency check that the shared read-back / scoring path + is configuration-stable. +- **The Repro / QuickRepro family is semantically competitive but weaker on + intensity:** median MI 0.530–0.549 (vs 0.723 for the `SyN` preset), + driven by the GC-based Rigid stage and the short 2-level SyN/BSplineSyN + refinement (50×20 or 50×0 iters, 2×1 pyramid, CC 2-bin metric) — a + workflow tuned for intra-modality T1w-T1w. `ants_synreprosr` still takes + the **best ImpactLoss median of all thirteen pipelines here** (0.151; + per-case wins on 4 of 9), so the semantic/intensity split of the main + table carries over to the stock ANTs workflows. +- **Time tracks budget, not quality:** the full Repro Rigid stage + (1000×500×250×100, 8×4×2×1, GC) dominates — `ants_synreprobr` is the + slowest ANTs configuration (370 s/case) at the same (low) MI as the + 84–96 s `QuickRepro` variants (which drop the 4th rigid level and the + 2nd SyN level, changing the MI by < 0.02). Only `ants_rigid_synonly` + challenges the `SyN` preset's quality (better MI at 6.6× the time); + every other variant is worse for more or for the same. +- **Production skull-ring QA (reference only):** acceptance 4/9 + `ants_syn` / `ants_synra` / `ants_rigid_synonly` (= baseline), 2/9 + `ants_synreprobr`, 1/9 the rest; rotational dispersion worst in the + `Repro` / `QuickRepro` family (24–29 mm std) vs 14–15 mm for + `ants_syn` / `ants_synra` / `ants_rigid_synonly` — consistent with + their alignment. + ## Transform conventions (deformable-specific) Unlike the rigid benchmark (where the elastix-style toolkits came out in the @@ -287,6 +423,16 @@ python scripts/bench_registration.py --deformable --workers 6 --threads 24 --qa- # single pipeline / case subset python scripts/bench_registration.py --deformable --pipelines simpleitk impact --n-cases 3 +# the ANTs-variant follow-up (direct SyNRA / SyN + Repro / QuickRepro +# sr/br; results go to the same results/bench_registration_deformable/ dir, +# 2026-10-03 run, 6 workers, ~28 min wall; `ants_syn` is the re-run of the +# stock SyN preset that replaced the initial SyNOnly row) +python scripts/bench_registration.py --deformable --pipelines "ants_synra ants_syn ants_synreprosr ants_synreprobr ants_synquickreprosr ants_synquickreprobr" --workers 6 --threads 24 --qa-images + +# the legacy two-call Rigid -> SyNOnly(100,100,100) workflow (2026-10-03, +# 6 workers, ~25 min wall; dominated by the full-resolution SyN iters) +python scripts/bench_registration.py --deformable --pipelines "ants_rigid_synonly" --workers 6 --threads 24 --qa-images + # 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 diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index 819afc9..d7e91c7 100644 --- a/scripts/bench_registration.py +++ b/scripts/bench_registration.py @@ -2,6 +2,9 @@ Pipelines (see src/registration_pipelines.py for the shared transform contract): simpleitk | elastix | antspyx | synthmorph | impact | niftyreg + + the ANTs variants (deformable only): ants_synra | ants_syn | + ants_synreprosr | ants_synreprobr | ants_synquickreprosr | + ants_synquickreprobr Each pipeline registers a head-and-neck CT to a tightly-cropped T1c MRI (cross-modality, inter-session tilt up to ~30 deg) and the GTV is warped onto @@ -298,6 +301,10 @@ def _call(pipeline, ct, t1c, tv, args): threads=args.threads) if pipeline == "antspyx": return RP.register_antspyx_deformable(ct, t1c) + if pipeline in RP.ANT_VARIANTS: + return RP.register_ants_variant_deformable(ct, t1c, variant=pipeline) + if pipeline == "ants_rigid_synonly": + return RP.register_ants_rigid_synonly_deformable(ct, t1c) if pipeline == "synthmorph": return RP.register_synthmorph_deformable(ct, t1c, gpu=args.synthmorph_gpu, @@ -445,7 +452,7 @@ def _fmt(v, spec="{}", default="-"): def print_table(rows): - order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)] + order = [p for p in RP.ALL_PIPELINE_NAMES if any(r["pipeline"] == p for r in rows)] print("\n" + "=" * 96) hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} " f"{'s/reg':>7s} {'ImpactLoss':>10s} {'MI':>7s}") @@ -510,7 +517,7 @@ def write_plot(rows, out_png, label="rigid"): def merge(args): out_dir = os.path.join(args.out_dir) rows = [] - for p in RP.ALL_PIPELINES: + for p in RP.ALL_PIPELINE_NAMES: fp = os.path.join(out_dir, p + ".json") if os.path.exists(fp): with open(fp) as f: @@ -651,7 +658,7 @@ def main(): args.qa_out = os.path.abspath(args.qa_out or os.path.join(args.out_dir, "qa")) os.makedirs(args.qa_out, exist_ok=True) pipelines = [p.lower() for p in args.pipelines.split()] - unknown = set(pipelines) - set(RP.ALL_PIPELINES) + unknown = set(pipelines) - set(RP.ALL_PIPELINE_NAMES) if unknown: raise SystemExit(f"unknown pipeline(s): {unknown}") diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index 76a8747..3f394e9 100644 --- a/src/registration_pipelines.py +++ b/src/registration_pipelines.py @@ -1402,6 +1402,247 @@ def register_antspyx_deformable(ct, t1c): return T, meta +# --------------------------------------------------------------------------- # +# ANTs variants: direct SyNRA / SyN presets + the reproducible +# Repro / QuickRepro workflows (sr = SyN stage, br = BSplineSyN stage), all +# driven through antspyx's stock presets (see +# REGISTRATION_BENCHMARK_DEFORMABLE.md, "ANTs variants"). +# --------------------------------------------------------------------------- # + +# pipeline name -> (antspyx type_of_transform, stock-preset description) +ANT_VARIANTS = { + "ants_synra": ( + "SyNRA", + "direct SyNRA preset: Rigid[0.25] + Affine[0.25] (both Mattes 32, " + "2100x1200x1200x0, 4x2x2x1, 3x2x1x0) + SyN[0.2,3,0] (Mattes 32, " + "40x20x0, 4x2x1); COM init"), + "ants_syn": ( + "SyN", + "direct SyN (affine + deformable): Affine[0.25] (Mattes 32, " + "2100x1200x1200x0, 4x2x2x1, 3x2x1x0) + SyN[0.2,3,0] (Mattes 32, " + "40x20x0, 4x2x1, 2x1x0); COM init -- the stock preset, same " + "configuration as the existing antspyx row"), + "ants_synreprosr": ( + "antsRegistrationSyNRepro[sr]", + "Repro (sr): Rigid[0.1] (GC 1/1/Regular/0.25, 1000x500x250x100, " + "8x4x2x1, 3x2x1x0vox) + SyN[0.2,3,0] (CC 1/2, 50x20, 2x1, 1x0vox); " + "COM init"), + "ants_synreprobr": ( + "antsRegistrationSyNRepro[br]", + "Repro (br): Rigid[0.1] (GC 1/1/Regular/0.25, 1000x500x250x100, " + "8x4x2x1, 3x2x1x0vox) + BSplineSyN[0.2,26,0,3] (CC 1/2, 50x20, 2x1, " + "1x0vox); COM init"), + "ants_synquickreprosr": ( + "antsRegistrationSyNQuickRepro[sr]", + "QuickRepro (sr): Rigid[0.1] (GC 1/1/Regular/0.25, 1000x500x250x0, " + "8x4x2x1, 3x2x1x0vox) + SyN[0.2,3,0] (CC 1/2, 50x0, 2x1, 1x0vox); " + "COM init"), + "ants_synquickreprobr": ( + "antsRegistrationSyNQuickRepro[br]", + "QuickRepro (br): Rigid[0.1] (GC 1/1/Regular/0.25, 1000x500x250x0, " + "8x4x2x1, 3x2x1x0vox) + BSplineSyN[0.2,26,0,3] (CC 1/2, 50x0, 2x1, " + "1x0vox); COM init"), +} + + +def _ants_stage_to_sitk(path, ref_img): + """Read one ANTs stage file into a SimpleITK transform. A Warp.nii.gz + vector field becomes a DisplacementFieldTransform on ``ref_img``'s grid; + a .mat becomes its affine. When two linear stages (Rigid + Affine) are + run in one registration, ANTs writes them as a single *binary* + '0GenericAffine.mat' (the merged rigid-then-affine mapping); that format + is not SimpleITK-readable, so it is recovered by point-probing the ANTs + transform (exact for affine maps: verified to 2e-4 mm against ANTs' own + point mapping and NCC 1.0 against its warpedmovout).""" + base = os.path.basename(path) + if "Warp.nii.gz" in base and "Inverse" not in base: + w = sitk.ReadImage(path) + arr = sitk.GetArrayFromImage(w) + if arr.ndim == 4 and arr.shape[0] == 3: + arr = np.transpose(arr, (1, 2, 3, 0)) + return _field_transform(arr, ref_img) + try: + return sitk.ReadTransform(path) + except Exception: # noqa: binary ANTs transform file + pass + import ants + tx = ants.read_transform(path) + o = np.array(ants.apply_ants_transform_to_point(tx, [0.0, 0.0, 0.0]), + dtype=float) + M = np.empty((3, 3), dtype=float) + for i in range(3): + p = [0.0, 0.0, 0.0] + p[i] = 1.0 + M[:, i] = np.array(ants.apply_ants_transform_to_point(tx, p), + dtype=float) - o + t = sitk.AffineTransform(3) + t.SetMatrix(M.ravel()) + t.SetTranslation(o) + return t + + +def register_ants_variant_deformable(ct, t1c, variant): + """Run one stock antspyx ANTs preset (:data:`ANT_VARIANTS`) as a + deformable pipeline under the unified contract (fixed = T1c, moving = CT, + full volumes, no clip). The stage files ANTs writes (in numeric order = + the order it applies them) are composed and rasterised to the T1c-grid + displacement field of the contract, exactly like + :func:`register_antspyx_deformable`. A per-case cross-check (NCC of the + composite's own resampling against ANTs' ``warpedmovout`` over the MR + head) is recorded in the meta. + + OMP threads are left at the system default (uncapped) to match the + antspyx baseline row, which runs ``ants.registration`` the same way.""" + if variant not in ANT_VARIANTS: + raise KeyError(f"unknown ANTs variant {variant!r} " + f"(one of {list(ANT_VARIANTS)})") + tx_type, desc = ANT_VARIANTS[variant] + try: + import ants + except Exception as e: # noqa + raise PipelineUnavailable(f"antspyx: {e!r}") + import glob + import time + fx, mv = fixed_moving(ct, t1c) + af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)) + am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) + tmp = tempfile.mkdtemp(prefix=f"ants_{variant}_") + dt = 0.0 + stage_names = [] + crosscheck = None + try: + t0 = time.time() + r = ants.registration(fixed=af, moving=am, type_of_transform=tx_type, + outprefix=tmp + "/") + dt = time.time() - t0 + files = glob.glob(os.path.join(tmp, "*")) + stage_files = [f for f in files if f.endswith(".mat") or + ("Warp.nii.gz" in f and "Inverse" not in f)] + if not stage_files: + raise RuntimeError(f"{variant}: no stage transform files in {tmp}") + stage_files.sort(key=lambda p: int(re.match( + r"(\d+)", os.path.basename(p)).group(1))) + stage_names = [os.path.basename(p) for p in stage_files] + stages = [_ants_stage_to_sitk(p, fx) for p in stage_files] + T = stages[0] if len(stages) == 1 else sitk.CompositeTransform(stages) + T = _rasterize_to_tdf(T, fx) + # read-back cross-check against ANTs' own resampling (best effort) + try: + wo = os.path.join(tmp, "warpedmovout_check.nii.gz") + ants.image_write(r["warpedmovout"], wo) + own = sitk.GetArrayFromImage(sitk.ReadImage(wo)) + mine = sitk.GetArrayFromImage(sitk.Resample(mv, fx, T, sitk.sitkLinear)) + if own.shape == mine.shape: + from src.common import head_mask_from_image + head = head_mask_from_image( + sitk.GetArrayFromImage(fx).astype("float32")) + if head.sum() > 1000: + crosscheck = round(float(np.corrcoef( + mine[head].ravel(), own[head].ravel())[0, 1]), 4) + except Exception: # noqa + crosscheck = None + finally: + shutil.rmtree(tmp, ignore_errors=True) + ha = sitk.GetArrayFromImage(t1c) + gap = _skull_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": variant, + "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", + "metric": "ANTs stock preset metrics (Mattes MI / GC / CC per stage, " + "see transform description)", + "transform": f"Deformable (antspyx preset '{tx_type}'): {desc}", + "convention": "T1c->CT (verified vs ANTs warpedmovout)", + "auto_inverted": False, + "fixed": _modality_name(fx), "moving": _modality_name(mv), + "seconds_reg": round(dt, 1), + "stage_files": stage_names, + "own_warp_ncc_head": crosscheck, + "skull_ring_gap_mm": gap, + } + return T, meta + + +def register_ants_rigid_synonly_deformable(ct, t1c): + """The legacy two-call antspyx workflow (src/registration_pipelines_ + deformable.py's original `antspyx` row, kept for comparison): a `Rigid` + stage (Mattes 32, 500x500x200 iters, shrink 4x2x1, smoothing 3x2x0, COM + init) followed by a `SyNOnly` stage (SyN[0.2,3,0], Mattes 32, 100x100x0 + iters, shrink 4x2x1, smoothing 2x1x0) **initialised from the Rigid + stage's transform** (`initial_transform`). The SyNOnly run is asked for + its `Composite.h5` (write_composite_transform=True), which folds the + init in, and is rasterised to the T1c-grid displacement field of the + contract like the other antspyx pipelines. Cross-checked per case + against ANTs' own resampling (NCC in the meta). + + OMP threads are left at the system default (uncapped), matching the + antspyx baseline row.""" + try: + import ants + except Exception as e: # noqa + raise PipelineUnavailable(f"antspyx: {e!r}") + import time + fx, mv = fixed_moving(ct, t1c) + af = ants.from_sitk(sitk.Cast(fx, sitk.sitkFloat32)) + am = ants.from_sitk(sitk.Cast(mv, sitk.sitkFloat32)) + tmp = tempfile.mkdtemp(prefix="ants_rigid_synonly_") + dt = 0.0 + crosscheck = None + try: + 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), + outprefix=tmp + "/rigid_") + rig_path = r1["fwdtransforms"][0] if isinstance(r1["fwdtransforms"], + list) \ + else r1["fwdtransforms"] + r2 = ants.registration(fixed=af, moving=am, type_of_transform="SyNOnly", + initial_transform=rig_path, + reg_iterations=(100, 100, 100), + write_composite_transform=True, + outprefix=tmp + "/syn_") + dt = time.time() - t0 + # Composite.h5 = the -r init + the SyN stage (fixed -> moving chain) + T = sitk.ReadTransform(r2["fwdtransforms"]) + T = _rasterize_to_tdf(T, fx) + try: + wo = os.path.join(tmp, "warpedmovout_check.nii.gz") + ants.image_write(r2["warpedmovout"], wo) + own = sitk.GetArrayFromImage(sitk.ReadImage(wo)) + mine = sitk.GetArrayFromImage(sitk.Resample(mv, fx, T, sitk.sitkLinear)) + if own.shape == mine.shape: + from src.common import head_mask_from_image + head = head_mask_from_image( + sitk.GetArrayFromImage(fx).astype("float32")) + if head.sum() > 1000: + crosscheck = round(float(np.corrcoef( + mine[head].ravel(), own[head].ravel())[0, 1]), 4) + except Exception: # noqa + crosscheck = None + finally: + shutil.rmtree(tmp, ignore_errors=True) + ha = sitk.GetArrayFromImage(t1c) + gap = _skull_gap(_warp_ct(ct, t1c, T), ha) + meta = { + "pipeline": "ants_rigid_synonly", + "toolkit": f"ANTsPy {getattr(ants, '__version__', '?')}", + "metric": "Mattes MI (rigid 500x500x200, regular 20% sampling; " + "SyNOnly 100x100x0, full volume)", + "transform": ("Deformable (legacy two-call: Rigid[0.25] 500x500x200 " + "4x2x1 3x2x0 -> SyNOnly SyN[0.2,3,0] 100x100x0 4x2x1 " + "2x1x0, init = rigid stage; Composite.h5 read-back)"), + "convention": "T1c->CT (verified vs ANTs warpedmovout)", + "auto_inverted": False, + "fixed": _modality_name(fx), "moving": _modality_name(mv), + "seconds_reg": round(dt, 1), + "own_warp_ncc_head": crosscheck, + "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, @@ -2052,6 +2293,12 @@ def register_synthmorph_deformable(ct, t1c, gpu=None, threads=None, extent=256, return T, meta +def _ants_variant_register(name): + def _reg(ct, t1c): + return register_ants_variant_deformable(ct, t1c, variant=name) + return _reg + + REGISTERS_DEFORMABLE = { "antspyx": register_antspyx_deformable, "elastix": register_elastix_deformable, @@ -2060,7 +2307,12 @@ REGISTERS_DEFORMABLE = { "simpleitk": register_simpleitk_deformable, "synthmorph": register_synthmorph_deformable, } +REGISTERS_DEFORMABLE.update({n: _ants_variant_register(n) for n in ANT_VARIANTS}) +REGISTERS_DEFORMABLE["ants_rigid_synonly"] = register_ants_rigid_synonly_deformable ALL_PIPELINES_DEFORMABLE = list(REGISTERS_DEFORMABLE) +# every pipeline name the harness may dispatch / merge on (rigid + variants) +ALL_PIPELINE_NAMES = (ALL_PIPELINES + list(ANT_VARIANTS) + + ["ants_rigid_synonly"]) def availability(): @@ -2121,4 +2373,7 @@ def availability(): f"static config)") except Exception as e: # noqa out[name] = (False, repr(e)) + # the ANTs variants all run through antspyx (same availability as antspyx) + for name in list(ANT_VARIANTS) + ["ants_rigid_synonly"]: + out[name] = out.get("antspyx", (False, "antspyx not checked")) return out \ No newline at end of file