From a9a4d9bba023c50e327914550309544fa396aab3 Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Mon, 28 Sep 2026 19:01:37 +0800 Subject: [PATCH] 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 | 2 + REGISTRATION_BENCHMARK.md | 131 ++++++++++++--------- scripts/bench_registration.py | 15 +-- src/registration_pipelines.py | 207 +++++++++++++++++++++------------- 4 files changed, 219 insertions(+), 136 deletions(-) diff --git a/.gitignore b/.gitignore index 3d000e0..4664ead 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,5 @@ runs_bench/ runs_nnu/ runs_monai/ __pycache__/ + +elastix/ \ No newline at end of file diff --git a/REGISTRATION_BENCHMARK.md b/REGISTRATION_BENCHMARK.md index 41971b1..d37bd52 100644 --- a/REGISTRATION_BENCHMARK.md +++ b/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) | @@ -84,58 +84,71 @@ Production QA, crop-independent: | 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.1 | 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.6 | 0.93 | **0.643 / 0.630** | | **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 | — | — | — | — | — | — | — | — | Per-case (mid / std, `OK` = accepted): -| case | simpleitk | antspyx | -|---|---|---| -| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | -| m6_22IRRKU6_20240115 | 25.7 / 29.3 | 26.3 / 29.2 | -| m6_23XIMO4M_20230222 | 22.2 / 23.9 | 23.8 / 25.0 | -| m6_24V62QQM_20230223 | 10.7 / 16.9 | 9.4 / 18.2 | -| 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** | -| m6_26TYQI3C_20221027 | 12.2 / 3.7 **OK** | 6.1 / 37.2 | -| m6_2BBN5XA4_20250411 | 18.6 / 8.0 **OK** | 14.3 / 12.7 | -| m6_2C4TKI6H_20240506 | 12.6 / 15.4 | 13.6 / 14.7 | -| m6_2CDNN7BH_20240102 | 19.4 / 30.2 | 18.2 / 33.6 | +| case | simpleitk | antspyx | elastix | +|---|---|---|---| +| m6_22GD5VL2_20230724 | incomplete-CT | incomplete-CT | incomplete-CT | +| 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 | 21.1 / 23.3 | +| 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 | 11.0 / 6.1 **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 | 11.4 / 3.2 **OK** | +| 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 | 12.0 / 6.8 **OK** | +| 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 - `≤ 20 mm` band). `tv-in-CT` is identical (0.93): the warped GTV lands inside - the CT head for both, so the target is in the right place. +- **Translation is recovered well by all three** (mid gaps mostly 6–20 mm, + within the `≤ 20 mm` band). `tv-in-CT` is identical (~0.93): the warped GTV + 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 - std 15.4 vs 18.2): its RSGD + physical-step-scales + 32-bin Mattes-MI - landscape handles the tilt better than ANTs' Rigid. Note `m6_26TYQI3C`: ANTs - nailed the position (mid 6.1) but missed the orientation (std 37.2) — a - textbook tilted local minimum. -- **Speed**: SimpleITK ~8.1 s/case, ANTs ~12.7 s/case (ANTs includes the - pre-shift resample + its affine-MI init stages). Both are CPU-bound and ran - 4 cases in parallel. +- **Elastix is the most accurate of the three standard single-init configs** + (5/10 accepted, lowest median std 9.9 vs SimpleITK 15.4 and ANTs 24.8) and + the fastest (3.2 s/case median wall, ~1.8 s pure registration vs 8.6 s + SimpleITK / 12.7 s ANTs): ASGD's adaptive step sizes climb out of small + misalignments better than RSGD/ANTs on these boxes. It accepts `m6_24V62QQM_20230223` + (std 8.2) and `m6_2C4TKI6H` (std 6.8), which the other two tools fail on + (std 16.9 / 15.4), and misses `m6_24V62QQM_20240820` (std 14.0) that both + 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**: - SimpleITK median 0.643 nats vs ANTs 0.600, 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. -- **ANTs MI is stochastic**: 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. + 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) @@ -173,29 +186,18 @@ 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 two available pipelines with: +visuals agree. Generate for the available pipelines with: ```bash /home/xfr/.conda/envs/longitudinal/bin/python scripts/bench_registration.py \ - --pipelines simpleitk --n-cases 10 --qa-images + --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 ``` (9 images each; the incomplete-CT case is skipped.) -## Enabling the two blocked pipelines - -**ITK-Elastix** — not on PyPI; needs conda-forge (which this host cannot reach). -Once conda-forge is available: - -```bash -conda install -c conda-forge elastix # provides the elastix/transformix CLI -python scripts/bench_registration.py --pipelines elastix --n-cases 10 -``` - -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 @@ -217,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 - "not installed", not measured. Their code paths are implemented and - unit-verified for import/convention; run them where the toolkits are present. +- `deepreg` could not be executed here, so its row is "not installed", not + measured; its code path is implemented and unit-verified for import/convention. - `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. \ No newline at end of file + 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. \ No newline at end of file diff --git a/scripts/bench_registration.py b/scripts/bench_registration.py index d62cd04..e3d98cd 100644 --- a/scripts/bench_registration.py +++ b/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 -`longitudinal` env (elastix needs conda-forge; antspyx has no cp314 wheel; -PyPI deepreg is a broken stub), the benchmark is env-portable: run it under -whichever conda env has a given toolkit installed, and it records only the -pipelines that are importable there. The `--merge` mode then combines the -per-pipeline JSON files (accumulated across envs) into one summary + plot. +Because two of the four toolkits are not pip-installable into the main +`longitudinal` env (antspyx has no cp314 wheel; PyPI deepreg is a broken +stub), the benchmark is env-portable: run it under whichever conda env has a +given toolkit installed, and it records only the pipelines that are +importable there. elastix runs in the main env through the ITK-Elastix +python API (pip ``itk`` + ``itk-elastix``). 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 diff --git a/src/registration_pipelines.py b/src/registration_pipelines.py index df93223..a3cba8c 100644 --- a/src/registration_pipelines.py +++ b/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): - """Rigid3DTransform + Mattes MI parameter file.""" - lines = [ - "(ImageDimension 3)", - "(FixedInternalImagePixelType float)", - "(MovingInternalImagePixelType float)", - "(Transform \"Rigid3DTransform\")", - "(TransformBase \"ElastixTransformBase\")", - "(FixedCenteredTransform 1)", - "(TransformParameterReadoutFormat 1)", - "(Optimizer \"AdaptiveStochasticGradientDescent\")", - "(MaximumNumberOfIterations 6000)", - "(StoppingValue 0.0005)", - "(StoppingIterationCount 2500)", - "(ImageMetric \"MattesMutualInformation\")", - "(NumberOfHistogramBins 10)", - "(ImageSampler \"RandomSampler\")", - "(ImageSamplerPatchSize 25 25 25)", - "(FixedSubstresamplerClass \"RandomSubstresampler\")", - "(FixedSubstresamplerScale 2)", - "(ResampleInterpolator \"FinalBSplineInterpolator\")", - "(FinalBSplineInterpolationOrder 1)", - "(ResultStationaryImage \"result.mha\")", - "(WriteResultImage 0)", - "(CompressResultImages 1)", - f"(FixedImage {fixed})", - f"(FixedImageMask {fixed_mask})", - f"(MovingImage {moving})", - f"(MovingImageMask {moving_mask})", - ] - p = os.path.join(workdir, "elastix_params.txt") - with open(p, "w") as f: - f.write("\n".join(lines) + "\n") - return p +def _elastix_api(): + """Import the ITK-Elastix python API, or raise PipelineUnavailable.""" + try: + import itk + except Exception as e: # noqa + raise PipelineUnavailable( + f"itk not importable ({e!r}). Install with `pip install itk itk-elastix` " + "(ITK-Elastix >= 0.25 exposes the elastix components in the `itk` " + "python namespace; no CLI / conda-forge needed).") + try: + _ = itk.ElastixRegistrationMethod # loads the elastix module + components + except Exception as e: # noqa + raise PipelineUnavailable( + f"ITK-Elastix components not registered in the `itk` python namespace " + f"({type(e).__name__}: {e}). `pip install itk-elastix` (with a " + "matching `itk`) provides `itk.ElastixRegistrationMethod`.") + return itk + + +def _sitk_to_itk(img, itk): + """SimpleITK image -> itk python image (float32, geometry preserved).""" + out = itk.GetImageFromArray(sitk.GetArrayFromImage(sitk.Cast(img, sitk.sitkFloat32))) + out.SetSpacing([float(s) for s in img.GetSpacing()]) + out.SetOrigin([float(o) for o in img.GetOrigin()]) + out.SetDirection(itk.GetMatrixFromArray( + np.array(img.GetDirection(), dtype=float).reshape(3, 3))) + return out + + +def _itk_transform_to_sitk_rigid(itk_tf): + """Recover the affine (R, t) of an ITK transform by point mapping and build + a SimpleITK transform from it. The mapping is p -> R p + t in world + coordinates, so a plain AffineTransform represents the rigid result exactly + (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.""" - if shutil.which("elastix") is None: - raise PipelineUnavailable( - "elastix/transformix binaries not found. Install with " - "`conda install -c conda-forge elastix` (conda-forge must be " - "reachable) or build from https://elastix.org. Not on PyPI.") - workdir = workdir or d("data/m6reg/elastix") - workdir = os.path.join(workdir, f"pid{os.getpid()}") # per-worker scratch - os.makedirs(workdir, exist_ok=True) + """ITK-Elastix rigid CT->T1c (fixed=T1c, moving=CT), Mattes MI. + + Runs in-process through the ITK-Elastix python API (pip ``itk`` + + ``itk-elastix``, elastix 5.x). Rigid 6-DOF via ``EulerTransform`` (the + 5.x replacement for ``Rigid3DTransform``) + Mattes MI + (``AdvancedMattesMutualInformation``) + ASGD, 3 resolutions, on the common + CT-space head box, from the *shared* centroid-aligned start + (``centroid_offset`` as an initial TranslationTransform, same init the + 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 - # choice) so MI is computed on head tissue, not neck/air. + # Resample both into a common CT-space head box (mirrors the other pipelines) + # 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") - sitk.WriteImage(sitk.Cast(fixed, sitk.sitkFloat32), fixed_f, True) - sitk.WriteImage(sitk.Cast(moving, sitk.sitkFloat32), moving_f, True) - # all-ones masks sized to the box (registration on the whole head box) - ones = sitk.GetImageFromArray(np.ones(shape, np.uint8)) - ones.CopyInformation(ref) - fm, mm = os.path.join(workdir, "fixedmask.nii.gz"), os.path.join(workdir, "movingmask.nii.gz") - sitk.WriteImage(ones, fm, True); sitk.WriteImage(ones, mm, True) - - params = _elastix_param_file(workdir, fixed_f, moving_f, fm, mm) - outdir = os.path.join(workdir, "out") - shutil.rmtree(outdir, ignore_errors=True); os.makedirs(outdir) - env = dict(os.environ) - env["OMP_NUM_THREADS"] = str(threads) - cmd = ["elastix", "-f", fixed_f, "-m", moving_f, "-out", outdir, - "-t", params, "-threads", str(threads), - "-loglevel", "0", "-logtoconsole", "false"] - subprocess.run(cmd, check=True, env=env, capture_output=True, text=True) - # elastix fixed(T1c) -> moving(CT) transform, saved as an ITK transform file - res = os.path.join(outdir, "Result", "0GenericAnatTransform.txt") - if not os.path.exists(res): - raise PipelineUnavailable(f"elastix did not write a result transform: {res}") - cand = sitk.ReadTransform(res) - # cand maps fixed(T1c) -> moving(CT) == the contract (T1c->CT); still verify. + # shared centered start: elastix transforms map fixed(T1c) -> moving(CT), + # i.e. the inverse direction of centroid_offset, hence the negation. + off = centroid_offset(ct, t1c, ct_box=moving) or (0.0, 0.0, 0.0) + init_tp0 = { + "Transform": "TranslationTransform", + "TransformBase": "ElastixTransformBase", + "TransformParameterReadoutFormat": 1, + "NumberOfParameters": 3, + "TransformParameters": [-float(off[0]), -float(off[1]), -float(off[2])], + } + params = { + "ImageDimension": 3, + "FixedInternalImagePixelType": "float", + "MovingInternalImagePixelType": "float", + "Transform": "EulerTransform", + "TransformBase": "ElastixTransformBase", + "FixedCenteredTransform": True, + "TransformParameterReadoutFormat": 1, + "HowToCombineTransforms": "Compose", + "Optimizer": "AdaptiveStochasticGradientDescent", + "MaximumNumberOfIterations": [500, 500, 250], + "Metric": "AdvancedMattesMutualInformation", + "NumberOfHistogramBins": 10, + "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)", - "auto_inverted": inverted, "fixed": "T1c", "moving": "CT", + "transform": "Rigid (EulerTransform, ASGD, 3 resolutions, centroid init)", + "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 - out[name] = (ok, "binaries present" if ok else "elastix/transformix missing") + try: + 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