diff --git a/.github/workflows/nightly-health.yml b/.github/workflows/nightly-health.yml index c96ee24..ce10a05 100644 --- a/.github/workflows/nightly-health.yml +++ b/.github/workflows/nightly-health.yml @@ -58,7 +58,7 @@ jobs: # Python 3.11 is required because the [physicsnemo] extra pulls in # nvidia-physicsnemo, which requires Python >= 3.11. run: | - & "C:\Users\saylward\AppData\Local\Programs\Python\Python311\python.exe" -m venv "$env:RUNNER_TEMP\physiomotion4d-venv" + & py -m venv "$env:RUNNER_TEMP\physiomotion4d-venv" echo "$env:RUNNER_TEMP\physiomotion4d-venv\Scripts" >> $env:GITHUB_PATH - name: Cache uv packages diff --git a/README.md b/README.md index df74e12..2d33e7f 100644 --- a/README.md +++ b/README.md @@ -159,7 +159,7 @@ print(WorkflowConvertImageToUSD.__name__) ### Key Dependencies -- **Medical Imaging**: ITK, TubeTK, MONAI, nibabel, PyVista +- **Medical Imaging**: ITK, MONAI, nibabel, PyVista - **AI/ML**: PyTorch, CuPy (CUDA 13), transformers, MONAI - **Registration**: icon-registration, unigradicon - **Visualization**: USD-core, PyVista diff --git a/docs/conf.py b/docs/conf.py index 9a7cfaf..91f9326 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -26,7 +26,6 @@ def __getattr__(cls, name): # Mock modules that need special handling -sys.modules["itk.TubeTK"] = Mock() sys.modules["icon_registration.losses"] = Mock() sys.modules["icon_registration.network_wrappers"] = Mock() diff --git a/docs/installation.rst b/docs/installation.rst index f63ce96..962f867 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -21,7 +21,7 @@ Software Dependencies PhysioMotion4D relies on several key packages: -* **Medical Imaging**: ITK, TubeTK, MONAI, nibabel, PyVista +* **Medical Imaging**: ITK, MONAI, nibabel, PyVista * **AI/ML**: PyTorch, CuPy (CUDA 13), transformers, MONAI * **Registration**: icon-registration, unigradicon * **Visualization**: USD-core, PyVista diff --git a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py index 92c2027..dc53964 100644 --- a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py +++ b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py @@ -22,7 +22,6 @@ import itk import numpy as np import pyvista as pv -from itk import TubeTK as ttk # Import from PhysioMotion4D package from physiomotion4d import ( @@ -30,6 +29,7 @@ RegisterModelsICPITK, TransformTools, ) +from physiomotion4d.image_tools import ImageTools from physiomotion4d.test_tools import TestTools # %% [markdown] @@ -64,10 +64,7 @@ # Resample to 1mm isotropic spacing print("Resampling to sotropic...") -resampler = ttk.ResampleImage.New(Input=patient_image) -resampler.SetMakeHighResIso(True) -resampler.Update() -patient_image = resampler.GetOutput() +patient_image = ImageTools().make_isotropic_image(patient_image) print(f" Resampled size: {itk.size(patient_image)}") print(f" Resampled spacing: {itk.spacing(patient_image)}") diff --git a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py index 18e4b1e..a428264 100644 --- a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py +++ b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py @@ -24,7 +24,6 @@ import itk import numpy as np import pyvista as pv -from itk import TubeTK as ttk # Import from PhysioMotion4D package from physiomotion4d import ( @@ -33,6 +32,7 @@ RegisterModelsPCA, TransformTools, ) +from physiomotion4d.image_tools import ImageTools from physiomotion4d.test_tools import TestTools # %% [markdown] @@ -75,10 +75,7 @@ # Resample to 1mm isotropic spacing print("Resampling to sotropic...") -resampler = ttk.ResampleImage.New(Input=patient_image) -resampler.SetMakeHighResIso(True) -resampler.Update() -patient_image = resampler.GetOutput() +patient_image = ImageTools().make_isotropic_image(patient_image) print(f" Resampled size: {itk.size(patient_image)}") print(f" Resampled spacing: {itk.spacing(patient_image)}") diff --git a/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py b/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py index ef3cbb5..fa4a641 100644 --- a/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py +++ b/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py @@ -6,8 +6,8 @@ import itk import numpy as np from data_dirlab_4d_ct import DataDirLab4DCT -from itk import TubeTK as tube +from physiomotion4d.image_tools import ImageTools from physiomotion4d.register_images_icon import RegisterImagesICON from physiomotion4d.segment_chest_total_segmentator import SegmentChestTotalSegmentator from physiomotion4d.transform_tools import TransformTools @@ -33,10 +33,7 @@ # %% def dilate_mask(mask: Optional[itk.image], dilation: int) -> Optional[itk.image]: if mask is not None: - im_math = tube.ImageMath.New(mask) - im_math.Dilate(dilation, 1, 0) - dilated_mask = im_math.GetOutputShort() - return dilated_mask + return ImageTools().binary_dilate_image(mask, dilation, 1, 0) return None def register_image( diff --git a/pyproject.toml b/pyproject.toml index 1c95126..0aa4b0a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,7 +55,6 @@ keywords = [ dependencies = [ # Core medical imaging "itk>=5.3.0", - "itk-tubetk>=1.4.0", "nibabel>=4.0.0", "numpy>=1.21.0", "pydicom>=2.4.0", diff --git a/src/physiomotion4d/image_tools.py b/src/physiomotion4d/image_tools.py index 2d880a0..3eaebf6 100644 --- a/src/physiomotion4d/image_tools.py +++ b/src/physiomotion4d/image_tools.py @@ -234,9 +234,6 @@ def convert_array_to_image_of_vectors( def make_isotropic_image(self, image: itk.Image) -> itk.Image: """Resample a 3-D *image* to isotropic spacing using the finest voxel pitch. - Equivalent to TubeTK's ResampleImage.SetMakeHighResIso(True), implemented - with standard ITK ResampleImageFilter so that TubeTK is not needed here. - Args: image: 3-D ITK image to resample. @@ -273,6 +270,118 @@ def make_isotropic_image(self, image: itk.Image) -> itk.Image: result.DisconnectPipeline() return result + def binary_dilate_image( + self, + image: itk.Image, + radius: int, + foreground_value: int = 1, + background_value: int = 0, + ) -> itk.Image: + """Binary-dilate *image* with a ball structuring element. + + Args: + image: Binary (or label) image to dilate. + radius: Radius, in voxels, of the ball structuring element. + foreground_value: Pixel value treated as foreground (default: 1). + background_value: Pixel value written for background voxels + (default: 0). + + Returns: + Dilated image with the same pixel type as *image*. + """ + ImageType = type(image) + dimension = image.GetImageDimension() + StructuringElementType = itk.FlatStructuringElement[dimension] + structuring_element = StructuringElementType.Ball(int(radius)) + + dilate_filter = itk.BinaryDilateImageFilter[ + ImageType, ImageType, StructuringElementType + ].New() + dilate_filter.SetInput(image) + dilate_filter.SetKernel(structuring_element) + dilate_filter.SetForegroundValue(foreground_value) + dilate_filter.SetBackgroundValue(background_value) + dilate_filter.Update() + result = dilate_filter.GetOutput() + result.DisconnectPipeline() + return result + + def binary_erode_image( + self, + image: itk.Image, + radius: int, + foreground_value: int = 1, + background_value: int = 0, + ) -> itk.Image: + """Binary-erode *image* with a ball structuring element. + + Args: + image: Binary (or label) image to erode. + radius: Radius, in voxels, of the ball structuring element. + foreground_value: Pixel value treated as foreground (default: 1). + background_value: Pixel value written for eroded-away voxels + (default: 0). + + Returns: + Eroded image with the same pixel type as *image*. + """ + ImageType = type(image) + dimension = image.GetImageDimension() + StructuringElementType = itk.FlatStructuringElement[dimension] + structuring_element = StructuringElementType.Ball(int(radius)) + + erode_filter = itk.BinaryErodeImageFilter[ + ImageType, ImageType, StructuringElementType + ].New() + erode_filter.SetInput(image) + erode_filter.SetKernel(structuring_element) + erode_filter.SetForegroundValue(foreground_value) + erode_filter.SetBackgroundValue(background_value) + erode_filter.Update() + result = erode_filter.GetOutput() + result.DisconnectPipeline() + return result + + def keep_largest_connected_component( + self, + image: itk.Image, + foreground_value: int = 1, + fully_connected: bool = False, + ) -> itk.Image: + """Keep only the largest connected component of a binary image. + + Args: + image: Binary (non-zero = foreground) image. + foreground_value: Value written for the retained component's + voxels (default: 1). + fully_connected: Whether diagonally-adjacent voxels are + considered connected (default: False, i.e. face connectivity + only). + + Returns: + Binary image, same pixel type as *image*, containing only the + largest connected component with value *foreground_value* + (background is 0). + """ + ImageType = type(image) + + connected_component_image = itk.connected_component_image_filter( + image, fully_connected=fully_connected + ) + relabeled_image = itk.relabel_component_image_filter( + connected_component_image, sort_by_object_size=True + ) + result = itk.binary_threshold_image_filter( + Input=relabeled_image, + LowerThreshold=1, + UpperThreshold=1, + InsideValue=foreground_value, + OutsideValue=0, + ) + if type(result) is not ImageType: + result = itk.cast_image_filter(result, ttype=(type(result), ImageType)) + return result + @overload def flip_image( self, diff --git a/src/physiomotion4d/register_images_icon.py b/src/physiomotion4d/register_images_icon.py index 3eb82b0..948aa3b 100644 --- a/src/physiomotion4d/register_images_icon.py +++ b/src/physiomotion4d/register_images_icon.py @@ -24,8 +24,7 @@ def _load_icon(): """Lazy-load icon_registration, torch, and unigradicon to avoid - initializing GPU/CUDA resources at import time, which interferes with - TubeTK's memory allocator on Windows.""" + initializing GPU/CUDA resources at import time.""" import icon_registration as icon import icon_registration.itk_wrapper import torch diff --git a/src/physiomotion4d/segment_anatomy_base.py b/src/physiomotion4d/segment_anatomy_base.py index 38497a1..3c0a564 100644 --- a/src/physiomotion4d/segment_anatomy_base.py +++ b/src/physiomotion4d/segment_anatomy_base.py @@ -10,9 +10,9 @@ import itk import numpy as np -from itk import TubeTK as tube from .anatomy_taxonomy import AnatomyTaxonomy +from .image_tools import ImageTools from .physiomotion4d_base import PhysioMotion4DBase @@ -427,10 +427,13 @@ def segment_connected_component( InsideValue=1, OutsideValue=0, ) - imMath = tube.ImageMath.New(connected_component_image) - imMath.Dilate(hole_fill, 1, 0) - imMath.Erode(hole_fill, 1, 0) - connected_component_image = imMath.GetOutputUChar() + image_tools = ImageTools() + connected_component_image = image_tools.binary_dilate_image( + connected_component_image, hole_fill, 1, 0 + ) + connected_component_image = image_tools.binary_erode_image( + connected_component_image, hole_fill, 1, 0 + ) labelmap_arr = itk.GetArrayFromImage(labelmap_image) connected_component_arr = itk.GetArrayFromImage(connected_component_image) diff --git a/src/physiomotion4d/segment_heart_simpleware.py b/src/physiomotion4d/segment_heart_simpleware.py index b7982e8..0673581 100644 --- a/src/physiomotion4d/segment_heart_simpleware.py +++ b/src/physiomotion4d/segment_heart_simpleware.py @@ -15,8 +15,8 @@ import itk import numpy as np -from itk import TubeTK as tube +from .image_tools import ImageTools from .segment_anatomy_base import SegmentAnatomyBase @@ -304,11 +304,14 @@ def segmentation_method(self, preprocessed_image: itk.image) -> itk.image: # Dilate the interior regions to simulate 3mm myocardium (heart) interior_image = itk.GetImageFromArray(interior_array.astype(np.uint8)) interior_image.CopyInformation(preprocessed_image) - imMath = tube.ImageMath.New(interior_image) + imMath = ImageTools() spacing = interior_image.GetSpacing() - imMath.Dilate(round(7 / spacing[0]), 1, 0) - imMath.Erode(round(4 / spacing[0]), 1, 0) - exterior_image = imMath.GetOutputUChar() + exterior_image = imMath.binary_dilate_image( + interior_image, round(7 / spacing[0]), 1, 0 + ) + exterior_image = imMath.binary_erode_image( + exterior_image, round(4 / spacing[0]), 1, 0 + ) exterior_array = itk.GetArrayFromImage(exterior_image) mask_id = 6 # Heart mask id exterior_array = exterior_array * mask_id diff --git a/src/physiomotion4d/segment_heart_simpleware_trimmed_branches.py b/src/physiomotion4d/segment_heart_simpleware_trimmed_branches.py index 9d8e6e2..1a4efc1 100644 --- a/src/physiomotion4d/segment_heart_simpleware_trimmed_branches.py +++ b/src/physiomotion4d/segment_heart_simpleware_trimmed_branches.py @@ -9,8 +9,8 @@ import itk import numpy as np -from itk import TubeTK as tube +from .image_tools import ImageTools from .segment_heart_simpleware import SegmentHeartSimpleware @@ -85,63 +85,47 @@ def trim_branches(self, labelmap_image: itk.image) -> itk.image: heart_arr[heart_arr == 6] = 0 heart_arr[heart_arr == 5] = 0 - img = itk.image_from_array(heart_arr) - img.CopyInformation(labelmap_image) - imMath = tube.ImageMath.New(img) - - # 2) Erode then Dilate Left Atrium label to clip vessels + image_tools = ImageTools() spacing = labelmap_image.GetSpacing() - imMath.Erode(round(7 / spacing[0]), 3, 0) - imMath.Dilate(round(7 / spacing[0]), 3, 0) + # 2) Erode then Dilate Left Atrium label to clip vessels # 3) Erode then Dilate Right Atrium label to clip vessels - imMath.Erode(round(7 / spacing[0]), 4, 0) - imMath.Dilate(round(7 / spacing[0]), 4, 0) - simple_img = imMath.GetOutput() - simple_arr = itk.array_from_image(simple_img) - - # Keep the largest component of the left atrium - simple_arr_3 = simple_arr.copy() - simple_arr_3[simple_arr_3 != 3] = 0 - simple_arr_3[simple_arr_3 == 3] = 1 - simple_img_3 = itk.image_from_array(simple_arr_3) - connComp = tube.SegmentConnectedComponents.New(simple_img_3) - connComp.SetKeepOnlyLargestComponent(True) - connComp.Update() - mask_img_3 = connComp.GetOutput() - mask_arr_3 = itk.array_from_image(mask_img_3) - simple_arr_3[mask_arr_3 == 0] = 0 - - # Keep the largest component of the right atrium - simple_arr_4 = simple_arr.copy() - simple_arr_4[simple_arr_4 != 4] = 0 - simple_arr_4[simple_arr_4 == 4] = 1 - simple_img_4 = itk.image_from_array(simple_arr_4) - connComp = tube.SegmentConnectedComponents.New(simple_img_4) - connComp.SetKeepOnlyLargestComponent(True) - connComp.Update() - mask_img_4 = connComp.GetOutput() - mask_arr_4 = itk.array_from_image(mask_img_4) - simple_arr_4[mask_arr_4 == 0] = 0 - - # Replace the left and right atrium labels with the largest components - simple_arr[simple_arr == 3] = 0 - simple_arr[simple_arr == 4] = 0 - simple_arr[simple_arr_3 > 0] = 3 - simple_arr[simple_arr_4 > 0] = 4 - simple_img = itk.image_from_array(simple_arr) - simple_img.CopyInformation(labelmap_image) + # + # Each label is isolated into its own binary mask before the + # open (erode-then-dilate) operation so that the opening of one + # label can never bleed into a neighboring label. + simple_arr = heart_arr.copy() + for label_id in (3, 4): + label_mask_arr = (heart_arr == label_id).astype(np.uint8) + label_mask_img = itk.image_from_array(label_mask_arr) + label_mask_img.CopyInformation(labelmap_image) + radius = round(7 / spacing[0]) + label_mask_img = image_tools.binary_erode_image( + label_mask_img, radius, 1, 0 + ) + label_mask_img = image_tools.binary_dilate_image( + label_mask_img, radius, 1, 0 + ) + + # Keep only the largest connected component of this label + label_mask_img = image_tools.keep_largest_connected_component( + label_mask_img, foreground_value=1 + ) + label_mask_arr = itk.array_from_image(label_mask_img) + + simple_arr[simple_arr == label_id] = 0 + simple_arr[label_mask_arr > 0] = label_id # 4) Dilate all others = keep_mask keep_mask_arr = heart_arr.copy() keep_mask_arr[keep_mask_arr == 2] = 1 keep_mask_arr[keep_mask_arr == 5] = 1 keep_mask_arr[keep_mask_arr != 1] = 0 - keep_mask = itk.image_from_array(keep_mask_arr) + keep_mask = itk.image_from_array(keep_mask_arr.astype(np.uint8)) keep_mask.CopyInformation(labelmap_image) - imMath.SetInput(keep_mask) - imMath.Dilate(round(7 / spacing[0]), 1, 0) - keep_mask = imMath.GetOutput() + keep_mask = image_tools.binary_dilate_image( + keep_mask, round(7 / spacing[0]), 1, 0 + ) keep_mask_arr = itk.array_from_image(keep_mask) # Add the left and right atrium labels to the keep_mask @@ -155,12 +139,14 @@ def trim_branches(self, labelmap_image: itk.image) -> itk.image: keep_mask_arr = heart_arr.copy() keep_mask_arr[keep_mask_arr == 1] = 0 keep_mask_arr[keep_mask_arr > 0] = 1 - keep_mask = itk.image_from_array(keep_mask_arr) + keep_mask = itk.image_from_array(keep_mask_arr.astype(np.uint8)) keep_mask.CopyInformation(labelmap_image) - imMath.SetInput(keep_mask) - imMath.Dilate(round(5 / spacing[0]), 1, 0) - imMath.Erode(round(2 / spacing[0]), 1, 0) - heart_mask = imMath.GetOutput() + keep_mask = image_tools.binary_dilate_image( + keep_mask, round(5 / spacing[0]), 1, 0 + ) + heart_mask = image_tools.binary_erode_image( + keep_mask, round(2 / spacing[0]), 1, 0 + ) # Insert the heart and myo labels back into the labelmap heart_mask_arr = itk.array_from_image(heart_mask) @@ -175,11 +161,9 @@ def trim_branches(self, labelmap_image: itk.image) -> itk.image: # Add in missing pieces / gaps of the myocardium lv_arr = heart_arr.copy() lv_arr[lv_arr != 1] = 0 - lv_img = itk.image_from_array(lv_arr) + lv_img = itk.image_from_array(lv_arr.astype(np.uint8)) lv_img.CopyInformation(labelmap_image) - imMath.SetInput(lv_img) - imMath.Dilate(round(2 / spacing[0]), 1, 0) - lv_img = imMath.GetOutput() + lv_img = image_tools.binary_dilate_image(lv_img, round(2 / spacing[0]), 1, 0) lv_arr = itk.array_from_image(lv_img) lv_arr = lv_arr * 5 # Myocardium label is 5 diff --git a/statistics.md b/statistics.md index e1c3412..db3ead9 100644 --- a/statistics.md +++ b/statistics.md @@ -105,7 +105,7 @@ PhysioMotion4D operates across several technically demanding domains: | Category | Key Packages | | --------------------- | ----------------------------------------------- | -| **Medical Imaging** | ITK, TubeTK, MONAI, nibabel, pynrrd | +| **Medical Imaging** | ITK, MONAI, nibabel, pynrrd | | **Deep Learning** | PyTorch, CuPy (CUDA 13), transformers | | **Registration** | ANTs, icon-registration, UniGradICON | | **3D Graphics / USD** | VTK, PyVista, USD-core | diff --git a/tests/test_segment_heart_simpleware_trimmed_branches.py b/tests/test_segment_heart_simpleware_trimmed_branches.py index 54bab88..99fa760 100644 --- a/tests/test_segment_heart_simpleware_trimmed_branches.py +++ b/tests/test_segment_heart_simpleware_trimmed_branches.py @@ -2,17 +2,11 @@ """ Tests for SegmentHeartSimplewareTrimmedBranches. -Structural tests (subclass identity, method placement) run everywhere, -including the Python 3.11 baseline. Tests that actually execute -trim_branches() are gated behind Python >= 3.12: trim_branches() uses -itk.TubeTK's SegmentConnectedComponents, whose native module segfaults -when loaded under CPython 3.11 (it is stable on 3.12+). The end-to-end -segmentation test additionally requires Simpleware Medical with ASCardio, -matching tests/test_segment_heart_simpleware.py's gating. +trim_branches() is implemented entirely with standard ITK filters, so all +tests below run on every supported Python version. """ import os -import sys import itk import numpy as np @@ -23,16 +17,6 @@ SegmentHeartSimplewareTrimmedBranches, ) -# itk.TubeTK's SegmentConnectedComponents (used by trim_branches) segfaults -# when its native module loads under CPython 3.11; it is stable on 3.12+. -# A skipped segfault-guard is the only option here - a C-level segfault -# cannot be caught with try/except, so any test that reaches TubeTK on 3.11 -# would crash the whole pytest process. -_tubetk_segfaults_on_this_python = sys.version_info < (3, 12) -_TUBETK_SKIP_REASON = ( - "itk.TubeTK SegmentConnectedComponents segfaults on CPython < 3.12" -) - def _simpleware_available(segmenter: SegmentHeartSimpleware) -> bool: """Return True if Simpleware Medical executable and script exist.""" @@ -61,7 +45,7 @@ def _synthetic_heart_labelmap() -> itk.Image: def test_subclass_owns_trim_branches() -> None: """SegmentHeartSimplewareTrimmedBranches is a SegmentHeartSimpleware that owns trim_branches(); the base class no longer defines trim_branches() or - set_trim_branches(). This runs on every supported Python (no TubeTK).""" + set_trim_branches().""" segmenter = SegmentHeartSimplewareTrimmedBranches() assert isinstance(segmenter, SegmentHeartSimpleware) assert hasattr(segmenter, "trim_branches") @@ -69,7 +53,6 @@ def test_subclass_owns_trim_branches() -> None: assert not hasattr(SegmentHeartSimpleware, "set_trim_branches") -@pytest.mark.skipif(_tubetk_segfaults_on_this_python, reason=_TUBETK_SKIP_REASON) def test_trim_branches_runs_on_synthetic_labelmap() -> None: """trim_branches() is pure post-processing - it must run without Simpleware and return a labelmap matching the input's geometry.""" @@ -89,7 +72,6 @@ def test_trim_branches_runs_on_synthetic_labelmap() -> None: @pytest.mark.requires_gpu @pytest.mark.requires_simpleware @pytest.mark.slow -@pytest.mark.skipif(_tubetk_segfaults_on_this_python, reason=_TUBETK_SKIP_REASON) def test_segment_trims_branches_relative_to_plain_segmenter( test_images: list, ) -> None: