From aa7ac26f4c851b7b7cc611880d0edc468c78cdc0 Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Sat, 3 Oct 2026 13:25:43 +0800 Subject: [PATCH] Refactor code structure for improved readability and maintainability --- scripts/bench_registration_deformable.py | 439 --------- src/registration_pipelines_deformable.py | 1039 ---------------------- 2 files changed, 1478 deletions(-) delete mode 100644 scripts/bench_registration_deformable.py delete mode 100644 src/registration_pipelines_deformable.py diff --git a/scripts/bench_registration_deformable.py b/scripts/bench_registration_deformable.py deleted file mode 100644 index 420b288..0000000 --- a/scripts/bench_registration_deformable.py +++ /dev/null @@ -1,439 +0,0 @@ -"""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/src/registration_pipelines_deformable.py b/src/registration_pipelines_deformable.py deleted file mode 100644 index 60b93d6..0000000 --- a/src/registration_pipelines_deformable.py +++ /dev/null @@ -1,1039 +0,0 @@ -"""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