xfr/longitudinal
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

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

Browse source 
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.
This commit is contained in:
Furen Xiao 2026-09-28 19:01:37 +08:00
parent 214afd643b
commit a9a4d9bba0
4 changed files with 219 additions and 136 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
.gitignore vendored
View file

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

131
REGISTRATION_BENCHMARK.md
View file

@ -32,7 +32,7 @@ handling — not on who finds the ~30–60 mm gross offset first.
| # | pipeline | toolkit | metric | transform | status in this env | | # | 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 | | 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) | | 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) | | 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) | | 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** | | **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): 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 ### 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 binding failure is rotation** — `std > 12 mm` on the repositioned cases
(the ~30° inter-session tilt produces MI local minima that a single centered (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` start can lock onto). This is exactly why production `02_build_m6_dataset.py`
uses a **multi-start tilt grid**. 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.1 s/case, ANTs ~12.7 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**: - **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 medians 0.643 (SimpleITK) / 0.600 (ANTs) / 0.562 (Elastix), and it diverges
tilted cases — ANTs `m6_26TYQI3C` drops to 0.094 nats (rotation ~37 mm off) sharply on the tilted cases — ANTs `m6_26TYQI3C` drops to 0.094 nats
while SimpleITK holds 0.670. A high MI means the CT skull/tissue intensities (rotation ~37 mm off) while SimpleITK holds 0.670. A high MI means the CT
fall on the matching MR structures after the inverse transform. skull/tissue intensities fall on the matching MR structures after the inverse
- **ANTs MI is stochastic**: ANTs' Mattes-MI subsamples the image at random, so transform.
its per-case numbers carry run-to-run variance (this table is one draw); - **Stochasticity**: ANTs' Mattes-MI subsamples the image at random, so its
SimpleITK's MI (and registration) is deterministic. 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 ## Reproducing
```bash ```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 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 --workers 4 --threads 24 --simpleitk-mode benchmark
# isolated py3.12 env with antspyx — ANTs (+ optional SimpleITK cross-check) # 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 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 (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 ```bash
/home/xfr/.conda/envs/longitudinal/bin/python scripts/bench_registration.py \ /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 \ /home/xfr/.conda/envs/ants312/bin/python scripts/bench_registration.py \
--pipelines antspyx --n-cases 10 --qa-images --pipelines antspyx --n-cases 10 --qa-images
``` ```
(9 images each; the incomplete-CT case is skipped.) (9 images each; the incomplete-CT case is skipped.)
## Enabling the two blocked pipelines ## Enabling the blocked pipeline
**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.
**DeepReg (Oxford)** — the PyPI `deepreg` (0.1.2) is an unrelated stub that pins **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 `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, 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 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. 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 - `fix_orientation` (elastix/ants/deepreg) chooses between a transform and its
inverse by skull-ring concentricity — a correctness gate that makes the result 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.

15
scripts/bench_registration.py
View file

@ -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). - accepted : production rule (nvox>0 and tv_in_ct and concentric).
- reg_wall_s : wall-clock registration time (excludes I/O + scoring). - 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): 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 source /opt/conda/etc/profile.d/conda.sh && conda activate longitudinal
python scripts/bench_registration.py --pipelines simpleitk elastix deepreg python scripts/bench_registration.py --pipelines simpleitk elastix deepreg

207
src/registration_pipelines.py
View file

@ -26,7 +26,8 @@ Pipelines
- ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a - ``register_simpleitk`` : SimpleITK ImageRegistrationMethod, Mattes MI on a
CT-space head box, multi-start tilt grid. Reuses the *proven* production 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_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_antspyx`` : ANTs / ANTsPy ``ants.registration`` (Rigid).
- ``register_deepreg`` : Oxford DeepReg (CycleGAN deformation-field - ``register_deepreg`` : Oxford DeepReg (CycleGAN deformation-field
registration; needs a pretrained checkpoint, see module docstring). registration; needs a pretrained checkpoint, see module docstring).
@ -49,7 +50,6 @@ SimpleITK path.
""" """
import os import os
import json import json
import shutil
import subprocess import subprocess
import importlib.util import importlib.util
@ -242,57 +242,74 @@ def register_simpleitk(ct, t1c, threads=12, mode="benchmark", iters_full=150):
# 1. ITK-Elastix # 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): 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) lo, hi = _ct_head_box_mm(ct)
if lo is None: if lo is None:
raise PipelineUnavailable("elastix: no CT head box found (corrupt CT?)") 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()) tsp = np.array(t1c.GetSpacing())
shape = tuple(int(np.ceil((hi[i] - lo[i]) / tsp[i])) + 1 for i in range(3)) 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)) 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)) ref.SetOrigin(tuple(float(x) for x in lo))
fixed = sitk.Resample(t1c, ref, sitk.Transform(), sitk.sitkLinear, 0.0) # T1c = reference 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 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) 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 = { meta = {
"pipeline": "elastix", "toolkit": "ITK-Elastix", "pipeline": "elastix", "toolkit": version,
"metric": "MattesMutualInformation (head box, CT space)", "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 return T, meta
@ -547,8 +597,13 @@ def availability():
if name == "simpleitk": if name == "simpleitk":
out[name] = (True, f"SimpleITK {sitk.Version_VersionString()}") out[name] = (True, f"SimpleITK {sitk.Version_VersionString()}")
elif name == "elastix": 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": elif name == "antspyx":
try: try:
import ants # noqa import ants # noqa