diff --git a/docs/volume.rst b/docs/volume.rst index 6150f998..10859b67 100644 --- a/docs/volume.rst +++ b/docs/volume.rst @@ -394,8 +394,8 @@ spatial metadata in the output object is correct. """This is a stand-in for a generic segmentation tool. We assume that the tool has certain requirements on the input array, in - this case that it has patient orientation "PRF" and a shape of (400, 400, - 2). + this case that it has patient orientation "FLP" and a shape of (2, 400, + 400). Further, we assume that the tool takes in a numpy array and returns a binary segmentation that is pixel-for-pixel aligned with its input array @@ -414,8 +414,8 @@ spatial metadata in the output object is correct. # Manipulate the original volume to give a suitable input for the tool input_volume = ( original_volume - .to_patient_orientation("PRF") - .crop_to_spatial_shape((400, 400, 2)) + .to_patient_orientation("FLP") + .crop_to_spatial_shape((2, 400, 400)) ) # Run the "complex segmentation tool" diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index f741e207..9d34c97a 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -974,8 +974,10 @@ def __init__( 'orientation.' ) - source_plane_orientation = deepcopy( - src_sfg.PlaneOrientationSequence + source_plane_orientation = ( + PlaneOrientationSequence.from_sequence( + src_sfg.PlaneOrientationSequence + ) ) else: iop = src_img.ImageOrientationPatient diff --git a/src/highdicom/volume.py b/src/highdicom/volume.py index 41b5123e..d69e3a9f 100644 --- a/src/highdicom/volume.py +++ b/src/highdicom/volume.py @@ -1097,14 +1097,14 @@ def random_permute_spatial_axes( "Argument 'axes' should contain unique values." ) - if set(axes) <= {0, 1, 2}: + if not set(axes) <= {0, 1, 2}: raise ValueError( "Argument 'axes' should contain only 0, 1, and 2." ) indices = np.random.permutation(axes).tolist() if len(indices) == 2: - missing_index = {0, 1, 2} - set(indices) + missing_index = list({0, 1, 2} - set(indices))[0] indices.insert(missing_index, missing_index) return self.permute_spatial_axes(indices) @@ -1289,7 +1289,7 @@ def random_flip_spatial(self, axes: Sequence[int] = (0, 1, 2)) -> Self: "Argument 'axes' should contain unique values." ) - if set(axes) <= {0, 1, 2}: + if not set(axes) <= {0, 1, 2}: raise ValueError( "Argument 'axes' should contain only 0, 1, and 2." ) @@ -1692,9 +1692,9 @@ def match_geometry( permute_indices = [] step_sizes = [] - for u, s in zip(self.unit_vectors(), self.spacing): + for u, s in zip(other.unit_vectors(), other.spacing): for j, (v, t) in enumerate( - zip(other.unit_vectors(), other.spacing) + zip(self.unit_vectors(), self.spacing) ): dot_product = u @ v if ( @@ -1703,7 +1703,7 @@ def match_geometry( ): permute_indices.append(j) - scale_factor = t / s + scale_factor = s / t step = int(np.round(scale_factor)) if abs(scale_factor - step) > tol: raise RuntimeError( @@ -1724,7 +1724,6 @@ def match_geometry( requires_permute = permute_indices != [0, 1, 2] if requires_permute: new_volume = self.permute_spatial_axes(permute_indices) - step_sizes = [step_sizes[i] for i in permute_indices] else: new_volume = self @@ -2789,10 +2788,13 @@ def permute_spatial_axes(self, indices: Sequence[int]) -> Self: """ new_affine = self._permute_affine(indices) - if self._array.ndim == 3: - new_array = np.transpose(self._array, indices) - else: - new_array = np.transpose(self._array, [*indices, 3]) + new_array = np.transpose( + self._array, + [ + *indices, + *[d + 3 for d in range(self.number_of_channel_dimensions)] + ] + ) return self.__class__( array=new_array, diff --git a/tests/test_seg.py b/tests/test_seg.py index bcafd696..2c1db771 100644 --- a/tests/test_seg.py +++ b/tests/test_seg.py @@ -1952,6 +1952,58 @@ def test_construction_volume(self): pp[0].ImagePositionPatient ) + def test_construction_volume_multiframe(self): + # Construction with a multiiframe source image and non-spatially + # aligned volume + arr = np.zeros((50, 50, 10), np.uint8) + arr[40:45, 34:39, 2:9] = 1 + volume = Volume( + arr, + np.eye(4), + coordinate_system="PATIENT", + frame_of_reference_uid=self._ct_multiframe.FrameOfReferenceUID, + ) + + instance = Segmentation( + [self._ct_multiframe], + volume, + SegmentationTypeValues.BINARY.value, + self._segment_descriptions, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + omit_empty_frames=False + ) + assert np.array_equal( + instance.pixel_array, + arr, + ) + + self.check_dimension_index_vals(instance) + assert instance.DimensionOrganizationType == '3D' + shared_item = instance.SharedFunctionalGroupsSequence[0] + assert len(shared_item.PixelMeasuresSequence) == 1 + pm_item = shared_item.PixelMeasuresSequence[0] + assert pm_item.PixelSpacing == [1.0, 1.0] + assert pm_item.SliceThickness == 1.0 + assert len(shared_item.PlaneOrientationSequence) == 1 + po_item = shared_item.PlaneOrientationSequence[0] + assert po_item.ImageOrientationPatient == \ + [0.0, 0.0, 1.0, 0.0, 1.0, 0.0] + for plane_item, pp in zip( + instance.PerFrameFunctionalGroupsSequence, + volume.get_plane_positions(), + ): + assert ( + plane_item.PlanePositionSequence[0].ImagePositionPatient == + pp[0].ImagePositionPatient + ) + def test_construction_volume_channels(self): # Segmentation instance from a series of single-frame CT images # with empty frames kept in, as volume with channels diff --git a/tests/test_volume.py b/tests/test_volume.py index 02c1e57f..093e5094 100644 --- a/tests/test_volume.py +++ b/tests/test_volume.py @@ -1,22 +1,40 @@ import numpy as np import pydicom +from pydicom.sr.codedict import codes from pydicom.data import get_testdata_file +from pydicom.datadict import tag_for_keyword +from pydicom.tag import BaseTag import pytest from highdicom import ( + AlgorithmIdentificationSequence, ChannelDescriptor, + PadModes, + PatientOrientationValuesBiped, + PlaneOrientationSequence, + PlanePositionSequence, + RGBColorChannels, + RGB_COLOR_CHANNEL_DESCRIPTOR, + UID, Volume, VolumeGeometry, VolumeToVolumeTransformer, - imread, get_volume_from_series, + imread, +) +from highdicom.seg import ( + Segmentation, + SegmentDescription, + SegmentAlgorithmTypeValues, + SegmentationTypeValues, ) -from highdicom.enum import PatientOrientationValuesBiped from highdicom.spatial import ( _normalize_patient_orientation, _translate_affine_matrix, ) +from tests.utils import write_and_read_dataset + def read_multiframe_ct_volume(): im = imread(get_testdata_file('eCT_Supplemental.dcm')) @@ -35,26 +53,49 @@ def read_ct_series_volume(): def test_transforms(): array = np.zeros((25, 50, 50)) + orientation = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0] volume = Volume.from_attributes( array=array, image_position=[0.0, 0.0, 0.0], - image_orientation=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + image_orientation=orientation, pixel_spacing=[1.0, 1.0], spacing_between_slices=10.0, coordinate_system="PATIENT", ) plane_positions = volume.get_plane_positions() for i, pos in enumerate(plane_positions): + assert isinstance(pos, PlanePositionSequence) assert np.array_equal( pos[0].ImagePositionPatient, [0.0, 0.0, -10.0 * i] ) + # Same thing but retrieve plane position individually + pos_2 = volume.get_plane_position(i) + assert isinstance(pos_2, PlanePositionSequence) + assert np.array_equal( + pos_2[0].ImagePositionPatient, + [0.0, 0.0, -10.0 * i] + ) + + ori = volume.get_plane_orientation() + assert isinstance(ori, PlaneOrientationSequence) + assert np.array_equal( + ori[0].ImageOrientationPatient, + orientation, + ) + indices = np.array([[1, 2, 3]]) coords = volume.map_indices_to_reference(indices) assert np.array_equal(coords, np.array([[3.0, 2.0, -10.0]])) round_trip = volume.map_reference_to_indices(coords) assert np.array_equal(round_trip, indices) + round_trip = volume.map_reference_to_indices( + coords, + check_bounds=True, + round_output=True, + ) + assert np.array_equal(round_trip, indices) index_center = volume.center_indices assert np.array_equal(index_center, [12.0, 24.5, 24.5]) index_center = volume.nearest_center_indices @@ -62,6 +103,20 @@ def test_transforms(): coord_center = volume.center_position assert np.array_equal(coord_center, [24.5, 24.5, -120]) + ras_affine = volume.get_affine("RAH") + expected = np.array( + [ + [0.0, 0.0, -1.0, 0.0], + [0.0, -1.0, 0.0, 0.0], + [-10.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + ], + ) + assert np.array_equal(ras_affine, expected) + + geom = volume.get_geometry() + assert isinstance(geom, VolumeGeometry) + @pytest.mark.parametrize( 'image_position,image_orientation,pixel_spacing,spacing_between_slices', @@ -118,6 +173,9 @@ def test_volume_from_attributes( assert volume.spatial_shape == (10, 10, 10) assert volume.channel_shape == () assert volume.channel_descriptors == () + assert volume.physical_extent == tuple( + [n * s for n, s in zip(volume.spatial_shape, volume.spacing)] + ) def test_volume_from_components(): @@ -531,6 +589,10 @@ def test_geometry_from_attributes(): assert geom.pixel_spacing == pixel_spacing assert geom.spatial_shape == (number_of_frames, rows, columns) + array = np.zeros((number_of_frames, rows, columns)) + vol = geom.with_array(array) + assert isinstance(vol, Volume) + def test_geometry_from_components(): @@ -815,3 +877,408 @@ def test_match_geometry_failure_rotation(): with pytest.raises(RuntimeError): vol.match_geometry(geometry) + + +def test_swap_axes(): + vol, _ = read_multiframe_ct_volume() + + swapped = vol.swap_spatial_axes(1, 2) + assert swapped.spatial_shape == (2, 512, 512) + + u1, u2, u3 = vol.spacing_vectors() + v1, v2, v3 = swapped.spacing_vectors() + + assert np.array_equal(u1, v1) + assert np.array_equal(u2, v3) + assert np.array_equal(u3, v2) + + +def test_random_operations(): + vol = Volume.from_components( + direction=np.eye(3), + center_position=[98.1, 78.4, 23.1], + spacing=[0.5, 0.5, 2.0], + coordinate_system="PATIENT", + array=np.random.randint(0, 10, size=(20, 20, 50)), + ) + + vol_2 = ( + vol + .random_flip_spatial() + .random_permute_spatial_axes() + .random_spatial_crop([10, 10, 10]) + ) + + matched = vol.match_geometry(vol_2) + assert np.allclose(matched.affine, vol_2.affine) + assert np.array_equal(matched.array, vol_2.array) + + +def test_random_operations_subset_axes(): + vol = Volume.from_components( + direction=np.eye(3), + center_position=[98.1, 78.4, 23.1], + spacing=[0.5, 0.5, 2.0], + coordinate_system="PATIENT", + array=np.random.randint(0, 10, size=(20, 20, 50)), + ) + + vol_2 = ( + vol + .random_flip_spatial([0, 1]) + .random_permute_spatial_axes([0, 1]) + .random_spatial_crop([10, 10, 10]) + ) + + matched = vol.match_geometry(vol_2) + assert np.allclose(matched.affine, vol_2.affine) + assert np.array_equal(matched.array, vol_2.array) + + +def test_random_operations_geometry(): + vol = VolumeGeometry.from_components( + direction=np.eye(3), + center_position=[98.1, 78.4, 23.1], + spacing=[0.5, 0.5, 2.0], + coordinate_system="PATIENT", + spatial_shape=[20, 20, 50], + ) + + vol_2 = ( + vol + .random_flip_spatial() + .random_permute_spatial_axes() + .random_spatial_crop([10, 10, 10]) + ) + + matched = vol.match_geometry(vol_2) + assert np.allclose(matched.affine, vol_2.affine) + + +@pytest.mark.parametrize('mode', list(PadModes)) +def test_pad(mode): + vol = Volume.from_components( + direction=np.eye(3), + center_position=[98.1, 78.4, 23.1], + spacing=[0.5, 0.5, 2.0], + coordinate_system="PATIENT", + array=np.random.randint(0, 10, size=(20, 20, 50)), + ) + + padded = vol.pad(10, mode=mode) + assert padded.spatial_shape == (40, 40, 70) + + padded = vol.pad([5, 10], mode=mode) + assert padded.spatial_shape == (35, 35, 65) + + padded = vol.pad([[5, 10], [10, 10], [25, 15]], mode=mode) + assert padded.spatial_shape == (35, 40, 90) + + +@pytest.mark.parametrize('mode', list(PadModes)) +@pytest.mark.parametrize('per_channel', [False, True]) +def test_pad_with_channels(mode, per_channel): + vol = Volume.from_components( + direction=np.eye(3), + center_position=[98.1, 78.4, 23.1], + spacing=[0.5, 0.5, 2.0], + coordinate_system="PATIENT", + array=np.random.randint(0, 10, size=(20, 20, 50, 3)), + channels={RGB_COLOR_CHANNEL_DESCRIPTOR: ['R', 'G', 'B']}, + ) + + padded = vol.pad(10, mode=mode, per_channel=per_channel) + assert padded.spatial_shape == (40, 40, 70) + assert padded.match_geometry(vol).geometry_equal(vol) + + padded = vol.pad([5, 10], mode=mode, per_channel=per_channel) + assert padded.spatial_shape == (35, 35, 65) + assert padded.match_geometry(vol).geometry_equal(vol) + + padded = vol.pad( + [[5, 10], [10, 10], [25, 15]], + mode=mode, + per_channel=per_channel + ) + assert padded.spatial_shape == (35, 40, 90) + assert padded.match_geometry(vol).geometry_equal(vol) + + +def test_pad_to_spatial_shape(): + vol, _ = read_multiframe_ct_volume() + + shape = (10, 600, 600) + padded = vol.pad_to_spatial_shape(shape) + assert padded.spatial_shape == shape + + assert padded.match_geometry(vol).geometry_equal(vol) + + +def test_pad_or_crop_to_spatial_shape(): + vol, _ = read_multiframe_ct_volume() + + shape = (10, 240, 240) + padded = vol.pad_or_crop_to_spatial_shape(shape) + assert padded.spatial_shape == shape + + assert padded.match_geometry(vol).geometry_equal(vol) + + +def test_normalize(): + vol, _ = read_multiframe_ct_volume() + + normed = vol.normalize_mean_std() + assert np.isclose(normed.array.mean(), 0.0) + assert np.isclose(np.std(normed.array), 1.0) + + normed = vol.normalize_min_max() + assert np.isclose(normed.array.min(), 0.0) + assert np.isclose(normed.array.max(), 1.0) + + +@pytest.mark.parametrize( + 'kw,pytype', + [ + ('DiffusionBValue', float), + ('SegmentNumber', int), + ('SegmentLabel', str), + ] +) +def test_channel_descriptors(kw, pytype): + tag = tag_for_keyword(kw) + + d1 = ChannelDescriptor(kw) + d2 = ChannelDescriptor(tag) + d3 = ChannelDescriptor(BaseTag(tag)) + d4 = ChannelDescriptor(d1) + + for d in [d1, d2, d3, d4]: + assert d.tag == tag + assert isinstance(d.tag, BaseTag) + assert d.keyword == kw + assert not d.is_custom + assert not d.is_enumerated + assert d.value_type is pytype + assert hash(d) == hash(kw) + assert str(d) == kw + assert repr(d) == kw + + +def test_channel_descriptors_custom(): + d = ChannelDescriptor('name', is_custom=True, value_type=int) + + assert d.tag is None + assert d.keyword == 'name' + assert d.is_custom + assert not d.is_enumerated + assert d.value_type is int + + +def test_channel_descriptors_enum(): + d = ChannelDescriptor( + 'name', + is_custom=True, + value_type=SegmentationTypeValues + ) + + assert d.tag is None + assert d.keyword == 'name' + assert d.is_custom + assert d.is_enumerated + assert d.value_type is SegmentationTypeValues + + +def test_multi_channels(): + + optical_path_desc = ChannelDescriptor('OpticalPathIdentifier') + path_names = ['path1', 'path2', 'path3', 'path4'] + + vol = Volume.from_components( + direction=np.eye(3), + center_position=[98.1, 78.4, 23.1], + spacing=[0.5, 0.5, 2.0], + coordinate_system="PATIENT", + array=np.random.randint(0, 10, size=(20, 20, 50, 3, 4)), + channels={ + RGB_COLOR_CHANNEL_DESCRIPTOR: ['R', 'G', 'B'], + optical_path_desc: path_names + }, + ) + + assert vol.shape == (20, 20, 50, 3, 4) + assert vol.spatial_shape == (20, 20, 50) + assert vol.channel_shape == (3, 4) + assert vol.number_of_channel_dimensions == 2 + assert vol.channel_descriptors == ( + RGB_COLOR_CHANNEL_DESCRIPTOR, + optical_path_desc, + ) + assert vol.get_channel_values(optical_path_desc) == path_names + assert vol.get_channel_values('OpticalPathIdentifier') == path_names + assert vol.get_channel_values(optical_path_desc.tag) == path_names + rgb = [RGBColorChannels.R, RGBColorChannels.G, RGBColorChannels.B] + assert vol.get_channel_values(RGB_COLOR_CHANNEL_DESCRIPTOR) == rgb + assert vol.get_channel_values( + RGB_COLOR_CHANNEL_DESCRIPTOR.keyword + ) == rgb + + permuted = vol.permute_channel_axes( + [optical_path_desc, RGB_COLOR_CHANNEL_DESCRIPTOR] + ) + assert permuted.shape == (20, 20, 50, 4, 3) + assert permuted.channel_shape == (4, 3) + assert permuted.channel_descriptors == ( + optical_path_desc, + RGB_COLOR_CHANNEL_DESCRIPTOR, + ) + + permuted = vol.permute_channel_axes_by_index([1, 0]) + assert permuted.shape == (20, 20, 50, 4, 3) + assert permuted.channel_shape == (4, 3) + assert permuted.channel_descriptors == ( + optical_path_desc, + RGB_COLOR_CHANNEL_DESCRIPTOR, + ) + + path1 = vol.get_channel(OpticalPathIdentifier='path1') + assert path1.channel_shape == (3, ) + assert path1.number_of_channel_dimensions == 1 + assert path1.channel_descriptors == (RGB_COLOR_CHANNEL_DESCRIPTOR, ) + assert np.array_equal(path1.array, vol.array[:, :, :, :, 0]) + + path2 = vol.get_channel(OpticalPathIdentifier='path2', keepdims=True) + assert path2.channel_shape == (3, 1) + assert path2.number_of_channel_dimensions == 2 + assert path2.channel_descriptors == ( + RGB_COLOR_CHANNEL_DESCRIPTOR, + optical_path_desc + ) + assert np.array_equal(path2.array, vol.array[:, :, :, :, 1:2]) + squeezed = path2.squeeze_channel() + assert squeezed.channel_shape == (3, ) + assert squeezed.number_of_channel_dimensions == 1 + assert squeezed.channel_descriptors == (RGB_COLOR_CHANNEL_DESCRIPTOR, ) + assert np.array_equal(squeezed.array, vol.array[:, :, :, :, 1]) + + red_channel = vol.get_channel(RGBColorChannel='R') + assert red_channel.channel_shape == (4, ) + assert red_channel.number_of_channel_dimensions == 1 + assert red_channel.channel_descriptors == (optical_path_desc, ) + assert np.array_equal(red_channel.array, vol.array[:, :, :, 0, :]) + + red_channel = vol.get_channel(RGBColorChannel=RGBColorChannels.R) + assert red_channel.channel_shape == (4, ) + assert red_channel.number_of_channel_dimensions == 1 + assert red_channel.channel_descriptors == (optical_path_desc, ) + assert np.array_equal(red_channel.array, vol.array[:, :, :, 0, :]) + + +def test_match_geometry_segmentation(): + # Test that creates a segmentation from a manipulated volume and ensures + # that the result can be matched back to the input image + + # Load an enhanced (multiframe) CT image + im = imread(get_testdata_file('eCT_Supplemental.dcm')) + + # Load the input volume + original_volume = im.get_volume() + + # Manipulate the original volume + input_volume = ( + original_volume + .to_patient_orientation("PRF") + .crop_to_spatial_shape((400, 400, 2)) + ) + + # Form a segmentation from the manpulated array + seg_array = input_volume.array > 0 + + # Since the seg array shares its geometry with the inupt array, we can + # combine the two to create a volume of the segmentation array + seg_volume = input_volume.with_array(seg_array) + + algorithm_identification = AlgorithmIdentificationSequence( + name='Complex Segmentation Tool', + version='v1.0', + family=codes.cid7162.ArtificialIntelligence + ) + + # metadata needed for a segmentation + brain_description = SegmentDescription( + segment_number=1, + segment_label='brain', + segmented_property_category=codes.SCT.Organ, + segmented_property_type=codes.SCT.Brain, + algorithm_type=SegmentAlgorithmTypeValues.AUTOMATIC, + algorithm_identification=algorithm_identification, + ) + + # Use the segmentation volume as input to create a DICOM Segmentation + seg_dataset = Segmentation( + pixel_array=seg_volume, + source_images=[im], + segmentation_type=SegmentationTypeValues.LABELMAP, + segment_descriptions=[brain_description], + series_instance_uid=UID(), + series_number=1, + sop_instance_uid=UID(), + instance_number=1, + manufacturer='Complex Segmentations Plc.', + manufacturer_model_name='Complex Segmentation Tool', + software_versions='0.0.1', + device_serial_number='1234567890', + omit_empty_frames=False, + ) + + seg_dataset = write_and_read_dataset(seg_dataset) + seg_dataset = Segmentation.from_dataset(seg_dataset, copy=False) + + # The reconstructed volume should be the same as the input volume, but may + # have a different handedness + seg_volume_recon = ( + seg_dataset + .get_volume(combine_segments=True) + .ensure_handedness(seg_volume.handedness, flip_axis=0) + ) + + assert np.array_equal(seg_volume_recon.affine, seg_volume.affine) + assert np.array_equal(seg_volume_recon.array, seg_volume.array) + + # Alternatively, it may be desirable to match the geometry of the output + # segmentation image to that of the input image + seg_volume_matched = seg_volume.match_geometry(original_volume) + + assert np.array_equal(original_volume.affine, seg_volume_matched.affine) + + # Use the segmentation volume as input to create a DICOM Segmentation + seg_dataset_matched = Segmentation( + pixel_array=seg_volume_matched, + source_images=[im], + segmentation_type=SegmentationTypeValues.LABELMAP, + segment_descriptions=[brain_description], + series_instance_uid=UID(), + series_number=1, + sop_instance_uid=UID(), + instance_number=1, + manufacturer='Complex Segmentations Plc.', + manufacturer_model_name='Complex Segmentation Tool', + software_versions='0.0.1', + device_serial_number='1234567890', + ) + seg_dataset_matched = write_and_read_dataset(seg_dataset_matched) + seg_dataset_matched = Segmentation.from_dataset( + seg_dataset_matched, + copy=False + ) + seg_vol_from_matched_dataset = seg_dataset_matched.get_volume( + combine_segments=True + ) + + assert np.array_equal( + seg_vol_from_matched_dataset.affine, + original_volume.affine + ) + assert np.array_equal( + seg_vol_from_matched_dataset.array, + seg_volume_matched.array + )