diff --git a/src/highdicom/__init__.py b/src/highdicom/__init__.py index d3918f62..a482d02d 100644 --- a/src/highdicom/__init__.py +++ b/src/highdicom/__init__.py @@ -42,11 +42,13 @@ CoordinateSystemNames, ContentQualificationValues, DimensionOrganizationTypeValues, + InterpolationMethods, LateralityValues, PadModes, PatientSexValues, PhotometricInterpretationValues, PixelIndexDirections, + PixelDataKeywords, PixelRepresentationValues, PlanarConfigurationValues, PatientOrientationValuesBiped, @@ -59,6 +61,7 @@ ) from highdicom import frame from highdicom.image import ( + DimensionIndexSequence, Image, imread, get_volume_from_series, @@ -79,7 +82,6 @@ __all__ = [ - 'RGB_COLOR_CHANNEL_DESCRIPTOR', 'AlgorithmIdentificationSequence', 'AnatomicalOrientationTypeValues', 'AxisHandedness', @@ -88,8 +90,10 @@ 'ContentQualificationValues', 'ContributingEquipment', 'CoordinateSystemNames', + 'DimensionIndexSequence', 'DimensionOrganizationTypeValues', 'Image', + 'InterpolationMethods', 'IssuerOfIdentifier', 'LUT', 'LateralityValues', @@ -102,8 +106,9 @@ 'PatientOrientationValuesQuadruped', 'PatientSexValues', 'PhotometricInterpretationValues', - 'PixelMeasuresSequence', + 'PixelDataKeywords', 'PixelIndexDirections', + 'PixelMeasuresSequence', 'PixelRepresentationValues', 'PlanarConfigurationValues', 'PlaneOrientationSequence', @@ -111,9 +116,10 @@ 'PresentationLUT', 'PresentationLUTShapeValues', 'PresentationLUTTransformation', + 'RGBColorChannels', + 'RGB_COLOR_CHANNEL_DESCRIPTOR', 'ReferencedImageSequence', 'RescaleTypeValues', - 'RGBColorChannels', 'SOPClass', 'SegmentedPaletteColorLUT', 'SpecimenCollection', @@ -134,6 +140,7 @@ 'ann', 'color', 'frame', + 'get_volume_from_series', 'imread', 'io', 'ko', @@ -146,5 +153,4 @@ 'spatial', 'sr', 'utils', - 'get_volume_from_series', ] diff --git a/src/highdicom/_module_utils.py b/src/highdicom/_module_utils.py index 887ed5a8..50785c3d 100644 --- a/src/highdicom/_module_utils.py +++ b/src/highdicom/_module_utils.py @@ -1,5 +1,5 @@ from enum import Enum -from typing import Any, Dict, List, Optional, Union +from typing import Any from collections.abc import Sequence from pydicom import Dataset diff --git a/src/highdicom/_pyramid.py b/src/highdicom/_pyramid.py new file mode 100644 index 00000000..8ce434ea --- /dev/null +++ b/src/highdicom/_pyramid.py @@ -0,0 +1,409 @@ +"""Tools for constructing multi-resolution pyramids.""" +from typing import Any +from collections.abc import Generator, Sequence + +import numpy as np +from PIL import Image +from pydicom import Dataset +from pydicom.uid import VLWholeSlideMicroscopyImageStorage + +from highdicom.content import PixelMeasuresSequence +from highdicom.enum import InterpolationMethods +from highdicom.uid import UID + + +def iter_derived_pyramid_levels( + source_images: Sequence[Dataset], + pixel_arrays: Sequence[np.ndarray], + *, + interpolator: InterpolationMethods | str = InterpolationMethods.NEAREST, + downsample_factors: Sequence[float] | None = None, + sop_instance_uids: list[str] | None = None, +) -> Generator[ + tuple[ + Dataset, + np.ndarray, + PixelMeasuresSequence, + int, + str, + ], + None, + None, + ]: + """Create levels of derived multi-resolution pyramid with downsampling. + + A multi-resolution pyramid represents the same derived image array at + multiple resolutions. This function is a general-purpose function for + iterating through downsampled versions of an image for inclusion in a + derived pyramid series. + + This function handles multiple related scenarios: + + * Constructing a derived pyramid of a source image pyramid given a + derived pixel array of the highest resolution source image. + Highdicom performs the downsampling automatically to match the + resolution of the other source images. For this case, pass multiple + ``source_images`` and a single item in ``pixel_arrays``. + * Constructing a derived pyramid of a source image pyramid given + user-provided derived pixel arrays for each level in the source pyramid. + For this case, pass multiple ``source_images`` and a matching number of + ``pixel_arrays``. + * Constructing a derived pyramid of a single source image given multiple + user-provided downsampled derived pixel arrays. For this case, pass a + single item in ``source_images``, and multiple items in + ``pixel_arrays``). + * Constructing a derived pyramid of a single source image and a single + derived pixel array by downsampling by a given list of + ``downsample_factors``. For this case, pass a single item in + ``source_images``, a single item in ``pixel_arrays``, and a list of one + or more desired ``downsample_factors``. + + In all cases, the items in both ``source_images`` and ``pixel_arrays`` + should be sorted in pyramid order from highest resolution (smallest + spacing) to lowest resolution (largest spacing), and the pixel array + in ``pixel_arrays[0]`` must be derived from the source image in + ``source_images[0]`` with spatial locations preserved (a one-to-one + correspondence between pixels in the source image's total pixel matrix and + the provided derived pixel array). + + In all cases, the provided pixel arrays should be total pixel matrices. + Tiling is performed automatically. + + Parameters + ---------- + source_images: Sequence[pydicom.Dataset] + List of source images. If there are multiple source images, they should + be from the same series and pyramid. + pixel_arrays: Sequence[numpy.ndarray] + List of derived pixel arrays. Each should be a total pixel matrix, i.e. + have shape (rows, columns), (1, rows, columns), or (1, rows, columns, + segments/channels). + interpolator: highdicom.InterpolationMethods | str, optional + Interpolation method used for the downsampling operations. + downsample_factors: Optional[Sequence[float]], optional + Factors by which to downsample the pixel array to create each of the + output levels. This should be provided if and only if a single source + image and single pixel array are provided. Note that the original array + is always used to create the first output, so the number of created + levels is one greater than the number of items in this list. Items must + be numbers greater than 1 and sorted in ascending order. A downsampling + factor of *n* implies that the output array is *1/n* time the size of + input pixel array. For example a list ``[2, 4, 8]`` would be produce 4 + output levels. The first is the same size as the original pixel array, + the next is half the size, the next is a quarter of the size of the + original, and the last is one eighth the size of the original. Output + sizes are rounded to the nearest integer. + sop_instance_uids: Optional[List[str]], optional + SOP instance UIDS of the output instances. If not specified, UIDs are + generated automatically using highdicom's prefix. + + Yields + ------ + source_image: + Dataset to use as the source image at this level. This will be one of + the inputs provided to ``source_images``. + pixel_array: numpy.ndarray + Pixel array at this level. This will either be one of the inputs + provided to ``pixel_arrays`` or an array created by downsampling the + single provided pixel array. + pixel_measures: highdicom.PixelMeasuresSequence + Pixel measures object describing the pixel spacing of the pixel array + at this level. + instance_number: int + Instance number to include in the instance at this level. + sop_instance_uid: str + SOP Instance UID to include in the instance at this level. + + """ + n_sources = len(source_images) + n_pix_arrays = len(pixel_arrays) + + if n_sources == 0: + raise ValueError( + 'Argument "source_images" must not be empty.' + ) + if n_pix_arrays == 0: + raise ValueError( + 'Argument "pixel_arrays" must not be empty.' + ) + + if n_sources == 1 and n_pix_arrays == 1: + if downsample_factors is None: + raise TypeError( + 'Argument "downsample_factors" must be provided when providing ' + 'only a single source image and pixel array.' + ) + if len(downsample_factors) < 1: + raise ValueError('Argument "downsample_factors" may not be empty.') + if any(f <= 1.0 for f in downsample_factors): + raise ValueError( + 'All items in "downsample_factors" must be greater than 1.' + ) + if len(downsample_factors) > 1: + if any( + z1 > z2 for z1, z2 in zip( + downsample_factors[:-1], + downsample_factors[1:] + ) + ): + raise ValueError( + 'Items in argument "downsample_factors" must be sorted in ' + 'ascending order.' + ) + n_outputs = len(downsample_factors) + 1 # including original + else: + if downsample_factors is not None: + raise TypeError( + 'Argument "downsample_factors" must not be provided when ' + 'multiple source images or pixel arrays are provided.' + ) + if n_sources > 1 and n_pix_arrays > 1: + if n_sources != n_pix_arrays: + raise ValueError( + 'If providing multiple source images and multiple pixel ' + 'arrays, the number of items in the two lists must match.' + ) + n_outputs = n_sources + else: + # Either n_sources > 1 or n_pix_arrays > 1 but not both + n_outputs = max(n_sources, n_pix_arrays) + + if sop_instance_uids is not None: + if len(sop_instance_uids) != n_outputs: + raise ValueError( + 'Number of specified SOP Instance UIDs does not match number ' + 'of output images.' + ) + + # Check that source images are WSI + for im in source_images: + if im.SOPClassUID != VLWholeSlideMicroscopyImageStorage: + raise ValueError( + 'Source images must have IOD VLWholeSlideMicroscopyImage' + ) + + # Check the source images are appropriately ordered + for index in range(1, len(source_images)): + r0 = source_images[index - 1].TotalPixelMatrixRows + c0 = source_images[index - 1].TotalPixelMatrixColumns + r1 = source_images[index].TotalPixelMatrixRows + c1 = source_images[index].TotalPixelMatrixColumns + + if r0 <= r1 or c0 <= c1: + raise ValueError( + 'Items in argument "source_images" must be strictly ordered in ' + 'decreasing resolution.' + ) + + # Check that the source images are from the same series and pyramid + if len(source_images) > 1: + series_uid = source_images[0].SeriesInstanceUID + if not all( + dcm.SeriesInstanceUID == series_uid + for dcm in source_images[1:] + ): + raise ValueError( + 'All source images should belong to the same series.' + ) + if not all(hasattr(dcm, 'PyramidUID') for dcm in source_images): + raise ValueError( + 'All source images should belong to the same pyramid ' + '(share a Pyramid UID).' + ) + pyramid_uid = source_images[0].PyramidUID + if not all( + dcm.PyramidUID == pyramid_uid + for dcm in source_images[1:] + ): + raise ValueError( + 'All source images should belong to the same pyramid ' + '(share a Pyramid UID).' + ) + + # Check that pixel arrays have an appropriate shape + if len({p.ndim for p in pixel_arrays}) != 1: + raise ValueError( + 'Each item of argument "pixel_arrays" must have the same number of ' + 'dimensions.' + ) + if pixel_arrays[0].ndim == 4: + n_channels = pixel_arrays[0].shape[3] + else: + n_channels = None + + # Map the highdicom interpolation methods enum to value used by Pillow + interpolator = InterpolationMethods(interpolator) + resampler = { + InterpolationMethods.NEAREST: Image.Resampling.NEAREST, + InterpolationMethods.LINEAR: Image.Resampling.BILINEAR, + }[interpolator] + + # Checks on consistency of the pixel arrays + dtype = pixel_arrays[0].dtype + for pixel_array in pixel_arrays: + if pixel_array.ndim not in (2, 3, 4): + raise ValueError( + 'Each item of argument "pixel_arrays" must be a NumPy array ' + 'with 2, 3, or 4 dimensions.' + ) + if pixel_array.ndim > 2 and pixel_array.shape[0] != 1: + raise ValueError( + 'Each item of argument "pixel_arrays" must contain a single ' + 'frame, with a size of 1 along dimension 0.' + ) + if pixel_array.dtype != dtype: + raise TypeError( + 'Each item of argument "pixel_arrays" must have ' + 'the same datatype.' + ) + if pixel_array.ndim == 4: + if pixel_array.shape[3] != n_channels: + raise ValueError( + 'Each item of argument "pixel_arrays" must have ' + 'the same shape down axis 3.' + ) + + # Check the pixel arrays are appropriately ordered + for index in range(1, len(pixel_arrays)): + arr0 = pixel_arrays[index - 1] + arr1 = pixel_arrays[index] + + if arr0.ndim == 2: + r0 = arr0.shape[:2] + c0 = arr0.shape[:2] + else: + r0 = arr0.shape[1:3] + c0 = arr0.shape[1:3] + + if arr1.ndim == 2: + r1 = arr1.shape[:2] + c1 = arr1.shape[:2] + else: + r1 = arr1.shape[1:3] + c1 = arr1.shape[1:3] + + if r0 <= r1 or c0 <= c1: + raise ValueError( + 'Items in argument "pixel_arrays" must be strictly ordered in ' + 'decreasing size.' + ) + + # Check that input dimensions match + for index, (source_image, pixel_array) in enumerate( + zip(source_images, pixel_arrays) + ): + src_shape = ( + source_image.TotalPixelMatrixRows, + source_image.TotalPixelMatrixColumns + ) + pix_shape = ( + pixel_array.shape[1:3] if pixel_array.ndim > 2 + else pixel_array.shape + ) + if pix_shape != src_shape: + raise ValueError( + "The shape of each provided pixel array must match the shape " + "of the total pixel matrix of the corresponding source image. " + f"Got pixel array of shape {pix_shape} for a source image of " + f"shape {src_shape} at index {index}." + ) + + if n_pix_arrays == 1: + # Create a pillow image for use later with resizing + if pixel_arrays[0].ndim == 2: + pil_images = [Image.fromarray(pixel_arrays[0])] + elif pixel_arrays[0].ndim == 3: + # Remove frame dimension before casting + pil_images = [Image.fromarray(pixel_arrays[0][0])] + else: # ndim = 4 + # One "Image" for each channel + pil_images = [ + Image.fromarray(pixel_arrays[0][0, :, :, i]) + for i in range(pixel_arrays[0].shape[3]) + ] + + # Work "up" pyramid from high to low resolution + for output_level in range(n_outputs): + if n_sources > 1: + source_image = source_images[output_level] + else: + source_image = source_images[0] + + if n_pix_arrays > 1: + pixel_array = pixel_arrays[output_level] + else: + if output_level == 0: + pixel_array = pixel_arrays[0] + else: + if n_sources > 1: + output_size = ( + source_image.TotalPixelMatrixColumns, + source_image.TotalPixelMatrixRows + ) + else: + f = downsample_factors[output_level - 1] + output_size = ( + int(source_images[0].TotalPixelMatrixColumns / f), + int(source_images[0].TotalPixelMatrixRows / f) + ) + + # Resize each channel individually + resized_images = [ + np.array(im.resize(output_size, resampler)) + for im in pil_images + ] + if len(resized_images) > 1: + pixel_array = np.stack(resized_images, axis=-1)[None] + else: + pixel_array = resized_images[0] + + # Standardize shape of pixel array to include singleton frame dimension + if pixel_array.ndim == 2: + pixel_array = pixel_array[None] + + if n_sources == 1: + source_pixel_measures = ( + source_image + .SharedFunctionalGroupsSequence[0] + .PixelMeasuresSequence[0] + ) + src_pixel_spacing = source_pixel_measures.PixelSpacing + src_slice_thickness = source_pixel_measures.SliceThickness + + if pixel_arrays[0].ndim == 2: + # No frame dimension + orig_n_rows = pixel_arrays[0].shape[0] + orig_n_cols = pixel_arrays[0].shape[1] + else: + # Ignore 0th frame dimension + orig_n_rows = pixel_arrays[0].shape[1] + orig_n_cols = pixel_arrays[0].shape[2] + + row_spacing = ( + src_pixel_spacing[0] * + (orig_n_rows / pixel_array.shape[1]) + ) + column_spacing = ( + src_pixel_spacing[1] * + (orig_n_cols / pixel_array.shape[2]) + ) + pixel_measures = PixelMeasuresSequence( + pixel_spacing=(row_spacing, column_spacing), + slice_thickness=src_slice_thickness + ) + else: + # This will be copied from the source image + pixel_measures = None + + if sop_instance_uids is None: + sop_instance_uid = UID() + else: + sop_instance_uid = sop_instance_uids[output_level] + + yield ( + source_image, + pixel_array, + pixel_measures, + output_level + 1, # instance number + sop_instance_uid, + ) diff --git a/src/highdicom/ann/sop.py b/src/highdicom/ann/sop.py index 1607a83b..41982e21 100644 --- a/src/highdicom/ann/sop.py +++ b/src/highdicom/ann/sop.py @@ -28,11 +28,14 @@ PixelOriginInterpretationValues, ) from highdicom.ann.content import AnnotationGroup +from highdicom.content import ( + _add_content_information, + ContentCreatorIdentificationCodeSequence, +) from highdicom.base import SOPClass, _check_little_endian from highdicom.base_content import ContributingEquipment from highdicom.io import _wrapped_dcmread from highdicom.sr.coding import CodedConcept -from highdicom.valuerep import check_person_name, _check_code_string class MicroscopyBulkSimpleAnnotations(SOPClass): @@ -67,6 +70,9 @@ def __init__( contributing_equipment: Sequence[ ContributingEquipment ] | None = None, + content_creator_identification: None | ( + ContentCreatorIdentificationCodeSequence + ) = None, **kwargs: Any ) -> None: """ @@ -116,6 +122,9 @@ def __init__( contributing_equipment: Sequence[highdicom.ContributingEquipment] | None, optional Additional equipment that has contributed to the acquisition, creation or modification of this instance. + content_creator_identification: Union[highdicom.ContentCreatorIdentificationCodeSequence, None], optional + Identifying information for the person who created the content of + this microscopy bulk simple annotations instance. **kwargs: Any, optional Additional keyword arguments that will be passed to the constructor of `highdicom.base.SOPClass` @@ -178,15 +187,16 @@ def __init__( self._add_contributing_equipment(contributing_equipment, src_img) # Microscopy Bulk Simple Annotations - if content_label is not None: - _check_code_string(content_label) - self.ContentLabel = content_label - else: - self.ContentLabel = f'{src_img.Modality}_ANN' - self.ContentDescription = content_description - if content_creator_name is not None: - check_person_name(content_creator_name) - self.ContentCreatorName = content_creator_name + _add_content_information( + dataset=self, + content_label=( + content_label if content_label is not None + else f'{src_img.Modality}_ANN' + ), + content_description=content_description, + content_creator_name=content_creator_name, + content_creator_identification=content_creator_identification, + ) self.AnnotationCoordinateType = coordinate_type.value if coordinate_type == AnnotationCoordinateTypeValues.SCOORD: diff --git a/src/highdicom/content.py b/src/highdicom/content.py index 75e40f70..e28e3a9b 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -2,18 +2,24 @@ from collections import Counter import datetime from copy import deepcopy +from io import BytesIO from typing import cast from collections.abc import Sequence from typing_extensions import Self import numpy as np from PIL import ImageColor +from PIL.ImageCms import ImageCmsProfile from pydicom.dataset import Dataset from pydicom import DataElement from pydicom.sequence import Sequence as DataElementSequence from pydicom.sr.coding import Code from pydicom.sr.codedict import codes -from pydicom.valuerep import DS, format_number_as_ds +from pydicom.valuerep import ( + DS, + format_number_as_ds, + PersonName, +) from pydicom.uid import SegmentationStorage from highdicom.enum import ( @@ -43,9 +49,11 @@ ) from highdicom.uid import UID from highdicom.valuerep import ( + _check_code_string, _check_long_string, _check_long_text, _check_short_text, + check_person_name, ) from highdicom._module_utils import ( check_required_attributes, @@ -2320,6 +2328,11 @@ def __init__( raise TypeError( 'Argument "window_width" must not be an empty sequence.' ) + for x in window_width: + if x <= 0: + raise ValueError( + 'Window width must be greater than zero.' + ) self.WindowWidth = [ format_number_as_ds(float(x)) for x in window_width ] @@ -2329,6 +2342,8 @@ def __init__( 'Length of "window_width" must match length of ' '"window_center".' ) + if window_width <= 0: + raise ValueError('Window width must be greater than zero.') self.WindowWidth = format_number_as_ds(float(window_width)) if window_explanation is not None: if window_center is None: @@ -3663,3 +3678,103 @@ def apply( first_mapped_value=self.first_mapped_value, dtype=dtype, ) + + +def _add_content_information( + dataset, + content_label: str, + content_description: str | None, + content_creator_name: str | PersonName | None = None, + content_creator_identification: None | ( + ContentCreatorIdentificationCodeSequence + ) = None, +): + _check_code_string(content_label) + dataset.ContentLabel = content_label + if content_description is not None: + _check_long_string(content_description) + dataset.ContentDescription = content_description + if content_creator_name is not None: + check_person_name(content_creator_name) + dataset.ContentCreatorName = content_creator_name + if content_creator_identification is not None: + if not isinstance( + content_creator_identification, + ContentCreatorIdentificationCodeSequence + ): + raise TypeError( + 'Argument "content_creator_identification" must be of type ' + 'ContentCreatorIdentificationCodeSequence.' + ) + dataset.ContentCreatorIdentificationCodeSequence = \ + content_creator_identification + + +def _add_icc_profile_attributes( + dataset: Dataset, + icc_profile: bytes +) -> None: + """Add attributes of module ICC Profile. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset to which attributes should be added + icc_profile: bytes + ICC color profile to include in the presentation state. + The profile must follow the constraints listed in :dcm:`C.11.15 + `. + + """ + if icc_profile is None: + raise TypeError('Argument "icc_profile" is required.') + + cms_profile = ImageCmsProfile(BytesIO(icc_profile)) + device_class = cms_profile.profile.device_class.strip() + if device_class not in ('scnr', 'spac'): + raise ValueError( + 'The device class of the ICC Profile must be "scnr" or "spac", ' + f'got "{device_class}".' + ) + color_space = cms_profile.profile.xcolor_space.strip() + if color_space != 'RGB': + raise ValueError( + 'The color space of the ICC Profile must be "RGB", ' + f'got "{color_space}".' + ) + pcs = cms_profile.profile.connection_space.strip() + if pcs not in ('Lab', 'XYZ'): + raise ValueError( + 'The profile connection space of the ICC Profile must ' + f'be "Lab" or "XYZ", got "{pcs}".' + ) + + dataset.ICCProfile = icc_profile + + +def _add_palette_color_lookup_table_attributes( + dataset: Dataset, + palette_color_lut_transformation: PaletteColorLUTTransformation +) -> None: + """Add attributes from the Palette Color Lookup Table module. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset to which attributes should be added + palette_color_lut_transformation: highdicom.PaletteColorLUTTransformation + Description of the Palette Color LUT Transformation for transforming + grayscale into RGB color pixel values + + """ # noqa: E501 + if not isinstance( + palette_color_lut_transformation, + PaletteColorLUTTransformation + ): + raise TypeError( + 'Argument "palette_color_lut_transformation" must be of type ' + 'PaletteColorLUTTransformation.' + ) + + for element in palette_color_lut_transformation: + dataset[element.tag] = element diff --git a/src/highdicom/enum.py b/src/highdicom/enum.py index 0ca66faa..b1f25a51 100644 --- a/src/highdicom/enum.py +++ b/src/highdicom/enum.py @@ -350,7 +350,7 @@ class AxisHandedness(Enum): """ - LEFT_HANDED = "LEFT_HANDED" + LEFT_HANDED = 'LEFT_HANDED' """ The unit vectors of the first, second and third axes form a left hand when @@ -360,7 +360,7 @@ class AxisHandedness(Enum): """ - RIGHT_HANDED = "RIGHT_HANDED" + RIGHT_HANDED = 'RIGHT_HANDED' """ The unit vectors of the first, second and third axes form a right hand when @@ -369,3 +369,28 @@ class AxisHandedness(Enum): vector, and the middle finger representing the third vector. """ + + +class InterpolationMethods(Enum): + + """Interpolation methods that may be used for resampling arrays.""" + + NEAREST = 'NEAREST' + """Nearest-neighbor interpolator.""" + + LINEAR = 'LINEAR' + """Linear (or bi-linear or tri-linear) interpolator.""" + + +class PixelDataKeywords(Enum): + + """Keywords used to store pixel data.""" + + PIXEL_DATA = 'PixelData' + """Integer-valued pixel data of any size (Image Pixel Module).""" + + FLOAT_PIXEL_DATA = 'FloatPixelData' + """32-bit-float-valued pixel data (Floating Point Image Pixel Module).""" + + DOUBLE_FLOAT_PIXEL_DATA = 'DoubleFloatPixelData' + """64-bit-float-valued pixel data (Double Floating Point Image Pixel Module).""" # noqa: E501 diff --git a/src/highdicom/frame.py b/src/highdicom/frame.py index c030aab4..38eb24cf 100644 --- a/src/highdicom/frame.py +++ b/src/highdicom/frame.py @@ -24,6 +24,7 @@ PhotometricInterpretationValues, PixelRepresentationValues, PlanarConfigurationValues, + PixelDataKeywords, ) logger = logging.getLogger(__name__) @@ -345,6 +346,7 @@ def encode_frame( planar_configuration=planar_configuration, **kwargs, ) + return data @@ -357,10 +359,10 @@ def decode_frame( bits_allocated: int, bits_stored: int, photometric_interpretation: PhotometricInterpretationValues | str, - pixel_representation: PixelRepresentationValues | int = 0, + pixel_representation: PixelRepresentationValues | int | None = 0, planar_configuration: PlanarConfigurationValues | int | None = None, index: int = 0, - pixel_keyword: str = 'PixelData', + pixel_keyword: PixelDataKeywords | str = PixelDataKeywords.PIXEL_DATA, ) -> np.ndarray: """Decode pixel data of an individual frame. @@ -385,7 +387,7 @@ def decode_frame( Photometric interpretation pixel_representation: Union[highdicom.PixelRepresentationValues, int, None], optional Whether pixel samples are represented as unsigned integers or - 2's complements + 2's complements. May only be None for floating point pixel data. planar_configuration: Union[highdicom.PlanarConfigurationValues, int, None], optional Whether color samples are encoded by pixel (``R1G1B1R2G2B2...``) or by plane (``R1R2...G1G2...B1B2...``). @@ -398,7 +400,7 @@ def decode_frame( the byte array. In all other situations, this parameter is not required and will have no effect (since decoding a frame does not depend on the index of the frame). - pixel_keyword: str, optional + pixel_keyword: highdicom.enum.PixelDataKeywords | str, optional Keyword of the attribute where the pixel data value was stored. Should be ``'PixelData'``, ``'FloatPixelData'``, or ``'DoubleFloatPixelData'``. @@ -445,9 +447,11 @@ def decode_frame( pixel_array = unpacked_frame[pixel_offset:pixel_offset + n_pixels] return pixel_array.reshape(rows, columns) - pixel_representation = PixelRepresentationValues( - pixel_representation - ).value + pixel_keyword = PixelDataKeywords(pixel_keyword) + if pixel_keyword == PixelDataKeywords.PIXEL_DATA: + pixel_representation = PixelRepresentationValues( + pixel_representation + ).value photometric_interpretation = PhotometricInterpretationValues( photometric_interpretation ).value @@ -480,7 +484,7 @@ def decode_frame( bits_stored=bits_stored, photometric_interpretation=photometric_interpretation, pixel_representation=pixel_representation, - pixel_keyword=pixel_keyword, + pixel_keyword=pixel_keyword.value, **kwargs, ) return array diff --git a/src/highdicom/image.py b/src/highdicom/image.py index 07894048..80a54f2b 100644 --- a/src/highdicom/image.py +++ b/src/highdicom/image.py @@ -1,5 +1,6 @@ """Tools for working with general DICOM images.""" from collections import Counter, defaultdict +from concurrent.futures import Executor, Future, ProcessPoolExecutor from contextlib import contextmanager, nullcontext from copy import deepcopy from enum import Enum @@ -12,13 +13,25 @@ BinaryIO, cast, ) -from collections.abc import Iterable, Iterator, Generator, Sequence +from collections.abc import ( + Callable, + Generator, + Iterable, + Iterator, + Sequence, +) from typing_extensions import Self +import warnings import numpy as np from pydicom import Dataset from pydicom.dataelem import DataElement -from pydicom.encaps import get_frame +from pydicom.encaps import ( + get_frame, + encapsulate, + encapsulate_extended, +) +from pydicom.sequence import Sequence as DataElementSequence from pydicom.tag import BaseTag from pydicom.datadict import ( get_entry, @@ -27,29 +40,43 @@ ) from pydicom.filebase import DicomIO, DicomBytesIO from pydicom.multival import MultiValue +from pydicom.pixels.utils import pack_bits from pydicom.sr.coding import Code from pydicom.sr.codedict import codes -from pydicom.uid import ParametricMapStorage +from pydicom.uid import ParametricMapStorage, SegmentationStorage from pydicom.valuerep import format_number_as_ds from highdicom._module_utils import ( + ModuleUsageValues, + construct_module_tree, does_iod_have_pixel_data, + get_module_usage, is_multiframe_image, ) +from highdicom._modules import MODULE_ATTRIBUTE_MAP from highdicom.base import SOPClass, _check_little_endian from highdicom.color import ColorManager from highdicom.content import ( + _add_icc_profile_attributes, + _add_palette_color_lookup_table_attributes, LUT, + PaletteColorLUTTransformation, PixelMeasuresSequence, PlaneOrientationSequence, PlanePositionSequence, VOILUTTransformation ) from highdicom.enum import ( + AxisHandedness, CoordinateSystemNames, DimensionOrganizationTypeValues, + PhotometricInterpretationValues, + PixelDataKeywords, + PixelIndexDirections, + PixelRepresentationValues, + PlanarConfigurationValues, ) -from highdicom.frame import decode_frame +from highdicom.frame import decode_frame, encode_frame from highdicom.io import ImageFileReader, _wrapped_dcmread from highdicom.pixels import ( _check_rescale_dtype, @@ -63,23 +90,31 @@ from highdicom.spatial import ( _are_orientations_coplanar, _are_orientations_equal, + _get_slice_distances, _get_spatial_information, ImageToReferenceTransformer, compute_tile_positions_per_frame, get_image_coordinate_system, get_normal_vector, get_series_volume_positions, + get_tile_array, get_volume_positions, is_tiled_image, + map_pixel_into_coordinate_system, ) from highdicom.sr.coding import CodedConcept from highdicom.uid import UID as UID from highdicom.utils import ( iter_tiled_full_frame_data, are_plane_positions_tiled_full, + compute_plane_position_slide_per_frame, +) +from highdicom.valuerep import ( + _check_long_string, ) from highdicom.volume import ( _DCM_PYTHON_TYPE_MAP, + ChannelDescriptor, VolumeGeometry, Volume, RGB_COLOR_CHANNEL_DESCRIPTOR, @@ -104,6 +139,7 @@ 'SL': 'INTEGER', 'SS': 'INTEGER', 'ST': 'TEXT', + 'SQ': 'TEXT', # only for codes 'UI': 'TEXT', 'UL': 'INTEGER', 'UR': 'TEXT', @@ -112,10 +148,24 @@ 'UT': 'TEXT', } -# This code is needed many times in loops so we precompute it +# These codes are needed many times in loops so we precompute them _PURPOSE_CODE = CodedConcept.from_code( codes.cid7202.SourceImageForImageProcessingOperation ) +_DERIVATION_CODE = CodedConcept.from_code( + codes.cid7203.SegmentationImageDerivation +) + + +# Tags used for spatial dimensions +_SPATIAL_TAGS = { + 0x0020_0032, # ImagePositionPatient + 0x0040_072a, # XOffsetInSlideCoordinateSystem + 0x0040_073a, # YOffsetInSlideCoordinateSystem + 0x0040_074a, # ZOffsetInSlideCoordinateSystem + 0x0048_021e, # ColumnPositionInTotalImagePixelMatrix + 0x0048_021f, # RowPositionInTotalImagePixelMatrix +} class _ImageColorType(Enum): @@ -458,17 +508,38 @@ def __init__( # Determine input type and range of values input_range = None - if ( - image.SOPClassUID == ParametricMapStorage and - image.BitsAllocated > 16 - ): - # Parametric Maps are the only SOP Class (currently) that allows - # floating point pixels - if image.BitsAllocated == 32: - self.input_dtype = np.dtype(np.float32) - elif image.BitsAllocated == 64: - self.input_dtype = np.dtype(np.float64) + + # Determine what pixel data keyword is present in the image + for kw in [ + PixelDataKeywords.PIXEL_DATA, + PixelDataKeywords.FLOAT_PIXEL_DATA, + PixelDataKeywords.DOUBLE_FLOAT_PIXEL_DATA, + ]: + if kw.value in image: + self.pixel_keyword = kw + break else: + # No pixel data keyword, which may be due to just the header being + # loaded. Parametric Maps are the only SOP Class (currently) that + # allows floating point pixels. BitsStored is only present for + # PixelData (i.e. when the ImagePixel module is used) and as such + # can be used to decude the pixel type + if ( + image.SOPClassUID == ParametricMapStorage and + 'BitsStored' not in image + ): + if image.BitsAllocated == 32: + self.input_dtype = np.dtype(np.float32) + self.pixel_keyword = PixelDataKeywords.FLOAT_PIXEL_DATA + elif image.BitsAllocated == 64: + self.input_dtype = np.dtype(np.float64) + self.pixel_keyword = ( + PixelDataKeywords.DOUBLE_FLOAT_PIXEL_DATA + ) + else: + self.pixel_keyword = PixelDataKeywords.PIXEL_DATA + + if self.pixel_keyword == PixelDataKeywords.PIXEL_DATA: if image.PixelRepresentation == 1: if image.BitsAllocated == 8: self.input_dtype = np.dtype(np.int8) @@ -969,7 +1040,7 @@ def __init__( self.bits_allocated = image.BitsAllocated self.bits_stored = image.get('BitsAllocated', image.BitsAllocated) self.photometric_interpretation = image.PhotometricInterpretation - self.pixel_representation = image.PixelRepresentation + self.pixel_representation = image.get('PixelRepresentation') self.planar_configuration = image.get('PlanarConfiguration') def __call__( @@ -1009,6 +1080,7 @@ def __call__( pixel_representation=self.pixel_representation, planar_configuration=self.planar_configuration, index=frame_index, + pixel_keyword=self.pixel_keyword, ) elif isinstance(frame, np.ndarray): frame_out = frame @@ -1085,118 +1157,1504 @@ def __init__( ): """ - Parameters - ---------- - table_name: str, - Name of the temporary table. - column_data: Iterable[Sequence[Any]] - Column data to place into the table. - column_defs: Sequence[str] - SQL syntax strings defining each column in the temporary table, one - string per column. + Parameters + ---------- + table_name: str, + Name of the temporary table. + column_data: Iterable[Sequence[Any]] + Column data to place into the table. + column_defs: Sequence[str] + SQL syntax strings defining each column in the temporary table, one + string per column. + + """ + self.table_name = table_name + self.column_defs = list(column_defs) + + # It is important to convert numpy arrays to lists, otherwise the + # values get inserted into the table as binary values and this leads to + # very difficult to detect failures + if isinstance(column_data, np.ndarray): + column_data = column_data.tolist() + + sanitized_column_data = [] + + for col in column_data: + if isinstance(col, np.ndarray): + col = col.tolist() + + sanitized_column_data.append(col) + + self.column_data = sanitized_column_data + + +class DimensionIndexSequence(DataElementSequence): + + """Sequence of data elements describing dimension indices for the patient + or slide coordinate system based on the Dimension Index functional + group macro. + + Note + ---- + The order of indices is fixed. + + """ + + def __init__( + self, + coordinate_system: str | CoordinateSystemNames | None, + functional_groups_module: str, + channel_dimensions: Sequence[ChannelDescriptor] | None = None, + ) -> None: + """ + Parameters + ---------- + coordinate_system: Union[str, highdicom.CoordinateSystemNames, None] + Subject (``"PATIENT"`` or ``"SLIDE"``) that was the target of + imaging. If None, the imaging does not belong within a frame of + reference. + + """ + super().__init__() + if coordinate_system is None: + self._coordinate_system = None + else: + self._coordinate_system = CoordinateSystemNames(coordinate_system) + + if functional_groups_module not in MODULE_ATTRIBUTE_MAP: + raise KeyError( + f"Value of '{functional_groups_module}' for parameter " + "'functional_groups_module' is not an existing module name." + ) + if ( + ( + functional_groups_module == + 'sparse-multi-frame-functional-groups' + ) or ( + not functional_groups_module.endswith( + '-multi-frame-functional-groups' + ) + ) + ): + raise ValueError( + f"Value of '{functional_groups_module}' for parameter " + "'functional_groups_module' is not a multi-frame " + "functional groups module." + ) + + module_tree = construct_module_tree(functional_groups_module) + per_frame_tree = ( + module_tree + ['attributes'] + ['PerFrameFunctionalGroupsSequence'] + ['attributes'] + ) + + dim_uid = UID() + + if channel_dimensions is not None: + + for desc in channel_dimensions: + + if not isinstance(desc, ChannelDescriptor): + raise TypeError( + "Items of 'channel_dimensions' must have type " + 'highdicom.ChannelDescriptor.' + ) + + if desc.is_custom: + raise ValueError( + "Channels descriptors used in 'channel_dimensions' may " + 'not be custom.' + ) + + channel_index = Dataset() + channel_index.DimensionIndexPointer = desc.tag + for seq_name, seq_info in per_frame_tree.items(): + if desc.keyword in seq_info['attributes']: + break + else: + raise ValueError( + f"The specified channel descriptor ('{desc.keyword}') " + 'may not be used as an index within the specified ' + f"module ('{functional_groups_module}')." + ) + channel_index.FunctionalGroupPointer = tag_for_keyword(seq_name) + channel_index.DimensionOrganizationUID = dim_uid + channel_index.DimensionDescriptionLabel = desc.description + self.append(channel_index) + + if self._coordinate_system == CoordinateSystemNames.SLIDE: + + x_axis_index = Dataset() + x_axis_index.DimensionIndexPointer = tag_for_keyword( + 'XOffsetInSlideCoordinateSystem' + ) + x_axis_index.FunctionalGroupPointer = tag_for_keyword( + 'PlanePositionSlideSequence' + ) + x_axis_index.DimensionOrganizationUID = dim_uid + x_axis_index.DimensionDescriptionLabel = ( + 'X Offset in Slide Coordinate System' + ) + + y_axis_index = Dataset() + y_axis_index.DimensionIndexPointer = tag_for_keyword( + 'YOffsetInSlideCoordinateSystem' + ) + y_axis_index.FunctionalGroupPointer = tag_for_keyword( + 'PlanePositionSlideSequence' + ) + y_axis_index.DimensionOrganizationUID = dim_uid + y_axis_index.DimensionDescriptionLabel = ( + 'Y Offset in Slide Coordinate System' + ) + + z_axis_index = Dataset() + z_axis_index.DimensionIndexPointer = tag_for_keyword( + 'ZOffsetInSlideCoordinateSystem' + ) + z_axis_index.FunctionalGroupPointer = tag_for_keyword( + 'PlanePositionSlideSequence' + ) + z_axis_index.DimensionOrganizationUID = dim_uid + z_axis_index.DimensionDescriptionLabel = ( + 'Z Offset in Slide Coordinate System' + ) + + column_dimension_index = Dataset() + column_dimension_index.DimensionIndexPointer = tag_for_keyword( + 'ColumnPositionInTotalImagePixelMatrix' + ) + column_dimension_index.FunctionalGroupPointer = tag_for_keyword( + 'PlanePositionSlideSequence' + ) + column_dimension_index.DimensionOrganizationUID = dim_uid + column_dimension_index.DimensionDescriptionLabel = ( + 'Column Position In Total Image Pixel Matrix' + ) + + row_dimension_index = Dataset() + row_dimension_index.DimensionIndexPointer = tag_for_keyword( + 'RowPositionInTotalImagePixelMatrix' + ) + row_dimension_index.FunctionalGroupPointer = tag_for_keyword( + 'PlanePositionSlideSequence' + ) + row_dimension_index.DimensionOrganizationUID = dim_uid + row_dimension_index.DimensionDescriptionLabel = ( + 'Row Position In Total Image Pixel Matrix' + ) + + # Organize frames for each segment similar to TILED_FULL, with + # segment position changing least frequently, followed by position + # of the row (from top to bottom) and then position of the column + # (from left to right) changing most frequently + self.extend([ + row_dimension_index, + column_dimension_index, + x_axis_index, + y_axis_index, + z_axis_index, + ]) + + elif self._coordinate_system == CoordinateSystemNames.PATIENT: + + image_position_index = Dataset() + image_position_index.DimensionIndexPointer = tag_for_keyword( + 'ImagePositionPatient' + ) + image_position_index.FunctionalGroupPointer = tag_for_keyword( + 'PlanePositionSequence' + ) + image_position_index.DimensionOrganizationUID = dim_uid + image_position_index.DimensionDescriptionLabel = ( + 'Image Position Patient' + ) + + self.append(image_position_index) + + elif self._coordinate_system is None: + if channel_dimensions is None or len(channel_dimensions) == 0: + # Use frame label here just for the sake of using something + frame_label_index = Dataset() + frame_label_index.DimensionIndexPointer = tag_for_keyword( + 'FrameLabel' + ) + frame_label_index.FunctionalGroupPointer = tag_for_keyword( + 'FrameContentSequence' + ) + frame_label_index.DimensionOrganizationUID = dim_uid + frame_label_index.DimensionDescriptionLabel = 'Frame Label' + self.append(frame_label_index) + else: + raise ValueError( + f'Unknown coordinate system "{self._coordinate_system}"' + ) + + # Check that the resulting sequence does not contain duplicate + # attributes + ptrs = {item.DimensionIndexPointer for item in self} + if len(ptrs) < len(self): + raise ValueError( + 'Specified channel dimensions lead to duplicate dimension ' + 'indices.' + ) + + def get_plane_positions_of_image( + self, + image: Dataset + ) -> list[PlanePositionSequence]: + """Gets plane positions of frames in multi-frame image. + + Parameters + ---------- + image: Dataset + Multi-frame image + + Returns + ------- + List[highdicom.PlanePositionSequence] + Plane position of each frame in the image + + """ + is_multiframe = is_multiframe_image(image) + if not is_multiframe: + raise ValueError('Argument "image" must be a multi-frame image.') + + if self._coordinate_system is None: + raise ValueError( + 'Cannot calculate plane positions when images do not exist ' + 'within a frame of reference.' + ) + elif self._coordinate_system == CoordinateSystemNames.SLIDE: + if hasattr(image, 'PerFrameFunctionalGroupsSequence'): + plane_positions = [PlanePositionSequence.from_sequence( + item.PlanePositionSlideSequence + ) + for item in image.PerFrameFunctionalGroupsSequence + ] + else: + # If Dimension Organization Type is TILED_FULL, plane + # positions are implicit and need to be computed. + plane_positions = compute_plane_position_slide_per_frame(image) + else: + plane_positions = [ + PlanePositionSequence.from_sequence(item.PlanePositionSequence) + for item in image.PerFrameFunctionalGroupsSequence + ] + + return plane_positions + + def get_plane_positions_of_series( + self, + images: Sequence[Dataset] + ) -> list[PlanePositionSequence]: + """Gets plane positions for series of single-frame images. + + Parameters + ---------- + images: Sequence[Dataset] + Series of single-frame images + + Returns + ------- + List[highdicom.PlanePositionSequence] + Plane position of each frame in the image + + """ + is_multiframe = any([is_multiframe_image(img) for img in images]) + if is_multiframe: + raise ValueError( + 'Argument "images" must be a series of single-frame images.' + ) + + if self._coordinate_system is None: + raise ValueError( + 'Cannot calculate plane positions when images do not exist ' + 'within a frame of reference.' + ) + elif self._coordinate_system == CoordinateSystemNames.SLIDE: + plane_positions = [] + for img in images: + # Unfortunately, the image position is not specified relative to + # the top left corner but to the center of the image. + # Therefore, we need to compute the offset and subtract it. + center_item = img.ImageCenterPointCoordinatesSequence[0] + x_center = center_item.XOffsetInSlideCoordinateSystem + y_center = center_item.YOffsetInSlideCoordinateSystem + z_center = center_item.ZOffsetInSlideCoordinateSystem + offset_coordinate = map_pixel_into_coordinate_system( + index=((img.Columns / 2, img.Rows / 2)), + image_position=(x_center, y_center, z_center), + image_orientation=img.ImageOrientationSlide, + pixel_spacing=img.PixelSpacing + ) + center_coordinate = np.array((0., 0., 0.), dtype=float) + origin_coordinate = center_coordinate - offset_coordinate + plane_positions.append( + PlanePositionSequence( + coordinate_system=CoordinateSystemNames.SLIDE, + image_position=origin_coordinate, + pixel_matrix_position=(1, 1) + ) + ) + else: + plane_positions = [ + PlanePositionSequence( + coordinate_system=CoordinateSystemNames.PATIENT, + image_position=img.ImagePositionPatient + ) + for img in images + ] + + return plane_positions + + def get_index_position(self, pointer: str) -> int: + """Get relative position of a given dimension in the dimension index. + + Parameters + ---------- + pointer: str + Name of the dimension (keyword of the attribute), + e.g., ``"ReferencedSegmentNumber"`` + + Returns + ------- + int + Zero-based relative position + + Examples + -------- + >>> dimension_index = DimensionIndexSequence("SLIDE") + >>> i = dimension_index.get_index_position("ReferencedSegmentNumber") + >>> dimension_description = dimension_index[i] + >>> dimension_description + (0020, 9164) Dimension Organization UID ... + (0020, 9165) Dimension Index Pointer AT: (0062, 000b) + (0020, 9167) Functional Group Pointer AT: (0062, 000a) + (0020, 9421) Dimension Description Label LO: 'Segment Number' + + """ + indices = [ + i + for i, indexer in enumerate(self) + if indexer.DimensionIndexPointer == tag_for_keyword(pointer) + ] + if len(indices) == 0: + raise ValueError( + f'Dimension index does not contain a dimension "{pointer}".' + ) + return indices[0] + + def get_index_values( + self, + plane_positions: Sequence[PlanePositionSequence], + image_orientation: Sequence[float] | None = None, + index_convention: ( + str | + Sequence[PixelIndexDirections | str] + ) = ( + PixelIndexDirections.R, + PixelIndexDirections.D, + ), + handedness: AxisHandedness | str = AxisHandedness.RIGHT_HANDED, + ) -> tuple[np.ndarray, np.ndarray]: + """Get values of indexed attributes that specify position of planes. + + Parameters + ---------- + plane_positions: Sequence[highdicom.PlanePositionSequence] + Plane position of frames in a multi-frame image or in a series of + single-frame images. + image_orientation: Union[Sequence[float], None], optional + An image orientation to use to order frames within a 3D coordinate + system. By default (if ``image_orientation`` is ``None``), the + plane positions are ordered using their raw numerical values and + not along any particular spatial vector. If ``image_orientation`` + is provided, planes are ordered along the positive direction of the + vector normal to the specified. Should be a sequence of 6 floats. + This is only valid when plane position inputs contain only the + ImagePositionPatient. + index_convention: Sequence[Union[highdicom.enum.PixelIndexDirections, str]], optional + Convention used to determine how to order frames if + ``image_orientation`` is specified. Should be a sequence of two + :class:`highdicom.enum.PixelIndexDirections` or their string + representations, giving in order, the indexing conventions used for + specifying pixel indices. For example ``('R', 'D')`` means that the + first pixel index indexes the columns from left to right, and the + second pixel index indexes the rows from top to bottom (this is the + convention typically used within DICOM). As another example ``('D', + 'R')`` would switch the order of the indices to give the convention + typically used within NumPy. + + Alternatively, a single shorthand string may be passed that combines + the string representations of the two directions. So for example, + passing ``'RD'`` is equivalent to passing ``('R', 'D')``. + + This is used in combination with the ``handedness`` to determine + the positive direction used to order frames. + handedness: Union[highdicom.enum.AxisHandedness, str], optional + Choose the frame order such that the frame axis creates a + coordinate system with this handedness in the when combined with + the within-frame convention given by ``index_convention``. + + Returns + ------- + dimension_index_values: numpy.ndarray + Array of dimension index values. The first dimension corresponds + to the items in the input plane_positions sequence. The second + dimension corresponds to the dimensions of the dimension index. + The third dimension (if any) corresponds to the multiplicity + of the values, and is omitted if this is 1 for all dimensions. + plane_indices: numpy.ndarray + 1D array of planes indices for sorting frames according to their + spatial position specified by the dimension index + + Note + ---- + Includes only values of indexed attributes that specify the spatial + position of planes relative to the total pixel matrix or the frame of + reference, and excludes values of the Referenced Segment Number + attribute. + + """ # noqa: E501 + if self._coordinate_system is None: + raise RuntimeError( + 'Cannot calculate index values for multiple plane ' + 'positions when images do not exist within a frame of ' + 'reference.' + ) + + # For each spatial dimension obtain the value of the attribute that the + # Dimension Index Pointer points to in the element of the Plane + # Position Sequence or Plane Position Slide Sequence. Per definition, + # this is the Image Position Patient attribute in case of the patient + # coordinate system, or the X/Y/Z Offset In Slide Coordinate System and + # the Column/Row Position in Total Image Pixel Matrix attributes in + # case of the the slide coordinate system. + indexers = [ + dim_ind for dim_ind in self + if dim_ind.DimensionIndexPointer in _SPATIAL_TAGS + ] + plane_position_values = np.array([ + [ + np.array(p[0][indexer.DimensionIndexPointer].value) + for indexer in indexers + ] + for p in plane_positions + ]) + + if image_orientation is not None: + if not hasattr(plane_positions[0][0], 'ImagePositionPatient'): + raise ValueError( + 'Provided "image_orientation" is only valid when ' + 'plane_positions contain the ImagePositionPatient.' + ) + normal_vector = get_normal_vector( + image_orientation, + index_convention=index_convention, + handedness=handedness, + ) + origin_distances = _get_slice_distances( + plane_position_values[:, 0, :], + normal_vector, + ) + _, plane_sort_indices = np.unique( + origin_distances, + return_index=True, + ) + else: + # Build an array that can be used to sort planes according to the + # Dimension Index Value based on the order of the items in the + # Dimension Index Sequence. + _, plane_sort_indices = np.unique( + plane_position_values, + axis=0, + return_index=True + ) + + if len(plane_sort_indices) != len(plane_positions): + raise ValueError( + 'Input image/frame positions are not unique according to the ' + 'Dimension Index Pointers. The generated segmentation would be ' + 'ambiguous. Ensure that source images/frames have distinct ' + 'locations.' + ) + + return (plane_position_values, plane_sort_indices) + + def get_index_keywords(self) -> list[str]: + """Get keywords of attributes that specify the position of planes. + + Returns + ------- + List[str] + Keywords of indexed attributes + + Note + ---- + Includes only keywords of indexed attributes that specify the spatial + position of planes relative to the total pixel matrix or the frame of + reference, and excludes the keyword of the Referenced Segment Number + attribute. + + Examples + -------- + >>> dimension_index = DimensionIndexSequence('SLIDE') + >>> plane_positions = [ + ... PlanePositionSequence('SLIDE', [10.0, 0.0, 0.0], [1, 1]), + ... PlanePositionSequence('SLIDE', [30.0, 0.0, 0.0], [1, 2]), + ... PlanePositionSequence('SLIDE', [50.0, 0.0, 0.0], [1, 3]) + ... ] + >>> values, indices = dimension_index.get_index_values(plane_positions) + >>> names = dimension_index.get_index_keywords() + >>> for name in names: + ... print(name) + RowPositionInTotalImagePixelMatrix + ColumnPositionInTotalImagePixelMatrix + XOffsetInSlideCoordinateSystem + YOffsetInSlideCoordinateSystem + ZOffsetInSlideCoordinateSystem + >>> index = names.index("XOffsetInSlideCoordinateSystem") + >>> print(values[:, index]) + [10. 30. 50.] + + """ + return [ + keyword_for_tag(indexer.DimensionIndexPointer) + for indexer in self + if indexer.DimensionIndexPointer in _SPATIAL_TAGS + ] + + +class _Image(SOPClass): + + """Base class representing a general DICOM image. + + An "image" is any object representing an Image Information Entity. + + This class serves as a base class for specific image types, including + Segmentations and Parametric Maps, as well as the general Image base class. + + """ + + _coordinate_system: CoordinateSystemNames | None + _is_tiled_full: bool + _has_frame_references: bool + _dim_ind_pointers: list[BaseTag] + # Mapping of tag value to (index column name, val column name(s)) + _dim_ind_col_names: dict[int, tuple[str, str | tuple[str, ...] | None]] + _db_con: sqlite3.Connection + _file_reader: ImageFileReader | None + + @classmethod + def from_dataset( + cls, + dataset: Dataset, + copy: bool = True, + ) -> Self: + """Create an Image from an existing pydicom Dataset. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + copy: bool + If True, the underlying dataset is deep-copied such that the + original dataset remains intact. If False, this operation will + alter the original dataset in place. + + Returns + ------- + Self: + Image object from the input dataset. + + """ + if not isinstance(dataset, Dataset): + raise TypeError( + 'Dataset must be of type pydicom.dataset.Dataset.' + ) + _check_little_endian(dataset) + + # Checks on integrity of input dataset + if copy: + im = deepcopy(dataset) + else: + im = dataset + im.__class__ = cls + im = cast(Self, im) + + im._build_luts() + return im + + @property + def number_of_frames(self) -> int: + """int: Number of frames in the image.""" + return self.get('NumberOfFrames', 1) + + def _get_color_type(self) -> _ImageColorType: + """_ImageColorType: Color type of the image.""" + return _deduce_color_type(self) + + @property + def is_tiled(self) -> bool: + """bool: Whether the image is a tiled multi-frame image.""" + return is_tiled_image(self) + + @property + def coordinate_system(self) -> CoordinateSystemNames | None: + """highdicom.CoordinateSystemNames | None: Frame-of-reference + coordinate system, if any, within which the image exists. + + """ + return self._coordinate_system + + def _init_multiframe_image( + self, + pixel_array: np.ndarray | Volume, + *, + source_images: Sequence[Dataset], + image_type: Sequence[str], + photometric_interpretation: PhotometricInterpretationValues | str, + bits_allocated: int, + functional_groups_module: str, + bits_stored: int | None = None, + samples_per_pixel: int = 1, + planar_configuration: ( + PlanarConfigurationValues | int + ) = PlanarConfigurationValues.COLOR_BY_PIXEL, + pixel_representation: ( + PixelRepresentationValues | int + ) = PixelRepresentationValues.UNSIGNED_INTEGER, + contains_recognizable_visual_features: bool | None = None, + burned_in_annotation: bool | None = None, + use_default_pixel_value_transformation: bool = False, + shared_voi_lut_transformations: ( + Sequence[VOILUTTransformation] | None + ) = None, + palette_color_lut_transformation: PaletteColorLUTTransformation | None, + icc_profile: bytes | None = None, + pixel_measures: PixelMeasuresSequence | None = None, + plane_orientation: PlaneOrientationSequence | None = None, + plane_positions: Sequence[PlanePositionSequence] | None = None, + omit_empty_frames: bool = False, + workers: int | Executor = 0, + dimension_organization_type: ( + DimensionOrganizationTypeValues | + str | + None + ) = None, + tile_pixel_array: bool = False, + tile_size: Sequence[int] | None = None, + further_source_images: Sequence[Dataset] | None = None, + pyramid_label: str | None = None, + pyramid_uid: str | None = None, + use_extended_offset_table: bool = False, + pixel_data_keyword: PixelDataKeywords = PixelDataKeywords.PIXEL_DATA, + channel_values: Sequence[Any] | None = None, # TODO generalize + add_channel_callback: ( + Callable[[Dataset, Any], Dataset] | None + ) = None, # TODO generalize + channel_dimension_index: ChannelDescriptor | None = None, + preprocess_channel_callback: Callable[ + [np.ndarray, Any, int], np.ndarray + ] | None = None, + ): + """ + + Parameters + ---------- + channel_dimension_index: highdicom.ChannelDescriptor | None + If there is a channel dimension and it should be included in the + dimension index sequence, provide the descriptor of the relevant + attribute here. Otherwise, pass ``None``. It is not required that + the channel, if present, be included in the DimensionIndexSequence. + + """ + uniqueness_criteria = { + ( + image.StudyInstanceUID, + image.SeriesInstanceUID, + image.Rows, + image.Columns, + getattr(image, 'FrameOfReferenceUID', None), + ) + for image in source_images + } + if len(uniqueness_criteria) > 1: + raise ValueError( + 'Source images must all be part of the same series and must ' + 'have the same image dimensions (number of rows/columns).' + ) + + src_img = source_images[0] + is_multiframe = is_multiframe_image(src_img) + is_tiled = is_tiled_image(src_img) + if is_multiframe and len(source_images) > 1: + raise ValueError( + 'Only one source image should be provided in case images ' + 'are multi-frame images.' + ) + + # Frame of Reference + has_ref_frame_uid = hasattr(src_img, 'FrameOfReferenceUID') + if has_ref_frame_uid: + self.FrameOfReferenceUID = src_img.FrameOfReferenceUID + self.PositionReferenceIndicator = getattr( + src_img, + 'PositionReferenceIndicator', + None + ) + else: + # Only allow missing FrameOfReferenceUID if it is not required + # for this IOD + usage = get_module_usage('frame-of-reference', src_img.SOPClassUID) + if usage == ModuleUsageValues.MANDATORY: + raise ValueError( + "Source images have no Frame Of Reference UID, but it is " + "required by the IOD." + ) + + self._coordinate_system = get_image_coordinate_system(src_img) + + if self._coordinate_system is None: + # It may be possible to generalize this, but for now only a single + # source frame is permitted when there is no coordinate system + if ( + len(source_images) > 1 or + (is_multiframe and src_img.NumberOfFrames > 1) + ): + raise ValueError( + "Only a single frame is supported when the source " + "image has no Frame of Reference UID." + ) + if plane_positions is not None: + raise TypeError( + "If source images have no Frame Of Reference UID, the " + 'argument "plane_positions" may not be specified since the ' + "segmentation pixel array must be spatially aligned with " + "the source images." + ) + if plane_orientation is not None: + raise TypeError( + "If source images have no Frame Of Reference UID, the " + 'argument "plane_orientation" may not be specified since ' + "the provided pixel array must be spatially aligned " + "with the source images." + ) + + # (Float/DoubleFloat) Image Pixel module + self.BitsAllocated = bits_allocated + self.SamplesPerPixel = samples_per_pixel + if samples_per_pixel > 1: + self.PlanarConfiguration = PlanarConfigurationValues( + planar_configuration + ).value + photometric_interpretation = PhotometricInterpretationValues( + photometric_interpretation + ) + self.PhotometricInterpretation = photometric_interpretation.value + + pixel_data_keyword = PixelDataKeywords(pixel_data_keyword) + + if pixel_data_keyword == PixelDataKeywords.PIXEL_DATA: + # Attributes in the Image Pixel Module but not the + # Float/DoubleFloat Image Pixel modules + if bits_stored is None: + bits_stored = bits_allocated + self.BitsStored = bits_stored + self.HighBit = bits_allocated - 1 + self.PixelRepresentation = PixelRepresentationValues( + pixel_representation + ).value + + # General Image module + self.ImageType = list(image_type) + if burned_in_annotation is not None: + self.BurnedInAnnotation = ( + 'YES' if burned_in_annotation else 'NO' + ) + if contains_recognizable_visual_features is not None: + self.RecognizableVisualFeatures = ( + 'YES' if contains_recognizable_visual_features else 'NO' + ) + self.PresentationLUTShape = ( + 'INVERSE' + if photometric_interpretation == + PhotometricInterpretationValues.MONOCHROME1 + else 'IDENTITY' + ) + + if palette_color_lut_transformation is not None: + _add_palette_color_lookup_table_attributes( + self, + palette_color_lut_transformation, + ) + + if icc_profile is not None: + _add_icc_profile_attributes(self, icc_profile) + + self._add_source_image_references( + source_images=source_images, + further_source_images=further_source_images, + ) + + self._init_pyramid( + pyramid_label=pyramid_label, + pyramid_uid=pyramid_uid, + is_tiled=is_tiled, + ) + + self.LossyImageCompression = getattr( + src_img, + 'LossyImageCompression', + '00' + ) + if self.LossyImageCompression == '01': + if 'LossyImageCompressionRatio' in src_img: + self.LossyImageCompressionRatio = \ + src_img.LossyImageCompressionRatio + if 'LossyImageCompressionMethod' in src_img: + self.LossyImageCompressionMethod = \ + src_img.LossyImageCompressionMethod + + if isinstance(pixel_array, Volume): + ( + pixel_array, + plane_positions, + plane_orientation, + pixel_measures, + ) = self._get_spatial_data_from_volume( + volume=pixel_array, + coordinate_system=get_image_coordinate_system(src_img), + frame_of_reference_uid=src_img.FrameOfReferenceUID, + plane_positions=plane_positions, + plane_orientation=plane_orientation, + pixel_measures=pixel_measures, + ) + + pixel_array = cast(np.ndarray, pixel_array) + + if pixel_array.ndim == 2: + pixel_array = pixel_array[np.newaxis, ...] + if pixel_array.ndim not in [3, 4]: + raise ValueError('Pixel array must be a 2D, 3D, or 4D array.') + + if (channel_values is None) != (add_channel_callback is None): + raise TypeError( + "Argument 'add_channel_callback' should be provided if and " + "only if 'channel_values' is provided." + ) + + if tile_pixel_array and not is_tiled: + raise ValueError( + 'When argument "tile_pixel_array" is True, the source image ' + 'must be a tiled image.' + ) + if tile_pixel_array and pixel_array.shape[0] != 1: + raise ValueError( + 'When argument "tile_pixel_array" is True, the input pixel ' + 'array must contain only one "frame" representing the ' + 'entire pixel matrix.' + ) + + # Image Pixel + if tile_pixel_array: + # By default use the same tile size as the source image (even if + # they are not spatially aligned) + tile_size = tile_size or (src_img.Rows, src_img.Columns) + self.Rows, self.Columns = (tile_size) + else: + self.Rows = pixel_array.shape[1] + self.Columns = pixel_array.shape[2] + + # Multi-frame Dimension + channel_is_indexed = channel_dimension_index is not None + self.DimensionIndexSequence = DimensionIndexSequence( + coordinate_system=self._coordinate_system, + functional_groups_module=functional_groups_module, + channel_dimensions=( + [channel_dimension_index] if channel_is_indexed else None + ), + ) + dimension_organization = Dataset() + dimension_organization.DimensionOrganizationUID = ( + self.DimensionIndexSequence[0].DimensionOrganizationUID + ) + self.DimensionOrganizationSequence = [dimension_organization] + + ( + plane_positions, + plane_orientation, + pixel_measures, + plane_position_values, + plane_sort_index, + derivation_source_image_items, + ) = self._prepare_spatial_metadata( + plane_positions=plane_positions, + plane_orientation=plane_orientation, + pixel_measures=pixel_measures, + source_images=source_images, + further_source_images=further_source_images or [], + tile_pixel_array=tile_pixel_array, + tile_size=tile_size, + frame_shape=pixel_array.shape[1:3], + number_of_planes=pixel_array.shape[0], + dimension_organization_type=dimension_organization_type, + ) + + # Shared functional groops + sffg_item = Dataset() + if ( + self._coordinate_system is not None and + self._coordinate_system == CoordinateSystemNames.PATIENT + ): + sffg_item.PlaneOrientationSequence = plane_orientation + + # Automatically populate the spacing between slices in the + # pixel measures if it was not provided. This is done on the + # initial plane positions, before any removals, to give the + # receiver more information about how to reconstruct a volume + # from the frames in the case that slices are omitted + if 'SpacingBetweenSlices' not in pixel_measures[0]: + ori = plane_orientation[0].ImageOrientationPatient + slice_spacing, _ = get_volume_positions( + image_positions=plane_position_values[:, 0, :], + image_orientation=ori, + ) + if slice_spacing is not None: + pixel_measures[0].SpacingBetweenSlices = ( + format_number_as_ds(slice_spacing) + ) + + if pixel_measures is not None: + sffg_item.PixelMeasuresSequence = pixel_measures + + # If there is only a single channel value, add it to the shared + # functional groups + if channel_values is not None and len(channel_values) == 1: + sffg_item = add_channel_callback(sffg_item, channel_values[0]) + + # Set to None so this is not used for per-frame items later + add_channel_callback = None + + if use_default_pixel_value_transformation: + # Identity Pixel Value Transformation + transformation_item = Dataset() + transformation_item.RescaleIntercept = 0 + transformation_item.RescaleSlope = 1 + transformation_item.RescaleType = 'US' + sffg_item.PixelValueTransformationSequence = [transformation_item] + + if shared_voi_lut_transformations is not None: + sffg_item.FrameVOILUTSequence = shared_voi_lut_transformations + + self.SharedFunctionalGroupsSequence = [sffg_item] + + # Find indices such that empty planes are removed + if omit_empty_frames: + if tile_pixel_array: + included_plane_indices, is_empty = \ + self._get_nonempty_tile_indices( + pixel_array, + plane_positions=plane_positions, + rows=self.Rows, + columns=self.Columns, + ) + else: + included_plane_indices, is_empty = \ + self._get_nonempty_plane_indices(pixel_array) + if is_empty: + # Cannot omit empty frames when all frames are empty + omit_empty_frames = False + included_plane_indices = range(len(plane_positions)) + else: + # Remove all empty plane positions from the list of sorted + # plane position indices + included_plane_indices_set = set(included_plane_indices) + plane_sort_index = [ + ind for ind in plane_sort_index + if ind in included_plane_indices_set + ] + else: + included_plane_indices = range(len(plane_positions)) + + # Dimension Organization Type + dimension_organization_type = self._check_tiled_dimension_organization( + dimension_organization_type=dimension_organization_type, + is_tiled=is_tiled, + omit_empty_frames=omit_empty_frames, + plane_positions=plane_positions, + tile_pixel_array=tile_pixel_array, + rows=self.Rows, + columns=self.Columns, + num_channels=1 if channel_values is None else len(channel_values), + ) + + if ( + self._coordinate_system is not None and + dimension_organization_type != + DimensionOrganizationTypeValues.TILED_FULL + ): + # Get unique values of attributes in the Plane Position Sequence or + # Plane Position Slide Sequence, which define the position of the + # plane with respect to the three dimensional patient or slide + # coordinate system, respectively. These can subsequently be used + # to look up the relative position of a plane relative to the + # indexed dimension. + unique_dimension_values = [ + np.unique( + plane_position_values[included_plane_indices, index], + axis=0 + ) + for index in range(plane_position_values.shape[1]) + ] + else: + unique_dimension_values = [None] + + if self._coordinate_system == CoordinateSystemNames.PATIENT: + inferred_dim_org_type = None + + # To be considered "3D", an segmentation should have frames that + # are differentiated only by location. This rules out any image + # with channels where there is more than a single channel value + # Further, only images with multiple spatial positions in the final + # image should be considered to have 3D dimension organization type + if ( + len(included_plane_indices) > 1 and + ( + channel_values is None or + len(channel_values) == 1 + ) + ): + # Calculate the spacing using only the included planes, and + # enforce ordering + ori = plane_orientation[0].ImageOrientationPatient + spacing, _ = get_volume_positions( + image_positions=plane_position_values[ + included_plane_indices, 0, : + ], + image_orientation=ori, + sort=False, + ) + if spacing is not None and spacing > 0.0: + inferred_dim_org_type = ( + DimensionOrganizationTypeValues.THREE_DIMENSIONAL + ) + + if ( + dimension_organization_type == + DimensionOrganizationTypeValues.THREE_DIMENSIONAL + ) and inferred_dim_org_type is None: + raise ValueError( + 'Dimension organization "3D" has been specified, ' + 'but the source image is not a regularly-spaced 3D ' + 'volume.' + ) + dimension_organization_type = inferred_dim_org_type + + if dimension_organization_type is not None: + self.DimensionOrganizationType = dimension_organization_type.value + + if ( + self._coordinate_system is not None and + self._coordinate_system == CoordinateSystemNames.SLIDE + ): + plane_position_names = ( + self.DimensionIndexSequence.get_index_keywords() + ) + row_dim_index = plane_position_names.index( + 'RowPositionInTotalImagePixelMatrix' + ) + col_dim_index = plane_position_names.index( + 'ColumnPositionInTotalImagePixelMatrix' + ) + + is_encaps = self.file_meta.TransferSyntaxUID.is_encapsulated + process_pool: Executor | None = None + + if not isinstance(workers, (int, Executor)): + raise TypeError( + 'Argument "workers" must be of type int or ' + 'concurrent.futures.Executor (or a derived class).' + ) + using_multiprocessing = ( + isinstance(workers, Executor) or workers != 0 + ) + + # List of frames. In the case of native transfer syntaxes, we will + # collect a list of frames as flattened NumPy arrays for bitpacking at + # the end. In the case of encapsulated transfer syntaxes with no + # workers, we will accumulate a list of encoded frames to encapsulate + # at the end + frames: list[bytes] | list[np.ndarray] = [] + + # In the case of native encoding when the number pixels in a frame is + # not a multiple of 8. This array carries "leftover" pixels that + # couldn't be encoded in previous iterations, to future iterations. This + # saves having to keep the entire un-endoded array in memory, which can + # get extremely heavy on memory in the case of very large arrays + remainder_pixels = np.empty((0, ), dtype=np.uint8) + + if is_encaps: + if using_multiprocessing: + # In the case of encapsulated transfer syntaxes with multiple + # workers, we will accumulate a list of encoded frames to + # encapsulate at the end + frame_futures: list[Future] = [] + + # Use the existing executor or create one + if isinstance(workers, Executor): + process_pool = workers + else: + # If workers is negative, pass None to use all processors + process_pool = ProcessPoolExecutor( + workers if workers > 0 else None + ) + + # Parameters to use when calling the encode_frame function in + # either of the above two cases + encode_frame_kwargs = dict( + transfer_syntax_uid=self.file_meta.TransferSyntaxUID, + bits_allocated=self.BitsAllocated, + bits_stored=self.BitsStored, + photometric_interpretation=self.PhotometricInterpretation, + pixel_representation=self.PixelRepresentation + ) + else: + if using_multiprocessing: + warnings.warn( + "Setting workers != 0 or passing an instance of " + "concurrent.futures.Executor when using a non-encapsulated " + "transfer syntax has no effect.", + UserWarning, + stacklevel=2, + ) + using_multiprocessing = False + + # Information about individual frames is placed into the + # PerFrameFunctionalGroupsSequence. Note that a *very* significant + # efficiency gain is observed when building this as a Python list + # rather than a pydicom sequence, and then converting to a pydicom + # sequence at the end + pffg_sequence: list[Dataset] = [] + + if channel_values is None: + if self._coordinate_system is None: + # If there are no channels nor any spatial indoces, we need to + # create some dimension to index down + channel_values = ["Frame"] + + # Function to use as callback to add the segment identification + def _add_frame_label( + pffg_item: Dataset, + frame_label: str, + ): + pffg_item.FrameContentSequence[0].add( + DataElement( + 0x00209453, # FrameLabel + 'LO', + frame_label, + ) + ) + return pffg_item + + add_channel_callback = _add_frame_label + + def preprocess_channel_callback( + arr: np.ndarray, + _: Any, # value (unused0) + __: int, # index (unused) + ): + return arr + + channel_is_indexed = True + else: + # Placeholder used just to have something to iterate over to + # get a single loop over plane positions + channel_values = [None] + + if preprocess_channel_callback is None: + # Default channel callback just indexes down the final dimension + if pixel_array.ndim == 4: + def preprocess_channel_callback( + arr: np.ndarray, + _: Any, + ind: int + ): + return arr[:, :, ind] + else: + def preprocess_channel_callback( + arr: np.ndarray, + _: Any, + __: int, + ): + return arr + + # Main frame loop: encode frames and create per-frame functional groups + for channel_index, channel_value in enumerate(channel_values): + + for plane_dim_ind, plane_index in enumerate(plane_sort_index, 1): + + if tile_pixel_array: + if ( + dimension_organization_type == + DimensionOrganizationTypeValues.TILED_FULL + ): + row_offset = int( + plane_position_values[plane_index, row_dim_index] + ) + column_offset = int( + plane_position_values[plane_index, col_dim_index] + ) + else: + pos = plane_positions[plane_index][0] + row_offset = pos.RowPositionInTotalImagePixelMatrix + column_offset = ( + pos.ColumnPositionInTotalImagePixelMatrix + ) + + plane_array = get_tile_array( + pixel_array[0], + row_offset=row_offset, + column_offset=column_offset, + tile_rows=self.Rows, + tile_columns=self.Columns, + ) + else: + # Select the relevant existing frame + plane_array = pixel_array[plane_index] + + if channel_value is None: + # Deal with all segments at once + channel_array = plane_array + else: + channel_array = preprocess_channel_callback( + plane_array, + channel_value, + channel_index, + ) + + # Even though completely empty planes were removed earlier, + # there may still be planes in which this specific segment is + # absent. Such frames should be removed + if omit_empty_frames and not np.any(channel_array): + logger.debug( + f'skip empty plane {plane_index} of channel ' + f'#{channel_index}' + ) + continue + + # Log a debug message + if channel_value is None: + msg = f'add plane #{plane_index}' + else: + msg = ( + f'add plane #{plane_index} for channel ' + f'#{channel_value}' + ) + logger.debug(msg) + + if ( + dimension_organization_type != + DimensionOrganizationTypeValues.TILED_FULL + ): + # No per-frame functional group for TILED FULL - """ - self.table_name = table_name - self.column_defs = list(column_defs) + # Get the item of the PerFrameFunctionalGroupsSequence for + # this segmentation frame + if self._coordinate_system is not None: + plane_pos_val = plane_position_values[plane_index] + if ( + self._coordinate_system == + CoordinateSystemNames.SLIDE + ): + try: + dimension_index_values = [ + int( + np.where( + unique_dimension_values[idx] == pos + )[0][0] + 1 + ) + for idx, pos in enumerate(plane_pos_val) + ] + except IndexError as error: + raise IndexError( + 'Could not determine position of plane ' + f'#{plane_index} in three dimensional ' + 'coordinate system based on dimension ' + f'index values: {error}' + ) from error + else: + dimension_index_values = [plane_dim_ind] + else: + dimension_index_values = [] - # It is important to convert numpy arrays to lists, otherwise the - # values get inserted into the table as binary values and this leads to - # very difficult to detect failures - if isinstance(column_data, np.ndarray): - column_data = column_data.tolist() + plane_derivation_items = ( + derivation_source_image_items[plane_index] + if derivation_source_image_items is not None else None + ) + if channel_value is not None and channel_is_indexed: + dimension_index_values = ( + [int(channel_index) + 1] + + dimension_index_values + ) - sanitized_column_data = [] + pffg_item = self._get_pffg_item( + dimension_index_values=dimension_index_values, + plane_position=plane_positions[plane_index], + derivation_source_image_items=plane_derivation_items, + coordinate_system=self._coordinate_system, + ) - for col in column_data: - if isinstance(col, np.ndarray): - col = col.tolist() + if add_channel_callback is not None: + # Add the channel information to the per-frame item + pffg_item = add_channel_callback( + pffg_item, + channel_value, + ) - sanitized_column_data.append(col) + pffg_sequence.append(pffg_item) + + # Add the segmentation pixel array for this frame to the list + if is_encaps: + if process_pool is None: + # Encode this frame and add resulting bytes to the list + # for encapsulation at the end + frames.append( + encode_frame( + channel_array, + **encode_frame_kwargs, + ) + ) + else: + # Submit this frame for encoding this frame and add the + # future to the list for encapsulation at the end + future = process_pool.submit( + encode_frame, + array=channel_array, + **encode_frame_kwargs, + ) + frame_futures.append(future) + else: + flat_array = channel_array.flatten() + if ( + self.BitsAllocated == 1 and + (self.Rows * self.Columns) // 8 != 0 + ): + # Need to encode a multiple of 8 pixels at a time + full_array = np.concatenate( + [remainder_pixels, flat_array] + ) + # Round down to closest multiple of 8 + n_pixels_to_take = 8 * (len(full_array) // 8) + to_encode = full_array[:n_pixels_to_take] + remainder_pixels = full_array[n_pixels_to_take:] + else: + # Simple - each frame can be individually encoded + to_encode = flat_array - self.column_data = sanitized_column_data + frames.append(self._encode_pixels_native(to_encode)) + if ( + dimension_organization_type != + DimensionOrganizationTypeValues.TILED_FULL + ): + self.PerFrameFunctionalGroupsSequence = pffg_sequence -class _Image(SOPClass): + if is_encaps: + if process_pool is not None: + frames = [ + fut.result() for fut in frame_futures + ] - """Base class representing a general DICOM image. + # Shutdown the pool if we created it, otherwise it is the + # caller's responsibility + if process_pool is not workers: + process_pool.shutdown() - An "image" is any object representing an Image Information Entity. + # Encapsulate all pre-compressed frames + self.NumberOfFrames = len(frames) + if use_extended_offset_table: + ( + self.PixelData, + self.ExtendedOffsetTable, + self.ExtendedOffsetTableLengths, + ) = encapsulate_extended(frames) + else: + # Encapsulated pixels are always put in PixelData attribute + self.PixelData = encapsulate(frames) + else: + self.NumberOfFrames = len(frames) - This class serves as a base class for specific image types, including - Segmentations and Parametric Maps, as well as the general Image base class. + # May need to add in a final set of pixels + if len(remainder_pixels) > 0: + frames.append(self._encode_pixels_native(remainder_pixels)) - """ + setattr(self, pixel_data_keyword.value, b''.join(frames)) - _coordinate_system: CoordinateSystemNames | None - _is_tiled_full: bool - _has_frame_references: bool - _dim_ind_pointers: list[BaseTag] - # Mapping of tag value to (index column name, val column name(s)) - _dim_ind_col_names: dict[int, tuple[str, str | tuple[str, ...] | None]] - _db_con: sqlite3.Connection - _file_reader: ImageFileReader | None + # Add a null trailing byte if required (can't happen for floating pixel + # data) + if 'PixelData' in self and len(self.PixelData) % 2 == 1: + self.PixelData += b'0' - @classmethod - def from_dataset( - cls, - dataset: Dataset, - copy: bool = True, - ) -> Self: - """Create an Image from an existing pydicom Dataset. + # Build lookup tables for efficient decoding + self._build_luts() + + def _encode_pixels_native(self, planes: np.ndarray) -> bytes: + """Encode pixel planes using a native transfer syntax. Parameters ---------- - dataset: pydicom.Dataset - Dataset representing an image. - copy: bool - If True, the underlying dataset is deep-copied such that the - original dataset remains intact. If False, this operation will - alter the original dataset in place. + planes: numpy.ndarray + Array representing one or more segmentation image planes. If + multiple image planes, planes stacked down the first dimension + (index 0). Returns ------- - Self: - Image object from the input dataset. + bytes + Encoded pixels """ - if not isinstance(dataset, Dataset): - raise TypeError( - 'Dataset must be of type pydicom.dataset.Dataset.' - ) - _check_little_endian(dataset) - - # Checks on integrity of input dataset - if copy: - im = deepcopy(dataset) + if self.BitsAllocated == 1: + return pack_bits(planes, pad=False) else: - im = dataset - im.__class__ = cls - im = cast(Self, im) - - im._build_luts() - return im - - @property - def number_of_frames(self) -> int: - """int: Number of frames in the image.""" - return self.get('NumberOfFrames', 1) - - def _get_color_type(self) -> _ImageColorType: - """_ImageColorType: Color type of the image.""" - return _deduce_color_type(self) + return planes.tobytes() - @property - def is_tiled(self) -> bool: - """bool: Whether the image is a tiled multi-frame image.""" - return is_tiled_image(self) + def _init_pyramid( + self, + pyramid_label: str | None, + pyramid_uid: str | None, + is_tiled: bool, + ) -> None: + # Multi-Resolution Pyramid + if pyramid_uid is not None: + if not is_tiled: + raise TypeError( + 'Argument "pyramid_uid" should only be specified ' + 'for tiled images.' + ) + if ( + self._coordinate_system is None or + self._coordinate_system != CoordinateSystemNames.SLIDE + ): + raise TypeError( + 'Argument "pyramid_uid" should only be specified for ' + 'segmentations in the SLIDE coordinate system.' + ) + self.PyramidUID = pyramid_uid - @property - def coordinate_system(self) -> CoordinateSystemNames | None: - """highdicom.CoordinateSystemNames | None: Frame-of-reference - coordinate system, if any, within which the image exists. + if pyramid_label is not None: + _check_long_string(pyramid_label) + self.PyramidLabel = pyramid_label - """ - return self._coordinate_system + elif pyramid_label is not None: + raise TypeError( + 'Argument "pyramid_label" should not be specified if ' + '"pyramid_uid" is not specified.' + ) @staticmethod def _get_pixel_measures_sequence( @@ -2188,6 +3646,110 @@ def _prepare_tiled_spatial_metadata( plane_sort_index, ) + @staticmethod + def _get_nonempty_plane_indices( + pixel_array: np.ndarray + ) -> tuple[list[int], bool]: + """Get a list of all indices of original planes that are non-empty. + + Empty planes (without any positive pixels in any of the channels) do + not need to be included in the image. This method finds a + list of indices of the input frames that are non-empty, and therefore + should be included in the image. + + Parameters + ---------- + pixel_array: numpy.ndarray + Segmentation pixel array + + Returns + ------- + included_plane_indices : List[int] + List giving for each plane position in the resulting segmentation + image the index of the corresponding frame in the original pixel + array. + is_empty: bool + Whether the entire image is empty. If so, empty frames should not + be omitted. + + """ + # This list tracks which source image each non-empty frame came from + source_image_indices = [ + i for i, frm in enumerate(pixel_array) + if np.any(frm) + ] + + if len(source_image_indices) == 0: + logger.warning( + 'Encoding an empty image with "omit_empty_frames" ' + 'set to True. Reverting to encoding all frames since omitting ' + 'all frames is not possible.' + ) + return (list(range(pixel_array.shape[0])), True) + + return (source_image_indices, False) + + @staticmethod + def _get_nonempty_tile_indices( + pixel_array: np.ndarray, + plane_positions: Sequence[PlanePositionSequence], + rows: int, + columns: int, + ) -> tuple[list[int], bool]: + """Get a list of all indices of tile locations that are non-empty. + + This is similar to _get_nonempty_plane_indices, but works on a total + pixel matrix rather than a set of frames. Empty planes (without any + positive pixels in any of the channels) do not need to be included in + the image. This method finds a list of indices of the input frames that + are non-empty, and therefore should be included in the image. + + Parameters + ---------- + pixel_array: numpy.ndarray + Pixel array + plane_positions: Sequence[highdicom.PlanePositionSequence] + Plane positions of each tile. + rows: int + Number of rows in each tile. + columns: int + Number of columns in each tile. + + Returns + ------- + included_plane_indices : List[int] + List giving for each plane position in the resulting segmentation + image the index of the corresponding frame in the original pixel + array. + is_empty: bool + Whether the entire image is empty. If so, empty frames should not + be omitted. + + """ + # This list tracks which source image each non-empty frame came from + source_image_indices = [ + i for i, pos in enumerate(plane_positions) + if np.any( + get_tile_array( + pixel_array[0], + row_offset=pos[0].RowPositionInTotalImagePixelMatrix, + column_offset=pos[0].ColumnPositionInTotalImagePixelMatrix, + tile_rows=rows, + tile_columns=columns, + ) + ) + ] + + if len(source_image_indices) == 0: + logger.warning( + 'Encoding an empty image with "omit_empty_frames" ' + 'set to True. Reverting to encoding all frames since omitting ' + 'all frames is not possible.' + ) + return (list(range(len(plane_positions))), True) + + return (source_image_indices, False) + @staticmethod def _get_derivation_source_image_item( source_image: Dataset, @@ -2261,6 +3823,118 @@ def _get_derivation_source_image_item( ) return derivation_src_img_item + @staticmethod + def _get_pffg_item( + dimension_index_values: list[int], + plane_position: PlanePositionSequence, + derivation_source_image_items: Sequence[Dataset] | None, + coordinate_system: CoordinateSystemNames | None, + ) -> Dataset: + """Get a single item of the Per Frame Functional Groups Sequence. + + This is a helper method used in the common multiframe constructor. + + Parameters + ---------- + dimension_index_values: List[int] + Dimension index values for this frame. + plane_position: highdicom.seg.PlanePositionSequence + Plane position of this frame. + derivation_source_image_items: Sequence[pydicom.Dataset] | None + Items of the Source Image Sequence, to place into the Derivation + Image Sequence, if any, for this frame. + coordinate_system: Optional[highdicom.CoordinateSystemNames] + Coordinate system used, if any. + + Returns + ------- + pydicom.Dataset + Dataset representing the item of the Per Frame Functional Groups + Sequence for this frame frame. + + """ + # NB this function is called many times in a loop when there are a + # large number of frames, and has been observed to dominate the + # creation time of some segmentations. Therefore we use low-level + # pydicom primitives to improve performance as much as possible + pffg_item = Dataset() + frame_content_item = Dataset() + + frame_content_item.add( + DataElement( + 0x00209157, # DimensionIndexValues + 'UL', + dimension_index_values, + ) + ) + + pffg_item.add( + DataElement( + 0x00209111, # FrameContentSequence + 'SQ', + [frame_content_item] + ) + ) + + if coordinate_system is not None: + if coordinate_system == CoordinateSystemNames.SLIDE: + pffg_item.add( + DataElement( + 0x0048021a, # PlanePositionSlideSequence + 'SQ', + plane_position + ) + ) + else: + pffg_item.add( + DataElement( + 0x00209113, # PlanePositionSequence + 'SQ', + plane_position + ) + ) + + if ( + derivation_source_image_items is not None and + len(derivation_source_image_items) > 0 + ): + derivation_image_item = Dataset() + derivation_image_item.add( + DataElement( + 0x00089215, # DerivationCodeSequence + 'SQ', + [_DERIVATION_CODE] + ) + ) + + derivation_image_item.add( + DataElement( + 0x00082112, # SourceImageSequence + 'SQ', + derivation_source_image_items, + ) + ) + pffg_item.add( + DataElement( + 0x00089124, # DerivationImageSequence + 'SQ', + [derivation_image_item] + ) + ) + else: + # Determining the source images that map to the frame is not + # always trivial. Since DerivationImageSequence is a type 2 + # attribute, we leave its value empty. + pffg_item.add( + DataElement( + 0x00089124, # DerivationImageSequence + 'SQ', + [] + ) + ) + + return pffg_item + def _add_source_image_references( self, source_images: Sequence[Dataset], @@ -2450,8 +4124,8 @@ def _add_slide_coordinate_metadata( format_number_as_ds(z_center) self.ImageCenterPointCoordinatesSequence = [center_item] - @staticmethod def _check_tiled_dimension_organization( + self, dimension_organization_type: ( DimensionOrganizationTypeValues | str | @@ -2463,6 +4137,7 @@ def _check_tiled_dimension_organization( tile_pixel_array: bool, rows: int, columns: int, + num_channels: int, ) -> DimensionOrganizationTypeValues | None: """Checks that the specified Dimension Organization Type is valid. @@ -2482,6 +4157,8 @@ def _check_tiled_dimension_organization( Number of rows in each frame of the segmentation image. columns: int Number of columns in each frame of the segmentation image. + num_channels: + Number of channels included in this image. Returns ------- @@ -2505,10 +4182,10 @@ def _check_tiled_dimension_organization( dimension_organization_type = DimensionOrganizationTypeValues( dimension_organization_type ) - tiled_dimension_organization_types = [ + tiled_dimension_organization_types = ( DimensionOrganizationTypeValues.TILED_SPARSE, - DimensionOrganizationTypeValues.TILED_FULL - ] + DimensionOrganizationTypeValues.TILED_FULL, + ) if ( dimension_organization_type in @@ -2541,8 +4218,24 @@ def _check_tiled_dimension_organization( 'because the "plane_positions" of the segmentation ' 'do not follow the relevant requirements. See ' 'https://dicom.nema.org/medical/dicom/current/output/' - 'chtml/part03/sect_C.7.6.17.3.html#sect_C.7.6.17.3 .' + 'chtml/part03/sect_C.7.6.17.3.html#sect_C.7.6.17.3.' + ) + + # A TILED_FULL image can only contain nultiple channels if they + # use certain values. Specifically multiple segments or + # multiple optical paths (but multiple real world value maps + # for a parametric map are NOT allowed). This check needs to be + # generalized but the following works for the types of channels + # currently supported (segments and RWVMs) + if num_channels > 1 and self.SOPClassUID != SegmentationStorage: + raise ValueError( + 'A value of "TILED_FULL" for parameter ' + '"dimension_organization_type" is not permitted ' + 'because the image contains multiple channels. See ' + 'https://dicom.nema.org/medical/dicom/current/output/' + 'chtml/part03/sect_C.7.6.17.3.html#sect_C.7.6.17.3.' ) + if omit_empty_frames: raise ValueError( 'Parameter "omit_empty_frames" should be False if ' @@ -2663,7 +4356,14 @@ def get_raw_frame( start = frame_index * frame_length end = start + frame_length - return self.PixelData[start:end] + if 'PixelData' in self: + return self.PixelData[start:end] + elif 'FloatPixelData' in self: + return self.FloatPixelData[start:end] + elif 'DoubleFloatPixelData' in self: + return self.DoubleFloatPixelData[start:end] + else: + raise AttributeError('Dataset contains no pixel data.') def get_stored_frame( self, @@ -3284,6 +4984,32 @@ def _serialize_db(self) -> bytes: ] ) + @staticmethod + def _standardize_lut_value(elem: DataElement): + """Standardize a value of a data element to store in a LUT. + + Most values are returned unchanged, however sequences are assumed to + represent codes and the code is returned as a single string. + + Parameters + ---------- + elem: pydicom.DataElement + Data Element + + Returns + ------- + Any: + Value from the element suitable for storing in SQL table. + + """ + v = elem.value + + # Interpret a sequence as a code + if isinstance(v, DataElementSequence): + return v[0].CodingSchemeDesignator + v[0].CodeValue + + return v + def _build_luts(self) -> None: """Build lookup tables for efficient querying. @@ -3464,6 +5190,8 @@ def _build_luts_multiframe(self) -> None: (0x0020_9111, 0x0020_9056), # FrameContentSequence/InStackPositionNumber (0x0020_9111, 0x0020_9057), + # RealWorldValueMappingSequence/LUTLabel + (0x0040_9096, 0x0040_9210), ]: if ptr in self._dim_ind_pointers: # Skip if this attribute is already indexed due to being a @@ -3487,7 +5215,7 @@ def _build_luts_multiframe(self) -> None: found = True # Get the shared value - dim_val = grp_seq[ptr].value + dim_val = self._standardize_lut_value(grp_seq[ptr]) # Check whether the attribute is in the first per-frame functional # group. If so, assume that it is there for all per-frame functional @@ -3619,9 +5347,13 @@ def _build_luts_multiframe(self) -> None: grp_ptr = func_grp_pointers[ptr] if grp_ptr is not None: - dim_val = frame_item[grp_ptr][0][ptr].value + dim_val = self._standardize_lut_value( + frame_item[grp_ptr][0][ptr] + ) else: - dim_val = frame_item[ptr].value + dim_val = self._standardize_lut_value( + frame_item[ptr] + ) dim_values[ptr].append(dim_val) for ptr in extra_collection_pointers: @@ -3635,9 +5367,13 @@ def _build_luts_multiframe(self) -> None: grp_ptr = extra_collection_func_pointers[ptr] if grp_ptr is not None: - dim_val = frame_item[grp_ptr][0][ptr].value + dim_val = self._standardize_lut_value( + frame_item[grp_ptr][0][ptr] + ) else: - dim_val = frame_item[ptr].value + dim_val = self._standardize_lut_value( + frame_item[ptr] + ) extra_collection_values[ptr].append(dim_val) for der_im in getattr( @@ -4892,7 +6628,7 @@ def _generate_temp_tables( for tdef in table_defs: # Clean up the tables - cmd = (f'DROP TABLE {tdef.table_name}') + cmd = f'DROP TABLE {tdef.table_name}' with self._db_con: self._db_con.execute(cmd) @@ -4930,7 +6666,7 @@ def _generate_temp_view( yield # Clean up the view - cmd = (f'DROP VIEW {view_name}') + cmd = f'DROP VIEW {view_name}' with self._db_con: self._db_con.execute(cmd) diff --git a/src/highdicom/pm/__init__.py b/src/highdicom/pm/__init__.py index 8c32115b..8162d09c 100644 --- a/src/highdicom/pm/__init__.py +++ b/src/highdicom/pm/__init__.py @@ -2,6 +2,7 @@ from highdicom.pm.content import DimensionIndexSequence, RealWorldValueMapping from highdicom.pm.enum import DerivedPixelContrastValues, ImageFlavorValues +from highdicom.pm.pyramid import create_parametric_map_pyramid from highdicom.pm.sop import ParametricMap SOP_CLASS_UIDS = { @@ -14,4 +15,5 @@ 'ImageFlavorValues', 'ParametricMap', 'RealWorldValueMapping', + 'create_parametric_map_pyramid', ] diff --git a/src/highdicom/pm/content.py b/src/highdicom/pm/content.py index 6c5d46d6..0f632288 100644 --- a/src/highdicom/pm/content.py +++ b/src/highdicom/pm/content.py @@ -1,19 +1,15 @@ from collections.abc import Sequence +import warnings import numpy as np -from pydicom.datadict import keyword_for_tag, tag_for_keyword from pydicom.dataset import Dataset -from pydicom.sequence import Sequence as DataElementSequence from pydicom.sr.coding import Code -from highdicom._module_utils import is_multiframe_image -from highdicom.content import PlanePositionSequence from highdicom.enum import CoordinateSystemNames from highdicom.pixels import apply_lut from highdicom.sr.coding import CodedConcept from highdicom.sr.value_types import CodeContentItem -from highdicom.uid import UID -from highdicom.utils import compute_plane_position_slide_per_frame +from highdicom.image import DimensionIndexSequence as BaseDimensionIndexSequence from highdicom.valuerep import ( _check_long_string, _check_short_string, @@ -51,22 +47,22 @@ def __init__( should be restricted. For example, values may be stored as floating-point values with double precision, but limited to the range ``(-1.0, 1.0)`` or ``(0.0, 1.0)`` or stored as 16-bit - unsigned integer values but limited to range ``(0, 4094). - Note that the type of the values in `value_range` is significant + unsigned integer values but limited to range ``(0, 4094)``. + Note that the type of the values in ``value_range`` is significant and is used to determine whether values are stored as integers or floating-point values. Therefore, use ``(0.0, 1.0)`` instead of ``(0, 1)`` to specify a range of floating-point values. slope: Union[int, float, None], optional Slope of the linear mapping function applied to values in - `value_range`. + ``value_range``. intercept: Union[int, float, None], optional Intercept of the linear mapping function applied to values in - `value_range`. + ``value_range``. lut_data: Union[Sequence[int], Sequence[float], None], optional Sequence of values to serve as a lookup table for mapping stored values into real-world values in case of a non-linear relationship. The sequence should contain an entry for each value in the specified - `value_range` such that + ``value_range`` such that ``len(sequence) == value_range[1] - value_range[0] + 1``. For example, in case of a value range of ``(0, 255)``, the sequence shall have ``256`` entries - one for each value in the given range. @@ -77,9 +73,9 @@ def __init__( Note ---- - Either `slope` and `intercept` or `lut_data` must be specified. - Specify `slope` and `intercept` if the mapping can be described by a - linear function. Specify `lut_data` if the relationship between stored + Either ``slope`` and ``intercept`` or ``lut_data`` must be specified. + Specify ``slope` and ``intercept``` if the mapping can be described by a + linear function. Specify ``lut_data`` if the relationship between stored and real-world values is non-linear. Note, however, that a non-linear relationship can only be described for values that are stored as integers. Values stored as floating-point numbers must map linearly to @@ -173,10 +169,15 @@ def lut_data(self) -> np.ndarray | None: return np.array(self.RealWorldValueLUTData) return None + def is_floating_point(self) -> bool: + """bool: Whether the value range is defined with floating point + values.""" + return 'DoubleFloatRealWorldValueFirstValueMapped' in self + @property def value_range(self) -> tuple[float, float]: """Tuple[float, float]: Range of valid input values.""" - if 'DoubleFloatRealWorldValueFirstValueMapped' in self: + if self.is_floating_point(): return ( self.DoubleFloatRealWorldValueFirstValueMapped, self.DoubleFloatRealWorldValueLastValueMapped, @@ -238,281 +239,51 @@ def apply( return array * slope + intercept -class DimensionIndexSequence(DataElementSequence): +class DimensionIndexSequence(BaseDimensionIndexSequence): """Sequence of data elements describing dimension indices for the patient or slide coordinate system based on the Dimension Index functional group macro. + Note ---- The order of indices is fixed. + + Note + ---- + This class is deprecated and will be removed in a future version of + highdicom. User code should generally avoid this class, and if necessary, + the more general :class:`highdicom.DimensionIndexSequence` should be used + instead. + """ def __init__( self, - coordinate_system: str | CoordinateSystemNames + coordinate_system: str | CoordinateSystemNames | None, ) -> None: """ Parameters ---------- - coordinate_system: Union[str, highdicom.CoordinateSystemNames] + coordinate_system: Union[str, highdicom.CoordinateSystemNames, None] Subject (``"PATIENT"`` or ``"SLIDE"``) that was the target of - imaging - """ - super().__init__() - dim_uid = UID() - - self._coordinate_system = CoordinateSystemNames(coordinate_system) - if self._coordinate_system == CoordinateSystemNames.SLIDE: - x_axis_index = Dataset() - x_axis_index.DimensionIndexPointer = tag_for_keyword( - 'XOffsetInSlideCoordinateSystem' - ) - x_axis_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - x_axis_index.DimensionOrganizationUID = dim_uid - x_axis_index.DimensionDescriptionLabel = \ - 'X Offset in Slide Coordinate System' - - y_axis_index = Dataset() - y_axis_index.DimensionIndexPointer = tag_for_keyword( - 'YOffsetInSlideCoordinateSystem' - ) - y_axis_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - y_axis_index.DimensionOrganizationUID = dim_uid - y_axis_index.DimensionDescriptionLabel = \ - 'Y Offset in Slide Coordinate System' - - z_axis_index = Dataset() - z_axis_index.DimensionIndexPointer = tag_for_keyword( - 'ZOffsetInSlideCoordinateSystem' - ) - z_axis_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - z_axis_index.DimensionOrganizationUID = dim_uid - z_axis_index.DimensionDescriptionLabel = \ - 'Z Offset in Slide Coordinate System' - - row_dimension_index = Dataset() - row_dimension_index.DimensionIndexPointer = tag_for_keyword( - 'ColumnPositionInTotalImagePixelMatrix' - ) - row_dimension_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - row_dimension_index.DimensionOrganizationUID = dim_uid - row_dimension_index.DimensionDescriptionLabel = \ - 'Column Position In Total Image Pixel Matrix' - - column_dimension_index = Dataset() - column_dimension_index.DimensionIndexPointer = tag_for_keyword( - 'RowPositionInTotalImagePixelMatrix' - ) - column_dimension_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - column_dimension_index.DimensionOrganizationUID = dim_uid - column_dimension_index.DimensionDescriptionLabel = \ - 'Row Position In Total Image Pixel Matrix' - - # Organize frames for each segment similar to TILED_FULL, first - # along the row dimension (column indices from left to right) and - # then along the column dimension (row indices from top to bottom) - # of the Total Pixel Matrix. - self.extend([ - row_dimension_index, - column_dimension_index, - x_axis_index, - y_axis_index, - z_axis_index, - ]) - - elif self._coordinate_system == CoordinateSystemNames.PATIENT: - image_position_index = Dataset() - image_position_index.DimensionIndexPointer = tag_for_keyword( - 'ImagePositionPatient' - ) - image_position_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSequence' - ) - image_position_index.DimensionOrganizationUID = dim_uid - image_position_index.DimensionDescriptionLabel = \ - 'Image Position Patient' - - self.append(image_position_index) - - else: - raise ValueError( - f'Unknown coordinate system "{self._coordinate_system}"' - ) - - def get_plane_positions_of_image( - self, - image: Dataset - ) -> list[PlanePositionSequence]: - """Get plane positions of frames in multi-frame image. - - Parameters - ---------- - image: Dataset - Multi-frame image - - Returns - ------- - List[highdicom.PlanePositionSequence] - Plane position of each frame in the image + imaging. If None, the imaging does not belong within a frame of + reference. + include_segment_number: bool + Include the segment number as a dimension index. """ - is_multiframe = is_multiframe_image(image) - if not is_multiframe: - raise ValueError('Argument "image" must be a multi-frame image.') - - if self._coordinate_system == CoordinateSystemNames.SLIDE: - if hasattr(image, 'PerFrameFunctionalGroupsSequence'): - plane_positions = [ - item.PlanePositionSlideSequence - for item in image.PerFrameFunctionalGroupsSequence - ] - else: - # If Dimension Organization Type is TILED_FULL, plane - # positions are implicit and need to be computed. - plane_positions = compute_plane_position_slide_per_frame( - image - ) - else: - plane_positions = [ - item.PlanePositionSequence - for item in image.PerFrameFunctionalGroupsSequence - ] - - return plane_positions - - def get_plane_positions_of_series( - self, - images: Sequence[Dataset] - ) -> list[PlanePositionSequence]: - """Gets plane positions for series of single-frame images. - - Parameters - ---------- - images: Sequence[Dataset] - Series of single-frame images - - Returns - ------- - List[highdicom.PlanePositionSequence] - Plane position of each frame in the image - - """ - is_multiframe = any([is_multiframe_image(img) for img in images]) - if is_multiframe: - raise ValueError( - 'Argument "images" must be a series of single-frame images.' - ) - - plane_positions = [ - PlanePositionSequence( - coordinate_system=CoordinateSystemNames.PATIENT, - image_position=img.ImagePositionPatient - ) - for img in images - ] - - return plane_positions - - def get_index_values( - self, - plane_positions: Sequence[PlanePositionSequence] - ) -> tuple[np.ndarray, np.ndarray]: - """Get the values of indexed attributes. - - Parameters - ---------- - plane_positions: Sequence[highdicom.PlanePositionSequence] - Plane position of frames in a multi-frame image or in a series of - single-frame images - - Returns - ------- - dimension_index_values: numpy.ndarray - 2D array of spatial dimension index values - plane_indices: numpy.ndarray - 1D array of planes indices for sorting frames according to their - spatial position specified by the dimension index. - - """ - # For each indexed spatial dimension obtain the value of the attribute - # that the Dimension Index Pointer points to in the element of the - # Plane Position Sequence or Plane Position Slide Sequence. - # In case of the patient coordinate system, this is the Image Position - # Patient attribute. In case of the slide coordinate system, these are - # X/Y/Z Offset In Slide Coordinate System and the Column/Row - # Position in Total Image Pixel Matrix attributes. - plane_position_values = np.array([ - [ - np.array(p[0][indexer.DimensionIndexPointer].value) - for indexer in self - ] - for p in plane_positions - ]) - - # Build an array that can be used to sort planes according to the - # Dimension Index Value based on the order of the items in the - # Dimension Index Sequence. - _, plane_sort_indices = np.unique( - plane_position_values, - axis=0, - return_index=True + warnings.warn( + "The highdicom.pm.DimensionIndexSequence class is deprecated and " + "will be removed in a future version of the library. User code " + "should typically avoid this class, or, if required, use the more " + "general highdicom.DimensionIndexSequence instead.", + UserWarning, + stacklevel=2, + ) + super().__init__( + coordinate_system=coordinate_system, + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), ) - - return (plane_position_values, plane_sort_indices) - - def get_index_position(self, pointer: str) -> int: - """Get relative position of a given dimension in the dimension index. - - Parameters - ---------- - pointer: str - Name of the dimension (keyword of the attribute), - e.g., ``"XOffsetInSlideCoordinateSystem"`` - - Returns - ------- - int - Zero-based relative position - - Examples - -------- - >>> dimension_index = DimensionIndexSequence("SLIDE") - >>> i = dimension_index.get_index_position("XOffsetInSlideCoordinateSystem") - >>> x_offsets = dimension_index[i] - - """ # noqa: E501 - indices = [ - i - for i, indexer in enumerate(self) - if indexer.DimensionIndexPointer == tag_for_keyword(pointer) - ] - if len(indices) == 0: - raise ValueError( - f'Dimension index does not contain a dimension "{pointer}".' - ) - return indices[0] - - def get_index_keywords(self) -> list[str]: - """Get keywords of attributes that specify the position of planes. - - Returns - ------- - List[str] - Keywords of indexed attributes - - """ - return [ - keyword_for_tag(indexer.DimensionIndexPointer) - for indexer in self - ] diff --git a/src/highdicom/pm/pyramid.py b/src/highdicom/pm/pyramid.py new file mode 100644 index 00000000..55d1287a --- /dev/null +++ b/src/highdicom/pm/pyramid.py @@ -0,0 +1,226 @@ +"""Tools for constructing multi-resolution parametric map pyramids.""" +from typing import Any +from collections.abc import Sequence + +import numpy as np +from pydicom import Dataset + +from highdicom.content import VOILUTTransformation +from highdicom.enum import InterpolationMethods +from highdicom._pyramid import iter_derived_pyramid_levels +from highdicom.pm.sop import ParametricMap +from highdicom.pm.content import ( + RealWorldValueMapping, +) +from highdicom.uid import UID + + +def create_parametric_map_pyramid( + source_images: Sequence[Dataset], + pixel_arrays: Sequence[np.ndarray], + interpolator: InterpolationMethods | str, + series_instance_uid: str | None, + series_number: int, + manufacturer: str, + manufacturer_model_name: str, + software_versions: str | tuple[str], + device_serial_number: str, + contains_recognizable_visual_features: bool, + real_world_value_mappings: ( + Sequence[RealWorldValueMapping] | + Sequence[Sequence[RealWorldValueMapping]] + ), + voi_lut_transformations: ( + Sequence[VOILUTTransformation] | None + ) = None, + downsample_factors: Sequence[float] | None = None, + sop_instance_uids: list[str] | None = None, + pyramid_uid: str | None = None, + pyramid_label: str | None = None, + **kwargs: Any +) -> list[ParametricMap]: + """Construct a multi-resolution parametric map pyramid series. + + A multi-resolution pyramid represents the same parametric map array at + multiple resolutions. + + This function handles multiple related scenarios: + + * Constructing a parametric map of a source image pyramid given a + parametric map pixel array of the highest resolution source image. + Highdicom performs the downsampling automatically to match the resolution + of the other source images. For this case, pass multiple + ``source_images`` and a single item in ``pixel_arrays``. + * Constructing a parametric map of a source image pyramid given + user-provided parametric map pixel arrays for each level in the source + pyramid. For this case, pass multiple ``source_images`` and a matching + number of ``pixel_arrays``. + * Constructing a parametric map of a single source image given multiple + user-provided downsampled parametric map pixel arrays. For this case, + pass a single item in ``source_images``, and multiple items in + ``pixel_arrays``). + * Constructing a parametric map of a single source image and a single + parametric map pixel array by downsampling by a given list of + ``downsample_factors``. For this case, pass a single item in + ``source_images``, a single item in ``pixel_arrays``, and a list of one + or more desired ``downsample_factors``. + + In all cases, the items in both ``source_images`` and ``pixel_arrays`` + should be sorted in pyramid order from highest resolution (smallest + spacing) to lowest resolution (largest spacing), and the pixel array + in ``pixel_arrays[0]`` must be the parametric map of the source image in + ``source_images[0]`` with spatial locations preserved (a one-to-one + correspondence between pixels in the source image's total pixel matrix and + the provided parametric map pixel array). + + In all cases, the provided pixel arrays should be total pixel matrices. + Tiling is performed automatically. + + Parameters + ---------- + source_images: Sequence[pydicom.Dataset] + List of source images. If there are multiple source images, they should + be from the same series and pyramid. + pixel_arrays: Sequence[numpy.ndarray] + List of parametric maps pixel arrays. Each should be a total pixel + matrix, i.e. have shape (rows, columns), (1, rows, columns), or (1, + rows, columns, channels). Otherwise all options supported by the + constructor of :class:`highdicom.pm.ParametricMap` are permitted. + interpolator: highdicom.InterpolationMethods | str, optional + Interpolation method used for the downsampling operations. + series_instance_uid: str | None + UID for the output parametric series. If ``None``, a UID will be + generated using highdicom's prefix. + series_number: int + Number of the output parametric map series. + manufacturer: str + Name of the manufacturer of the device (developer of the software) + that creates the instance + manufacturer_model_name: str + Name of the device model (name of the software library or + application) that creates the instance + software_versions: Union[str, Tuple[str]] + Version(s) of the software that creates the instance. + device_serial_number: str + Manufacturer's serial number of the device + contains_recognizable_visual_features: bool + Whether the image contains recognizable visible features of the + patient + real_world_value_mappings: Union[Sequence[highdicom.pm.RealWorldValueMapping], Sequence[Sequence[highdicom.pm.RealWorldValueMapping]] + Descriptions of how stored values map to real-world values. Each + channel encoded in each item of ``pixel_arrays`` shall be described + with one or more real-world value mappings. Multiple mappings might be + used for different representations such as log versus linear scales or + for different representations in different units. If `pixel_array` is a + 2D or 3D array and only one channel exists at each spatial image + position, then one or more real-world value mappings shall be provided + in a flat sequence. If `pixel_array` is a 4D array and multiple + channels exist at each spatial image position, then one or more + mappings shall be provided for each channel in a nested sequence of + length ``m``, where ``m`` shall match the channel dimension of each + item of ``pixel_arrays``. + + In some situations the mapping may be difficult to describe (e.g., in + case of a transformation performed by a deep convolutional neural + network). The real-world value mapping may then simply describe an + identity function that maps stored values to unit-less real-world + values. + voi_lut_transformations: Sequence[highdicom.VOILUTTransformation] | None, optional + One or more VOI transformations that describe a pixel transformation to + apply to frames. + downsample_factors: Optional[Sequence[float]], optional + Factors by which to downsample the pixel array to create each of the + output parametric map objects. This should be provided if and only if a + single source image and single pixel array are provided. Note that the + original array is always used to create the first parametric map + output, so the number of created segmententation instances is one + greater than the number of items in this list. Items must be numbers + greater than 1 and sorted in ascending order. A downsampling factor of + *n* implies that the output array is *1/n* time the size of input pixel + array. For example a list ``[2, 4, 8]`` would be produce 4 output + parametric map instances. The first is the same size as the original + pixel array, the next is half the size, the next is a quarter of the + size of the original, and the last is one eighth the size of the + original. Output sizes are rounded to the nearest integer. + sop_instance_uids: Optional[List[str]], optional + SOP instance UIDS of the output instances. If not specified, UIDs are + generated automatically using highdicom's prefix. + pyramid_uid: Optional[str], optional + UID for the output imaging pyramid. If not specified, a UID is generated + using highdicom's prefix. + pyramid_label: Optional[str], optional + A human readable label for the output pyramid. + **kwargs: Any + Any further parameters are passed directly to the constructor of the + :class:`highdicom.pm.ParametricMap` object. However the following + parameters are disallowed: ``instance_number``, ``sop_instance_uid``, + ``plane_orientation``, ``plane_positions``, ``pixel_measures``, + ``pixel_array``, ``tile_pixel_array``. + + """ # noqa: E501 + # Disallow duplicate items in kwargs + kwarg_keys = set(kwargs.keys()) + disallowed_keys = { + 'instance_number', + 'sop_instance_uid', + 'plane_orientation', + 'plane_positions', + 'pixel_array', + 'tile_pixel_array', + } + error_keys = kwarg_keys & disallowed_keys + if len(error_keys) > 0: + raise TypeError( + f'kwargs supplied to the create_parametric_map_pyramid function ' + f'should not contain a value for parameter {list(error_keys)[0]}.' + ) + + if pyramid_uid is None: + pyramid_uid = UID() + if series_instance_uid is None: + series_instance_uid = UID() + + all_pmaps = [] + + for ( + source_image, + pixel_array, + pixel_measures, + instance_number, + sop_instance_uid, + ) in iter_derived_pyramid_levels( + source_images=source_images, + pixel_arrays=pixel_arrays, + interpolator=interpolator, + downsample_factors=downsample_factors, + sop_instance_uids=sop_instance_uids, + ): + # Create the output segmentation + pmap = ParametricMap( + source_images=[source_image], + pixel_array=pixel_array, + series_instance_uid=series_instance_uid, + series_number=series_number, + sop_instance_uid=sop_instance_uid, + instance_number=instance_number, + manufacturer=manufacturer, + manufacturer_model_name=manufacturer_model_name, + software_versions=software_versions, + device_serial_number=device_serial_number, + contains_recognizable_visual_features=( + contains_recognizable_visual_features + ), + real_world_value_mappings=real_world_value_mappings, + voi_lut_transformations=voi_lut_transformations, + pyramid_uid=pyramid_uid, + pyramid_label=pyramid_label, + tile_pixel_array=True, + plane_orientation=None, + plane_positions=None, + pixel_measures=pixel_measures, + **kwargs, + ) + + all_pmaps.append(pmap) + + return all_pmaps diff --git a/src/highdicom/pm/sop.py b/src/highdicom/pm/sop.py index 645cb0d3..4e12f097 100644 --- a/src/highdicom/pm/sop.py +++ b/src/highdicom/pm/sop.py @@ -1,31 +1,33 @@ -from collections import defaultdict -from copy import deepcopy -from typing import cast +"""Module for SOP classes of the PM modality.""" from collections.abc import Sequence -from enum import Enum +from concurrent.futures import Executor +from os import PathLike +from typing import cast, BinaryIO +import warnings import numpy as np -from pydicom.encaps import encapsulate, encapsulate_extended -from highdicom.base import SOPClass from highdicom.base_content import ContributingEquipment from highdicom.content import ( + _add_content_information, ContentCreatorIdentificationCodeSequence, PaletteColorLUTTransformation, PixelMeasuresSequence, PlaneOrientationSequence, PlanePositionSequence, + VOILUTTransformation, ) -from highdicom.enum import ContentQualificationValues, CoordinateSystemNames -from highdicom.frame import encode_frame -from highdicom.pm.content import DimensionIndexSequence, RealWorldValueMapping -from highdicom.pm.enum import DerivedPixelContrastValues, ImageFlavorValues -from highdicom.valuerep import ( - check_person_name, - _check_code_string, - _check_long_string, +from highdicom.enum import ( + ContentQualificationValues, + DimensionOrganizationTypeValues, + PhotometricInterpretationValues, + PixelDataKeywords, ) -from highdicom._module_utils import is_multiframe_image +from highdicom.image import _Image, Image +from highdicom.pm.content import RealWorldValueMapping +from highdicom.pm.enum import DerivedPixelContrastValues, ImageFlavorValues +from highdicom.volume import ChannelDescriptor, Volume from pydicom import Dataset +from pydicom.dataelem import DataElement from pydicom.uid import ( UID, ExplicitVRLittleEndian, @@ -34,18 +36,10 @@ JPEGLSLossless, RLELossless, ) -from pydicom.valuerep import format_number_as_ds - +from typing_extensions import Self -class _PixelDataType(Enum): - """Helper class for tracking the type of the pixel data.""" - USHORT = 1 - SINGLE = 2 - DOUBLE = 3 - - -class ParametricMap(SOPClass): +class ParametricMap(Image): """SOP class for a Parametric Map. @@ -60,7 +54,7 @@ class ParametricMap(SOPClass): def __init__( self, source_images: Sequence[Dataset], - pixel_array: np.ndarray, + pixel_array: np.ndarray | Volume, series_instance_uid: str, series_number: int, sop_instance_uid: str, @@ -74,8 +68,11 @@ def __init__( Sequence[RealWorldValueMapping] | Sequence[Sequence[RealWorldValueMapping]] ), - window_center: float, - window_width: float, + window_center: float | None = None, + window_width: float | None = None, + voi_lut_transformations: ( + Sequence[VOILUTTransformation] | None + ) = None, transfer_syntax_uid: str | UID = ExplicitVRLittleEndian, content_description: str | None = None, content_creator_name: str | None = None, @@ -102,6 +99,18 @@ def __init__( ContributingEquipment ] | None = None, use_extended_offset_table: bool = False, + icc_profile: bytes | None = None, + workers: int | Executor = 0, + dimension_organization_type: ( + DimensionOrganizationTypeValues | + str | + None + ) = None, + tile_pixel_array: bool = False, + tile_size: Sequence[int] | None = None, + pyramid_uid: str | None = None, + pyramid_label: str | None = None, + further_source_images: Sequence[Dataset] | None = None, **kwargs, ): """ @@ -111,7 +120,7 @@ def __init__( source_images: Sequence[pydicom.dataset.Dataset] One or more single- or multi-frame images (or metadata of images) from which the parametric map was derived - pixel_array: numpy.ndarray + pixel_array: numpy.ndarray | highdicom.Volume 2D, 3D, or 4D array of unsigned integer or floating-point data type representing one or more channels (images derived from source images via an image transformation) for one or more spatial image @@ -157,19 +166,19 @@ def __init__( contains_recognizable_visual_features: bool Whether the image contains recognizable visible features of the patient - real_world_value_mappings: Union[Sequence[highdicom.map.RealWorldValueMapping], Sequence[Sequence[highdicom.map.RealWorldValueMapping]] - Descriptions of how stored values map to real-world values. - Each channel encoded in `pixel_array` shall be described with one - or more real-world value mappings. Multiple mappings might be - used for different representations such as log versus linear scales - or for different representations in different units. - If `pixel_array` is a 2D or 3D array and only one channel exists - at each spatial image position, then one or more real-world value - mappings shall be provided in a flat sequence. - If `pixel_array` is a 4D array and multiple channels exist at each - spatial image position, then one or more mappings shall be provided - for each channel in a nested sequence of length ``m``, where ``m`` - shall match the channel dimension of the `pixel_array``. + real_world_value_mappings: Union[Sequence[highdicom.pm.RealWorldValueMapping], Sequence[Sequence[highdicom.pm.RealWorldValueMapping]] + Descriptions of how stored values map to real-world values. Each + channel encoded in ``pixel_array`` shall be described with one or + more real-world value mappings. Multiple mappings might be used for + different representations such as log versus linear scales or for + different representations in different units. If ``pixel_array`` is + a 2D or 3D array and only one channel exists at each spatial image + position, then one or more real-world value mappings shall be + provided in a flat sequence. If `pixel_array` is a 4D array and + multiple channels exist at each spatial image position, then one or + more mappings shall be provided for each channel in a nested + sequence of length ``m``, where ``m`` shall match the channel + dimension of the ``pixel_array``. In some situations the mapping may be difficult to describe (e.g., in case of a transformation performed by a deep convolutional neural @@ -177,23 +186,17 @@ def __init__( identity function that maps stored values to unit-less real-world values. window_center: Union[int, float, None], optional - Window center (intensity) for rescaling stored values for display - purposes by applying a linear transformation function. For example, - in case of floating-point values in the range ``[0.0, 1.0]``, the - window center may be ``0.5``, in case of floating-point values in - the range ``[-1.0, 1.0]`` the window center may be ``0.0``, in case - of unsigned integer values in the range ``[0, 255]`` the window - center may be ``128``. + This argument has been deprecated and will be removed in a future + release. Use the more flexible ``voi_lut_transformations`` argument + instead. window_width: Union[int, float, None], optional - Window width (contrast) for rescaling stored values for display - purposes by applying a linear transformation function. For example, - in case of floating-point values in the range ``[0.0, 1.0]``, the - window width may be ``1.0``, in case of floating-point values in the - range ``[-1.0, 1.0]`` the window width may be ``2.0``, and in - case of unsigned integer values in the range ``[0, 255]`` the - window width may be ``256``. In case of unbounded floating-point - values, a sensible window width should be chosen to allow for - stored values to be displayed on 8-bit monitors. + This argument has been deprecated and will be removed in a future + release. Use the more flexible ``voi_lut_transformations`` argument + instead. + voi_lut_transformations: Sequence[highdicom.VOILUTTransformation] | None, optional + One or more VOI transformations that describe a pixel + transformation to apply to frames. This will become a required + argument in a future release. transfer_syntax_uid: Union[str, None], optional UID of transfer syntax that should be used for encoding of data elements. Defaults to Explicit VR Little Endian @@ -207,7 +210,7 @@ def __init__( If ``None``, it will be assumed that the parametric map image has the same pixel measures as the source image(s). plane_orientation: Union[highdicom.PlaneOrientationSequence, None], optional - Orientation of planes in `pixel_array` relative to axes of + Orientation of planes in ``pixel_array`` relative to axes of three-dimensional patient or slide coordinate space. If ``None``, it will be assumed that the parametric map image as the same plane orientation as the source image(s). @@ -255,61 +258,35 @@ def __init__( ------ ValueError When - * Length of `source_images` is zero. - * Items of `source_images` are not all part of the same study + + * Length of ``source_images`` is zero. + * Items of ``source_images`` are not all part of the same study and series. - * Items of `source_images` have different number of rows and + * Items of ``source_images`` have different number of rows and columns. - * Length of `plane_positions` does not match number of 2D planes + * Length of ``plane_positions`` does not match number of 2D planes in `pixel_array` (size of first array dimension). - * Transfer Syntax specified by `transfer_syntax_uid` is not + * Transfer Syntax specified by ``transfer_syntax_uid`` is not supported for data type of `pixel_array`. Note ---- - The assumption is made that planes in `pixel_array` are defined in - the same frame of reference as `source_images`. It is further assumed - that all image frame have the same type (i.e., the same `image_flavor` - and `derived_pixel_contrast`). + The assumption is made that planes in ``pixel_array`` are defined in + the same frame of reference as ``source_images``. It is further assumed + that all image frame have the same type (i.e., the same ``image_flavor`` + and ``derived_pixel_contrast``). """ # noqa if len(source_images) == 0: - raise ValueError('At least one source image is required') - self._source_images = source_images - - uniqueness_criteria = { - ( - image.StudyInstanceUID, - image.SeriesInstanceUID, # TODO: Might be overly restrictive - image.Rows, - image.Columns, - image.FrameOfReferenceUID, - ) - for image in self._source_images - } - if len(uniqueness_criteria) > 1: - raise ValueError( - 'Source images must all be part of the same series and must' - 'have the same image dimensions (number of rows/columns).' - ) + raise ValueError('At least one source image is required.') - src_img = self._source_images[0] - is_multiframe = is_multiframe_image(src_img) - # TODO: Revisit, may be overly restrictive - # Check Source Image Sequence attribute in General Reference module - if is_multiframe: - if len(self._source_images) > 1: - raise ValueError( - 'Only one source image should be provided in case images ' - 'are multi-frame images.' - ) - self._src_num_frames = src_img.NumberOfFrames + src_img = source_images[0] supported_transfer_syntaxes = { ImplicitVRLittleEndian, ExplicitVRLittleEndian, } - if pixel_array.dtype.kind == 'u': + if np.dtype(pixel_array.dtype).kind == 'u': # If pixel data has unsigned or signed integer data type, then it # can be lossless compressed. The standard does not specify any # compression codecs for floating-point data types. @@ -325,22 +302,50 @@ def __init__( f'Transfer syntax "{transfer_syntax_uid}" is not supported.' ) - if window_width <= 0: - raise ValueError('Window width must be greater than zero.') - - if pixel_array.ndim == 2: - pixel_array = pixel_array[np.newaxis, ...] + if (window_center is None) != (window_width is None): + raise TypeError( + "Arguments 'window_center' and 'window_width' should both " + "be None, or neither should be None." + ) + if window_center is not None: + if voi_lut_transformations is not None: + raise TypeError( + "Arguments 'window_center' and 'window_width' must be " + "omitted if 'voi_lut_transformations' is provided." + ) + warnings.warn( + "Arguments 'window_center' and 'window_width' are deprecated " + "and will be removed in a future version of the library. " + "Use the more flexible 'voi_lut_transformations' argument " + "instead.", + UserWarning, + stacklevel=2, + ) + voi_lut_transformations = [ + VOILUTTransformation( + window_center=window_center, + window_width=window_width, + ) + ] + else: + if voi_lut_transformations is None: + raise TypeError( + "Argument 'voi_lut_transformations' is required." + ) + if len(voi_lut_transformations) < 1: + raise TypeError( + "Argument 'voi_lut_transformations' must contain at least " + 'one item.' + ) - # There are different DICOM Attributes in the SOP instance depending - # on what type of data is being saved. This lets us keep track of that - # a bit easier - self._pixel_data_type_map = { - _PixelDataType.USHORT: 'PixelData', - _PixelDataType.SINGLE: 'FloatPixelData', - _PixelDataType.DOUBLE: 'DoubleFloatPixelData', - } + for v in voi_lut_transformations: + if not isinstance(v, VOILUTTransformation): + raise TypeError( + "Argument 'voi_lut_transformations' must be a " + 'sequence of highdicom.VOILUTTransformation objects.' + ) - super().__init__( + super(_Image, self).__init__( study_instance_uid=src_img.StudyInstanceUID, series_instance_uid=series_instance_uid, series_number=series_number, @@ -367,106 +372,9 @@ def __init__( **kwargs, ) - if hasattr(src_img, 'ImageOrientationSlide') or hasattr( - src_img, 'ImageCenterPointCoordinatesSequence' - ): - coordinate_system = CoordinateSystemNames.SLIDE - else: - coordinate_system = CoordinateSystemNames.PATIENT - - # Frame of Reference - self.FrameOfReferenceUID = src_img.FrameOfReferenceUID - self.PositionReferenceIndicator = getattr( - src_img, 'PositionReferenceIndicator', None - ) - - # General Reference - self.SourceImageSequence: list[Dataset] = [] - referenced_series: dict[str, list[Dataset]] = defaultdict(list) - for s_img in self._source_images: - ref = Dataset() - ref.ReferencedSOPClassUID = s_img.SOPClassUID - ref.ReferencedSOPInstanceUID = s_img.SOPInstanceUID - self.SourceImageSequence.append(ref) - referenced_series[s_img.SeriesInstanceUID].append(ref) - - # Common Instance Reference - self.ReferencedSeriesSequence: list[Dataset] = [] - for ( - series_instance_uid, - referenced_images, - ) in referenced_series.items(): - ref = Dataset() - ref.SeriesInstanceUID = series_instance_uid - ref.ReferencedInstanceSequence = referenced_images - self.ReferencedSeriesSequence.append(ref) - - # Parametric Map Image - image_flavor = ImageFlavorValues(image_flavor) - derived_pixel_contrast = DerivedPixelContrastValues( - derived_pixel_contrast - ) - self.ImageType = [ - "DERIVED", - "PRIMARY", - image_flavor.value, - derived_pixel_contrast.value, - ] - content_qualification = ContentQualificationValues( - content_qualification - ) - self.ContentQualification = content_qualification.value - self.LossyImageCompression = getattr( - src_img, 'LossyImageCompression', '00' - ) - if self.LossyImageCompression == "01": - self.LossyImageCompressionRatio = ( - src_img.LossyImageCompressionRatio - ) - self.LossyImageCompressionMethod = ( - src_img.LossyImageCompressionMethod - ) - self.SamplesPerPixel = 1 - self.PhotometricInterpretation = 'MONOCHROME2' - self.BurnedInAnnotation = 'NO' - if contains_recognizable_visual_features: - self.RecognizableVisualFeatures = 'YES' - else: - self.RecognizableVisualFeatures = 'NO' - - if content_label is not None: - _check_code_string(content_label) - self.ContentLabel = content_label - else: - self.ContentLabel = 'MAP' - if content_description is not None: - _check_long_string(content_description) - self.ContentDescription = content_description - if content_creator_name is not None: - check_person_name(content_creator_name) - self.ContentCreatorName = content_creator_name - if content_creator_identification is not None: - if not isinstance( - content_creator_identification, - ContentCreatorIdentificationCodeSequence - ): - raise TypeError( - 'Argument "content_creator_identification" must be of type ' - 'ContentCreatorIdentificationCodeSequence.' - ) - self.ContentCreatorIdentificationCodeSequence = \ - content_creator_identification - - self.PresentationLUTShape = 'IDENTITY' - - self.DimensionIndexSequence = DimensionIndexSequence(coordinate_system) - dimension_organization = Dataset() - dimension_organization.DimensionOrganizationUID = ( - self.DimensionIndexSequence[0].DimensionOrganizationUID - ) - self.DimensionOrganizationSequence = [dimension_organization] - - sffg_item = Dataset() + self.copy_specimen_information(src_img) + self.copy_patient_and_study_information(src_img) + self._add_contributing_equipment(contributing_equipment, src_img) # If the same set of mappings applies to all frames, the information # is stored in the Shared Functional Groups Sequence. Otherwise, it @@ -478,7 +386,6 @@ def __init__( '"real_world_value_mappings" must be a flat sequence ' 'of one or more RealWorldValueMapping items.' ) - sffg_item.RealWorldValueMappingSequence = real_world_value_mappings try: real_world_value_mappings[0] except IndexError as e: @@ -488,10 +395,6 @@ def __init__( RealWorldValueMapping ): raise TypeError(error_message) - if pixel_array.ndim == 2: - pixel_array = pixel_array[np.newaxis, ..., np.newaxis] - elif pixel_array.ndim == 3: - pixel_array = pixel_array[..., np.newaxis] real_world_value_mappings = cast( Sequence[Sequence[RealWorldValueMapping]], [real_world_value_mappings] @@ -502,10 +405,11 @@ def __init__( '"real_world_value_mappings" must be a nested sequence ' 'of one or more RealWorldValueMapping items.' ) - try: - real_world_value_mappings[0][0] - except IndexError as e: - raise TypeError(error_message) from e + if isinstance( + real_world_value_mappings[0], + RealWorldValueMapping + ): + raise TypeError(error_message) if not isinstance( real_world_value_mappings[0][0], RealWorldValueMapping @@ -518,292 +422,211 @@ def __init__( else: raise ValueError('Pixel array must be a 2D, 3D, or 4D array.') - # Acquisition Context - self.AcquisitionContextSequence: list[Dataset] = [] + if isinstance(pixel_array, Volume): + if pixel_array.number_of_channel_dimensions == 1: + if pixel_array.channel_descriptors != ( + ChannelDescriptor('LUTLabel'), + ): + raise ValueError( + "Input volume should have no channels other than " + "'LUTLabel'." + ) - # Image Pixel - self.Rows = pixel_array.shape[1] - self.Columns = pixel_array.shape[2] + for channel_mappings, vol_lut_label in zip( + real_world_value_mappings, + pixel_array.get_channel_values('LUTLabel') + ): + if len(channel_mappings) != 1: + raise ValueError( + 'Only a single mapping should be provided in ' + "each item within 'real_world_value_mappings' " + "when a Volume is passed as the 'pixel_array'." + ) + mapping = channel_mappings[0] + + if vol_lut_label != mapping.LUTLabel: + raise ValueError( + "The LUTLabels of the 'real_world_value_mappings' " + "must match those within the channel indentifiers " + "of the 'pixel_array'." + ) + + elif pixel_array.number_of_channel_dimensions != 0: + raise ValueError( + "If 'pixel_array' is a highdicom.Volume, it should have " + "0 or 1 channel dimensions." + ) - if len(real_world_value_mappings) != pixel_array.shape[3]: + n_channels = pixel_array.shape[3] if pixel_array.ndim == 4 else 1 + if len(real_world_value_mappings) != n_channels: raise ValueError( 'Number of RealWorldValueMapping items provided via ' '"real_world_value_mappings" argument does not match size of ' 'last dimension of "pixel_array" argument.' ) - if is_multiframe: - source_plane_positions = \ - self.DimensionIndexSequence.get_plane_positions_of_image( - self._source_images[0] - ) - if coordinate_system == CoordinateSystemNames.SLIDE: - source_plane_orientation = PlaneOrientationSequence( - coordinate_system=coordinate_system, - image_orientation=[ - float(v) for v in src_img.ImageOrientationSlide - ], - ) - else: - src_sfg = src_img.SharedFunctionalGroupsSequence[0] - source_plane_orientation = src_sfg.PlaneOrientationSequence - else: - source_plane_positions = \ - self.DimensionIndexSequence.get_plane_positions_of_series( - self._source_images - ) - source_plane_orientation = PlaneOrientationSequence( - coordinate_system=coordinate_system, - image_orientation=[ - float(v) for v in src_img.ImageOrientationPatient - ], - ) + # Parametric Map Image + image_flavor = ImageFlavorValues(image_flavor) + derived_pixel_contrast = DerivedPixelContrastValues( + derived_pixel_contrast + ) + image_type = [ + "DERIVED", + "PRIMARY", + image_flavor.value, + derived_pixel_contrast.value, + ] + content_qualification = ContentQualificationValues( + content_qualification + ) + self.ContentQualification = content_qualification.value - if plane_positions is None: - if pixel_array.shape[0] != len(source_plane_positions): - raise ValueError( - 'Number of plane positions in source image(s) does not ' - 'match size of first dimension of "pixel_array" argument.' - ) - plane_positions = source_plane_positions - else: - if len(plane_positions) != pixel_array.shape[0]: - raise ValueError( - 'Number of PlanePositionSequence items provided via ' - '"plane_positions" argument does not match size of ' - 'first dimension of "pixel_array" argument.' - ) + _add_content_information( + dataset=self, + content_label=( + content_label if content_label is not None else 'MAP' + ), + content_description=content_description, + content_creator_name=content_creator_name, + content_creator_identification=content_creator_identification, + ) - if plane_orientation is None: - plane_orientation = source_plane_orientation + # Acquisition Context + self.AcquisitionContextSequence: list[Dataset] = [] - are_spatial_locations_preserved = ( - all( - plane_positions[i] == source_plane_positions[i] - for i in range(len(plane_positions)) - ) and - plane_orientation == source_plane_orientation + # Get the correct pixel data attribute + plain_array = ( + pixel_array.array + if isinstance(pixel_array, Volume) + else pixel_array ) + pixel_data_keyword = self._get_pixel_data_keyword(plain_array) + bits_allocated = { + PixelDataKeywords.PIXEL_DATA: int(plain_array.itemsize * 8), + PixelDataKeywords.FLOAT_PIXEL_DATA: 32, + PixelDataKeywords.DOUBLE_FLOAT_PIXEL_DATA: 64, + }[pixel_data_keyword] - plane_position_names = self.DimensionIndexSequence.get_index_keywords() - plane_position_values, plane_sort_index = \ - self.DimensionIndexSequence.get_index_values(plane_positions) - - if coordinate_system == CoordinateSystemNames.SLIDE: - self.ImageOrientationSlide = \ - plane_orientation[0].ImageOrientationSlide - if are_spatial_locations_preserved: - self.TotalPixelMatrixOriginSequence = \ - source_images[0].TotalPixelMatrixOriginSequence - self.TotalPixelMatrixRows = \ - source_images[0].TotalPixelMatrixRows - self.TotalPixelMatrixColumns = \ - source_images[0].TotalPixelMatrixColumns - else: - row_index = plane_position_names.index( - 'RowPositionInTotalImagePixelMatrix' - ) - row_offsets = plane_position_values[:, row_index] - col_index = plane_position_names.index( - 'ColumnPositionInTotalImagePixelMatrix' - ) - col_offsets = plane_position_values[:, col_index] - frame_indices = np.lexsort([row_offsets, col_offsets]) - first_frame_index = frame_indices[0] - last_frame_index = frame_indices[-1] - x_index = plane_position_names.index( - 'XOffsetInSlideCoordinateSystem' - ) - x_offset = plane_position_values[first_frame_index, x_index] - y_index = plane_position_names.index( - 'YOffsetInSlideCoordinateSystem' - ) - y_offset = plane_position_values[first_frame_index, y_index] - origin_item = Dataset() - origin_item.XOffsetInSlideCoordinateSystem = x_offset - origin_item.YOffsetInSlideCoordinateSystem = y_offset - self.TotalPixelMatrixOriginSequence = [origin_item] - self.TotalPixelMatrixRows = int( - plane_position_values[last_frame_index, row_index] + - self.Rows - ) - self.TotalPixelMatrixColumns = int( - plane_position_values[last_frame_index, col_index] + - self.Columns + # Palette color lookup table + self._configure_color( + palette_color_lut_transformation=palette_color_lut_transformation, + icc_profile=icc_profile, + pixel_data_keyword=pixel_data_keyword, + ) + + # Check that the real world value maps are consistent with the provided + # data type + if pixel_data_keyword == PixelDataKeywords.PIXEL_DATA: + if any( + any(m.is_floating_point() for m in m_list) + for m_list in real_world_value_mappings + ): + raise ValueError( + "When using an integer-valued 'pixel_array', all items " + "in 'real_world_value_mappings' must have their value " + "range specified with integers." ) - self.ImageOrientationSlide = deepcopy( - plane_orientation[0].ImageOrientationSlide + else: + if not all( + all(m.is_floating_point() for m in m_list) + for m_list in real_world_value_mappings + ): + raise ValueError( + "When using a floating point-valued pixel_array, " + "all items in 'real_world_value_mappings' must have " + "their value range specified with floats." ) - # Multi-Frame Functional Groups and Multi-Frame Dimensions - if pixel_measures is None: - if is_multiframe: - src_shared_fg = src_img.SharedFunctionalGroupsSequence[0] - pixel_measures = src_shared_fg.PixelMeasuresSequence - else: - pixel_measures = PixelMeasuresSequence( - pixel_spacing=[float(v) for v in src_img.PixelSpacing], - slice_thickness=float(src_img.SliceThickness), - spacing_between_slices=src_img.get( - 'SpacingBetweenSlices', None - ), + def add_channel_callback( + item: Dataset, + mappings: Sequence[RealWorldValueMapping], + ): + # Mappings may contain multiple mappings. Directly add the whole + # list as a sequence + item.add( + DataElement( + 0x0040_9096, # RealWorldValueMappingSequence + 'SQ', + mappings, ) + ) - sffg_item.PixelMeasuresSequence = pixel_measures - if coordinate_system == CoordinateSystemNames.PATIENT: - sffg_item.PlaneOrientationSequence = plane_orientation - - # Identity Pixel Value Transformation - transformation_item = Dataset() - transformation_item.RescaleIntercept = 0 - transformation_item.RescaleSlope = 1 - transformation_item.RescaleType = 'US' - sffg_item.PixelValueTransformationSequence = [transformation_item] + return item - # Frame VOI LUT With LUT - voi_lut_item = Dataset() - voi_lut_item.WindowCenter = format_number_as_ds(float(window_center)) - voi_lut_item.WindowWidth = format_number_as_ds(float(window_width)) - voi_lut_item.VOILUTFunction = 'LINEAR' - sffg_item.FrameVOILUTSequence = [voi_lut_item] + self._init_multiframe_image( + source_images=source_images, + pixel_array=pixel_array, + functional_groups_module=( + 'parametric-map-multi-frame-functional-groups' + ), + image_type=image_type, + photometric_interpretation=( + PhotometricInterpretationValues.MONOCHROME2 + ), + bits_allocated=bits_allocated, + samples_per_pixel=1, + use_default_pixel_value_transformation=True, + shared_voi_lut_transformations=voi_lut_transformations, + palette_color_lut_transformation=palette_color_lut_transformation, + icc_profile=icc_profile, + contains_recognizable_visual_features=( + contains_recognizable_visual_features + ), + burned_in_annotation=False, + pixel_measures=pixel_measures, + plane_orientation=plane_orientation, + plane_positions=plane_positions, + omit_empty_frames=False, + workers=workers, + dimension_organization_type=dimension_organization_type, + tile_pixel_array=tile_pixel_array, + tile_size=tile_size, + pyramid_label=pyramid_label, + pyramid_uid=pyramid_uid, + further_source_images=further_source_images, + use_extended_offset_table=use_extended_offset_table, + channel_values=real_world_value_mappings, + add_channel_callback=add_channel_callback, + pixel_data_keyword=pixel_data_keyword, + ) # Parametric Map Frame Type frame_type_item = Dataset() frame_type_item.FrameType = self.ImageType - sffg_item.ParametricMapFrameTypeSequence = [frame_type_item] - - has_multiple_mappings = pixel_array.shape[3] > 1 - if not has_multiple_mappings: - # Only if there is only a single set of mappings. Otherwise, - # the information will be stored in the Per-Frame Functional - # Groups Sequence. - sffg_item.RealWorldValueMappingSequence = \ - real_world_value_mappings[0] - - self.SharedFunctionalGroupsSequence = [sffg_item] - - # Get the correct pixel data attribute - pixel_data_type, pixel_data_attr = self._get_pixel_data_type_and_attr( - pixel_array - ) - if pixel_data_type == _PixelDataType.USHORT: - self.BitsAllocated = int(pixel_array.itemsize * 8) - self.BitsStored = self.BitsAllocated - self.HighBit = self.BitsStored - 1 - self.PixelRepresentation = 0 - elif pixel_data_type == _PixelDataType.SINGLE: - self.BitsAllocated = 32 - elif pixel_data_type == _PixelDataType.DOUBLE: - self.BitsAllocated = 64 - else: - raise ValueError('Encountered unexpected pixel data type.') + ( + self + .SharedFunctionalGroupsSequence[0] + .ParametricMapFrameTypeSequence + ) = [frame_type_item] - # Palette color lookup table + def _configure_color( + self, + palette_color_lut_transformation: PaletteColorLUTTransformation | None, + icc_profile: bytes | None, + pixel_data_keyword: PixelDataKeywords, + ) -> None: if palette_color_lut_transformation is not None: - if pixel_data_type != _PixelDataType.USHORT: + if pixel_data_keyword != PixelDataKeywords.PIXEL_DATA: raise ValueError( 'Use of palette_color_lut is only supported with integer-' 'valued pixel data.' ) - if not isinstance( - palette_color_lut_transformation, - PaletteColorLUTTransformation - ): - raise TypeError( - 'Argument "palette_color_lut_transformation" should be of ' - 'type PaletteColorLookupTable.' - ) - self.PixelPresentation = 'COLOR_RANGE' - - colors = ['Red', 'Green', 'Blue'] - for color in colors: - desc_kw = f'{color}PaletteColorLookupTableDescriptor' - data_kw = f'{color}PaletteColorLookupTableData' - desc = getattr(palette_color_lut_transformation, desc_kw) - lut = getattr(palette_color_lut_transformation, data_kw) - - setattr(self, desc_kw, desc) - setattr(self, data_kw, lut) - if hasattr( - palette_color_lut_transformation, - 'PaletteColorLookupTableUID' - ): - self.PaletteColorLookupTableUID = ( - palette_color_lut_transformation.PaletteColorLookupTableUID - ) + self.PixelPresentation = 'COLOR_RANGE' else: + if icc_profile is not None: + raise TypeError( + "Argument 'icc_profile' should " + "not be provided when 'palette_color_lut_transformation' " + "is not provided." + ) self.PixelPresentation = 'MONOCHROME' - self.copy_specimen_information(src_img) - self.copy_patient_and_study_information(src_img) - self._add_contributing_equipment(contributing_equipment, src_img) - - dimension_position_values = [ - np.unique(plane_position_values[:, index], axis=0) - for index in range(plane_position_values.shape[1]) - ] - - frames = [] - self.PerFrameFunctionalGroupsSequence = [] - for i in range(pixel_array.shape[0]): - for j in range(pixel_array.shape[3]): - pffg_item = Dataset() - - # Derivation Image - pffg_item.DerivationImageSequence = [] - - # Plane Position (Patient/Slide) - if coordinate_system == CoordinateSystemNames.SLIDE: - pffg_item.PlanePositionSlideSequence = plane_positions[i] - else: - pffg_item.PlanePositionSequence = plane_positions[i] - - # Frame Content - frame_content_item = Dataset() - frame_content_item.DimensionIndexValues = [ - int( - np.where( - dimension_position_values[idx] == pos - )[0][0] + 1 - ) - for idx, pos in enumerate(plane_position_values[i]) - ] - - pffg_item.FrameContentSequence = [frame_content_item] - - # Real World Value Mapping - if has_multiple_mappings: - # Only if there are multiple sets of mappings. Otherwise, - # the information will be stored in the Shared Functional - # Groups Sequence. - pffg_item.RealWorldValueMappingSequence = \ - real_world_value_mappings[j] - - self.PerFrameFunctionalGroupsSequence.append(pffg_item) - - plane = pixel_array[i, :, :, j] - frames.append(self._encode_frame(plane)) - - self.NumberOfFrames = len(frames) - if self.file_meta.TransferSyntaxUID.is_encapsulated: - if use_extended_offset_table: - ( - pixel_data, - self.ExtendedOffsetTable, - self.ExtendedOffsetTableLengths, - ) = encapsulate_extended(frames) - else: - pixel_data = encapsulate(frames) - else: - pixel_data = b''.join(frames) - setattr(self, pixel_data_attr, pixel_data) - - def _get_pixel_data_type_and_attr( + def _get_pixel_data_keyword( self, pixel_array: np.ndarray - ) -> tuple[_PixelDataType, str]: - """Get the data type and name of pixel data attribute. + ) -> PixelDataKeywords: + """Get the pixel data keyword to use. Parameters ---------- @@ -812,11 +635,8 @@ def _get_pixel_data_type_and_attr( Returns ------- - Tuple[highdicom.map.sop._PixelDataType, str] - A tuple where the first element is the enum value and the second - value is the DICOM pixel data attribute for the given datatype. - One of (``"PixelData"``, ``"FloatPixelData"``, - ``"DoubleFloatPixelData"``) + highdicom.enum.PixelDataKeywords: + Pixel data keyword where this pixel data should be stored. Raises ------ @@ -827,15 +647,9 @@ def _get_pixel_data_type_and_attr( """ if pixel_array.dtype.kind == 'f': if pixel_array.dtype.name == 'float32': - return ( - _PixelDataType.SINGLE, - self._pixel_data_type_map[_PixelDataType.SINGLE], - ) + return PixelDataKeywords.FLOAT_PIXEL_DATA elif pixel_array.dtype.name == 'float64': - return ( - _PixelDataType.DOUBLE, - self._pixel_data_type_map[_PixelDataType.DOUBLE], - ) + return PixelDataKeywords.DOUBLE_FLOAT_PIXEL_DATA else: raise ValueError( 'Unsupported floating-point type for pixel data: ' @@ -846,12 +660,9 @@ def _get_pixel_data_type_and_attr( if pixel_array.dtype not in (np.uint8, np.uint16): raise ValueError( 'Unsupported unsigned integer type for pixel data: ' - '16-bit unsigned integer types are supported.' + '8- and 16-bit unsigned integer types are supported.' ) - return ( - _PixelDataType.USHORT, - self._pixel_data_type_map[_PixelDataType.USHORT], - ) + return PixelDataKeywords.PIXEL_DATA raise ValueError( 'Unsupported data type for pixel data. ' 'Supported are 8-bit or 16-bit unsigned integer types as well as ' @@ -859,38 +670,66 @@ def _get_pixel_data_type_and_attr( 'floating-point types.' ) - def _encode_frame(self, pixel_array: np.ndarray) -> bytes: - """Encode a given pixel array. + @classmethod + def from_dataset( + cls, + dataset: Dataset, + copy: bool = True, + ) -> Self: + """Create instance from an existing dataset. Parameters ---------- - pixel_array: np.ndarray - The pixel array to encode + dataset: pydicom.dataset.Dataset + Dataset representing a Parametric Map. + copy: bool + If True, the underlying dataset is deep-copied such that the + original dataset remains intact. If False, this operation will + alter the original dataset in place. Returns ------- - bytes - Encoded frame - - Raises - ------ - ValueError - If the `pixel_array` is not exactly two-dimensional. + highdicom.pm.ParametricMap + Representation of the supplied dataset as a highdicom + ParametricMap. """ - if pixel_array.ndim != 2: - raise ValueError( - 'Only single frame can be encoded at at time ' - 'in case of encapsulated format encoding.' - ) - if self.file_meta.TransferSyntaxUID.is_encapsulated: - return encode_frame( - pixel_array, - transfer_syntax_uid=self.file_meta.TransferSyntaxUID, - bits_allocated=self.BitsAllocated, - bits_stored=self.BitsStored, - photometric_interpretation=self.PhotometricInterpretation, - pixel_representation=self.PixelRepresentation - ) - else: - return pixel_array.flatten().tobytes() + if dataset.SOPClassUID != '1.2.840.10008.5.1.4.1.1.30': + raise ValueError('Dataset is not a Parametric Map.') + + pm = super().from_dataset(dataset, copy=copy) + + return cast(Self, pm) + + +def pmread( + fp: str | bytes | PathLike | BinaryIO, + lazy_frame_retrieval: bool = False, +) -> ParametricMap: + """Read a parametric map image stored in DICOM File Format. + + Parameters + ---------- + fp: Union[str, bytes, os.PathLike] + Any file-like object representing a DICOM file containing a + Parametric Map image. + lazy_frame_retrieval: bool + If True, the returned parametric map will retrieve frames from the file + as requested, rather than loading in the entire object to memory + initially. This may be a good idea if file reading is slow and you are + likely to need only a subset of the frames in the parametric map. + + Returns + ------- + highdicom.pm.ParametricMap + Parametric Map image read from the file. + + """ + # This is essentially a convenience alias for the classmethod (which is + # used so that it is inherited correctly by subclasses). It is used + # because it follows the format of other similar functions around the + # library + return ParametricMap.from_file( + fp, + lazy_frame_retrieval=lazy_frame_retrieval, + ) diff --git a/src/highdicom/pr/content.py b/src/highdicom/pr/content.py index f47c6a6b..b40b84e6 100644 --- a/src/highdicom/pr/content.py +++ b/src/highdicom/pr/content.py @@ -2,10 +2,8 @@ import datetime import logging from collections import defaultdict -from io import BytesIO import numpy as np -from PIL.ImageCms import ImageCmsProfile from pydicom.dataset import Dataset from pydicom.sr.coding import Code from pydicom.multival import MultiValue @@ -14,6 +12,8 @@ from highdicom.color import CIELabColor from highdicom.content import ( + _add_content_information, + _add_palette_color_lookup_table_attributes, ContentCreatorIdentificationCodeSequence, ModalityLUTTransformation, PaletteColorLUTTransformation, @@ -36,7 +36,6 @@ from highdicom.uid import UID from highdicom.spatial import is_tiled_image from highdicom.valuerep import ( - check_person_name, _check_code_string, _check_long_string, _check_short_text @@ -801,14 +800,13 @@ def _add_presentation_state_identification_attributes( this presentation state. """ # noqa: E501 - _check_code_string(content_label) - dataset.ContentLabel = content_label - if content_description is not None: - if len(content_description) > 64: - raise ValueError( - 'Argument "content_description" must not exceed 64 characters.' - ) - dataset.ContentDescription = content_description + _add_content_information( + dataset=dataset, + content_label=content_label, + content_description=content_description, + content_creator_name=content_creator_name, + content_creator_identification=content_creator_identification, + ) now = datetime.datetime.now() dataset.PresentationCreationDate = DA(now.date()) dataset.PresentationCreationTime = TM(now.time()) @@ -829,22 +827,6 @@ def _add_presentation_state_identification_attributes( ) ] - if content_creator_name is not None: - check_person_name(content_creator_name) - dataset.ContentCreatorName = content_creator_name - - if content_creator_identification is not None: - if not isinstance( - content_creator_identification, - ContentCreatorIdentificationCodeSequence - ): - raise TypeError( - 'Argument "content_creator_identification" must be of type ' - 'ContentCreatorIdentificationCodeSequence.' - ) - dataset.ContentCreatorIdentificationCodeSequence = \ - content_creator_identification - # Not technically part of PR IODs, but we include anyway now = datetime.datetime.now() dataset.ContentDate = DA(now.date()) @@ -1587,76 +1569,6 @@ def _get_icc_profile(referenced_images: Sequence[Dataset]) -> bytes: return icc_profiles[0] -def _add_icc_profile_attributes( - dataset: Dataset, - icc_profile: bytes -) -> None: - """Add attributes of module ICC Profile. - - Parameters - ---------- - dataset: pydicom.Dataset - Dataset to which attributes should be added - icc_profile: bytes - ICC color profile to include in the presentation state. - The profile must follow the constraints listed in :dcm:`C.11.15 - `. - - """ - if icc_profile is None: - raise TypeError('Argument "icc_profile" is required.') - - cms_profile = ImageCmsProfile(BytesIO(icc_profile)) - device_class = cms_profile.profile.device_class.strip() - if device_class not in ('scnr', 'spac'): - raise ValueError( - 'The device class of the ICC Profile must be "scnr" or "spac", ' - f'got "{device_class}".' - ) - color_space = cms_profile.profile.xcolor_space.strip() - if color_space != 'RGB': - raise ValueError( - 'The color space of the ICC Profile must be "RGB", ' - f'got "{color_space}".' - ) - pcs = cms_profile.profile.connection_space.strip() - if pcs not in ('Lab', 'XYZ'): - raise ValueError( - 'The profile connection space of the ICC Profile must ' - f'be "Lab" or "XYZ", got "{pcs}".' - ) - - dataset.ICCProfile = icc_profile - - -def _add_palette_color_lookup_table_attributes( - dataset: Dataset, - palette_color_lut_transformation: PaletteColorLUTTransformation -) -> None: - """Add attributes from the Palette Color Lookup Table module. - - Parameters - ---------- - dataset: pydicom.Dataset - Dataset to which attributes should be added - palette_color_lut_transformation: highdicom.PaletteColorLUTTransformation - Description of the Palette Color LUT Transformation for transforming - grayscale into RGB color pixel values - - """ # noqa: E501 - if not isinstance( - palette_color_lut_transformation, - PaletteColorLUTTransformation - ): - raise TypeError( - 'Argument "palette_color_lut_transformation" must be of type ' - 'PaletteColorLUTTransformation.' - ) - - for element in palette_color_lut_transformation: - dataset[element.tag] = element - - def _add_softcopy_presentation_lut_attributes( dataset: Dataset, presentation_lut_transformation: PresentationLUTTransformation, diff --git a/src/highdicom/pr/sop.py b/src/highdicom/pr/sop.py index 1428e73f..90d64366 100644 --- a/src/highdicom/pr/sop.py +++ b/src/highdicom/pr/sop.py @@ -19,6 +19,8 @@ from highdicom.base import SOPClass from highdicom.base_content import ContributingEquipment from highdicom.content import ( + _add_icc_profile_attributes, + _add_palette_color_lookup_table_attributes, ContentCreatorIdentificationCodeSequence, ModalityLUTTransformation, PaletteColorLUTTransformation, @@ -28,9 +30,7 @@ from highdicom.pr.content import ( _add_displayed_area_attributes, _add_graphic_group_annotation_layer_attributes, - _add_icc_profile_attributes, _add_modality_lut_attributes, - _add_palette_color_lookup_table_attributes, _add_presentation_state_identification_attributes, _add_presentation_state_relationship_attributes, _add_softcopy_presentation_lut_attributes, diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index baf956f4..cdbe3a2a 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -3,36 +3,25 @@ from typing import cast from collections.abc import Sequence from typing_extensions import Self +import warnings -import numpy as np -from pydicom.datadict import keyword_for_tag, tag_for_keyword from pydicom.dataset import Dataset -from pydicom.sequence import Sequence as DataElementSequence from pydicom.sr.coding import Code from highdicom.color import CIELabColor from highdicom.content import ( AlgorithmIdentificationSequence, - PlanePositionSequence, ) from highdicom.enum import ( - AxisHandedness, CoordinateSystemNames, - PixelIndexDirections, ) +from highdicom.image import DimensionIndexSequence as BaseDimensionIndexSequence from highdicom.seg.enum import SegmentAlgorithmTypeValues -from highdicom.spatial import ( - _get_slice_distances, - get_normal_vector, - map_pixel_into_coordinate_system, -) from highdicom.sr.coding import CodedConcept -from highdicom.uid import UID -from highdicom.utils import compute_plane_position_slide_per_frame from highdicom._module_utils import ( check_required_attributes, - is_multiframe_image, ) +from highdicom.volume import ChannelDescriptor class SegmentDescription(Dataset): @@ -326,7 +315,7 @@ def primary_anatomic_structures(self) -> list[CodedConcept]: return list(self.PrimaryAnatomicStructureSequence) -class DimensionIndexSequence(DataElementSequence): +class DimensionIndexSequence(BaseDimensionIndexSequence): """Sequence of data elements describing dimension indices for the patient or slide coordinate system based on the Dimension Index functional @@ -336,6 +325,13 @@ class DimensionIndexSequence(DataElementSequence): ---- The order of indices is fixed. + Note + ---- + This class is deprecated and will be removed in a future version of + highdicom. User code should generally avoid this class, and if necessary, + the more general :class:`highdicom.DimensionIndexSequence` should be used + instead. + """ def __init__( @@ -354,454 +350,22 @@ def __init__( Include the segment number as a dimension index. """ - super().__init__() - if coordinate_system is None: - self._coordinate_system = None - else: - self._coordinate_system = CoordinateSystemNames(coordinate_system) - - dim_uid = UID() - - if include_segment_number: - segment_number_index = Dataset() - segment_number_index.DimensionIndexPointer = tag_for_keyword( - 'ReferencedSegmentNumber' - ) - segment_number_index.FunctionalGroupPointer = tag_for_keyword( - 'SegmentIdentificationSequence' - ) - segment_number_index.DimensionOrganizationUID = dim_uid - segment_number_index.DimensionDescriptionLabel = 'Segment Number' - self.append(segment_number_index) - - if self._coordinate_system == CoordinateSystemNames.SLIDE: - - x_axis_index = Dataset() - x_axis_index.DimensionIndexPointer = tag_for_keyword( - 'XOffsetInSlideCoordinateSystem' - ) - x_axis_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - x_axis_index.DimensionOrganizationUID = dim_uid - x_axis_index.DimensionDescriptionLabel = \ - 'X Offset in Slide Coordinate System' - - y_axis_index = Dataset() - y_axis_index.DimensionIndexPointer = tag_for_keyword( - 'YOffsetInSlideCoordinateSystem' - ) - y_axis_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - y_axis_index.DimensionOrganizationUID = dim_uid - y_axis_index.DimensionDescriptionLabel = \ - 'Y Offset in Slide Coordinate System' - - z_axis_index = Dataset() - z_axis_index.DimensionIndexPointer = tag_for_keyword( - 'ZOffsetInSlideCoordinateSystem' - ) - z_axis_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - z_axis_index.DimensionOrganizationUID = dim_uid - z_axis_index.DimensionDescriptionLabel = \ - 'Z Offset in Slide Coordinate System' - - column_dimension_index = Dataset() - column_dimension_index.DimensionIndexPointer = tag_for_keyword( - 'ColumnPositionInTotalImagePixelMatrix' - ) - column_dimension_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - column_dimension_index.DimensionOrganizationUID = dim_uid - column_dimension_index.DimensionDescriptionLabel = \ - 'Column Position In Total Image Pixel Matrix' - - row_dimension_index = Dataset() - row_dimension_index.DimensionIndexPointer = tag_for_keyword( - 'RowPositionInTotalImagePixelMatrix' - ) - row_dimension_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSlideSequence' - ) - row_dimension_index.DimensionOrganizationUID = dim_uid - row_dimension_index.DimensionDescriptionLabel = \ - 'Row Position In Total Image Pixel Matrix' - - # Organize frames for each segment similar to TILED_FULL, with - # segment position changing least frequently, followed by position - # of the row (from top to bottom) and then position of the column - # (from left to right) changing most frequently - self.extend([ - row_dimension_index, - column_dimension_index, - x_axis_index, - y_axis_index, - z_axis_index, - ]) - - elif self._coordinate_system == CoordinateSystemNames.PATIENT: - - image_position_index = Dataset() - image_position_index.DimensionIndexPointer = tag_for_keyword( - 'ImagePositionPatient' - ) - image_position_index.FunctionalGroupPointer = tag_for_keyword( - 'PlanePositionSequence' - ) - image_position_index.DimensionOrganizationUID = dim_uid - image_position_index.DimensionDescriptionLabel = \ - 'Image Position Patient' - - self.append(image_position_index) - - elif self._coordinate_system is None: - if not include_segment_number: - # Use frame label here just for the sake of using something - frame_label_index = Dataset() - frame_label_index.DimensionIndexPointer = tag_for_keyword( - 'FrameLabel' - ) - frame_label_index.FunctionalGroupPointer = tag_for_keyword( - 'FrameContentSequence' - ) - frame_label_index.DimensionOrganizationUID = dim_uid - frame_label_index.DimensionDescriptionLabel = 'Frame Label' - self.append(frame_label_index) - else: - raise ValueError( - f'Unknown coordinate system "{self._coordinate_system}"' - ) - - def get_plane_positions_of_image( - self, - image: Dataset - ) -> list[PlanePositionSequence]: - """Gets plane positions of frames in multi-frame image. - - Parameters - ---------- - image: Dataset - Multi-frame image - - Returns - ------- - List[highdicom.PlanePositionSequence] - Plane position of each frame in the image - - """ - is_multiframe = is_multiframe_image(image) - if not is_multiframe: - raise ValueError('Argument "image" must be a multi-frame image.') - - if self._coordinate_system is None: - raise ValueError( - 'Cannot calculate plane positions when images do not exist ' - 'within a frame of reference.' - ) - elif self._coordinate_system == CoordinateSystemNames.SLIDE: - if hasattr(image, 'PerFrameFunctionalGroupsSequence'): - plane_positions = [PlanePositionSequence.from_sequence( - item.PlanePositionSlideSequence - ) - for item in image.PerFrameFunctionalGroupsSequence - ] - else: - # If Dimension Organization Type is TILED_FULL, plane - # positions are implicit and need to be computed. - plane_positions = compute_plane_position_slide_per_frame(image) - else: - plane_positions = [ - PlanePositionSequence.from_sequence(item.PlanePositionSequence) - for item in image.PerFrameFunctionalGroupsSequence - ] - - return plane_positions - - def get_plane_positions_of_series( - self, - images: Sequence[Dataset] - ) -> list[PlanePositionSequence]: - """Gets plane positions for series of single-frame images. - - Parameters - ---------- - images: Sequence[Dataset] - Series of single-frame images - - Returns - ------- - List[highdicom.PlanePositionSequence] - Plane position of each frame in the image - - """ - is_multiframe = any([is_multiframe_image(img) for img in images]) - if is_multiframe: - raise ValueError( - 'Argument "images" must be a series of single-frame images.' - ) - - if self._coordinate_system is None: - raise ValueError( - 'Cannot calculate plane positions when images do not exist ' - 'within a frame of reference.' - ) - elif self._coordinate_system == CoordinateSystemNames.SLIDE: - plane_positions = [] - for img in images: - # Unfortunately, the image position is not specified relative to - # the top left corner but to the center of the image. - # Therefore, we need to compute the offset and subtract it. - center_item = img.ImageCenterPointCoordinatesSequence[0] - x_center = center_item.XOffsetInSlideCoordinateSystem - y_center = center_item.YOffsetInSlideCoordinateSystem - z_center = center_item.ZOffsetInSlideCoordinateSystem - offset_coordinate = map_pixel_into_coordinate_system( - index=((img.Columns / 2, img.Rows / 2)), - image_position=(x_center, y_center, z_center), - image_orientation=img.ImageOrientationSlide, - pixel_spacing=img.PixelSpacing - ) - center_coordinate = np.array((0., 0., 0.), dtype=float) - origin_coordinate = center_coordinate - offset_coordinate - plane_positions.append( - PlanePositionSequence( - coordinate_system=CoordinateSystemNames.SLIDE, - image_position=origin_coordinate, - pixel_matrix_position=(1, 1) - ) - ) - else: - plane_positions = [ - PlanePositionSequence( - coordinate_system=CoordinateSystemNames.PATIENT, - image_position=img.ImagePositionPatient - ) - for img in images - ] - - return plane_positions - - def get_index_position(self, pointer: str) -> int: - """Get relative position of a given dimension in the dimension index. - - Parameters - ---------- - pointer: str - Name of the dimension (keyword of the attribute), - e.g., ``"ReferencedSegmentNumber"`` - - Returns - ------- - int - Zero-based relative position - - Examples - -------- - >>> dimension_index = DimensionIndexSequence("SLIDE") - >>> i = dimension_index.get_index_position("ReferencedSegmentNumber") - >>> dimension_description = dimension_index[i] - >>> dimension_description - (0020, 9164) Dimension Organization UID ... - (0020, 9165) Dimension Index Pointer AT: (0062, 000b) - (0020, 9167) Functional Group Pointer AT: (0062, 000a) - (0020, 9421) Dimension Description Label LO: 'Segment Number' - - """ - indices = [ - i - for i, indexer in enumerate(self) - if indexer.DimensionIndexPointer == tag_for_keyword(pointer) - ] - if len(indices) == 0: - raise ValueError( - f'Dimension index does not contain a dimension "{pointer}".' - ) - return indices[0] - - def get_index_values( - self, - plane_positions: Sequence[PlanePositionSequence], - image_orientation: Sequence[float] | None = None, - index_convention: ( - str | - Sequence[PixelIndexDirections | str] - ) = ( - PixelIndexDirections.R, - PixelIndexDirections.D, - ), - handedness: AxisHandedness | str = AxisHandedness.RIGHT_HANDED, - ) -> tuple[np.ndarray, np.ndarray]: - """Get values of indexed attributes that specify position of planes. - - Parameters - ---------- - plane_positions: Sequence[highdicom.PlanePositionSequence] - Plane position of frames in a multi-frame image or in a series of - single-frame images. - image_orientation: Union[Sequence[float], None], optional - An image orientation to use to order frames within a 3D coordinate - system. By default (if ``image_orientation`` is ``None``), the - plane positions are ordered using their raw numerical values and - not along any particular spatial vector. If ``image_orientation`` - is provided, planes are ordered along the positive direction of the - vector normal to the specified. Should be a sequence of 6 floats. - This is only valid when plane position inputs contain only the - ImagePositionPatient. - index_convention: Sequence[Union[highdicom.enum.PixelIndexDirections, str]], optional - Convention used to determine how to order frames if - ``image_orientation`` is specified. Should be a sequence of two - :class:`highdicom.enum.PixelIndexDirections` or their string - representations, giving in order, the indexing conventions used for - specifying pixel indices. For example ``('R', 'D')`` means that the - first pixel index indexes the columns from left to right, and the - second pixel index indexes the rows from top to bottom (this is the - convention typically used within DICOM). As another example ``('D', - 'R')`` would switch the order of the indices to give the convention - typically used within NumPy. - - Alternatively, a single shorthand string may be passed that combines - the string representations of the two directions. So for example, - passing ``'RD'`` is equivalent to passing ``('R', 'D')``. - - This is used in combination with the ``handedness`` to determine - the positive direction used to order frames. - handedness: Union[highdicom.enum.AxisHandedness, str], optional - Choose the frame order such that the frame axis creates a - coordinate system with this handedness in the when combined with - the within-frame convention given by ``index_convention``. - - Returns - ------- - dimension_index_values: numpy.ndarray - Array of dimension index values. The first dimension corresponds - to the items in the input plane_positions sequence. The second - dimension corresponds to the dimensions of the dimension index. - The third dimension (if any) corresponds to the multiplicity - of the values, and is omitted if this is 1 for all dimensions. - plane_indices: numpy.ndarray - 1D array of planes indices for sorting frames according to their - spatial position specified by the dimension index - - Note - ---- - Includes only values of indexed attributes that specify the spatial - position of planes relative to the total pixel matrix or the frame of - reference, and excludes values of the Referenced Segment Number - attribute. - - """ # noqa: E501 - if self._coordinate_system is None: - raise RuntimeError( - 'Cannot calculate index values for multiple plane ' - 'positions when images do not exist within a frame of ' - 'reference.' - ) - - # For each dimension other than the Referenced Segment Number, - # obtain the value of the attribute that the Dimension Index Pointer - # points to in the element of the Plane Position Sequence or - # Plane Position Slide Sequence. - # Per definition, this is the Image Position Patient attribute - # in case of the patient coordinate system, or the - # X/Y/Z Offset In Slide Coordinate System and the Column/Row - # Position in Total Image Pixel Matrix attributes in case of the - # the slide coordinate system. - ref_seg_tag = tag_for_keyword("ReferencedSegmentNumber") - indexers = [ - dim_ind for dim_ind in self - if dim_ind.DimensionIndexPointer != ref_seg_tag - ] - plane_position_values = np.array([ - [ - np.array(p[0][indexer.DimensionIndexPointer].value) - for indexer in indexers - ] - for p in plane_positions - ]) - - if image_orientation is not None: - if not hasattr(plane_positions[0][0], 'ImagePositionPatient'): - raise ValueError( - 'Provided "image_orientation" is only valid when ' - 'plane_positions contain the ImagePositionPatient.' - ) - normal_vector = get_normal_vector( - image_orientation, - index_convention=index_convention, - handedness=handedness, - ) - origin_distances = _get_slice_distances( - plane_position_values[:, 0, :], - normal_vector, - ) - _, plane_sort_indices = np.unique( - origin_distances, - return_index=True, - ) - else: - # Build an array that can be used to sort planes according to the - # Dimension Index Value based on the order of the items in the - # Dimension Index Sequence. - _, plane_sort_indices = np.unique( - plane_position_values, - axis=0, - return_index=True - ) - - if len(plane_sort_indices) != len(plane_positions): - raise ValueError( - 'Input image/frame positions are not unique according to the ' - 'Dimension Index Pointers. The generated segmentation would be ' - 'ambiguous. Ensure that source images/frames have distinct ' - 'locations.' - ) - - return (plane_position_values, plane_sort_indices) - - def get_index_keywords(self) -> list[str]: - """Get keywords of attributes that specify the position of planes. - - Returns - ------- - List[str] - Keywords of indexed attributes - - Note - ---- - Includes only keywords of indexed attributes that specify the spatial - position of planes relative to the total pixel matrix or the frame of - reference, and excludes the keyword of the Referenced Segment Number - attribute. - - Examples - -------- - >>> dimension_index = DimensionIndexSequence('SLIDE') - >>> plane_positions = [ - ... PlanePositionSequence('SLIDE', [10.0, 0.0, 0.0], [1, 1]), - ... PlanePositionSequence('SLIDE', [30.0, 0.0, 0.0], [1, 2]), - ... PlanePositionSequence('SLIDE', [50.0, 0.0, 0.0], [1, 3]) - ... ] - >>> values, indices = dimension_index.get_index_values(plane_positions) - >>> names = dimension_index.get_index_keywords() - >>> for name in names: - ... print(name) - RowPositionInTotalImagePixelMatrix - ColumnPositionInTotalImagePixelMatrix - XOffsetInSlideCoordinateSystem - YOffsetInSlideCoordinateSystem - ZOffsetInSlideCoordinateSystem - >>> index = names.index("XOffsetInSlideCoordinateSystem") - >>> print(values[:, index]) - [10. 30. 50.] - - """ - referenced_segment_tag = keyword_for_tag('ReferencedSegmentNumber') - return [ - keyword_for_tag(indexer.DimensionIndexPointer) - for indexer in self - if indexer.DimensionIndexPointer != referenced_segment_tag - ] + warnings.warn( + "The highdicom.seg.DimensionIndexSequence class is deprecated and " + "will be removed in a future version of the library. User code " + "should typically avoid this class, or, if required, use the more " + "general highdicom.DimensionIndexSequence instead.", + UserWarning, + stacklevel=2, + ) + channel_dimensions = ( + [ChannelDescriptor(0x0062_000b)] # ReferencedSegmentNumber + if include_segment_number else None + ) + super().__init__( + coordinate_system=coordinate_system, + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), + channel_dimensions=channel_dimensions, + ) diff --git a/src/highdicom/seg/pyramid.py b/src/highdicom/seg/pyramid.py index f65bc799..4482875f 100644 --- a/src/highdicom/seg/pyramid.py +++ b/src/highdicom/seg/pyramid.py @@ -3,11 +3,10 @@ from collections.abc import Sequence import numpy as np -from PIL import Image from pydicom import Dataset -from pydicom.uid import VLWholeSlideMicroscopyImageStorage -from highdicom.content import PixelMeasuresSequence +from highdicom.enum import InterpolationMethods +from highdicom._pyramid import iter_derived_pyramid_levels from highdicom.seg.sop import Segmentation from highdicom.seg.enum import ( SegmentationTypeValues, @@ -23,7 +22,7 @@ def create_segmentation_pyramid( pixel_arrays: Sequence[np.ndarray], segmentation_type: str | SegmentationTypeValues, segment_descriptions: Sequence[SegmentDescription], - series_instance_uid: str, + series_instance_uid: str | None, series_number: int, manufacturer: str, manufacturer_model_name: str, @@ -88,6 +87,9 @@ def create_segmentation_pyramid( Description of each segment encoded in `pixel_array`. In the case of pixel arrays with multiple integer values, the segment description with the corresponding segment number is used to describe each segment. + series_instance_uid: str | None + UID for the output segmentation series. If ``None``, a UID will be + generated using highdicom's prefix. series_number: int Number of the output segmentation series. manufacturer: str @@ -114,9 +116,6 @@ def create_segmentation_pyramid( next is half the size, the next is a quarter of the size of the original, and the last is one eighth the size of the original. Output sizes are rounded to the nearest integer. - series_instance_uid: Optional[str], optional - UID of the output segmentation series. If not specified, UIDs are - generated automatically using highdicom's prefix. sop_instance_uids: Optional[List[str]], optional SOP instance UIDS of the output instances. If not specified, UIDs are generated automatically using highdicom's prefix. @@ -164,298 +163,34 @@ def create_segmentation_pyramid( if series_instance_uid is None: series_instance_uid = UID() - n_sources = len(source_images) - n_pix_arrays = len(pixel_arrays) - - if n_sources == 0: - raise ValueError( - 'Argument "source_images" must not be empty.' - ) - if n_pix_arrays == 0: - raise ValueError( - 'Argument "pixel_arrays" must not be empty.' - ) - - if n_sources == 1 and n_pix_arrays == 1: - if downsample_factors is None: - raise TypeError( - 'Argument "downsample_factors" must be provided when providing ' - 'only a single source image and pixel array.' - ) - if len(downsample_factors) < 1: - raise ValueError('Argument "downsample_factors" may not be empty.') - if any(f <= 1.0 for f in downsample_factors): - raise ValueError( - 'All items in "downsample_factors" must be greater than 1.' - ) - if len(downsample_factors) > 1: - if any( - z1 > z2 for z1, z2 in zip( - downsample_factors[:-1], - downsample_factors[1:] - ) - ): - raise ValueError( - 'Items in argument "downsample_factors" must be sorted in ' - 'ascending order.' - ) - n_outputs = len(downsample_factors) + 1 # including original - else: - if downsample_factors is not None: - raise TypeError( - 'Argument "downsample_factors" must not be provided when ' - 'multiple source images or pixel arrays are provided.' - ) - if n_sources > 1 and n_pix_arrays > 1: - if n_sources != n_pix_arrays: - raise ValueError( - 'If providing multiple source images and multiple pixel ' - 'arrays, the number of items in the two lists must match.' - ) - n_outputs = n_sources - else: - # Either n_sources > 1 or n_pix_arrays > 1 but not both - n_outputs = max(n_sources, n_pix_arrays) - - if sop_instance_uids is not None: - if len(sop_instance_uids) != n_outputs: - raise ValueError( - 'Number of specified SOP Instance UIDs does not match number ' - 'of output images.' - ) - - # Check the source images are appropriately ordered - for index in range(1, len(source_images)): - r0 = source_images[index - 1].TotalPixelMatrixRows - c0 = source_images[index - 1].TotalPixelMatrixColumns - r1 = source_images[index].TotalPixelMatrixRows - c1 = source_images[index].TotalPixelMatrixColumns - - if r0 <= r1 or c0 <= c1: - raise ValueError( - 'Items in argument "source_images" must be strictly ordered in ' - 'decreasing resolution.' - ) - - # Check that source images are WSI - for im in source_images: - if im.SOPClassUID != VLWholeSlideMicroscopyImageStorage: - raise ValueError( - 'Source images must have IOD VLWholeSlideMicroscopyImageStorage' - ) - - # Check that the source images are from the same series and pyramid - if len(source_images) > 1: - series_uid = source_images[0].SeriesInstanceUID - if not all( - dcm.SeriesInstanceUID == series_uid - for dcm in source_images[1:] - ): - raise ValueError( - 'All source images should belong to the same series.' - ) - if not all(hasattr(dcm, 'PyramidUID') for dcm in source_images): - raise ValueError( - 'All source images should belong to the same pyramid ' - '(share a Pyramid UID).' - ) - pyramid_uid = source_images[0].PyramidUID - if not all( - dcm.PyramidUID == pyramid_uid - for dcm in source_images[1:] - ): - raise ValueError( - 'All source images should belong to the same pyramid ' - '(share a Pyramid UID).' - ) - - # Check that pixel arrays have an appropriate shape - if len({p.ndim for p in pixel_arrays}) != 1: - raise ValueError( - 'Each item of argument "pixel_arrays" must have the same number of ' - 'dimensions.' - ) - if pixel_arrays[0].ndim == 4: - n_segment_channels = pixel_arrays[0].shape[3] - else: - n_segment_channels = None - dtype = pixel_arrays[0].dtype if dtype in (np.bool_, np.uint8, np.uint16): - resampler = Image.Resampling.NEAREST + interpolator = InterpolationMethods.NEAREST elif dtype in (np.float32, np.float64): if segmentation_type == SegmentationTypeValues.FRACTIONAL: - resampler = Image.Resampling.BILINEAR + interpolator = InterpolationMethods.LINEAR else: # This is a floating point image that will ultimately be treated as # binary - resampler = Image.Resampling.NEAREST + interpolator = InterpolationMethods.NEAREST else: raise TypeError('Pixel array has an invalid data type.') - # Checks on consistency of the pixel arrays - for pixel_array in pixel_arrays: - if pixel_array.ndim not in (2, 3, 4): - raise ValueError( - 'Each item of argument "pixel_arrays" must be a NumPy array ' - 'with 2, 3, or 4 dimensions.' - ) - if pixel_array.ndim > 2 and pixel_array.shape[0] != 1: - raise ValueError( - 'Each item of argument "pixel_arrays" must contain a single ' - 'frame, with a size of 1 along dimension 0.' - ) - if pixel_array.dtype != dtype: - raise TypeError( - 'Each item of argument "pixel_arrays" must have ' - 'the same datatype.' - ) - if pixel_array.ndim == 4: - if pixel_array.shape[3] != n_segment_channels: - raise ValueError( - 'Each item of argument "pixel_arrays" must have ' - 'the same shape down axis 3.' - ) - - # Check the pixel arrays are appropriately ordered - for index in range(1, len(pixel_arrays)): - arr0 = pixel_arrays[index - 1] - arr1 = pixel_arrays[index] - - if arr0.ndim == 2: - r0 = arr0.shape[:2] - c0 = arr0.shape[:2] - else: - r0 = arr0.shape[1:3] - c0 = arr0.shape[1:3] - - if arr1.ndim == 2: - r1 = arr1.shape[:2] - c1 = arr1.shape[:2] - else: - r1 = arr1.shape[1:3] - c1 = arr1.shape[1:3] - - if r0 <= r1 or c0 <= c1: - raise ValueError( - 'Items in argument "pixel_arrays" must be strictly ordered in ' - 'decreasing resolution.' - ) - - # Check that input dimensions match - for index, (source_image, pixel_array) in enumerate( - zip(source_images, pixel_arrays) - ): - src_shape = ( - source_image.TotalPixelMatrixRows, - source_image.TotalPixelMatrixColumns - ) - pix_shape = ( - pixel_array.shape[1:3] if pixel_array.ndim > 2 - else pixel_array.shape - ) - if pix_shape != src_shape: - raise ValueError( - "The shape of each provided pixel array must match the shape " - "of the total pixel matrix of the corresponding source image. " - f"Got pixel array of shape {pix_shape} for a source image of " - f"shape {src_shape} at index {index}." - ) - - if n_pix_arrays == 1: - # Create a pillow image for use later with resizing - if pixel_arrays[0].ndim == 2: - mask_images = [Image.fromarray(pixel_arrays[0])] - elif pixel_arrays[0].ndim == 3: - # Remove frame dimension before casting - mask_images = [Image.fromarray(pixel_arrays[0][0])] - else: # ndim = 4 - # One "Image" for each segment - mask_images = [ - Image.fromarray(pixel_arrays[0][0, :, :, i]) - for i in range(pixel_arrays[0].shape[3]) - ] - all_segs = [] - # Work "up" pyramid from high to low resolution - for output_level in range(n_outputs): - if n_sources > 1: - source_image = source_images[output_level] - else: - source_image = source_images[0] - - if n_pix_arrays > 1: - pixel_array = pixel_arrays[output_level] - else: - if output_level == 0: - pixel_array = pixel_arrays[0] - else: - if n_sources > 1: - output_size = ( - source_image.TotalPixelMatrixColumns, - source_image.TotalPixelMatrixRows - ) - else: - f = downsample_factors[output_level - 1] - output_size = ( - int(source_images[0].TotalPixelMatrixColumns / f), - int(source_images[0].TotalPixelMatrixRows / f) - ) - - # Resize each segment individually - resized_masks = [ - np.array(im.resize(output_size, resampler)) - for im in mask_images - ] - if len(resized_masks) > 1: - pixel_array = np.stack(resized_masks, axis=-1)[None] - else: - pixel_array = resized_masks[0] - - # Standardize shape of pixel array to include singleton frame dimension - if pixel_array.ndim == 2: - pixel_array = pixel_array[None] - - if n_sources == 1: - source_pixel_measures = ( - source_image - .SharedFunctionalGroupsSequence[0] - .PixelMeasuresSequence[0] - ) - src_pixel_spacing = source_pixel_measures.PixelSpacing - src_slice_thickness = source_pixel_measures.SliceThickness - - if pixel_arrays[0].ndim == 2: - # No frame dimension - orig_n_rows = pixel_arrays[0].shape[0] - orig_n_cols = pixel_arrays[0].shape[1] - else: - # Ignore 0th frame dimension - orig_n_rows = pixel_arrays[0].shape[1] - orig_n_cols = pixel_arrays[0].shape[2] - - row_spacing = ( - src_pixel_spacing[0] * - (orig_n_rows / pixel_array.shape[1]) - ) - column_spacing = ( - src_pixel_spacing[1] * - (orig_n_cols / pixel_array.shape[2]) - ) - pixel_measures = PixelMeasuresSequence( - pixel_spacing=(row_spacing, column_spacing), - slice_thickness=src_slice_thickness - ) - else: - # This will be copied from the source image - pixel_measures = None - - if sop_instance_uids is None: - sop_instance_uid = UID() - else: - sop_instance_uid = sop_instance_uids[output_level] - + for ( + source_image, + pixel_array, + pixel_measures, + instance_number, + sop_instance_uid, + ) in iter_derived_pyramid_levels( + source_images=source_images, + pixel_arrays=pixel_arrays, + interpolator=interpolator, + downsample_factors=downsample_factors, + sop_instance_uids=sop_instance_uids, + ): # Create the output segmentation seg = Segmentation( source_images=[source_image], @@ -465,7 +200,7 @@ def create_segmentation_pyramid( series_instance_uid=series_instance_uid, series_number=series_number, sop_instance_uid=sop_instance_uid, - instance_number=output_level + 1, + instance_number=instance_number, manufacturer=manufacturer, manufacturer_model_name=manufacturer_model_name, software_versions=software_versions, diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 6399856a..d7e93b25 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1,9 +1,8 @@ """Module for SOP classes of the SEG modality.""" +from concurrent.futures import Executor import logging -from concurrent.futures import Executor, Future, ProcessPoolExecutor from copy import deepcopy from os import PathLike -import pkgutil from typing import ( Any, BinaryIO, @@ -11,14 +10,11 @@ ) from collections.abc import Iterator, Sequence from typing_extensions import Self -import warnings import numpy as np from pydicom.dataelem import DataElement from pydicom.dataset import Dataset from pydicom.datadict import keyword_for_tag, tag_for_keyword -from pydicom.encaps import encapsulate, encapsulate_extended -from pydicom.pixels.utils import pack_bits from pydicom.tag import BaseTag, Tag from pydicom.uid import ( ExplicitVRLittleEndian, @@ -29,38 +25,28 @@ UID, ) from pydicom.sr.codedict import codes -from pydicom.valuerep import PersonName, format_number_as_ds +from pydicom.valuerep import PersonName from pydicom.sr.coding import Code -from highdicom._module_utils import ( - ModuleUsageValues, - get_module_usage, - is_multiframe_image, -) from highdicom.image import _Image from highdicom.base import _check_little_endian from highdicom.base_content import ContributingEquipment from highdicom.color import CIELabColor from highdicom.content import ( + _add_content_information, ContentCreatorIdentificationCodeSequence, PaletteColorLUTTransformation, PlaneOrientationSequence, PlanePositionSequence, - PixelMeasuresSequence + PixelMeasuresSequence, ) from highdicom.enum import ( CoordinateSystemNames, DimensionOrganizationTypeValues, + PhotometricInterpretationValues, + PixelRepresentationValues, ) -from highdicom.frame import encode_frame -from highdicom.pr.content import ( - _add_icc_profile_attributes, - _add_palette_color_lookup_table_attributes, -) -from highdicom.seg.content import ( - DimensionIndexSequence, - SegmentDescription, -) +from highdicom.seg.content import SegmentDescription from highdicom.seg.enum import ( SegmentationFractionalTypeValues, SegmentationTypeValues, @@ -68,17 +54,7 @@ SegmentAlgorithmTypeValues, ) from highdicom.seg.utils import iter_segments -from highdicom.spatial import ( - get_image_coordinate_system, - get_tile_array, - get_volume_positions, -) from highdicom.sr.coding import CodedConcept -from highdicom.valuerep import ( - check_person_name, - _check_code_string, - _check_long_string, -) from highdicom.volume import ( ChannelDescriptor, Volume, @@ -90,12 +66,6 @@ logger = logging.getLogger(__name__) -# This code is needed many times in loops so we precompute it -_DERIVATION_CODE = CodedConcept.from_code( - codes.cid7203.SegmentationImageDerivation -) - - def _get_unsigned_dtype(max_val: int | np.integer) -> type: """Get the smallest unsigned NumPy datatype to accommodate a value. @@ -508,29 +478,8 @@ def __init__( if len(source_images) == 0: raise ValueError('At least one source image is required.') - uniqueness_criteria = { - ( - image.StudyInstanceUID, - image.SeriesInstanceUID, - image.Rows, - image.Columns, - getattr(image, 'FrameOfReferenceUID', None), - ) - for image in source_images - } - if len(uniqueness_criteria) > 1: - raise ValueError( - 'Source images must all be part of the same series and must ' - 'have the same image dimensions (number of rows/columns).' - ) - src_img = source_images[0] - is_multiframe = is_multiframe_image(src_img) - if is_multiframe and len(source_images) > 1: - raise ValueError( - 'Only one source image should be provided in case images ' - 'are multi-frame images.' - ) + supported_transfer_syntaxes = { ImplicitVRLittleEndian, ExplicitVRLittleEndian, @@ -576,52 +525,9 @@ def __init__( **kwargs ) - # Frame of Reference - has_ref_frame_uid = hasattr(src_img, 'FrameOfReferenceUID') - if has_ref_frame_uid: - self.FrameOfReferenceUID = src_img.FrameOfReferenceUID - self.PositionReferenceIndicator = getattr( - src_img, - 'PositionReferenceIndicator', - None - ) - else: - # Only allow missing FrameOfReferenceUID if it is not required - # for this IOD - usage = get_module_usage('frame-of-reference', src_img.SOPClassUID) - if usage == ModuleUsageValues.MANDATORY: - raise ValueError( - "Source images have no Frame Of Reference UID, but it is " - "required by the IOD." - ) - - self._coordinate_system = get_image_coordinate_system(src_img) - - if self._coordinate_system is None: - # It may be possible to generalize this, but for now only a single - # source frame is permitted when there is no coordinate system - if ( - len(source_images) > 1 or - (is_multiframe and src_img.NumberOfFrames > 1) - ): - raise ValueError( - "Only a single frame is supported when the source " - "image has no Frame of Reference UID." - ) - if plane_positions is not None: - raise TypeError( - "If source images have no Frame Of Reference UID, the " - 'argument "plane_positions" may not be specified since the ' - "segmentation pixel array must be spatially aligned with " - "the source images." - ) - if plane_orientation is not None: - raise TypeError( - "If source images have no Frame Of Reference UID, the " - 'argument "plane_orientation" may not be specified since ' - "the segmentation pixel array must be spatially aligned " - "with the source images." - ) + self.copy_specimen_information(src_img) + self.copy_patient_and_study_information(src_img) + self._add_contributing_equipment(contributing_equipment, src_img) # Check segment numbers described_segment_numbers = np.array([ @@ -656,86 +562,23 @@ def __init__( "0 or 1 channel dimensions." ) - ( - pixel_array, - plane_positions, - plane_orientation, - pixel_measures, - ) = self._get_spatial_data_from_volume( - volume=pixel_array, - coordinate_system=self._coordinate_system, - frame_of_reference_uid=src_img.FrameOfReferenceUID, - plane_positions=plane_positions, - plane_orientation=plane_orientation, - pixel_measures=pixel_measures, - ) - pixel_array = cast(np.ndarray, pixel_array) - - if pixel_array.ndim == 2: - pixel_array = pixel_array[np.newaxis, ...] - if pixel_array.ndim not in [3, 4]: - raise ValueError('Pixel array must be a 2D, 3D, or 4D array.') - - is_tiled = hasattr(src_img, 'TotalPixelMatrixRows') - if tile_pixel_array and not is_tiled: - raise ValueError( - 'When argument "tile_pixel_array" is True, the source image ' - 'must be a tiled image.' - ) - if tile_pixel_array and pixel_array.shape[0] != 1: - raise ValueError( - 'When argument "tile_pixel_array" is True, the input pixel ' - 'array must contain only one "frame" representing the ' - 'entire pixel matrix.' - ) - - # Image Pixel - if tile_pixel_array: - # By default use the same tile size as the source image (even if - # they are not spatially aligned) - tile_size = tile_size or (src_img.Rows, src_img.Columns) - self.Rows, self.Columns = (tile_size) - else: - self.Rows = pixel_array.shape[1] - self.Columns = pixel_array.shape[2] - - # General Reference - self._add_source_image_references( - source_images, further_source_images, - ) - # Segmentation Image - self.ImageType = ['DERIVED', 'PRIMARY'] - self.SamplesPerPixel = 1 - self.PhotometricInterpretation = 'MONOCHROME2' - self.PixelRepresentation = 0 self.SegmentationType = segmentation_type.value - if content_label is not None: - _check_code_string(content_label) - self.ContentLabel = content_label - else: - self.ContentLabel = f'{src_img.Modality}_SEG' - self.ContentDescription = content_description - if content_creator_name is not None: - check_person_name(content_creator_name) - self.ContentCreatorName = content_creator_name - if content_creator_identification is not None: - if not isinstance( - content_creator_identification, - ContentCreatorIdentificationCodeSequence - ): - raise TypeError( - 'Argument "content_creator_identification" must be of type ' - 'ContentCreatorIdentificationCodeSequence.' - ) - self.ContentCreatorIdentificationCodeSequence = \ - content_creator_identification + _add_content_information( + dataset=self, + content_label=( + content_label if content_label is not None + else f'{src_img.Modality}_SEG' + ), + content_description=content_description, + content_creator_name=content_creator_name, + content_creator_identification=content_creator_identification, + ) if segmentation_type == SegmentationTypeValues.BINARY: dtype = np.uint8 - self.BitsAllocated = 1 - self.HighBit = 0 + bits_allocated = 1 if ( self.file_meta.TransferSyntaxUID != JPEG2000Lossless and self.file_meta.TransferSyntaxUID.is_encapsulated @@ -747,8 +590,7 @@ def __init__( ) elif segmentation_type == SegmentationTypeValues.FRACTIONAL: dtype = np.uint8 - self.BitsAllocated = 8 - self.HighBit = 7 + bits_allocated = 8 segmentation_fractional_type = SegmentationFractionalTypeValues( fractional_type ) @@ -766,26 +608,10 @@ def __init__( raise ValueError( "Too many segments to represent with a 16 bit integer." ) - self.BitsAllocated = np.iinfo(dtype).bits - self.HighBit = self.BitsAllocated - 1 - self.BitsStored = self.BitsAllocated + bits_allocated = np.iinfo(dtype).bits self.PixelPaddingValue = 0 - self.BitsStored = self.BitsAllocated - self.LossyImageCompression = getattr( - src_img, - 'LossyImageCompression', - '00' - ) - if self.LossyImageCompression == '01': - if 'LossyImageCompressionRatio' in src_img: - self.LossyImageCompressionRatio = \ - src_img.LossyImageCompressionRatio - if 'LossyImageCompressionMethod' in src_img: - self.LossyImageCompressionMethod = \ - src_img.LossyImageCompressionMethod - - self._configure_color( + photometric_interpretation = self._configure_color( segmentation_type=segmentation_type, palette_color_lut_transformation=palette_color_lut_transformation, icc_profile=icc_profile, @@ -793,50 +619,47 @@ def __init__( max_described_segment=int(described_segment_numbers.max()), ) - # Multi-Resolution Pyramid - if pyramid_uid is not None: - if not is_tiled: - raise TypeError( - 'Argument "pyramid_uid" should only be specified ' - 'for tiled images.' - ) - if ( - self._coordinate_system is None or - self._coordinate_system != CoordinateSystemNames.SLIDE - ): - raise TypeError( - 'Argument "pyramid_uid" should only be specified for ' - 'segmentations in the SLIDE coordinate system.' - ) - self.PyramidUID = pyramid_uid - - if pyramid_label is not None: - _check_long_string(pyramid_label) - self.PyramidLabel = pyramid_label - - elif pyramid_label is not None: - raise TypeError( - 'Argument "pyramid_label" should not be specified if ' - '"pyramid_uid" is not specified.' - ) + self._add_segment_descriptions( + segment_descriptions, + segmentation_type, + ) # Multi-Frame Functional Groups and Multi-Frame Dimensions include_segment_number = ( segmentation_type != SegmentationTypeValues.LABELMAP ) - self.DimensionIndexSequence = DimensionIndexSequence( - coordinate_system=self._coordinate_system, - include_segment_number=include_segment_number, - ) - dimension_organization = Dataset() - dimension_organization.DimensionOrganizationUID = \ - self.DimensionIndexSequence[0].DimensionOrganizationUID - self.DimensionOrganizationSequence = [dimension_organization] + if include_segment_number: + # Include ReferencedSegmentNumber as an indexed channel + channel_dimension_index = ChannelDescriptor(0x0062_000b) + channel_values = described_segment_numbers.tolist() + + # Function to use as callback to add the segment identification + def _add_segment_identification( + pffg_item: Dataset, + segment_number: int, + ): + identification = Dataset() + identification.add( + DataElement( + 0x0062000b, # ReferencedSegmentNumber + 'US', + int(segment_number) + ) + ) + pffg_item.add( + DataElement( + 0x0062000a, # SegmentIdentificationSequence + 'SQ', + [identification] + ) + ) + return pffg_item - self._add_segment_descriptions( - segment_descriptions, - segmentation_type, - ) + add_channel_callback = _add_segment_identification + else: + channel_dimension_index = None + channel_values = None + add_channel_callback = None # Checks on pixels and overlap pixel_array, segments_overlap = self._check_and_cast_pixel_array( @@ -847,461 +670,52 @@ def __init__( ) self.SegmentsOverlap = segments_overlap.value - ( - plane_positions, - plane_orientation, - pixel_measures, - plane_position_values, - plane_sort_index, - derivation_source_image_items, - ) = self._prepare_spatial_metadata( - plane_positions=plane_positions, - plane_orientation=plane_orientation, - pixel_measures=pixel_measures, - source_images=source_images, - further_source_images=further_source_images or [], - tile_pixel_array=tile_pixel_array, - tile_size=tile_size, - frame_shape=pixel_array.shape[1:3], - number_of_planes=pixel_array.shape[0], - dimension_organization_type=dimension_organization_type, - ) - - # Shared functional groops - sffg_item = Dataset() - if ( - self._coordinate_system is not None and - self._coordinate_system == CoordinateSystemNames.PATIENT + def preprocess_channel_callback( + plane_array: np.ndarray, + segment_number: int, + segment_index: int, # unused ): - sffg_item.PlaneOrientationSequence = plane_orientation - - # Automatically populate the spacing between slices in the - # pixel measures if it was not provided. This is done on the - # initial plane positions, before any removals, to give the - # receiver more information about how to reconstruct a volume - # from the frames in the case that slices are omitted - if 'SpacingBetweenSlices' not in pixel_measures[0]: - ori = plane_orientation[0].ImageOrientationPatient - slice_spacing, _ = get_volume_positions( - image_positions=plane_position_values[:, 0, :], - image_orientation=ori, - ) - if slice_spacing is not None: - pixel_measures[0].SpacingBetweenSlices = ( - format_number_as_ds(slice_spacing) - ) - - if pixel_measures is not None: - sffg_item.PixelMeasuresSequence = pixel_measures - self.SharedFunctionalGroupsSequence = [sffg_item] - - # Find indices such that empty planes are removed - if omit_empty_frames: - if tile_pixel_array: - included_plane_indices, is_empty = \ - self._get_nonempty_tile_indices( - pixel_array, - plane_positions=plane_positions, - rows=self.Rows, - columns=self.Columns, - ) - else: - included_plane_indices, is_empty = \ - self._get_nonempty_plane_indices(pixel_array) - if is_empty: - # Cannot omit empty frames when all frames are empty - omit_empty_frames = False - included_plane_indices = range(len(plane_positions)) - else: - # Remove all empty plane positions from the list of sorted - # plane position indices - included_plane_indices_set = set(included_plane_indices) - plane_sort_index = [ - ind for ind in plane_sort_index - if ind in included_plane_indices_set - ] - else: - included_plane_indices = range(len(plane_positions)) + return self._get_segment_pixel_array( + plane_array, + segment_number=segment_number, + described_segment_numbers=described_segment_numbers, + segmentation_type=segmentation_type, + max_fractional_value=max_fractional_value, + dtype=dtype, + ) - # Dimension Organization Type - dimension_organization_type = self._check_tiled_dimension_organization( - dimension_organization_type=dimension_organization_type, - is_tiled=is_tiled, - omit_empty_frames=omit_empty_frames, + self._init_multiframe_image( + source_images=source_images, + pixel_array=pixel_array, + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), + photometric_interpretation=photometric_interpretation, + bits_allocated=bits_allocated, + samples_per_pixel=1, + image_type=['DERIVED', 'PRIMARY'], + pixel_representation=PixelRepresentationValues.UNSIGNED_INTEGER, + use_default_pixel_value_transformation=False, + palette_color_lut_transformation=palette_color_lut_transformation, + icc_profile=icc_profile, + pixel_measures=pixel_measures, + plane_orientation=plane_orientation, plane_positions=plane_positions, + omit_empty_frames=omit_empty_frames, + workers=workers, + dimension_organization_type=dimension_organization_type, tile_pixel_array=tile_pixel_array, - rows=self.Rows, - columns=self.Columns, - ) - - if ( - self._coordinate_system is not None and - dimension_organization_type != - DimensionOrganizationTypeValues.TILED_FULL - ): - # Get unique values of attributes in the Plane Position Sequence or - # Plane Position Slide Sequence, which define the position of the - # plane with respect to the three dimensional patient or slide - # coordinate system, respectively. These can subsequently be used - # to look up the relative position of a plane relative to the - # indexed dimension. - unique_dimension_values = [ - np.unique( - plane_position_values[included_plane_indices, index], - axis=0 - ) - for index in range(plane_position_values.shape[1]) - ] - else: - unique_dimension_values = [None] - - if self._coordinate_system == CoordinateSystemNames.PATIENT: - inferred_dim_org_type = None - - # To be considered "3D", a segmentation should have frames that are - # differentiated only by location. This rules out any non-labelmap - # segmentations with more than a single segment. - # Further, only segmentation with multiple spatial positions in the - # final segmentation should be considered to have 3D dimension - # organization type - if ( - len(included_plane_indices) > 1 and - ( - segmentation_type == SegmentationTypeValues.LABELMAP or - len(described_segment_numbers) == 1 - ) - ): - # Calculate the spacing using only the included planes, and - # enforce ordering - ori = plane_orientation[0].ImageOrientationPatient - spacing, _ = get_volume_positions( - image_positions=plane_position_values[ - included_plane_indices, 0, : - ], - image_orientation=ori, - sort=False, - ) - if spacing is not None and spacing > 0.0: - inferred_dim_org_type = ( - DimensionOrganizationTypeValues.THREE_DIMENSIONAL - ) - - if ( - dimension_organization_type == - DimensionOrganizationTypeValues.THREE_DIMENSIONAL - ) and inferred_dim_org_type is None: - raise ValueError( - 'Dimension organization "3D" has been specified, ' - 'but the source image is not a regularly-spaced 3D ' - 'volume.' - ) - dimension_organization_type = inferred_dim_org_type - - if dimension_organization_type is not None: - self.DimensionOrganizationType = dimension_organization_type.value - - if ( - self._coordinate_system is not None and - self._coordinate_system == CoordinateSystemNames.SLIDE - ): - plane_position_names = ( - self.DimensionIndexSequence.get_index_keywords() - ) - row_dim_index = plane_position_names.index( - 'RowPositionInTotalImagePixelMatrix' - ) - col_dim_index = plane_position_names.index( - 'ColumnPositionInTotalImagePixelMatrix' - ) - - is_encaps = self.file_meta.TransferSyntaxUID.is_encapsulated - process_pool: Executor | None = None - - if not isinstance(workers, (int, Executor)): - raise TypeError( - 'Argument "workers" must be of type int or ' - 'concurrent.futures.Executor (or a derived class).' - ) - using_multiprocessing = ( - isinstance(workers, Executor) or workers != 0 - ) - - # List of frames. In the case of native transfer syntaxes, we will - # collect a list of frames as flattened NumPy arrays for bitpacking at - # the end. In the case of encapsulated transfer syntaxes with no - # workers, we will accumulate a list of encoded frames to encapsulate - # at the end - frames: list[bytes] | list[np.ndarray] = [] - - # In the case of native encoding when the number pixels in a frame is - # not a multiple of 8. This array carries "leftover" pixels that - # couldn't be encoded in previous iterations, to future iterations. This - # saves having to keep the entire un-endoded array in memory, which can - # get extremely heavy on memory in the case of very large arrays - remainder_pixels = np.empty((0, ), dtype=np.uint8) - - if is_encaps: - if using_multiprocessing: - # In the case of encapsulated transfer syntaxes with multiple - # workers, we will accumulate a list of encoded frames to - # encapsulate at the end - frame_futures: list[Future] = [] - - # Use the existing executor or create one - if isinstance(workers, Executor): - process_pool = workers - else: - # If workers is negative, pass None to use all processors - process_pool = ProcessPoolExecutor( - workers if workers > 0 else None - ) - - # Parameters to use when calling the encode_frame function in - # either of the above two cases - encode_frame_kwargs = dict( - transfer_syntax_uid=self.file_meta.TransferSyntaxUID, - bits_allocated=self.BitsAllocated, - bits_stored=self.BitsStored, - photometric_interpretation=self.PhotometricInterpretation, - pixel_representation=self.PixelRepresentation - ) - else: - if using_multiprocessing: - warnings.warn( - "Setting workers != 0 or passing an instance of " - "concurrent.futures.Executor when using a non-encapsulated " - "transfer syntax has no effect.", - UserWarning, - stacklevel=2, - ) - using_multiprocessing = False - - # Information about individual frames is placed into the - # PerFrameFunctionalGroupsSequence. Note that a *very* significant - # efficiency gain is observed when building this as a Python list - # rather than a pydicom sequence, and then converting to a pydicom - # sequence at the end - pffg_sequence: list[Dataset] = [] - - # We want the larger loop to work in the labelmap cases (where segments - # are dealt with together) and the other cases (where segments are - # dealt with separately). So we define a suitable iterable here for - # each case - segments_iterable = ( - [None] if segmentation_type == SegmentationTypeValues.LABELMAP - else described_segment_numbers + tile_size=tile_size, + pyramid_label=pyramid_label, + pyramid_uid=pyramid_uid, + further_source_images=further_source_images, + use_extended_offset_table=use_extended_offset_table, + channel_dimension_index=channel_dimension_index, + channel_values=channel_values, + add_channel_callback=add_channel_callback, + preprocess_channel_callback=preprocess_channel_callback, ) - for segment_number in segments_iterable: - - for plane_dim_ind, plane_index in enumerate(plane_sort_index, 1): - - if tile_pixel_array: - if ( - dimension_organization_type == - DimensionOrganizationTypeValues.TILED_FULL - ): - row_offset = int( - plane_position_values[plane_index, row_dim_index] - ) - column_offset = int( - plane_position_values[plane_index, col_dim_index] - ) - else: - pos = plane_positions[plane_index][0] - row_offset = pos.RowPositionInTotalImagePixelMatrix - column_offset = ( - pos.ColumnPositionInTotalImagePixelMatrix - ) - - plane_array = get_tile_array( - pixel_array[0], - row_offset=row_offset, - column_offset=column_offset, - tile_rows=self.Rows, - tile_columns=self.Columns, - ) - else: - # Select the relevant existing frame - plane_array = pixel_array[plane_index] - - if segment_number is None: - # Deal with all segments at once - segment_array = plane_array - else: - # Pixel array for just this segment and this position - segment_array = self._get_segment_pixel_array( - plane_array, - segment_number=segment_number, - described_segment_numbers=described_segment_numbers, - segmentation_type=segmentation_type, - max_fractional_value=max_fractional_value, - dtype=dtype, - ) - - # Even though completely empty planes were removed earlier, - # there may still be planes in which this specific segment is - # absent. Such frames should be removed - if segment_number is not None: - if omit_empty_frames and not np.any(segment_array): - logger.debug( - f'skip empty plane {plane_index} of segment ' - f'#{segment_number}' - ) - continue - - # Log a debug message - if segment_number is None: - msg = f'add plane #{plane_index}' - else: - msg = ( - f'add plane #{plane_index} for segment ' - f'#{segment_number}' - ) - logger.debug(msg) - - if ( - dimension_organization_type != - DimensionOrganizationTypeValues.TILED_FULL - ): - # No per-frame functional group for TILED FULL - - # Get the item of the PerFrameFunctionalGroupsSequence for - # this segmentation frame - if self._coordinate_system is not None: - plane_pos_val = plane_position_values[plane_index] - if ( - self._coordinate_system == - CoordinateSystemNames.SLIDE - ): - try: - dimension_index_values = [ - int( - np.where( - unique_dimension_values[idx] == pos - )[0][0] + 1 - ) - for idx, pos in enumerate(plane_pos_val) - ] - except IndexError as error: - raise IndexError( - 'Could not determine position of plane ' - f'#{plane_index} in three dimensional ' - 'coordinate system based on dimension ' - f'index values: {error}' - ) from error - else: - dimension_index_values = [plane_dim_ind] - else: - if segmentation_type == SegmentationTypeValues.LABELMAP: - # Here we have to use the "Frame Label" dimension - # value (which is used just to have one index since - # Referenced Segment cannot be used) - dimension_index_values = [1] - else: - dimension_index_values = [] - - plane_derivation_items = ( - derivation_source_image_items[plane_index] - if derivation_source_image_items is not None else None - ) - pffg_item = self._get_pffg_item( - segment_number=segment_number, - dimension_index_values=dimension_index_values, - plane_position=plane_positions[plane_index], - derivation_source_image_items=plane_derivation_items, - coordinate_system=self._coordinate_system, - ) - pffg_sequence.append(pffg_item) - - # Add the segmentation pixel array for this frame to the list - if is_encaps: - if process_pool is None: - # Encode this frame and add resulting bytes to the list - # for encapsulation at the end - frames.append( - encode_frame( - segment_array, - **encode_frame_kwargs, - ) - ) - else: - # Submit this frame for encoding this frame and add the - # future to the list for encapsulation at the end - future = process_pool.submit( - encode_frame, - array=segment_array, - **encode_frame_kwargs, - ) - frame_futures.append(future) - else: - flat_array = segment_array.flatten() - if ( - self.SegmentationType == - SegmentationTypeValues.BINARY.value and - (self.Rows * self.Columns) // 8 != 0 - ): - # Need to encode a multiple of 8 pixels at a time - full_array = np.concatenate( - [remainder_pixels, flat_array] - ) - # Round down to closest multiple of 8 - n_pixels_to_take = 8 * (len(full_array) // 8) - to_encode = full_array[:n_pixels_to_take] - remainder_pixels = full_array[n_pixels_to_take:] - else: - # Simple - each frame can be individually encoded - to_encode = flat_array - - frames.append(self._encode_pixels_native(to_encode)) - - if ( - dimension_organization_type != - DimensionOrganizationTypeValues.TILED_FULL - ): - self.PerFrameFunctionalGroupsSequence = pffg_sequence - - if is_encaps: - if process_pool is not None: - frames = [ - fut.result() for fut in frame_futures - ] - - # Shutdown the pool if we created it, otherwise it is the - # caller's responsibility - if process_pool is not workers: - process_pool.shutdown() - - # Encapsulate all pre-compressed frames - self.NumberOfFrames = len(frames) - if use_extended_offset_table: - ( - self.PixelData, - self.ExtendedOffsetTable, - self.ExtendedOffsetTableLengths, - ) = encapsulate_extended(frames) - else: - self.PixelData = encapsulate(frames) - else: - self.NumberOfFrames = len(frames) - - # May need to add in a final set of pixels - if len(remainder_pixels) > 0: - frames.append(self._encode_pixels_native(remainder_pixels)) - - self.PixelData = b''.join(frames) - - # Add a null trailing byte if required - if len(self.PixelData) % 2 == 1: - self.PixelData += b'0' - - self.copy_specimen_information(src_img) - self.copy_patient_and_study_information(src_img) - self._add_contributing_equipment(contributing_equipment, src_img) - - # Build lookup tables for efficient decoding - self._build_luts() - def add_segments( self, pixel_array: np.ndarray, @@ -1458,7 +872,7 @@ def _configure_color( icc_profile: bytes | None, segment_descriptions: Sequence[SegmentDescription], max_described_segment: int, - ) -> None: + ) -> PhotometricInterpretationValues: # Use PALETTE COLOR photometric interpretation in the case # of a labelmap segmentation with a provided LUT, MONOCHROME2 # otherwise @@ -1472,10 +886,10 @@ def _configure_color( "'palette_color_lut_transformation' " "is not specified." ) + return PhotometricInterpretationValues.MONOCHROME2 else: # Using photometric interpretation "PALETTE COLOR" # need to specify the LUT in this case - self.PhotometricInterpretation = 'PALETTE COLOR' # Checks on the validity of the LUT if not isinstance( @@ -1514,25 +928,9 @@ def _configure_color( 'color when using a palette color LUT.' ) - # Add the LUT to this instance - _add_palette_color_lookup_table_attributes( - self, - palette_color_lut_transformation, - ) - - if icc_profile is None: - # Use default sRGB profile - icc_profile = pkgutil.get_data( - 'highdicom', - '_icc_profiles/sRGB_v4_ICC_preference.icc' - ) - _add_icc_profile_attributes( - self, - icc_profile=icc_profile - ) + return PhotometricInterpretationValues.PALETTE_COLOR else: - self.PhotometricInterpretation = 'MONOCHROME2' if palette_color_lut_transformation is not None: raise TypeError( "Argument 'palette_color_lut_transformation' should " @@ -1545,22 +943,23 @@ def _configure_color( "not be provided when 'segmentation_type' is " f"'{segmentation_type.value}'." ) + return PhotometricInterpretationValues.MONOCHROME2 @classmethod def _check_and_cast_pixel_array( cls, - pixel_array: np.ndarray, + pixel_array: np.ndarray | Volume, segment_numbers: np.ndarray, segmentation_type: SegmentationTypeValues, dtype: type, - ) -> tuple[np.ndarray, SegmentsOverlapValues]: + ) -> tuple[np.ndarray | Volume, SegmentsOverlapValues]: """Checks on the shape and data type of the pixel array. Also checks for overlapping segments and returns the result. Parameters ---------- - pixel_array: numpy.ndarray + pixel_array: numpy.ndarray | highdicom.Volume The segmentation pixel array. segment_numbers: numpy.ndarray The segment numbers from the segment descriptions, in the order @@ -1572,7 +971,7 @@ def _check_and_cast_pixel_array( Returns ------- - pixel_array: numpyp.ndarray + pixel_array: numpy.ndarray | highdicom.Volume Input pixel array with the data type simplified if possible. segments_overlap: highdicom.seg.SegmentationOverlaps The value for the SegmentationOverlaps attribute, inferred from the @@ -1587,6 +986,11 @@ def _check_and_cast_pixel_array( # results are reused wherever possible number_of_segments = len(segment_numbers) + plain_array = ( + pixel_array.array if isinstance(pixel_array, Volume) + else pixel_array + ) + if pixel_array.ndim == 4: # Check that the number of segments in the array matches if pixel_array.shape[-1] != number_of_segments: @@ -1599,7 +1003,7 @@ def _check_and_cast_pixel_array( if pixel_array.dtype in (np.bool_, np.uint8, np.uint16): - if pixel_array.ndim == 3: + if pixel_array.ndim in (2, 3): # A label-map style array where pixel values represent # segment associations @@ -1617,7 +1021,7 @@ def _check_and_cast_pixel_array( # to check the max pixel value, which is MUCH more # efficient than calculating the set of unique values has_undescribed_segments = ( - pixel_array.max() > number_of_segments + plain_array.max() > number_of_segments ) else: # The general case, much slower @@ -1625,7 +1029,7 @@ def _check_and_cast_pixel_array( [np.array([0]), segment_numbers] ) has_undescribed_segments = len( - np.setdiff1d(pixel_array, numbers_with_bg) + np.setdiff1d(plain_array, numbers_with_bg) ) != 0 if has_undescribed_segments: @@ -1638,7 +1042,7 @@ def _check_and_cast_pixel_array( # cannot overlap segments_overlap = SegmentsOverlapValues.NO else: - max_pixel = pixel_array.max() + max_pixel = plain_array.max() # Pixel array is 4D where each segment is stacked down # the last dimension @@ -1658,14 +1062,14 @@ def _check_and_cast_pixel_array( # A single segment does not overlap segments_overlap = SegmentsOverlapValues.NO else: - sum_over_segments = pixel_array.sum(axis=-1) + sum_over_segments = plain_array.sum(axis=-1) if np.any(sum_over_segments > 1): segments_overlap = SegmentsOverlapValues.YES else: segments_overlap = SegmentsOverlapValues.NO elif pixel_array.dtype in (np.float32, np.float64): - unique_values = np.unique(pixel_array) + unique_values = np.unique(plain_array) if np.min(unique_values) < 0.0 or np.max(unique_values) > 1.0: raise ValueError( 'Floating point pixel array values must be in the ' @@ -1691,15 +1095,18 @@ def _check_and_cast_pixel_array( if len(unique_values) == 1 and unique_values[0] == 0.0: # All pixels are zero: there can be no overlap segments_overlap = SegmentsOverlapValues.NO - elif pixel_array.ndim == 3 or pixel_array.shape[-1] == 1: + elif pixel_array.ndim in (2, 3) or pixel_array.shape[-1] == 1: # A single segment does not overlap segments_overlap = SegmentsOverlapValues.NO - elif pixel_array.sum(axis=-1).max() > 1: + elif plain_array.sum(axis=-1).max() > 1: segments_overlap = SegmentsOverlapValues.YES else: segments_overlap = SegmentsOverlapValues.NO else: - if (pixel_array.ndim == 3) or (pixel_array.shape[-1] == 1): + if ( + (pixel_array.ndim in (2, 3)) or + (pixel_array.shape[-1] == 1) + ): # A single segment does not overlap segments_overlap = SegmentsOverlapValues.NO else: @@ -1718,58 +1125,20 @@ def _check_and_cast_pixel_array( ) if pixel_array.ndim == 4: - pixel_array = cls._combine_segments( - pixel_array, - labelmap_dtype=dtype - ) + if isinstance(pixel_array, Volume): + pixel_array = pixel_array.with_array( + cls._combine_segments(plain_array, dtype) + ) + else: + pixel_array = cls._combine_segments( + pixel_array, + labelmap_dtype=dtype + ) else: pixel_array = pixel_array.astype(dtype) return pixel_array, segments_overlap - @staticmethod - def _get_nonempty_plane_indices( - pixel_array: np.ndarray - ) -> tuple[list[int], bool]: - """Get a list of all indices of original planes that are non-empty. - - Empty planes (without any positive pixels in any of the segments) do - not need to be included in the segmentation image. This method finds a - list of indices of the input frames that are non-empty, and therefore - should be included in the segmentation image. - - Parameters - ---------- - pixel_array: numpy.ndarray - Segmentation pixel array - - Returns - ------- - included_plane_indices : List[int] - List giving for each plane position in the resulting segmentation - image the index of the corresponding frame in the original pixel - array. - is_empty: bool - Whether the entire image is empty. If so, empty frames should not - be omitted. - - """ - # This list tracks which source image each non-empty frame came from - source_image_indices = [ - i for i, frm in enumerate(pixel_array) - if np.any(frm) - ] - - if len(source_image_indices) == 0: - logger.warning( - 'Encoding an empty segmentation with "omit_empty_frames" ' - 'set to True. Reverting to encoding all frames since omitting ' - 'all frames is not possible.' - ) - return (list(range(pixel_array.shape[0])), True) - - return (source_image_indices, False) - @staticmethod def _combine_segments( pixel_array: np.ndarray, @@ -1808,68 +1177,6 @@ def _combine_segments( return pixel_array - @staticmethod - def _get_nonempty_tile_indices( - pixel_array: np.ndarray, - plane_positions: Sequence[PlanePositionSequence], - rows: int, - columns: int, - ) -> tuple[list[int], bool]: - """Get a list of all indices of tile locations that are non-empty. - - This is similar to _get_nonempty_plane_indices, but works on a total - pixel matrix rather than a set of frames. Empty planes (without any - positive pixels in any of the segments) do not need to be included in - the segmentation image. This method finds a list of indices of the - input frames that are non-empty, and therefore should be included in - the segmentation image. - - Parameters - ---------- - pixel_array: numpy.ndarray - Segmentation pixel array - plane_positions: Sequence[highdicom.PlanePositionSequence] - Plane positions of each tile. - rows: int - Number of rows in each tile. - columns: int - Number of columns in each tile. - - Returns - ------- - included_plane_indices : List[int] - List giving for each plane position in the resulting segmentation - image the index of the corresponding frame in the original pixel - array. - is_empty: bool - Whether the entire image is empty. If so, empty frames should not - be omitted. - - """ - # This list tracks which source image each non-empty frame came from - source_image_indices = [ - i for i, pos in enumerate(plane_positions) - if np.any( - get_tile_array( - pixel_array[0], - row_offset=pos[0].RowPositionInTotalImagePixelMatrix, - column_offset=pos[0].ColumnPositionInTotalImagePixelMatrix, - tile_rows=rows, - tile_columns=columns, - ) - ) - ] - - if len(source_image_indices) == 0: - logger.warning( - 'Encoding an empty segmentation with "omit_empty_frames" ' - 'set to True. Reverting to encoding all frames since omitting ' - 'all frames is not possible.' - ) - return (list(range(len(plane_positions))), True) - - return (source_image_indices, False) - @staticmethod def _get_segment_pixel_array( pixel_array: np.ndarray, @@ -1957,178 +1264,6 @@ def _get_segment_pixel_array( return segment_array - @staticmethod - def _get_pffg_item( - segment_number: int | None, - dimension_index_values: list[int], - plane_position: PlanePositionSequence, - derivation_source_image_items: Sequence[Dataset] | None, - coordinate_system: CoordinateSystemNames | None, - ) -> Dataset: - """Get a single item of the Per Frame Functional Groups Sequence. - - This is a helper method used in the constructor. - - Parameters - ---------- - segment_number: Optional[int] - Segment number of this segmentation frame. If None, this is a - LABELMAP segmentation in which each frame has no segment number. - dimension_index_values: List[int] - Dimension index values (except segment number) for this frame. - plane_position: highdicom.seg.PlanePositionSequence - Plane position of this frame. - derivation_source_image_items: Sequence[pydicom.Dataset] | None - Items of the Source Image Sequence, to place into the Derivation - Image Sequence, if any, for this frame. - coordinate_system: Optional[highdicom.CoordinateSystemNames] - Coordinate system used, if any. - - Returns - ------- - pydicom.Dataset - Dataset representing the item of the Per Frame Functional Groups - Sequence for this segmentation frame. - - """ - # NB this function is called many times in a loop when there are a - # large number of frames, and has been observed to dominate the - # creation time of some segmentations. Therefore we use low-level - # pydicom primitives to improve performance as much as possible - pffg_item = Dataset() - frame_content_item = Dataset() - - if segment_number is None: - all_index_values = dimension_index_values - else: - all_index_values = [int(segment_number)] + dimension_index_values - - frame_content_item.add( - DataElement( - 0x00209157, # DimensionIndexValues - 'UL', - all_index_values, - ) - ) - - if segment_number is None and coordinate_system is None: - # If this is an labelmap segmentation of an image that has no frame - # of reference, we need to create a dummy frame label to be pointed - # to as a dimension index because there is nothing else appropriate - # to use - frame_content_item.add( - DataElement( - 0x00209453, # FrameLabel - 'LO', - "Segmentation Frame", - ) - ) - - pffg_item.add( - DataElement( - 0x00209111, # FrameContentSequence - 'SQ', - [frame_content_item] - ) - ) - - if coordinate_system is not None: - if coordinate_system == CoordinateSystemNames.SLIDE: - pffg_item.add( - DataElement( - 0x0048021a, # PlanePositionSlideSequence - 'SQ', - plane_position - ) - ) - else: - pffg_item.add( - DataElement( - 0x00209113, # PlanePositionSequence - 'SQ', - plane_position - ) - ) - - if ( - derivation_source_image_items is not None and - len(derivation_source_image_items) > 0 - ): - derivation_image_item = Dataset() - derivation_image_item.add( - DataElement( - 0x00089215, # DerivationCodeSequence - 'SQ', - [_DERIVATION_CODE] - ) - ) - - derivation_image_item.add( - DataElement( - 0x00082112, # SourceImageSequence - 'SQ', - derivation_source_image_items, - ) - ) - pffg_item.add( - DataElement( - 0x00089124, # DerivationImageSequence - 'SQ', - [derivation_image_item] - ) - ) - else: - # Determining the source images that map to the frame is not - # always trivial. Since DerivationImageSequence is a type 2 - # attribute, we leave its value empty. - pffg_item.add( - DataElement( - 0x00089124, # DerivationImageSequence - 'SQ', - [] - ) - ) - - if segment_number is not None: - identification = Dataset() - identification.add( - DataElement( - 0x0062000b, # ReferencedSegmentNumber - 'US', - int(segment_number) - ) - ) - pffg_item.add( - DataElement( - 0x0062000a, # SegmentIdentificationSequence - 'SQ', - [identification] - ) - ) - - return pffg_item - - def _encode_pixels_native(self, planes: np.ndarray) -> bytes: - """Encode pixel planes using a native transfer syntax. - - Parameters - ---------- - planes: numpy.ndarray - Array representing one or more segmentation image planes. If - multiple image planes, planes stacked down the first dimension - (index 0). - - Returns - ------- - bytes - Encoded pixels - - """ - if self.SegmentationType == SegmentationTypeValues.BINARY.value: - return pack_bits(planes, pad=False) - else: - return planes.tobytes() - @classmethod def from_dataset( cls, diff --git a/src/highdicom/spatial.py b/src/highdicom/spatial.py index 4685e97b..325935ab 100644 --- a/src/highdicom/spatial.py +++ b/src/highdicom/spatial.py @@ -304,8 +304,8 @@ def iter_tiled_full_frame_data( Parameters ---------- dataset: pydicom.dataset.Dataset - VL Whole Slide Microscopy Image or Segmentation Image using the - "TILED_FULL" DimensionOrganizationType. + VL Whole Slide Microscopy Image, Segmentation Image, or Parametric Map + using the "TILED_FULL" DimensionOrganizationType. Returns ------- @@ -334,15 +334,33 @@ def iter_tiled_full_frame_data( units. """ - allowed_sop_class_uids = { - '1.2.840.10008.5.1.4.1.1.77.1.6', # VL Whole Slide Microscopy Image - '1.2.840.10008.5.1.4.1.1.66.4', # Segmentation Image - '1.2.840.10008.5.1.4.1.1.66.7', # Label Map Segmentation Image - } - if dataset.SOPClassUID not in allowed_sop_class_uids: - raise ValueError( - 'Expected a VL Whole Slide Microscopy Image or Segmentation Image.' - ) + # The "channels" output is either segment for segmentations, or optical + # path for other images + match dataset.SOPClassUID: + case ( + '1.2.840.10008.5.1.4.1.1.66.7' | # Label Map Segmentation + '1.2.840.10008.5.1.4.1.1.30' # Parametric Map + ): + # No "channel" in this case -> return None + channels = [None] + + case '1.2.840.10008.5.1.4.1.1.66.4': # Segmentation Image + channels = range(1, len(dataset.SegmentSequence) + 1) + + case '1.2.840.10008.5.1.4.1.1.77.1.6': # VL Whole Slide Microscopy + num_optical_paths = getattr( + dataset, + 'NumberOfOpticalPaths', + len(dataset.OpticalPathSequence) + ) + channels = range(1, num_optical_paths + 1) + + case _: + raise ValueError( + 'Expected a VL Whole Slide Microscopy Image, Segmentation ' + 'Image, or Parametric Map.' + ) + if ( not hasattr(dataset, "DimensionOrganizationType") or dataset.DimensionOrganizationType != "TILED_FULL" @@ -366,27 +384,6 @@ def iter_tiled_full_frame_data( 1 ) - is_segmentation = dataset.SOPClassUID in ( - '1.2.840.10008.5.1.4.1.1.66.4', - '1.2.840.10008.5.1.4.1.1.66.7', - ) - - # The "channels" output is either segment for segmentations, or optical - # path for other images - if is_segmentation: - if dataset.SegmentationType == "LABELMAP": - # No "channel" in this case -> return None - channels = [None] - else: - channels = range(1, len(dataset.SegmentSequence) + 1) - else: - num_optical_paths = getattr( - dataset, - 'NumberOfOpticalPaths', - len(dataset.OpticalPathSequence) - ) - channels = range(1, num_optical_paths + 1) - shared_fg = dataset.SharedFunctionalGroupsSequence[0] pixel_measures = shared_fg.PixelMeasuresSequence[0] pixel_spacing = ( diff --git a/src/highdicom/utils.py b/src/highdicom/utils.py index 4fdfaa24..a1ee94a3 100644 --- a/src/highdicom/utils.py +++ b/src/highdicom/utils.py @@ -106,7 +106,8 @@ def compute_plane_position_tiled_full( warnings.warn( "Passing a slice_thickness other than None has no effect and " "will be deprecated in a future version of the library.", - UserWarning + UserWarning, + stacklevel=2, ) if row_index < 1 or column_index < 1: raise ValueError("Row and column indices must be positive integers.") diff --git a/src/highdicom/volume.py b/src/highdicom/volume.py index 36ea4a79..af364949 100644 --- a/src/highdicom/volume.py +++ b/src/highdicom/volume.py @@ -39,6 +39,7 @@ from highdicom.uid import UID from pydicom.datadict import ( + dictionary_description, get_entry, tag_for_keyword, keyword_for_tag, @@ -188,6 +189,19 @@ def keyword(self) -> str: """str: The DICOM keyword or custom string for the descriptor.""" return self._keyword + @property + def description(self) -> str | None: + """str: The DICOM description for the descriptor. + + The description is a short textual description of the attribute taken + from the standard. + + ``None`` for custom descriptors. + + """ + if not self.is_custom: + return dictionary_description(self._tag) + @property def tag(self) -> BaseTag | None: """str: The DICOM tag for the attribute. @@ -2236,6 +2250,16 @@ def shape(self) -> tuple[int, ...]: """ return self.spatial_shape + @property + def ndim(self) -> int: + """int: Number of dimensions. + + For objects of type :class:`highdicom.VolumeGeometry`, this is + always 3. + + """ + return 3 + def __getitem__( self, index: int | slice | tuple[int | slice], @@ -2790,9 +2814,9 @@ def get_geometry(self) -> VolumeGeometry: ) @property - def dtype(self) -> type: + def dtype(self) -> np.dtype: """type: Datatype of the array.""" - return self._array.dtype.type + return self._array.dtype @property def shape(self) -> tuple[int, ...]: @@ -2812,6 +2836,15 @@ def spatial_shape(self) -> tuple[int, int, int]: """ return tuple(self._array.shape[:3]) + @property + def ndim(self) -> int: + """int: Number of dimensions. + + This includes spatial and channel dimensions. + + """ + return self._array.ndim + @property def number_of_channel_dimensions(self) -> int: """int: Number of channel dimensions.""" diff --git a/tests/test_image.py b/tests/test_image.py index 4f40e980..d63de354 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -12,6 +12,7 @@ import re from highdicom import ( + DimensionIndexSequence, Image, Volume, imread, @@ -19,25 +20,95 @@ from highdicom._module_utils import ( does_iod_have_pixel_data, ) -from highdicom.content import VOILUTTransformation +from highdicom.content import ( + _add_icc_profile_attributes, + VOILUTTransformation, +) from highdicom.image import ( _CombinedPixelTransform, ) from highdicom.pixels import ( apply_voi_window, ) -from highdicom.pr.content import ( - _add_icc_profile_attributes, -) from highdicom.pm import ( RealWorldValueMapping, ParametricMap, ) from highdicom.sr.coding import CodedConcept from highdicom.uid import UID +from highdicom.volume import ChannelDescriptor from tests.utils import find_readable_images +class TestDimensionIndexSequence(): + + def test_construction_1(self): + seq = DimensionIndexSequence( + coordinate_system='PATIENT', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), + ) + assert len(seq) == 1 + assert seq[0].DimensionIndexPointer == 0x00200032 + assert seq[0].FunctionalGroupPointer == 0x00209113 + + def test_construction_2(self): + seq = DimensionIndexSequence( + coordinate_system='PATIENT', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), + channel_dimensions=[ChannelDescriptor("ReferencedSegmentNumber")], + ) + assert len(seq) == 2 + assert seq[0].DimensionIndexPointer == 0x0062000B + assert seq[0].FunctionalGroupPointer == 0x0062000A + assert seq[1].DimensionIndexPointer == 0x00200032 + assert seq[1].FunctionalGroupPointer == 0x00209113 + + def test_construction_3(self): + seq = DimensionIndexSequence( + coordinate_system='SLIDE', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), + ) + assert len(seq) == 5 + assert seq[0].DimensionIndexPointer == 0x0048021F + assert seq[0].FunctionalGroupPointer == 0x0048021A + assert seq[1].DimensionIndexPointer == 0x0048021E + assert seq[1].FunctionalGroupPointer == 0x0048021A + assert seq[2].DimensionIndexPointer == 0x0040072A + assert seq[2].FunctionalGroupPointer == 0x0048021A + assert seq[3].DimensionIndexPointer == 0x0040073A + assert seq[3].FunctionalGroupPointer == 0x0048021A + assert seq[4].DimensionIndexPointer == 0x0040074A + assert seq[4].FunctionalGroupPointer == 0x0048021A + + def test_construction_4(self): + seq = DimensionIndexSequence( + coordinate_system='SLIDE', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), + channel_dimensions=[ChannelDescriptor("ReferencedSegmentNumber")], + ) + assert len(seq) == 6 + assert seq[0].DimensionIndexPointer == 0x0062000B + assert seq[0].FunctionalGroupPointer == 0x0062000A + assert seq[1].DimensionIndexPointer == 0x0048021F + assert seq[1].FunctionalGroupPointer == 0x0048021A + assert seq[2].DimensionIndexPointer == 0x0048021E + assert seq[2].FunctionalGroupPointer == 0x0048021A + assert seq[3].DimensionIndexPointer == 0x0040072A + assert seq[3].FunctionalGroupPointer == 0x0048021A + assert seq[4].DimensionIndexPointer == 0x0040073A + assert seq[4].FunctionalGroupPointer == 0x0048021A + assert seq[5].DimensionIndexPointer == 0x0040074A + assert seq[5].FunctionalGroupPointer == 0x0048021A + + def test_slice_spacing(): ct_multiframe = pydicom.dcmread( get_testdata_file('eCT_Supplemental.dcm') diff --git a/tests/test_pm.py b/tests/test_pm.py index da9f4dc6..fb431cde 100644 --- a/tests/test_pm.py +++ b/tests/test_pm.py @@ -1,11 +1,9 @@ -import unittest from collections import defaultdict from pathlib import Path from typing import Sequence import numpy as np import pytest -from pydicom import dcmread from pydicom.data import get_testdata_file, get_testdata_files from pydicom.sr.codedict import codes from pydicom.sr.coding import Code @@ -15,26 +13,36 @@ ) from highdicom import ( + PaletteColorLUT, PaletteColorLUTTransformation, PlanePositionSequence, PixelMeasuresSequence, PlaneOrientationSequence, + Volume, + imread, +) +from highdicom.content import VOILUTTransformation +from highdicom.enum import ( + ContentQualificationValues, + CoordinateSystemNames, + DimensionOrganizationTypeValues, + InterpolationMethods, ) -from highdicom.enum import ContentQualificationValues, CoordinateSystemNames from highdicom.pm.content import RealWorldValueMapping from highdicom.pm.enum import ( DerivedPixelContrastValues, ImageFlavorValues, ) from highdicom.pm.sop import ParametricMap +from highdicom.pm.pyramid import create_parametric_map_pyramid +from highdicom.spatial import sort_datasets from highdicom.uid import UID +from .utils import write_and_read_dataset -class TestRealWorldValueMapping(unittest.TestCase): - def setUp(self): - super().setUp() +class TestRealWorldValueMapping(): def test_failed_construction_missing_or_unnecessary_parameters(self): lut_label = '1' @@ -284,10 +292,10 @@ def test_failed_construction_floating_point_nonlinear_relationship(self): ) -class TestParametricMap(unittest.TestCase): +class TestParametricMap(): + @pytest.fixture(autouse=True) def setUp(self): - super().setUp() file_path = Path(__file__) data_dir = file_path.parent.parent.joinpath('data') @@ -302,39 +310,73 @@ def setUp(self): self._content_description = 'Test Parametric Map' self._content_creator_name = 'Will^I^Am' - self._ct_image = dcmread( + self._integer_mapping = RealWorldValueMapping( + lut_label='1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=[0, 255], + intercept=0, + slope=1 + ) + self._float_mapping = RealWorldValueMapping( + lut_label='1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=[-10000.0, 10000.0], + intercept=0, + slope=1 + ) + + self._voi_transformation = VOILUTTransformation( + window_width=128, + window_center=120, + ) + + self._ct_image = imread( str(data_dir.joinpath('test_files', 'ct_image.dcm')) ) - self._ct_multiframe_image = dcmread( + self._ct_multiframe_image = imread( get_testdata_file('eCT_Supplemental.dcm') ) - self._sm_image = dcmread( + self._sm_image = imread( str(data_dir.joinpath('test_files', 'sm_image.dcm')) ) ct_series = [ - dcmread(f) + imread(f) for f in get_testdata_files('dicomdirtests/77654033/CT2/*') ] - self._ct_series = sorted( - ct_series, - key=lambda x: x.ImagePositionPatient[2] + self._ct_series = sort_datasets(ct_series) + + self._ct_volume = ( + imread(get_testdata_file('eCT_Supplemental.dcm')) + .get_volume() ) @staticmethod - def check_dimension_index_vals(seg): + @pytest.fixture( + params=[ + DimensionOrganizationTypeValues.TILED_FULL, + DimensionOrganizationTypeValues.TILED_SPARSE + ] + ) + def tiled_dimension_organization(request): + return request.param + + @staticmethod + def check_dimension_index_vals(pm): # Function to apply some checks (necessary but not sufficient for # correctness) to ensure that the dimension indices are correct is_patient_coord_system = hasattr( - seg.PerFrameFunctionalGroupsSequence[0], + pm.PerFrameFunctionalGroupsSequence[0], 'PlanePositionSequence' ) if is_patient_coord_system: # Build up the mapping from index to value index_mapping = defaultdict(list) - for f in seg.PerFrameFunctionalGroupsSequence: + for f in pm.PerFrameFunctionalGroupsSequence: posn_index = f.FrameContentSequence[0].DimensionIndexValues[1] # This is not general, but all the tests run here use axial # images so just check the z coordinate @@ -361,7 +403,7 @@ def check_dimension_index_vals(seg): 'RowPositionInTotalImagePixelMatrix' ], [1, 2]): index_mapping = defaultdict(list) - for f in seg.PerFrameFunctionalGroupsSequence: + for f in pm.PerFrameFunctionalGroupsSequence: content_item = f.FrameContentSequence[0] posn_index = content_item.DimensionIndexValues[dim_ind] # This is not general, but all the tests run here use axial @@ -389,6 +431,14 @@ def test_multi_frame_sm_image_single_native(self): pixel_array = pixel_array.astype(np.float32) window_center = 0.5 window_width = 1.0 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -411,15 +461,10 @@ def test_multi_frame_sm_image_single_native(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, content_label=content_label ) - # Work around pydicom 3 decoding issue (should be able to remove this - # soon) - pmap.pixel_array_options(use_v2_backend=True) - assert pmap.SOPClassUID == '1.2.840.10008.5.1.4.1.1.30' assert pmap.SOPInstanceUID == self._sop_instance_uid assert pmap.SeriesInstanceUID == self._series_instance_uid @@ -448,6 +493,13 @@ def test_multi_frame_sm_image_single_native(self): voi_lut_item = sffg_item.FrameVOILUTSequence[0] assert voi_lut_item.WindowCenter == str(window_center) assert voi_lut_item.WindowWidth == str(window_width) + assert not hasattr(pmap, 'ICCProfile') + assert not hasattr(pmap, 'RedPaletteColorLookupTableDescriptor') + assert not hasattr(pmap, 'RedPaletteColorLookupTableData') + assert not hasattr(pmap, 'GreenPaletteColorLookupTableDescriptor') + assert not hasattr(pmap, 'GreenPaletteColorLookupTableData') + assert not hasattr(pmap, 'BluePaletteColorLookupTableDescriptor') + assert not hasattr(pmap, 'BluePaletteColorLookupTableData') def test_multi_frame_sm_image_ushort_native(self): pixel_array = np.random.randint( @@ -458,6 +510,14 @@ def test_multi_frame_sm_image_ushort_native(self): ) window_center = 128 window_width = 256 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -479,8 +539,7 @@ def test_multi_frame_sm_image_ushort_native(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, content_qualification=ContentQualificationValues.SERVICE, image_flavor=ImageFlavorValues.WHOLE_BODY, derived_pixel_contrast=DerivedPixelContrastValues.NONE @@ -488,6 +547,8 @@ def test_multi_frame_sm_image_ushort_native(self): sffg_item = instance.SharedFunctionalGroupsSequence[0] assert hasattr(sffg_item, 'RealWorldValueMappingSequence') assert len(sffg_item.RealWorldValueMappingSequence) == 1 + assert hasattr(sffg_item, 'PixelValueTransformationSequence') + assert hasattr(sffg_item, 'FrameVOILUTSequence') pffg_item = instance.PerFrameFunctionalGroupsSequence[0] assert not hasattr(pffg_item, 'RealWorldValueMappingSequence') assert instance.BitsAllocated == 8 @@ -499,6 +560,13 @@ def test_multi_frame_sm_image_ushort_native(self): assert instance.ImageType[2] == 'WHOLE_BODY' assert instance.ImageType[3] == 'NONE' assert instance.PixelPresentation == 'MONOCHROME' + assert not hasattr(instance, 'ICCProfile') + assert not hasattr(instance, 'RedPaletteColorLookupTableDescriptor') + assert not hasattr(instance, 'RedPaletteColorLookupTableData') + assert not hasattr(instance, 'GreenPaletteColorLookupTableDescriptor') + assert not hasattr(instance, 'GreenPaletteColorLookupTableData') + assert not hasattr(instance, 'BluePaletteColorLookupTableDescriptor') + assert not hasattr(instance, 'BluePaletteColorLookupTableData') def test_multi_frame_palette_lut(self): pixel_array = np.random.randint( @@ -509,6 +577,14 @@ def test_multi_frame_palette_lut(self): ) window_center = 128 window_width = 256 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -546,8 +622,7 @@ def test_multi_frame_palette_lut(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, content_qualification=ContentQualificationValues.SERVICE, image_flavor=ImageFlavorValues.WHOLE_BODY, derived_pixel_contrast=DerivedPixelContrastValues.NONE, @@ -586,6 +661,14 @@ def test_multi_frame_sm_image_ushort_palette_lut(self): ) window_center = 128 window_width = 256 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -607,8 +690,7 @@ def test_multi_frame_sm_image_ushort_palette_lut(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, content_qualification=ContentQualificationValues.SERVICE, image_flavor=ImageFlavorValues.WHOLE_BODY, derived_pixel_contrast=DerivedPixelContrastValues.NONE @@ -616,6 +698,8 @@ def test_multi_frame_sm_image_ushort_palette_lut(self): sffg_item = instance.SharedFunctionalGroupsSequence[0] assert hasattr(sffg_item, 'RealWorldValueMappingSequence') assert len(sffg_item.RealWorldValueMappingSequence) == 1 + assert hasattr(sffg_item, 'PixelValueTransformationSequence') + assert hasattr(sffg_item, 'FrameVOILUTSequence') pffg_item = instance.PerFrameFunctionalGroupsSequence[0] assert not hasattr(pffg_item, 'RealWorldValueMappingSequence') assert instance.BitsAllocated == 8 @@ -627,6 +711,13 @@ def test_multi_frame_sm_image_ushort_palette_lut(self): assert instance.ImageType[2] == 'WHOLE_BODY' assert instance.ImageType[3] == 'NONE' assert instance.PixelPresentation == 'MONOCHROME' + assert not hasattr(instance, 'ICCProfile') + assert not hasattr(instance, 'RedPaletteColorLookupTableDescriptor') + assert not hasattr(instance, 'RedPaletteColorLookupTableData') + assert not hasattr(instance, 'GreenPaletteColorLookupTableDescriptor') + assert not hasattr(instance, 'GreenPaletteColorLookupTableData') + assert not hasattr(instance, 'BluePaletteColorLookupTableDescriptor') + assert not hasattr(instance, 'BluePaletteColorLookupTableData') def test_multi_frame_sm_image_ushort_encapsulated_jpeg2000(self): pytest.importorskip("openjpeg") @@ -639,6 +730,13 @@ def test_multi_frame_sm_image_ushort_encapsulated_jpeg2000(self): window_center = 128 window_width = 256 + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -660,8 +758,7 @@ def test_multi_frame_sm_image_ushort_encapsulated_jpeg2000(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, transfer_syntax_uid=JPEG2000Lossless ) assert pmap.BitsAllocated == 8 @@ -678,6 +775,13 @@ def test_multi_frame_sm_image_ushort_encapsulated_jpegls(self): window_center = 128 window_width = 256 + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -699,8 +803,7 @@ def test_multi_frame_sm_image_ushort_encapsulated_jpegls(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, transfer_syntax_uid=JPEGLSLossless, use_extended_offset_table=True ) @@ -713,11 +816,19 @@ def test_single_frame_ct_image_double(self): pixel_array = np.random.uniform(-1, 1, self._ct_image.pixel_array.shape) window_center = 0.0 window_width = 2.0 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', unit=codes.UCUM.NoUnits, - value_range=[-1, 1], + value_range=(-1.0, 1.0), intercept=0, slope=1 ) @@ -734,16 +845,53 @@ def test_single_frame_ct_image_double(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, ) assert pmap.BitsAllocated == 64 + assert np.array_equal(pmap.pixel_array, pixel_array) - # Work around pydicom 3 decoding issue (should be able to remove this - # soon) - pmap.pixel_array_options(use_v2_backend=True) + def test_single_frame_ct_image_double_wrong_value_range(self): + pixel_array = np.random.uniform(-1, 1, self._ct_image.pixel_array.shape) + window_center = 0.0 + window_width = 2.0 - assert np.array_equal(pmap.pixel_array, pixel_array) + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + + real_world_value_mapping = RealWorldValueMapping( + lut_label='1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=(-1, 1), # should be float + intercept=0, + slope=1 + ) + + msg = ( + "When using a floating point-valued pixel_array, " + "all items in 'real_world_value_mappings' must have " + "their value range specified with floats." + ) + with pytest.raises(ValueError, match=msg): + ParametricMap( + [self._ct_image], + pixel_array, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[real_world_value_mapping], + voi_lut_transformations=voi_transformations, + ) def test_single_frame_ct_image_ushort_native(self): pixel_array = np.random.randint( @@ -754,6 +902,14 @@ def test_single_frame_ct_image_ushort_native(self): ) window_center = 2**12 / 2.0 window_width = 2**12 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -775,12 +931,59 @@ def test_single_frame_ct_image_ushort_native(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, ) assert pmap.BitsAllocated == 16 assert np.array_equal(pmap.pixel_array, pixel_array) + def test_single_frame_ct_image_ushort_wrong_value_range(self): + pixel_array = np.random.randint( + low=0, + high=2**12, + size=self._ct_image.pixel_array.shape, + dtype=np.uint16 + ) + window_center = 2**12 / 2.0 + window_width = 2**12 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + + real_world_value_mapping = RealWorldValueMapping( + lut_label='1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=(0.0, 4095.0), # should be int + intercept=0, + slope=1 + ) + + msg = ( + "When using an integer-valued 'pixel_array', all items " + "in 'real_world_value_mappings' must have their value " + "range specified with integers." + ) + with pytest.raises(ValueError, match=msg): + ParametricMap( + [self._ct_image], + pixel_array, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[real_world_value_mapping], + voi_lut_transformations=voi_transformations, + ) + def test_single_frame_ct_image_ushort(self): pixel_array = np.random.randint( low=120, @@ -790,6 +993,14 @@ def test_single_frame_ct_image_ushort(self): ) window_center = 2**16 / 2 window_width = 2**16 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', @@ -811,8 +1022,7 @@ def test_single_frame_ct_image_ushort(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, ) retrieved_pixel_array = pmap.pixel_array @@ -833,11 +1043,19 @@ def test_series_single_frame_ct_image_single(self): pixel_array = pixel_array.astype(np.float32) window_center = 0.0 window_width = 2.0 + + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + real_world_value_mapping = RealWorldValueMapping( lut_label='1', lut_explanation='feature_001', unit=codes.UCUM.NoUnits, - value_range=[-1, 1], + value_range=(-1.0, 1.0), intercept=0, slope=1 ) @@ -854,14 +1072,9 @@ def test_series_single_frame_ct_image_single(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, ) - # Work around pydicom 3 decoding issue (should be able to remove this - # soon) - pmap.pixel_array_options(use_v2_backend=True) - assert np.array_equal(pmap.pixel_array, pixel_array) def test_multi_frame_sm_image_with_spatial_positions_not_preserved(self): @@ -874,6 +1087,13 @@ def test_multi_frame_sm_image_with_spatial_positions_not_preserved(self): window_center = 128 window_width = 256 + voi_transformations = [ + VOILUTTransformation( + window_width=window_width, + window_center=window_center, + ) + ] + pixel_spacing = (0.5, 0.5) slice_thickness = 0.3 pixel_measures = PixelMeasuresSequence( @@ -916,8 +1136,7 @@ def test_multi_frame_sm_image_with_spatial_positions_not_preserved(self): self._device_serial_number, contains_recognizable_visual_features=False, real_world_value_mappings=[real_world_value_mapping], - window_center=window_center, - window_width=window_width, + voi_lut_transformations=voi_transformations, pixel_measures=pixel_measures, plane_orientation=plane_orientation, plane_positions=plane_positions @@ -932,3 +1151,305 @@ def test_multi_frame_sm_image_with_spatial_positions_not_preserved(self): assert not hasattr(shared_item, 'PlaneOrientationSequence') assert instance.ImageOrientationSlide == list(image_orientation) self.check_dimension_index_vals(instance) + + def test_from_volume(self): + # Creating a parametric map from a volume aligned with the source + # images + instance = ParametricMap( + [self._ct_multiframe_image], + self._ct_volume, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[self._float_mapping], + voi_lut_transformations=[self._voi_transformation], + ) + for frame_item in instance.PerFrameFunctionalGroupsSequence: + assert len(frame_item.DerivationImageSequence) == 1 + + instance = ParametricMap.from_dataset( + write_and_read_dataset(instance) + ) + + vol = instance.get_volume() + assert vol.geometry_equal(self._ct_volume) + assert np.array_equal(vol.array, self._ct_volume.array) + + def test_from_volume_multichannel(self): + # Creating a parametric map from a volume with multiple channels + # aligned with the source images + channels = {'LUTLabel': ['LUT1', 'LUT2']} + + channel_2 = 1000.0 - self._ct_volume.array + + array = np.stack([self._ct_volume.array, channel_2], -1) + volume = self._ct_volume.with_array(array, channels=channels) + + mapping1 = RealWorldValueMapping( + lut_label='LUT1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=[-10000.0, 10000.0], + intercept=0, + slope=1 + ) + mapping2 = RealWorldValueMapping( + lut_label='LUT2', + lut_explanation='feature_002', + unit=codes.UCUM.NoUnits, + value_range=[-10000.0, 10000.0], + intercept=0, + slope=1 + ) + + # Should fail if the LUTLabels are mimatched between the volume and the + # real_world_value_mappings parameter + msg = ( + "The LUTLabels of the 'real_world_value_mappings' " + "must match those within the channel indentifiers " + "of the 'pixel_array'." + ) + with pytest.raises(ValueError, match=msg): + ParametricMap( + [self._ct_multiframe_image], + volume, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[[mapping2], [mapping1]], + voi_lut_transformations=[self._voi_transformation], + ) + + instance = ParametricMap( + [self._ct_multiframe_image], + volume, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[[mapping1], [mapping2]], + voi_lut_transformations=[self._voi_transformation], + ) + for frame_item in instance.PerFrameFunctionalGroupsSequence: + assert len(frame_item.DerivationImageSequence) == 1 + + instance = ParametricMap.from_dataset( + write_and_read_dataset(instance) + ) + + vol = instance.get_volume() + assert vol.geometry_equal(self._ct_volume) + assert np.array_equal(vol.array, array) + + def test_from_volume_non_aligned(self): + # Creating a parametric map from a volume that is not aligned with the + # source images + volume = Volume( + array=self._ct_volume.array, + affine=np.eye(4), + coordinate_system='PATIENT', + ) + + instance = ParametricMap( + [self._ct_multiframe_image], + volume, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[self._float_mapping], + voi_lut_transformations=[self._voi_transformation], + ) + for frame_item in instance.PerFrameFunctionalGroupsSequence: + assert len(frame_item.DerivationImageSequence) == 0 + + instance = ParametricMap.from_dataset( + write_and_read_dataset(instance) + ) + + vol = instance.get_volume() + assert vol.geometry_equal(volume) + assert np.array_equal(vol.array, volume.array) + + def test_autotile(self, tiled_dimension_organization): + # Creating a parametric map from a total pixel matrix + tpm = self._sm_image.get_total_pixel_matrix().mean(axis=-1) + + instance = ParametricMap( + [self._sm_image], + tpm, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[self._float_mapping], + voi_lut_transformations=[self._voi_transformation], + tile_pixel_array=True, + dimension_organization_type=tiled_dimension_organization, + ) + + instance = ParametricMap.from_dataset( + write_and_read_dataset(instance) + ) + + new_tpm = instance.get_total_pixel_matrix() + assert np.array_equal(new_tpm, tpm) + + volume = instance.get_volume() + assert np.array_equal(volume.array[0], tpm) + + def test_invalid_tiled_full(self): + # Cannot use TILED_FULL when there are multiple channels + tpm = self._sm_image.get_total_pixel_matrix()[None, :, :, :2] + + mapping1 = RealWorldValueMapping( + lut_label='LUT1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=[-10000.0, 10000.0], + intercept=0, + slope=1 + ) + mapping2 = RealWorldValueMapping( + lut_label='LUT2', + lut_explanation='feature_002', + unit=codes.UCUM.NoUnits, + value_range=[-10000.0, 10000.0], + intercept=0, + slope=1 + ) + + msg = ( + 'A value of "TILED_FULL" for parameter ' + '"dimension_organization_type" is not permitted ' + 'because the image contains multiple channels. See ' + 'https://dicom.nema.org/medical/dicom/current/output/' + 'chtml/part03/sect_C.7.6.17.3.html#sect_C.7.6.17.3.' + ) + + with pytest.raises(ValueError, match=msg): + ParametricMap( + [self._sm_image], + tpm, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[[mapping1], [mapping2]], + voi_lut_transformations=[self._voi_transformation], + tile_pixel_array=True, + dimension_organization_type="TILED_FULL", + ) + + +class TestParametricMapPyramid(): + + @pytest.fixture(autouse=True) + def setUp(self): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data') + self._sm_image = imread( + str(data_dir.joinpath('test_files', 'sm_image.dcm')) + ) + self._pixel_array = ( + self._sm_image + .get_total_pixel_matrix().mean(axis=-1) + [None] + ) + + self._series_instance_uid = UID() + self._series_number = 76 + self._manufacturer = 'MyManufacturer' + self._manufacturer_model_name = 'MyModel' + self._software_versions = 'v1.0' + self._device_serial_number = '1-2-3' + self._content_description = 'test parametric map' + self._content_creator_name = 'will^i^am' + self._pyramid_uid = UID() + self._pyramid_label = 'Giza001' + self._series_description = 'A test pyramid' + + self._float_mapping = RealWorldValueMapping( + lut_label='1', + lut_explanation='feature_001', + unit=codes.UCUM.NoUnits, + value_range=(-10000.0, 10000.0), + intercept=0, + slope=1 + ) + + self._voi_transformation = VOILUTTransformation( + window_width=128, + window_center=120, + ) + + def test_contruct_pm_pyramid(self): + pmaps = create_parametric_map_pyramid( + source_images=[self._sm_image], + pixel_arrays=[self._pixel_array], + interpolator=InterpolationMethods.LINEAR, + series_instance_uid=self._series_instance_uid, + series_number=self._series_number, + manufacturer=self._manufacturer, + manufacturer_model_name=self._manufacturer_model_name, + software_versions=self._software_versions, + device_serial_number=self._device_serial_number, + contains_recognizable_visual_features=False, + real_world_value_mappings=[self._float_mapping], + voi_lut_transformations=[self._voi_transformation], + downsample_factors=[5.0], + pyramid_uid=self._pyramid_uid, + pyramid_label=self._pyramid_label, + series_description=self._series_description, + ) + + assert len(pmaps) == 2 + + for pm in pmaps: + assert isinstance(pm, ParametricMap) + assert pm.SeriesInstanceUID == self._series_instance_uid + assert pm.SeriesNumber == self._series_number + assert pm.Manufacturer == self._manufacturer + assert ( + pm.ManufacturerModelName == + self._manufacturer_model_name + ) + assert pm.SoftwareVersions == self._software_versions + assert pm.DeviceSerialNumber == self._device_serial_number + assert pm.PyramidLabel == self._pyramid_label + assert pm.PyramidUID == self._pyramid_uid + + assert pmaps[1].TotalPixelMatrixRows == 10 + assert pmaps[1].TotalPixelMatrixColumns == 10 diff --git a/tests/test_seg.py b/tests/test_seg.py index 8ac46656..6a023301 100644 --- a/tests/test_seg.py +++ b/tests/test_seg.py @@ -31,6 +31,7 @@ PaletteColorLUT, PaletteColorLUTTransformation, ) +from highdicom.image import DimensionIndexSequence from highdicom.base_content import ContributingEquipment from highdicom.color import CIELabColor from highdicom.content import ( @@ -48,7 +49,6 @@ from highdicom.seg import ( create_segmentation_pyramid, segread, - DimensionIndexSequence, SegmentationTypeValues, SegmentAlgorithmTypeValues, Segmentation, @@ -549,40 +549,6 @@ def test_construction_wrong_attribute_enumated_value(self): ) -class TestDimensionIndexSequence(unittest.TestCase): - - def setUp(self): - super().setUp() - - def test_construction(self): - seq = DimensionIndexSequence( - coordinate_system='PATIENT' - ) - assert len(seq) == 2 - assert seq[0].DimensionIndexPointer == 0x0062000B - assert seq[0].FunctionalGroupPointer == 0x0062000A - assert seq[1].DimensionIndexPointer == 0x00200032 - assert seq[1].FunctionalGroupPointer == 0x00209113 - - def test_construction_2(self): - seq = DimensionIndexSequence( - coordinate_system='SLIDE' - ) - assert len(seq) == 6 - assert seq[0].DimensionIndexPointer == 0x0062000B - assert seq[0].FunctionalGroupPointer == 0x0062000A - assert seq[1].DimensionIndexPointer == 0x0048021F - assert seq[1].FunctionalGroupPointer == 0x0048021A - assert seq[2].DimensionIndexPointer == 0x0048021E - assert seq[2].FunctionalGroupPointer == 0x0048021A - assert seq[3].DimensionIndexPointer == 0x0040072A - assert seq[3].FunctionalGroupPointer == 0x0048021A - assert seq[4].DimensionIndexPointer == 0x0040073A - assert seq[4].FunctionalGroupPointer == 0x0048021A - assert seq[5].DimensionIndexPointer == 0x0040074A - assert seq[5].FunctionalGroupPointer == 0x0048021A - - class TestSegmentation: @pytest.fixture(autouse=True) @@ -868,7 +834,10 @@ def sort_frames(sources, mask): else: orientation = src.ImageOrientationPatient - dim_index = DimensionIndexSequence(coordinate_system) + dim_index = DimensionIndexSequence( + coordinate_system, + 'segmentation-multi-frame-functional-groups' + ) if hasattr(src, 'NumberOfFrames'): plane_positions = dim_index.get_plane_positions_of_image(src) else: @@ -897,6 +866,8 @@ def get_array_after_writing(instance): def check_dimension_index_vals(seg): # Function to apply some checks (necessary but not sufficient for # correctness) to ensure that the dimension indices are correct + sfgs = seg.SharedFunctionalGroupsSequence[0] + if seg.SegmentationType != "LABELMAP": all_segment_numbers = [] for f in seg.PerFrameFunctionalGroupsSequence: @@ -908,9 +879,18 @@ def check_dimension_index_vals(seg): # VM=1 so this is an int posn_index = dim_ind_vals - seg_number = ( - f.SegmentIdentificationSequence[0].ReferencedSegmentNumber - ) + if 'SegmentIdentificationSequence' in f: + seg_number = ( + f + .SegmentIdentificationSequence[0] + .ReferencedSegmentNumber + ) + else: + seg_number = ( + sfgs + .SegmentIdentificationSequence[0] + .ReferencedSegmentNumber + ) assert seg_number == posn_index all_segment_numbers.append(seg_number) @@ -1045,7 +1025,7 @@ def test_construction(self): assert instance.MaximumFractionalValue == 255 assert instance.ContentLabel == self._content_label assert instance.ContentDescription is None - assert instance.ContentCreatorName is None + assert 'ContentCreatorName' not in instance with pytest.raises(AttributeError): instance.LossyImageCompressionRatio # noqa: B018 with pytest.raises(AttributeError): @@ -1077,12 +1057,15 @@ def test_construction(self): po_item = shared_item.PlaneOrientationSequence[0] assert po_item.ImageOrientationPatient == \ self._ct_image.ImageOrientationPatient + assert len(shared_item.SegmentIdentificationSequence) == 1 + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert len(instance.DimensionOrganizationSequence) == 1 assert len(instance.DimensionIndexSequence) == 2 assert instance.NumberOfFrames == 1 assert len(instance.PerFrameFunctionalGroupsSequence) == 1 frame_item = instance.PerFrameFunctionalGroupsSequence[0] - assert len(frame_item.SegmentIdentificationSequence) == 1 + assert not hasattr(frame_item, 'SegmentIdentificationSequence') assert len(frame_item.FrameContentSequence) == 1 assert len(frame_item.DerivationImageSequence) == 1 assert len(frame_item.PlanePositionSequence) == 1 @@ -1161,6 +1144,9 @@ def test_construction_2(self): assert pm_item.PixelSpacing == src_pm_item.PixelSpacing assert pm_item.SliceThickness == src_pm_item.SliceThickness assert not hasattr(shared_item, "PlaneOrientationSequence") + assert len(shared_item.SegmentIdentificationSequence) == 1 + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert instance.ImageOrientationSlide == \ self._sm_image.ImageOrientationSlide assert instance.TotalPixelMatrixOriginSequence == \ @@ -1174,7 +1160,7 @@ def test_construction_2(self): assert instance.NumberOfFrames == num_frames assert len(instance.PerFrameFunctionalGroupsSequence) == num_frames frame_item = instance.PerFrameFunctionalGroupsSequence[0] - assert len(frame_item.SegmentIdentificationSequence) == 1 + assert not hasattr(frame_item, 'SegmentIdentificationSequence') assert len(frame_item.DerivationImageSequence) == 1 assert len(frame_item.FrameContentSequence) == 1 assert len(frame_item.PlanePositionSlideSequence) == 1 @@ -1242,6 +1228,9 @@ def test_construction_3(self): po_item = shared_item.PlaneOrientationSequence[0] assert po_item.ImageOrientationPatient == \ src_im.ImageOrientationPatient + assert len(shared_item.SegmentIdentificationSequence) == 1 + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert len(instance.DimensionOrganizationSequence) == 1 assert len(instance.DimensionIndexSequence) == 2 n_frames = len(self._ct_series_nonempty) @@ -1254,7 +1243,7 @@ def test_construction_3(self): ), 1 ): - assert len(frame_item.SegmentIdentificationSequence) == 1 + assert not hasattr(frame_item, 'SegmentIdentificationSequence') assert len(frame_item.FrameContentSequence) == 1 assert len(frame_item.DerivationImageSequence) == 1 assert len(frame_item.PlanePositionSequence) == 1 @@ -1342,12 +1331,15 @@ def test_construction_4(self): po_item = shared_item.PlaneOrientationSequence[0] assert po_item.ImageOrientationPatient == \ src_shared_item.PlaneOrientationSequence[0].ImageOrientationPatient + assert len(shared_item.SegmentIdentificationSequence) == 1 + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert len(instance.DimensionOrganizationSequence) == 1 assert len(instance.DimensionIndexSequence) == 2 assert len(instance.PerFrameFunctionalGroupsSequence) == \ self._ct_multiframe.NumberOfFrames frame_item = instance.PerFrameFunctionalGroupsSequence[0] - assert len(frame_item.SegmentIdentificationSequence) == 1 + assert not hasattr(frame_item, 'SegmentIdentificationSequence') assert len(frame_item.FrameContentSequence) == 1 assert len(frame_item.DerivationImageSequence) == 1 assert len(frame_item.PlanePositionSequence) == 1 @@ -1429,12 +1421,15 @@ def test_construction_5(self): po_item = shared_item.PlaneOrientationSequence[0] assert po_item.ImageOrientationPatient == \ src_im.ImageOrientationPatient + assert len(shared_item.SegmentIdentificationSequence) == 1 + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert len(instance.DimensionOrganizationSequence) == 1 assert len(instance.DimensionIndexSequence) == 2 assert instance.NumberOfFrames == 4 assert len(instance.PerFrameFunctionalGroupsSequence) == 4 frame_item = instance.PerFrameFunctionalGroupsSequence[0] - assert len(frame_item.SegmentIdentificationSequence) == 1 + assert not hasattr(frame_item, 'SegmentIdentificationSequence') assert len(frame_item.FrameContentSequence) == 1 assert len(frame_item.DerivationImageSequence) == 1 assert len(frame_item.PlanePositionSequence) == 1 @@ -1527,7 +1522,7 @@ def test_construction_6(self): assert instance.MaximumFractionalValue == 255 assert instance.ContentLabel == self._content_label assert instance.ContentDescription is None - assert instance.ContentCreatorName is None + assert 'ContentCreatorName' not in instance with pytest.raises(AttributeError): instance.LossyImageCompressionRatio # noqa: B018 with pytest.raises(AttributeError): @@ -1551,14 +1546,17 @@ def test_construction_6(self): assert instance.Columns == self._cr_image.pixel_array.shape[1] assert len(instance.SharedFunctionalGroupsSequence) == 1 shared_item = instance.SharedFunctionalGroupsSequence[0] + assert len(shared_item.SegmentIdentificationSequence) == 1 assert not hasattr(shared_item, 'PixelMeasuresSequence') assert not hasattr(shared_item, 'PlaneOrientationSequence') + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert len(instance.DimensionOrganizationSequence) == 1 assert len(instance.DimensionIndexSequence) == 1 assert instance.NumberOfFrames == 1 assert len(instance.PerFrameFunctionalGroupsSequence) == 1 frame_item = instance.PerFrameFunctionalGroupsSequence[0] - assert len(frame_item.SegmentIdentificationSequence) == 1 + assert not hasattr(frame_item, 'SegmentIdentificationSequence') assert len(frame_item.FrameContentSequence) == 1 assert len(frame_item.DerivationImageSequence) == 1 assert not hasattr(frame_item, 'PlanePositionSequence') @@ -1623,7 +1621,7 @@ def test_construction_7(self): assert instance.MaximumFractionalValue == 255 assert instance.ContentLabel == self._content_label assert instance.ContentDescription is None - assert instance.ContentCreatorName is None + assert 'ContentCreatorName' not in instance with pytest.raises(AttributeError): instance.LossyImageCompressionRatio # noqa: B018 with pytest.raises(AttributeError): @@ -1650,6 +1648,8 @@ def test_construction_7(self): shared_item = instance.SharedFunctionalGroupsSequence[0] assert not hasattr(shared_item, 'PixelMeasuresSequence') assert not hasattr(shared_item, 'PlaneOrientationSequence') + assert not hasattr(shared_item, 'PixelValueTransformationSequence') + assert not hasattr(shared_item, 'FrameVOILUTSequence') assert len(instance.DimensionOrganizationSequence) == 1 assert len(instance.DimensionIndexSequence) == 1 assert instance.NumberOfFrames == 2 @@ -1743,7 +1743,7 @@ def test_construction_9(self): ) assert instance.SOPClassUID == '1.2.840.10008.5.1.4.1.1.66.7' assert instance.PhotometricInterpretation == 'PALETTE COLOR' - assert hasattr(instance, 'ICCProfile') + assert not hasattr(instance, 'ICCProfile') assert hasattr(instance, 'RedPaletteColorLookupTableDescriptor') assert hasattr(instance, 'RedPaletteColorLookupTableData') assert hasattr(instance, 'GreenPaletteColorLookupTableDescriptor') @@ -4509,7 +4509,10 @@ def test_spatial_positions_not_preserved(self): def test_get_plane_positions_of_image_patient(self): seq = DimensionIndexSequence( - coordinate_system='PATIENT' + coordinate_system='PATIENT', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), ) plane_positions = seq.get_plane_positions_of_image(self._ct_multiframe) for position in plane_positions: @@ -4517,7 +4520,10 @@ def test_get_plane_positions_of_image_patient(self): def test_get_plane_positions_of_image_slide(self): seq = DimensionIndexSequence( - coordinate_system='SLIDE' + coordinate_system='SLIDE', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), ) plane_positions = seq.get_plane_positions_of_image(self._sm_image) for position in plane_positions: @@ -4525,7 +4531,10 @@ def test_get_plane_positions_of_image_slide(self): def test_get_plane_positions_of_series(self): seq = DimensionIndexSequence( - coordinate_system='PATIENT' + coordinate_system='PATIENT', + functional_groups_module=( + 'segmentation-multi-frame-functional-groups' + ), ) plane_positions = seq.get_plane_positions_of_series(self._ct_series) for position in plane_positions: diff --git a/tests/utils.py b/tests/utils.py index 79d5184f..371dc122 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -28,7 +28,9 @@ def write_and_read_dataset(dataset: Dataset): implicit_vr=implicit_vr, little_endian=little_endian, ) - return dcmread(fp, force=True) + dcm = dcmread(fp, force=True) + dcm.buffer = None # avoid warnings when deepcopied + return dcm def find_readable_images() -> list[tuple[str, str | None]]: