refactor(benchmark): use ITK Python API for elastix registration

Replace the conda-forge-dependent Elastix CLI path with the pip-installable
ITK-Elastix Python components so the benchmark can run elastix in the main
environment. Add SimpleITK/ITK image and transform conversion helpers, and
use the elastix 5.x EulerTransform rigid pipeline with Mattes MI and ASGD.

Update REGISTRATION_BENCHMARK.md and bench_registration.py documentation to
reflect the new elastix availability, results, and multi-environment merge
workflow. Ignore local elastix output directories in .gitignore.

.gitignore

@@ -7,3 +7,5 @@ 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 | Mattes MI | `Rigid3DTransform` (ASGD) | **not installed** (conda-forge unreachable) |
+| 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) |
 | 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,56 +71,84 @@ 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 |
+| pipeline | toolkit | registered | accepted | gap mid (med/mean) | gap std (med/mean) | s/reg (med) | tv-in-CT | MI (med/mean, nats) |
-|---|---|---|---|---|---|---|---|
+|---|---|---|---|---|---|---|---|---|
-| **simpleitk** | SimpleITK 2.5.6 | 9 | **4/10 (0.40)** | 12.6 / 16.1 | 15.4 / 16.0 | 8.3 | 0.93 |
+| **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** |
-| **antspyx** | ANTsPy 0.6.3 | 9 | **1/10 (0.10)** | 14.3 / 15.4 | 18.2 / 21.7 | 13.4 | 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** |
-| elastix | — | — | — | — | — | — | — |
+| **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** |
-| deepreg | — | — | — | — | — | — | — |
+| deepreg | — | — | — | — | — | — | — | — |
 
 Per-case (mid / std, `OK` = accepted):
 
-| case | simpleitk | antspyx |
+| case | simpleitk | antspyx | elastix |
-|---|---|---|
+|---|---|---|---|
-| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT |
+| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT |
-| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 26.3 / 29.2 |
+| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 26.3 / 29.2 | 26.4 / 30.3 |
-| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.8 / 25.0 |
+| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.8 / 25.0 | 21.1 / 23.3 |
-| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.4 / 18.2 |
+| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.4 / 18.2 | 10.6 / 8.2 **OK** |
-| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 12.5 / 12.8 |
+| m6_24V62QQM_20231017 | 11.5 / 7.3 **OK** | 12.5 / 12.8 | 11.0 / 6.1 **OK** |
-| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.8 / 11.7 **OK** |
+| m6_24V62QQM_20240820 | 11.6 / 9.6 **OK** | 14.8 / 11.7 **OK** | 15.3 / 14.0 |
-| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 6.1 / 37.2 |
+| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 6.1 / 37.2 | 11.4 / 3.2 **OK** |
-| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 |
+| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | 18.1 / 9.9 **OK** |
-| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.6 / 14.7 |
+| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.6 / 14.7 | 12.0 / 6.8 **OK** |
-| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 18.2 / 33.6 |
+| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 18.2 / 33.6 | 18.9 / 22.7 |
 
 ### Interpretation
 
-- **Translation is recovered well by both** (mid gaps mostly 6–20 mm, within the
+- **Translation is recovered well by all three** (mid gaps mostly 6–20 mm,
-  `≤ 20 mm` band). `tv-in-CT` is identical (0.93): the warped GTV lands inside
+  within the `≤ 20 mm` band). `tv-in-CT` is identical (~0.93): the warped GTV
-  the CT head for both, so the target is in the right place.
+  lands inside the CT head for every pipeline, so the target is in the right
+  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**.
-- **SimpleITK is the more robust of the two** (4 vs 1 accepted, lower median
+- **Elastix is the most accurate of the three standard single-init configs**
-  std 15.4 vs 18.2): its RSGD + physical-step-scales + 32-bin Mattes-MI
+  (5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 24.8) and
-  landscape handles the tilt better than ANTs' Rigid. Note `m6_26TYQI3C`: ANTs
+  the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s
-  nailed the position (mid 6.1) but missed the orientation (std 37.2) — a
+  SimpleITK / 12.7 s ANTs): ASGD's adaptive step sizes climb out of small
-  textbook tilted local minimum.
+  misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223`
-- **Speed**: SimpleITK ~8.3 s/case, ANTs ~13.4 s/case (ANTs includes the
+  (std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on
-  pre-shift resample + its affine-MI init stages). Both are CPU-bound and ran
+  (std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both
-  4 cases in parallel.
+  others accept — different optimizers settle in different local minima on the
+  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 baseline (elastix/deepreg report unavailable here)
+# main env — SimpleITK + Elastix (itk-elastix is pip-installed 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 --pipelines simpleitk --n-cases 10 \
+python scripts/bench_registration.py --n-cases 10 \
     --workers 4 --threads 24 --simpleitk-mode benchmark
 
 # isolated py3.12 env with antspyx — ANTs (+ optional SimpleITK cross-check)
@@ -134,20 +162,42 @@ 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).
 
-## Enabling the two blocked pipelines
+## Visual QA images (checkerboard + color overlay)
 
-**ITK-Elastix** — not on PyPI; needs conda-forge (which this host cannot reach).
+`--qa-images` writes a `checkreg`-style PNG per (pipeline, case) to
-Once conda-forge is available:
+`<out-dir>/qa/<pipeline>__<case>.png` for eyeballing the alignment.  Each image
+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
-conda install -c conda-forge elastix        # provides the elastix/transformix CLI
+/home/xfr/.conda/envs/longitudinal/bin/python scripts/bench_registration.py \
-python scripts/bench_registration.py --pipelines elastix --n-cases 10
+    --pipelines simpleitk elastix --n-cases 10 --qa-images
+/home/xfr/.conda/envs/ants312/bin/python scripts/bench_registration.py \
+    --pipelines antspyx --n-cases 10 --qa-images
 ```
 
-The adapter registers `fixed=T1c, moving=CT` with `Rigid3DTransform` + Mattes MI
+(9 images each; the incomplete-CT case is skipped.)
-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
@@ -169,9 +219,32 @@ 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.
-- `elastix` / `deepreg` could not be executed here, so their rows are
+- `deepreg` could not be executed here, so its row is "not installed", not
-  "not installed", not measured. Their code paths are implemented and
+  measured; its code path is implemented and unit-verified for import/convention.
-  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.
+  independent of each toolkit's transform-direction convention. In practice it
+  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,15 +17,16 @@ 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 three of the four toolkits are not pip-installable into the main
+Because two of the four toolkits are not pip-installable into the main
-`longitudinal` env (elastix needs conda-forge; antspyx has no cp314 wheel;
+`longitudinal` env (antspyx has no cp314 wheel; PyPI deepreg is a broken
-PyPI deepreg is a broken stub), the benchmark is env-portable: run it under
+stub), the benchmark is env-portable: run it under whichever conda env has a
-whichever conda env has a given toolkit installed, and it records only the
+given toolkit installed, and it records only the pipelines that are
-pipelines that are importable there.  The `--merge` mode then combines the
+importable there.  elastix runs in the main env through the ITK-Elastix
-per-pipeline JSON files (accumulated across envs) into one summary + plot.
+python API (pip ``itk`` + ``itk-elastix``).  The `--merge` mode then combines
+the per-pipeline JSON files (accumulated across envs) into one summary + plot.
 
 Usage (repo root):
-  # main env (SimpleITK baseline; elastix/deepreg report as unavailable)
+  # main env (SimpleITK + elastix; deepreg reports as unavailable)
   source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal
   python scripts/bench_registration.py --pipelines simpleitk elastix deepreg
 
@@ -52,10 +53,34 @@ import numpy as np
 import SimpleITK as sitk
 from scipy import ndimage
 
-from src.common import DATA, d, read_nii_img
+from src.common import DATA, d, read_nii_img, head_mask_from_image
 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()
 
@@ -106,6 +131,10 @@ 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)
@@ -119,6 +148,7 @@ 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",
@@ -126,6 +156,64 @@ 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)
@@ -158,6 +246,18 @@ 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
 
 
@@ -184,6 +284,7 @@ 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)})
@@ -197,6 +298,7 @@ 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"),
@@ -212,6 +314,8 @@ 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"),
     }
@@ -238,23 +342,24 @@ 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}")
+            f"{'accept':>7s} {'gapMid':>7s} {'gapStd':>7s} {'s/reg':>7s} {'tvInCT':>7s} {'MI':>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}   [{r['reason']}]")
+                  f"{'-':>7s} {'-':>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}")
+              f"{_fmt(tv, '{:.2f}'):>7s} {_fmt(mi, '{:.3f}'):>7s}")
     print("=" * 100)
 
 
@@ -275,13 +380,15 @@ 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]
-    fig, axes = plt.subplots(1, 4, figsize=(18, 4))
+    mi = [ok_rows[p].get("mi_ct_mr", {}).get("median") or 0 for p in labels]
-    colors = ["#4878CF", "#D64F5C", "#6ACC65", "#9467BD"]
+    fig, axes = plt.subplots(1, 5, figsize=(20, 4))
+    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()),
@@ -325,12 +432,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 |",
+    md = ["| pipeline | available | toolkit | cases | ok | err | accepted | gap mid (mm) | gap std (mm) | s/reg (median) | tv-in-CT | MI (nats) |",
-          "|---|---|---|---|---|---|---|---|---|---|---|"]
+          "|---|---|---|---|---|---|---|---|---|---|---|---|"]
     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)} | "
@@ -338,7 +445,8 @@ 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")
@@ -363,6 +471,11 @@ 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
@@ -375,6 +488,9 @@ 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,7 +26,8 @@ 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 CLI (Rigid3DTransform, Mattes MI).
+- ``register_elastix``   : ITK-Elastix Python API (pip ``itk`` + ``itk-elastix``,
+  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).
@@ -49,7 +50,6 @@ SimpleITK path.
 """
 import os
 import json
-import shutil
 import subprocess
 import importlib.util
 
@@ -242,57 +242,74 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
 # 1. ITK-Elastix
 # --------------------------------------------------------------------------- #
 
-def _elastix_param_file(workdir, fixed, moving, fixed_mask, moving_mask):
+def _elastix_api():
-    """Rigid3DTransform + Mattes MI parameter file."""
+    """Import the ITK-Elastix python API, or raise PipelineUnavailable."""
-    lines = [
+    try:
-        "(ImageDimension 3)",
+        import itk
-        "(FixedInternalImagePixelType float)",
+    except Exception as e:  # noqa
-        "(MovingInternalImagePixelType float)",
+        raise PipelineUnavailable(
-        "(Transform \"Rigid3DTransform\")",
+            f"itk not importable ({e!r}). Install with `pip install itk itk-elastix` "
-        "(TransformBase \"ElastixTransformBase\")",
+            "(ITK-Elastix >= 0.25 exposes the elastix components in the `itk` "
-        "(FixedCenteredTransform 1)",
+            "python namespace; no CLI / conda-forge needed).")
-        "(TransformParameterReadoutFormat 1)",
+    try:
-        "(Optimizer \"AdaptiveStochasticGradientDescent\")",
+        _ = itk.ElastixRegistrationMethod  # loads the elastix module + components
-        "(MaximumNumberOfIterations 6000)",
+    except Exception as e:  # noqa
-        "(StoppingValue 0.0005)",
+        raise PipelineUnavailable(
-        "(StoppingIterationCount 2500)",
+            f"ITK-Elastix components not registered in the `itk` python namespace "
-        "(ImageMetric \"MattesMutualInformation\")",
+            f"({type(e).__name__}: {e}). `pip install itk-elastix` (with a "
-        "(NumberOfHistogramBins 10)",
+            "matching `itk`) provides `itk.ElastixRegistrationMethod`.")
-        "(ImageSampler \"RandomSampler\")",
+    return itk
-        "(ImageSamplerPatchSize 25 25 25)",
+
-        "(FixedSubstresamplerClass \"RandomSubstresampler\")",
+
-        "(FixedSubstresamplerScale 2)",
+def _sitk_to_itk(img, itk):
-        "(ResampleInterpolator \"FinalBSplineInterpolator\")",
+    """SimpleITK image -> itk python image (float32, geometry preserved)."""
-        "(FinalBSplineInterpolationOrder 1)",
+    out = itk.GetImageFromArray(sitk.GetArrayFromImage(sitk.Cast(img, sitk.sitkFloat32)))
-        "(ResultStationaryImage \"result.mha\")",
+    out.SetSpacing([float(s) for s in img.GetSpacing()])
-        "(WriteResultImage 0)",
+    out.SetOrigin([float(o) for o in img.GetOrigin()])
-        "(CompressResultImages 1)",
+    out.SetDirection(itk.GetMatrixFromArray(
-        f"(FixedImage {fixed})",
+        np.array(img.GetDirection(), dtype=float).reshape(3, 3)))
-        f"(FixedImageMask {fixed_mask})",
+    return out
-        f"(MovingImage {moving})",
+
-        f"(MovingImageMask {moving_mask})",
+
-    ]
+def _itk_transform_to_sitk_rigid(itk_tf):
-    p = os.path.join(workdir, "elastix_params.txt")
+    """Recover the affine (R, t) of an ITK transform by point mapping and build
-    with open(p, "w") as f:
+    a SimpleITK transform from it.  The mapping is p -> R p + t in world
-        f.write("\n".join(lines) + "\n")
+    coordinates, so a plain AffineTransform represents the rigid result exactly
-    return p
+    (independent of the toolkit's internal angle parameterization)."""
+    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:
+
-        raise PipelineUnavailable(
+    Runs in-process through the ITK-Elastix python API (pip ``itk`` +
-            "elastix/transformix binaries not found. Install with "
+    ``itk-elastix``, elastix 5.x).  Rigid 6-DOF via ``EulerTransform`` (the
-            "`conda install -c conda-forge elastix` (conda-forge must be "
+    5.x replacement for ``Rigid3DTransform``) + Mattes MI
-            "reachable) or build from https://elastix.org. Not on PyPI.")
+    (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common
-    workdir = workdir or d("data/m6reg/elastix")
+    CT-space head box, from the *shared* centroid-aligned start
-    workdir = os.path.join(workdir, f"pid{os.getpid()}")  # per-worker scratch
+    (``centroid_offset`` as an initial TranslationTransform, same init the
-    os.makedirs(workdir, exist_ok=True)
+    SimpleITK/ANTs benchmark modes use).  Elastix 5.x component names: sampler
+    ``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 production grid
+    # Resample both into a common CT-space head box (mirrors the other pipelines)
-    # choice) so MI is computed on head tissue, not neck/air.
+    # 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))
@@ -301,36 +318,69 @@ 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
-    fixed_f, moving_f = os.path.join(workdir, "fixed.nii.gz"), os.path.join(workdir, "moving.nii.gz")
+    # shared centered start: elastix transforms map fixed(T1c) -> moving(CT),
-    sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), fixed_f, True)
+    # i.e. the inverse direction of centroid_offset, hence the negation.
-    sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), moving_f, True)
+    off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0)
-    # all-ones masks sized to the box (registration on the whole head box)
+    init_tp0 = {
-    ones = sitk.GetImageFromArray(np.ones(shape, np.uint8))
+        "Transform": "TranslationTransform",
-    ones.CopyInformation(ref)
+        "TransformBase": "ElastixTransformBase",
-    fm, mm = os.path.join(workdir, "fixedmask.nii.gz"), os.path.join(workdir, "movingmask.nii.gz")
+        "TransformParameterReadoutFormat": 1,
-    sitk.WriteImage(ones, fm, True); sitk.WriteImage(ones, mm, True)
+        "NumberOfParameters": 3,
-
+        "TransformParameters": [-float(off[0]), -float(off[1]), -float(off[2])],
-    params = _elastix_param_file(workdir, fixed_f, moving_f, fm, mm)
+    }
-    outdir = os.path.join(workdir, "out")
+    params = {
-    shutil.rmtree(outdir, ignore_errors=True); os.makedirs(outdir)
+        "ImageDimension": 3,
-    env = dict(os.environ)
+        "FixedInternalImagePixelType": "float",
-    env["OMP_NUM_THREADS"] = str(threads)
+        "MovingInternalImagePixelType": "float",
-    cmd = ["elastix", "-f", fixed_f, "-m", moving_f, "-out", outdir,
+        "Transform": "EulerTransform",
-           "-t", params, "-threads", str(threads),
+        "TransformBase": "ElastixTransformBase",
-           "-loglevel", "0", "-logtoconsole", "false"]
+        "FixedCenteredTransform": True,
-    subprocess.run(cmd, check=True, env=env, capture_output=True, text=True)
+        "TransformParameterReadoutFormat": 1,
-    # elastix fixed(T1c) -> moving(CT) transform, saved as an ITK transform file
+        "HowToCombineTransforms": "Compose",
-    res = os.path.join(outdir, "Result", "0GenericAnatTransform.txt")
+        "Optimizer": "AdaptiveStochasticGradientDescent",
-    if not os.path.exists(res):
+        "MaximumNumberOfIterations": [500, 500, 250],
-        raise PipelineUnavailable(f"elastix did not write a result transform: {res}")
+        "Metric": "AdvancedMattesMutualInformation",
-    cand = sitk.ReadTransform(res)
+        "NumberOfHistogramBins": 10,
-    # cand maps fixed(T1c) -> moving(CT) == the contract (T1c->CT); still verify.
+        "ImageSampler": "Random",
+        "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": "ITK-Elastix",
+        "pipeline": "elastix", "toolkit": version,
         "metric": "MattesMutualInformation (head box, CT space)",
-        "transform": "Rigid3DTransform (ASGD)", "convention": "T1c->CT (verified)",
+        "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, centroid init)",
-        "auto_inverted": inverted, "fixed": "T1c", "moving": "CT",
+        "convention": "T1c->CT (verified)", "auto_inverted": inverted,
+        "fixed": "T1c", "moving": "CT", "seconds_reg": round(dt, 1),
     }
     return T, meta
 
@@ -547,8 +597,13 @@ def availability():
             if name == "simpleitk":
                 out[name] = (True, f"SimpleITK {sitk.Version_VersionString()}")
             elif name == "elastix":
-                ok = shutil.which("elastix") is not None
+                try:
-                out[name] = (ok, "binaries present" if ok else "elastix/transformix missing")
+                    import importlib.metadata
+                    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