From 6d4fe128a78bc5b0b2758c7e3ae0b8ff5cb99576 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Thu, 15 Apr 2021 23:58:47 -0400 Subject: [PATCH 01/31] Initial attempt at seg parsing --- src/highdicom/seg/sop.py | 511 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 509 insertions(+), 2 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 6452b9cf..8efdbb43 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1,8 +1,9 @@ """Module for the SOP class of the Segmentation IOD.""" import logging import numpy as np +from operator import eq from collections import defaultdict -from typing import Any, Dict, List, Optional, Set, Sequence, Union, Tuple +from typing import Any, Dict, Iterable, List, Optional, Set, Sequence, Union, Tuple from pydicom.dataset import Dataset from pydicom.encaps import decode_data_sequence, encapsulate @@ -17,6 +18,7 @@ ) from pydicom.sr.codedict import codes from pydicom.valuerep import PersonName +from pydicom.sr.coding import Code from highdicom.base import SOPClass from highdicom.content import ( @@ -34,6 +36,7 @@ SegmentationFractionalTypeValues, SegmentationTypeValues, SegmentsOverlapValues, + SpatialLocationsPreservedValues ) from highdicom.sr.coding import CodedConcept from highdicom.valuerep import check_person_name @@ -79,7 +82,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 +539,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 +901,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 +950,500 @@ 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': + # 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.' + ) + dataset.__class__ = Segmentation + + dataset._source_images = None + dataset._source_plane_orientation = None + sf_groups = dataset.SharedFunctionalGroupsSequence[0] + plane_ori_seq = sf_groups.PlaneOrientationSequence[0] + if hasattr(plane_ori_seq, 'ImageOrientationSlide'): + dataset._coordinate_system = CoordinateSystemNames.SLIDE + dataset._plane_orientation = plane_ori_seq.ImageOrientationSlide + elif hasattr(plane_ori_seq, 'ImageOrientationPatient'): + dataset._coordinate_system = CoordinateSystemNames.PATIENT + dataset._plane_orientation = plane_ori_seq.ImageOrientationPatient + else: + raise ValueError( + 'Expected Plane Orientation Sequence to have either ' + 'ImageOrientationSlide or ImageOrientationPatient ' + 'attribute.' + ) + + for i, segment in enumerate(dataset.SegmentSequence, 1): + if segment.SegmentNumber != i: + raise ValueError( + 'Segments are expected to start at 1 and be consecutive ' + 'integers.' + ) + dataset._segment_inventory = { + s.SegmentNumber for s in dataset.SegmentSequence + } + + dataset._build_luts() + + return dataset + + def _build_luts(self) -> None: + + ref_series = self.ReferencedSeriesSequence[0] + self.source_sop_instance_uids = np.array([ + ref_ins.ReferencedSOPInstanceUID for ref_ins in + ref_series.ReferencedInstanceSequence + ]) + self.segment_description_lut = { + int(item.SegmentNumber): item + for item in self.SegmentSequence + } + + segnum_col_data = [] + source_frame_col_data = [] + + for frame_item in self.PerFrameFunctionalGroupsSequence: + # Get segment number for this frame + seg_no = frame_item.SegmentIdentificationSequence[0].\ + ReferencedSegmentNumber + segnum_col_data.append(int(seg_no)) + + source_images = [] + for der_im in frame_item.DerivationImageSequence: + for src_im in der_im.SourceImageSequence: + if 'SpatialLocationsPreserved' in src_im: + if eq( + src_im.SpatialLocationsPreserved, + SpatialLocationsPreservedValues.NO.value + ): + raise RuntimeError( + 'Segmentation dataset specifies that spatial ' + 'locations are not preserved' + ) + source_images.append(src_im.ReferencedSOPInstanceUID) + if len(set(source_images)) == 0: + raise RuntimeError( + 'Could not find source image information for ' + 'segmentation frames' + ) + if len(set(source_images)) > 1: + raise RuntimeError( + "Multiple source images found for segmentation frames" + ) + src_fm_ind = self._get_src_fm_index(source_images[0]) + source_frame_col_data.append(src_fm_ind) + + # Frame LUT is a 2D numpy array with two columns + # Row i represents frame i of the segmentation dataset + # The two columns represent the segment number and the source frame + # index in the source_sop_instance_uids list + # This allows for fairly efficient querying by any of the three + # variables: seg frame number, source frame number, segment number + self.frame_lut = np.array([segnum_col_data, source_frame_col_data]).T + self.lut_seg_col = 0 + self.lut_src_col = 1 + + # Build LUT from tracking ID -> segment numbers + self.tracking_id_lut: Dict[str, List[int]] = defaultdict(list) + for n, seg_info in self.segment_description_lut.items(): + if 'TrackingID' in seg_info: + self.tracking_id_lut[seg_info.TrackingID].append(n) + + # Build LUT from tracking UID -> segment numbers + self.tracking_uid_lut: Dict[str, int] = defaultdict(list) + for n, seg_info in self.segment_description_lut.items(): + if 'TrackingUID' in seg_info: + self.tracking_uid_lut[seg_info.TrackingUID].append(n) + + def is_binary(self) -> bool: + return ( + self.SegmentationType == + SegmentationTypeValues.BINARY.value + ) + + def iter_segments(self): + return iter_segments(self) + + @property + def number_of_segments(self) -> int: + return len(self.SegmentSequence) + + @property + def segment_numbers(self) -> List[int]: + # Do not assume segments are sorted (although they should be) + return sorted(list(self.segment_description_lut.keys())) + + def get_segment_description(self, segment_number: int) -> Dataset: + + if segment_number < 1: + raise ValueError(f'{segment_number} is an invalid segment number') + + if segment_number not in self.segment_description_lut: + raise KeyError( + f'No segment number {segment_number} found in dataset' + ) + return self.segment_description_lut[segment_number] + + @classmethod + def _coded_concept_sequence_to_code( + cls, + coded_concept: Dataset + ) -> Code: + scheme_version = ( + coded_concept.CodingSchemeVersion + if 'CodingSchemeVersion' in coded_concept else None + ) + return Code( + value=coded_concept.CodeValue, + meaning=coded_concept.CodeMeaning, + scheme_designator=coded_concept.CodingSchemeDesignator, + scheme_version=scheme_version + ) + + def get_segmented_property_category_code( + self, + segment_number: int + ) -> Code: + seg_desc = self.segment_description_lut[segment_number] + return self._coded_concept_sequence_to_code( + seg_desc.SegmentedPropertyCategoryCodeSequence[0] + ) + + def get_segmented_property_type_code(self, segment_number: int) -> Code: + seg_desc = self.segment_description_lut[segment_number] + return self._coded_concept_sequence_to_code( + seg_desc.SegmentedPropertyTypeCodeSequence[0] + ) + + def get_segment_number_property_mapping( + self + ) -> Dict[int, Tuple[Code, Code]]: + return { + n: + ( + self.get_segmented_property_category_code(n), + self.get_segmented_property_type_code(n) + ) + for n in self.segment_numbers + } + + def get_segment_number_for_segmented_property( + self, + segmented_property_category: Union[Code, CodedConcept], + segmented_property_type: Union[Code, CodedConcept], + ) -> int: + if isinstance(segmented_property_category, CodedConcept): + segmented_property_category = self._coded_concept_sequence_to_code( + segmented_property_category + ) + if isinstance(segmented_property_type, CodedConcept): + segmented_property_type = self._coded_concept_sequence_to_code( + segmented_property_type + ) + + segment_number: Optional[int] = None + + for n in self.segment_numbers: + category_match = ( + segmented_property_category == + self.get_segmented_property_category_code(n) + ) + type_match = ( + segmented_property_type == + self.get_segmented_property_type_code(n) + ) + + if category_match and type_match: + if segment_number is not None: + raise KeyError( + 'Multiple segments found with specified segment ' + 'property' + ) + segment_number = n + + if segment_number is None: + raise KeyError( + 'No segment found with specified segment property' + ) + + return segment_number + + def get_segment_tracking_id(self, segment_number: int) -> Optional[str]: + seg_info = self.segment_description_lut[segment_number] + if 'TrackingID' in seg_info: + return self.segment_description_lut[segment_number].TrackingID + return None + + def get_segment_tracking_uid(self, segment_number: int) -> Optional[str]: + seg_info = self.segment_description_lut[segment_number] + if 'TrackingUID' in seg_info: + return self.segment_description_lut[segment_number].TrackingUID + return None + + def get_segment_numbers_for_tracking_id( + self, + tracking_id: str + ) -> List[int]: + if tracking_id not in self.tracking_id_lut: + raise KeyError( + 'No segment found matching specified tracking ID' + ) + return self.tracking_id_lut[tracking_id] + + def get_segment_numbers_for_tracking_uid( + self, + tracking_uid: str + ) -> List[int]: + if tracking_uid not in self.tracking_uid_lut: + raise KeyError( + 'No segment found matching specified tracking UID' + ) + return self.tracking_uid_lut[tracking_uid] + + def get_tracking_ids(self) -> List[str]: + return list(self.tracking_id_lut.keys()) + + def get_tracking_uids(self) -> List[str]: + return list(self.tracking_uid_lut.keys()) + + def _get_src_fm_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 list_seg_frames_for_source_frame( + self, + source_sop_instance_uid: str + ) -> List[int]: + src_ind = self._get_src_fm_index(source_sop_instance_uid) + seg_frames = np.where( + self.frame_lut[:, self.lut_src_col] == src_ind + )[0] + return seg_frames.tolist() + + def list_segments_in_source_frame( + self, + source_sop_instance_uid: str + ) -> List[int]: + src_ind = self._get_src_fm_index(source_sop_instance_uid) + segments = self.frame_lut[ + self.frame_lut[:, self.lut_src_col] == src_ind, + self.lut_seg_col + ] + return segments.tolist() + + def list_seg_frames_for_segment_number( + self, + segment_number: int + ) -> List[int]: + if segment_number not in self.segment_numbers: + raise KeyError( + 'No segment found with specified segment number' + ) + seg_frames = np.where( + self.frame_lut[:, self.lut_seg_col] == segment_number + )[0].tolist() + return seg_frames + + def list_source_frames_for_segment_number( + self, + segment_number: int + ) -> List[str]: + if segment_number not in self.segment_numbers: + raise KeyError( + 'No segment found with specified segment number' + ) + source_frame_indices = self.frame_lut[ + self.frame_lut[:, self.lut_seg_col] == segment_number, + self.lut_src_col + ] + source_frame_uids = self.source_sop_instance_uids[source_frame_indices] + return source_frame_uids.tolist() + + def get_pixels( + self, + source_frames: Optional[Iterable[str]] = None, + segment_numbers: Optional[Iterable[int]] = None, + combine_segments: bool = False, + remap_segment_numbers: bool = False, + ) -> np.ndarray: + """Get a numpy array of the reconstructed segmentation. + Parameters + ---------- + source_frames: Optional[Iterable[str]] + Iterable containing SOP Instance UIDs of the source images to + include in the segmentation. If unspecified, all source images + are included. + segment_numbers: Optional[Sequence[int]] + Sequence containing segment numbers to include. If unspecified, + all segments are included. + 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. + remap_segment_numbers: bool + If True and ``combine_segments`` is ``True``, the pixel values in + the output array are remapped 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. + Returns + ------- + pixel_array: np.ndarray + Pixel array representing the segmentation + Note + ---- + 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. If ``source_frames`` + has been specified, then ``pixel_array[i, ...]`` represents the + segmentation of ``source_frames[i]``. If ``source_frames`` was not + specified, then ``pixel_array[i, ...]`` represents the segmentation of + ``parser.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 type of the segmentation is ``BINARY`` and + 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 ``remap_segment_numbers`` parameter. In both + cases, pixels that appear in no segments with have a value of ``0``. + If ``remap_segment_numbers`` 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 ``remap_segment_numbers`` 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. + """ + if combine_segments and not self.is_binary(): + raise ValueError( + 'Cannot combine segments if the segmentation is not binary' + ) + + # If no source frames were specified, use all source frames + if source_frames is None: + source_frames = self.source_sop_instance_uids + + # If no segments were specified, use all segments + if segment_numbers is None: + segment_numbers = self.segment_numbers + + # Check all provided segmentation numbers are valid + seg_num_arr = np.array(segment_numbers) + seg_nums_exist = np.in1d( + seg_num_arr, np.array(self.segment_numbers) + ) + if not seg_nums_exist.all(): + missing_seg_nums = seg_num_arr[np.logical_not(seg_nums_exist)] + raise KeyError( + 'Segment numbers contained non-existent segments: ' + f'{missing_seg_nums}' + ) + + # Check segment numbers are unique + if len(np.unique(seg_num_arr)) != len(seg_num_arr): + raise ValueError( + 'Segment numbers must not contain duplicates' + ) + + # Initialize empty pixel array + pixel_array = np.zeros( + ( + len(source_frames), + self.Rows, + self.Columns, + len(segment_numbers) + ), + self.pixel_array.dtype + ) + + # Loop through source frames + for out_ind, src_frame in enumerate(source_frames): + src_ind = self._get_src_fm_index(src_frame) + + # Find segmentation frames relating to this source frame + seg_frames = np.where( + self.frame_lut[:, self.lut_src_col] == src_ind + )[0] + seg_nums_for_frame = self.frame_lut[seg_frames, self.lut_seg_col] + + # Loop through segmentation frames + for seg_frame, seg_num in zip(seg_frames, seg_nums_for_frame): + # Find output index for this segmentation number (if any) + seg_ind = np.where(seg_num_arr == seg_num)[0] + + if len(seg_ind) > 0: + # Copy data to to output array + if self.pixel_array.ndim == 2: + # Special case with a single segmentation frame + pixel_array[out_ind, :, :, seg_ind.item()] = \ + self.pixel_array + else: + pixel_array[out_ind, :, :, seg_ind.item()] = \ + self.pixel_array[seg_frame, :, :] + + if combine_segments: + # 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 remap_segment_numbers: + pixel_value_map = np.arange(1, len(segment_numbers) + 1) + else: + pixel_value_map = seg_num_arr + 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 + + # TODO + # Shift segment information getters to SegmentDescription + # Multi-frame source image + # Segs with a single frame + # Allow binary fractional segs to combine segments + # Option to force standards compliance + # Integrity checks: + # Each source frame specified in pffgs exists + # No combination of source frame and segment number has + # multiple source frames + # Tracking IDs do not have multiple segments with different + # segmented properties + # Lazy decoding of segmentation frames + # Should functions raise error or return empty lists? + # Standardize method names + # List segmented property for tracking id and vice versa + # Correct for spatial ordering of source frames? + # Optimize combine_segments to exploit sparse nature of array + # Ensure output array cannot be used to change value of pixel_array + From f8ea54ae639e6569f4f0cfad1759957b74d537b0 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 9 May 2021 18:34:55 -0400 Subject: [PATCH 02/31] Searching for segments by property can return multiple segments --- src/highdicom/seg/sop.py | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 8efdbb43..41d10003 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1131,11 +1131,19 @@ def get_segment_number_property_mapping( for n in self.segment_numbers } - def get_segment_number_for_segmented_property( + def get_segment_numbers_for_segmented_property( self, - segmented_property_category: Union[Code, CodedConcept], - segmented_property_type: Union[Code, CodedConcept], - ) -> int: + segmented_property_category: Union[Code, CodedConcept, None] = None, + segmented_property_type: Union[Code, CodedConcept, None] = None, + ) -> List[int]: + if ( + segmented_property_category is None and + segmented_property_type is None + ): + raise TypeError( + "At least one of the segmented_property_category and " + "segmented_property_type must be provided." + ) if isinstance(segmented_property_category, CodedConcept): segmented_property_category = self._coded_concept_sequence_to_code( segmented_property_category @@ -1145,7 +1153,7 @@ def get_segment_number_for_segmented_property( segmented_property_type ) - segment_number: Optional[int] = None + matched_segment_numbers: List[int] = [] for n in self.segment_numbers: category_match = ( @@ -1158,19 +1166,9 @@ def get_segment_number_for_segmented_property( ) if category_match and type_match: - if segment_number is not None: - raise KeyError( - 'Multiple segments found with specified segment ' - 'property' - ) - segment_number = n - - if segment_number is None: - raise KeyError( - 'No segment found with specified segment property' - ) + matched_segment_numbers.append(n) - return segment_number + return matched_segment_numbers def get_segment_tracking_id(self, segment_number: int) -> Optional[str]: seg_info = self.segment_description_lut[segment_number] From d861a6252ff0d239dc910714bf2132175bad941a Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 9 May 2021 18:50:16 -0400 Subject: [PATCH 03/31] Rename parameter to 'relabel' --- src/highdicom/seg/sop.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 41d10003..c9f7af11 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1270,9 +1270,10 @@ def get_pixels( source_frames: Optional[Iterable[str]] = None, segment_numbers: Optional[Iterable[int]] = None, combine_segments: bool = False, - remap_segment_numbers: bool = False, + relabel: bool = False, ) -> np.ndarray: """Get a numpy array of the reconstructed segmentation. + Parameters ---------- source_frames: Optional[Iterable[str]] @@ -1287,16 +1288,18 @@ def get_pixels( 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. - remap_segment_numbers: bool + relabel: bool If True and ``combine_segments`` is ``True``, the pixel values in - the output array are remapped into the range ``0`` to + 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. + Returns ------- pixel_array: np.ndarray Pixel array representing the segmentation + Note ---- The output array will have 4 dimensions under the default behavior, @@ -1323,12 +1326,12 @@ def get_pixels( This is only possible if the type of the segmentation is ``BINARY`` and 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 ``remap_segment_numbers`` parameter. In both + depends upon the ``relabel`` parameter. In both cases, pixels that appear in no segments with have a value of ``0``. - If ``remap_segment_numbers`` is ``False``, a pixel that appears + 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 ``remap_segment_numbers`` is ``True``, the 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 @@ -1413,7 +1416,7 @@ def get_pixels( ) # Scale the array by the segment number using broadcasting - if remap_segment_numbers: + if relabel: pixel_value_map = np.arange(1, len(segment_numbers) + 1) else: pixel_value_map = seg_num_arr From e515918232393d04e925bc6f2ac13b20a5b21f7c Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 16 May 2021 01:04:49 -0400 Subject: [PATCH 04/31] Added attribute checking logic and applied to segment description --- src/highdicom/module_utils.py | 168 ++++++++++++++++++++++++++++++++++ src/highdicom/seg/content.py | 80 ++++++++++++++++ 2 files changed, 248 insertions(+) create mode 100644 src/highdicom/module_utils.py diff --git a/src/highdicom/module_utils.py b/src/highdicom/module_utils.py new file mode 100644 index 00000000..e657bc8e --- /dev/null +++ b/src/highdicom/module_utils.py @@ -0,0 +1,168 @@ +from enum import Enum +from typing import Any, Dict, Iterator, List, Optional, Sequence, Tuple, Union + +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) + + for p in base_path: + try: + tree = tree['attributes'][p] + except: + 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 = {'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/seg/content.py b/src/highdicom/seg/content.py index 99f5fc06..9287aeb1 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -15,6 +15,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 +153,85 @@ 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.content.SegmentDescription + Segment description. + + """ + check_required_attributes( + dataset, + module='segmentation-image', + base_path=['SegmentSequence'] + ) + dataset.__class__ = SegmentDescription + return dataset + + @property + def segmented_property_category(self)-> CodedConcept: + code_seq = self.SegmentedPropertyCategoryCodeSequence[0] + scheme_version = getattr(code_seq, 'CodingSchemeVersion', None) + return CodedConcept( + value=code_seq.CodeValue, + scheme_designator=code_seq.CodingSchemeDesignator, + meaning=code_seq.CodeMeaning, + scheme_version=scheme_version + ) + + @property + def segmented_property_type(self)-> CodedConcept: + code_seq = self.SegmentedPropertyTypeCodeSequence[0] + scheme_version = getattr(code_seq, 'CodingSchemeVersion', None) + return CodedConcept( + value=code_seq.CodeValue, + scheme_designator=code_seq.CodingSchemeDesignator, + meaning=code_seq.CodeMeaning, + scheme_version=scheme_version + ) + + @property + def algorithm_type(self) -> SegmentAlgorithmTypeValues: + return SegmentAlgorithmTypeValues(self.SegmentAlgorithmType) + + @property + def tracking_id(self) -> Optional[str]: + if 'TrackingID' in self: + return self.TrackingID + return None + + @property + def tracking_uid(self) -> Optional[str]: + if 'TrackingUID' in self: + return self.TrackingUID + return None + + @property + def segment_number(self) -> int: + return self.SegmentNumber + + #@property + #def anatomic_regions(self) -> List[CodedConcept]: + # pass + + #@property + #def primary_anatomic_structures(self) -> List[CodedConcept]: + # pass + + #@property + #def algorithm_identification(self) -> Optional[ + # AlgorithmIdentificationSequence + # ]: + # pass + class DimensionIndexSequence(DataElementSequence): From 7f5a54544c481d5ac56eea60bc09365c94b23c8d Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 16 May 2021 16:47:34 -0400 Subject: [PATCH 05/31] Implemented conversion and parsing for algorithm id and segment description --- src/highdicom/content.py | 84 +++++++++++++++++++++- src/highdicom/seg/content.py | 130 ++++++++++++++++++++++++++--------- src/highdicom/seg/sop.py | 35 ++++++---- src/highdicom/sr/coding.py | 40 +++++++++++ 4 files changed, 241 insertions(+), 48 deletions(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index ce878302..c9bc88fc 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]: + """Optional[str]: + 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]]: + """Optional[Dict[str, str]]: + 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 AttributeError('Malformed parameter string') + parameters[split[0]] = split[1] + return parameters + class PixelMeasuresSequence(DataElementSequence): diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index 9287aeb1..fdeb0d7a 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 @@ -168,69 +169,130 @@ def from_dataset(cls, dataset: Dataset) -> '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'] ) - dataset.__class__ = SegmentDescription - return dataset + 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 self.SegmentNumber + + @property + def segment_label(self) -> str: + """str: Label of the segment.""" + return self.SegmentLabel @property def segmented_property_category(self)-> CodedConcept: - code_seq = self.SegmentedPropertyCategoryCodeSequence[0] - scheme_version = getattr(code_seq, 'CodingSchemeVersion', None) - return CodedConcept( - value=code_seq.CodeValue, - scheme_designator=code_seq.CodingSchemeDesignator, - meaning=code_seq.CodeMeaning, - scheme_version=scheme_version - ) + """highdicom.sr.coding.CodedConcept: + Category of the property the segment represents. + + """ + return self.SegmentedPropertyCategoryCodeSequence[0] @property def segmented_property_type(self)-> CodedConcept: - code_seq = self.SegmentedPropertyTypeCodeSequence[0] - scheme_version = getattr(code_seq, 'CodingSchemeVersion', None) - return CodedConcept( - value=code_seq.CodeValue, - scheme_designator=code_seq.CodingSchemeDesignator, - meaning=code_seq.CodeMeaning, - scheme_version=scheme_version - ) + """highdicom.sr.coding.CodedConcept: + Type of the property the segment represents. + + """ + return self.SegmentedPropertyTypeCodeSequence[0] @property def algorithm_type(self) -> SegmentAlgorithmTypeValues: + """highdicom.seg.enum.SegmentAlgorithmTypeValues: + Type of algorithm used to create the segment. + + """ return SegmentAlgorithmTypeValues(self.SegmentAlgorithmType) @property - def tracking_id(self) -> Optional[str]: - if 'TrackingID' in self: - return self.TrackingID + def algorithm_identification( + self + ) -> Optional[AlgorithmIdentificationSequence]: + """Optional[highdicom.content.AlgorithmIdentificationSequence] + Information useful for identification of the algorithm, if any. + + """ + if hasattr(self, 'SegmentationAlgorithmIdentificationSequence'): + return self.SegmentationAlgorithmIdentificationSequence return None @property def tracking_uid(self) -> Optional[str]: + """Optional[str]: + Tracking unique identifier for the segment, if any. + + """ if 'TrackingUID' in self: return self.TrackingUID return None @property - def segment_number(self) -> int: - return self.SegmentNumber + def tracking_id(self) -> Optional[str]: + """Optional[str]: Tracking identifier for the segment, if any.""" + if 'TrackingID' in self: + return self.TrackingID + return None - #@property - #def anatomic_regions(self) -> List[CodedConcept]: - # pass + @property + def anatomic_regions(self) -> List[CodedConcept]: + """List[highdicom.sr.coding.CodedConcept]: + List of anatomic regions into which the segment falls. + May be empty. - #@property - #def primary_anatomic_structures(self) -> List[CodedConcept]: - # pass + """ + if not hasattr(self, 'AnatomicRegionSequence'): + return [] + return self.AnatomicRegionSequence - #@property - #def algorithm_identification(self) -> Optional[ - # AlgorithmIdentificationSequence - # ]: - # pass + @property + def primary_anatomic_structures(self) -> List[CodedConcept]: + """List[highdicom.sr.coding.CodedConcept]: + List of anatomic anatomic structures the segment represents. + May be empty. + + """ + if not hasattr(self, 'PrimaryAnatomicStructureSequence'): + return [] + return self.PrimaryAnatomicStructureSequence class DimensionIndexSequence(DataElementSequence): diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index c9f7af11..a985ca00 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1,4 +1,5 @@ """Module for the SOP class of the Segmentation IOD.""" +from copy import deepcopy import logging import numpy as np from operator import eq @@ -953,23 +954,28 @@ def _encode_pixels(self, planes: np.ndarray) -> bytes: @classmethod def from_dataset(cls, dataset: Dataset) -> '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.' ) - dataset.__class__ = Segmentation + seg = deepcopy(dataset) + seg.__class__ = cls - dataset._source_images = None - dataset._source_plane_orientation = None - sf_groups = dataset.SharedFunctionalGroupsSequence[0] + 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'): - dataset._coordinate_system = CoordinateSystemNames.SLIDE - dataset._plane_orientation = plane_ori_seq.ImageOrientationSlide + seg._coordinate_system = CoordinateSystemNames.SLIDE + seg._plane_orientation = plane_ori_seq.ImageOrientationSlide elif hasattr(plane_ori_seq, 'ImageOrientationPatient'): - dataset._coordinate_system = CoordinateSystemNames.PATIENT - dataset._plane_orientation = plane_ori_seq.ImageOrientationPatient + seg._coordinate_system = CoordinateSystemNames.PATIENT + seg._plane_orientation = plane_ori_seq.ImageOrientationPatient else: raise ValueError( 'Expected Plane Orientation Sequence to have either ' @@ -977,19 +983,22 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': 'attribute.' ) - for i, segment in enumerate(dataset.SegmentSequence, 1): + for i, segment in enumerate(seg.SegmentSequence, 1): if segment.SegmentNumber != i: raise ValueError( 'Segments are expected to start at 1 and be consecutive ' 'integers.' ) - dataset._segment_inventory = { - s.SegmentNumber for s in dataset.SegmentSequence + seg._segment_inventory = { + s.SegmentNumber for s in seg.SegmentSequence } + seg.SegmentSequence = [ + SegmentDescription.from_dataset(ds) for ds in seg.SegmentSequence + ] - dataset._build_luts() + seg._build_luts() - return dataset + return seg def _build_luts(self) -> None: diff --git a/src/highdicom/sr/coding.py b/src/highdicom/sr/coding.py index a64a2378..75d402a2 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: {attr}.' + ) + 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` From bf85baef6ce7b856248d00bcdf386037a82c7d8c Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 16 May 2021 22:46:40 -0400 Subject: [PATCH 06/31] Further implementation of segmentation from_dataset conversion and parsing. --- src/highdicom/content.py | 136 ++++++++++++++++++++++ src/highdicom/module_utils.py | 4 +- src/highdicom/seg/content.py | 6 +- src/highdicom/seg/sop.py | 208 ++++++++++++++++------------------ src/highdicom/sr/coding.py | 2 +- 5 files changed, 237 insertions(+), 119 deletions(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index c9bc88fc..dcdde048 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -191,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 + ------- + 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): @@ -294,6 +339,53 @@ 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 + ------- + 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): @@ -374,6 +466,50 @@ 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 + ------- + 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/module_utils.py b/src/highdicom/module_utils.py index e657bc8e..60b061b2 100644 --- a/src/highdicom/module_utils.py +++ b/src/highdicom/module_utils.py @@ -1,5 +1,5 @@ from enum import Enum -from typing import Any, Dict, Iterator, List, Optional, Sequence, Tuple, Union +from typing import Any, Dict, List, Optional, Sequence from pydicom import Dataset @@ -84,7 +84,7 @@ def check_required_attributes( for p in base_path: try: tree = tree['attributes'][p] - except: + except KeyError: raise AttributeError(f"Invalid base path: {base_path}.") # Define recursive function diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index fdeb0d7a..e65206ca 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -220,7 +220,7 @@ def segment_label(self) -> str: return self.SegmentLabel @property - def segmented_property_category(self)-> CodedConcept: + def segmented_property_category(self) -> CodedConcept: """highdicom.sr.coding.CodedConcept: Category of the property the segment represents. @@ -228,7 +228,7 @@ def segmented_property_category(self)-> CodedConcept: return self.SegmentedPropertyCategoryCodeSequence[0] @property - def segmented_property_type(self)-> CodedConcept: + def segmented_property_type(self) -> CodedConcept: """highdicom.sr.coding.CodedConcept: Type of the property the segment represents. @@ -433,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 a985ca00..75f9828e 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -4,7 +4,9 @@ import numpy as np from operator import eq from collections import defaultdict -from typing import Any, Dict, Iterable, List, Optional, Set, Sequence, Union, Tuple +from typing import ( + Any, Dict, Iterable, List, Optional, Set, Sequence, Union, Tuple +) from pydicom.dataset import Dataset from pydicom.encaps import decode_data_sequence, encapsulate @@ -37,8 +39,10 @@ SegmentationFractionalTypeValues, SegmentationTypeValues, SegmentsOverlapValues, - SpatialLocationsPreservedValues + SpatialLocationsPreservedValues, + SegmentAlgorithmTypeValues, ) +from highdicom.seg.utils import iter_segments from highdicom.sr.coding import CodedConcept from highdicom.valuerep import check_person_name @@ -992,10 +996,48 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': 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 @@ -1068,10 +1110,16 @@ def _build_luts(self) -> None: if 'TrackingUID' in seg_info: self.tracking_uid_lut[seg_info.TrackingUID].append(n) - def is_binary(self) -> bool: - return ( - self.SegmentationType == - SegmentationTypeValues.BINARY.value + @property + def segmentation_type(self) -> SegmentationTypeValues: + return SegmentationTypeValues(self.SegmentationType) + + @property + def segmentation_fractional_type( + self + ) -> SegmentationFractionalTypeValues: + return SegmentationFractionalTypeValues( + self.SegmentationFractionalType ) def iter_segments(self): @@ -1086,7 +1134,10 @@ def segment_numbers(self) -> List[int]: # Do not assume segments are sorted (although they should be) return sorted(list(self.segment_description_lut.keys())) - def get_segment_description(self, segment_number: int) -> Dataset: + def get_segment_description( + self, + segment_number: int + ) -> SegmentDescription: if segment_number < 1: raise ValueError(f'{segment_number} is an invalid segment number') @@ -1097,119 +1148,52 @@ def get_segment_description(self, segment_number: int) -> Dataset: ) return self.segment_description_lut[segment_number] - @classmethod - def _coded_concept_sequence_to_code( - cls, - coded_concept: Dataset - ) -> Code: - scheme_version = ( - coded_concept.CodingSchemeVersion - if 'CodingSchemeVersion' in coded_concept else None - ) - return Code( - value=coded_concept.CodeValue, - meaning=coded_concept.CodeMeaning, - scheme_designator=coded_concept.CodingSchemeDesignator, - scheme_version=scheme_version - ) - - def get_segmented_property_category_code( + def get_segment_numbers( self, - segment_number: int - ) -> Code: - seg_desc = self.segment_description_lut[segment_number] - return self._coded_concept_sequence_to_code( - seg_desc.SegmentedPropertyCategoryCodeSequence[0] - ) - - def get_segmented_property_type_code(self, segment_number: int) -> Code: - seg_desc = self.segment_description_lut[segment_number] - return self._coded_concept_sequence_to_code( - seg_desc.SegmentedPropertyTypeCodeSequence[0] - ) - - def get_segment_number_property_mapping( - self - ) -> Dict[int, Tuple[Code, Code]]: - return { - n: - ( - self.get_segmented_property_category_code(n), - self.get_segmented_property_type_code(n) - ) - for n in self.segment_numbers - } - - def get_segment_numbers_for_segmented_property( - self, - segmented_property_category: Union[Code, CodedConcept, None] = None, - segmented_property_type: Union[Code, CodedConcept, None] = None, + 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]: - if ( - segmented_property_category is None and - segmented_property_type is None - ): - raise TypeError( - "At least one of the segmented_property_category and " - "segmented_property_type must be provided." + filter_funcs = [] + if segment_label is not None: + filter_funcs.append( + lambda desc: desc.segment_label == segment_label ) - if isinstance(segmented_property_category, CodedConcept): - segmented_property_category = self._coded_concept_sequence_to_code( - segmented_property_category + if segmented_property_category is not None: + filter_funcs.append( + lambda desc: + desc.segmented_property_category == segmented_property_category ) - if isinstance(segmented_property_type, CodedConcept): - segmented_property_type = self._coded_concept_sequence_to_code( - segmented_property_type + if segmented_property_type is not None: + filter_funcs.append( + lambda desc: + desc.segmented_property_type == segmented_property_type ) - - matched_segment_numbers: List[int] = [] - - for n in self.segment_numbers: - category_match = ( - segmented_property_category == - self.get_segmented_property_category_code(n) + if algorithm_type is not None: + algo_type = SegmentAlgorithmTypeValues(algorithm_type) + filter_funcs.append( + lambda desc: + SegmentAlgorithmTypeValues(desc.algorithm_type) == algo_type ) - type_match = ( - segmented_property_type == - self.get_segmented_property_type_code(n) + if tracking_uid is not None: + filter_funcs.append( + lambda desc: desc.tracking_uid == tracking_uid ) - - if category_match and type_match: - matched_segment_numbers.append(n) - - return matched_segment_numbers - - def get_segment_tracking_id(self, segment_number: int) -> Optional[str]: - seg_info = self.segment_description_lut[segment_number] - if 'TrackingID' in seg_info: - return self.segment_description_lut[segment_number].TrackingID - return None - - def get_segment_tracking_uid(self, segment_number: int) -> Optional[str]: - seg_info = self.segment_description_lut[segment_number] - if 'TrackingUID' in seg_info: - return self.segment_description_lut[segment_number].TrackingUID - return None - - def get_segment_numbers_for_tracking_id( - self, - tracking_id: str - ) -> List[int]: - if tracking_id not in self.tracking_id_lut: - raise KeyError( - 'No segment found matching specified tracking ID' + if tracking_id is not None: + filter_funcs.append( + lambda desc: desc.tracking_id == tracking_id ) - return self.tracking_id_lut[tracking_id] - def get_segment_numbers_for_tracking_uid( - self, - tracking_uid: str - ) -> List[int]: - if tracking_uid not in self.tracking_uid_lut: - raise KeyError( - 'No segment found matching specified tracking UID' - ) - return self.tracking_uid_lut[tracking_uid] + matches = [ + desc.segment_number + for desc in self.SegmentSequence + if all(f(desc) for f in filter_funcs) + ] + + return matches def get_tracking_ids(self) -> List[str]: return list(self.tracking_id_lut.keys()) @@ -1438,7 +1422,6 @@ def get_pixels( return pixel_array # TODO - # Shift segment information getters to SegmentDescription # Multi-frame source image # Segs with a single frame # Allow binary fractional segs to combine segments @@ -1456,4 +1439,3 @@ def get_pixels( # Correct for spatial ordering of source frames? # Optimize combine_segments to exploit sparse nature of array # Ensure output array cannot be used to change value of pixel_array - diff --git a/src/highdicom/sr/coding.py b/src/highdicom/sr/coding.py index 75d402a2..97b7cd63 100644 --- a/src/highdicom/sr/coding.py +++ b/src/highdicom/sr/coding.py @@ -110,7 +110,7 @@ def from_dataset(cls, dataset: Dataset) -> 'CodedConcept': if not hasattr(dataset, kw): raise AttributeError( 'Dataset does not contain the following attribute ' - f'required for coded concepts: {attr}.' + f'required for coded concepts: {kw}.' ) return cls( value=dataset.CodeValue, From 1ec456f72a9fad72a3227122f2acc46f362d766e Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 16 May 2021 22:53:26 -0400 Subject: [PATCH 07/31] Fix outdated method call --- src/highdicom/seg/sop.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 75f9828e..47c67bcb 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1334,10 +1334,11 @@ def get_pixels( pixel array will always lie in the range ``0`` to ``len(segment_numbers)`` inclusive. """ - if combine_segments and not self.is_binary(): - raise ValueError( - 'Cannot combine segments if the segmentation is not binary' - ) + if combine_segments: + if self.segmentation_type != SegmentationTypeValues.BINARY: + raise ValueError( + 'Cannot combine segments if the segmentation is not binary' + ) # If no source frames were specified, use all source frames if source_frames is None: @@ -1438,4 +1439,4 @@ def get_pixels( # List segmented property for tracking id and vice versa # Correct for spatial ordering of source frames? # Optimize combine_segments to exploit sparse nature of array - # Ensure output array cannot be used to change value of pixel_array + # Implement from_dataset from DimensionIndexSequence From 326289137400f08925644d64819d449377835d53 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Wed, 19 May 2021 16:45:02 -0400 Subject: [PATCH 08/31] Docstrings and refactoring of segmentation object parsing --- src/highdicom/seg/sop.py | 135 +++++++++++++++++++++++++++++---------- 1 file changed, 102 insertions(+), 33 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 47c67bcb..7ea8781e 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -29,6 +29,7 @@ PlanePositionSequence, PixelMeasuresSequence ) +from highdicom.errors import DicomComplianceError from highdicom.enum import CoordinateSystemNames from highdicom.frame import encode_frame from highdicom.seg.content import ( @@ -958,6 +959,20 @@ def _encode_pixels(self, planes: np.ndarray) -> bytes: @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.' @@ -981,7 +996,7 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': seg._coordinate_system = CoordinateSystemNames.PATIENT seg._plane_orientation = plane_ori_seq.ImageOrientationPatient else: - raise ValueError( + raise DicomComplianceError( 'Expected Plane Orientation Sequence to have either ' 'ImageOrientationSlide or ImageOrientationPatient ' 'attribute.' @@ -989,10 +1004,11 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': for i, segment in enumerate(seg.SegmentSequence, 1): if segment.SegmentNumber != i: - raise ValueError( + raise DicomComplianceError( 'Segments are expected to start at 1 and be consecutive ' 'integers.' ) + seg._segment_inventory = { s.SegmentNumber for s in seg.SegmentSequence } @@ -1049,10 +1065,6 @@ def _build_luts(self) -> None: ref_ins.ReferencedSOPInstanceUID for ref_ins in ref_series.ReferencedInstanceSequence ]) - self.segment_description_lut = { - int(item.SegmentNumber): item - for item in self.SegmentSequence - } segnum_col_data = [] source_frame_col_data = [] @@ -1098,26 +1110,22 @@ def _build_luts(self) -> None: self.lut_seg_col = 0 self.lut_src_col = 1 - # Build LUT from tracking ID -> segment numbers - self.tracking_id_lut: Dict[str, List[int]] = defaultdict(list) - for n, seg_info in self.segment_description_lut.items(): - if 'TrackingID' in seg_info: - self.tracking_id_lut[seg_info.TrackingID].append(n) - - # Build LUT from tracking UID -> segment numbers - self.tracking_uid_lut: Dict[str, int] = defaultdict(list) - for n, seg_info in self.segment_description_lut.items(): - if 'TrackingUID' in seg_info: - self.tracking_uid_lut[seg_info.TrackingUID].append(n) - @property def segmentation_type(self) -> SegmentationTypeValues: + """highdicom.seg.SegmentationTypeValues: Segmentation type.""" return SegmentationTypeValues(self.SegmentationType) @property def segmentation_fractional_type( self - ) -> SegmentationFractionalTypeValues: + ) -> Optional[SegmentationFractionalTypeValues]: + """ + highdicom.seg.SegmentationFractionalTypeValues: + Segmentation fractional type. + + """ + if not hasattr(self, 'SegmentationFractionalType'): + return None return SegmentationFractionalTypeValues( self.SegmentationFractionalType ) @@ -1127,28 +1135,39 @@ def 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) -> List[int]: - # Do not assume segments are sorted (although they should be) - return sorted(list(self.segment_description_lut.keys())) + 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. - if segment_number < 1: - raise ValueError(f'{segment_number} is an invalid segment number') + Parameters + ---------- + segment_number: int + Segment number for the segment, as a 1-based index. - if segment_number not in self.segment_description_lut: - raise KeyError( - f'No segment number {segment_number} found in dataset' + 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.segment_description_lut[segment_number] + return self.SegmentSequence[segment_number - 1] - def get_segment_numbers( + def search_segments( self, segment_label: Optional[str] = None, segmented_property_category: Optional[Union[Code, CodedConcept]] = None, @@ -1157,6 +1176,32 @@ def get_segment_numbers( 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: Optional[str] + Segment label filter to apply. + segmented_property_category: Optional[Union[Code, CodedConcept]] + Segmented property category filter to apply. + segmented_property_type: Optional[Union[Code, CodedConcept]] + Segmented property type filter to apply. + algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]] + Segmented property type filter to apply. + tracking_uid: Optional[str] + Tracking unique identifier filter to apply. + tracking_id: Optional[str] + Tracking identifier filter to apply. + + Returns + ------- + List[int] + List of all segment numbers matching the provided criteria. + + """ filter_funcs = [] if segment_label is not None: filter_funcs.append( @@ -1195,11 +1240,35 @@ def get_segment_numbers( return matches - def get_tracking_ids(self) -> List[str]: - return list(self.tracking_id_lut.keys()) + def all_tracking_ids(self) -> Set[str]: + """Get all unique tracking identifiers in this SEG image. - def get_tracking_uids(self) -> List[str]: - return list(self.tracking_uid_lut.keys()) + Returns + ------- + Set[str] + All unique tracking identifiers referenced in segment descriptions + in this SEG image. + + """ + return { + desc.tracking_id for desc in self.SegmentSequence + if desc.tracking_id is not None + } + + def all_tracking_uids(self) -> Set[str]: + """Get all unique tracking unique identifiers in this SEG image. + + Returns + ------- + Set[str] + All unique tracking unique identifiers referenced in segment + descriptions in this SEG image. + + """ + return { + desc.tracking_uid for desc in self.SegmentSequence + if desc.tracking_uid is not None + } def _get_src_fm_index(self, sop_instance_uid: str) -> int: ind = np.argwhere(self.source_sop_instance_uids == sop_instance_uid) From 9ac18ce68ae445840bc927bf3ba8a0e7cb123256 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Wed, 19 May 2021 23:40:35 -0400 Subject: [PATCH 09/31] Work in progress on altering indexing logic --- src/highdicom/errors.py | 10 ++ src/highdicom/seg/sop.py | 340 ++++++++++++++++++++++++++++++++++----- 2 files changed, 310 insertions(+), 40 deletions(-) create mode 100644 src/highdicom/errors.py diff --git a/src/highdicom/errors.py b/src/highdicom/errors.py new file mode 100644 index 00000000..0d4eef1e --- /dev/null +++ b/src/highdicom/errors.py @@ -0,0 +1,10 @@ +"""Errors raised by highdicom processes.""" + +class DicomComplianceError(Exception): + """DICOM standard compliance error. + + Exception indicating that a user-provided DICOM dataset is not in + compliance with the DICOM standard. + + """ + pass diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 7ea8781e..62cf74b1 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1,5 +1,6 @@ """Module for the SOP class of the Segmentation IOD.""" from copy import deepcopy +from collections import OrderedDict import logging import numpy as np from operator import eq @@ -10,6 +11,7 @@ from pydicom.dataset import Dataset 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.uid import ( @@ -51,6 +53,9 @@ logger = logging.getLogger(__name__) +_NO_FRAME_REF_VALUE = -1 + + class Segmentation(SOPClass): """SOP class for a Segmentation, which represents one or more @@ -1058,57 +1063,133 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': return seg + def _build_ref_instance_lut(self) -> None: + # 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] = ( + self.StudyInstanceUID, + ref_series.SeriesInstanceUID, + 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] = ( + ref_study.StudyInstanceUID, + ref_series.SeriesInstanceUID, + ref_ins.ReferencedSOPInstanceUID + ) + + self._source_sop_instance_uids = np.array(list(self._ref_ins_lut.keys())) + def _build_luts(self) -> None: - ref_series = self.ReferencedSeriesSequence[0] - self.source_sop_instance_uids = np.array([ - ref_ins.ReferencedSOPInstanceUID for ref_ins in - ref_series.ReferencedInstanceSequence - ]) + self._build_ref_instance_lut() segnum_col_data = [] + source_instance_col_data = [] source_frame_col_data = [] + # Create a list of source images and check for spatial locations + # preserved and that there is a single source frame per seg frame + locations_preserved = [] + self._single_source_frame_per_seg_frame = True for frame_item in self.PerFrameFunctionalGroupsSequence: # Get segment number for this frame - seg_no = frame_item.SegmentIdentificationSequence[0].\ - ReferencedSegmentNumber - segnum_col_data.append(int(seg_no)) + seg_id_seg = frame_item.SegmentIdentificationSequence[0] + seg_num = seg_id_seg.ReferencedSegmentNumber + segnum_col_data.append(int(seg_num)) - source_images = [] + frame_source_instances = [] + frame_source_frames = [] for der_im in frame_item.DerivationImageSequence: for src_im in der_im.SourceImageSequence: - if 'SpatialLocationsPreserved' in src_im: - if eq( - src_im.SpatialLocationsPreserved, - SpatialLocationsPreservedValues.NO.value + frame_source_instances.append( + src_im.ReferencedSOPInstanceUID + ) + if hasattr(src_im, 'SpatialLocationsPreserved'): + locations_preserved.append( + SpatialLocationsPreservedValues( + src_im.SpatialLocationsPreserved + ) + ) + else: + locations_preserved.append( + SpatialLocationsPreservedValues.UNKNOWN + ) + + if hasattr(src_im, 'ReferencedFrameNumber'): + if isinstance( + src_im.ReferencedFrameNumber, + MultiValue ): - raise RuntimeError( - 'Segmentation dataset specifies that spatial ' - 'locations are not preserved' + frame_source_frames.extend( + [ + int(f) + for f in src_im.ReferencedFrameNumber + ] ) - source_images.append(src_im.ReferencedSOPInstanceUID) - if len(set(source_images)) == 0: - raise RuntimeError( - 'Could not find source image information for ' - 'segmentation frames' - ) - if len(set(source_images)) > 1: - raise RuntimeError( - "Multiple source images found for segmentation frames" - ) - src_fm_ind = self._get_src_fm_index(source_images[0]) - source_frame_col_data.append(src_fm_ind) + 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 DicomComplianceError( + 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_fm_ind = self._get_src_fm_index(frame_source_instances[0]) + source_instance_col_data.append(src_fm_ind) + source_frame_col_data.append(frame_source_frames[0]) + + # Summarise + if any( + v == SpatialLocationsPreservedValues.NO + for v in locations_preserved + ): + self._locations_preserved = SpatialLocationsPreservedValues.NO + elif all( + v == SpatialLocationsPreservedValues.YES + for v in locations_preserved + ): + self._locations_preserved = SpatialLocationsPreservedValues.YES + else: + self._locations_preserved = SpatialLocationsPreservedValues.UNKNOWN - # Frame LUT is a 2D numpy array with two columns + # Frame LUT is a 2D numpy array with three columns # Row i represents frame i of the segmentation dataset - # The two columns represent the segment number and the source frame - # index in the source_sop_instance_uids list + # The three columns represent the segment number, the source frame + # index in the source_sop_instance_uids list, and the source frame + # number (if applicable) # This allows for fairly efficient querying by any of the three - # variables: seg frame number, source frame number, segment number - self.frame_lut = np.array([segnum_col_data, source_frame_col_data]).T - self.lut_seg_col = 0 - self.lut_src_col = 1 + # variables: seg frame number, source instance/frame number, segment + # number + if self._single_source_frame_per_seg_frame: + self._frame_lut = np.array([ + segnum_col_data, + source_instance_col_data, + source_frame_col_data + ]).T + self._lut_seg_col = 0 + self._lut_src_instance_col = 1 + self._lut_src_frame_col = 2 @property def segmentation_type(self) -> SegmentationTypeValues: @@ -1167,7 +1248,7 @@ def get_segment_description( ) return self.SegmentSequence[segment_number - 1] - def search_segments( + def search_for_segments( self, segment_label: Optional[str] = None, segmented_property_category: Optional[Union[Code, CodedConcept]] = None, @@ -1270,8 +1351,42 @@ def all_tracking_uids(self) -> Set[str]: if desc.tracking_uid is not None } + def all_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 + + def all_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_fm_index(self, sop_instance_uid: str) -> int: - ind = np.argwhere(self.source_sop_instance_uids == sop_instance_uid) + 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}' @@ -1324,9 +1439,151 @@ def list_source_frames_for_segment_number( self.frame_lut[:, self.lut_seg_col] == segment_number, self.lut_src_col ] - source_frame_uids = self.source_sop_instance_uids[source_frame_indices] + source_frame_uids = self._source_sop_instance_uids[source_frame_indices] return source_frame_uids.tolist() + def _get_pixels_by_seg_frame( + self, + seg_frames_matrix: np.ndarray, + segment_numbers: np.array, + combine_segments: bool = False, + relabel: bool = False, + ) -> np.ndarray: + """ + + 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 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. + + """ + if combine_segments: + if self.segmentation_type != SegmentationTypeValues.BINARY: + raise ValueError( + 'Cannot combine segments if the segmentation is not binary' + ) + + # 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.' + '-1.' + ) + + # 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 combine_segments: + # 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_pixels_by_source_frame( + source_sop_instance_uids: Sequence[str], + source_frame_numbers: Optional[Sequence[int]] = None, + segment_numbers: Optional[Iterable[int]] = None, + combine_segments: bool = False, + relabel: bool = False, + assert_locations_preserved: bool = False + ): + if self._locations_preserved == SpatialLocationsPreservedValues.NO: + 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 eq( + self._locations_preserved, + SpatialLocationsPreservedValues.UNKNOWN + ): + 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, you may override this behavior with the ' + "'assert_locations_preserved' 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( self, source_frames: Optional[Iterable[str]] = None, @@ -1370,7 +1627,7 @@ def get_pixels( has been specified, then ``pixel_array[i, ...]`` represents the segmentation of ``source_frames[i]``. If ``source_frames`` was not specified, then ``pixel_array[i, ...]`` represents the segmentation of - ``parser.source_sop_instance_uids[i]``. + ``parser._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 @@ -1411,7 +1668,7 @@ def get_pixels( # If no source frames were specified, use all source frames if source_frames is None: - source_frames = self.source_sop_instance_uids + source_frames = self._source_sop_instance_uids # If no segments were specified, use all segments if segment_numbers is None: @@ -1509,3 +1766,6 @@ def get_pixels( # Correct for spatial ordering of source frames? # Optimize combine_segments to exploit sparse nature of array # Implement from_dataset from DimensionIndexSequence + # Allowing indexing from multiple source series + # Index by spatial location + # Get fill pixel array with metadata From d10af53bd347aa1ddc23b57121cfd1cd900e6bb0 Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Fri, 21 May 2021 17:56:25 -0400 Subject: [PATCH 10/31] Got indexing by source frame working --- src/highdicom/errors.py | 5 +- src/highdicom/seg/sop.py | 200 +++++++++++++++++++++++++++++++++------ 2 files changed, 173 insertions(+), 32 deletions(-) diff --git a/src/highdicom/errors.py b/src/highdicom/errors.py index 0d4eef1e..841e7e79 100644 --- a/src/highdicom/errors.py +++ b/src/highdicom/errors.py @@ -1,10 +1,11 @@ """Errors raised by highdicom processes.""" -class DicomComplianceError(Exception): + +class DicomAttributeError(Exception): """DICOM standard compliance error. Exception indicating that a user-provided DICOM dataset is not in - compliance with the DICOM standard. + compliance with the DICOM standard due to a missing attribute. """ pass diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 62cf74b1..9df6dd68 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -3,7 +3,6 @@ from collections import OrderedDict import logging import numpy as np -from operator import eq from collections import defaultdict from typing import ( Any, Dict, Iterable, List, Optional, Set, Sequence, Union, Tuple @@ -31,7 +30,7 @@ PlanePositionSequence, PixelMeasuresSequence ) -from highdicom.errors import DicomComplianceError +from highdicom.errors import DicomAttributeError from highdicom.enum import CoordinateSystemNames from highdicom.frame import encode_frame from highdicom.seg.content import ( @@ -1001,7 +1000,7 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': seg._coordinate_system = CoordinateSystemNames.PATIENT seg._plane_orientation = plane_ori_seq.ImageOrientationPatient else: - raise DicomComplianceError( + raise DicomAttributeError( 'Expected Plane Orientation Sequence to have either ' 'ImageOrientationSlide or ImageOrientationPatient ' 'attribute.' @@ -1009,7 +1008,7 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': for i, segment in enumerate(seg.SegmentSequence, 1): if segment.SegmentNumber != i: - raise DicomComplianceError( + raise DicomAttributeError( 'Segments are expected to start at 1 and be consecutive ' 'integers.' ) @@ -1085,7 +1084,9 @@ def _build_ref_instance_lut(self) -> None: ref_ins.ReferencedSOPInstanceUID ) - self._source_sop_instance_uids = np.array(list(self._ref_ins_lut.keys())) + self._source_sop_instance_uids = np.array( + list(self._ref_ins_lut.keys()) + ) def _build_luts(self) -> None: @@ -1149,14 +1150,14 @@ def _build_luts(self) -> None: else: ref_instance_uid = frame_source_instances[0] if ref_instance_uid not in self._source_sop_instance_uids: - raise DicomComplianceError( + raise DicomAttributeError( 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_fm_ind = self._get_src_fm_index(frame_source_instances[0]) - source_instance_col_data.append(src_fm_ind) + 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 @@ -1164,18 +1165,19 @@ def _build_luts(self) -> None: v == SpatialLocationsPreservedValues.NO for v in locations_preserved ): - self._locations_preserved = SpatialLocationsPreservedValues.NO + Type = Optional[SpatialLocationsPreservedValues] + self._locations_preserved: Type = SpatialLocationsPreservedValues.NO elif all( v == SpatialLocationsPreservedValues.YES for v in locations_preserved ): self._locations_preserved = SpatialLocationsPreservedValues.YES else: - self._locations_preserved = SpatialLocationsPreservedValues.UNKNOWN + self._locations_preserved = None # Frame LUT is a 2D numpy array with three columns # Row i represents frame i of the segmentation dataset - # The three columns represent the segment number, the source frame + # The three columns represent the segment number, the source uid # index in the source_sop_instance_uids list, and the source frame # number (if applicable) # This allows for fairly efficient querying by any of the three @@ -1385,7 +1387,7 @@ def all_segmented_property_types(self) -> List[CodedConcept]: return types - def _get_src_fm_index(self, sop_instance_uid: str) -> int: + 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( @@ -1397,7 +1399,7 @@ def list_seg_frames_for_source_frame( self, source_sop_instance_uid: str ) -> List[int]: - src_ind = self._get_src_fm_index(source_sop_instance_uid) + src_ind = self._get_src_uid_index(source_sop_instance_uid) seg_frames = np.where( self.frame_lut[:, self.lut_src_col] == src_ind )[0] @@ -1407,7 +1409,7 @@ def list_segments_in_source_frame( self, source_sop_instance_uid: str ) -> List[int]: - src_ind = self._get_src_fm_index(source_sop_instance_uid) + src_ind = self._get_src_uid_index(source_sop_instance_uid) segments = self.frame_lut[ self.frame_lut[:, self.lut_src_col] == src_ind, self.lut_seg_col @@ -1550,40 +1552,179 @@ def _get_pixels_by_seg_frame( return pixel_array + def get_source_sop_instances(self): + return [self._ref_ins_lut[sop_uid] for sop_uid in self._ref_ins_lut] + def get_pixels_by_source_frame( + self, source_sop_instance_uids: Sequence[str], source_frame_numbers: Optional[Sequence[int]] = None, segment_numbers: Optional[Iterable[int]] = None, combine_segments: bool = False, relabel: bool = False, - assert_locations_preserved: bool = False + assert_locations_preserved: bool = False, + assert_missing_frames_are_empty: bool = True ): - if self._locations_preserved == SpatialLocationsPreservedValues.NO: + # Checks that it is possible to index using source frames in this + # dataset + if self._locations_preserved is None: + if not assert_locations_preserved: + 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, you may override this behavior ' + "with the 'assert_locations_preserved' parameter." + ) + elif self._locations_preserved == SpatialLocationsPreservedValues.NO: 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 eq( - self._locations_preserved, - SpatialLocationsPreservedValues.UNKNOWN - ): - 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, you may override this behavior with the ' - "'assert_locations_preserved' 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.' ) + # 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_sop_instance_uids) == 0: + raise ValueError( + 'Source SOP instance UIDs may not be empty.' + ) + + if source_frame_numbers is None: + # 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, assume this frame + # is empty and leave the entry in seg_frames as -1 + + else: + if len(source_sop_instance_uids) != 1: + raise ValueError( + 'Only a single source SOP instance UID may be specified ' + 'when frame numbers are specified.' + ) + 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_uids[0]) + 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[:, 1]: + 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 KeyError(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, 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, + ) + def get_pixels( self, source_frames: Optional[Iterable[str]] = None, @@ -1762,7 +1903,6 @@ def get_pixels( # Lazy decoding of segmentation frames # Should functions raise error or return empty lists? # Standardize method names - # List segmented property for tracking id and vice versa # Correct for spatial ordering of source frames? # Optimize combine_segments to exploit sparse nature of array # Implement from_dataset from DimensionIndexSequence From 70fe01c9d7ca8761e6cec5affda153feff68295a Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Thu, 27 May 2021 17:21:55 -0400 Subject: [PATCH 11/31] WIP --- src/highdicom/seg/sop.py | 877 +++++++++++++++++++++++++++------------ 1 file changed, 602 insertions(+), 275 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 9df6dd68..8d31322c 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -9,6 +9,7 @@ ) from pydicom.dataset import Dataset +from pydicom.datadict import 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 @@ -1096,6 +1097,16 @@ def _build_luts(self) -> None: 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_index_col_data = {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_preserved = [] @@ -1106,6 +1117,18 @@ def _build_luts(self) -> None: 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, ind in zip(self._dim_ind_pointers, indices): + dim_index_col_data[ptr].append(ind) + frame_source_instances = [] frame_source_frames = [] for der_im in frame_item.DerivationImageSequence: @@ -1175,23 +1198,36 @@ def _build_luts(self) -> None: else: self._locations_preserved = None - # Frame LUT is a 2D numpy array with three columns - # Row i represents frame i of the segmentation dataset - # The three columns represent the segment number, the source uid - # index in the source_sop_instance_uids list, and the source frame - # number (if applicable) + # 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 if self._single_source_frame_per_seg_frame: - self._frame_lut = np.array([ + # Include columns related + col_data = [ segnum_col_data, source_instance_col_data, source_frame_col_data - ]).T + ] self._lut_seg_col = 0 self._lut_src_instance_col = 1 self._lut_src_frame_col = 2 + else: + self._lut_seg_col = 0 + col_data = [segnum_col_data] + + self._lut_dim_ind_cols = { + ptr: i for ptr, i in enumerate(self._dim_ind_pointers, len(col_data)) + } + col_data += [ + indices for indices in dim_index_col_data.values() + ] + self._frame_lut = np.array(col_data).T @property def segmentation_type(self) -> SegmentationTypeValues: @@ -1451,27 +1487,43 @@ def _get_pixels_by_seg_frame( combine_segments: bool = False, relabel: bool = False, ) -> 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 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. + 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. - """ - if combine_segments: - if self.segmentation_type != SegmentationTypeValues.BINARY: - raise ValueError( - 'Cannot combine segments if the segmentation is not binary' - ) + 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.') @@ -1533,6 +1585,20 @@ def _get_pixels_by_seg_frame( self.pixel_array[seg_frame_ind, :, :] if combine_segments: + # Check whether segmentation is binary, or fractional with only + # binary values + if self.segmentation_type != SegmentationTypeValues.BINARY: + is_binary = np.isin( + pixel_array.unique(), + 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( @@ -1555,41 +1621,166 @@ def _get_pixels_by_seg_frame( def get_source_sop_instances(self): return [self._ref_ins_lut[sop_uid] for sop_uid in self._ref_ins_lut] - def get_pixels_by_source_frame( + def _check_indexing_with_source_frames( self, - source_sop_instance_uids: Sequence[str], - source_frame_numbers: Optional[Sequence[int]] = None, - segment_numbers: Optional[Iterable[int]] = None, - combine_segments: bool = False, - relabel: bool = False, - assert_locations_preserved: bool = False, - assert_missing_frames_are_empty: bool = True - ): + 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 assert_locations_preserved: + 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, you may override this behavior ' - "with the 'assert_locations_preserved' parameter." + """ + 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: - 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 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[Iterable[int]] = None, + combine_segments: bool = False, + relabel: bool = False, + ignore_spatial_locations: bool = False, + assert_missing_frames_are_empty: bool = False + ) -> np.ndarray: + """Get a pixel array for a list of source instances. + + This is intended for retrieving segmentation masks derived from a + 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: Optional[Sequence[int]] + Sequence containing segment numbers to include. If unspecified, + all segments are included. + 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. + ignore_spatial_locations: bool + 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 + 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. + + Returns + ------- + pixel_array: np.ndarray + Pixel array representing the segmentation. See notes for full + explanation. + + """ + # 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) @@ -1602,121 +1793,60 @@ def get_pixels_by_source_frame( 'Source SOP instance UIDs may not be empty.' ) - if source_frame_numbers is None: - # 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) + # 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] + # 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, assume this frame - # is empty and leave the entry in seg_frames as -1 - - else: - if len(source_sop_instance_uids) != 1: - raise ValueError( - 'Only a single source SOP instance UID may be specified ' - 'when frame numbers are specified.' - ) - 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_uids[0]) - 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.' - ) + # 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, src_frm_num in enumerate(source_frame_numbers): - # Check whether this source frame exists in the LUT - if src_frm_num not in lut[:, 1]: - 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 KeyError(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, assume this frame - # is empty and leave the entry in seg_frames as -1 + # 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, @@ -1725,21 +1855,69 @@ def get_pixels_by_source_frame( relabel=relabel, ) - def get_pixels( + def get_pixels_by_source_frame( self, - source_frames: Optional[Iterable[str]] = None, + source_sop_instance_uid: str, + source_frame_numbers: Sequence[int], segment_numbers: Optional[Iterable[int]] = None, combine_segments: bool = False, relabel: bool = False, - ) -> np.ndarray: - """Get a numpy array of the reconstructed segmentation. + ignore_spatial_locations: bool = False, + assert_missing_frames_are_empty: bool = False + ): + """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_frames: Optional[Iterable[str]] - Iterable containing SOP Instance UIDs of the source images to - include in the segmentation. If unspecified, all source images - are included. + 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: Optional[Sequence[int]] Sequence containing segment numbers to include. If unspecified, all segments are included. @@ -1754,23 +1932,142 @@ def get_pixels( ``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 + 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 + 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. Returns ------- pixel_array: np.ndarray - Pixel array representing the segmentation + Pixel array representing the segmentation. See notes for full + explanation. - Note - ---- - 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. If ``source_frames`` - has been specified, then ``pixel_array[i, ...]`` represents the - segmentation of ``source_frames[i]``. If ``source_frames`` was not - specified, then ``pixel_array[i, ...]`` represents the segmentation of - ``parser._source_sop_instance_uids[i]``. - The next two dimensions are the rows and columns of the frames, + """ + # 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[:, 1]: + 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 KeyError(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, + ) + + def get_pixels_by_dimension_index_values( + self, + dimension_index_values: Sequence[Sequence[int]], + dimension_index_pointers: Optional[Sequence[int]] = None, + segment_numbers: Optional[Iterable[int]] = None, + combine_segments: bool = False, + relabel: bool = False, + assert_missing_frames_are_empty: bool = False + ): + """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 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]`` @@ -1780,132 +2077,162 @@ def get_pixels( 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 type of the segmentation is ``BINARY`` and - 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. - """ - if combine_segments: - if self.segmentation_type != SegmentationTypeValues.BINARY: - raise ValueError( - 'Cannot combine segments if the segmentation is not binary' - ) + 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. - # If no source frames were specified, use all source frames - if source_frames is None: - source_frames = self._source_sop_instance_uids - - # If no segments were specified, use all segments - if segment_numbers is None: - segment_numbers = self.segment_numbers + Parameters + ---------- + dimension_index_values: Sequence[Sequence[int]] + TODO + dimension_index_pointers: Optional[Sequence[int] + TODO + segment_numbers: Optional[Sequence[int]] + Sequence containing segment numbers to include. If unspecified, + all segments are included. + 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. + assert_missing_frames_are_empty: bool + 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. - # Check all provided segmentation numbers are valid - seg_num_arr = np.array(segment_numbers) - seg_nums_exist = np.in1d( - seg_num_arr, np.array(self.segment_numbers) - ) - if not seg_nums_exist.all(): - missing_seg_nums = seg_num_arr[np.logical_not(seg_nums_exist)] - raise KeyError( - 'Segment numbers contained non-existent segments: ' - f'{missing_seg_nums}' - ) + Returns + ------- + pixel_array: np.ndarray + Pixel array representing the segmentation. See notes for full + explanation. - # Check segment numbers are unique - if len(np.unique(seg_num_arr)) != len(seg_num_arr): + """ + # 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 must not contain duplicates' + 'Segment numbers may not be empty.' ) - # Initialize empty pixel array - pixel_array = np.zeros( - ( - len(source_frames), - self.Rows, - self.Columns, - len(segment_numbers) - ), - self.pixel_array.dtype - ) - - # Loop through source frames - for out_ind, src_frame in enumerate(source_frames): - src_ind = self._get_src_fm_index(src_frame) - - # Find segmentation frames relating to this source frame - seg_frames = np.where( - self.frame_lut[:, self.lut_src_col] == src_ind - )[0] - seg_nums_for_frame = self.frame_lut[seg_frames, self.lut_seg_col] - - # Loop through segmentation frames - for seg_frame, seg_num in zip(seg_frames, seg_nums_for_frame): - # Find output index for this segmentation number (if any) - seg_ind = np.where(seg_num_arr == seg_num)[0] - - if len(seg_ind) > 0: - # Copy data to to output array - if self.pixel_array.ndim == 2: - # Special case with a single segmentation frame - pixel_array[out_ind, :, :, seg_ind.item()] = \ - self.pixel_array - else: - pixel_array[out_ind, :, :, seg_ind.item()] = \ - self.pixel_array[seg_frame, :, :] + 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 {ptr} ({kw}) is not used as a dimension index ' + 'in this image.' + ) - if combine_segments: - # 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' - ) + if len(dimension_index_values) == 0: + raise ValueError( + 'Dimension index values 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 + query_cols = [ + self._dim_ind_pointers.index(ptr) + for ptr in dimension_index_pointers + ] + lut = lut[:, query_cols] - # 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 = seg_num_arr - scaled_array = pixel_array * pixel_value_map.reshape(1, 1, 1, -1) + 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.' + ) - # Combine segments by taking maximum down final dimension - max_array = scaled_array.max(axis=-1) - pixel_array = max_array + # Build the segmentation frame matrix + for r, ind_vals in enumerate(dimension_index_values): + # Check whether this source frame exists in the LUT + if src_frm_num not in lut[:, 1]: + 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 KeyError(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 pixel_array + return self._get_pixels_by_seg_frame( + seg_frames_matrix=seg_frames, + segment_numbers=np.array(segment_numbers), + combine_segments=combine_segments, + relabel=relabel, + ) # TODO - # Multi-frame source image - # Segs with a single frame - # Allow binary fractional segs to combine segments - # Option to force standards compliance + # Errors? # Integrity checks: - # Each source frame specified in pffgs exists - # No combination of source frame and segment number has - # multiple source frames # Tracking IDs do not have multiple segments with different # segmented properties # Lazy decoding of segmentation frames - # Should functions raise error or return empty lists? # Standardize method names - # Correct for spatial ordering of source frames? # Optimize combine_segments to exploit sparse nature of array # Implement from_dataset from DimensionIndexSequence # Allowing indexing from multiple source series # Index by spatial location # Get fill pixel array with metadata + # method to access seg pffgs? From deaa590832a955ff016d2644fdb8db591da5e30d Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Thu, 27 May 2021 18:13:35 -0400 Subject: [PATCH 12/31] Implemented indexing by dimension index values. Needs docstrings --- src/highdicom/seg/sop.py | 60 ++++++++++++++++++++-------------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 8d31322c..ddae3dcb 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1207,26 +1207,28 @@ def _build_luts(self) -> None: # This allows for fairly efficient querying by any of the three # variables: seg frame number, source instance/frame number, segment # number - if self._single_source_frame_per_seg_frame: - # Include columns related - col_data = [ - segnum_col_data, - source_instance_col_data, - source_frame_col_data - ] - self._lut_seg_col = 0 - self._lut_src_instance_col = 1 - self._lut_src_frame_col = 2 - else: - self._lut_seg_col = 0 - col_data = [segnum_col_data] + # 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 = { ptr: i 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 @@ -2103,7 +2105,7 @@ def get_pixels_by_dimension_index_values( ---------- dimension_index_values: Sequence[Sequence[int]] TODO - dimension_index_pointers: Optional[Sequence[int] + dimension_index_pointers: Optional[Sequence[TODO]] TODO segment_numbers: Optional[Sequence[int]] Sequence containing segment numbers to include. If unspecified, @@ -2171,7 +2173,7 @@ def get_pixels_by_dimension_index_values( self._dim_ind_pointers.index(ptr) for ptr in dimension_index_pointers ] - lut = lut[:, query_cols] + lut = self._frame_lut[:, query_cols + [self._lut_seg_col]] if np.unique(lut, axis=0).shape[0] != lut.shape[0]: raise RuntimeError( @@ -2182,25 +2184,22 @@ def get_pixels_by_dimension_index_values( # Build the segmentation frame matrix for r, ind_vals in enumerate(dimension_index_values): - # Check whether this source frame exists in the LUT - if src_frm_num not in lut[:, 1]: - 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 KeyError(msg) + # TODO Check whether any frame exists in the LUT with these indices + 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.' + ) # 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]) + # the index values and segment number + qry = np.array(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 @@ -2216,6 +2215,7 @@ def get_pixels_by_dimension_index_values( # assume this frame is empty and leave the entry in seg_frames # as -1 + print(seg_frames) return self._get_pixels_by_seg_frame( seg_frames_matrix=seg_frames, segment_numbers=np.array(segment_numbers), From 4581715395cd07de395d0c64f478a3019b6a6d6b Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 30 May 2021 17:04:47 -0400 Subject: [PATCH 13/31] Completed documentation, fixed indexing bugs --- src/highdicom/content.py | 10 +- src/highdicom/module_utils.py | 17 +- src/highdicom/seg/sop.py | 458 ++++++++++++++++++++++++---------- src/highdicom/uid.py | 11 +- 4 files changed, 343 insertions(+), 153 deletions(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index dcdde048..872b16b9 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -340,7 +340,10 @@ def __eq__(self, other: Any) -> bool: ) @classmethod - def from_sequence(cls, sequence: DataElementSequence) -> 'PlanePositionSequence': + 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. @@ -467,7 +470,10 @@ def __eq__(self, other: Any) -> bool: return False @classmethod - def from_sequence(cls, sequence: DataElementSequence) -> 'PlaneOrientationSequence': + 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. diff --git a/src/highdicom/module_utils.py b/src/highdicom/module_utils.py index 60b061b2..57357f47 100644 --- a/src/highdicom/module_utils.py +++ b/src/highdicom/module_utils.py @@ -81,11 +81,12 @@ def check_required_attributes( # Construct tree once and re-use in all recursive calls tree = construct_module_tree(module) - for p in base_path: - try: - tree = tree['attributes'][p] - except KeyError: - raise AttributeError(f"Invalid base path: {base_path}.") + 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( @@ -112,8 +113,8 @@ def check( ) if recursive: sequence_exists = ( - 'attributes' in subtree['attributes'][kw] - and hasattr(dataset, kw) + '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 @@ -155,7 +156,7 @@ def construct_module_tree(module: str) -> Dict[str, Any]: """ if module not in MODULE_ATTRIBUTE_MAP: raise AttributeError(f"No such module found: '{module}'.") - tree = {'attributes': {}} + tree: Dict[str, Any] = {'attributes': {}} for item in MODULE_ATTRIBUTE_MAP[module]: location = tree['attributes'] for p in item['path']: diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index ddae3dcb..7edc44a4 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -5,15 +5,16 @@ import numpy as np from collections import defaultdict from typing import ( - Any, Dict, Iterable, List, Optional, Set, Sequence, Union, Tuple + Any, Dict, List, Optional, Set, Sequence, Union, Tuple ) from pydicom.dataset import Dataset -from pydicom.datadict import tag_for_keyword +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 from pydicom.uid import ( ExplicitVRLittleEndian, ImplicitVRLittleEndian, @@ -48,6 +49,7 @@ 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__) @@ -1064,15 +1066,22 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': 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] = ( - self.StudyInstanceUID, - ref_series.SeriesInstanceUID, - 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): @@ -1080,9 +1089,9 @@ def _build_ref_instance_lut(self) -> None: for ref_series in ref_study.ReferencedSeriesSequence: for ref_ins in ref_series.ReferencedInstanceSequence: self._ref_ins_lut[ref_ins.ReferencedSOPInstanceUID] = ( - ref_study.StudyInstanceUID, - ref_series.SeriesInstanceUID, - 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( @@ -1090,6 +1099,19 @@ def _build_ref_instance_lut(self) -> None: ) 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() @@ -1105,11 +1127,19 @@ def _build_luts(self) -> None: for dim_ind in self.DimensionIndexSequence if dim_ind.DimensionIndexPointer != seg_num_tag ] - dim_index_col_data = {ptr: [] for ptr in self._dim_ind_pointers} + 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_preserved = [] + 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 @@ -1126,8 +1156,8 @@ def _build_luts(self) -> None: 'per-frames functional groups sequence and items in the ' 'dimension index sequence.' ) - for ptr, ind in zip(self._dim_ind_pointers, indices): - dim_index_col_data[ptr].append(ind) + for ptr in self._dim_ind_pointers: + dim_index_col_data[ptr].append(indices[dim_ind_positions[ptr]]) frame_source_instances = [] frame_source_frames = [] @@ -1144,7 +1174,7 @@ def _build_luts(self) -> None: ) else: locations_preserved.append( - SpatialLocationsPreservedValues.UNKNOWN + None ) if hasattr(src_im, 'ReferencedFrameNumber'): @@ -1185,12 +1215,14 @@ def _build_luts(self) -> None: # 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 ): @@ -1213,7 +1245,8 @@ def _build_luts(self) -> None: # Columns for other dimension index values self._lut_dim_ind_cols = { - ptr: i for ptr, i in enumerate(self._dim_ind_pointers, len(col_data)) + 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() @@ -1304,17 +1337,17 @@ def search_for_segments( Parameters ---------- - segment_label: Optional[str] + segment_label: Optional[str], optional Segment label filter to apply. - segmented_property_category: Optional[Union[Code, CodedConcept]] + segmented_property_category: Optional[Union[Code, CodedConcept]], optional Segmented property category filter to apply. - segmented_property_type: Optional[Union[Code, CodedConcept]] + segmented_property_type: Optional[Union[Code, CodedConcept]], optional Segmented property type filter to apply. - algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]] + algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]], optional Segmented property type filter to apply. - tracking_uid: Optional[str] + tracking_uid: Optional[str], optional Tracking unique identifier filter to apply. - tracking_id: Optional[str] + tracking_id: Optional[str], optional Tracking identifier filter to apply. Returns @@ -1322,7 +1355,7 @@ def search_for_segments( List[int] List of all segment numbers matching the provided criteria. - """ + """ # noqa: E501 filter_funcs = [] if segment_label is not None: filter_funcs.append( @@ -1361,37 +1394,109 @@ def search_for_segments( return matches - def all_tracking_ids(self) -> Set[str]: + def search_for_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[str]: """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. + + Parameters + ---------- + segmented_property_category: Optional[Union[Code, CodedConcept]], optional + Segmented property category filter to apply. + segmented_property_type: Optional[Union[Code, CodedConcept]], optional + Segmented property type filter to apply. + algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]], optional + Segmented property type filter to apply. + Returns ------- - Set[str] + List[str] All unique tracking identifiers referenced in segment descriptions - in this SEG image. + in this SEG image that match all provided filters. - """ - return { + """ # 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 for desc in self.SegmentSequence - if desc.tracking_id is not None - } + if desc.tracking_id is not None and + all(f(desc) for f in filter_funcs) + }) - def all_tracking_uids(self) -> Set[str]: + def search_for_tracking_uids( + 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[str]: """Get all unique tracking unique 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. + + Parameters + ---------- + segmented_property_category: Optional[Union[Code, CodedConcept]], optional + Segmented property category filter to apply. + segmented_property_type: Optional[Union[Code, CodedConcept]], optional + Segmented property type filter to apply. + algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]], optional + Segmented property type filter to apply. + Returns ------- - Set[str] + List[str] All unique tracking unique identifiers referenced in segment - descriptions in this SEG image. + in this SEG image that match all provided filters. - """ - return { + """ # 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_uid for desc in self.SegmentSequence - if desc.tracking_uid is not None - } + if desc.tracking_uid is not None and + all(f(desc) for f in filter_funcs) + }) - def all_segmented_property_categories(self) -> List[CodedConcept]: + def segmented_property_categories(self) -> List[CodedConcept]: """Get all unique segmented property categories in this SEG image. Returns @@ -1408,7 +1513,7 @@ def all_segmented_property_categories(self) -> List[CodedConcept]: return categories - def all_segmented_property_types(self) -> List[CodedConcept]: + def segmented_property_types(self) -> List[CodedConcept]: """Get all unique segmented property types in this SEG image. Returns @@ -1433,59 +1538,10 @@ def _get_src_uid_index(self, sop_instance_uid: str) -> int: ) return ind.item() - def list_seg_frames_for_source_frame( - self, - source_sop_instance_uid: str - ) -> List[int]: - src_ind = self._get_src_uid_index(source_sop_instance_uid) - seg_frames = np.where( - self.frame_lut[:, self.lut_src_col] == src_ind - )[0] - return seg_frames.tolist() - - def list_segments_in_source_frame( - self, - source_sop_instance_uid: str - ) -> List[int]: - src_ind = self._get_src_uid_index(source_sop_instance_uid) - segments = self.frame_lut[ - self.frame_lut[:, self.lut_src_col] == src_ind, - self.lut_seg_col - ] - return segments.tolist() - - def list_seg_frames_for_segment_number( - self, - segment_number: int - ) -> List[int]: - if segment_number not in self.segment_numbers: - raise KeyError( - 'No segment found with specified segment number' - ) - seg_frames = np.where( - self.frame_lut[:, self.lut_seg_col] == segment_number - )[0].tolist() - return seg_frames - - def list_source_frames_for_segment_number( - self, - segment_number: int - ) -> List[str]: - if segment_number not in self.segment_numbers: - raise KeyError( - 'No segment found with specified segment number' - ) - source_frame_indices = self.frame_lut[ - self.frame_lut[:, self.lut_seg_col] == segment_number, - self.lut_src_col - ] - source_frame_uids = self._source_sop_instance_uids[source_frame_indices] - return source_frame_uids.tolist() - def _get_pixels_by_seg_frame( self, seg_frames_matrix: np.ndarray, - segment_numbers: np.array, + segment_numbers: np.ndarray, combine_segments: bool = False, relabel: bool = False, ) -> np.ndarray: @@ -1591,7 +1647,7 @@ def _get_pixels_by_seg_frame( # binary values if self.segmentation_type != SegmentationTypeValues.BINARY: is_binary = np.isin( - pixel_array.unique(), + np.unique(pixel_array), np.array([0.0, 1.0]), assume_unique=True ) @@ -1620,9 +1676,89 @@ def _get_pixels_by_seg_frame( return pixel_array - def get_source_sop_instances(self): + def get_source_sop_instances(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 dimension_indices_are_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 @@ -1681,7 +1817,7 @@ def _check_indexing_with_source_frames( def get_pixels_by_source_instance( self, source_sop_instance_uids: Sequence[str], - segment_numbers: Optional[Iterable[int]] = None, + segment_numbers: Optional[Sequence[int]] = None, combine_segments: bool = False, relabel: bool = False, ignore_spatial_locations: bool = False, @@ -1689,8 +1825,8 @@ def get_pixels_by_source_instance( ) -> np.ndarray: """Get a pixel array for a list of source instances. - This is intended for retrieving segmentation masks derived from a - series of single frame source images. + 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 @@ -1735,21 +1871,21 @@ def get_pixels_by_source_instance( source_sop_instance_uids: str SOP Instance UID of the source instances to for which segmentations are requested. - segment_numbers: Optional[Sequence[int]] + segment_numbers: Optional[Sequence[int]], optional Sequence containing segment numbers to include. If unspecified, all segments are included. - combine_segments: bool + 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 + 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 + 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 @@ -1763,7 +1899,7 @@ def get_pixels_by_source_instance( 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 + 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 @@ -1861,7 +1997,7 @@ def get_pixels_by_source_frame( self, source_sop_instance_uid: str, source_frame_numbers: Sequence[int], - segment_numbers: Optional[Iterable[int]] = None, + segment_numbers: Optional[Sequence[int]] = None, combine_segments: bool = False, relabel: bool = False, ignore_spatial_locations: bool = False, @@ -1920,21 +2056,21 @@ def get_pixels_by_source_frame( source_frame_numbers: Sequence[int] A sequence of frame numbers (1-based) within the source instance for which segmentations are requested. - segment_numbers: Optional[Sequence[int]] + segment_numbers: Optional[Sequence[int]], optional Sequence containing segment numbers to include. If unspecified, all segments are included. - combine_segments: bool + 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 + 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 + 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 @@ -1948,7 +2084,7 @@ def get_pixels_by_source_frame( 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 + 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 @@ -2052,7 +2188,7 @@ def get_pixels_by_dimension_index_values( self, dimension_index_values: Sequence[Sequence[int]], dimension_index_pointers: Optional[Sequence[int]] = None, - segment_numbers: Optional[Iterable[int]] = None, + segment_numbers: Optional[Sequence[int]] = None, combine_segments: bool = False, relabel: bool = False, assert_missing_frames_are_empty: bool = False @@ -2066,9 +2202,9 @@ def get_pixels_by_dimension_index_values( 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. + 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. @@ -2104,24 +2240,39 @@ def get_pixels_by_dimension_index_values( Parameters ---------- dimension_index_values: Sequence[Sequence[int]] - TODO - dimension_index_pointers: Optional[Sequence[TODO]] - TODO - segment_numbers: Optional[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: Optional[Sequence[Union[int, pydicom.tag.BaseTag]]], 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: Optional[Sequence[int]], optional Sequence containing segment numbers to include. If unspecified, all segments are included. - combine_segments: bool + 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 + 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 + 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 @@ -2137,7 +2288,38 @@ def get_pixels_by_dimension_index_values( Pixel array representing the segmentation. See notes for full explanation. - """ + Example + ------- + + >>> import highdicom as hd + >>> from pydicom.datadict import keyword_for_tag, tag_for_keyword + >>> from pydicom import dcmread + >>> + >>> # Read a test image of a segmentation of a slide microscopy image + >>> 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.dimension_indices_are_unique(tags) # True + >>> + >>> # It is therefore possible to index using just this subset of + >>> # dimension indices + >>> pixels = seg.get_pixels_by_dimension_index_values() + + """ # noqa: E501 # Checks on validity of the inputs if segment_numbers is None: segment_numbers = list(self.segment_numbers) @@ -2169,22 +2351,21 @@ def get_pixels_by_dimension_index_values( # Create the sub-matrix of the look up table that indexes # by the dimension index values - query_cols = [ - self._dim_ind_pointers.index(ptr) + dim_ind_cols = [ + self._lut_dim_ind_cols[ptr] for ptr in dimension_index_pointers ] - lut = self._frame_lut[:, query_cols + [self._lut_seg_col]] + 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 ' + '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): - # TODO Check whether any frame exists in the LUT with these indices if len(ind_vals) != len(dimension_index_pointers): raise ValueError( 'Number of provided indices does not match the expected' @@ -2192,14 +2373,30 @@ def get_pixels_by_dimension_index_values( ) if not all(v > 0 for v in ind_vals): raise ValueError( - 'Indices are 1 based and must be greater than 1.' + 'Indices are 1-based and must be greater than 1.' ) - # Iterate over segment numbers for this source frame + # 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(ind_vals + [seg_num]) + 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 @@ -2215,24 +2412,9 @@ def get_pixels_by_dimension_index_values( # assume this frame is empty and leave the entry in seg_frames # as -1 - print(seg_frames) return self._get_pixels_by_seg_frame( seg_frames_matrix=seg_frames, segment_numbers=np.array(segment_numbers), combine_segments=combine_segments, relabel=relabel, ) - - # TODO - # Errors? - # Integrity checks: - # Tracking IDs do not have multiple segments with different - # segmented properties - # Lazy decoding of segmentation frames - # Standardize method names - # Optimize combine_segments to exploit sparse nature of array - # Implement from_dataset from DimensionIndexSequence - # Allowing indexing from multiple source series - # Index by spatial location - # Get fill pixel array with metadata - # method to access seg pffgs? 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) From 0dcab5033a812a724c8177b89fcaf53f559b1a71 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 30 May 2021 17:14:30 -0400 Subject: [PATCH 14/31] removed highdicom errors --- src/highdicom/errors.py | 11 ----------- src/highdicom/seg/sop.py | 7 +++---- 2 files changed, 3 insertions(+), 15 deletions(-) delete mode 100644 src/highdicom/errors.py diff --git a/src/highdicom/errors.py b/src/highdicom/errors.py deleted file mode 100644 index 841e7e79..00000000 --- a/src/highdicom/errors.py +++ /dev/null @@ -1,11 +0,0 @@ -"""Errors raised by highdicom processes.""" - - -class DicomAttributeError(Exception): - """DICOM standard compliance error. - - Exception indicating that a user-provided DICOM dataset is not in - compliance with the DICOM standard due to a missing attribute. - - """ - pass diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 7edc44a4..6f1ccf9e 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -32,7 +32,6 @@ PlanePositionSequence, PixelMeasuresSequence ) -from highdicom.errors import DicomAttributeError from highdicom.enum import CoordinateSystemNames from highdicom.frame import encode_frame from highdicom.seg.content import ( @@ -1003,7 +1002,7 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': seg._coordinate_system = CoordinateSystemNames.PATIENT seg._plane_orientation = plane_ori_seq.ImageOrientationPatient else: - raise DicomAttributeError( + raise AttributeError( 'Expected Plane Orientation Sequence to have either ' 'ImageOrientationSlide or ImageOrientationPatient ' 'attribute.' @@ -1011,7 +1010,7 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': for i, segment in enumerate(seg.SegmentSequence, 1): if segment.SegmentNumber != i: - raise DicomAttributeError( + raise AttributeError( 'Segments are expected to start at 1 and be consecutive ' 'integers.' ) @@ -1203,7 +1202,7 @@ def _build_luts(self) -> None: else: ref_instance_uid = frame_source_instances[0] if ref_instance_uid not in self._source_sop_instance_uids: - raise DicomAttributeError( + 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 ' From 4890c14a0d84c0c3de6ad3241346310ce087d2f3 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 30 May 2021 17:17:59 -0400 Subject: [PATCH 15/31] made module_utils private --- src/highdicom/{module_utils.py => _module_utils.py} | 0 src/highdicom/content.py | 2 +- src/highdicom/seg/content.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename src/highdicom/{module_utils.py => _module_utils.py} (100%) diff --git a/src/highdicom/module_utils.py b/src/highdicom/_module_utils.py similarity index 100% rename from src/highdicom/module_utils.py rename to src/highdicom/_module_utils.py diff --git a/src/highdicom/content.py b/src/highdicom/content.py index 872b16b9..7a8e3f3e 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -21,7 +21,7 @@ NumContentItem, TextContentItem, ) -from highdicom.module_utils import check_required_attributes +from highdicom._module_utils import check_required_attributes class AlgorithmIdentificationSequence(DataElementSequence): diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index e65206ca..80c82d4c 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -16,7 +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 +from highdicom._module_utils import check_required_attributes class SegmentDescription(Dataset): From a224a1142e2ce580a44a41cb45e7cbb4bc0ab96a Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Wed, 16 Jun 2021 19:38:18 +0000 Subject: [PATCH 16/31] Apply hackermd's suggestions from code review Co-authored-by: Markus D. Herrmann --- src/highdicom/content.py | 6 +++--- src/highdicom/seg/content.py | 4 ++-- src/highdicom/seg/sop.py | 5 ++--- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index 7a8e3f3e..04caa151 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -205,7 +205,7 @@ def from_sequence( Returns ------- - PixelMeasuresSequence: + highdicom.PixelMeasuresSequence: Plane position sequence. Raises @@ -355,7 +355,7 @@ def from_sequence( Returns ------- - PlanePositionSequence: + highdicom.PlanePositionSequence: Plane position sequence. Raises @@ -485,7 +485,7 @@ def from_sequence( Returns ------- - PlaneOrientationSequence: + highdicom.PlaneOrientationSequence: Plane orientation sequence. Raises diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index 80c82d4c..a5f669fd 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -212,12 +212,12 @@ def from_dataset(cls, dataset: Dataset) -> 'SegmentDescription': @property def segment_number(self) -> int: """int: Number of the segment.""" - return self.SegmentNumber + return int(self.SegmentNumber) @property def segment_label(self) -> str: """str: Label of the segment.""" - return self.SegmentLabel + return str(self.SegmentLabel) @property def segmented_property_category(self) -> CodedConcept: diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 6f1ccf9e..f1ecdf0b 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1385,13 +1385,12 @@ def search_for_segments( lambda desc: desc.tracking_id == tracking_id ) - matches = [ + return [ desc.segment_number for desc in self.SegmentSequence if all(f(desc) for f in filter_funcs) ] - return matches def search_for_tracking_ids( self, @@ -1675,7 +1674,7 @@ def _get_pixels_by_seg_frame( return pixel_array - def get_source_sop_instances(self) -> List[Tuple[hd_UID, hd_UID, hd_UID]]: + def get_source_instance_uids(self) -> List[Tuple[hd_UID, hd_UID, hd_UID]]: """Get UIDs for all source SOP instances referenced in the dataset. Returns From 65548a0b150284789f2118f4b7b34eb7e114fc72 Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Fri, 25 Jun 2021 23:13:10 -0400 Subject: [PATCH 17/31] Updated some docstring types to new import pattern --- src/highdicom/seg/content.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index a5f669fd..192f5795 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -165,7 +165,7 @@ def from_dataset(cls, dataset: Dataset) -> 'SegmentDescription': Returns ------- - highdicom.seg.content.SegmentDescription + highdicom.seg.SegmentDescription Segment description. """ @@ -221,7 +221,7 @@ def segment_label(self) -> str: @property def segmented_property_category(self) -> CodedConcept: - """highdicom.sr.coding.CodedConcept: + """highdicom.sr.CodedConcept: Category of the property the segment represents. """ @@ -229,7 +229,7 @@ def segmented_property_category(self) -> CodedConcept: @property def segmented_property_type(self) -> CodedConcept: - """highdicom.sr.coding.CodedConcept: + """highdicom.sr.CodedConcept: Type of the property the segment represents. """ @@ -237,7 +237,7 @@ def segmented_property_type(self) -> CodedConcept: @property def algorithm_type(self) -> SegmentAlgorithmTypeValues: - """highdicom.seg.enum.SegmentAlgorithmTypeValues: + """highdicom.seg.SegmentAlgorithmTypeValues: Type of algorithm used to create the segment. """ @@ -247,7 +247,7 @@ def algorithm_type(self) -> SegmentAlgorithmTypeValues: def algorithm_identification( self ) -> Optional[AlgorithmIdentificationSequence]: - """Optional[highdicom.content.AlgorithmIdentificationSequence] + """Optional[highdicom.AlgorithmIdentificationSequence] Information useful for identification of the algorithm, if any. """ @@ -274,7 +274,7 @@ def tracking_id(self) -> Optional[str]: @property def anatomic_regions(self) -> List[CodedConcept]: - """List[highdicom.sr.coding.CodedConcept]: + """List[highdicom.sr.CodedConcept]: List of anatomic regions into which the segment falls. May be empty. @@ -285,7 +285,7 @@ def anatomic_regions(self) -> List[CodedConcept]: @property def primary_anatomic_structures(self) -> List[CodedConcept]: - """List[highdicom.sr.coding.CodedConcept]: + """List[highdicom.sr.CodedConcept]: List of anatomic anatomic structures the segment represents. May be empty. From 4d44cbbcd1a37eb831f07c20228ae40a45f48e8e Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Fri, 25 Jun 2021 23:19:04 -0400 Subject: [PATCH 18/31] Changed Optional return markings to Union[x, None] --- src/highdicom/seg/content.py | 6 +++--- src/highdicom/seg/sop.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index 192f5795..9a993194 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -246,7 +246,7 @@ def algorithm_type(self) -> SegmentAlgorithmTypeValues: @property def algorithm_identification( self - ) -> Optional[AlgorithmIdentificationSequence]: + ) -> Union[AlgorithmIdentificationSequence, None]: """Optional[highdicom.AlgorithmIdentificationSequence] Information useful for identification of the algorithm, if any. @@ -256,7 +256,7 @@ def algorithm_identification( return None @property - def tracking_uid(self) -> Optional[str]: + def tracking_uid(self) -> Union[str, None]: """Optional[str]: Tracking unique identifier for the segment, if any. @@ -266,7 +266,7 @@ def tracking_uid(self) -> Optional[str]: return None @property - def tracking_id(self) -> Optional[str]: + def tracking_id(self) -> Union[str, None]: """Optional[str]: Tracking identifier for the segment, if any.""" if 'TrackingID' in self: return self.TrackingID diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index f1ecdf0b..8ed89ad4 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1271,7 +1271,7 @@ def segmentation_type(self) -> SegmentationTypeValues: @property def segmentation_fractional_type( self - ) -> Optional[SegmentationFractionalTypeValues]: + ) -> Union[SegmentationFractionalTypeValues, None]: """ highdicom.seg.SegmentationFractionalTypeValues: Segmentation fractional type. From 81fe054fd1126253af1a626aea66f31306269cae Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Fri, 25 Jun 2021 23:32:08 -0400 Subject: [PATCH 19/31] Rename methods --- src/highdicom/seg/sop.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 8ed89ad4..35fb5b4b 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1320,7 +1320,7 @@ def get_segment_description( ) return self.SegmentSequence[segment_number - 1] - def search_for_segments( + def get_segment_numbers( self, segment_label: Optional[str] = None, segmented_property_category: Optional[Union[Code, CodedConcept]] = None, @@ -1708,7 +1708,7 @@ def get_default_dimension_index_pointers( """ return self._dim_ind_pointers[:] - def dimension_indices_are_unique( + def are_dimension_indices_unique( self, dimension_index_pointers: Sequence[Union[int, BaseTag]] ) -> bool: @@ -2311,7 +2311,7 @@ def get_pixels_by_dimension_index_values( >>> tag_for_keyword('ColumnPositionInTotalImagePixelMatrix'), >>> tag_for_keyword('RowPositionInTotalImagePixelMatrix') >>> ] - >>> assert seg.dimension_indices_are_unique(tags) # True + >>> assert seg.are_dimension_indices_unique(tags) # True >>> >>> # It is therefore possible to index using just this subset of >>> # dimension indices From e2d8bdb710a21dc4709c6191c946647a1a9ba7a8 Mon Sep 17 00:00:00 2001 From: Christopher Bridge Date: Fri, 25 Jun 2021 23:52:14 -0400 Subject: [PATCH 20/31] Rename tracking id methods --- src/highdicom/seg/sop.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 35fb5b4b..fd03cd7e 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1392,7 +1392,7 @@ def get_segment_numbers( ] - def search_for_tracking_ids( + def get_tracking_ids( self, segmented_property_category: Optional[Union[Code, CodedConcept]] = None, segmented_property_type: Optional[Union[Code, CodedConcept]] = None, @@ -1443,7 +1443,7 @@ def search_for_tracking_ids( all(f(desc) for f in filter_funcs) }) - def search_for_tracking_uids( + def get_tracking_uids( self, segmented_property_category: Optional[Union[Code, CodedConcept]] = None, segmented_property_type: Optional[Union[Code, CodedConcept]] = None, From 4ee6173d6285e95a7cbd109be99532a7d15f3a24 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sat, 26 Jun 2021 17:15:49 -0400 Subject: [PATCH 21/31] Added dosctring for iter_segments method --- src/highdicom/seg/sop.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index fd03cd7e..25526745 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1284,6 +1284,17 @@ def segmentation_fractional_type( ) 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 From 85050f31384e6bfe81cee6245a8b8ad6b143f75a Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sat, 26 Jun 2021 18:01:33 -0400 Subject: [PATCH 22/31] Added check for correct segment numbers --- src/highdicom/seg/sop.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 25526745..d83c1d46 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1015,6 +1015,13 @@ def from_dataset(cls, dataset: Dataset) -> 'Segmentation': '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 } From 7fe1b3db698f5e669123317bc811fb87bbcfe0f7 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sat, 26 Jun 2021 23:55:50 -0400 Subject: [PATCH 23/31] Added tests for seg parsing --- src/highdicom/content.py | 2 +- src/highdicom/seg/sop.py | 90 +++----- tests/test_seg.py | 463 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 486 insertions(+), 69 deletions(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index 04caa151..e236625f 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -151,7 +151,7 @@ def parameters(self) -> Optional[Dict[str, str]]: for param in self[0].AlgorithmParameters.split(','): split = param.split('=') if len(split) != 2: - raise AttributeError('Malformed parameter string') + raise ValueError('Malformed parameter string') parameters[split[0]] = split[1] return parameters diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index d83c1d46..964f5ec7 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -14,7 +14,7 @@ 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 +from pydicom.tag import BaseTag, Tag from pydicom.uid import ( ExplicitVRLittleEndian, ImplicitVRLittleEndian, @@ -1118,7 +1118,6 @@ def _build_luts(self) -> None: index values. """ - self._build_ref_instance_lut() segnum_col_data = [] @@ -1409,68 +1408,22 @@ def get_segment_numbers( 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[str]: + ) -> 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. - Parameters - ---------- - segmented_property_category: Optional[Union[Code, CodedConcept]], optional - Segmented property category filter to apply. - segmented_property_type: Optional[Union[Code, CodedConcept]], optional - Segmented property type filter to apply. - algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]], optional - Segmented property type filter to apply. + The tracking IDs and the accompanying tracking UIDs are returned + in a list of tuples. - Returns - ------- - List[str] - All unique tracking identifiers referenced in segment descriptions - in this SEG image that match all provided filters. - - """ # 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 for desc in self.SegmentSequence - if desc.tracking_id is not None and - all(f(desc) for f in filter_funcs) - }) - - def get_tracking_uids( - 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[str]: - """Get all unique tracking unique 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. + Note that the order of the returned list is not significant and will + not in general match the order of segments. Parameters ---------- @@ -1483,9 +1436,10 @@ def get_tracking_uids( Returns ------- - List[str] - All unique tracking unique identifiers referenced in segment - in this SEG image that match all provided filters. + 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. """ # noqa: E501 filter_funcs = [] @@ -1507,11 +1461,13 @@ def get_tracking_uids( ) return list({ - desc.tracking_uid for desc in self.SegmentSequence - if desc.tracking_uid is not None and + (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. @@ -1529,6 +1485,7 @@ def segmented_property_categories(self) -> List[CodedConcept]: return categories + @property def segmented_property_types(self) -> List[CodedConcept]: """Get all unique segmented property types in this SEG image. @@ -1629,7 +1586,6 @@ def _get_pixels_by_seg_frame( ): raise ValueError( 'Segment numbers array contains invalid values.' - '-1.' ) # Initialize empty pixel array @@ -1692,7 +1648,7 @@ def _get_pixels_by_seg_frame( return pixel_array - def get_source_instance_uids(self) -> List[Tuple[hd_UID, hd_UID, hd_UID]]: + 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 @@ -2160,7 +2116,7 @@ def get_pixels_by_source_frame( # 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[:, 1]: + if src_frm_num not in lut[:, 0]: if assert_missing_frames_are_empty: continue else: @@ -2171,7 +2127,7 @@ def get_pixels_by_source_frame( "use the 'assert_missing_frames_are_empty' " 'parameter.' ) - raise KeyError(msg) + raise ValueError(msg) # Iterate over segment numbers for this source frame for c, seg_num in enumerate(segment_numbers): @@ -2353,14 +2309,18 @@ def get_pixels_by_dimension_index_values( if kw == '': kw = '' raise KeyError( - f'Tag {ptr} ({kw}) is not used as a dimension index ' - 'in this image.' + 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) @@ -2384,7 +2344,7 @@ def get_pixels_by_dimension_index_values( 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 of provided indices does not match the expected ' 'number.' ) if not all(v > 0 for v in ind_vals): diff --git a/tests/test_seg.py b/tests/test_seg.py index 5fbe1945..0cbb51dd 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 ( @@ -33,6 +34,8 @@ ) 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,378 @@ 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 + ) + + def test_from_dataset(self): + assert isinstance(self._sm_control_seg, Segmentation) + + def test_properties(self): + 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 + + 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] + ) + + class TestSegUtilities(unittest.TestCase): def setUp(self): From 77ea1458fabf2785f131644e1b89092d0c122674 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 27 Jun 2021 15:37:46 -0400 Subject: [PATCH 24/31] Full testing of seg parsing --- data/test_files/README.md | 29 +++ data/test_files/seg_image_ct_binary.dcm | Bin 0 -> 4380 bytes .../seg_image_ct_binary_fractional.dcm | Bin 0 -> 5044 bytes .../seg_image_ct_binary_overlap.dcm | Bin 0 -> 6980 bytes .../seg_image_ct_true_fractional.dcm | Bin 0 -> 5044 bytes src/highdicom/seg/sop.py | 62 ++++- tests/test_seg.py | 241 ++++++++++++++++++ 7 files changed, 329 insertions(+), 3 deletions(-) create mode 100644 data/test_files/README.md create mode 100644 data/test_files/seg_image_ct_binary.dcm create mode 100644 data/test_files/seg_image_ct_binary_fractional.dcm create mode 100644 data/test_files/seg_image_ct_binary_overlap.dcm create mode 100644 data/test_files/seg_image_ct_true_fractional.dcm 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 0000000000000000000000000000000000000000..04a19ae009880359bf10518f14865c8dab6a4897 GIT binary patch literal 4380 zcmeHKOK)3M5FR^)n(}fSQAJf%ofL!uD!1oxpDZBG18Rir$Z=XyMSvgI&V|IbuANA% z3KEMh>j<%jCBK3N3l=O8OIH0C>{lN1BP{;FqWt0~f&sC1=xq-BmcG8rEp4=U_U>SSJVdla;Ed-Wl z8|b%Gnr3M2N|v^w_Es}$ba%v_v^+~@Ck5uWX16ou%{m{-CpE zXqFZxGXLOg84}n9Yz2WM18Fk!ax&IZE9q)4Lz156dCy3rGi2xm9x3VDaJR2)&%u6I zlc;ELAJZw!udTcbzDLi2FV;G{ty*KS8;5a*Zs7~Tg3O)+U8-IZS61@NmlrGfLP=aJ zi=}d@TFI|fA(>^n#cng)hOU>Ue@+qoen7W>1kKQU(>KTC>xsVAY1Vf?k7=w}pO%S) zv0Ma&QQ`!y=J)^@{3;c=j*-3tExu6!J_Uv?o~Q-;#4R`{Zow5m#qqt%QOV8yAIHCm z#eX4=S=1xW@8xdzVLonbMfKWtQ)I_A!t-3mmd0e0SQBIyAj>qH<*b*s-pElssQ?=* zYa&ODB;T8jTBqCD>~uxc3_F9U6@_ts^V2ZwiDuYu#8Gb$bvxrPJ4edJNG1*kQ5-f$ zgI*l%)dnFe+l`uiroC_gJv-fYT*==n!LvVw}Jr z$gV*hIyt-Y>(z2GUo9-*?9UvTeetf7vb6TJiz|kdu>04uvW>-V2PYN8$YmRo;~}aD zZnq_Hsp{$A#cC5itbnG)D_qbuWDXc?XnI1;&xSrybmRY+u{HA86xi*VfW= zeBZC*4#$2jm3#Dga=yY^`J+Nbs6hAV3(OGiJ6|T0DzEL|0>{z32`cF9@SXAatz_)M zU%V%r7tbbcANOPKk`ch>@F-SharXC-LhigPp;T#YVb?hO?9Nld}}$=?&{ZYpUkKrINoIxxgh zV$F<^C4?HgJA10^~* deK2){-cH^d#P7l6O&z2TrVbwR{PvLfKLF7Wakc;e literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..90ad2b778703f1d767eed947f4967ffc92d48b41 GIT binary patch literal 5044 zcmeHK&2L*p5Fe+7ngYdfL={z4wJ8V%R9<)AzE2J)&PSDSk}3lH@!~v4Z0p&H zxKShyJ=PK8&VRyz0|yR>Bd7i)gx{>6P5MFG6osOe7ezC(J3BKwUpppxw%=lLp)>&# z)~=nT-vB2lL*;pf_%aiKYsK@VEZYhkX(=gXV7ZpE9Q-TG^DJkAPJm}_$YXmJxGe1j zQaT#^ma9NoO8KsIq^E7|c%kP?+ty*AL(j9Fz_ml=LrxBI^7rY2X4sI&&V`3+2rb{y zkZ-98nxfUKIogU_Tg|x9*%5ov@-3OAEF=>xE>()xmliixD#cQvvM#c8md<1SQFqJI zG|i28{=wa{Bwz(}g`q1$sk8K!p>1g^>DhjkB)!c0zLwg~lBP3EX{{t2A39Fp!`>AW z6&u*+bQ1fkOCJI6(rds=wf632tJME2jN5svj-H$ipDCuo{8AV;ujCzeE?)KwOJ6+j1)GmfRNz{*%s5vCM zNxWCC z$(Gnc(Bv4KGMK_zvFXr45E46+MRkQFTdpAR~$Wrq|I_u?d5n!rBRg5wP40^GRuwHg%zVU(}(KbM8u|uuv_QmkZbG&~+;_jKo)fbxbQZ zIM@&G89kKYmpO2K*OL-2@kGYFtqWC=rIQ(c9jkG#w;OTA4R2MamrVU`#!0U)dZ~Zw zbc()6CD#Er=`ziy^ap^@Vx%csD}||KsAqAZWb~j~sE}~gHsad0D7O32_IBKg+9;h( zQqawHW64eO42}xx*CCT{TUe`>ONDB20r!99(C!-#+?2U(-}W$OND007-j#LCciOnA zFj9v249#KvVVOIg?n3Qhs>8tNds0c?(>{i~!1g@X^Sr=OfeP#J7<=yk5y+XGSgDlf z3-iTm#p*ij9@gG+xl6aqorEyf-3ZD8v7JLCZ#fxH?V&W~wKU7R1dzspi)(HZifrRCs;fr9iQf~Fz&e}@a% zuJ4g)S$cbL4v^`*q4P7p;3U&O0b-szo@vUNrJGYQg!pxsqwWvaQ?1KJ>nIDAwlQhm zn}wIavD(>98j*-iX>WFt6ze9PM%3%Y?Uv|tqXd}xp7+E888#6nJ2MmPl_yNo`4t*r8uX z4sn!NGUq(^M$W#BS$U?>ai53ioPaNn9v#QBz5jT0yyqc0T>D3lj*Fk_kB-jWuTiB4 z{lQ-!G##c5Iocl*`UtDAq;}_Bc Fe*(t}mmL5A literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..a12fc40cea79170f29dfa9890483a8b87d7c6274 GIT binary patch literal 6980 zcmeHLOK%%h6h3wfm;%LdN)=UAbyE-usGZ)KJ5LrUPMlIBY)7`!k}3l1@x&QOY-{X9 z>;Q>H7cA-svFD$3!GZ+~#4qSi=q~`iGh;^~=_{>F)2d@x=iJAfbMN`i%O-lZPocP2 zO5+GCpPi(?K+}|>$^uh7sWj+H@jU68p}D5ibSb5)S(>hy`0AQ%Yi62GV9$KA$FMc* z(v*wtKo{$m#I>b-;mSfaUo49o z6|q<;*Q@zj9XqqEx72L~n~?R+)SoG$zyG1Tzk_Dz)2Umd@y%G?@)XN^ki`tVu`(s& z9ZItS3VVs=*#_eS$Ka2yd$y&dYe9;ubPrRWVvQ5kU`||vHE|8LfL9#dx{OM0_9c#w z>5KmZj#=a*j_*P@_#z+q8)37t*%H}Nj&K~?GNn@47;D1bdF*AH&2rRB(^$<>GtL03 zOEr-rKaTg7-{^EZYn`qLTR~?Ku7^R?U;8EqdZHEd{V41W!ft0Yv$Zc>BzvM@5Jo{O z8T6uXr!fdv*ml_JGwlTf$l2<)qegEd^v9C1tH0LrTLXWiE&Po}ryYnEPl_2F`4rAu z2zEk05TCE!oFS3ry0W7Sdg}>1S%-yAKRO>9o@FUhx+pRV#+16LJFe-uQaLC;cr578 zj&iA>sXcBw=5!&b`UIylak}~VewC7-*l()B9#4yjHE0mJBh6e8j6RxCDhP_9ZAV$6ZJA zZfcK9#~fkx#2%NO<$Y2bq~f^H9Y-{wlk3?+v5XSmgq|}e5){5ov4u|H9?AGaP*mg; zS7q07EL%!c;dCmltBZA!p_3_=R1f?8?SP{#IaQNhiDP#?jQRu7AKA7^uhG{d$xYC8 zxxYd^QS1zY&CReKbPzV{q(e8yj3(E~ zj(J4apb42AUHO%IrIfE17jgI3_pQEt=%&2ea169tXyTyvzjtLF3*8QGDq<-KKTGqN zzn|Ccb#|^Y(8-y$t6Q#TTBavG6Rn?8(pHAyaDT@>_tVJ!2d%=Vf>jWkXc!Nj>I_!% zoyMTM75cDDZ4hn8EzTH^G%!u@=m&nc(;5SE9Pu!57q-MDQ1%DDxykJx;5dA9_Lf>N z>;%A7#6CS9Jq9Ld(@n$FP0vQ}<(USka*_X(Cv{tSrsEh*bjy5#m2zR}YO%aj5Iy=n zo+;KUUlyxE_lzF>fE|Pm^T(J{^@H8p;P@7|Kn0y0ULB1;j>ir>l>G%)!n1~c??KF= zJp$QaxIa2)(C>0|!uu?iGKZX;!&qvj>l&^O=?Q~$EqCuB&G5N?fuCjQz2P}f9_Qk5 zZuSqn$>ZNZ(MeCJkDa5N(;P`$m*=VX>-CY=a;$Zn19Sw&IR3^tfW4yD-H!Y~gmF?` z??xl87j^xh-w!)&(d`8hax79Udi7nlHQ$LV&0dDlm2~Wb(qWu&9Y|iMsUc?_eme^) z=*Ez34P(V$&hh)P;QW6PYnOGy8e;D2-upRyY literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..dbd114f3c41a3b65314a206e3fbdccf34ef4e8c7 GIT binary patch literal 5044 zcmeHK&2L*p5FfjRngYdfL={z4wJ8V%R9^S(`|{&}>^Lbks_n>eT2e)TKVF;%iETYQ z5jTv)p~pHx-1$#9aNxiJapcs$gz}s9vq_V-G)1APsiu4~vyIt8BGm`Aw=xD3n9 zN@-i*HyjhBVVb^W`B~32O(h-AQOe4Cj-{lplnZXlb&`~ZoXjCz&{G(1D5Cb&E=(ym1?n+uda&}ou`Xff85gFRc{4rD-LVDP85INYcx^?^)7PDYEDsQ%ZS`Yj~b?lnr}V zbW{|u&*?PwSC>Bm-lNxmm+GC}&3bdN8wF8{?qY?YL8dPN7HU_-_2vB1)x~PQSQfV` zqEIQvxEu)IfkOyR`ZjMN;K z;UXxUC64bZjt>%p|Bwl;W67)oDOt<(vE*B<@kBM)C$7OcaSg6ORvh2E9F@%c|8ac6 zT>KZ}m_|S1`d;b=pXQ_HR@kU-w?sPb5uWEdwzRCYjx~Wd2VREhG*`W})p~{+x&y2) zuZj#cb-lNm^-j05(dmk?6?6vSW*9{MjW2_sCt5+j8HK$;*zLqCI|tgum?sJbVHC8+ zL@x^W>VtrV?S`#B!(K3eoSkkvs`s|S=0q~M`WvlgYtYpU&yh&;T-h@Pz5W!LY{En5o}?eU7%r~oIkt;dgK1*}k*0&0gGuAt(iiY1 zTcU)ZnJG4_F%xUmr?Ae3bvb>%N=eZCz5)5%6W#$7)EZp^)SdIt9Em5GGAu!t_kWB} z)sAs07;#=S=vC-x$M{*`^WSlntp@EKh?lgjkE!sG?Y(Z)&ThQ!?WuJ$X;!Yj$ zUg`Ej&MmQ0AB0e~06(&K6SROV)zhg`t`x3UE2U*67SNA)*8(<&B3qei`@=xDzua_ zK2JHUKP+>nuRBO{vX=Bb%d;%g@@(7lJlpkbnax_tw;d%L@ECjVJ|d9TIk8f$Ean%B zH;T1&*gdYja-~N@eW&7T<+EZ{n7-=KH#kRl=X|RnRo~se3ydRp2T;)Y(Y2WVP}3fK z#1`X@_%`tBcocGkMe9$f%rI-}{#{4Y4k^iP0z&z;O`%3h$`voM5sI^?ML!_8Q0S!_9gVnGmVb(JVfUde0lun*oJcdw%08PgML!XIx=&mX3T*cpGZy4>H;!wEyof0T|Hru&EJVQk`L{|FsR PobJE4NbbeG`f&Isf*N2U literal 0 HcmV?d00001 diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 964f5ec7..293354a2 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1517,6 +1517,7 @@ def _get_pixels_by_seg_frame( 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. @@ -1548,6 +1549,13 @@ def _get_pixels_by_seg_frame( ``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 ------- @@ -1614,10 +1622,26 @@ def _get_pixels_by_seg_frame( 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]), @@ -1793,7 +1817,8 @@ def get_pixels_by_source_instance( combine_segments: bool = False, relabel: bool = False, ignore_spatial_locations: bool = False, - assert_missing_frames_are_empty: 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. @@ -1880,6 +1905,13 @@ def get_pixels_by_source_instance( 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 ------- @@ -1898,6 +1930,11 @@ def get_pixels_by_source_instance( 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.' @@ -1963,6 +2000,7 @@ def get_pixels_by_source_instance( segment_numbers=np.array(segment_numbers), combine_segments=combine_segments, relabel=relabel, + rescale_fractional=rescale_fractional ) def get_pixels_by_source_frame( @@ -1973,7 +2011,8 @@ def get_pixels_by_source_frame( combine_segments: bool = False, relabel: bool = False, ignore_spatial_locations: bool = False, - assert_missing_frames_are_empty: 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. @@ -2065,6 +2104,13 @@ def get_pixels_by_source_frame( 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 ------- @@ -2154,6 +2200,7 @@ def get_pixels_by_source_frame( segment_numbers=np.array(segment_numbers), combine_segments=combine_segments, relabel=relabel, + rescale_fractional=rescale_fractional ) def get_pixels_by_dimension_index_values( @@ -2163,7 +2210,8 @@ def get_pixels_by_dimension_index_values( segment_numbers: Optional[Sequence[int]] = None, combine_segments: bool = False, relabel: bool = False, - assert_missing_frames_are_empty: bool = False + assert_missing_frames_are_empty: bool = False, + rescale_fractional: bool = True ): """Get a pixel array for a list of dimension index values. @@ -2253,6 +2301,13 @@ def get_pixels_by_dimension_index_values( 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 ------- @@ -2393,4 +2448,5 @@ def get_pixels_by_dimension_index_values( segment_numbers=np.array(segment_numbers), combine_segments=combine_segments, relabel=relabel, + rescale_fractional=rescale_fractional ) diff --git a/tests/test_seg.py b/tests/test_seg.py index 0cbb51dd..0d648a8f 100644 --- a/tests/test_seg.py +++ b/tests/test_seg.py @@ -31,6 +31,7 @@ SegmentAlgorithmTypeValues, SegmentsOverlapValues, SegmentationTypeValues, + SegmentationFractionalTypeValues, ) from highdicom.seg import Segmentation from highdicom.seg.utils import iter_segments @@ -1680,10 +1681,44 @@ def setUp(self): 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_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 @@ -1696,6 +1731,30 @@ def test_properties(self): 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 @@ -2042,6 +2101,188 @@ def test_get_pixels_by_dimension_index_values_invalid(self): 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): From 98b5e2ff984c33407b87ef70e2673654e3cb1e76 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Fri, 2 Jul 2021 00:02:29 +0000 Subject: [PATCH 25/31] Update src/highdicom/content.py Co-authored-by: Markus D. Herrmann --- src/highdicom/content.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index e236625f..35f536df 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -205,7 +205,7 @@ def from_sequence( Returns ------- - highdicom.PixelMeasuresSequence: + highdicom.PixelMeasuresSequence Plane position sequence. Raises From 9f449dcbe65d8b9ac8dc026f77ccbfa2a69ed453 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sat, 3 Jul 2021 21:48:15 -0400 Subject: [PATCH 26/31] Changed Optional[x] -> Union[x, None] in docstrings --- src/highdicom/content.py | 4 ++-- src/highdicom/seg/content.py | 6 +++--- src/highdicom/seg/sop.py | 26 +++++++++++++------------- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index 35f536df..790cd7b3 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -131,7 +131,7 @@ def version(self) -> str: @property def source(self) -> Optional[str]: - """Optional[str]: + """Union[str, None]: Source of the algorithm, e.g. name of the algorithm manufacturer, if any @@ -140,7 +140,7 @@ def source(self) -> Optional[str]: @property def parameters(self) -> Optional[Dict[str, str]]: - """Optional[Dict[str, str]]: + """Union[Dict[str, str], None]: Dictionary mapping algorithm parameter names to values, if any diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index 9a993194..27541644 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -247,7 +247,7 @@ def algorithm_type(self) -> SegmentAlgorithmTypeValues: def algorithm_identification( self ) -> Union[AlgorithmIdentificationSequence, None]: - """Optional[highdicom.AlgorithmIdentificationSequence] + """Union[highdicom.AlgorithmIdentificationSequence, None] Information useful for identification of the algorithm, if any. """ @@ -257,7 +257,7 @@ def algorithm_identification( @property def tracking_uid(self) -> Union[str, None]: - """Optional[str]: + """Union[str, None]: Tracking unique identifier for the segment, if any. """ @@ -267,7 +267,7 @@ def tracking_uid(self) -> Union[str, None]: @property def tracking_id(self) -> Union[str, None]: - """Optional[str]: Tracking identifier for the segment, if any.""" + """Union[str, None]: Tracking identifier for the segment, if any.""" if 'TrackingID' in self: return self.TrackingID return None diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 293354a2..70a6c9ee 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1353,17 +1353,17 @@ def get_segment_numbers( Parameters ---------- - segment_label: Optional[str], optional + segment_label: Union[str, None], optional Segment label filter to apply. - segmented_property_category: Optional[Union[Code, CodedConcept]], optional + segmented_property_category: Union[Code, CodedConcept, None], optional Segmented property category filter to apply. - segmented_property_type: Optional[Union[Code, CodedConcept]], optional + segmented_property_type: Union[Code, CodedConcept, None], optional Segmented property type filter to apply. - algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]], optional + algorithm_type: Union[SegmentAlgorithmTypeValues, str, None], optional Segmented property type filter to apply. - tracking_uid: Optional[str], optional + tracking_uid: Union[str, None], optional Tracking unique identifier filter to apply. - tracking_id: Optional[str], optional + tracking_id: Union[str, None], optional Tracking identifier filter to apply. Returns @@ -1427,11 +1427,11 @@ def get_tracking_ids( Parameters ---------- - segmented_property_category: Optional[Union[Code, CodedConcept]], optional + segmented_property_category: Union[Code, CodedConcept, None], optional Segmented property category filter to apply. - segmented_property_type: Optional[Union[Code, CodedConcept]], optional + segmented_property_type: Union[Code, CodedConcept, None], optional Segmented property type filter to apply. - algorithm_type: Optional[Union[SegmentAlgorithmTypeValues, str]], optional + algorithm_type: Union[SegmentAlgorithmTypeValues, str, None], optional Segmented property type filter to apply. Returns @@ -1868,7 +1868,7 @@ def get_pixels_by_source_instance( source_sop_instance_uids: str SOP Instance UID of the source instances to for which segmentations are requested. - segment_numbers: Optional[Sequence[int]], optional + segment_numbers: Union[Sequence[int], None], optional Sequence containing segment numbers to include. If unspecified, all segments are included. combine_segments: bool, optional @@ -2067,7 +2067,7 @@ def get_pixels_by_source_frame( source_frame_numbers: Sequence[int] A sequence of frame numbers (1-based) within the source instance for which segmentations are requested. - segment_numbers: Optional[Sequence[int]], optional + segment_numbers: Sequence[int, None], optional Sequence containing segment numbers to include. If unspecified, all segments are included. combine_segments: bool, optional @@ -2267,7 +2267,7 @@ def get_pixels_by_dimension_index_values( 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: Optional[Sequence[Union[int, pydicom.tag.BaseTag]]], optional + 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 @@ -2278,7 +2278,7 @@ def get_pixels_by_dimension_index_values( 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: Optional[Sequence[int]], optional + segment_numbers: Union[Sequence[int], None], optional Sequence containing segment numbers to include. If unspecified, all segments are included. combine_segments: bool, optional From fe168eb1e83362f21e52cf7c02b2845a3744dfd4 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sat, 3 Jul 2021 23:37:00 -0400 Subject: [PATCH 27/31] Added more examples to seg docstrings --- src/highdicom/seg/sop.py | 193 +++++++++++++++++++++++++++++++++++---- 1 file changed, 177 insertions(+), 16 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 70a6c9ee..21dc83e2 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1371,6 +1371,36 @@ def get_segment_numbers( 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: @@ -1441,6 +1471,35 @@ def get_tracking_ids( 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 highdicom.seg.utils import segread + >>> from pydicom.sr.codedict import codes + >>> + >>> 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: @@ -1919,6 +1978,36 @@ def get_pixels_by_source_instance( Pixel array representing the segmentation. See notes for full explanation. + Examples + -------- + + Read in an example from the highdicom test data: + + >>> import highdicom as hd + >>> from highdicom.seg.utils import segread + >>> + >>> 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) @@ -2118,6 +2207,68 @@ def get_pixels_by_source_frame( 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 + >>> from highdicom.seg.utils import segread + >>> + >>> 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) @@ -2315,36 +2466,46 @@ def get_pixels_by_dimension_index_values( Pixel array representing the segmentation. See notes for full explanation. - Example - ------- + 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 >>> - >>> # Read a test image of a segmentation of a slide microscopy image >>> 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 + + 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 + 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() + + 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 From ed62f994d3f6dca21c8058f2d576f1a4ec3a80ae Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sat, 3 Jul 2021 23:58:46 -0400 Subject: [PATCH 28/31] Add segread function --- src/highdicom/seg/__init__.py | 3 ++- src/highdicom/seg/sop.py | 22 +++++++++++++++++++++- tests/test_seg.py | 15 +++++++++++---- 3 files changed, 34 insertions(+), 6 deletions(-) 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/sop.py b/src/highdicom/seg/sop.py index 21dc83e2..de07d874 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -2,10 +2,11 @@ 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 + Any, Dict, List, Optional, Set, Sequence, Union, Tuple, BinaryIO ) from pydicom.dataset import Dataset @@ -25,6 +26,7 @@ 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 ( @@ -2611,3 +2613,21 @@ def get_pixels_by_dimension_index_values( 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/tests/test_seg.py b/tests/test_seg.py index 0d648a8f..613097e2 100644 --- a/tests/test_seg.py +++ b/tests/test_seg.py @@ -24,16 +24,15 @@ ) 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 @@ -1717,6 +1716,14 @@ def setUp(self): 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 From 2c44e3171cc670c1166c4df2196fd10319e030b1 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 4 Jul 2021 00:10:31 -0400 Subject: [PATCH 29/31] Fix segread paths in docstring examples --- src/highdicom/seg/sop.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index de07d874..282c98ec 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1479,10 +1479,9 @@ def get_tracking_ids( Read in an example segmentation image in the highdicom test data: >>> import highdicom as hd - >>> from highdicom.seg.utils import segread >>> from pydicom.sr.codedict import codes >>> - >>> seg = segread('data/test_files/seg_image_ct_binary_overlap.dcm') + >>> seg = hd.seg.segread('data/test_files/seg_image_ct_binary_overlap.dcm') List the tracking IDs and UIDs present in the segmentation image: @@ -1986,9 +1985,8 @@ def get_pixels_by_source_instance( Read in an example from the highdicom test data: >>> import highdicom as hd - >>> from highdicom.seg.utils import segread >>> - >>> seg = segread('data/test_files/seg_image_ct_binary.dcm') + >>> seg = hd.seg.segread('data/test_files/seg_image_ct_binary.dcm') List the source images for this segmentation: @@ -2216,9 +2214,8 @@ def get_pixels_by_source_frame( multiframe slide microscopy image: >>> import highdicom as hd - >>> from highdicom.seg.utils import segread >>> - >>> seg = segread('data/test_files/seg_image_sm_control.dcm') + >>> seg = hd.seg.segread('data/test_files/seg_image_sm_control.dcm') List the source image SOP instance UID for this segmentation: From 2c33af3b75f043aa5ba9cc092ce2e3bed4fe9759 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 4 Jul 2021 11:22:11 -0400 Subject: [PATCH 30/31] Fix exception text formatting --- src/highdicom/seg/sop.py | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index 282c98ec..6935f6e3 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -1842,27 +1842,23 @@ def _check_indexing_with_source_frames( 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. - """ + '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. - """ + '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( From 94d7c9d85ce5cf46b5230170ea8abe08782b0ff6 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Sun, 4 Jul 2021 11:32:33 -0400 Subject: [PATCH 31/31] Ensure anatomic regions and anatomic structures properties return plain list --- src/highdicom/seg/content.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/highdicom/seg/content.py b/src/highdicom/seg/content.py index 27541644..bece510f 100644 --- a/src/highdicom/seg/content.py +++ b/src/highdicom/seg/content.py @@ -281,7 +281,7 @@ def anatomic_regions(self) -> List[CodedConcept]: """ if not hasattr(self, 'AnatomicRegionSequence'): return [] - return self.AnatomicRegionSequence + return list(self.AnatomicRegionSequence) @property def primary_anatomic_structures(self) -> List[CodedConcept]: @@ -292,7 +292,7 @@ def primary_anatomic_structures(self) -> List[CodedConcept]: """ if not hasattr(self, 'PrimaryAnatomicStructureSequence'): return [] - return self.PrimaryAnatomicStructureSequence + return list(self.PrimaryAnatomicStructureSequence) class DimensionIndexSequence(DataElementSequence):