From 202ab4af6fc6c7882b8a41dc527afc0edfb7847a Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Thu, 6 Nov 2025 20:42:43 -0500 Subject: [PATCH 1/6] Increment version to 0.28.0 --- pyproject.toml | 2 +- src/highdicom/version.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 80bd4fd2..7abfba2e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "highdicom" -version = "0.27.0" +version = "0.28.0" description = "High-level DICOM abstractions." readme = "README.md" requires-python = ">=3.10" diff --git a/src/highdicom/version.py b/src/highdicom/version.py index cf7b6d65..1bf36757 100644 --- a/src/highdicom/version.py +++ b/src/highdicom/version.py @@ -1 +1 @@ -__version__ = '0.27.0' +__version__ = '0.28.0' From 4b5a03877c349205891848cfd912369a5def8b18 Mon Sep 17 00:00:00 2001 From: "Mason C. Cleveland" <104479423+mccle@users.noreply.github.com> Date: Sat, 21 Feb 2026 13:19:57 -0500 Subject: [PATCH 2/6] Feature/optional dependencies (#387) Add import_optional_dependency util function --- pyproject.toml | 1 + src/highdicom/_dependency_utils.py | 64 ++++++++++++++++++++++++++++++ 2 files changed, 65 insertions(+) create mode 100644 src/highdicom/_dependency_utils.py diff --git a/pyproject.toml b/pyproject.toml index 74de21df..8eeacc4b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ dependencies = [ "pydicom>=3.0.1", "pyjpegls>=1.0.0", "typing-extensions>=4.0.0", + "packaging>=25.0" ] [project.optional-dependencies] diff --git a/src/highdicom/_dependency_utils.py b/src/highdicom/_dependency_utils.py new file mode 100644 index 00000000..9c7643c6 --- /dev/null +++ b/src/highdicom/_dependency_utils.py @@ -0,0 +1,64 @@ +from importlib import import_module, metadata +from packaging.requirements import Requirement +from types import ModuleType + + +def import_optional_dependency( + module_name: str, + feature: str +) -> ModuleType: + """Import an optional dependency. + + This function is designed to support interaction with other common + libraries that are not required for `highdicom` by default. + + Parameters + ---------- + module_name: str + Name of the module to be imported. + feature: str + Name or description of the feature that requires this dependency. + This is used for improving the clarity of error messages. + + Returns + ------- + ModuleType: + Imported module. + + Raises + ------ + ImportError: + When the specified module cannot be imported. + + """ + for req_str in metadata.requires('highdicom'): + req = Requirement(req_str) + if req.name == module_name: + break + + else: + raise ValueError( + f'`{module_name}` is not a requirement of highdicom ' + f'but is required for {feature}.' + ) + + try: + module = import_module(name=module_name) + + except ImportError as error: + raise ImportError( + f'Optional dependency `{module_name}` could not be imported' + f' but is required for {feature}.' + f' highdicom requires {module_name}{req.specifier}.' + ) from error + + installed_version = metadata.version(module_name) + + if installed_version not in req.specifier: + raise ImportError( + f'Optional dependency `{module_name}` has an unsuitable ' + f'version. Found {module_name}=={installed_version}, but ' + f'highdicom requires {module_name}{req.specifier}.' + ) + + return module From 41deabf7bdd09158f83f7a7a5f9c4bb7c6eface2 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Thu, 16 Apr 2026 20:36:21 +0000 Subject: [PATCH 3/6] Remove packaging version specifier (#399) * Remove packaging version specifier * Update packaging dependency version in pyproject.toml --------- Co-authored-by: Mason C. Cleveland <104479423+mccle@users.noreply.github.com> --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index ebb1ea49..6c259330 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,7 +39,7 @@ dependencies = [ "pydicom>=3.0.1", "pyjpegls>=1.0.0", "typing-extensions>=4.0.0", - "packaging>=25.0" + "packaging>=17.1" ] [project.optional-dependencies] From 58b6cda74e87d2243e4e6bfc5ebd56b636ea931d Mon Sep 17 00:00:00 2001 From: "Mason C. Cleveland" <104479423+mccle@users.noreply.github.com> Date: Fri, 15 May 2026 17:59:55 -0400 Subject: [PATCH 4/6] Feature/itk volume conversions (#394) * Add import_optional_dependency util function * Update optional import function. * Reformat code * Add module version compatibility check. * Move version incompatibility error raise location and update error message * Update src/highdicom/utils.py Co-authored-by: Chris Bridge * Update src/highdicom/utils.py Co-authored-by: Chris Bridge * Update src/highdicom/utils.py Co-authored-by: Chris Bridge * Update src/highdicom/utils.py Co-authored-by: Chris Bridge * Update src/highdicom/utils.py Co-authored-by: Chris Bridge * Move to new utils file * Add packaging to requirements * add itk and sitk conversions * Add frame_of_reference_uid and coordinate_system arguments; update docstrings; reformat strings * remove channel argument * add noqa to work with flake8 * update type hints, minor tweaks to bool retyping * sync changes * flake8 * Add tests for sitk * flake8 * skip sitk test on failed import * fix url syntax * Fix itk conversion and add tests * Update .github/workflows/run_unit_tests.yml Co-authored-by: Chris Bridge * Replace common url elements * Update docstrings and error messages * Update volume docs file * Update dtype tests with array equality and random values * Apply suggestion from @CPBridge Co-authored-by: Chris Bridge * Reorganize and update docs files * Update docstrings * Fix multichannel check and add tests --------- Co-authored-by: Chris Bridge --- .github/workflows/run_unit_tests.yml | 2 +- docs/general.rst | 1 + docs/itk_lib.rst | 49 ++ docs/other_libraries.rst | 24 + docs/sitk_lib.rst | 51 ++ docs/volume.rst | 60 +-- pyproject.toml | 2 + src/highdicom/volume.py | 309 ++++++++++++ tests/test_itk.py | 702 +++++++++++++++++++++++++++ tests/test_sitk.py | 650 +++++++++++++++++++++++++ 10 files changed, 1799 insertions(+), 51 deletions(-) create mode 100644 docs/itk_lib.rst create mode 100644 docs/other_libraries.rst create mode 100644 docs/sitk_lib.rst create mode 100644 tests/test_itk.py create mode 100644 tests/test_sitk.py diff --git a/.github/workflows/run_unit_tests.yml b/.github/workflows/run_unit_tests.yml index dea7ffbf..e22eb6fe 100644 --- a/.github/workflows/run_unit_tests.yml +++ b/.github/workflows/run_unit_tests.yml @@ -19,7 +19,7 @@ jobs: # and the likelihood of 3.12-specific bugs is considered low at this # stage python-version: ["3.10", "3.11", "3.13", "3.14"] - dependencies: [".", "'.[libjpeg]'"] + dependencies: [".", "'.[libjpeg,itk,sitk]'"] env: # Set this otherwise coverage on python 3.12 is absurdly slow diff --git a/docs/general.rst b/docs/general.rst index 13ad6fa3..9103f9a6 100644 --- a/docs/general.rst +++ b/docs/general.rst @@ -16,3 +16,4 @@ parts of the library. volume coding remote + other_libraries diff --git a/docs/itk_lib.rst b/docs/itk_lib.rst new file mode 100644 index 00000000..b6609fc8 --- /dev/null +++ b/docs/itk_lib.rst @@ -0,0 +1,49 @@ +.. _itk_lib: + +ITK +=== + +`ITK`_ is a widely-used library for volumetric image processing. In order to +use ITK with highdicom, the ``itk`` python package must be installed separately. +Version 5.4.0 or later is required. + +.. _itk_vol: + +Volume Conversions +------------------ + +Highdicom supports conversions with the ``itk.Image`` class through the +:meth:`highdicom.Volume.to_itk` and :meth:`highdicom.Volume.from_itk` methods. +Like highdicom, ITK uses the "LPS" convention. However, when converting to and +from NumPy arrays, ITK reverses the order of dimensions. This permutation is +handled automatically by highdicom and requires no intervention by the user. + +Creating an ITK Image from a Volume: + +.. code-block:: python + + import highdicom as hd + + + vol = hd.Volume(...) + + itk_im = vol.to_itk() + +Creating a volume from an ITK Image: + +.. code-block:: python + + import itk + import highdicom as hd + + + itk_im = itk.image(...) + + vol = hd.Volume.from_itk( + itk_im=itk_im, + coordinate_system='PATIENT', + frame_of_reference_uid=None + ) + + +.. _`ITK`: https://itk.org/ diff --git a/docs/other_libraries.rst b/docs/other_libraries.rst new file mode 100644 index 00000000..dcab858d --- /dev/null +++ b/docs/other_libraries.rst @@ -0,0 +1,24 @@ +.. _other_libraries: + +Interactions with Other Libraries +================================= + +Highdicom is able to interact with other imaging and scientific python +libraries. These integrations are intended to provide streamlined access to +complementary functionality available in external ecosystems that would +otherwise require custom user code. + +Because these libraries are not required for all workflows, they are +treated as optional dependencies and are not installed as part of the +default ``highdicom`` installation. Functionality requiring an external +library becomes available when an appropriate version of the corresponding +package is installed alongside ``highdicom``. The following sections describe +which libraries are currently supported and the additional features each +library enables. + +.. toctree:: + :maxdepth: 3 + :caption: Contents: + + itk_lib + sitk_lib diff --git a/docs/sitk_lib.rst b/docs/sitk_lib.rst new file mode 100644 index 00000000..d9321571 --- /dev/null +++ b/docs/sitk_lib.rst @@ -0,0 +1,51 @@ +.. _sitk_lib: + +SimpleITK +========= + +`SimpleITK`_ provides a simplified interface for the algorithms and data structures +of `ITK`_. In order to use ITK with highdicom, the ``SimpleITK`` python package must +be installed separately. Version 2.2.1 or later is required. + +.. _sitk_vol: + +Volume Conversions +------------------ + +Highdicom supports conversions with the ``SimpleITK.Image`` class through the +:meth:`highdicom.Volume.to_sitk` and :meth:`highdicom.Volume.from_sitk` methods. +Like highdicom, SimpleITK uses the "LPS" convention. However, when converting to +and from NumPy arrays, SimpleITK reverses the order of dimensions. This permutation +is handled automatically by highdicom and requires no intervention by the user. + + +Creating a SimpleITK Image from a Volume: + +.. code-block:: python + + import highdicom as hd + + + vol = hd.Volume(...) + + sitk_im = vol.to_sitk() + +Creating a volume from a SimpleITK Image: + +.. code-block:: python + + import SimpleITK as sitk + import highdicom as hd + + + sitk_im = sitk.image(...) + + vol = hd.Volume.from_sitk( + sitk_im=sitk_im, + coordinate_system='PATIENT', + frame_of_reference_uid=None + ) + + +.. _`ITK`: https://itk.org/ +.. _`SimpleITK`: https://simpleitk.org/ diff --git a/docs/volume.rst b/docs/volume.rst index ee9f99de..7bbff751 100644 --- a/docs/volume.rst +++ b/docs/volume.rst @@ -731,60 +731,20 @@ Writing a volume to a NIfTI file: nifti_path = '/path/to/nifti.nii' # or .nii.gz nib.save(nifti, nifti_path) -Volumes To/From ITK Images --------------------------- -`ITK`_ is a widely-used library for volumetric image processing. Its ``Image`` -class shares many similarities with our :class:`highdicom.Volume` class. Like -highdicom, ITK uses the "LPS" convention. However, when converting to and from -NumPy arrays, ITK reverses the order of dimensions. It is important to account -for this when performing conversions. +Volumes To/From Other Libraries +------------------------------- -We plan to add tools to handle this conversion in the near future, but for now -these snippets should correctly handle simple situations converting to and from -ITK Images. - -Creating a volume from an ITK Image: - -.. code-block:: python - - import itk - import numpy as np - import highdicom as hd - - - im = itk.image(...) - - # Reverse array dimension order - array = itk.array_from_image(im).transpose([2, 1, 0]) - - vol2 = hd.Volume.from_components( - array=array, - direction=np.asarray(im.GetDirection()), - spacing=np.asarray(im.GetSpacing()), - position=np.asarray(im.GetOrigin()), - coordinate_system="PATIENT" - ) - -Creating an ITK Image from a Volume: - -.. code-block:: python - - import itk - import highdicom as hd - - - vol = hd.Volume(...) - - # Reverse array dimension order - array = vol.array.transpose([2, 1, 0]) +Many other libraries have classes that share similarities with our +:class:`highdicom.Volume` class and allow for simple conversions between +formats. For commonly used libraries, highdicom Volumes provide +``to_()`` and ``from_()`` methods to simplify moving to and +from other libraries as much as possible. For more information about each library, +visit the following: - im = itk.image_from_array(array) - im.SetOrigin(vol.position) - im.SetDirection(vol.direction) - im.SetSpacing(vol.spacing) + - :ref:`ITK ` + - :ref:`SimpleITK ` .. _`NIfTI`: https://nifti.nimh.nih.gov/ -.. _`ITK`: https://itk.org/ .. _`nibabel`: https://nipy.org/nibabel/ diff --git a/pyproject.toml b/pyproject.toml index 6c259330..01ab8c05 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,6 +48,8 @@ libjpeg = [ "pylibjpeg-openjpeg>=2.0.0", "pylibjpeg>=2.0", ] +itk = ["itk>=5.4.0"] +sitk = ["SimpleITK>=2.2.1"] [dependency-groups] test = [ diff --git a/src/highdicom/volume.py b/src/highdicom/volume.py index d5c54d6a..0e9d949d 100644 --- a/src/highdicom/volume.py +++ b/src/highdicom/volume.py @@ -44,6 +44,8 @@ keyword_for_tag, ) +from highdicom._dependency_utils import import_optional_dependency +import warnings _DCM_PYTHON_TYPE_MAP = { 'CS': str, @@ -3583,6 +3585,313 @@ def pad_array(array: np.ndarray, cval: float) -> float: channels=self._channels, ) + def to_sitk(self) -> 'SimpleITK.Image': # noqa: F821 + """Convert the Volume to ``SimpleITK.Image`` format. + + The Volume is converted to a 3D ``SimpleITK.Image``. If + its array's current datatype is not supported by SimpleITK, + it is safely cast to a compatible type where possible. If + impossible to cast safely, a ``ValueError`` is raised. + Casting is performed on the following data types: + + - ``bool`` -> ``uint8`` + - ``float16`` -> ``float32`` (with warning) + - ``float128`` -> ``float64`` (with warning if possible, + else raises error) + + Spatial metadata (spacing, direction, origin) is preserved + with both highdicom and SimpletITK using "LPS" convention. + As SimpleITK uses column-major order and NumPy uses row-major, + this method automatically applies a transpose to the original + array. + + Returns + ------- + SimpleITK.Image: + Image constructed from the Volume. + + Raises + ------ + ValueError + When the volume is not 3D (multiple channels are unsupported). + ValueError + When the array's current datatype is not supported + and it is not possible to safely cast to a new datatype. + + """ + func = self.to_sitk + sitk = import_optional_dependency( + module_name='SimpleITK', + feature=f'{func.__module__}.{func.__qualname__}' + ) + + if self.array.ndim != 3: + raise ValueError( + 'SimpleITK conversion does not currently support' + ' volumes with multiple channels.' + ) + + array = self.array.transpose(2, 1, 0) + + if array.dtype == np.bool_: + array = array.astype(np.uint8) + + elif array.dtype == np.float16: + warnings.warn( + 'SimpleITK does not support float16 data.' + ' Safely casting to float32.' + ) + array = array.astype(np.float32) + + elif array.dtype == np.float128: + f64 = np.finfo(np.float64) + if array.min() >= f64.min and array.max() <= f64.max: + warnings.warn( + 'SimpleITK does not support float128 data.' + ' Casting to float64, precision may be lost.' + ) + array = array.astype(np.float64) + + else: + raise ValueError( + 'SimpleITK does not support float128 data.' + ' Casting to float64 is not possible.' + ) + + sitk_im = sitk.GetImageFromArray(array) + sitk_im.SetSpacing(self.spacing) + sitk_im.SetDirection(self.direction.flatten()) + sitk_im.SetOrigin(self.position) + + return sitk_im + + @classmethod + def from_sitk( + cls, + sitk_im: 'SimpleITK.Image', # noqa: F821 + coordinate_system: CoordinateSystemNames | str = 'PATIENT', + frame_of_reference_uid: str | None = None, + ) -> Self: + """Construct a Volume from a `SimpleITK.Image`. + + The ``SimpleITK.Image`` is converted to a 3D Volume. + Spatial metadata (spacing, direction, origin) is preserved + with both highdicom and SimpletITK using "LPS" convention. + As SimpleITK uses column-major order and NumPy uses row-major, + this method automatically applies a transpose to the original + array. + + Parameters + ---------- + sitk_im: SimpleITK.Image + A `SimpleITK.Image` to convert to a volume. + coordinate_system: highdicom.CoordinateSystemNames | str + Coordinate system (``"PATIENT"`` or ``"SLIDE"``) in which the volume + is defined. + frame_of_reference_uid: Union[str, None], optional + Frame of reference UID for the frame of reference, if known. + + Returns + ------- + highdicom.Volume: + Volume constructed from the `SimpleITK.Image`. + + Raises + ------ + ValueError + When the volume is not 3D (multiple channels are unsupported). + + """ + func = cls.from_sitk + sitk = import_optional_dependency( + module_name='SimpleITK', + feature=f'{func.__module__}.{func.__qualname__}' + ) + + array = sitk.GetArrayFromImage(sitk_im) + + if array.ndim != 3: + raise ValueError( + 'SimpleITK conversion does not currently support' + ' volumes with multiple channels.' + ) + + array = array.transpose(2, 1, 0) + + return cls.from_components( + array=array, + spacing=sitk_im.GetSpacing(), + coordinate_system=coordinate_system, + direction=np.reshape(sitk_im.GetDirection(), (3, 3)), + position=sitk_im.GetOrigin(), + frame_of_reference_uid=frame_of_reference_uid + ) + + def to_itk(self) -> 'itk.Image': # noqa: F821 + """Convert the volume to `itk.Image` format. + + The Volume is converted to a 3D ``itk.image``. If its array's + current datatype is not supported by ITK, it is safely cast to + a compatible type where possible. If impossible to cast safely, + a ``ValueError`` is raised. Casting is performed on the following + data types: + + - ``bool`` -> ``uint8`` + - ``int8`` -> ``int16`` (with warning) + - ``int64`` -> ``int32`` (with warning if possible, else + raises error) + - ``float16`` -> ``float32`` (with warning) + - ``float128`` -> ``float64`` (with warning if possible, else + raises error) + + Spatial metadata (spacing, direction, origin) is preserved + with both highdicom and ITK using "LPS" convention. As ITK uses + column-major order and NumPy uses row-major, this method automatically + applies a transpose to the original array. + + Returns + ------- + itk.Image: + Image constructed from the volume. + + Raises + ------ + ValueError + When the volume is not 3D (multiple channels are unsupported). + ValueError + When the array's current datatype is not supported + and it is not possible to safely cast to a new datatype. + + """ + func = self.to_itk + itk = import_optional_dependency( + module_name='itk', + feature=f'{func.__module__}.{func.__qualname__}' + ) + + if self.array.ndim != 3: + raise ValueError( + 'ITK conversion does not currently support' + ' volumes with multiple channels.' + ) + + array = self.array.transpose(2, 1, 0).copy() + + if array.dtype == np.bool_: + array = array.astype(np.uint8) + + elif array.dtype == np.int8: + warnings.warn( + 'ITK does not support int8 data.' + ' Safely casting to int16.' + ) + array = array.astype(np.int16) + + elif array.dtype == np.int64: + i32 = np.iinfo(np.int32) + if array.min() >= i32.min and array.max() <= i32.max: + warnings.warn( + 'ITK does not support int64 data.' + ' Safely casting to int32.' + ) + array = array.astype(np.int32) + + else: + raise ValueError( + 'ITK does not support int64 data.' + ' Safely casting to int32 is not possible.' + ) + + elif array.dtype == np.float16: + warnings.warn( + 'ITK does not support float16 data.' + ' Safely casting to float32.' + ) + array = array.astype(np.float32) + + elif array.dtype == np.float128: + f64 = np.finfo(np.float64) + if array.min() >= f64.min and array.max() <= f64.max: + warnings.warn( + 'ITK does not support float128 data.' + ' Casting to float64, precision may be lost.' + ) + array = array.astype(np.float64) + + else: + raise ValueError( + 'ITK does not support float128 data.' + ' Casting to float64 is not possible.' + ) + + itk_im = itk.GetImageFromArray(array) + itk_im.SetSpacing(self.spacing) + itk_im.SetDirection(self.direction) + itk_im.SetOrigin(self.position) + + return itk_im + + @classmethod + def from_itk( + cls, + itk_im: 'itk.Image', # noqa: F821 + coordinate_system: CoordinateSystemNames | str = 'PATIENT', + frame_of_reference_uid: str | None = None, + ) -> Self: + """Construct a Volume from an `itk.Image`. + + The ``itk.Image`` is converted to a 3D Volume. + Spatial metadata (spacing, direction, origin) is preserved + with both highdicom and ITK using "LPS" convention. As ITK uses + column-major order and NumPy uses row-major, this method automatically + applies a transpose to the original array. + + Parameters + ---------- + itk_im: itk.Image + A `itk.Image` to convert to a volume. + coordinate_system: highdicom.CoordinateSystemNames | str + Coordinate system (``"PATIENT"`` or ``"SLIDE"``) in which the volume + is defined. + frame_of_reference_uid: Union[str, None], optional + Frame of reference UID for the frame of reference, if known. + + Returns + ------- + highdicom.Volume: + Volume constructed from the `itk.Image`. + + Raises + ------ + ValueError + When the volume is not 3D (multiple channels are unsupported). + + """ + func = cls.from_itk + itk = import_optional_dependency( + module_name='itk', + feature=f'{func.__module__}.{func.__qualname__}' + ) + + array = itk.GetArrayFromImage(itk_im) + + if array.ndim != 3: + raise ValueError( + 'ITK conversion does not currently support' + ' volumes with multiple channels.' + ) + + array = array.transpose(2, 1, 0).copy() + + return cls.from_components( + array=array, + spacing=np.array(itk_im.GetSpacing()), + coordinate_system=coordinate_system, + direction=np.reshape(itk_im.GetDirection(), (3, 3)), + position=np.array(itk_im.GetOrigin()), + frame_of_reference_uid=frame_of_reference_uid + ) + class VolumeToVolumeTransformer: diff --git a/tests/test_itk.py b/tests/test_itk.py new file mode 100644 index 00000000..77da27d8 --- /dev/null +++ b/tests/test_itk.py @@ -0,0 +1,702 @@ +import numpy as np +import pydicom +import tempfile +import zipfile +import urllib +import pytest + +from pathlib import Path +from typing import Sequence +from pydicom.data import get_testdata_file +from highdicom import ( + Volume, + get_volume_from_series, + imread, +) +from highdicom.spatial import get_closest_patient_orientation +from highdicom._dependency_utils import import_optional_dependency + +try: + itk = import_optional_dependency('itk', feature='itk tests') + +except Exception: + pytest.skip("Optional dependency not available", allow_module_level=True) + +DCM_QA_MPRAGE = 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main' # noqa: E501 +DCM_QA_ME = 'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master' +DCM_QA_PDT2 = 'https://github.com/neurolabusc/dcm_qa_pdt2/raw/refs/heads/main' + + +def read_multiframe_ct_volume(): + im = imread(get_testdata_file('eCT_Supplemental.dcm')) + return im.get_volume() + + +def read_ct_series_volume(): + ct_files = [ + get_testdata_file('dicomdirtests/77654033/CT2/17136'), + get_testdata_file('dicomdirtests/77654033/CT2/17196'), + get_testdata_file('dicomdirtests/77654033/CT2/17166'), + ] + ct_series = [pydicom.dcmread(f) for f in ct_files] + return get_volume_from_series(ct_series) + + +def read_github_zip_volume(url: str): + with tempfile.TemporaryDirectory() as temp_dir: + zipfilename = Path(temp_dir) / Path(url).name + urllib.request.urlretrieve(url, zipfilename) + + with zipfile.ZipFile(zipfilename, 'r') as zf: + zf.extractall(temp_dir) + + series = [pydicom.dcmread(f) for f in Path(temp_dir).glob('**/*.dcm')] + + return get_volume_from_series(series), series + + +def read_github_series_volume(urls: Sequence[str]): + series = [] + with tempfile.TemporaryDirectory() as temp_dir: + for url in urls: + filename = Path(temp_dir) / Path(url).name + urllib.request.urlretrieve(url, filename) + + series.append(pydicom.dcmread(filename)) + + return get_volume_from_series(series), series + + +def read_github_itk(url: str): + with tempfile.TemporaryDirectory() as temp_dir: + filename = Path(temp_dir) / Path(url).name + urllib.request.urlretrieve(url, filename) + + itk_im = itk.imread(filename) + + return itk_im + + +@pytest.mark.parametrize( + 'vol', + [ + # testdata_files + read_multiframe_ct_volume(), + read_ct_series_volume(), + # different orientations + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RAF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RAH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RPF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RPH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LAF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LAH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LPF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LPH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('HLP'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('FPR'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('HRP'), + # isotropic + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[0.5, 0.5], + spacing_between_slices=0.5, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[2.0, 2.0], + spacing_between_slices=2.0, + coordinate_system='PATIENT' + ), + # anisotropic + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[0.5, 0.5], + spacing_between_slices=2.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[2.0, 0.5], + spacing_between_slices=0.5, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[0.5, 2.0], + spacing_between_slices=0.5, + coordinate_system='PATIENT' + ), + # non-square + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 32, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (64, 128, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # single-slice + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 1)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 1, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (1, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # random position offset + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[84.40363858, 105.04467386, 143.73326388], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[-21.03512292, 35.19549233, -184.42393696], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[-197.36060235, 86.22231644, -14.79874245], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # random orientation + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[ + -0.9267662161157189, + 0.32283606313442387, + -0.1920449348627007, + -0.3751482085550474, + -0.7693372329674889, + 0.5170919101937937 + ], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[ + -0.16411694786392106, + -0.7887835415036736, + 0.5923564400567902, + 0.9859501024515502, + -0.11222667947664411, + 0.12372375636644917 + ], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[ + -0.4161787354389772, + 0.440351232633788, + 0.7955413578729377, + -0.2771687765466385, + -0.8947097563255662, + 0.350245515665062 + ], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # entirely random + Volume.from_attributes( + array=np.random.randint(113, 257, (192, 249, 84)), + image_position=[156.03935104, -57.61106994, -108.37601079], + image_orientation=[ + -0.3056572521325831, + -0.9434667295206645, + -0.12823484123411946, + 0.9440777770580763, + -0.3177984972478506, + 0.08787073467366081 + ], + pixel_spacing=[3.34201481, 2.35548103], + spacing_between_slices=2.82618053, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(180, 214, (96, 121, 50)), + image_position=[93.43769804, -184.44672839, -153.64700033], + image_orientation=[ + -0.7034764816836532, + 0.05445328347346446, + -0.7086294374614612, + -0.24703727011677554, + -0.9536256538386763, + 0.1719613314498601 + ], + pixel_spacing=[2.84370598, 0.69898499], + spacing_between_slices=0.57265037, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(81, 214, (34, 123, 59)), + image_position=[-252.23051789, 146.90528128, 84.40363858], + image_orientation=[ + -0.7374895773326877, + 0.5716640328411087, + 0.35959610242811746, + -0.35709313761636013, + -0.7820074118766409, + 0.510831575802926 + ], + pixel_spacing=[3.52582689, 3.90364516], + spacing_between_slices=2.29457888, + coordinate_system='PATIENT' + ), + + ] +) +def test_roundtrip(vol: Volume): + itk_im = vol.to_itk() + + assert np.allclose(vol.position, itk_im.GetOrigin(), atol=1e-4) + assert np.allclose(vol.spacing, itk_im.GetSpacing(), atol=1e-4) + assert np.allclose( + vol.direction, + itk_im.GetDirection(), + atol=1e-4 + ) + assert ( + vol.array == itk.GetArrayFromImage(itk_im).transpose(2, 1, 0) + ).all() + + itk_roundtrip = Volume.from_itk(itk_im) + + assert np.allclose(vol.affine, itk_roundtrip.affine, atol=1e-4) + assert (vol.array == itk_roundtrip.array).all() + + +@pytest.mark.parametrize( + 'dtype', + [ + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.int8, + np.int16, + np.int32, + np.int64, + np.float16, + np.float32, + np.float64, + np.float128, + np.bool_ + ] +) +def test_dtype_itk(dtype): + rng = np.random.default_rng() + size = (10, 10, 10) + + volume = Volume.from_attributes( + array=np.zeros(size), + image_position=(0.0, 0.0, 0.0), + image_orientation=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing=(1.0, 1.0), + spacing_between_slices=1.0, + coordinate_system='PATIENT', + ) + + if dtype == np.bool_: + volume.array = np.round(rng.random(size=size)).astype(dtype) + + itk_roundtrip = Volume.from_itk(volume.to_itk()) + assert itk_roundtrip.dtype == np.uint8 + assert (itk_roundtrip.array == volume.array).all() + + elif dtype == np.int8: + i8 = np.iinfo(np.int8) + volume.array = rng.integers(i8.min, i8.max, size=size, dtype=dtype) + + with pytest.warns( + UserWarning, + match=( + 'ITK does not support int8 data.' + ' Safely casting to int16.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + assert itk_roundtrip.dtype == np.int16 + assert (itk_roundtrip.array == volume.array).all() + + elif dtype == np.int64: + i32 = np.iinfo(np.int32) + volume.array = rng.integers(i32.min, i32.max, size=size, dtype=dtype) + + with pytest.warns( + UserWarning, + match=( + 'ITK does not support int64 data.' + ' Safely casting to int32.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + assert itk_roundtrip.dtype == np.int32 + assert (itk_roundtrip.array == volume.array).all() + + volume.array[(0, 0, 0)] = np.int64(i32.max) + 1 + with pytest.raises( + ValueError, + match=( + 'ITK does not support int64 data.' + ' Safely casting to int32 is not possible.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + volume.array[(0, 0, 0)] = np.int64(i32.min) - 1 + with pytest.raises( + ValueError, + match=( + 'ITK does not support int64 data.' + ' Safely casting to int32 is not possible.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + elif dtype == np.float16: + f16 = np.finfo(np.float16) + volume.array = rng.uniform(f16.min, f16.max, size=size).astype(dtype) + + with pytest.warns( + UserWarning, + match=( + 'ITK does not support float16 data.' + ' Safely casting to float32.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + assert itk_roundtrip.dtype == np.float32 + assert (itk_roundtrip.array == volume.array).all() + + elif dtype == np.float128: + f64 = np.finfo(np.float64) + array = rng.random(size).astype(dtype) + volume.array = (2 * array - 1) * 0.9 * f64.max + + with pytest.warns( + UserWarning, + match=( + 'ITK does not support float128 data.' + ' Casting to float64, precision may be lost.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + assert itk_roundtrip.dtype == np.float64 + assert np.allclose(itk_roundtrip.array, volume.array, atol=1e-4) + + volume.array[(0, 0, 0)] = 1.1 * np.float128(f64.max) + with pytest.raises( + ValueError, + match=( + 'ITK does not support float128 data.' + ' Casting to float64 is not possible.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + volume.array[(0, 0, 0)] = 1.1 * np.float128(f64.min) + with pytest.raises( + ValueError, + match=( + 'ITK does not support float128 data.' + ' Casting to float64 is not possible.' + ) + ): + itk_roundtrip = Volume.from_itk(volume.to_itk()) + + elif np.issubdtype(dtype, np.integer): + ib = np.iinfo(dtype) + volume.array = rng.integers(ib.min, ib.max, size=size, dtype=dtype) + + itk_roundtrip = Volume.from_itk(volume.to_itk()) + assert itk_roundtrip.dtype == volume.dtype + assert (itk_roundtrip.array == volume.array).all() + + else: + fb = np.finfo(dtype) + array = rng.random(size).astype(dtype) + volume.array = (2 * array - 1) * 0.9 * fb.max + + itk_roundtrip = Volume.from_itk(volume.to_itk()) + assert itk_roundtrip.dtype == volume.dtype + assert (itk_roundtrip.array == volume.array).all() + + +@pytest.mark.parametrize( + 'zip_url,nifti_url', + [ + ( + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/In/2_t1_mp2rage_sag_p3_32.zip', # noqa: E501 + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/Ref/Si_2_t1_mp2rage_sag_p3_32_INV1.nii.gz' # noqa: E501 + ), + ( + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/In/3_t1_mp2rage_sag_p3_32.zip', # noqa: E501 + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/Ref/Si_3_t1_mp2rage_sag_p3_32_INV2.nii.gz' # noqa: E501 + ), + ( + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/In/4_t1_mp2rage_sag_p3_32.zip', # noqa: E501 + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/Ref/Si_4_t1_mp2rage_sag_p3_32_UNI_Images.nii.gz' # noqa: E501 + ), + ( + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/In/5_HCP_T1.zip', # noqa: E501 + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/Ref/Si_5_HCP_T1.nii.gz' # noqa: E501 + ), + ( + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/In/6_T1_mprage_ns_sag_p2.zip', # noqa: E501 + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/Ref/Si_6_T1_mprage_ns_sag_p2.nii.gz' # noqa: E501 + ), + ( + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/In/8_T1_memprage_rms.zip', # noqa: E501 + 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main/Ref/Si_8_T1_memprage_rms_RMS.nii.gz' # noqa: E501 + ), + ] +) +def test_nifti_equivalence_zip(zip_url: str, nifti_url: str): + vol, series = read_github_zip_volume(zip_url) + itk_im = read_github_itk(nifti_url) + + orientation = get_closest_patient_orientation( + np.reshape(itk_im.GetDirection(), (3, 3)) + ) + oriented_vol = vol.to_patient_orientation(orientation) + + assert np.allclose(oriented_vol.position, itk_im.GetOrigin(), atol=1e-4) + assert np.allclose(oriented_vol.spacing, itk_im.GetSpacing(), atol=1e-4) + assert np.allclose( + oriented_vol.direction, + itk_im.GetDirection(), + atol=1e-4 + ) + assert ( + oriented_vol.array == itk.GetArrayFromImage(itk_im).transpose(2, 1, 0) + ).all() + + +@pytest.mark.parametrize( + 'dcm_urls,nifti_url', + [ + ( + [ + f'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master/In/2_me_FieldMap_GRE/{i:04d}.dcm' # noqa: E501 + for i in range(1, 37) + ], + 'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master/Ref/me_FieldMap_GRE_2_e1.nii' # noqa: E501 + ), + ( + [ + f'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master/In/2_me_FieldMap_GRE/{i:04d}_e2.dcm' # noqa: E501 + for i in range(1, 37) + ], + 'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master/Ref/me_FieldMap_GRE_2_e2.nii' # noqa: E501 + ), + ( + [ + f'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master/In/3_me_FieldMap_GRE/{i:04d}_e2_ph.dcm' # noqa: E501 + for i in range(1, 37) + ], + 'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master/Ref/me_FieldMap_GRE_3_e2_ph.nii' # noqa: E501 + ), + ( + [ + f'https://github.com/neurolabusc/dcm_qa_pdt2/raw/refs/heads/main/In/Siemens/VE11/{i:04d}.dcm' # noqa: E501 + for i in range(1, 36) + ], + 'https://github.com/neurolabusc/dcm_qa_pdt2/raw/refs/heads/main/Ref/Siemens_pd+t2_tse_sag_ISO_1.8mm_3_e1.nii' # noqa: E501 + ), + ( + [ + f'https://github.com/neurolabusc/dcm_qa_pdt2/raw/refs/heads/main/In/Siemens/VE11/{i:04d}_e2.dcm' # noqa: E501 + for i in range(36, 71) + ], + 'https://github.com/neurolabusc/dcm_qa_pdt2/raw/refs/heads/main/Ref/Siemens_pd+t2_tse_sag_ISO_1.8mm_3_e2.nii' # noqa: E501 + ), + ] +) +def test_nifti_equivalence_series(dcm_urls: Sequence[str], nifti_url: str): + vol, series = read_github_series_volume(dcm_urls) + itk_im = read_github_itk(nifti_url) + + orientation = get_closest_patient_orientation( + np.reshape(itk_im.GetDirection(), (3, 3)) + ) + oriented_vol = vol.to_patient_orientation(orientation) + + assert np.allclose(oriented_vol.position, itk_im.GetOrigin(), atol=1e-4) + assert np.allclose(oriented_vol.spacing, itk_im.GetSpacing(), atol=1e-4) + assert np.allclose( + oriented_vol.direction, + itk_im.GetDirection(), + atol=1e-4 + ) + assert ( + oriented_vol.array == itk.GetArrayFromImage(itk_im).transpose(2, 1, 0) + ).all() + + +def test_multichannel_volume(): + array = np.zeros((10, 10, 10, 2)) + volume = Volume.from_attributes( + array=array, + image_position=(0.0, 0.0, 0.0), + image_orientation=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing=(1.0, 1.0), + spacing_between_slices=2.0, + channels={'OpticalPathIdentifier': ['path1', 'path2']}, + coordinate_system="PATIENT", + ) + + with pytest.raises( + ValueError, + match=( + 'ITK conversion does not currently support' + ' volumes with multiple channels.' + ) + ): + volume.to_itk() + + array = np.zeros((10, 10, 10, 2)) + itk_im = itk.GetImageFromArray(array) + + with pytest.raises( + ValueError, + match=( + 'ITK conversion does not currently support' + ' volumes with multiple channels.' + ) + ): + Volume.from_itk(itk_im) diff --git a/tests/test_sitk.py b/tests/test_sitk.py new file mode 100644 index 00000000..1d6d605d --- /dev/null +++ b/tests/test_sitk.py @@ -0,0 +1,650 @@ +import numpy as np +import pydicom +import tempfile +import zipfile +import urllib +import pytest + +from pathlib import Path +from typing import Sequence +from pydicom.data import get_testdata_file +from highdicom import ( + Volume, + get_volume_from_series, + imread, +) +from highdicom.spatial import get_closest_patient_orientation +from highdicom._dependency_utils import import_optional_dependency + +try: + sitk = import_optional_dependency('SimpleITK', feature='sitk tests') + +except Exception: + pytest.skip("Optional dependency not available", allow_module_level=True) + +DCM_QA_MPRAGE = 'https://github.com/neurolabusc/dcm_qa_mprage/raw/refs/heads/main' # noqa: E501 +DCM_QA_ME = 'https://github.com/neurolabusc/dcm_qa_me/raw/refs/heads/master' +DCM_QA_PDT2 = 'https://github.com/neurolabusc/dcm_qa_pdt2/raw/refs/heads/main' + + +def read_multiframe_ct_volume(): + im = imread(get_testdata_file('eCT_Supplemental.dcm')) + return im.get_volume() + + +def read_ct_series_volume(): + ct_files = [ + get_testdata_file('dicomdirtests/77654033/CT2/17136'), + get_testdata_file('dicomdirtests/77654033/CT2/17196'), + get_testdata_file('dicomdirtests/77654033/CT2/17166'), + ] + ct_series = [pydicom.dcmread(f) for f in ct_files] + return get_volume_from_series(ct_series) + + +def read_github_zip_volume(url: str): + with tempfile.TemporaryDirectory() as temp_dir: + zipfilename = Path(temp_dir) / Path(url).name + urllib.request.urlretrieve(url, zipfilename) + + with zipfile.ZipFile(zipfilename, 'r') as zf: + zf.extractall(temp_dir) + + series = [pydicom.dcmread(f) for f in Path(temp_dir).glob('**/*.dcm')] + + return get_volume_from_series(series), series + + +def read_github_series_volume(urls: Sequence[str]): + series = [] + with tempfile.TemporaryDirectory() as temp_dir: + for url in urls: + filename = Path(temp_dir) / Path(url).name + urllib.request.urlretrieve(url, filename) + + series.append(pydicom.dcmread(filename)) + + return get_volume_from_series(series), series + + +def read_github_sitk(url: str): + with tempfile.TemporaryDirectory() as temp_dir: + filename = Path(temp_dir) / Path(url).name + urllib.request.urlretrieve(url, filename) + + sitk_im = sitk.ReadImage(filename) + + return sitk_im + + +@pytest.mark.parametrize( + 'vol', + [ + # testdata_files + read_multiframe_ct_volume(), + read_ct_series_volume(), + # different orientations + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RAF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RAH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RPF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('RPH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LAF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LAH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LPF'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('LPH'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('HLP'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('FPR'), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ).to_patient_orientation('HRP'), + # isotropic + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[0.5, 0.5], + spacing_between_slices=0.5, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[2.0, 2.0], + spacing_between_slices=2.0, + coordinate_system='PATIENT' + ), + # anisotropic + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[0.5, 0.5], + spacing_between_slices=2.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[2.0, 0.5], + spacing_between_slices=0.5, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[0.5, 2.0], + spacing_between_slices=0.5, + coordinate_system='PATIENT' + ), + # non-square + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 32, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (64, 128, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # single-slice + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 1)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 1, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (1, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # random position offset + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[84.40363858, 105.04467386, 143.73326388], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[-21.03512292, 35.19549233, -184.42393696], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[-197.36060235, 86.22231644, -14.79874245], + image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # random orientation + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[ + -0.9267662161157189, + 0.32283606313442387, + -0.1920449348627007, + -0.3751482085550474, + -0.7693372329674889, + 0.5170919101937937 + ], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[ + -0.16411694786392106, + -0.7887835415036736, + 0.5923564400567902, + 0.9859501024515502, + -0.11222667947664411, + 0.12372375636644917 + ], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(0, 100, (16, 16, 16)), + image_position=[0.0, 0.0, 0.0], + image_orientation=[ + -0.4161787354389772, + 0.440351232633788, + 0.7955413578729377, + -0.2771687765466385, + -0.8947097563255662, + 0.350245515665062 + ], + pixel_spacing=[1.0, 1.0], + spacing_between_slices=1.0, + coordinate_system='PATIENT' + ), + # entirely random + Volume.from_attributes( + array=np.random.randint(113, 257, (192, 249, 84)), + image_position=[156.03935104, -57.61106994, -108.37601079], + image_orientation=[ + -0.3056572521325831, + -0.9434667295206645, + -0.12823484123411946, + 0.9440777770580763, + -0.3177984972478506, + 0.08787073467366081 + ], + pixel_spacing=[3.34201481, 2.35548103], + spacing_between_slices=2.82618053, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(180, 214, (96, 121, 50)), + image_position=[93.43769804, -184.44672839, -153.64700033], + image_orientation=[ + -0.7034764816836532, + 0.05445328347346446, + -0.7086294374614612, + -0.24703727011677554, + -0.9536256538386763, + 0.1719613314498601 + ], + pixel_spacing=[2.84370598, 0.69898499], + spacing_between_slices=0.57265037, + coordinate_system='PATIENT' + ), + Volume.from_attributes( + array=np.random.randint(81, 214, (34, 123, 59)), + image_position=[-252.23051789, 146.90528128, 84.40363858], + image_orientation=[ + -0.7374895773326877, + 0.5716640328411087, + 0.35959610242811746, + -0.35709313761636013, + -0.7820074118766409, + 0.510831575802926 + ], + pixel_spacing=[3.52582689, 3.90364516], + spacing_between_slices=2.29457888, + coordinate_system='PATIENT' + ), + + ] +) +def test_roundtrip(vol: Volume): + sitk_im = vol.to_sitk() + + assert np.allclose(vol.position, sitk_im.GetOrigin(), atol=1e-4) + assert np.allclose(vol.spacing, sitk_im.GetSpacing(), atol=1e-4) + assert np.allclose( + vol.direction.flatten(), + sitk_im.GetDirection(), + atol=1e-4 + ) + assert ( + vol.array == sitk.GetArrayFromImage(sitk_im).transpose(2, 1, 0) + ).all() + + sitk_roundtrip = Volume.from_sitk(sitk_im) + + assert np.allclose(vol.affine, sitk_roundtrip.affine, atol=1e-4) + assert (vol.array == sitk_roundtrip.array).all() + + +@pytest.mark.parametrize( + 'dtype', + [ + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.int8, + np.int16, + np.int32, + np.int64, + np.float16, + np.float32, + np.float64, + np.float128, + np.bool_ + ] +) +def test_dtype_sitk(dtype: np.dtype): + rng = np.random.default_rng() + size = (10, 10, 10) + + volume = Volume.from_attributes( + array=np.zeros(size), + image_position=(0.0, 0.0, 0.0), + image_orientation=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing=(1.0, 1.0), + spacing_between_slices=1.0, + coordinate_system='PATIENT', + ) + + if dtype == np.bool_: + volume.array = np.round(rng.random(size=size)).astype(dtype) + + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + assert sitk_roundtrip.dtype == np.uint8 + assert (sitk_roundtrip.array == volume.array).all() + + elif dtype == np.float16: + f16 = np.finfo(np.float16) + volume.array = rng.uniform(f16.min, f16.max, size=size).astype(dtype) + + with pytest.warns( + UserWarning, + match=( + 'SimpleITK does not support float16 data.' + ' Safely casting to float32.' + ) + ): + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + + assert sitk_roundtrip.dtype == np.float32 + assert (sitk_roundtrip.array == volume.array).all() + + elif dtype == np.float128: + f64 = np.finfo(np.float64) + array = rng.random(size).astype(dtype) + volume.array = (2 * array - 1) * 0.9 * f64.max + + with pytest.warns( + UserWarning, + match=( + 'SimpleITK does not support float128 data.' + ' Casting to float64, precision may be lost.' + ) + ): + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + + assert sitk_roundtrip.dtype == np.float64 + assert np.allclose(sitk_roundtrip.array, volume.array, atol=1e-4) + + volume.array[(0, 0, 0)] = 1.1 * np.float128(np.finfo(np.float64).max) + with pytest.raises( + ValueError, + match=( + 'SimpleITK does not support float128 data.' + ' Casting to float64 is not possible.' + ) + ): + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + + volume.array[(0, 0, 0)] = 1.1 * np.float128(np.finfo(np.float64).min) + with pytest.raises( + ValueError, + match=( + 'SimpleITK does not support float128 data.' + ' Casting to float64 is not possible.' + ) + ): + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + + elif np.issubdtype(dtype, np.integer): + ib = np.iinfo(dtype) + volume.array = rng.integers(ib.min, ib.max, size=size, dtype=dtype) + + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + assert sitk_roundtrip.dtype == volume.dtype + assert (sitk_roundtrip.array == volume.array).all() + + else: + fb = np.finfo(dtype) + array = rng.random(size).astype(dtype) + volume.array = (2 * array - 1) * 0.9 * fb.max + + sitk_roundtrip = Volume.from_sitk(volume.to_sitk()) + assert sitk_roundtrip.dtype == volume.dtype + assert (sitk_roundtrip.array == volume.array).all() + + +@pytest.mark.parametrize( + 'zip_url,nifti_url', + [ + ( + f'{DCM_QA_MPRAGE}/In/2_t1_mp2rage_sag_p3_32.zip', + f'{DCM_QA_MPRAGE}/Ref/Si_2_t1_mp2rage_sag_p3_32_INV1.nii.gz' + ), + ( + f'{DCM_QA_MPRAGE}/In/3_t1_mp2rage_sag_p3_32.zip', + f'{DCM_QA_MPRAGE}/Ref/Si_3_t1_mp2rage_sag_p3_32_INV2.nii.gz' + ), + ( + f'{DCM_QA_MPRAGE}/In/4_t1_mp2rage_sag_p3_32.zip', + f'{DCM_QA_MPRAGE}/Ref/Si_4_t1_mp2rage_sag_p3_32_UNI_Images.nii.gz' + ), + ( + f'{DCM_QA_MPRAGE}/In/5_HCP_T1.zip', + f'{DCM_QA_MPRAGE}/Ref/Si_5_HCP_T1.nii.gz' + ), + ( + f'{DCM_QA_MPRAGE}/In/6_T1_mprage_ns_sag_p2.zip', + f'{DCM_QA_MPRAGE}/Ref/Si_6_T1_mprage_ns_sag_p2.nii.gz' + ), + ( + f'{DCM_QA_MPRAGE}/In/8_T1_memprage_rms.zip', + f'{DCM_QA_MPRAGE}/Ref/Si_8_T1_memprage_rms_RMS.nii.gz' + ), + ] +) +def test_nifti_equivalence_zip(zip_url: str, nifti_url: str): + vol, series = read_github_zip_volume(zip_url) + sitk_im = read_github_sitk(nifti_url) + + orientation = get_closest_patient_orientation( + np.reshape(sitk_im.GetDirection(), (3, 3)) + ) + oriented_vol = vol.to_patient_orientation(orientation) + + assert np.allclose(oriented_vol.position, sitk_im.GetOrigin(), atol=1e-4) + assert np.allclose(oriented_vol.spacing, sitk_im.GetSpacing(), atol=1e-4) + assert np.allclose( + oriented_vol.direction.flatten(), + sitk_im.GetDirection(), + atol=1e-4 + ) + assert ( + oriented_vol.array == sitk.GetArrayFromImage(sitk_im).transpose(2, 1, 0) + ).all() + + +@pytest.mark.parametrize( + 'dcm_urls,nifti_url', + [ + ( + [ + f'{DCM_QA_ME}/In/2_me_FieldMap_GRE/{i:04d}.dcm' + for i in range(1, 37) + ], + f'{DCM_QA_ME}/Ref/me_FieldMap_GRE_2_e1.nii' + ), + ( + [ + f'{DCM_QA_ME}/In/2_me_FieldMap_GRE/{i:04d}_e2.dcm' + for i in range(1, 37) + ], + f'{DCM_QA_ME}/Ref/me_FieldMap_GRE_2_e2.nii' + ), + ( + [ + f'{DCM_QA_ME}/In/3_me_FieldMap_GRE/{i:04d}_e2_ph.dcm' + for i in range(1, 37) + ], + f'{DCM_QA_ME}/Ref/me_FieldMap_GRE_3_e2_ph.nii' + ), + ( + [ + f'{DCM_QA_PDT2}/In/Siemens/VE11/{i:04d}.dcm' + for i in range(1, 36) + ], + f'{DCM_QA_PDT2}/Ref/Siemens_pd+t2_tse_sag_ISO_1.8mm_3_e1.nii' + ), + ( + [ + f'{DCM_QA_PDT2}/In/Siemens/VE11/{i:04d}_e2.dcm' + for i in range(36, 71) + ], + f'{DCM_QA_PDT2}/Ref/Siemens_pd+t2_tse_sag_ISO_1.8mm_3_e2.nii' + ), + ] +) +def test_nifti_equivalence_series(dcm_urls: Sequence[str], nifti_url: str): + vol, series = read_github_series_volume(dcm_urls) + sitk_im = read_github_sitk(nifti_url) + + orientation = get_closest_patient_orientation( + np.reshape(sitk_im.GetDirection(), (3, 3)) + ) + oriented_vol = vol.to_patient_orientation(orientation) + + assert np.allclose(oriented_vol.position, sitk_im.GetOrigin(), atol=1e-4) + assert np.allclose(oriented_vol.spacing, sitk_im.GetSpacing(), atol=1e-4) + assert np.allclose( + oriented_vol.direction.flatten(), + sitk_im.GetDirection(), + atol=1e-4 + ) + assert ( + oriented_vol.array == sitk.GetArrayFromImage(sitk_im).transpose(2, 1, 0) + ).all() + + +def test_multichannel_volume(): + array = np.zeros((10, 10, 10, 2)) + volume = Volume.from_attributes( + array=array, + image_position=(0.0, 0.0, 0.0), + image_orientation=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing=(1.0, 1.0), + spacing_between_slices=2.0, + channels={'OpticalPathIdentifier': ['path1', 'path2']}, + coordinate_system="PATIENT", + ) + + with pytest.raises( + ValueError, + match=( + 'SimpleITK conversion does not currently support' + ' volumes with multiple channels.' + ) + ): + volume.to_sitk() + + array = np.zeros((10, 10, 10, 2)) + sitk_im = sitk.GetImageFromArray(array) + + with pytest.raises( + ValueError, + match=( + 'SimpleITK conversion does not currently support' + ' volumes with multiple channels.' + ) + ): + Volume.from_sitk(sitk_im) From 522b13e959d7cb661fd8803ad09102e33b7669b5 Mon Sep 17 00:00:00 2001 From: Mason Cleveland Date: Mon, 18 May 2026 11:57:32 -0400 Subject: [PATCH 5/6] Enforce volume datatypes must be a subtype of float, integer, or bool --- src/highdicom/volume.py | 19 +++++++++++ tests/test_volume.py | 74 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) diff --git a/src/highdicom/volume.py b/src/highdicom/volume.py index 0e9d949d..c7fd9f27 100644 --- a/src/highdicom/volume.py +++ b/src/highdicom/volume.py @@ -2504,6 +2504,15 @@ def __init__( "Argument 'array' must be at least three dimensional." ) + if not any([ + np.issubdtype(array.dtype, dtype) for dtype in + [np.floating, np.integer, np.bool_] + ]): + raise ValueError( + "Argument 'array' must have a dtype of float, integer," + f" or bool, received '{array.dtype}'." + ) + if channels is None: channels = {} @@ -2856,6 +2865,16 @@ def array(self, value: np.ndarray) -> None: raise ValueError( "Array must match the shape of the existing array." ) + + if not any([ + np.issubdtype(value.dtype, dtype) for dtype in + [np.floating, np.integer, np.bool_] + ]): + raise ValueError( + "Array must have a dtype of float, integer," + f" or bool, received '{value.dtype}'." + ) + self._array = value def astype(self, dtype: type) -> Self: diff --git a/tests/test_volume.py b/tests/test_volume.py index 0b9d17e8..55e633c0 100644 --- a/tests/test_volume.py +++ b/tests/test_volume.py @@ -1341,3 +1341,77 @@ def test_match_geometry_segmentation(): seg_vol_from_matched_dataset.array, seg_volume_matched.array ) + + +@pytest.mark.parametrize( + 'dtype', + [ + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.int8, + np.int16, + np.int32, + np.int64, + int, + np.float16, + np.float32, + np.float64, + float, + np.float128, + np.complex64, + complex, + np.complex128, + np.complex256, + np.bool_, + bool, + str + ] +) +def test_dtype_restriction(dtype): + if any([ + np.issubdtype(dtype, t) for t in + [np.floating, np.integer, np.bool_] + ]): + volume = Volume.from_attributes( + array=np.zeros((10, 10, 10), dtype=dtype), + image_position=(0.0, 0.0, 0.0), + image_orientation=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing=(1.0, 1.0), + spacing_between_slices=1.0, + coordinate_system='PATIENT', + ) + + for disallowed_dtype in [ + np.complex64, + complex, + np.complex128, + np.complex256, + str + ]: + with pytest.raises( + ValueError, + match=( + "Array must have a dtype of float, integer," + " or bool, received" + ) + ): + volume.array = np.zeros((10, 10, 10), dtype=disallowed_dtype) + + else: + with pytest.raises( + ValueError, + match=( + "Argument 'array' must have a dtype of float, integer," + " or bool, received" + ) + ): + volume = Volume.from_attributes( + array=np.zeros((10, 10, 10), dtype=dtype), + image_position=(0.0, 0.0, 0.0), + image_orientation=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing=(1.0, 1.0), + spacing_between_slices=1.0, + coordinate_system='PATIENT', + ) From 98e07679cd5fd9dc0a0ed57aeaae0424bc354dfc Mon Sep 17 00:00:00 2001 From: Mason Cleveland Date: Mon, 18 May 2026 12:46:16 -0400 Subject: [PATCH 6/6] Update docstrings to indicate optional dependency versions --- src/highdicom/volume.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/highdicom/volume.py b/src/highdicom/volume.py index c7fd9f27..158d3998 100644 --- a/src/highdicom/volume.py +++ b/src/highdicom/volume.py @@ -3607,6 +3607,9 @@ def pad_array(array: np.ndarray, cval: float) -> float: def to_sitk(self) -> 'SimpleITK.Image': # noqa: F821 """Convert the Volume to ``SimpleITK.Image`` format. + This method requires an optional dependency to be installed + separately from highdicom, specifically ``SimpleITK>=2.2.1``. + The Volume is converted to a 3D ``SimpleITK.Image``. If its array's current datatype is not supported by SimpleITK, it is safely cast to a compatible type where possible. If @@ -3693,6 +3696,9 @@ def from_sitk( ) -> Self: """Construct a Volume from a `SimpleITK.Image`. + This method requires an optional dependency to be installed + separately from highdicom, specifically ``SimpleITK>=2.2.1``. + The ``SimpleITK.Image`` is converted to a 3D Volume. Spatial metadata (spacing, direction, origin) is preserved with both highdicom and SimpletITK using "LPS" convention. @@ -3749,6 +3755,9 @@ def from_sitk( def to_itk(self) -> 'itk.Image': # noqa: F821 """Convert the volume to `itk.Image` format. + This method requires an optional dependency to be installed + separately from highdicom, specifically ``itk>=5.4.0``. + The Volume is converted to a 3D ``itk.image``. If its array's current datatype is not supported by ITK, it is safely cast to a compatible type where possible. If impossible to cast safely, @@ -3859,6 +3868,9 @@ def from_itk( ) -> Self: """Construct a Volume from an `itk.Image`. + This method requires an optional dependency to be installed + separately from highdicom, specifically ``itk>=5.4.0``. + The ``itk.Image`` is converted to a 3D Volume. Spatial metadata (spacing, direction, origin) is preserved with both highdicom and ITK using "LPS" convention. As ITK uses