Merge pull request #319 from ImagingDataCommons/bug/nonaligned_volume

0.24.0 bugfixes

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"

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

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,

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