BUG:ACCESS_VIOLATION crashes and NumPy shape bug in 4D CT pipeline (#68)

* BUG:ACCESS_VIOLATION crashes and NumPy shape bug in 4D CT pipeline

Four bugs fixed:

1. NumPy array shape wrong in create_distance_map (contour_tools.py)
   ITK GetSize() returns (Nx, Ny, Nz) but NumPy needs (Nz, Ny, Nx)
   ordering. np.zeros(size) was creating the wrong shape, causing an
   IndexError when ix >= size[2].

2. TubeTK ResampleImage crashes ~80% on Windows (workflow_fit_statistical_model_to_patient.py)
   ITK's lazy SWIG DLL loader for TubeTK is intermittently unstable on
   Windows. Replace TubeTK.ResampleImage.SetMakeHighResIso() with a new
   _make_isotropic_image() helper using standard ITK ResampleImageFilter
   + LinearInterpolateImageFunction. Equivalent output, no TubeTK DLL
   needed for this path.

3. GPU libraries initialized at import time conflict with TubeTK DLL loading
   Several modules imported CUDA-initializing libraries at module level,
   which interfered with TubeTK's DLL initialization on Windows:
   - register_images_icon.py: torch/icon_registration/unigradicon moved
     into a _load_icon() lazy loader function
   - segment_chest_total_segmentator.py: totalsegmentator import moved
     inside segmentation_method() where it is used
   - __init__.py: `import cupy` (which initializes CUDA at import time)
     replaced with importlib.util.find_spec("cupy") existence check
   - transform_tools.py: `import cupy` moved inside the one method
     that uses it for GPU-accelerated norm computation

4. ITK lazy-load reference dropped before .New() (workflow_fit_statistical_model_to_patient.py)
   ITK 6.0b2's __getattr__ returns a temporary that can be GC'd before
   .New() completes. Fixed by storing the class reference before calling
   .New() (now superseded by the TubeTK replacement in bug #2).

* BUG: Remove test of removed function

* BUG: Suggestions from coderabbit

* BUG: CodeRabbit

* ENH: Update env for cupy commit hooks

* BUG: Force commit of cupy methods until commit hooks updated

Author: aylward

docs/architecture.rst

@@ -55,8 +55,8 @@ Primary Workflows
 
 ``WorkflowFitStatisticalModelToPatient``
    Fits a template/statistical model to patient-specific surfaces with ICP,
-   optional PCA fitting, mask-to-mask registration, and optional image
-   refinement.
+   optional PCA fitting, labelmap-to-labelmap registration, and optional
+   labelmap-to-image refinement.
 
 ``WorkflowReconstructHighres4DCT``
    Reconstructs higher-resolution 4D CT frames from a time series and a fixed

docs/cli_scripts/fit_statistical_model_to_patient.rst

@@ -16,8 +16,8 @@ The registration pipeline consists of four stages:
 
 1. **ICP Alignment**: Rigid/affine alignment using surface matching
 2. **PCA Registration** (optional): Statistical shape model fitting
-3. **Mask-to-Mask Registration**: Greedy affine + ICON deformable registration using distance maps
-4. **Mask-to-Image Refinement** (optional): Final intensity-based refinement
+3. **Labelmap-to-Labelmap Registration**: Greedy affine + ICON deformable registration using distance maps
+4. **Labelmap-to-Image Refinement** (optional): Final intensity-based refinement
 
 Installation
 ============
@@ -84,7 +84,7 @@ Optional inputs:
 
 ``--template-labelmap PATH``
    Path to template labelmap image (.nii.gz, .nrrd, .mha). Required only when
-   ``--mask-to-image`` is set.
+   ``--labelmap-to-image`` is set.
 
 See :class:`physiomotion4d.WorkflowFitStatisticalModelToPatient` for API documentation.
 
@@ -112,16 +112,16 @@ PCA Registration Options
 Registration Configuration
 ---------------------------
 
-``--no-mask-to-mask``
-   Disable mask-to-mask deformable registration (default: enabled)
+``--no-labelmap-to-labelmap``
+   Disable labelmap-to-labelmap deformable registration (default: enabled)
 
-``--mask-to-image``
-   Enable mask-to-image refinement registration. Requires
+``--labelmap-to-image``
+   Enable labelmap-to-image refinement registration. Requires
    ``--template-labelmap`` and template label IDs. Disabled by default.
 
 ``--use-ICON-refinement``
-   Enable ICON deep learning refinement in the mask-to-image stage (Stage 4).
-   The mask-to-mask stage always uses Greedy affine + ICON deformable.
+   Enable ICON deep learning refinement in the labelmap-to-image stage (Stage 4).
+   The labelmap-to-labelmap stage always uses Greedy affine + ICON deformable.
    Default: disabled
 
 Output Options
@@ -180,7 +180,7 @@ Intermediate Results
 
 * ``{prefix}_icp_surface.vtp`` - Result after ICP alignment
 * ``{prefix}_pca_surface.vtp`` - Result after PCA fitting (if used)
-* ``{prefix}_m2m_surface.vtp`` - Result after mask-to-mask registration
+* ``{prefix}_l2l_surface.vtp`` - Result after labelmap-to-labelmap registration
 
 See Also
 ========

experiments/Heart-Create_Statistical_Model/3-registration_based_correspondence.py

@@ -107,7 +107,7 @@
         moving_model=moving_model,
         fixed_model=fixed_model,
         reference_image=reference_image,
-        roi_dilation_mm=20.0,  # Dilation for ROI mask
+        mask_dilation_mm=20.0,  # Dilation for binary registration mask
     )
 
     # Perform Greedy affine + ICON deformable registration
@@ -338,5 +338,5 @@
 # - The `RegisterModelsDistanceMaps` class uses a two-stage pipeline:
 #   1. **Greedy affine** registration (fast CPU-based alignment)
 #   2. **ICON deformable** registration on the affine-pre-aligned masks (deep learning)
-# - The `roi_dilation_mm` parameter controls the dilation of the ROI mask (default 20mm)
+# - The `mask_dilation_mm` parameter controls the dilation of the binary registration mask (default 20mm)
 # - Composed Greedy + ICON transforms provide smooth, invertible deformation fields for anatomical correspondence

experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient-CHOPValve.py

@@ -21,21 +21,21 @@
 
 # %%
 # Patient CT image (defines coordinate frame)
-patient_data_dir = Path.cwd().parent.parent / "data" / "CHOP-Valve4D" / "CT"
+patient_data_dir = Path(__file__).parent.parent.parent / "data" / "CHOP-Valve4D" / "CT"
 patient_ct_path = patient_data_dir / "RVOT28-Dias.mha"
 
 # Template model (moving)
-model_data_dir = Path.cwd().parent.parent / "data" / "KCL-Heart-Model"
+model_data_dir = Path(__file__).parent.parent.parent / "data" / "KCL-Heart-Model"
 model_labelmap_path = model_data_dir / "labelmap" / "average_labelmap_with_bkg.mha"
 model_pca_data_dir = (
-    Path.cwd().parent / "Heart-Create_Statistical_Model" / "kcl-heart-model"
+    Path(__file__).parent.parent / "Heart-Create_Statistical_Model" / "kcl-heart-model"
 )
 model_pca_json_path = model_pca_data_dir / "pca_model.json"
 model_mesh_path = model_pca_data_dir / "pca_mean.vtp"
 model_pca_n_modes = 10
 
 # Output directory
-output_dir = Path.cwd() / "results-chop"
+output_dir = Path(__file__).parent / "results-chop"
 
 # %%
 patient_image = itk.imread(str(patient_ct_path))
@@ -58,7 +58,7 @@
     True, pca_model=model_pca_data, pca_number_of_modes=model_pca_n_modes
 )
 
-registrar.set_use_mask_to_mask_registration(True)
+# registrar.set_use_labelmap_to_labelmap_registration(True)
 
 # %%
 patient_image = registrar.patient_image
@@ -75,10 +75,6 @@
 
 registered_model.save(str(output_dir / "registered_model.vtp"))
 registered_model_surface.save(str(output_dir / "registered_model_surface.vtp"))
-registered_labelmap = results["registered_template_labelmap"]
-itk.imwrite(
-    registered_labelmap, str(output_dir / "registered_labelmap.mha"), compression=True
-)
 
 # %%
 pca_model = registrar.pca_template_model

experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient.py

@@ -145,16 +145,16 @@
     registrar.set_use_pca_registration(
         True, pca_model=pca_model, pca_number_of_modes=pca_n_modes
     )
-    registrar.set_use_mask_to_image_registration(
+    registrar.set_use_labelmap_to_image_registration(
         True,
         template_labelmap=template_labelmap,
         template_labelmap_organ_mesh_ids=[1],
         template_labelmap_organ_extra_ids=[2, 3, 4, 5],
         template_labelmap_background_ids=[6],
     )
 
-    registrar.set_mask_dilation_mm(0)
-    registrar.set_roi_dilation_mm(25)
+    registrar.set_labelmap_dilation_mm(0)
+    registrar.set_mask_dilation_mm(25)
 
     patient_image = registrar.patient_image
     itk.imwrite(
@@ -184,37 +184,37 @@
     itk.imwrite(pca_labelmap, str(output_dir / "pca_labelmap.mha"), compression=True)
 
     # %% [markdown]
-    # ## Mask Alignment
+    # ## Labelmap Alignment
 
     # %%
     # Perform deformable registration
-    print("Starting deformable mask-to-mask registration...")
+    print("Starting deformable labelmap-to-labelmap registration...")
 
-    m2m_results = registrar.register_mask_to_mask()
-    m2m_inverse_transform = m2m_results["inverse_transform"]
-    m2m_forward_transform = m2m_results["forward_transform"]
-    m2m_model_surface = m2m_results["registered_template_model_surface"]
-    m2m_labelmap = m2m_results["registered_template_labelmap"]
+    l2l_results = registrar.register_labelmap_to_labelmap()
+    l2l_inverse_transform = l2l_results["inverse_transform"]
+    l2l_forward_transform = l2l_results["forward_transform"]
+    l2l_model_surface = l2l_results["registered_template_model_surface"]
+    l2l_labelmap = l2l_results["registered_template_labelmap"]
 
     print("Registration complete!")
 
-    m2m_model_surface.save(str(output_dir / "m2m_model_surface.vtp"))
-    itk.imwrite(m2m_labelmap, str(output_dir / "m2m_labelmap.mha"), compression=True)
+    l2l_model_surface.save(str(output_dir / "l2l_model_surface.vtp"))
+    itk.imwrite(l2l_labelmap, str(output_dir / "l2l_labelmap.mha"), compression=True)
 
     # %%
     print("Starting deformable registration...")
     print("This may take several minutes depending on GPU availability.")
 
-    m2i_results = registrar.register_labelmap_to_image()
-    m2i_inverse_transform = m2i_results["inverse_transform"]
-    m2i_forward_transform = m2i_results["forward_transform"]
-    m2i_surface = m2i_results["registered_template_model_surface"]
-    m2i_labelmap = m2i_results["registered_template_labelmap"]
+    l2i_results = registrar.register_labelmap_to_image()
+    l2i_inverse_transform = l2i_results["inverse_transform"]
+    l2i_forward_transform = l2i_results["forward_transform"]
+    l2i_surface = l2i_results["registered_template_model_surface"]
+    l2i_labelmap = l2i_results["registered_template_labelmap"]
     print("\nRegistration complete!")
 
     # Save registration results to output folder
-    m2i_surface.save(str(output_dir / "m2i_model_surface.vtp"))
-    itk.imwrite(m2i_labelmap, str(output_dir / "m2i_labelmap.mha"), compression=True)
+    l2i_surface.save(str(output_dir / "l2i_model_surface.vtp"))
+    itk.imwrite(l2i_labelmap, str(output_dir / "l2i_labelmap.mha"), compression=True)
 
     # %%
     tmp_p = itk.Point[itk.D, 3]()
@@ -241,16 +241,16 @@
     print(f"PCA transform time: {time.time() - start_time} seconds", flush=True)
 
     start_time = time.time()
-    tmp_p = registrar.m2m_inverse_transform.TransformPoint(tmp_p)
-    print(f"M2M inverse transform time: {time.time() - start_time} seconds", flush=True)
+    tmp_p = registrar.l2l_inverse_transform.TransformPoint(tmp_p)
+    print(f"L2L inverse transform time: {time.time() - start_time} seconds", flush=True)
 
     start_time = time.time()
-    tmp_p = registrar.m2i_inverse_transform.TransformPoint(tmp_p)
-    print(f"M2I inverse transform time: {time.time() - start_time} seconds", flush=True)
+    tmp_p = registrar.l2i_inverse_transform.TransformPoint(tmp_p)
+    print(f"L2I inverse transform time: {time.time() - start_time} seconds", flush=True)
 
     # %%
     # Verify registration using the transform member function
-    surface_transformed = registrar.m2i_template_model_surface
+    surface_transformed = registrar.l2i_template_model_surface
     surface_transformed.save(str(output_dir / "registered_template_surface.vtp"))
 
     model_transformed = registrar.transform_model()
@@ -265,8 +265,8 @@
     registered_surface = registrar.registered_template_model_surface
     icp_surface = registrar.icp_template_model_surface
     pca_surface = registrar.pca_template_model_surface
-    m2m_surface = registrar.m2m_template_model_surface
-    m2i_surface = registrar.m2i_template_model_surface
+    l2l_surface = registrar.l2l_template_model_surface
+    l2i_surface = registrar.l2i_template_model_surface
 
     # Create side-by-side comparison
     plotter = pv.Plotter(shape=(1, 2))
@@ -280,7 +280,7 @@
     # After deformable registration
     plotter.subplot(0, 1)
     plotter.add_mesh(patient_surface, color="red", opacity=0.5, label="Patient")
-    plotter.add_mesh(m2i_surface, color="blue", opacity=1.0, label="Registered")
+    plotter.add_mesh(l2i_surface, color="blue", opacity=1.0, label="Registered")
     plotter.add_title("Final Registration")
 
     plotter.link_views()

experiments/README.md

@@ -149,9 +149,9 @@ to handle diverse cases.
 
 **Technologies:**
 - ICP (Iterative Closest Point) registration for initial alignment
-- Mask-based deformable registration for anatomical fitting
+- Labelmap-based deformable registration for anatomical fitting
 - PCA (Principal Component Analysis) shape modeling for shape constraints
-- Three-stage registration pipeline (ICP → Mask-to-Mask → Mask-to-Image)
+- Three-stage registration pipeline (ICP → Labelmap-to-Labelmap → Labelmap-to-Image)
 - Computationally intensive (>1 hour on typical PC)
 
 **Prerequisites:**

pyproject.toml

@@ -221,6 +221,8 @@ show_error_codes = true
 module = [
     "ants",
     "cupy",
+    "cupy_backends",
+    "cupy_backends.*",
     "icon_registration",
     "icon_registration.*",
     "itk",
@@ -270,6 +272,12 @@ module = [
 ]
 disable_error_code = ["import-not-found", "import-untyped"]
 
+[[tool.mypy.overrides]]
+# torch/icon_registration/unigradicon are lazy-loaded at runtime; string
+# forward-reference annotations like "torch.Size" are intentional.
+module = ["physiomotion4d.register_images_icon"]
+ignore_errors = true
+
 [tool.pyright]
 # Third-party packages (e.g. pyvista) are in dependencies but may have no stubs;
 # do not report import-not-found so analysis matches mypy overrides above.
@@ -398,6 +406,8 @@ ignore = [
 "test_*.py" = ["S101", "PLR2004", "ARG001", "ARG002"]
 "tests/*.py" = ["S101", "PLR2004", "ARG001", "ARG002"]
 "experiments/**/*.py" = ["T201", "S101"]
+# torch/icon_registration are lazy-loaded at runtime; forward-reference strings are intentional
+"src/physiomotion4d/register_images_icon.py" = ["F821"]
 
 [tool.ruff.lint.isort]
 known-first-party = ["physiomotion4d"]

src/physiomotion4d/__init__.py

@@ -19,11 +19,10 @@
 
 __version__ = "2026.05.9"
 
+import importlib.util as _importlib_util
 import warnings as _warnings
 
-try:
-    import cupy as _cupy  # noqa: F401
-except ImportError:
+if _importlib_util.find_spec("cupy") is None:
     _warnings.warn(
         "CuPy is not installed — GPU acceleration is unavailable and processing "
         "will be slow. Re-install with uv to get CuPy and CUDA-enabled PyTorch "