diff --git a/data/test_files/README.md b/data/test_files/README.md new file mode 100644 index 00000000..7a161a06 --- /dev/null +++ b/data/test_files/README.md @@ -0,0 +1,29 @@ +# Test Files + +These files are used for automated tests of the highdicom library. Note that +many more test files are available as part of the pydicom package. + +### Images +* `ct_image.dcm` - A small CT image. +* `sm_image.dcm` - A slide microscopy image. + +### Segmentation Images +* `seg_image_ct_binary.dcm` - Segmentation image of `BINARY` type derived + from the series of small CT images in the pydicom test files + (`dicomdirtests/77654033/CT2/*`) using the highdicom library. +* `seg_image_ct_binary_overlap.dcm` - Segmentation image of `BINARY` type derived + from the series of small CT images in the pydicom test files + (`dicomdirtests/77654033/CT2/*`) using the highdicom library. This example + contains 2 overlapping segments. +* `seg_image_ct_binary_fractional.dcm` - Segmentation image of `FRACTIONAL` + type but with binary values (i.e. only 0.0 and 1.0) derived from the series + of small CT images in the pydicom test files (`dicomdirtests/77654033/CT2/*`) + using the highdicom library. +* `seg_image_ct_true_fractional.dcm` - Segmentation image of `FRACTIONAL` + type with true fractional values derived from the series + of small CT images in the pydicom test files (`dicomdirtests/77654033/CT2/*`) + using the highdicom library. + +### Structured Report Documents +* `sr_document.dcm` - A simple SR document following TID 1500 describing + measurements of a CT image. diff --git a/data/test_files/seg_image_ct_binary.dcm b/data/test_files/seg_image_ct_binary.dcm new file mode 100644 index 00000000..04a19ae0 Binary files /dev/null and b/data/test_files/seg_image_ct_binary.dcm differ diff --git a/data/test_files/seg_image_ct_binary_fractional.dcm b/data/test_files/seg_image_ct_binary_fractional.dcm new file mode 100644 index 00000000..90ad2b77 Binary files /dev/null and b/data/test_files/seg_image_ct_binary_fractional.dcm differ diff --git a/data/test_files/seg_image_ct_binary_overlap.dcm b/data/test_files/seg_image_ct_binary_overlap.dcm new file mode 100644 index 00000000..a12fc40c Binary files /dev/null and b/data/test_files/seg_image_ct_binary_overlap.dcm differ diff --git a/data/test_files/seg_image_ct_true_fractional.dcm b/data/test_files/seg_image_ct_true_fractional.dcm new file mode 100644 index 00000000..dbd114f3 Binary files /dev/null and b/data/test_files/seg_image_ct_true_fractional.dcm differ diff --git a/src/highdicom/_module_utils.py b/src/highdicom/_module_utils.py new file mode 100644 index 00000000..57357f47 --- /dev/null +++ b/src/highdicom/_module_utils.py @@ -0,0 +1,169 @@ +from enum import Enum +from typing import Any, Dict, List, Optional, Sequence + +from pydicom import Dataset + +from highdicom._modules import MODULE_ATTRIBUTE_MAP + + +# Allowed values for the type of an attribute +class AttributeTypeValues(Enum): + + """Enumerated values for the type of an attribute.""" + + REQUIRED = '1' + CONDITIONALLY_REQUIRED = '1C' + REQUIRED_EMPTY_IF_UNKNOWN = '2' + CONDITIONALLY_REQUIRED_EMPTY_IF_UNKNOWN = '2C' + OPTIONAL = '3' + + +def check_required_attributes( + dataset: Dataset, + module: str, + base_path: Optional[Sequence[str]] = None, + recursive: bool = True, + check_optional_sequences: bool = True +) -> None: + """Check that a dataset contains a module's required attributes. + + This may be used to check a top-level dataset, or for a dataset + representing an element of a sequence at an arbitrary level of + nesting, as specified by passing the path parameter. + + The function does not return anything, but throws an AttributeError + if the checks fail. + + Parameters + ---------- + dataset: pydicom.dataset.Dataset + The dataset to be checked. + module: str + Name of the module whose attributes should be checked. + base_path: Optional[Sequence[str]] + Path within the module that the dataset is intended to occupy, + represented as a sequence of strings with each string representing + the name of a data element sequence. If omitted, the dataset is + assumed to represent the top-level dataset. + recursive: bool + If True (the default), the attributes within nested data element + sequences will also be checked. If False, only attributes expected + in the top level of the passed dataset will be checked. + check_optional_sequences: bool + If True, the required attributes of an optional sequence will be + checked if the optional sequence is present in the dataset or + any of its nested sequences. This is ignored if recursive is + False. + + Raises + ------ + AttributeError + If any of the required (type 1 or 2) attributes are not present + in the dataset for the given module. + + Notes + ----- + This function merely checks for the presence of required attributes. + It does not check whether the data elements are empty or not, whether + there are additional, invalid attributes, or whether the values of the + data elements are allowed. Furthermore, it does not attempt to + check for conditionally required attributes. Therefore, this check + represents a necessary but not sufficient condition for a dataset + to be valid according to the DICOM standard. + + """ + # Only check for type 1 and type 2 attributes + types_to_check = [ + AttributeTypeValues.REQUIRED, + AttributeTypeValues.REQUIRED_EMPTY_IF_UNKNOWN + ] + + # Construct tree once and re-use in all recursive calls + tree = construct_module_tree(module) + + if base_path is not None: + for p in base_path: + try: + tree = tree['attributes'][p] + except KeyError: + raise AttributeError(f"Invalid base path: {base_path}.") + + # Define recursive function + def check( + dataset: Dataset, + subtree: Dict[str, Any], + path: List[str] + ) -> None: + for kw, item in subtree['attributes'].items(): + required = item['type'] in types_to_check + if required: + if not hasattr(dataset, kw): + if len(path) > 0: + msg = ( + "Dataset does not have required attribute " + f"'{kw}' at path {path}" + ) + else: + msg = ( + "Dataset does not have required attribute " + f"'{kw}'." + ) + raise AttributeError( + msg + ) + if recursive: + sequence_exists = ( + 'attributes' in subtree['attributes'][kw] and + hasattr(dataset, kw) + ) + if required or (sequence_exists and check_optional_sequences): + # Recurse down to the next level of the tree, if it exists + new_subtree = subtree['attributes'][kw] + if 'attributes' in new_subtree: + # Need to perform the check on all elements of the + # sequence + for elem in dataset[kw]: + check( + elem, + subtree=new_subtree, + path=path + [kw] + ) + + # Kick off recursion + check(dataset, tree, []) + + +def construct_module_tree(module: str) -> Dict[str, Any]: + """Return module attributes arranged in a tree structure. + + Parameters + ---------- + module: str + Name of the module. + + Returns + ------- + Dict[str, Any] + Tree-structured representation of the module attributes in the form + of nested dictionaries. Each level of the tree consists of a dictionary + with up to two keys: 'type' and 'attributes'. 'type' maps to a + :class:AttributeTypeValues that describes the type of the attribute + (1, 1C, 2, 3), and is present at every level except the root. + 'attributes' is present for any attribute that is a sequence + containing other attributes, and maps attribute keywords to a + dictionary that forms an item in the next level of the tree structure. + + """ + if module not in MODULE_ATTRIBUTE_MAP: + raise AttributeError(f"No such module found: '{module}'.") + tree: Dict[str, Any] = {'attributes': {}} + for item in MODULE_ATTRIBUTE_MAP[module]: + location = tree['attributes'] + for p in item['path']: + if 'attributes' not in location[p]: + location[p]['attributes'] = {} + location = location[p]['attributes'] + location[item['keyword']] = { + 'type': AttributeTypeValues(item['type']) + } + return tree diff --git a/src/highdicom/content.py b/src/highdicom/content.py index ce878302..790cd7b3 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -1,5 +1,6 @@ """Generic Data Elements that can be included in a variety of IODs.""" import datetime +from copy import deepcopy from typing import Any, Dict, List, Optional, Union, Sequence, Tuple import numpy as np @@ -20,6 +21,7 @@ NumContentItem, TextContentItem, ) +from highdicom._module_utils import check_required_attributes class AlgorithmIdentificationSequence(DataElementSequence): @@ -47,7 +49,7 @@ def __init__( Version of the algorithm source: str, optional Source of the algorithm, e.g. name of the algorithm manufacturer - parameters: Dict[str: str], optional + parameters: Dict[str, str], optional Name and actual value of the parameters with which the algorithm was invoked @@ -73,6 +75,86 @@ def __init__( ]) self.append(item) + @classmethod + def from_sequence( + cls, + sequence: DataElementSequence + ) -> 'AlgorithmIdentificationSequence': + """Construct instance from an existing data element sequence. + + Parameters + ---------- + sequence: pydicom.sequence.Sequence + Data element sequence representing the + AlgorithmIdentificationSequence Sequence. + + Returns + ------- + highdicom.seg.content.AlgorithmIdentificationSequence + Algorithm identification sequence. + + """ + if not isinstance(sequence, DataElementSequence): + raise TypeError( + 'Sequence should be of type pydicom.sequence.Sequence.' + ) + if len(sequence) != 1: + raise ValueError('Sequence should contain a single item.') + check_required_attributes( + sequence[0], + module='segmentation-image', + base_path=[ + 'SegmentSequence', + 'SegmentationAlgorithmIdentificationSequence' + ] + ) + algo_id_sequence = deepcopy(sequence) + algo_id_sequence.__class__ = cls + return algo_id_sequence + + @property + def name(self) -> str: + """str: Name of the algorithm.""" + return self[0].AlgorithmName + + @property + def family(self) -> CodedConcept: + """highdicom.sr.coding.CodedConcept: Kind of the algorithm family.""" + return CodedConcept.from_dataset( + self[0].AlgorithmFamilyCodeSequence[0] + ) + + @property + def version(self) -> str: + """str: Version of the algorithm.""" + return self[0].AlgorithmVersion + + @property + def source(self) -> Optional[str]: + """Union[str, None]: + Source of the algorithm, e.g. name of the algorithm + manufacturer, if any + + """ + return getattr(self[0], 'AlgorithmSource', None) + + @property + def parameters(self) -> Optional[Dict[str, str]]: + """Union[Dict[str, str], None]: + Dictionary mapping algorithm parameter names to values, + if any + + """ + if not hasattr(self[0], 'AlgorithmParameters'): + return None + parameters = {} + for param in self[0].AlgorithmParameters.split(','): + split = param.split('=') + if len(split) != 2: + raise ValueError('Malformed parameter string') + parameters[split[0]] = split[1] + return parameters + class PixelMeasuresSequence(DataElementSequence): @@ -109,6 +191,51 @@ def __init__( item.SpacingBetweenSlices = spacing_between_slices self.append(item) + @classmethod + def from_sequence( + cls, + sequence: DataElementSequence + ) -> 'PixelMeasuresSequence': + """Create a PixelMeasuresSequence from an existing Sequence. + + Parameters + ---------- + sequence: pydicom.sequence.Sequence + Sequence to be converted. + + Returns + ------- + highdicom.PixelMeasuresSequence + Plane position sequence. + + Raises + ------ + TypeError: + If sequence is not of the correct type. + ValueError: + If sequence does not contain exactly one item. + AttributeError: + If sequence does not contain the attributes required for a + pixel measures sequence. + + """ + if not isinstance(sequence, DataElementSequence): + raise TypeError( + 'Sequence must be of type pydicom.sequence.Sequence' + ) + if len(sequence) != 1: + raise ValueError('Sequence must contain a single item.') + req_kws = ['SliceThickness', 'PixelSpacing'] + if not all(hasattr(sequence[0], kw) for kw in req_kws): + raise AttributeError( + 'Sequence does not have the required attributes for ' + 'a Pixel Measures Sequence.' + ) + + pixel_measures = deepcopy(sequence) + pixel_measures.__class__ = cls + return pixel_measures + class PlanePositionSequence(DataElementSequence): @@ -212,6 +339,56 @@ def __eq__(self, other: Any) -> bool: ]), ) + @classmethod + def from_sequence( + cls, + sequence: DataElementSequence + ) -> 'PlanePositionSequence': + """Create a PlanePositionSequence from an existing Sequence. + + The coordinate system is inferred from the attributes in the sequence. + + Parameters + ---------- + sequence: pydicom.sequence.Sequence + Sequence to be converted. + + Returns + ------- + highdicom.PlanePositionSequence: + Plane position sequence. + + Raises + ------ + TypeError: + If sequence is not of the correct type. + ValueError: + If sequence does not contain exactly one item. + AttributeError: + If sequence does not contain the attributes required for a + plane position sequence. + + """ + if not isinstance(sequence, DataElementSequence): + raise TypeError( + 'Sequence must be of type pydicom.sequence.Sequence' + ) + if len(sequence) != 1: + raise ValueError('Sequence must contain a single item.') + if not hasattr(sequence[0], 'ImagePositionPatient'): + check_required_attributes( + dataset=sequence[0], + module='segmentation-multi-frame-functional-groups', + base_path=[ + 'PerFrameFunctionalGroupsSequence', + 'PlanePositionSlideSequence' + ] + ) + + plane_position = deepcopy(sequence) + plane_position.__class__ = cls + return plane_position + class PlaneOrientationSequence(DataElementSequence): @@ -292,6 +469,53 @@ def __eq__(self, other: Any) -> bool: else: return False + @classmethod + def from_sequence( + cls, + sequence: DataElementSequence + ) -> 'PlaneOrientationSequence': + """Create a PlaneOrientationSequence from an existing Sequence. + + The coordinate system is inferred from the attributes in the sequence. + + Parameters + ---------- + sequence: pydicom.sequence.Sequence + Sequence to be converted. + + Returns + ------- + highdicom.PlaneOrientationSequence: + Plane orientation sequence. + + Raises + ------ + TypeError: + If sequence is not of the correct type. + ValueError: + If sequence does not contain exactly one item. + AttributeError: + If sequence does not contain the attributes required for a + plane orientation sequence. + + """ + if not isinstance(sequence, DataElementSequence): + raise TypeError( + 'Sequence must be of type pydicom.sequence.Sequence' + ) + if len(sequence) != 1: + raise ValueError('Sequence must contain a single item.') + if not hasattr(sequence[0], 'ImageOrientationPatient'): + if not hasattr(sequence[0], 'ImageOrientationSlide'): + raise AttributeError( + 'The sequence does not contain required attributes for ' + 'either the PATIENT or SLIDE coordinate system.' + ) + + plane_orientation = deepcopy(sequence) + plane_orientation.__class__ = cls + return plane_orientation + class IssuerOfIdentifier(Dataset): diff --git a/src/highdicom/seg/__init__.py b/src/highdicom/seg/__init__.py index a71c68c4..cd15d75d 100644 --- a/src/highdicom/seg/__init__.py +++ b/src/highdicom/seg/__init__.py @@ -1,5 +1,5 @@ """Package for creation of Segmentation (SEG) instances.""" -from highdicom.seg.sop import Segmentation +from highdicom.seg.sop import Segmentation, segread from highdicom.seg.enum import ( SegmentAlgorithmTypeValues, SegmentationTypeValues, @@ -21,6 +21,7 @@ __all__ = [ 'DimensionIndexSequence', 'Segmentation', + 'segread', 'SegmentAlgorithmTypeValues', 'SegmentationFractionalTypeValues', 'SegmentationTypeValues', diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index 99f5fc06..bece510f 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -1,4 +1,5 @@ """Data Elements that are specific to the Segmentation IOD.""" +from copy import deepcopy from typing import List, Optional, Sequence, Tuple, Union import numpy as np @@ -15,6 +16,7 @@ from highdicom.seg.enum import SegmentAlgorithmTypeValues from highdicom.sr.coding import CodedConcept from highdicom.utils import compute_plane_position_slide_per_frame +from highdicom._module_utils import check_required_attributes class SegmentDescription(Dataset): @@ -152,6 +154,146 @@ def __init__( for structure in primary_anatomic_structures ] + @classmethod + def from_dataset(cls, dataset: Dataset) -> 'SegmentDescription': + """Construct instance from an existing dataset. + + Parameters + ---------- + dataset: pydicom.dataset.Dataset + Dataset representing an item of the Segment Sequence. + + Returns + ------- + highdicom.seg.SegmentDescription + Segment description. + + """ + if not isinstance(dataset, Dataset): + raise TypeError( + 'Dataset must be of type pydicom.dataset.Dataset.' + ) + check_required_attributes( + dataset, + module='segmentation-image', + base_path=['SegmentSequence'] + ) + desc = deepcopy(dataset) + desc.__class__ = cls + + # Convert sub sequences to highdicom types + desc.SegmentedPropertyCategoryCodeSequence = [ + CodedConcept.from_dataset( + desc.SegmentedPropertyCategoryCodeSequence[0] + ) + ] + desc.SegmentedPropertyTypeCodeSequence = [ + CodedConcept.from_dataset( + desc.SegmentedPropertyTypeCodeSequence[0] + ) + ] + if hasattr(desc, 'SegmentationAlgorithmIdentificationSequence'): + desc.SegmentationAlgorithmIdentificationSequence = \ + AlgorithmIdentificationSequence.from_sequence( + desc.SegmentationAlgorithmIdentificationSequence + ) + if hasattr(desc, 'AnatomicRegionSequence'): + desc.AnatomicRegionSequence = [ + CodedConcept.from_dataset(ds) + for ds in desc.AnatomicRegionSequence + ] + if hasattr(desc, 'PrimaryAnatomicStructureSequence'): + desc.PrimaryAnatomicStructureSequence = [ + CodedConcept.from_dataset(ds) + for ds in desc.PrimaryAnatomicStructureSequence + ] + return desc + + @property + def segment_number(self) -> int: + """int: Number of the segment.""" + return int(self.SegmentNumber) + + @property + def segment_label(self) -> str: + """str: Label of the segment.""" + return str(self.SegmentLabel) + + @property + def segmented_property_category(self) -> CodedConcept: + """highdicom.sr.CodedConcept: + Category of the property the segment represents. + + """ + return self.SegmentedPropertyCategoryCodeSequence[0] + + @property + def segmented_property_type(self) -> CodedConcept: + """highdicom.sr.CodedConcept: + Type of the property the segment represents. + + """ + return self.SegmentedPropertyTypeCodeSequence[0] + + @property + def algorithm_type(self) -> SegmentAlgorithmTypeValues: + """highdicom.seg.SegmentAlgorithmTypeValues: + Type of algorithm used to create the segment. + + """ + return SegmentAlgorithmTypeValues(self.SegmentAlgorithmType) + + @property + def algorithm_identification( + self + ) -> Union[AlgorithmIdentificationSequence, None]: + """Union[highdicom.AlgorithmIdentificationSequence, None] + Information useful for identification of the algorithm, if any. + + """ + if hasattr(self, 'SegmentationAlgorithmIdentificationSequence'): + return self.SegmentationAlgorithmIdentificationSequence + return None + + @property + def tracking_uid(self) -> Union[str, None]: + """Union[str, None]: + Tracking unique identifier for the segment, if any. + + """ + if 'TrackingUID' in self: + return self.TrackingUID + return None + + @property + def tracking_id(self) -> Union[str, None]: + """Union[str, None]: Tracking identifier for the segment, if any.""" + if 'TrackingID' in self: + return self.TrackingID + return None + + @property + def anatomic_regions(self) -> List[CodedConcept]: + """List[highdicom.sr.CodedConcept]: + List of anatomic regions into which the segment falls. + May be empty. + + """ + if not hasattr(self, 'AnatomicRegionSequence'): + return [] + return list(self.AnatomicRegionSequence) + + @property + def primary_anatomic_structures(self) -> List[CodedConcept]: + """List[highdicom.sr.CodedConcept]: + List of anatomic anatomic structures the segment represents. + May be empty. + + """ + if not hasattr(self, 'PrimaryAnatomicStructureSequence'): + return [] + return list(self.PrimaryAnatomicStructureSequence) + class DimensionIndexSequence(DataElementSequence): @@ -291,7 +433,7 @@ def __init__( else: raise ValueError( - f'Unknown coordinate system "{self._coordinat_system}"' + f'Unknown coordinate system "{self._coordinate_system}"' ) def get_plane_positions_of_image( diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 6452b9cf..6935f6e3 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1,13 +1,21 @@ """Module for the SOP class of the Segmentation IOD.""" +from copy import deepcopy +from collections import OrderedDict import logging +from os import PathLike import numpy as np from collections import defaultdict -from typing import Any, Dict, List, Optional, Set, Sequence, Union, Tuple +from typing import ( + Any, Dict, List, Optional, Set, Sequence, Union, Tuple, BinaryIO +) from pydicom.dataset import Dataset +from pydicom.datadict import keyword_for_tag, tag_for_keyword from pydicom.encaps import decode_data_sequence, encapsulate +from pydicom.multival import MultiValue from pydicom.pixel_data_handlers.numpy_handler import pack_bits from pydicom.pixel_data_handlers.util import get_expected_length +from pydicom.tag import BaseTag, Tag from pydicom.uid import ( ExplicitVRLittleEndian, ImplicitVRLittleEndian, @@ -17,6 +25,8 @@ ) from pydicom.sr.codedict import codes from pydicom.valuerep import PersonName +from pydicom.sr.coding import Code +from pydicom.filereader import dcmread from highdicom.base import SOPClass from highdicom.content import ( @@ -34,14 +44,21 @@ SegmentationFractionalTypeValues, SegmentationTypeValues, SegmentsOverlapValues, + SpatialLocationsPreservedValues, + SegmentAlgorithmTypeValues, ) +from highdicom.seg.utils import iter_segments from highdicom.sr.coding import CodedConcept from highdicom.valuerep import check_person_name +from highdicom.uid import UID as hd_UID logger = logging.getLogger(__name__) +_NO_FRAME_REF_VALUE = -1 + + class Segmentation(SOPClass): """SOP class for a Segmentation, which represents one or more @@ -79,7 +96,7 @@ def __init__( """ Parameters ---------- - source_images: Sequence[pydicom.dataset.Dataset] + source_images: Sequence[Dataset] One or more single- or multi-frame images (or metadata of images) from which the segmentation was derived pixel_array: numpy.ndarray @@ -536,6 +553,11 @@ def add_segments( permitted in `segment_descriptions`. """ # noqa + if self._source_images is None: + raise AttributeError( + 'Further segments may not be added to Segmentation objects ' + 'created from existing datasets.' + ) if pixel_array.ndim == 2: pixel_array = pixel_array[np.newaxis, ...] if pixel_array.ndim != 3: @@ -893,6 +915,8 @@ def add_segments( if len(self.PixelData) % 2 == 1: self.PixelData += b'0' + self._build_luts() + def _encode_pixels(self, planes: np.ndarray) -> bytes: """Encodes pixel planes. @@ -940,3 +964,1663 @@ def _encode_pixels(self, planes: np.ndarray) -> bytes: return pack_bits(planes.flatten()) else: return planes.flatten().tobytes() + + @classmethod + def from_dataset(cls, dataset: Dataset) -> 'Segmentation': + """Create a Segmentation object from an existing pydicom dataset. + + Parameters + ---------- + dataset: pydicom.dataset.Dataset + Pydicom dataset representing a SEG image. + + Returns + ------- + highdicom.seg.Segmentation + Representation of the supplied dataset as a highdicom + Segmentation. + + """ + if not isinstance(dataset, Dataset): + raise TypeError( + 'Dataset must be of type pydicom.dataset.Dataset.' + ) + # Checks on integrity of input dataset + if dataset.SOPClassUID != '1.2.840.10008.5.1.4.1.1.66.4': + raise ValueError( + 'Dataset is not a Segmentation.' + ) + seg = deepcopy(dataset) + seg.__class__ = cls + + seg._source_images = None + seg._source_plane_orientation = None + sf_groups = seg.SharedFunctionalGroupsSequence[0] + plane_ori_seq = sf_groups.PlaneOrientationSequence[0] + if hasattr(plane_ori_seq, 'ImageOrientationSlide'): + seg._coordinate_system = CoordinateSystemNames.SLIDE + seg._plane_orientation = plane_ori_seq.ImageOrientationSlide + elif hasattr(plane_ori_seq, 'ImageOrientationPatient'): + seg._coordinate_system = CoordinateSystemNames.PATIENT + seg._plane_orientation = plane_ori_seq.ImageOrientationPatient + else: + raise AttributeError( + 'Expected Plane Orientation Sequence to have either ' + 'ImageOrientationSlide or ImageOrientationPatient ' + 'attribute.' + ) + + for i, segment in enumerate(seg.SegmentSequence, 1): + if segment.SegmentNumber != i: + raise AttributeError( + 'Segments are expected to start at 1 and be consecutive ' + 'integers.' + ) + + for i, s in enumerate(seg.SegmentSequence, 1): + if s.SegmentNumber != i: + raise ValueError( + 'Segment numbers in the segmentation image must start at ' + '1 and increase by 1 with the segments sequence.' + ) + # Needed for compatibility with add_segments + seg._segment_inventory = { + s.SegmentNumber for s in seg.SegmentSequence + } + + # Convert contained items to highdicom types + # Segment descriptions + seg.SegmentSequence = [ + SegmentDescription.from_dataset(ds) for ds in seg.SegmentSequence + ] + + # Shared functional group elements + if hasattr(sf_groups, 'PlanePositionSequence'): + plane_pos = PlanePositionSequence.from_sequence( + sf_groups.PlanePositionSequence + ) + sf_groups.PlanePositionSequence = plane_pos + if hasattr(sf_groups, 'PlaneOrientationSequence'): + plane_ori = PlaneOrientationSequence.from_sequence( + sf_groups.PlaneOrientationSequence + ) + sf_groups.PlaneOrientationSequence = plane_ori + if hasattr(sf_groups, 'PixelMeasuresSequence'): + pixel_measures = PixelMeasuresSequence.from_sequence( + sf_groups.PixelMeasuresSequence + ) + sf_groups.PixelMeasuresSequence = pixel_measures + + # Per-frame functional group items + for pffg_item in seg.PerFrameFunctionalGroupsSequence: + if hasattr(pffg_item, 'PlanePositionSequence'): + plane_pos = PlanePositionSequence.from_sequence( + pffg_item.PlanePositionSequence + ) + pffg_item.PlanePositionSequence = plane_pos + if hasattr(pffg_item, 'PlaneOrientationSequence'): + plane_ori = PlaneOrientationSequence.from_sequence( + pffg_item.PlaneOrientationSequence + ) + pffg_item.PlaneOrientationSequence = plane_ori + if hasattr(pffg_item, 'PixelMeasuresSequence'): + pixel_measures = PixelMeasuresSequence.from_sequence( + pffg_item.PixelMeasuresSequence + ) + pffg_item.PixelMeasuresSequence = pixel_measures + + seg._build_luts() + + return seg + + def _build_ref_instance_lut(self) -> None: + """Build lookup table for all instance referenced in the segmentation. + + Builds a lookup table mapping the SOPInstanceUIDs of all datasets + referenced in the segmentation to a tuple containing the + StudyInstanceUID, SeriesInstanceUID and SOPInstanceUID. + + """ + # Map sop uid to tuple (study uid, series uid, sop uid) + self._ref_ins_lut = OrderedDict() + if hasattr(self, 'ReferencedSeriesSequence'): + for ref_series in self.ReferencedSeriesSequence: + for ref_ins in ref_series.ReferencedInstanceSequence: + self._ref_ins_lut[ref_ins.ReferencedSOPInstanceUID] = ( + hd_UID(self.StudyInstanceUID), + hd_UID(ref_series.SeriesInstanceUID), + hd_UID(ref_ins.ReferencedSOPInstanceUID) + ) + other_studies_kw = 'StudiesContainingOtherReferencedInstancesSequence' + if hasattr(self, other_studies_kw): + for ref_study in getattr(self, other_studies_kw): + for ref_series in ref_study.ReferencedSeriesSequence: + for ref_ins in ref_series.ReferencedInstanceSequence: + self._ref_ins_lut[ref_ins.ReferencedSOPInstanceUID] = ( + hd_UID(ref_study.StudyInstanceUID), + hd_UID(ref_series.SeriesInstanceUID), + hd_UID(ref_ins.ReferencedSOPInstanceUID) + ) + + self._source_sop_instance_uids = np.array( + list(self._ref_ins_lut.keys()) + ) + + def _build_luts(self) -> None: + """Build lookup tables for efficient querying. + + Two lookup tables are currently constructed. The first maps the + SOPInstanceUIDs of all datasets referenced in the segmentation to a + tuple containing the StudyInstanceUID, SeriesInstanceUID and + SOPInstanceUID. + + The second look-up table contains information about each frame of the + segmentation, including the segment it contains, the instance and frame + from which it was derived (if these are unique), and its dimension + index values. + + """ + self._build_ref_instance_lut() + + segnum_col_data = [] + source_instance_col_data = [] + source_frame_col_data = [] + + # Get list of all dimension index pointers, excluding the segment + # number, since this is treated differently + seg_num_tag = tag_for_keyword('ReferencedSegmentNumber') + self._dim_ind_pointers = [ + dim_ind.DimensionIndexPointer + for dim_ind in self.DimensionIndexSequence + if dim_ind.DimensionIndexPointer != seg_num_tag + ] + dim_ind_positions = { + dim_ind.DimensionIndexPointer: i + for i, dim_ind in enumerate(self.DimensionIndexSequence) + if dim_ind.DimensionIndexPointer != seg_num_tag + } + dim_index_col_data: Dict[int, List[int]] = { + ptr: [] for ptr in self._dim_ind_pointers + } + + # Create a list of source images and check for spatial locations + # preserved and that there is a single source frame per seg frame + locations_list_type = List[Optional[SpatialLocationsPreservedValues]] + locations_preserved: locations_list_type = [] + self._single_source_frame_per_seg_frame = True + for frame_item in self.PerFrameFunctionalGroupsSequence: + # Get segment number for this frame + seg_id_seg = frame_item.SegmentIdentificationSequence[0] + seg_num = seg_id_seg.ReferencedSegmentNumber + segnum_col_data.append(int(seg_num)) + + # Get dimension indices for this frame + indices = frame_item.FrameContentSequence[0].DimensionIndexValues + if len(indices) != len(self._dim_ind_pointers) + 1: + # (+1 because referenced segment number is ignored) + raise RuntimeError( + 'Unexpected mismatch between dimension index values in ' + 'per-frames functional groups sequence and items in the ' + 'dimension index sequence.' + ) + for ptr in self._dim_ind_pointers: + dim_index_col_data[ptr].append(indices[dim_ind_positions[ptr]]) + + frame_source_instances = [] + frame_source_frames = [] + for der_im in frame_item.DerivationImageSequence: + for src_im in der_im.SourceImageSequence: + frame_source_instances.append( + src_im.ReferencedSOPInstanceUID + ) + if hasattr(src_im, 'SpatialLocationsPreserved'): + locations_preserved.append( + SpatialLocationsPreservedValues( + src_im.SpatialLocationsPreserved + ) + ) + else: + locations_preserved.append( + None + ) + + if hasattr(src_im, 'ReferencedFrameNumber'): + if isinstance( + src_im.ReferencedFrameNumber, + MultiValue + ): + frame_source_frames.extend( + [ + int(f) + for f in src_im.ReferencedFrameNumber + ] + ) + else: + frame_source_frames.append( + int(src_im.ReferencedFrameNumber) + ) + else: + frame_source_frames.append(_NO_FRAME_REF_VALUE) + + if ( + len(set(frame_source_instances)) != 1 or + len(set(frame_source_frames)) != 1 + ): + self._single_source_frame_per_seg_frame = False + else: + ref_instance_uid = frame_source_instances[0] + if ref_instance_uid not in self._source_sop_instance_uids: + raise AttributeError( + f'SOP instance {ref_instance_uid} referenced in the ' + 'source image sequence is not included in the ' + 'Referenced Series Sequence or Studies Containing ' + 'Other Referenced Instances Sequence.' + ) + src_uid_ind = self._get_src_uid_index(frame_source_instances[0]) + source_instance_col_data.append(src_uid_ind) + source_frame_col_data.append(frame_source_frames[0]) + + # Summarise + if any( + isinstance(v, SpatialLocationsPreservedValues) and + v == SpatialLocationsPreservedValues.NO + for v in locations_preserved + ): + Type = Optional[SpatialLocationsPreservedValues] + self._locations_preserved: Type = SpatialLocationsPreservedValues.NO + elif all( + isinstance(v, SpatialLocationsPreservedValues) and + v == SpatialLocationsPreservedValues.YES + for v in locations_preserved + ): + self._locations_preserved = SpatialLocationsPreservedValues.YES + else: + self._locations_preserved = None + + # Frame LUT is a 2D numpy array. Columns represent different ways to + # index segmentation frames. Row i represents frame i of the + # segmentation dataset Possible columns include the segment number, the + # source uid index in the source_sop_instance_uids list, the source + # frame number (if applicable), and the dimension index values of the + # segmentation frames + # This allows for fairly efficient querying by any of the three + # variables: seg frame number, source instance/frame number, segment + # number + # Column for segment number + self._lut_seg_col = 0 + col_data = [segnum_col_data] + + # Columns for other dimension index values + self._lut_dim_ind_cols = { + i: ptr for ptr, i in + enumerate(self._dim_ind_pointers, len(col_data)) + } + col_data += [ + indices for indices in dim_index_col_data.values() + ] + + # Columns related to source frames, if they are usable for indexing + if self._single_source_frame_per_seg_frame: + self._lut_src_instance_col = len(col_data) + self._lut_src_frame_col = len(col_data) + 1 + col_data += [ + source_instance_col_data, + source_frame_col_data + ] + + # Build LUT from columns + self._frame_lut = np.array(col_data).T + + @property + def segmentation_type(self) -> SegmentationTypeValues: + """highdicom.seg.SegmentationTypeValues: Segmentation type.""" + return SegmentationTypeValues(self.SegmentationType) + + @property + def segmentation_fractional_type( + self + ) -> Union[SegmentationFractionalTypeValues, None]: + """ + highdicom.seg.SegmentationFractionalTypeValues: + Segmentation fractional type. + + """ + if not hasattr(self, 'SegmentationFractionalType'): + return None + return SegmentationFractionalTypeValues( + self.SegmentationFractionalType + ) + + def iter_segments(self): + """Iterates over segments in this segmentation image. + + Returns + ------- + Iterator[Tuple[numpy.ndarray, Tuple[pydicom.dataset.Dataset, ...], pydicom.dataset.Dataset]] + For each segment in the Segmentation image instance, provides the + Pixel Data frames representing the segment, items of the Per-Frame + Functional Groups Sequence describing the individual frames, and + the item of the Segment Sequence describing the segment + + """ # noqa + return iter_segments(self) + + @property + def number_of_segments(self) -> int: + """int: The number of segments in this SEG image.""" + return len(self.SegmentSequence) + + @property + def segment_numbers(self) -> range: + """range: The segment numbers present in the SEG image as a range.""" + return range(1, self.number_of_segments + 1) + + def get_segment_description( + self, + segment_number: int + ) -> SegmentDescription: + """Get segment description for a segment. + + Parameters + ---------- + segment_number: int + Segment number for the segment, as a 1-based index. + + Returns + ------- + highdicom.seg.SegmentDescription + Description of the given segment. + + """ + if segment_number < 1 or segment_number > self.number_of_segments: + raise IndexError( + f'{segment_number} is an invalid segment number for this ' + 'dataset.' + ) + return self.SegmentSequence[segment_number - 1] + + def get_segment_numbers( + self, + segment_label: Optional[str] = None, + segmented_property_category: Optional[Union[Code, CodedConcept]] = None, + segmented_property_type: Optional[Union[Code, CodedConcept]] = None, + algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]] = None, + tracking_uid: Optional[str] = None, + tracking_id: Optional[str] = None, + ) -> List[int]: + """Get a list of segment numbers matching provided criteria. + + Any number of optional filters may be provided. A segment must match + all provided filters to be included in the returned list. + + Parameters + ---------- + segment_label: Union[str, None], optional + Segment label filter to apply. + segmented_property_category: Union[Code, CodedConcept, None], optional + Segmented property category filter to apply. + segmented_property_type: Union[Code, CodedConcept, None], optional + Segmented property type filter to apply. + algorithm_type: Union[SegmentAlgorithmTypeValues, str, None], optional + Segmented property type filter to apply. + tracking_uid: Union[str, None], optional + Tracking unique identifier filter to apply. + tracking_id: Union[str, None], optional + Tracking identifier filter to apply. + + Returns + ------- + List[int] + List of all segment numbers matching the provided criteria. + + Examples + -------- + + Get segment numbers of all segments that both represent tumors and were + generated by an automatic algorithm from a segmentation object ``seg``: + + >>> from pydicom.sr.codedict import codes + >>> from highdicom.seg import SegmentAlgorithmTypeValues + >>> + >>> segment_numbers = seg.get_segment_numbers( + >>> segmented_property_type=codes.SCT.Tumor, + >>> algorithm_type=SegmentAlgorithmTypeValues.AUTOMATIC + >>> ) + [1, 2, 3] + + Get segment numbers of all segments identified by a given + institution-specific tracking ID: + + >>> segment_numbers = seg.get_segment_numbers( + >>> tracking_id='Tumor #1234567' + >>> ) + [13] + + Get segment numbers of all segments identified a globally unique + tracking UID: + + >>> uid = '1.2.826.0.1.3680043.10.511.3.73025483745501512180439199223117347' + >>> segment_numbers = seg.get_segment_numbers(tracking_uid=uid) + [5] + + """ # noqa: E501 + filter_funcs = [] + if segment_label is not None: + filter_funcs.append( + lambda desc: desc.segment_label == segment_label + ) + if segmented_property_category is not None: + filter_funcs.append( + lambda desc: + desc.segmented_property_category == segmented_property_category + ) + if segmented_property_type is not None: + filter_funcs.append( + lambda desc: + desc.segmented_property_type == segmented_property_type + ) + if algorithm_type is not None: + algo_type = SegmentAlgorithmTypeValues(algorithm_type) + filter_funcs.append( + lambda desc: + SegmentAlgorithmTypeValues(desc.algorithm_type) == algo_type + ) + if tracking_uid is not None: + filter_funcs.append( + lambda desc: desc.tracking_uid == tracking_uid + ) + if tracking_id is not None: + filter_funcs.append( + lambda desc: desc.tracking_id == tracking_id + ) + + return [ + desc.segment_number + for desc in self.SegmentSequence + if all(f(desc) for f in filter_funcs) + ] + + def get_tracking_ids( + self, + segmented_property_category: Optional[Union[Code, CodedConcept]] = None, + segmented_property_type: Optional[Union[Code, CodedConcept]] = None, + algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]] = None + ) -> List[Tuple[str, UID]]: + """Get all unique tracking identifiers in this SEG image. + + Any number of optional filters may be provided. A segment must match + all provided filters to be included in the returned list. + + The tracking IDs and the accompanying tracking UIDs are returned + in a list of tuples. + + Note that the order of the returned list is not significant and will + not in general match the order of segments. + + Parameters + ---------- + segmented_property_category: Union[Code, CodedConcept, None], optional + Segmented property category filter to apply. + segmented_property_type: Union[Code, CodedConcept, None], optional + Segmented property type filter to apply. + algorithm_type: Union[SegmentAlgorithmTypeValues, str, None], optional + Segmented property type filter to apply. + + Returns + ------- + List[Tuple[str, UID]] + List of all unique (tracking_id, tracking_uid) tuples that are + referenced in segment descriptions in this SEG image that match + all provided filters. + + Examples + -------- + + Read in an example segmentation image in the highdicom test data: + + >>> import highdicom as hd + >>> from pydicom.sr.codedict import codes + >>> + >>> seg = hd.seg.segread('data/test_files/seg_image_ct_binary_overlap.dcm') + + List the tracking IDs and UIDs present in the segmentation image: + + >>> seg.get_tracking_ids() + [('Spine', '1.2.826.0.1.3680043.10.511.3.10042414969629429693880339016394772'), + ('Bone', '1.2.826.0.1.3680043.10.511.3.83271046815894549094043330632275067')] + + >>> for seg_num in seg.segment_numbers: + >>> desc = seg.get_segment_description(seg_num) + >>> print(desc.segmented_property_type.meaning) + Bone + Spine + + List tracking IDs only for those segments with a segmented property + category of 'Spine': + + >>> seg.get_tracking_ids(segmented_property_type=codes.SCT.Spine) + [('Spine', '1.2.826.0.1.3680043.10.511.3.10042414969629429693880339016394772')] + + """ # noqa: E501 + filter_funcs = [] + if segmented_property_category is not None: + filter_funcs.append( + lambda desc: + desc.segmented_property_category == segmented_property_category + ) + if segmented_property_type is not None: + filter_funcs.append( + lambda desc: + desc.segmented_property_type == segmented_property_type + ) + if algorithm_type is not None: + algo_type = SegmentAlgorithmTypeValues(algorithm_type) + filter_funcs.append( + lambda desc: + SegmentAlgorithmTypeValues(desc.algorithm_type) == algo_type + ) + + return list({ + (desc.tracking_id, UID(desc.tracking_uid)) + for desc in self.SegmentSequence + if desc.tracking_id is not None and + all(f(desc) for f in filter_funcs) + }) + + @property + def segmented_property_categories(self) -> List[CodedConcept]: + """Get all unique segmented property categories in this SEG image. + + Returns + ------- + List[CodedConcept] + All unique segmented property categories referenced in segment + descriptions in this SEG image. + + """ + categories = [] + for desc in self.SegmentSequence: + if desc.segmented_property_category not in categories: + categories.append(desc.segmented_property_category) + + return categories + + @property + def segmented_property_types(self) -> List[CodedConcept]: + """Get all unique segmented property types in this SEG image. + + Returns + ------- + List[CodedConcept] + All unique segmented property types referenced in segment + descriptions in this SEG image. + + """ + types = [] + for desc in self.SegmentSequence: + if desc.segmented_property_type not in types: + types.append(desc.segmented_property_type) + + return types + + def _get_src_uid_index(self, sop_instance_uid: str) -> int: + ind = np.argwhere(self._source_sop_instance_uids == sop_instance_uid) + if len(ind) == 0: + raise KeyError( + f'No such source frame: {sop_instance_uid}' + ) + return ind.item() + + def _get_pixels_by_seg_frame( + self, + seg_frames_matrix: np.ndarray, + segment_numbers: np.ndarray, + combine_segments: bool = False, + relabel: bool = False, + rescale_fractional: bool = True + ) -> np.ndarray: + """Construct a segmentation array given an array of frame numbers. + + The output array is either 4D (combine_segments=False) or 3D + (combine_segments=True), where dimensions are frames x rows x columes x + segments. + + Parameters + ---------- + seg_frames_matrix: np.ndarray + Two dimensional numpy array containing integers. Each row of the + array corresponds to a frame (0th dimension) of the output array. + Each column of the array corresponds to a segment of the output + array. Elements specify the (1-based) frame numbers of the + segmentation image that contain the pixel data for that output + frame and segment. A value of -1 signifies that there is no pixel + data for that frame/segment combination. + segment_numbers: np.ndarray + One dimensional numpy array containing segment numbers + corresponding to the columns of the seg frames matrix. + combine_segments: bool + If True, combine the different segments into a single label + map in which the value of a pixel represents its segment. + If False (the default), segments are binary and stacked down the + last dimension of the output array. + relabel: bool + If True and ``combine_segments`` is ``True``, the pixel values in + the output array are relabelled into the range ``0`` to + ``len(segment_numbers)`` (inclusive) accoring to the position of + the original segment numbers in ``segment_numbers`` parameter. If + ``combine_segments`` is ``False``, this has no effect. + rescale_fractional: bool + If this is a FRACTIONAL segmentation and ``rescale_fractional`` is + True, the raw integer-valued array stored in the segmentation image + output will be rescaled by the MaximumFractionalValue such that + each pixel lies in the range 0.0 to 1.0. If False, the raw integer + values are returned. If the segmentation has BINARY type, this + parameter has no effect. + + Returns + ------- + pixel_array: np.ndarray + Segmentation pixel array + + """ + # Checks on input array + if seg_frames_matrix.ndim != 2: + raise ValueError('Seg frames matrix must be a 2D array.') + if not issubclass(seg_frames_matrix.dtype.type, np.integer): + raise TypeError( + 'Seg frames matrix must have an integer data type.' + ) + + if seg_frames_matrix.min() < -1: + raise ValueError( + 'Seg frames matrix may not contain negative values other than ' + '-1.' + ) + if seg_frames_matrix.max() > self.NumberOfFrames: + raise ValueError( + 'Seg frames matrix contains values outside the range of ' + 'segmentation frames in the image.' + ) + + if segment_numbers.shape != (seg_frames_matrix.shape[1], ): + raise ValueError( + 'Segment numbers array does not match the shape of the ' + 'seg frames matrix.' + ) + + if ( + segment_numbers.min() < 1 or + segment_numbers.max() > self.number_of_segments + ): + raise ValueError( + 'Segment numbers array contains invalid values.' + ) + + # Initialize empty pixel array + pixel_array = np.zeros( + ( + seg_frames_matrix.shape[0], + self.Rows, + self.Columns, + seg_frames_matrix.shape[1] + ), + self.pixel_array.dtype + ) + + # Loop through output frames + for out_frm, frm_row in enumerate(seg_frames_matrix): + # Loop through segmentation frames + for out_seg, seg_frame_num in enumerate(frm_row): + if seg_frame_num >= 1: + seg_frame_ind = seg_frame_num.item() - 1 + # Copy data to to output array + if self.pixel_array.ndim == 2: + # Special case with a single segmentation frame + pixel_array[out_frm, :, :, out_seg] = \ + self.pixel_array + else: + pixel_array[out_frm, :, :, out_seg] = \ + self.pixel_array[seg_frame_ind, :, :] + + if rescale_fractional: + if self.segmentation_type == SegmentationTypeValues.FRACTIONAL: + if pixel_array.max() > self.MaximumFractionalValue: + raise RuntimeError( + 'Segmentation image contains values greater than the ' + 'MaximumFractionalValue recorded in the dataset.' + ) + max_val = self.MaximumFractionalValue + pixel_array = pixel_array.astype(np.float32) / max_val + + if combine_segments: + # Check whether segmentation is binary, or fractional with only + # binary values + if self.segmentation_type != SegmentationTypeValues.BINARY: + if not rescale_fractional: + raise ValueError( + 'In order to combine segments of a FRACTIONAL ' + 'segmentation image, rescale_fractional must be ' + 'set to True.' + ) + is_binary = np.isin( + np.unique(pixel_array), + np.array([0.0, 1.0]), + assume_unique=True + ) + if not is_binary.all(): + raise ValueError( + 'Cannot combine segments if the segmentation is not' + 'binary' + ) + + # Check for overlap by summing over the segments dimension + if np.any(pixel_array.sum(axis=-1) > 1): + raise RuntimeError( + 'Segments cannot be combined because they overlap' + ) + + # Scale the array by the segment number using broadcasting + if relabel: + pixel_value_map = np.arange(1, len(segment_numbers) + 1) + else: + pixel_value_map = segment_numbers + scaled_array = pixel_array * pixel_value_map.reshape(1, 1, 1, -1) + + # Combine segments by taking maximum down final dimension + max_array = scaled_array.max(axis=-1) + pixel_array = max_array + + return pixel_array + + def get_source_image_uids(self) -> List[Tuple[hd_UID, hd_UID, hd_UID]]: + """Get UIDs for all source SOP instances referenced in the dataset. + + Returns + ------- + List[Tuple[highdicom.UID, highdicom.UID, highdicom.UID]] + List of tuples containing Study Instance UID, Series Instance UID + and SOP Instance UID for every SOP Instance referenced in the + dataset. + + """ + return [self._ref_ins_lut[sop_uid] for sop_uid in self._ref_ins_lut] + + def get_default_dimension_index_pointers( + self + ) -> List[BaseTag]: + """Get the default list of tags used to index frames. + + The list of tags used to index dimensions depends upon how the + segmentation image was constructed, and is stored in the + DimensionIndexPointer attribute within the DimensionIndexSequence. The + list returned by this method matches the order of items in the + DimensionIndexSequence, but omits the ReferencedSegmentNumber + attribute, since this is handled differently to other tags when + indexing frames in highdicom. + + Returns + ------- + List[pydicom.tag.BaseTag] + List of tags used as the default dimension index pointers. + + """ + return self._dim_ind_pointers[:] + + def are_dimension_indices_unique( + self, + dimension_index_pointers: Sequence[Union[int, BaseTag]] + ) -> bool: + """Check if a list of index pointers uniquely identifies frames. + + For a given list of dimension index pointers, check whether every + combination of index values for these pointers identifies a unique + frame per segment in the segmentation image. This is a pre-requisite + for indexing using this list of dimension index pointers in the + :meth:`Segmentation.get_pixels_by_dimension_index_values()` method. + + Parameters + ---------- + dimension_index_pointers: Sequence[Union[int, pydicom.tag.BaseTag]] + Sequence of tags serving as dimension index pointers. + + Returns + ------- + bool + True if the specified list of dimension index pointers uniquely + identifies frames in the segmentation image. False otherwise. + + Raises + ------ + KeyError + If any of the elements of the ``dimension_index_pointers`` are not + valid dimension index pointers in this segmentation image. + + """ + for ptr in dimension_index_pointers: + if ptr not in self._dim_ind_pointers: + kw = keyword_for_tag(ptr) + if kw == '': + kw = '' + raise KeyError( + f'Tag {ptr} ({kw}) is not used as a dimension index ' + 'in this image.' + ) + # Create the sub-matrix of the look up table that indexes + # by the dimension index values + dim_ind_cols = [ + self._lut_dim_ind_cols[ptr] + for ptr in dimension_index_pointers + ] + lut = self._frame_lut[:, dim_ind_cols + [self._lut_seg_col]] + + return np.unique(lut, axis=0).shape[0] == lut.shape[0] + + def _check_indexing_with_source_frames( + self, + ignore_spatial_locations: bool = False + ) -> None: + """Check if indexing by source frames is possible. + + Raise exceptions with useful messages otherwise. + + Possible problems include: + * Spatial locations are not preserved. + * The dataset does not specify that spatial locations are preserved + and the user has not asserted that they are. + * At least one frame in the segmentation lists multiple + source frames. + + Parameters + ---------- + ignore_spatial_locations: bool + Allows the user to ignore whether spatial locations are preserved + in the frames. + + """ + # Checks that it is possible to index using source frames in this + # dataset + if self._locations_preserved is None: + if not ignore_spatial_locations: + raise RuntimeError( + 'Indexing via source frames is not permissible since this ' + 'image does not specify that spatial locations are ' + 'preserved in the course of deriving the segmentation ' + 'from the source image. If you are confident that spatial ' + 'locations are preserved, or do not require that spatial ' + 'locations are preserved, you may override this behavior ' + "with the 'ignore_spatial_locations' parameter." + ) + elif self._locations_preserved == SpatialLocationsPreservedValues.NO: + if not ignore_spatial_locations: + raise RuntimeError( + 'Indexing via source frames is not permissible since this ' + 'image specifies that spatial locations are not preserved ' + 'in the course of deriving the segmentation from the ' + 'source image. If you do not require that spatial ' + ' locations are preserved you may override this behavior ' + "with the 'ignore_spatial_locations' parameter." + ) + if not self._single_source_frame_per_seg_frame: + raise RuntimeError( + 'Indexing via source frames is not permissible since some ' + 'frames in the segmentation specify multiple source frames.' + ) + + def get_pixels_by_source_instance( + self, + source_sop_instance_uids: Sequence[str], + segment_numbers: Optional[Sequence[int]] = None, + combine_segments: bool = False, + relabel: bool = False, + ignore_spatial_locations: bool = False, + assert_missing_frames_are_empty: bool = False, + rescale_fractional: bool = True + ) -> np.ndarray: + """Get a pixel array for a list of source instances. + + This is intended for retrieving segmentation masks derived from + (series of) single frame source images. + + The output array will have 4 dimensions under the default behavior, and + 3 dimensions if ``combine_segments`` is set to ``True``. The first + dimension represents the source instances. ``pixel_array[i, ...]`` + represents the segmentation of ``source_sop_instance_uids[i]``. The + next two dimensions are the rows and columns of the frames, + respectively. + + When ``combine_segments`` is ``False`` (the default behavior), the + segments are stacked down the final (4th) dimension of the pixel array. + If ``segment_numbers`` was specified, then ``pixel_array[:, :, :, i]`` + represents the data for segment ``segment_numbers[i]``. If + ``segment_numbers`` was unspecified, then ``pixel_array[:, :, :, i]`` + represents the data for segment ``parser.segment_numbers[i]``. Note + that in neither case does ``pixel_array[:, :, :, i]`` represent + the segmentation data for the segment with segment number ``i``, since + segment numbers begin at 1 in DICOM. + + When ``combine_segments`` is ``True``, then the segmentation data from + all specified segments is combined into a multi-class array in which + pixel value is used to denote the segment to which a pixel belongs. + This is only possible if the segments do not overlap and either the + type of the segmentation is ``BINARY`` or the type of the segmentation + is ``FRACTIONAL`` but all values are exactly 0.0 or 1.0. the segments + do not overlap. If the segments do overlap, a ``RuntimeError`` will be + raised. After combining, the value of a pixel depends upon the + ``relabel`` parameter. In both cases, pixels that appear in no segments + with have a value of ``0``. If ``relabel`` is ``False``, a pixel that + appears in the segment with segment number ``i`` (according to the + original segment numbering of the segmentation object) will have a + value of ``i``. If ``relabel`` is ``True``, the value of a pixel in + segment ``i`` is related not to the original segment number, but to the + index of that segment number in the ``segment_numbers`` parameter of + this method. Specifically, pixels belonging to the segment with segment + number ``segment_numbers[i]`` is given the value ``i + 1`` in the + output pixel array (since 0 is reserved for pixels that belong to no + segments). In this case, the values in the output pixel array will + always lie in the range ``0`` to ``len(segment_numbers)`` inclusive. + + Parameters + ---------- + source_sop_instance_uids: str + SOP Instance UID of the source instances to for which segmentations + are requested. + segment_numbers: Union[Sequence[int], None], optional + Sequence containing segment numbers to include. If unspecified, + all segments are included. + combine_segments: bool, optional + If True, combine the different segments into a single label + map in which the value of a pixel represents its segment. + If False (the default), segments are binary and stacked down the + last dimension of the output array. + relabel: bool, optional + If True and ``combine_segments`` is ``True``, the pixel values in + the output array are relabelled into the range ``0`` to + ``len(segment_numbers)`` (inclusive) accoring to the position of + the original segment numbers in ``segment_numbers`` parameter. If + ``combine_segments`` is ``False``, this has no effect. + ignore_spatial_locations: bool, optional + Ignore whether or not spatial locations were preserved in the + derivation of the segmentation frames from the source frames. In + some segmentation images, the pixel locations in the segmentation + frames may not correspond to pixel locations in the frames of the + source image from which they were derived. The segmentation image + may or may not specify whether or not spatial locations are + preserved in this way through use of the optional (0028,135A) + SpatialLocationsPreserved attribute. If this attribute specifies + that spatial locations are not preserved, or is absent from the + segmentation image, highdicom's default behavior is to disallow + indexing by source frames. To override this behavior and retrieve + segmentation pixels regardless of the presence or value of the + spatial locations preserved attribute, set this parameter to True. + assert_missing_frames_are_empty: bool, optional + Assert that requested source frame numbers that are not referenced + by the segmentation image contain no segments. If a source frame + number is not referenced by the segmentation image, highdicom is + unable to check that the frame number is valid in the source image. + By default, highdicom will raise an error if any of the requested + source frames are not referenced in the source image. To override + this behavior and return a segmentation frame of all zeros for such + frames, set this parameter to True. + rescale_fractional: bool + If this is a FRACTIONAL segmentation and ``rescale_fractional`` is + True, the raw integer-valued array stored in the segmentation image + output will be rescaled by the MaximumFractionalValue such that + each pixel lies in the range 0.0 to 1.0. If False, the raw integer + values are returned. If the segmentation has BINARY type, this + parameter has no effect. + + Returns + ------- + pixel_array: np.ndarray + Pixel array representing the segmentation. See notes for full + explanation. + + Examples + -------- + + Read in an example from the highdicom test data: + + >>> import highdicom as hd + >>> + >>> seg = hd.seg.segread('data/test_files/seg_image_ct_binary.dcm') + + List the source images for this segmentation: + + >>> for study_uid, series_uid, sop_uid in seg.get_source_image_uids(): + >>> print(sop_uid) + 1.3.6.1.4.1.5962.1.1.0.0.0.1196530851.28319.0.93 + 1.3.6.1.4.1.5962.1.1.0.0.0.1196530851.28319.0.94 + 1.3.6.1.4.1.5962.1.1.0.0.0.1196530851.28319.0.95 + 1.3.6.1.4.1.5962.1.1.0.0.0.1196530851.28319.0.96 + + Get the segmentation array for a subset of these images: + + >>> pixels = seg.get_pixels_by_source_instance( + >>> source_sop_instance_uids=[ + >>> '1.3.6.1.4.1.5962.1.1.0.0.0.1196530851.28319.0.93', + >>> '1.3.6.1.4.1.5962.1.1.0.0.0.1196530851.28319.0.94' + >>> ] + >>> ) + >>> pixels.shape + (2, 16, 16, 1) + + """ + # Check that indexing in this way is possible + self._check_indexing_with_source_frames(ignore_spatial_locations) + + # Checks on validity of the inputs + if segment_numbers is None: + segment_numbers = list(self.segment_numbers) + if len(segment_numbers) == 0: + raise ValueError( + 'Segment numbers may not be empty.' + ) + if isinstance(source_sop_instance_uids, str): + raise TypeError( + 'source_sop_instance_uids should be a sequence of UIDs, not a ' + 'single UID' + ) + if len(source_sop_instance_uids) == 0: + raise ValueError( + 'Source SOP instance UIDs may not be empty.' + ) + + # Initialize seg frame numbers matrix with value signifying + # "empty" (-1) + rows = len(source_sop_instance_uids) + cols = len(segment_numbers) + seg_frames = -np.ones(shape=(rows, cols), dtype=np.int32) + + # Sub-matrix of the LUT used for indexing by instance and + # segment number + query_cols = [self._lut_src_instance_col, self._lut_seg_col] + lut = self._frame_lut[:, query_cols] + + # Check for uniqueness + if np.unique(lut, axis=0).shape[0] != lut.shape[0]: + raise RuntimeError( + 'Source SOP instance UIDs and segment numbers do not ' + 'uniquely identify frames of the segmentation image.' + ) + + # Build the segmentation frame matrix + for r, sop_uid in enumerate(source_sop_instance_uids): + # Check whether this source UID exists in the LUT + try: + src_uid_ind = self._get_src_uid_index(sop_uid) + except KeyError as e: + if assert_missing_frames_are_empty: + continue + else: + msg = ( + f'SOP Instance UID {sop_uid} does not match any ' + 'referenced source frames. To return an empty ' + 'segmentation mask in this situation, use the ' + "'assert_missing_frames_are_empty' parameter." + ) + raise KeyError(msg) from e + + # Iterate over segment numbers for this source instance + for c, seg_num in enumerate(segment_numbers): + # Use LUT to find the segmentation frame containing + # the source frame and segment number + qry = np.array([src_uid_ind, seg_num]) + seg_frm_indices = np.where(np.all(lut == qry, axis=1))[0] + if len(seg_frm_indices) == 1: + seg_frames[r, c] = seg_frm_indices[0] + 1 + elif len(seg_frm_indices) > 0: + # This should never happen due to earlier checks + # on unique rows + raise RuntimeError( + 'An unexpected error was encountered during ' + 'indexing of the segmentation image. Please ' + 'file a bug report with the highdicom repository.' + ) + # else no segmentation frame found for this segment number, + # assume this frame is empty and leave the entry in seg_frames + # as -1 + + return self._get_pixels_by_seg_frame( + seg_frames_matrix=seg_frames, + segment_numbers=np.array(segment_numbers), + combine_segments=combine_segments, + relabel=relabel, + rescale_fractional=rescale_fractional + ) + + def get_pixels_by_source_frame( + self, + source_sop_instance_uid: str, + source_frame_numbers: Sequence[int], + segment_numbers: Optional[Sequence[int]] = None, + combine_segments: bool = False, + relabel: bool = False, + ignore_spatial_locations: bool = False, + assert_missing_frames_are_empty: bool = False, + rescale_fractional: bool = True + ): + """Get a pixel array for a list of frames within a source instance. + + This is intended for retrieving segmentation masks derived from + multi-frame (enhanced) source images. All source frames for + which segmentations are requested must belong within the same + SOP Instance UID. + + The output array will have 4 dimensions under the default behavior, and + 3 dimensions if ``combine_segments`` is set to ``True``. The first + dimension represents the source frames. ``pixel_array[i, ...]`` + represents the segmentation of ``source_frame_numbers[i]``. The + next two dimensions are the rows and columns of the frames, + respectively. + + When ``combine_segments`` is ``False`` (the default behavior), the + segments are stacked down the final (4th) dimension of the pixel array. + If ``segment_numbers`` was specified, then ``pixel_array[:, :, :, i]`` + represents the data for segment ``segment_numbers[i]``. If + ``segment_numbers`` was unspecified, then ``pixel_array[:, :, :, i]`` + represents the data for segment ``parser.segment_numbers[i]``. Note + that in neither case does ``pixel_array[:, :, :, i]`` represent + the segmentation data for the segment with segment number ``i``, since + segment numbers begin at 1 in DICOM. + + When ``combine_segments`` is ``True``, then the segmentation data from + all specified segments is combined into a multi-class array in which + pixel value is used to denote the segment to which a pixel belongs. + This is only possible if the segments do not overlap and either the + type of the segmentation is ``BINARY`` or the type of the segmentation + is ``FRACTIONAL`` but all values are exactly 0.0 or 1.0. the segments + do not overlap. If the segments do overlap, a ``RuntimeError`` will be + raised. After combining, the value of a pixel depends upon the + ``relabel`` parameter. In both cases, pixels that appear in no segments + with have a value of ``0``. If ``relabel`` is ``False``, a pixel that + appears in the segment with segment number ``i`` (according to the + original segment numbering of the segmentation object) will have a + value of ``i``. If ``relabel`` is ``True``, the value of a pixel in + segment ``i`` is related not to the original segment number, but to the + index of that segment number in the ``segment_numbers`` parameter of + this method. Specifically, pixels belonging to the segment with segment + number ``segment_numbers[i]`` is given the value ``i + 1`` in the + output pixel array (since 0 is reserved for pixels that belong to no + segments). In this case, the values in the output pixel array will + always lie in the range ``0`` to ``len(segment_numbers)`` inclusive. + + Parameters + ---------- + source_sop_instance_uid: str + SOP Instance UID of the source instance that contains the source + frames. + source_frame_numbers: Sequence[int] + A sequence of frame numbers (1-based) within the source instance + for which segmentations are requested. + segment_numbers: Sequence[int, None], optional + Sequence containing segment numbers to include. If unspecified, + all segments are included. + combine_segments: bool, optional + If True, combine the different segments into a single label + map in which the value of a pixel represents its segment. + If False (the default), segments are binary and stacked down the + last dimension of the output array. + relabel: bool, optional + If True and ``combine_segments`` is ``True``, the pixel values in + the output array are relabelled into the range ``0`` to + ``len(segment_numbers)`` (inclusive) accoring to the position of + the original segment numbers in ``segment_numbers`` parameter. If + ``combine_segments`` is ``False``, this has no effect. + ignore_spatial_locations: bool, optional + Ignore whether or not spatial locations were preserved in the + derivation of the segmentation frames from the source frames. In + some segmentation images, the pixel locations in the segmentation + frames may not correspond to pixel locations in the frames of the + source image from which they were derived. The segmentation image + may or may not specify whether or not spatial locations are + preserved in this way through use of the optional (0028,135A) + SpatialLocationsPreserved attribute. If this attribute specifies + that spatial locations are not preserved, or is absent from the + segmentation image, highdicom's default behavior is to disallow + indexing by source frames. To override this behavior and retrieve + segmentation pixels regardless of the presence or value of the + spatial locations preserved attribute, set this parameter to True. + assert_missing_frames_are_empty: bool, optional + Assert that requested source frame numbers that are not referenced + by the segmentation image contain no segments. If a source frame + number is not referenced by the segmentation image, highdicom is + unable to check that the frame number is valid in the source image. + By default, highdicom will raise an error if any of the requested + source frames are not referenced in the source image. To override + this behavior and return a segmentation frame of all zeros for such + frames, set this parameter to True. + rescale_fractional: bool + If this is a FRACTIONAL segmentation and ``rescale_fractional`` is + True, the raw integer-valued array stored in the segmentation image + output will be rescaled by the MaximumFractionalValue such that + each pixel lies in the range 0.0 to 1.0. If False, the raw integer + values are returned. If the segmentation has BINARY type, this + parameter has no effect. + + Returns + ------- + pixel_array: np.ndarray + Pixel array representing the segmentation. See notes for full + explanation. + + Examples + -------- + + Read in an example from the highdicom test data derived from a + multiframe slide microscopy image: + + >>> import highdicom as hd + >>> + >>> seg = hd.seg.segread('data/test_files/seg_image_sm_control.dcm') + + List the source image SOP instance UID for this segmentation: + + >>> sop_uid = seg.get_source_image_uids()[0][2] + '1.2.826.0.1.3680043.9.7433.3.12857516184849951143044513877282227' + + Get the segmentation array for 3 of the frames in the multiframe source + image. The resulting segmentation array has 3 10 x 10 frames, one for + each source frame. The final dimension contains the 20 different + segments present in this segmentation. + + >>> pixels = seg.get_pixels_by_source_frame( + >>> source_sop_instance_uid=sop_uid, + >>> source_frame_numbers=[4, 5, 6] + >>> ) + >>> pixels.shape + (3, 10, 10, 20) + + This time, select only 4 of the 20 segments: + + >>> pixels = seg.get_pixels_by_source_frame( + >>> source_sop_instance_uid=sop_uid, + >>> source_frame_numbers=[4, 5, 6], + >>> segment_numbers=[10, 11, 12, 13] + >>> ) + >>> pixels.shape + (3, 10, 10, 4) + + Instead create a multiclass label map for each source frame. Note + that segments 6, 8, and 10 are present in the three chosen frames. + + >>> pixels = seg.get_pixels_by_source_frame( + >>> source_sop_instance_uid=sop_uid, + >>> source_frame_numbers=[4, 5, 6], + >>> combine_segments=True + >>> ) + >>> pixels.shape, np.unique(pixels) + (3, 10, 10), array([0, 6, 8, 10]) + + Now relabel the segments to give a pixel map with values between 0 + and 3 (inclusive): + + >>> pixels = seg.get_pixels_by_source_frame( + >>> source_sop_instance_uid=sop_uid, + >>> source_frame_numbers=[4, 5, 6], + >>> segment_numbers=[6, 8, 10] + >>> combine_segments=True, + >>> relabel=True + >>> ) + >>> pixels.shape, np.unique(pixels) + (3, 10, 10), array([0, 1, 2, 3]) + + """ + # Check that indexing in this way is possible + self._check_indexing_with_source_frames(ignore_spatial_locations) + + # Checks on validity of the inputs + if segment_numbers is None: + segment_numbers = list(self.segment_numbers) + if len(segment_numbers) == 0: + raise ValueError( + 'Segment numbers may not be empty.' + ) + + if len(source_frame_numbers) == 0: + raise ValueError( + 'Source frame numbers should not be empty.' + ) + if not all(f > 0 for f in source_frame_numbers): + raise ValueError( + 'Frame numbers are 1-based indices and must be > 0.' + ) + + rows = len(source_frame_numbers) + cols = len(segment_numbers) + seg_frames = -np.ones(shape=(rows, cols), dtype=np.int32) + + # Create the sub-matrix of the look up table that indexes + # by frame number and segment number within this + # instance + src_uid_ind = self._get_src_uid_index(source_sop_instance_uid) + col = self._lut_src_instance_col + lut_instance_mask = self._frame_lut[:, col] == src_uid_ind + lut = self._frame_lut[lut_instance_mask, :] + query_cols = [self._lut_src_frame_col, self._lut_seg_col] + lut = lut[:, query_cols] + + if np.unique(lut, axis=0).shape[0] != lut.shape[0]: + raise RuntimeError( + 'Source frame numbers and segment numbers do not ' + 'uniquely identify frames of the segmentation image.' + ) + + # Build the segmentation frame matrix + for r, src_frm_num in enumerate(source_frame_numbers): + # Check whether this source frame exists in the LUT + if src_frm_num not in lut[:, 0]: + if assert_missing_frames_are_empty: + continue + else: + msg = ( + f'Source frame number {src_frm_num} does not ' + 'match any referenced source frame. To return ' + 'an empty segmentation mask in this situation, ' + "use the 'assert_missing_frames_are_empty' " + 'parameter.' + ) + raise ValueError(msg) + + # Iterate over segment numbers for this source frame + for c, seg_num in enumerate(segment_numbers): + # Use LUT to find the segmentation frame containing + # the source frame and segment number + qry = np.array([src_frm_num, seg_num]) + seg_frm_indices = np.where(np.all(lut == qry, axis=1))[0] + if len(seg_frm_indices) == 1: + seg_frames[r, c] = seg_frm_indices[0] + 1 + elif len(seg_frm_indices) > 0: + # This should never happen due to earlier checks + # on unique rows + raise RuntimeError( + 'An unexpected error was encountered during ' + 'indexing of the segmentation image. Please ' + 'file a bug report with the highdicom repository.' + ) + # else no segmentation frame found for this segment number, + # assume this frame is empty and leave the entry in seg_frames + # as -1 + + return self._get_pixels_by_seg_frame( + seg_frames_matrix=seg_frames, + segment_numbers=np.array(segment_numbers), + combine_segments=combine_segments, + relabel=relabel, + rescale_fractional=rescale_fractional + ) + + def get_pixels_by_dimension_index_values( + self, + dimension_index_values: Sequence[Sequence[int]], + dimension_index_pointers: Optional[Sequence[int]] = None, + segment_numbers: Optional[Sequence[int]] = None, + combine_segments: bool = False, + relabel: bool = False, + assert_missing_frames_are_empty: bool = False, + rescale_fractional: bool = True + ): + """Get a pixel array for a list of dimension index values. + + This is intended for retrieving segmentation masks using the index + values within the segmentation object, without referring to the + source images from which the segmentation as derived. + + The output array will have 4 dimensions under the default behavior, and + 3 dimensions if ``combine_segments`` is set to ``True``. The first + dimension represents the source frames. ``pixel_array[i, ...]`` + represents the segmentation frame with index + ``dimension_index_values[i]``. The next two dimensions are the rows + and columns of the frames, respectively. + + When ``combine_segments`` is ``False`` (the default behavior), the + segments are stacked down the final (4th) dimension of the pixel array. + If ``segment_numbers`` was specified, then ``pixel_array[:, :, :, i]`` + represents the data for segment ``segment_numbers[i]``. If + ``segment_numbers`` was unspecified, then ``pixel_array[:, :, :, i]`` + represents the data for segment ``parser.segment_numbers[i]``. Note + that in neither case does ``pixel_array[:, :, :, i]`` represent + the segmentation data for the segment with segment number ``i``, since + segment numbers begin at 1 in DICOM. + + When ``combine_segments`` is ``True``, then the segmentation data from + all specified segments is combined into a multi-class array in which + pixel value is used to denote the segment to which a pixel belongs. + This is only possible if the segments do not overlap and either the + type of the segmentation is ``BINARY`` or the type of the segmentation + is ``FRACTIONAL`` but all values are exactly 0.0 or 1.0. the segments + do not overlap. If the segments do overlap, a ``RuntimeError`` will be + raised. After combining, the value of a pixel depends upon the + ``relabel`` parameter. In both cases, pixels that appear in no segments + with have a value of ``0``. If ``relabel`` is ``False``, a pixel that + appears in the segment with segment number ``i`` (according to the + original segment numbering of the segmentation object) will have a + value of ``i``. If ``relabel`` is ``True``, the value of a pixel in + segment ``i`` is related not to the original segment number, but to the + index of that segment number in the ``segment_numbers`` parameter of + this method. Specifically, pixels belonging to the segment with segment + number ``segment_numbers[i]`` is given the value ``i + 1`` in the + output pixel array (since 0 is reserved for pixels that belong to no + segments). In this case, the values in the output pixel array will + always lie in the range ``0`` to ``len(segment_numbers)`` inclusive. + + Parameters + ---------- + dimension_index_values: Sequence[Sequence[int]] + Dimension index values for the requested frames. Each element of + the sequence is a sequence of 1-based index values representing the + dimension index values for a single frame of the output + segmentation. The order of the index values within the inner + sequence is determined by the ``dimension_index_pointers`` + parameter, and as such the length of each inner sequence must + match the length of ``dimension_index_pointers`` parameter. + dimension_index_pointers: Union[Sequence[Union[int, pydicom.tag.BaseTag]], None], optional + The data element tags that identify the indices used in the + ``dimension_index_values`` parameter. Each element identifies a + data element tag by which frames are ordered in the segmentation + image dataset. If this parameter is set to ``None`` (the default), + the value of + :meth:`Segmentation.get_default_dimension_index_pointers()` is + used. Valid values of this parameter are are determined by + the construction of the segmentation image and include any + permutation of any subset of elements in the + :meth:`Segmentation.get_default_dimension_index_pointers()` list. + segment_numbers: Union[Sequence[int], None], optional + Sequence containing segment numbers to include. If unspecified, + all segments are included. + combine_segments: bool, optional + If True, combine the different segments into a single label + map in which the value of a pixel represents its segment. + If False (the default), segments are binary and stacked down the + last dimension of the output array. + relabel: bool, optional + If True and ``combine_segments`` is ``True``, the pixel values in + the output array are relabelled into the range ``0`` to + ``len(segment_numbers)`` (inclusive) accoring to the position of + the original segment numbers in ``segment_numbers`` parameter. If + ``combine_segments`` is ``False``, this has no effect. + assert_missing_frames_are_empty: bool, optional + Assert that requested source frame numbers that are not referenced + by the segmentation image contain no segments. If a source frame + number is not referenced by the segmentation image, highdicom is + unable to check that the frame number is valid in the source image. + By default, highdicom will raise an error if any of the requested + source frames are not referenced in the source image. To override + this behavior and return a segmentation frame of all zeros for such + frames, set this parameter to True. + rescale_fractional: bool + If this is a FRACTIONAL segmentation and ``rescale_fractional`` is + True, the raw integer-valued array stored in the segmentation image + output will be rescaled by the MaximumFractionalValue such that + each pixel lies in the range 0.0 to 1.0. If False, the raw integer + values are returned. If the segmentation has BINARY type, this + parameter has no effect. + + Returns + ------- + pixel_array: np.ndarray + Pixel array representing the segmentation. See notes for full + explanation. + + Examples + -------- + + Read a test image of a segmentation of a slide microscopy image + + >>> import highdicom as hd + >>> from pydicom.datadict import keyword_for_tag, tag_for_keyword + >>> from pydicom import dcmread + >>> + >>> ds = dcmread('data/test_files/seg_image_sm_control.dcm') + >>> seg = hd.seg.Segmentation.from_dataset(ds) + + Get the default list of dimension index values + + >>> for tag in seg.get_default_dimension_index_pointers(): + >>> print(keyword_for_tag(tag)) + ColumnPositionInTotalImagePixelMatrix + RowPositionInTotalImagePixelMatrix + XOffsetInSlideCoordinateSystem + YOffsetInSlideCoordinateSystem + ZOffsetInSlideCoordinateSystem + + + Use a subset of these index pointers to index the image + + >>> tags = [ + >>> tag_for_keyword('ColumnPositionInTotalImagePixelMatrix'), + >>> tag_for_keyword('RowPositionInTotalImagePixelMatrix') + >>> ] + >>> assert seg.are_dimension_indices_unique(tags) # True + + It is therefore possible to index using just this subset of + dimension indices + + >>> pixels = seg.get_pixels_by_dimension_index_values( + >>> dimension_index_pointers=tags, + >>> dimension_index_values=[[1, 1], [1, 2]] + >>> ) + >>> pixels.shape + (2, 10, 10, 20) + + """ # noqa: E501 + # Checks on validity of the inputs + if segment_numbers is None: + segment_numbers = list(self.segment_numbers) + if len(segment_numbers) == 0: + raise ValueError( + 'Segment numbers may not be empty.' + ) + + if dimension_index_pointers is None: + dimension_index_pointers = self._dim_ind_pointers + else: + for ptr in dimension_index_pointers: + if ptr not in self._dim_ind_pointers: + kw = keyword_for_tag(ptr) + if kw == '': + kw = '' + raise KeyError( + f'Tag {Tag(ptr)} ({kw}) is not used as a dimension ' + 'index in this image.' + ) + + if len(dimension_index_values) == 0: + raise ValueError( + 'Dimension index values should not be empty.' + ) + if len(dimension_index_pointers) == 0: + raise ValueError( + 'Dimension index pointers should not be empty.' + ) + rows = len(dimension_index_values) + cols = len(segment_numbers) + seg_frames = -np.ones(shape=(rows, cols), dtype=np.int32) + + # Create the sub-matrix of the look up table that indexes + # by the dimension index values + dim_ind_cols = [ + self._lut_dim_ind_cols[ptr] + for ptr in dimension_index_pointers + ] + lut = self._frame_lut[:, dim_ind_cols + [self._lut_seg_col]] + + if np.unique(lut, axis=0).shape[0] != lut.shape[0]: + raise RuntimeError( + 'The chosen dimension indices do not uniquely identify ' + 'frames of the segmentation image. You may need to provide ' + 'further indices to disambiguate.' + ) + + # Build the segmentation frame matrix + for r, ind_vals in enumerate(dimension_index_values): + if len(ind_vals) != len(dimension_index_pointers): + raise ValueError( + 'Number of provided indices does not match the expected ' + 'number.' + ) + if not all(v > 0 for v in ind_vals): + raise ValueError( + 'Indices are 1-based and must be greater than 1.' + ) + + # Check whether this frame exists in the LUT at all, ignoring the + # segment number column + qry = np.array(ind_vals) + seg_frm_indices = np.where(np.all(lut[:, :-1] == qry, axis=1))[0] + if len(seg_frm_indices) == 0: + if assert_missing_frames_are_empty: + # Nothing more to do + continue + else: + raise RuntimeError( + f'No frame with dimension index values {ind_vals} ' + 'found in the segmentation image. To return a frame of ' + 'zeros in this situation, set the ' + "'assert_missing_frames_are_empty' parameter to True." + ) + + # Iterate over requested segment numbers for this source frame + for c, seg_num in enumerate(segment_numbers): + # Use LUT to find the segmentation frame containing + # the index values and segment number + qry = np.array(list(ind_vals) + [seg_num]) + seg_frm_indices = np.where(np.all(lut == qry, axis=1))[0] + if len(seg_frm_indices) == 1: + seg_frames[r, c] = seg_frm_indices[0] + 1 + elif len(seg_frm_indices) > 0: + # This should never happen due to earlier checks + # on unique rows + raise RuntimeError( + 'An unexpected error was encountered during ' + 'indexing of the segmentation image. Please ' + 'file a bug report with the highdicom repository.' + ) + # else no segmentation frame found for this segment number, + # assume this frame is empty and leave the entry in seg_frames + # as -1 + + return self._get_pixels_by_seg_frame( + seg_frames_matrix=seg_frames, + segment_numbers=np.array(segment_numbers), + combine_segments=combine_segments, + relabel=relabel, + rescale_fractional=rescale_fractional + ) + + +def segread(fp: Union[str, bytes, PathLike, BinaryIO]) -> Segmentation: + """Read a segmentation image stored in DICOM File Format. + + Parameters + ---------- + fp: Union[str, bytes, os.PathLike] + Any file-like object representing a DICOM file containing a + Segmentation image. + + Returns + ------- + highdicom.seg.Segmentation + Segmentation image read from the file. + + """ + return Segmentation.from_dataset(dcmread(fp)) diff --git a/src/highdicom/sr/coding.py b/src/highdicom/sr/coding.py index a64a2378..97b7cd63 100644 --- a/src/highdicom/sr/coding.py +++ b/src/highdicom/sr/coding.py @@ -79,6 +79,46 @@ def __ne__(self, other: Any) -> bool: """ return not (self == other) + @classmethod + def from_dataset(cls, dataset: Dataset) -> 'CodedConcept': + """Construct a CodedConcept from an existing dataset. + + Parameters + ---------- + dataset: pydicom.dataset.Dataset + Dataset representing a coded concept. + + Returns + ------- + highdicom.sr.coding.CodedConcept: + Coded concept representation of the dataset. + + Raises + ------ + TypeError: + If the passed dataset is not a pydicom dataset. + AttributeError: + If the dataset does not contain the required elements for a + coded concept. + + """ + if not isinstance(dataset, Dataset): + raise TypeError( + 'Dataset must be a pydicom.dataset.Dataset.' + ) + for kw in ['CodeValue', 'CodeMeaning', 'CodingSchemeDesignator']: + if not hasattr(dataset, kw): + raise AttributeError( + 'Dataset does not contain the following attribute ' + f'required for coded concepts: {kw}.' + ) + return cls( + value=dataset.CodeValue, + scheme_designator=dataset.CodingSchemeDesignator, + meaning=dataset.CodeMeaning, + scheme_version=getattr(dataset, 'CodingSchemeVersion', None) + ) + @property def value(self) -> str: """str: value of either `CodeValue`, `LongCodeValue` or `URNCodeValue` diff --git a/src/highdicom/uid.py b/src/highdicom/uid.py index ccf4a08c..051e91f0 100644 --- a/src/highdicom/uid.py +++ b/src/highdicom/uid.py @@ -1,5 +1,5 @@ import logging -from typing import Type, TypeVar +from typing import Optional, Type, TypeVar import pydicom @@ -13,7 +13,8 @@ class UID(pydicom.uid.UID): """Unique DICOM identifier with a highdicom-specific UID prefix.""" - def __new__(cls: Type[T]) -> T: - prefix = '1.2.826.0.1.3680043.10.511.3.' - identifier = pydicom.uid.generate_uid(prefix=prefix) - return super().__new__(cls, identifier) + def __new__(cls: Type[T], val: Optional[str] = None) -> T: + if val is None: + prefix = '1.2.826.0.1.3680043.10.511.3.' + val = pydicom.uid.generate_uid(prefix=prefix) + return super().__new__(cls, val) diff --git a/tests/test_seg.py b/tests/test_seg.py index 5fbe1945..613097e2 100644 --- a/tests/test_seg.py +++ b/tests/test_seg.py @@ -6,6 +6,7 @@ import pytest from pydicom.data import get_testdata_file, get_testdata_files +from pydicom.datadict import tag_for_keyword from pydicom.filereader import dcmread from pydicom.sr.codedict import codes from pydicom.uid import ( @@ -23,16 +24,18 @@ ) from highdicom.enum import CoordinateSystemNames from highdicom.seg import ( + Segmentation, + segread, DimensionIndexSequence, SegmentDescription, -) -from highdicom.seg import ( + SegmentationTypeValues, SegmentAlgorithmTypeValues, SegmentsOverlapValues, - SegmentationTypeValues, + SegmentationFractionalTypeValues, ) -from highdicom.seg import Segmentation from highdicom.seg.utils import iter_segments +from highdicom.sr.coding import CodedConcept +from highdicom.uid import UID class TestAlgorithmIdentificationSequence(unittest.TestCase): @@ -59,8 +62,15 @@ def test_construction(self): assert item.AlgorithmFamilyCodeSequence[0] == self._family with pytest.raises(AttributeError): item.AlgorithmSource + with pytest.raises(AttributeError): item.AlgorithmParameters + assert seq.name == self._name + assert seq.version == self._version + assert seq.family == self._family + assert seq.source is None + assert seq.parameters is None + def test_construction_missing_required_argument(self): with pytest.raises(TypeError): AlgorithmIdentificationSequence( @@ -95,6 +105,8 @@ def test_construction_optional_argument(self): with pytest.raises(AttributeError): item.AlgorithmParameters + assert seq.source == self._source + def test_construction_optional_argument_2(self): seq = AlgorithmIdentificationSequence( name=self._name, @@ -109,9 +121,20 @@ def test_construction_optional_argument_2(self): for key, value in self._parameters.items() ]) assert item.AlgorithmParameters == parsed_params + assert seq.parameters == self._parameters with pytest.raises(AttributeError): item.AlgorithmSource + def test_malformed_params(self): + seq = AlgorithmIdentificationSequence( + self._name, + self._family, + self._version + ) + seq[0].AlgorithmParameters = 'some invalid parameters' + with pytest.raises(ValueError): + seq.parameters + class TestSegmentDescription(unittest.TestCase): @@ -124,7 +147,7 @@ def setUp(self): codes.SCT.MorphologicallyAbnormalStructure self._segmented_property_type = codes.SCT.Neoplasm self._segment_algorithm_type = \ - SegmentAlgorithmTypeValues.AUTOMATIC.value + SegmentAlgorithmTypeValues.AUTOMATIC self._algorithm_identification = AlgorithmIdentificationSequence( name='bla', family=codes.DCM.ArtificialIntelligence, @@ -150,7 +173,7 @@ def test_construction(self): self._segmented_property_category assert item.SegmentedPropertyTypeCodeSequence[0] == \ self._segmented_property_type - assert item.SegmentAlgorithmType == self._segment_algorithm_type + assert item.SegmentAlgorithmType == self._segment_algorithm_type.value assert item.SegmentAlgorithmName == \ self._algorithm_identification[0].AlgorithmName assert len(item.SegmentationAlgorithmIdentificationSequence) == 1 @@ -160,6 +183,26 @@ def test_construction(self): item.AnatomicRegionSequence item.PrimaryAnatomicStructureSequence + assert item.segment_number == self._segment_number + assert item.segment_label == self._segment_label + assert isinstance(item.segmented_property_category, CodedConcept) + property_category = item.segmented_property_category + assert property_category == self._segmented_property_category + assert isinstance(item.segmented_property_type, CodedConcept) + assert item.segmented_property_type == self._segmented_property_type + assert isinstance(item.algorithm_type, SegmentAlgorithmTypeValues) + algo_type = item.algorithm_type + assert algo_type == SegmentAlgorithmTypeValues( + self._segment_algorithm_type + ) + algo_id = item.algorithm_identification + assert isinstance(algo_id, AlgorithmIdentificationSequence) + + assert item.tracking_id is None + assert item.tracking_uid is None + assert len(item.anatomic_regions) == 0 + assert len(item.primary_anatomic_structures) == 0 + def test_construction_invalid_segment_number(self): with pytest.raises(ValueError): SegmentDescription( @@ -233,13 +276,14 @@ def test_construction_missing_required_argument_6(self): def test_construction_no_algo_id_manual_seg(self): # Omitting the algo id should not give an error if the segmentation # type is MANUAL - SegmentDescription( + item = SegmentDescription( segment_number=self._segment_number, segment_label=self._segment_label, segmented_property_category=self._segmented_property_category, segmented_property_type=self._segmented_property_type, algorithm_type=SegmentAlgorithmTypeValues.MANUAL ) + assert item.algorithm_identification is None def test_construction_optional_argument(self): item = SegmentDescription( @@ -254,8 +298,11 @@ def test_construction_optional_argument(self): ) assert item.TrackingID == self._tracking_id assert item.TrackingUID == self._tracking_uid + assert item.tracking_id == self._tracking_id + assert item.tracking_uid == self._tracking_uid with pytest.raises(AttributeError): item.AnatomicRegionSequence + with pytest.raises(AttributeError): item.PrimaryAnatomicStructureSequence def test_construction_optional_argument_2(self): @@ -271,13 +318,51 @@ def test_construction_optional_argument_2(self): ) assert len(item.AnatomicRegionSequence) == 1 assert item.AnatomicRegionSequence[0] == self._anatomic_region + assert len(item.anatomic_regions) == 1 + assert all( + isinstance(el, CodedConcept) for el in item.anatomic_regions + ) + assert item.anatomic_regions[0] == self._anatomic_region + assert len(item.PrimaryAnatomicStructureSequence) == 1 assert item.PrimaryAnatomicStructureSequence[0] == \ self._anatomic_structure + assert len(item.primary_anatomic_structures) == 1 + assert all( + isinstance(el, CodedConcept) + for el in item.primary_anatomic_structures + ) + assert item.primary_anatomic_structures[0] == self._anatomic_structure + with pytest.raises(AttributeError): item.TrackingID + with pytest.raises(AttributeError): item.TrackingUID + def test_construction_mismatched_ids(self): + with pytest.raises(TypeError): + SegmentDescription( + self._segment_number, + self._segment_label, + self._segmented_property_category, + self._segmented_property_type, + self._segment_algorithm_type, + self._algorithm_identification, + tracking_id=self._tracking_id, + ) + + def test_construction_mismatched_ids_2(self): + with pytest.raises(TypeError): + SegmentDescription( + self._segment_number, + self._segment_label, + self._segmented_property_category, + self._segmented_property_type, + self._segment_algorithm_type, + self._algorithm_identification, + tracking_uid=self._tracking_uid, + ) + class TestPixelMeasuresSequence(unittest.TestCase): @@ -1586,6 +1671,626 @@ def test_construction_optional_arguments_3(self): assert po_item.ImageOrientationSlide == list(image_orientation) +class TestSegmentationParsing(unittest.TestCase): + def setUp(self): + self._sm_control_seg_ds = dcmread( + 'data/test_files/seg_image_sm_control.dcm' + ) + self._sm_control_seg = Segmentation.from_dataset( + self._sm_control_seg_ds + ) + + self._ct_binary_seg_ds = dcmread( + 'data/test_files/seg_image_ct_binary.dcm' + ) + self._ct_binary_seg = Segmentation.from_dataset( + self._ct_binary_seg_ds + ) + + self._ct_binary_overlap_seg_ds = dcmread( + 'data/test_files/seg_image_ct_binary_overlap.dcm' + ) + self._ct_binary_overlap_seg = Segmentation.from_dataset( + self._ct_binary_overlap_seg_ds + ) + + self._ct_binary_fractional_seg_ds = dcmread( + 'data/test_files/seg_image_ct_binary_fractional.dcm' + ) + self._ct_binary_fractional_seg = Segmentation.from_dataset( + self._ct_binary_fractional_seg_ds + ) + + self._ct_true_fractional_seg_ds = dcmread( + 'data/test_files/seg_image_ct_true_fractional.dcm' + ) + self._ct_true_fractional_seg = Segmentation.from_dataset( + self._ct_true_fractional_seg_ds + ) + self._ct_segs = [ + self._ct_binary_seg, + self._ct_binary_fractional_seg, + self._ct_true_fractional_seg + ] + + def test_from_dataset(self): + assert isinstance(self._sm_control_seg, Segmentation) + + def test_segread(self): + seg = segread('data/test_files/seg_image_ct_true_fractional.dcm') + assert isinstance(seg, Segmentation) + seg = segread('data/test_files/seg_image_ct_binary_overlap.dcm') + assert isinstance(seg, Segmentation) + seg = segread('data/test_files/seg_image_sm_numbers.dcm') + assert isinstance(seg, Segmentation) + + def test_properties(self): + # SM segs + seg_type = self._sm_control_seg.segmentation_type + assert seg_type == SegmentationTypeValues.BINARY + assert self._sm_control_seg.segmentation_fractional_type is None + assert self._sm_control_seg.number_of_segments == 20 + assert self._sm_control_seg.segment_numbers == range(1, 21) + + assert len(self._sm_control_seg.segmented_property_categories) == 1 + seg_category = self._sm_control_seg.segmented_property_categories[0] + assert seg_category == codes.SCT.Tissue + seg_property = self._sm_control_seg.segmented_property_types[0] + assert seg_property == codes.SCT.ConnectiveTissue + + # CT segs + for seg in self._ct_segs: + seg_type = seg.segmentation_type + assert seg.number_of_segments == 1 + assert seg.segment_numbers == range(1, 2) + + assert len(seg.segmented_property_categories) == 1 + seg_category = seg.segmented_property_categories[0] + assert seg_category == codes.SCT.Tissue + seg_property = seg.segmented_property_types[0] + assert seg_property == codes.SCT.Bone + + seg_type = self._ct_binary_seg.segmentation_type + assert seg_type == SegmentationTypeValues.BINARY + seg_type = self._ct_binary_fractional_seg.segmentation_type + assert seg_type == SegmentationTypeValues.FRACTIONAL + seg_type = self._ct_true_fractional_seg.segmentation_type + assert seg_type == SegmentationTypeValues.FRACTIONAL + + frac_type = self._ct_binary_fractional_seg.segmentation_fractional_type + assert frac_type == SegmentationFractionalTypeValues.PROBABILITY + frac_type = self._ct_true_fractional_seg.segmentation_fractional_type + assert frac_type == SegmentationFractionalTypeValues.PROBABILITY + + def test_get_source_image_uids(self): + uids = self._sm_control_seg.get_source_image_uids() + assert len(uids) == 1 + ins_uids = uids[0] + assert len(ins_uids) == 3 + assert all(isinstance(uid, UID) for uid in ins_uids) + + def test_get_segment_description(self): + desc1 = self._sm_control_seg.get_segment_description(1) + desc20 = self._sm_control_seg.get_segment_description(20) + assert isinstance(desc1, SegmentDescription) + assert desc1.segment_number == 1 + assert isinstance(desc20, SegmentDescription) + assert desc20.segment_number == 20 + + def test_get_segment_numbers_no_filters(self): + seg_nums = self._sm_control_seg.get_segment_numbers() + assert seg_nums == list(self._sm_control_seg.segment_numbers) + + def test_get_segment_numbers_with_filters(self): + desc1 = self._sm_control_seg.get_segment_description(1) + + seg_nums = self._sm_control_seg.get_segment_numbers( + tracking_id=desc1.tracking_id + ) + assert seg_nums == [1] + + seg_nums = self._sm_control_seg.get_segment_numbers( + tracking_uid=desc1.tracking_uid + ) + assert seg_nums == [1] + + # All segments match these filters + seg_nums = self._sm_control_seg.get_segment_numbers( + segmented_property_category=codes.SCT.Tissue, + segmented_property_type=codes.SCT.ConnectiveTissue, + algorithm_type=SegmentAlgorithmTypeValues.AUTOMATIC + ) + assert seg_nums == list(self._sm_control_seg.segment_numbers) + + def test_get_tracking_ids(self): + desc1 = self._sm_control_seg.get_segment_description(1) + + tracking_id_tuples = self._sm_control_seg.get_tracking_ids() + n_segs = self._sm_control_seg.number_of_segments + assert len(tracking_id_tuples) == n_segs + ids, uids = zip(*tracking_id_tuples) + assert desc1.tracking_id in ids + assert desc1.tracking_uid in uids + + def test_get_tracking_ids_with_filters(self): + desc1 = self._sm_control_seg.get_segment_description(1) + + # All segments in this test image match these filters + tracking_id_tuples = self._sm_control_seg.get_tracking_ids( + segmented_property_category=codes.SCT.Tissue, + segmented_property_type=codes.SCT.ConnectiveTissue, + algorithm_type=SegmentAlgorithmTypeValues.AUTOMATIC + ) + n_segs = self._sm_control_seg.number_of_segments + assert len(tracking_id_tuples) == n_segs + ids, uids = zip(*tracking_id_tuples) + assert desc1.tracking_id in ids + assert desc1.tracking_uid in uids + + def test_get_tracking_ids_with_filters_2(self): + # No segments in this test image match these filters + tracking_id_tuples = self._sm_control_seg.get_tracking_ids( + segmented_property_category=codes.SCT.Tissue, + segmented_property_type=codes.SCT.Lung, + ) + assert len(tracking_id_tuples) == 0 + + def test_get_pixels_by_source_frames(self): + source_sop_uid = self._sm_control_seg.get_source_image_uids()[0][-1] + + source_frames_valid = [1, 2, 4, 5] + pixels = self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_valid + ) + + out_shape = ( + len(source_frames_valid), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + self._sm_control_seg.number_of_segments + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_invalid_source_frames(self): + source_sop_uid = self._sm_control_seg.get_source_image_uids()[0][-1] + + # (frame 3 has no segment) + source_frames_invalid = [1, 3, 4, 5] + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_invalid + ) + + def test_get_pixels_by_invalid_source_frames_with_assert(self): + source_sop_uid = self._sm_control_seg.get_source_image_uids()[0][-1] + + # (frame 3 has no segment) + source_frames_invalid = [1, 3, 4, 5] + pixels = self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_invalid, + assert_missing_frames_are_empty=True + ) + + out_shape = ( + len(source_frames_invalid), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + self._sm_control_seg.number_of_segments + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_source_frames_with_segments(self): + source_sop_uid = self._sm_control_seg.get_source_image_uids()[0][-1] + + source_frames_valid = [1, 2, 4, 5] + segments_valid = [1, 20] + pixels = self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_valid, + segment_numbers=segments_valid + ) + + out_shape = ( + len(source_frames_valid), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + len(segments_valid) + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_source_frames_with_invalid_segments(self): + source_sop_uid = self._sm_control_seg.get_source_image_uids()[0][-1] + + source_frames_valid = [1, 2, 4, 5] + segments_invalid = [1, 21] # 21 > 20 + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_valid, + segment_numbers=segments_invalid + ) + + def test_get_pixels_by_source_frames_combine(self): + source_sop_uid = self._sm_control_seg.get_source_image_uids()[0][-1] + + source_frames_valid = [1, 2, 4, 5] + # These segments match the above frames for this test image + segments_valid = [6, 7, 8, 9] + pixels = self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_valid, + segment_numbers=segments_valid, + combine_segments=True + ) + + out_shape = ( + len(source_frames_valid), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns + ) + assert pixels.shape == out_shape + assert np.all(np.unique(pixels) == np.array([0] + segments_valid)) + + pixels = self._sm_control_seg.get_pixels_by_source_frame( + source_sop_instance_uid=source_sop_uid, + source_frame_numbers=source_frames_valid, + segment_numbers=segments_valid, + combine_segments=True, + relabel=True + ) + assert pixels.shape == out_shape + assert np.all(np.unique(pixels) == np.arange(len(segments_valid) + 1)) + + def test_get_default_dimension_index_pointers(self): + ptrs = self._sm_control_seg.get_default_dimension_index_pointers() + assert len(ptrs) == 5 + + def test_are_dimension_indices_unique(self): + ptrs = self._sm_control_seg.get_default_dimension_index_pointers() + assert self._sm_control_seg.are_dimension_indices_unique(ptrs) + + ptr_kws = [ + 'ColumnPositionInTotalImagePixelMatrix', + 'RowPositionInTotalImagePixelMatrix' + ] + ptrs = [tag_for_keyword(kw) for kw in ptr_kws] + assert self._sm_control_seg.are_dimension_indices_unique(ptrs) + + ptr_kws = [ + 'XOffsetInSlideCoordinateSystem', + 'YOffsetInSlideCoordinateSystem' + ] + ptrs = [tag_for_keyword(kw) for kw in ptr_kws] + assert self._sm_control_seg.are_dimension_indices_unique(ptrs) + + ptr_kws = [ + 'ZOffsetInSlideCoordinateSystem' + ] + ptrs = [tag_for_keyword(kw) for kw in ptr_kws] + assert not self._sm_control_seg.are_dimension_indices_unique(ptrs) + + def test_are_dimension_indices_unique_invalid_ptrs(self): + ptr_kws = [ + 'ImagePositionPatient' + ] + ptrs = [tag_for_keyword(kw) for kw in ptr_kws] + with pytest.raises(KeyError): + self._sm_control_seg.are_dimension_indices_unique(ptrs) + + def test_get_pixels_by_dimension_index_values(self): + ind_values = [ + (1, 1, 5, 5, 1), + (2, 1, 4, 5, 1), + (3, 1, 3, 5, 1) + ] + pixels = self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + ) + + out_shape = ( + len(ind_values), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + self._sm_control_seg.number_of_segments + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_dimension_index_values_subset(self): + ptr_kws = [ + 'ColumnPositionInTotalImagePixelMatrix', + 'RowPositionInTotalImagePixelMatrix' + ] + ptrs = [tag_for_keyword(kw) for kw in ptr_kws] + + ind_values = [ + (1, 1), + (2, 1), + (3, 1) + ] + pixels = self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + dimension_index_pointers=ptrs + ) + + out_shape = ( + len(ind_values), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + self._sm_control_seg.number_of_segments + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_dimension_index_values_missing(self): + ind_values = [ + (1, 1, 4, 5, 1), + ] + with pytest.raises(RuntimeError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + ) + + pixels = self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + assert_missing_frames_are_empty=True + ) + + out_shape = ( + len(ind_values), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + self._sm_control_seg.number_of_segments + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_dimension_index_values_with_segments(self): + ind_values = [ + (1, 1, 5, 5, 1), + (2, 1, 4, 5, 1), + (3, 1, 3, 5, 1) + ] + segments = [1, 6, 11] + pixels = self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + segment_numbers=segments + ) + + out_shape = ( + len(ind_values), + self._sm_control_seg.Rows, + self._sm_control_seg.Columns, + len(segments) + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_dimension_index_values_invalid(self): + ind_values = [ + (1, 1, 5, 5, 1), + (2, 1, 4, 5, 1), + (3, 1, 3, 5, 1) + ] + ptrs = [tag_for_keyword('ImagePositionPatient')] + + # Invalid pointers + with pytest.raises(KeyError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + dimension_index_pointers=ptrs + ) + # Invalid values + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=[(-1, 1, 1, 1, 1)], + ) + # Empty values + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=[], + ) + # Empty pointers + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + dimension_index_pointers=[] + ) + # Empty segment numbers + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + segment_numbers=[] + ) + # Invalid segment numbers + with pytest.raises(ValueError): + self._sm_control_seg.get_pixels_by_dimension_index_values( + dimension_index_values=ind_values, + segment_numbers=[-1] + ) + + def test_get_pixels_by_source_instances(self): + all_source_sop_uids = [ + tup[-1] for tup in self._ct_binary_seg.get_source_image_uids() + ] + source_sop_uids = all_source_sop_uids[1:3] + + pixels = self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + ) + + out_shape = ( + len(source_sop_uids), + self._ct_binary_seg.Rows, + self._ct_binary_seg.Columns, + self._ct_binary_seg.number_of_segments + ) + assert pixels.shape == out_shape + + pixels = self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + combine_segments=True + ) + + out_shape = ( + len(source_sop_uids), + self._ct_binary_seg.Rows, + self._ct_binary_seg.Columns, + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_source_instances_with_segments(self): + all_source_sop_uids = [ + tup[-1] for tup in self._ct_binary_seg.get_source_image_uids() + ] + source_sop_uids = all_source_sop_uids[1:3] + segment_numbers = [1] + + pixels = self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + segment_numbers=segment_numbers + ) + + out_shape = ( + len(source_sop_uids), + self._ct_binary_seg.Rows, + self._ct_binary_seg.Columns, + len(segment_numbers) + ) + assert pixels.shape == out_shape + + def test_get_pixels_by_source_instances_invalid(self): + all_source_sop_uids = [ + tup[-1] for tup in self._ct_binary_seg.get_source_image_uids() + ] + source_sop_uids = all_source_sop_uids[1:3] + + # Empty SOP uids + with pytest.raises(ValueError): + self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=[], + ) + # Empty SOP uids + with pytest.raises(KeyError): + self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=['1.2.3.4'], + ) + # Empty segments + with pytest.raises(ValueError): + self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + segment_numbers=[] + ) + # Invalid segments + with pytest.raises(ValueError): + self._ct_binary_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + segment_numbers=[0] + ) + + def test_get_pixels_by_source_instances_binary_fractional(self): + all_source_sop_uids = [ + tup[-1] for tup in + self._ct_binary_fractional_seg.get_source_image_uids() + ] + source_sop_uids = all_source_sop_uids[1:3] + + pixels = self._ct_binary_fractional_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + ) + + out_shape = ( + len(source_sop_uids), + self._ct_binary_fractional_seg.Rows, + self._ct_binary_fractional_seg.Columns, + self._ct_binary_fractional_seg.number_of_segments + ) + assert pixels.shape == out_shape + assert np.all(np.unique(pixels) == np.array([0.0, 1.0])) + + pixels = self._ct_binary_fractional_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + combine_segments=True + ) + + out_shape = ( + len(source_sop_uids), + self._ct_binary_fractional_seg.Rows, + self._ct_binary_fractional_seg.Columns, + ) + assert pixels.shape == out_shape + assert np.all(np.unique(pixels) == np.array([0.0, 1.0])) + + def test_get_pixels_by_source_instances_true_fractional(self): + all_source_sop_uids = [ + tup[-1] for tup in + self._ct_true_fractional_seg.get_source_image_uids() + ] + source_sop_uids = all_source_sop_uids[1:3] + + pixels = self._ct_true_fractional_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + ) + + out_shape = ( + len(source_sop_uids), + self._ct_true_fractional_seg.Rows, + self._ct_true_fractional_seg.Columns, + self._ct_true_fractional_seg.number_of_segments + ) + assert pixels.shape == out_shape + assert pixels.max() <= 1.0 + assert pixels.min() >= 0.0 + assert len(np.unique(pixels)) > 2 + + # Without fractional rescaling + pixels = self._ct_true_fractional_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + rescale_fractional=False + ) + + out_shape = ( + len(source_sop_uids), + self._ct_true_fractional_seg.Rows, + self._ct_true_fractional_seg.Columns, + self._ct_true_fractional_seg.number_of_segments + ) + assert pixels.shape == out_shape + assert pixels.max() == 128 + assert len(np.unique(pixels)) > 2 + + # Can't combine segments with a true fractional segmentation + with pytest.raises(ValueError): + self._ct_true_fractional_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + combine_segments=True + ) + + def test_get_pixels_by_source_instances_overlap(self): + all_source_sop_uids = [ + tup[-1] for tup in + self._ct_binary_overlap_seg.get_source_image_uids() + ] + source_sop_uids = all_source_sop_uids + + pixels = self._ct_binary_overlap_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + ) + + out_shape = ( + len(source_sop_uids), + self._ct_binary_overlap_seg.Rows, + self._ct_binary_overlap_seg.Columns, + self._ct_binary_overlap_seg.number_of_segments + ) + assert pixels.shape == out_shape + + with pytest.raises(RuntimeError): + self._ct_binary_overlap_seg.get_pixels_by_source_instance( + source_sop_instance_uids=source_sop_uids, + combine_segments=True + ) + + class TestSegUtilities(unittest.TestCase): def setUp(self):