.gitignore

@@ -7,5 +7,3 @@ runs_bench/
 runs_nnu/
 runs_monai/
 __pycache__/
-
-elastix/

REGISTRATION_BENCHMARK.md

@@ -32,7 +32,7 @@ handling — not on who finds the ~30–60 mm gross offset first.
 
 | # | pipeline | toolkit | metric | transform | status in this env |
 |---|----------|---------|--------|-----------|--------------------|
-| 1 | `elastix`   | ITK-Elastix 0.25.4 (pip `itk` + `itk-elastix`, elastix 5.3.1) | Mattes MI (10 bins) | `EulerTransform` rigid (5.x replacement for `Rigid3DTransform`), ASGD, 3 resolutions | ✅ ran (main env, in-process python API — no CLI) |
+| 1 | `elastix`   | ITK-Elastix | Mattes MI | `Rigid3DTransform` (ASGD) | **not installed** (conda-forge unreachable) |
 | 2 | `simpleitk` | SimpleITK 2.5.6 | Mattes MI (32 bins) | `Euler3D` rigid, pyramid `[4,2,1]` | ✅ ran |
 | 3 | `antspyx`   | ANTsPy / ANTs 0.6.3 | Mattes MI (multi-res) | ANTs `Rigid` | ✅ ran (isolated py3.12 env) |
 | 4 | `deepreg`   | Oxford DeepReg (CycleGAN) | CycleGAN + reg. loss (pretrained) | 3D displacement field | **not installed** (see below) |
@@ -71,84 +71,56 @@ Production QA, crop-independent:
 - **`tv_nvox`** — warped GTV volume.
 - **`accepted`** — `tv_nvox > 0 and tv_in_ct_frac ≥ 0.5 and aligned` (the
   production acceptance rule).
-- **`mi_ct_mr`** — mean Mattes MI (nats) between the **CT inverse-transformed
-  into MR coordinates** (the warped CT on the T1c grid) and the **original MR**,
-  restricted to the MR head and each image percentile-binned independently
-  (SimpleITK Mattes-MI style). This is the registration's *achieved* metric read
-  off as a standalone alignment score: higher = the warped-CT and MR intensities
-  are better in-sync = better aligned; it collapses toward 0 when the CT is
-  rotated/translated away from the MR.
 - **`reg_wall_s`** — wall-clock registration time (excludes I/O + scoring).
 
 ## Results (10 cases)
 
-| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
+| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT |
-|---|---|---|---|---|---|---|---|---|
+|---|---|---|---|---|---|---|---|
-| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.6 | 0.93 | **0.643 / 0.630** |
+| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.3 | 0.93 |
-| **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 14.2 / 17.3 | 24.8 / 22.4 | 12.7 | 0.93 | **0.600 / 0.523** |
+| **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 14.3 / 15.4 | 18.2 / 21.7 | 13.4 | 0.93 |
-| **elastix** | ITK-Elastix 0.25.4 | 9 | **5/10 (0.50)** | 15.3 / 16.1 | 9.9 / 13.8 | 3.2 | 0.93 | **0.562 / 0.549** |
+| elastix | — | — | — | — | — | — | — |
-| deepreg | — | — | — | — | — | — | — | — |
+| deepreg | — | — | — | — | — | — | — |
 
 Per-case (mid / std, `OK` = accepted):
 
-| case | simpleitk | antspyx | elastix |
+| case | simpleitk | antspyx |
-|---|---|---|---|
+|---|---|---|
-| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT |
+| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT |
-| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 26.3 / 29.2 | 26.4 / 30.3 |
+| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 26.3 / 29.2 |
-| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.8 / 25.0 | 21.1 / 23.3 |
+| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.8 / 25.0 |
-| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.4 / 18.2 | 10.6 / 8.2 **OK** |
+| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.4 / 18.2 |
-| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 12.5 / 12.8 | 11.0 / 6.1 **OK** |
+| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 12.5 / 12.8 |
-| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.8 / 11.7 **OK** | 15.3 / 14.0 |
+| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.8 / 11.7 **OK** |
-| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 6.1 / 37.2 | 11.4 / 3.2 **OK** |
+| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 6.1 / 37.2 |
-| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** |
+| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 |
-| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.6 / 14.7 | 12.0 / 6.8 **OK** |
+| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.6 / 14.7 |
-| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 18.2 / 33.6 | 18.9 / 22.7 |
+| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 18.2 / 33.6 |
 
 ### Interpretation
 
-- **Translation is recovered well by all three** (mid gaps mostly 6–20 mm,
+- **Translation is recovered well by both** (mid gaps mostly 6–20 mm, within the
-  within the `≤ 20 mm` band). `tv-in-CT` is identical (~0.93): the warped GTV
+  `≤ 20 mm` band). `tv-in-CT` is identical (0.93): the warped GTV lands inside
-  lands inside the CT head for every pipeline, so the target is in the right
+  the CT head for both, so the target is in the right place.
-  place.
 - **The binding failure is rotation** — `std > 12 mm` on the repositioned cases
   (the ~30° inter-session tilt produces MI local minima that a single centered
   start can lock onto). This is exactly why production `02_build_m6_dataset.py`
   uses a **multi-start tilt grid**.
-- **Elastix is the most accurate of the three standard single-init configs**
+- **SimpleITK is the more robust of the two** (4 vs 1 accepted, lower median
-  (5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 24.8) and
+  std 15.4 vs 18.2): its RSGD + physical-step-scales + 32-bin Mattes-MI
-  the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s
+  landscape handles the tilt better than ANTs' Rigid. Note `m6_26TYQI3C`: ANTs
-  SimpleITK / 12.7 s ANTs): ASGD's adaptive step sizes climb out of small
+  nailed the position (mid 6.1) but missed the orientation (std 37.2) — a
-  misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223`
+  textbook tilted local minimum.
-  (std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on
+- **Speed**: SimpleITK ~8.3 s/case, ANTs ~13.4 s/case (ANTs includes the
-  (std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both
+  pre-shift resample + its affine-MI init stages). Both are CPU-bound and ran
-  others accept — different optimizers settle in different local minima on the
+  4 cases in parallel.
-  hard tilted cases.
-- **ANTs' Rigid is the weakest** (1/10 accepted): it frequently locks onto the
-  tilted local minimum. Note `m6_26TYQI3C`: ANTs nailed the position (mid 6.1)
-  but missed the orientation (std 37.2) — a textbook tilted local minimum.
-- **Speed** (median wall, 4 cases in parallel, 24 threads): Elastix ~3.2 s
-  (warm workers; a worker's first case pays a ~10 s warm-up), SimpleITK ~8.6 s,
-  ANTs ~12.7 s (includes the pre-shift resample + its affine-MI init stages).
-  All CPU-bound.
-- **Mattes MI (warped-CT ~ MR) is an independent, intensity-based confirmation**:
-  medians 0.643 (SimpleITK) / 0.600 (ANTs) / 0.562 (Elastix), and it diverges
-  sharply on the tilted cases — ANTs `m6_26TYQI3C` drops to 0.094 nats
-  (rotation ~37 mm off) while SimpleITK holds 0.670. A high MI means the CT
-  skull/tissue intensities fall on the matching MR structures after the inverse
-  transform.
-- **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
-  per-case numbers carry run-to-run variance (this table is one draw);
-  SimpleITK's MI (and registration) is deterministic at fixed thread count
-  (RSGD wall time is thread-sensitive, convergence is not); Elastix's ASGD is
-  deterministic given its fixed seed.
 
 ## Reproducing
 
 ```bash
-# main env — SimpleITK + Elastix (itk-elastix is pip-installed here;
+# main env — SimpleITK baseline (elastix/deepreg report unavailable here)
-# antspyx/deepreg report unavailable and are recorded as such)
 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal
-python scripts/bench_registration.py --n-cases 10 \
+python scripts/bench_registration.py --pipelines simpleitk --n-cases 10 \
     --workers 4 --threads 24 --simpleitk-mode benchmark
 
 # isolated py3.12 env with antspyx — ANTs (+ optional SimpleITK cross-check)
@@ -162,42 +134,20 @@ python scripts/bench_registration.py --merge
 Other flags: `--sample random --seed 0` (reproducible random 10),
 `--keys k1,k2,...` (explicit cases), `--n-cases N`, `--workers W`,
 `--threads T`, `--simpleitk-mode {benchmark,production}`,
-`--deepreg-ckpt <.ckpt>` (DeepReg), `--out-dir <dir>`,
+`--deepreg-ckpt <.ckpt>` (DeepReg), `--out-dir <dir>`.
-`--qa-images` / `--qa-out <dir>` (visual QA, below).
 
-## Visual QA images (checkerboard + color overlay)
+## Enabling the two blocked pipelines
 
-`--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to
+**ITK-Elastix** — not on PyPI; needs conda-forge (which this host cannot reach).
-`<out-dir>/qa/<pipeline>__<case>.png` for eyeballing the alignment.  Each image
+Once conda-forge is available:
-is a 3×2 grid: **rows** are the three orthogonal mid-slices through the head
-center (**mid-axial, mid-coronal, mid-sagittal**), and **columns** are:
-
-- **left — checkerboard blend:** alternating blocks of the CT inverse-transformed
-  into MR coordinates and the original MR; a clean registration gives a seamless
-  mosaic, a misregistration shows hard edges at the block boundaries.
-- **right — additive color overlay:** CT in red, MR in green. Structures that
-  coincide blend to yellow; rotation/translation error appears as red/green
-  fringing at the skull and brain edges.
-
-The CT is shown in a **brain window (L40/W80 = 0–80 HU)** so the brain
-parenchyma is visible and comparable to the MR (the skull, ~800+ HU, clips to
-white); the MR uses a 2–98 % head window. The coronal and sagittal views are
-oriented superior-up to match the axial view.
-
-The title carries the pipeline, case key, Mattes MI, and the skull-ring gap
-(mid/std), so a glance tells you both the achieved-metric score and whether the
-visuals agree. Generate for the available pipelines with:
 
 ```bash
-/home/xfr/.conda/envs/longitudinal/bin/python scripts/bench_registration.py \
+conda install -c conda-forge elastix        # provides the elastix/transformix CLI
-    --pipelines simpleitk elastix --n-cases 10 --qa-images
+python scripts/bench_registration.py --pipelines elastix --n-cases 10
-/home/xfr/.conda/envs/ants312/bin/python scripts/bench_registration.py \
-    --pipelines antspyx --n-cases 10 --qa-images
 ```
 
-(9 images each; the incomplete-CT case is skipped.)
+The adapter registers `fixed=T1c, moving=CT` with `Rigid3DTransform` + Mattes MI
-
+on the head box and reads the result transform.
-## Enabling the blocked pipeline
 
 **DeepReg (Oxford)** — the PyPI `deepreg` (0.1.2) is an unrelated stub that pins
 `numpy==1.18.5`; do not install it. DeepReg is not a one-shot rigid tool: it is a
@@ -219,32 +169,9 @@ it to a SimpleITK `DisplacementFieldTransform`, and scores it like the others.
   start). Production uses `--simpleitk-mode production`'s multi-start tilt grid,
   which raises SimpleITK's acceptance at ~20× the cost; that is out of scope for
   the head-to-head but is the config the dataset actually ships with.
-- `deepreg` could not be executed here, so its row is "not installed", not
+- `elastix` / `deepreg` could not be executed here, so their rows are
-  measured; its code path is implemented and unit-verified for import/convention.
+  "not installed", not measured. Their code paths are implemented and
+  unit-verified for import/convention; run them where the toolkits are present.
 - `fix_orientation` (elastix/ants/deepreg) chooses between a transform and its
   inverse by skull-ring concentricity — a correctness gate that makes the result
-  independent of each toolkit's transform-direction convention. In practice it
+  independent of each toolkit's transform-direction convention.
-  never fired (all `auto_inverted: False`): elastix/ants returned the T1c→CT
-  direction as expected.
-- **Elastix specifics** (ITK-Elastix 0.25.4, elastix 5.3.1, in-process python
-  API — no CLI binary involved):
-  - Elastix 5.x renamed components vs the classic parameter files: the rigid
-    3D transform is `EulerTransform` (6-DOF, `Rigid3DTransform` is gone from
-    the component set), the metric key is `Metric` (not `ImageMetric`), the
-    random sampler is `Random` (not `RandomSampler`). The adapter uses the 5.x
-    names and is not compatible with elastix ≤ 4.x parameter files.
-  - The initial transform (the shared `centroid_offset` start) is passed as an
-    `InitialTransformParameterObject` (a `TranslationTransform` with
-    `NumberOfParameters 3`); the C++ filter only wires it up when an
-    **output directory is set** (`SetOutputDirectory`), so the adapter uses a
-    scratch dir under `data/m6reg/elastix/itkelastix/`.
-  - `WriteResultImage` must be the string `"false"` (Python `False`); the
-    default (absent) makes the filter raise "Could not read result image"
-    after a *successful* registration.
-  - The result transform is recovered by point-mapping
-    (`GetCombinationTransform().TransformPoint`) into an affine — independent
-    of elastix's Euler-angle parameterization — and wrapped in a SimpleITK
-    `AffineTransform(3)`.
-  - Registration wall time is ~1.8 s/case at full budget (ASGD 500/500/250
-    over 3 resolutions); each worker's first case pays a ~10 s warm-up, so the
-    median wall is 3.2 s and the mean is inflated by warm-up.

scripts/bench_registration.py

@@ -17,16 +17,15 @@ geometric QA used in production (scripts/02_build_m6_dataset.py):
   - accepted           : production rule (nvox>0 and tv_in_ct and concentric).
   - reg_wall_s         : wall-clock registration time (excludes I/O + scoring).
 
-Because two of the four toolkits are not pip-installable into the main
+Because three of the four toolkits are not pip-installable into the main
-`longitudinal` env (antspyx has no cp314 wheel; PyPI deepreg is a broken
+`longitudinal` env (elastix needs conda-forge; antspyx has no cp314 wheel;
-stub), the benchmark is env-portable: run it under whichever conda env has a
+PyPI deepreg is a broken stub), the benchmark is env-portable: run it under
-given toolkit installed, and it records only the pipelines that are
+whichever conda env has a given toolkit installed, and it records only the
-importable there.  elastix runs in the main env through the ITK-Elastix
+pipelines that are importable there.  The `--merge` mode then combines the
-python API (pip ``itk`` + ``itk-elastix``).  The `--merge` mode then combines
+per-pipeline JSON files (accumulated across envs) into one summary + plot.
-the per-pipeline JSON files (accumulated across envs) into one summary + plot.
 
 Usage (repo root):
-  # main env (SimpleITK + elastix; deepreg reports as unavailable)
+  # main env (SimpleITK baseline; elastix/deepreg report as unavailable)
   source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal
   python scripts/bench_registration.py --pipelines simpleitk elastix deepreg
 
@@ -53,34 +52,10 @@ import numpy as np
 import SimpleITK as sitk
 from scipy import ndimage
 
-from src.common import DATA, d, read_nii_img, head_mask_from_image
+from src.common import DATA, d, read_nii_img
 from src import registration_pipelines as RP
 
 
-def _mi_masked(fixed, moving, mask, bins=32):
-    """Mattes MI (nats) between two co-registered volumes restricted to ``mask``,
-    each image percentile-binned independently (SimpleITK Mattes-MI style).
-    Higher = the warped-CT and MR intensities are better in-sync = better
-    alignment.  ``fixed`` is the CT inverse-transformed into MR coordinates."""
-    f = np.asarray(fixed)[mask]
-    m = np.asarray(moving)[mask]
-    if f.size < 1000:
-        return None
-
-    def q(a):
-        lo, hi = np.percentile(a, (1, 99))
-        return np.clip(np.round((a - lo) * (bins - 1) / max(hi - lo, 1e-9)).astype(int),
-                       0, bins - 1)
-
-    A, B = q(f), q(m)
-    h = np.bincount(A.ravel() * bins + B.ravel(), minlength=bins * bins)
-    p = h.reshape(bins, bins) / max(h.sum(), 1)
-    px, py = p.sum(axis=1, keepdims=True), p.sum(axis=0, keepdims=True)
-    with np.errstate(divide="ignore", invalid="ignore"):
-        ll = np.log(p / (px * py))
-    return float(np.where((p > 0) & (px > 0) & (py > 0), p * ll, 0.0).sum())
-
-
 def _m6():
     return RP._load_m6()
 
@@ -131,10 +106,6 @@ def evaluate(ct, t1c, tv, T):
     warped_ct = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
     warped_tv = sitk.GetArrayFromImage(sitk.Resample(tv, t1c, T, sitk.sitkNearestNeighbor, 0.0))
     gap = m6._skull_ring_gap(warped_ct, ha)
-    # achieved-metric alignment score: Mattes MI between the CT inverse-transformed
-    # into MR coordinates (warped_ct) and the original MR (ha), over the MR head.
-    ha_f = ha.astype("float32", copy=False)
-    mi = _mi_masked(warped_ct, ha_f, head_mask_from_image(ha_f))
     posw = warped_ct[warped_ct > 0]
     tissue_w = (warped_ct > max(2.0, float(np.percentile(posw, 2)))) if posw.size \
         else np.zeros_like(warped_ct, bool)
@@ -148,7 +119,6 @@ def evaluate(ct, t1c, tv, T):
     accepted = nvox > 0 and tv_in_ct_frac >= 0.5 and ok_gap
     return {
         "skull_ring_gap_mm": gap,
-        "mi_ct_mr": round(mi, 4) if mi is not None else None,
         "tv_nvox": nvox,
         "tv_in_ct_frac": round(tv_in_ct_frac, 3),
         "alignment": "concentric" if ok_gap else "not_concentric",
@@ -156,64 +126,6 @@ def evaluate(ct, t1c, tv, T):
     }
 
 
-def make_qa_png(warped_ct, mr, out_png, title, block=24, ct_brain_window=(0.0, 80.0)):
-    """checkreg-style visual QA of (inverse-warped CT) vs (original MR):
-    rows = the three orthogonal mid-slices through the head center
-    (mid-axial, mid-coronal, mid-sagittal); columns = left checkerboard blend,
-    right additive color overlay (CT=red, MR=green -> yellow where aligned,
-    red/green fringes where not).  The CT is shown in a brain window
-    (default L40/W80 = 0..80 HU).  Writes ``out_png``."""
-    import matplotlib
-    matplotlib.use("Agg")
-    import matplotlib.pyplot as plt
-    mr = np.asarray(mr, dtype="float32")
-    wc = np.asarray(warped_ct, dtype="float32")
-    pos = mr > 0
-    thr = np.percentile(mr[pos], 2) if pos.any() else 0.0
-    head = mr > max(thr, 1e-6)
-    cen = np.argwhere(head).mean(0) if head.any() else np.array(mr.shape) / 2.0
-    ia, ic, isg = (int(round(cen[0])), int(round(cen[1])), int(round(cen[2])))
-    # orthogonal mid-slices (arrays are [z, y, x]); CT and MR sliced identically
-    # so the two stay registered in each view. coronal/sagittal flip the z-axis
-    # (superior up) to match the axial view.
-    def _flip0(a):
-        return np.ascontiguousarray(a[::-1, :])
-
-    views = [
-        ("Mid-axial", np.ascontiguousarray(mr[ia]), np.ascontiguousarray(wc[ia])),
-        ("Mid-coronal", _flip0(mr[:, ic, :]), _flip0(wc[:, ic, :])),
-        ("Mid-sagittal", _flip0(mr[:, :, isg]), _flip0(wc[:, :, isg])),
-    ]
-    mr_lo, mr_hi = np.percentile(mr[head], (2, 98)) if head.any() else (0.0, 1.0)
-    ct_lo, ct_hi = ct_brain_window  # CT brain window (0..80 HU default)
-
-    def w(a, lo, hi):
-        return np.clip((a - lo) / max(hi - lo, 1e-9), 0.0, 1.0)
-
-    n = len(views)
-    fig, axes = plt.subplots(nrows=n, ncols=2, figsize=(7.6, 2.15 * n))
-    for r, (vname, mrs_raw, wcs_raw) in enumerate(views):
-        mrs, wcs = w(mrs_raw, mr_lo, mr_hi), w(wcs_raw, ct_lo, ct_hi)
-        yy, xx = np.indices(mrs.shape)
-        cb = mrs.copy()
-        cb[((yy // block) + (xx // block)) % 2 == 0] = wcs[((yy // block) + (xx // block)) % 2 == 0]
-        axes[r, 0].imshow(cb, cmap="gray", vmin=0, vmax=1); axes[r, 0].axis("off")
-        axes[r, 0].text(0.01, 0.02, vname, transform=axes[r, 0].transAxes, color="white",
-                        bbox=dict(facecolor="black", alpha=0.55, pad=1.2), fontsize=8)
-        ov = np.zeros((mrs.shape[0], mrs.shape[1], 3))
-        ov[..., 0] = wcs; ov[..., 1] = mrs
-        ov[..., 2] = 0.28 * np.minimum(wcs, mrs)
-        axes[r, 1].imshow(ov); axes[r, 1].axis("off")
-        axes[r, 1].text(0.01, 0.02, vname, transform=axes[r, 1].transAxes, color="white",
-                        bbox=dict(facecolor="black", alpha=0.55, pad=1.2), fontsize=8)
-    axes[0, 0].set_title("Checkerboard blend (CT | MR)", fontsize=9)
-    axes[0, 1].set_title("Color overlay (CT=red, MR=green)", fontsize=9)
-    fig.suptitle(title, fontsize=10)
-    fig.tight_layout(rect=(0, 0, 1, 0.97))
-    fig.savefig(out_png, dpi=110)
-    plt.close(fig)
-
-
 def _call(pipeline, ct, t1c, tv, args):
     if pipeline == "simpleitk":
         return RP.register_simpleitk(ct, t1c, threads=args.threads, mode=args.simpleitk_mode)
@@ -246,18 +158,6 @@ def _run_case(pipeline, case, args):
     out["reg_wall_s"] = round(time.time() - t0, 1)
     out.update(meta)
     out.update(evaluate(ct, t1c, tv, T))
-    if getattr(args, "qa_images", False):
-        g = out.get("skull_ring_gap_mm")
-        title = (f"{pipeline} | {key} | MI={out.get('mi_ct_mr')} | "
-                 f"gap mid {g['mid_mm']}mm std {g['std_mm']}mm" if g
-                 else f"{pipeline} | {key}")
-        qp = os.path.join(args.qa_out, f"{pipeline}__{key}.png")
-        try:
-            wc = sitk.GetArrayFromImage(sitk.Resample(ct, t1c, T, sitk.sitkLinear))
-            make_qa_png(wc, sitk.GetArrayFromImage(t1c), qp, title)
-            out["qa_png"] = qp
-        except Exception as e:  # noqa
-            out["qa_png_error"] = repr(e)
     return out
 
 
@@ -284,7 +184,6 @@ def run_pipeline(pipeline, cases, args):
                 print(f"   {c['key']}: accepted={r['accepted']} "
                       f"gap={g['mid_mm'] if g else 'None'}mm "
                       f"std={g['std_mm'] if g else 'None'}mm "
-                      f"mi={r.get('mi_ct_mr')} "
                       f"tv_in_ct={r.get('tv_in_ct_frac')} {reg}", flush=True)
             except Exception as e:  # noqa
                 errors.append({"key": c["key"], "error": repr(e)})
@@ -298,7 +197,6 @@ def run_pipeline(pipeline, cases, args):
     gaps_std = [r["skull_ring_gap_mm"]["std_mm"] for r in registered if r.get("skull_ring_gap_mm")]
     regt = [r["reg_wall_s"] for r in registered]
     tvfr = [r["tv_in_ct_frac"] for r in registered if r.get("tv_in_ct_frac") is not None]
-    mifr = [r["mi_ct_mr"] for r in registered if r.get("mi_ct_mr") is not None]
     first = (registered or results)[0] if (registered or results) else {}
     row = {
         "pipeline": pipeline, "available": True, "toolkit": first.get("toolkit"),
@@ -314,8 +212,6 @@ def run_pipeline(pipeline, cases, args):
         "skull_ring_gap_std_mm": {"mean": round(float(np.mean(gaps_std)), 1),
                                   "median": round(float(np.median(gaps_std)), 1)} if gaps_std else None,
         "tv_in_ct_frac": {"mean": round(float(np.mean(tvfr)), 3)} if tvfr else None,
-        "mi_ct_mr": {"mean": round(float(np.mean(mifr)), 4),
-                     "median": round(float(np.median(mifr)), 4)} if mifr else None,
         "cases": results, "errors": errors,
         "timestamp": datetime.now().isoformat(timespec="seconds"),
     }
@@ -342,24 +238,23 @@ def print_table(rows):
     order = [p for p in RP.ALL_PIPELINES if any(r["pipeline"] == p for r in rows)]
     print("\n" + "=" * 100)
     hdr = (f"{'pipeline':10s} {'avail':6s} {'toolkit':22s} {'n':>3s} {'ok':>3s} {'err':>4s} "
-            f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s} {'MI':>7s}")
+           f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s}")
     print(hdr)
     print("-" * 100)
     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'):6s} {'-':22s} {r['n_cases']:>3d} {0:>3d} {'-':>4s} "
-                  f"{'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s}   [{r['reason']}]")
+                  f"{'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s} {'-':>7s}   [{r['reason']}]")
             continue
         gm = (r.get("skull_ring_gap_mid_mm") or {}).get("median")
         gs = (r.get("skull_ring_gap_std_mm") or {}).get("median")
         rt = (r.get("reg_wall_s") or {}).get("median")
         tv = (r.get("tv_in_ct_frac") or {}).get("mean")
-        mi = (r.get("mi_ct_mr") or {}).get("median")
         print(f"{p:10s} {'yes':6s} {(r.get('toolkit') or '')[:22]:22s} {r['n_cases']:>3d} "
               f"{r.get('n_registered', 0):>3d} {r['n_error']:>4d} {r['accepted_rate']:>7.2f} "
               f"{_fmt(gm, '{:.1f}'):>7s} {_fmt(gs, '{:.1f}'):>7s} {_fmt(rt, '{:.1f}'):>7s} "
-              f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
+              f"{_fmt(tv, '{:.2f}'):>7s}")
     print("=" * 100)
 
 
@@ -380,15 +275,13 @@ def write_plot(rows, out_png):
     gap_std = [ok_rows[p].get("skull_ring_gap_std_mm", {}).get("median") or 0 for p in labels]
     regt = [ok_rows[p].get("reg_wall_s", {}).get("median") or 0 for p in labels]
     acc = [ok_rows[p].get("accepted_rate") or 0 for p in labels]
-    mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels]
+    fig, axes = plt.subplots(1, 4, figsize=(18, 4))
-    fig, axes = plt.subplots(1, 5, figsize=(20, 4))
+    colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD"]
-    colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD", "#E08E45"]
     for ax, (title, vals, fmt, ylab, inv) in zip(axes, [
         ("Skull-ring gap mid (mm, median)", gap_mid, "{:.1f}", "mm", True),
         ("Skull-ring gap std (mm, median)", gap_std, "{:.1f}", "mm", True),
         ("Registration time (s, median)", regt, "{:.0f}", "s", True),
-        ("Acceptance rate", acc, "{:.2f}", "frac", False),
+        ("Acceptance rate", acc, "{:.2f}", "frac", False)]):
-        ("Mean Mattes MI: warped-CT ~ MR (nats, median)", mi, "{:.3f}", "nats", False)]):
         bars = ax.bar(range(len(vals)), vals, color=colors[:len(vals)])
         for b, v in zip(bars, vals):
             ax.annotate(fmt.format(v), (b.get_x() + b.get_width() / 2, b.get_height()),
@@ -432,12 +325,12 @@ def merge(args):
     print("summary:", sf)
     write_plot(rows, os.path.join(out_dir, "bench_registration.png"))
     # markdown table for the report
-    md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT | MI (nats) |",
+    md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT |",
-          "|---|---|---|---|---|---|---|---|---|---|---|---|"]
+          "|---|---|---|---|---|---|---|---|---|---|---|"]
     for p in [x["pipeline"] for x in rows]:
         r = next(x for x in rows if x["pipeline"] == p)
         if not r.get("available"):
-            md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | - | - | - | ({r['reason']}) |")
+            md.append(f"| {p} | NO | - | {r['n_cases']} | 0 | - | - | - | - | - | - | ({r['reason']}) |")
             continue
         md.append(
             f"| {p} | yes | {(r.get('toolkit') or '')} | {r['n_cases']} | {r.get('n_registered', 0)} | "
@@ -445,8 +338,7 @@ def merge(args):
             f"{(r.get('skull_ring_gap_mid_mm') or {}).get('median')} | "
             f"{(r.get('skull_ring_gap_std_mm') or {}).get('median')} | "
             f"{(r.get('reg_wall_s') or {}).get('median')} | "
-            f"{(r.get('tv_in_ct_frac') or {}).get('mean')} | "
+            f"{(r.get('tv_in_ct_frac') or {}).get('mean')} |")
-            f"{(r.get('mi_ct_mr') or {}).get('median')} |")
     mf = os.path.join(out_dir, "summary.md")
     with open(mf, "w") as f:
         f.write("\n".join(md) + "\n")
@@ -471,11 +363,6 @@ def parse_args():
     ap.add_argument("--elastix-dir", default=None, help="elastix scratch dir (default data/m6reg/elastix)")
     ap.add_argument("--deepreg-dir", default=None, help="deepreg scratch dir")
     ap.add_argument("--deepreg-ckpt", default=None, help="path to a DeepReg 3D pretrained checkpoint")
-    ap.add_argument("--qa-images", action="store_true",
-                    help="also write a checkreg-style PNG (checkerboard blend + color "
-                         "overlay of the inverse-warped CT vs the original MR) per case")
-    ap.add_argument("--qa-out", default=None,
-                    help="QA image dir (default <out-dir>/qa)")
     ap.add_argument("--merge", action="store_true", help="combine existing per-pipeline JSONs")
     args = ap.parse_args()
     return args
@@ -488,9 +375,6 @@ def main():
         return
     args.elastix_dir = args.elastix_dir or d("data/m6reg/elastix")
     args.deepreg_dir = args.deepreg_dir or d("data/m6reg/deepreg")
-    if args.qa_images:
-        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)
     if unknown:

src/registration_pipelines.py

@@ -26,8 +26,7 @@ Pipelines
 - ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a
   CT-space head box, multi-start tilt grid.  Reuses the *proven* production
   ``register_ct_to_t1c`` so the baseline is bit-identical to 02_build_m6_dataset.
-- ``register_elastix``   : ITK-Elastix Python API (pip ``itk`` + ``itk-elastix``,
+- ``register_elastix``   : ITK-Elastix CLI (Rigid3DTransform, Mattes MI).
-  elastix 5.x: ``EulerTransform`` rigid + Mattes MI, ASGD).
 - ``register_antspyx``   : ANTs / ANTsPy ``ants.registration`` (Rigid).
 - ``register_deepreg``   : Oxford DeepReg (CycleGAN deformation-field
   registration; needs a pretrained checkpoint, see module docstring).
@@ -50,6 +49,7 @@ SimpleITK path.
 """
 import os
 import json
+import shutil
 import subprocess
 import importlib.util
 
@@ -242,74 +242,57 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
 # 1. ITK-Elastix
 # --------------------------------------------------------------------------- #
 
-def _elastix_api():
+def _elastix_param_file(workdir, fixed, moving, fixed_mask, moving_mask):
-    """Import the ITK-Elastix python API, or raise PipelineUnavailable."""
+    """Rigid3DTransform + Mattes MI parameter file."""
-    try:
+    lines = [
-        import itk
+        "(ImageDimension 3)",
-    except Exception as e:  # noqa
+        "(FixedInternalImagePixelType float)",
-        raise PipelineUnavailable(
+        "(MovingInternalImagePixelType float)",
-            f"itk not importable ({e!r}). Install with `pip install itk itk-elastix` "
+        "(Transform \"Rigid3DTransform\")",
-            "(ITK-Elastix >= 0.25 exposes the elastix components in the `itk` "
+        "(TransformBase \"ElastixTransformBase\")",
-            "python namespace; no CLI / conda-forge needed).")
+        "(FixedCenteredTransform 1)",
-    try:
+        "(TransformParameterReadoutFormat 1)",
-        _ = itk.ElastixRegistrationMethod  # loads the elastix module + components
+        "(Optimizer \"AdaptiveStochasticGradientDescent\")",
-    except Exception as e:  # noqa
+        "(MaximumNumberOfIterations 6000)",
-        raise PipelineUnavailable(
+        "(StoppingValue 0.0005)",
-            f"ITK-Elastix components not registered in the `itk` python namespace "
+        "(StoppingIterationCount 2500)",
-            f"({type(e).__name__}: {e}). `pip install itk-elastix` (with a "
+        "(ImageMetric \"MattesMutualInformation\")",
-            "matching `itk`) provides `itk.ElastixRegistrationMethod`.")
+        "(NumberOfHistogramBins 10)",
-    return itk
+        "(ImageSampler \"RandomSampler\")",
-
+        "(ImageSamplerPatchSize 25 25 25)",
-
+        "(FixedSubstresamplerClass \"RandomSubstresampler\")",
-def _sitk_to_itk(img, itk):
+        "(FixedSubstresamplerScale 2)",
-    """SimpleITK image -> itk python image (float32, geometry preserved)."""
+        "(ResampleInterpolator \"FinalBSplineInterpolator\")",
-    out = itk.GetImageFromArray(sitk.GetArrayFromImage(sitk.Cast(img, sitk.sitkFloat32)))
+        "(FinalBSplineInterpolationOrder 1)",
-    out.SetSpacing([float(s) for s in img.GetSpacing()])
+        "(ResultStationaryImage \"result.mha\")",
-    out.SetOrigin([float(o) for o in img.GetOrigin()])
+        "(WriteResultImage 0)",
-    out.SetDirection(itk.GetMatrixFromArray(
+        "(CompressResultImages 1)",
-        np.array(img.GetDirection(), dtype=float).reshape(3, 3)))
+        f"(FixedImage {fixed})",
-    return out
+        f"(FixedImageMask {fixed_mask})",
-
+        f"(MovingImage {moving})",
-
+        f"(MovingImageMask {moving_mask})",
-def _itk_transform_to_sitk_rigid(itk_tf):
+    ]
-    """Recover the affine (R, t) of an ITK transform by point mapping and build
+    p = os.path.join(workdir, "elastix_params.txt")
-    a SimpleITK transform from it.  The mapping is p -> R p + t in world
+    with open(p, "w") as f:
-    coordinates, so a plain AffineTransform represents the rigid result exactly
+        f.write("\n".join(lines) + "\n")
-    (independent of the toolkit's internal angle parameterization)."""
+    return p
-    o = np.array(itk_tf.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(itk_tf.TransformPoint(p), dtype=float) - o
-    t = sitk.AffineTransform(3)
-    t.SetMatrix([float(v) for v in R.ravel()])
-    t.SetTranslation([float(v) for v in o])
-    return t
 
 
 def register_elastix(ct, t1c, workdir=None, threads=8):
-    """ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI.
+    """ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI."""
-
+    if shutil.which("elastix") is None:
-    Runs in-process through the ITK-Elastix python API (pip ``itk`` +
+        raise PipelineUnavailable(
-    ``itk-elastix``, elastix 5.x).  Rigid 6-DOF via ``EulerTransform`` (the
+            "elastix/transformix binaries not found. Install with "
-    5.x replacement for ``Rigid3DTransform``) + Mattes MI
+            "`conda install -c conda-forge elastix` (conda-forge must be "
-    (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common
+            "reachable) or build from https://elastix.org. Not on PyPI.")
-    CT-space head box, from the *shared* centroid-aligned start
+    workdir = workdir or d("data/m6reg/elastix")
-    (``centroid_offset`` as an initial TranslationTransform, same init the
+    workdir = os.path.join(workdir, f"pid{os.getpid()}")  # per-worker scratch
-    SimpleITK/ANTs benchmark modes use).  Elastix 5.x component names: sampler
+    os.makedirs(workdir, exist_ok=True)
-    ``Random`` (not ``RandomSampler``), metric key ``Metric`` (not
-    ``ImageMetric``).  An output directory is required so the initial
-    transform parameter object is honored (the C++ filter only wires it up
-    when ``SetOutputDirectory`` is set), and ``WriteResultImage`` must be the
-    string ``"false"`` (a Python bool False) or the filter raises "Could not
-    read result image" after a successful registration."""
-    itk = _elastix_api()
     lo, hi = _ct_head_box_mm(ct)
     if lo is None:
         raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)")
-    # Resample both into a common CT-space head box (mirrors the other pipelines)
+    # Resample both into a common CT-space head box (mirrors the production grid
-    # so MI is computed on head tissue, not neck/air.
+    # choice) so MI is computed on head tissue, not neck/air.
     tsp = np.array(t1c.GetSpacing())
     shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3))
     ref = sitk.GetImageFromArray(np.zeros(shape, np.float32))
@@ -318,69 +301,36 @@ def register_elastix(ct, t1c, workdir=None, threads=8):
     ref.SetOrigin(tuple(float(x) for x in lo))
     fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0)   # T1c = reference
     moving = sitk.Resample(ct, ref, sitk.Transform(), sitk.sitkLinear, 0.0)    # CT  = moving
-    # shared centered start: elastix transforms map fixed(T1c) -> moving(CT),
+    fixed_f, moving_f = os.path.join(workdir, "fixed.nii.gz"), os.path.join(workdir, "moving.nii.gz")
-    # i.e. the inverse direction of centroid_offset, hence the negation.
+    sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), fixed_f, True)
-    off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0)
+    sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), moving_f, True)
-    init_tp0 = {
+    # all-ones masks sized to the box (registration on the whole head box)
-        "Transform": "TranslationTransform",
+    ones = sitk.GetImageFromArray(np.ones(shape, np.uint8))
-        "TransformBase": "ElastixTransformBase",
+    ones.CopyInformation(ref)
-        "TransformParameterReadoutFormat": 1,
+    fm, mm = os.path.join(workdir, "fixedmask.nii.gz"), os.path.join(workdir, "movingmask.nii.gz")
-        "NumberOfParameters": 3,
+    sitk.WriteImage(ones, fm, True); sitk.WriteImage(ones, mm, True)
-        "TransformParameters": [-float(off[0]), -float(off[1]), -float(off[2])],
+
-    }
+    params = _elastix_param_file(workdir, fixed_f, moving_f, fm, mm)
-    params = {
+    outdir = os.path.join(workdir, "out")
-        "ImageDimension": 3,
+    shutil.rmtree(outdir, ignore_errors=True); os.makedirs(outdir)
-        "FixedInternalImagePixelType": "float",
+    env = dict(os.environ)
-        "MovingInternalImagePixelType": "float",
+    env["OMP_NUM_THREADS"] = str(threads)
-        "Transform": "EulerTransform",
+    cmd = ["elastix", "-f", fixed_f, "-m", moving_f, "-out", outdir,
-        "TransformBase": "ElastixTransformBase",
+           "-t", params, "-threads", str(threads),
-        "FixedCenteredTransform": True,
+           "-loglevel", "0", "-logtoconsole", "false"]
-        "TransformParameterReadoutFormat": 1,
+    subprocess.run(cmd, check=True, env=env, capture_output=True, text=True)
-        "HowToCombineTransforms": "Compose",
+    # elastix fixed(T1c) -> moving(CT) transform, saved as an ITK transform file
-        "Optimizer": "AdaptiveStochasticGradientDescent",
+    res = os.path.join(outdir, "Result", "0GenericAnatTransform.txt")
-        "MaximumNumberOfIterations": [500, 500, 250],
+    if not os.path.exists(res):
-        "Metric": "AdvancedMattesMutualInformation",
+        raise PipelineUnavailable(f"elastix did not write a result transform: {res}")
-        "NumberOfHistogramBins": 10,
+    cand = sitk.ReadTransform(res)
-        "ImageSampler": "Random",
+    # cand maps fixed(T1c) -> moving(CT) == the contract (T1c->CT); still verify.
-        "ImageSamplerPatchSize": [25, 25, 25],
-        "WriteResultImage": False,
-    }
-    outdir = os.path.join(workdir or d("data/m6reg/elastix"), "itkelastix")
-    os.makedirs(outdir, exist_ok=True)
-    import time, 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(_sitk_to_itk(fixed, itk))
-    reg.SetMovingImage(_sitk_to_itk(moving, itk))
-    reg.SetParameterObject(itk.dict_to_parameter_object(params))
-    reg.SetInitialTransformParameterObject(itk.dict_to_parameter_object(init_tp0))
-    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 registration failed: {e}\n{tail}") from e
-    dt = time.time() - t0
-    # The combination transform maps fixed(T1c) -> moving(CT) == the contract
-    # (T1c->CT); fix_orientation still verifies it against the skull ring.
-    cand = _itk_transform_to_sitk_rigid(reg.GetCombinationTransform())
     T, inverted = fix_orientation(ct, t1c, cand)
-    try:
-        import importlib.metadata
-        version = f"ITK-Elastix {importlib.metadata.version('itk-elastix')}"
-    except Exception:  # noqa
-        version = "ITK-Elastix (python API)"
     meta = {
-        "pipeline": "elastix", "toolkit": version,
+        "pipeline": "elastix", "toolkit": "ITK-Elastix",
         "metric": "MattesMutualInformation (head box, CT space)",
-        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, centroid init)",
+        "transform": "Rigid3DTransform (ASGD)", "convention": "T1c->CT (verified)",
-        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
+        "auto_inverted": inverted, "fixed": "T1c", "moving": "CT",
-        "fixed": "T1c", "moving": "CT", "seconds_reg": round(dt, 1),
     }
     return T, meta
 
@@ -597,13 +547,8 @@ def availability():
             if name == "simpleitk":
                 out[name] = (True, f"SimpleITK {sitk.Version_VersionString()}")
             elif name == "elastix":
-                try:
+                ok = shutil.which("elastix") is not None
-                    import importlib.metadata
+                out[name] = (ok, "binaries present" if ok else "elastix/transformix missing")
-                    import itk
-                    _ = itk.ElastixRegistrationMethod
-                    out[name] = (True, f"ITK-Elastix {importlib.metadata.version('itk-elastix')} (python API)")
-                except Exception as e:  # noqa
-                    out[name] = (False, f"itk-elastix python API missing ({type(e).__name__})")
             elif name == "antspyx":
                 try:
                     import ants  # noqa