BF+TST: Test and fix a bunch of multiframe fixes

Corrects issue where order of slice indices was assumed to match
the order needed to move along the direction of the slice normal,
which resulted in slice orientation flips.

Ignores indices that don't evenly divide data, and at the end will
try to combine those indices (if needed) into a single tuple index.

Author: moloney

nibabel/nicom/dicomwrappers.py

@@ -467,6 +467,25 @@ def __init__(self, dcm_data):
             self.shared = dcm_data.get('SharedFunctionalGroupsSequence')[0]
         except TypeError:
             raise WrapperError('SharedFunctionalGroupsSequence is empty.')
+        # Try to determine slice order and minimal image position patient
+        self._frame_slc_ord = self._ipp = None
+        try:
+            frame_ipps = [self.shared.PlanePositionSequence[0].ImagePositionPatient]
+        except AttributeError:
+            try:
+                frame_ipps = [f.PlanePositionSequence[0].ImagePositionPatient for f in self.frames]
+            except AttributeError:
+                frame_ipps = None
+        if frame_ipps is not None and all(ipp is not None for ipp in frame_ipps):
+            frame_ipps = [np.array(list(map(float, ipp))) for ipp in frame_ipps]
+            frame_slc_pos = [np.inner(ipp, self.slice_normal) for ipp in frame_ipps]
+            rnd_slc_pos = np.round(frame_slc_pos, 4)
+            uniq_slc_pos = np.unique(rnd_slc_pos)
+            pos_ord_map = {
+                val: order for val, order in zip(uniq_slc_pos, np.argsort(uniq_slc_pos))
+            }
+            self._frame_slc_ord = [pos_ord_map[pos] for pos in rnd_slc_pos]
+            self._ipp = frame_ipps[np.argmin(frame_slc_pos)]
         self._shape = None
 
     @cached_property
@@ -509,14 +528,16 @@ def image_shape(self):
         if hasattr(first_frame, 'get') and first_frame.get([0x18, 0x9117]):
             # DWI image may include derived isotropic, ADC or trace volume
             try:
-                anisotropic = pydicom.Sequence(
-                    frame
-                    for frame in self.frames
-                    if frame.MRDiffusionSequence[0].DiffusionDirectionality != 'ISOTROPIC'
-                )
+                aniso_frames = pydicom.Sequence()
+                aniso_slc_ord = []
+                for slc_ord, frame in zip(self._frame_slc_ord, self.frames):
+                    if frame.MRDiffusionSequence[0].DiffusionDirectionality != 'ISOTROPIC':
+                        aniso_frames.append(frame)
+                        aniso_slc_ord.append(slc_ord)
                 # Image contains DWI volumes followed by derived images; remove derived images
-                if len(anisotropic) != 0:
-                    self.frames = anisotropic
+                if len(aniso_frames) != 0:
+                    self.frames = aniso_frames
+                    self._frame_slc_ord = aniso_slc_ord
             except IndexError:
                 # Sequence tag is found but missing items!
                 raise WrapperError('Diffusion file missing information')
@@ -554,20 +575,70 @@ def image_shape(self):
                 raise WrapperError('Missing information, cannot remove indices with confidence.')
             derived_dim_idx = dim_seq.index(derived_tag)
             frame_indices = np.delete(frame_indices, derived_dim_idx, axis=1)
+        # Determine the shape and which indices to use
+        shape = [rows, cols]
+        curr_parts = n_frames
+        frames_per_part = 1
+        del_indices = {}
+        for row_idx, row in enumerate(frame_indices.T):
+            if curr_parts == 1:
+                break
+            unique = np.unique(row)
+            count = len(unique)
+            if count == 1:
+                continue
+            # Replace slice indices with order determined from slice positions along normal
+            if len(shape) == 2:
+                row = self._frame_slc_ord
+                frame_indices.T[row_idx, :] = row
+                unique = np.unique(row)
+                if len(unique) != count:
+                    raise WrapperError("Number of slice indices and positions don't match")
+            new_parts, leftover = divmod(curr_parts, count)
+            allowed_val_counts = [new_parts * frames_per_part]
+            if len(shape) > 2:
+                # Except for the slice dim, having a unique value for each frame is valid
+                allowed_val_counts.append(n_frames)
+            if leftover != 0 or any(
+                np.count_nonzero(row == val) not in allowed_val_counts for val in unique
+            ):
+                if len(shape) == 2:
+                    raise WrapperError('Missing slices from multiframe')
+                del_indices[row_idx] = count
+                continue
+            frames_per_part *= count
+            shape.append(count)
+            curr_parts = new_parts
+        if del_indices:
+            if curr_parts > 1:
+                ns_failed = [k for k, v in del_indices.items() if v != 1]
+                if len(ns_failed) > 1:
+                    # If some indices weren't used yet but we still have unaccounted for
+                    # partitions, try combining indices into single tuple and using that
+                    tup_dtype = np.dtype(','.join(['I'] * len(ns_failed)))
+                    row = [tuple(x for x in vals) for vals in frame_indices[:, ns_failed]]
+                    row = np.array(row, dtype=tup_dtype)
+            frame_indices = np.delete(frame_indices, np.array(list(del_indices.keys())), axis=1)
+            if curr_parts > 1 and len(ns_failed) > 1:
+                unique = np.unique(row, axis=0)
+                count = len(unique)
+                new_parts, rem = divmod(curr_parts, count)
+                allowed_val_counts = [new_parts * frames_per_part, n_frames]
+                if rem == 0 and all(
+                    np.count_nonzero(row == val) in allowed_val_counts for val in unique
+                ):
+                    shape.append(count)
+                    curr_parts = new_parts
+                    ord_vals = np.argsort(unique)
+                    order = {tuple(unique[i]): ord_vals[i] for i in range(count)}
+                    ord_row = np.array([order[tuple(v)] for v in row])
+                    frame_indices = np.hstack(
+                        [frame_indices, np.array(ord_row).reshape((n_frames, 1))]
+                    )
+        if curr_parts > 1:
+            raise WrapperError('Unable to determine sorting of final dimension(s)')
         # Store frame indices
         self._frame_indices = frame_indices
-        # Determine size of any extra-spatial dimensions
-        ns_unique = [len(np.unique(row)) for row in self._frame_indices.T]
-        shape = (rows, cols) + tuple(x for i, x in enumerate(ns_unique) if i == 0 or x != 1)
-        n_dim = len(shape)
-        if n_dim > 3:
-            n_vols = np.prod(shape[3:])
-            n_frames_calc = n_vols * shape[2]
-            if n_frames != n_frames_calc:
-                raise WrapperError(
-                    f'Calculated # of frames ({n_frames_calc}={n_vols}*{shape[2]}) '
-                    f'of shape {shape} does not match NumberOfFrames {n_frames}.'
-                )
         return tuple(shape)
 
     @cached_property
@@ -607,18 +678,11 @@ def voxel_sizes(self):
         # Ensure values are float rather than Decimal
         return tuple(map(float, list(pix_space) + [zs]))
 
-    @cached_property
+    @property
     def image_position(self):
-        try:
-            ipp = self.shared.PlanePositionSequence[0].ImagePositionPatient
-        except AttributeError:
-            try:
-                ipp = self.frames[0].PlanePositionSequence[0].ImagePositionPatient
-            except AttributeError:
-                raise WrapperError('Cannot get image position from dicom')
-        if ipp is None:
-            return None
-        return np.array(list(map(float, ipp)))
+        if self._ipp is None:
+            raise WrapperError('Not enough information for image_position_patient')
+        return self._ipp
 
     @cached_property
     def series_signature(self):

nibabel/nicom/tests/test_dicomwrappers.py

@@ -364,7 +364,7 @@ def test_decimal_rescale():
     assert dw.get_data().dtype != np.dtype(object)
 
 
-def fake_frames(seq_name, field_name, value_seq):
+def fake_frames(seq_name, field_name, value_seq, frame_seq=None):
     """Make fake frames for multiframe testing
 
     Parameters
@@ -375,6 +375,8 @@ def fake_frames(seq_name, field_name, value_seq):
         name of field within sequence
     value_seq : length N sequence
         sequence of values
+    frame_seq : length N list
+        previous result from this function to update
 
     Returns
     -------
@@ -386,19 +388,28 @@ def fake_frames(seq_name, field_name, value_seq):
     class Fake:
         pass
 
-    frames = []
-    for value in value_seq:
-        fake_frame = Fake()
+    if frame_seq == None:
+        frame_seq = [Fake() for _ in range(len(value_seq))]
+    for value, fake_frame in zip(value_seq, frame_seq):
         fake_element = Fake()
         setattr(fake_element, field_name, value)
         setattr(fake_frame, seq_name, [fake_element])
-        frames.append(fake_frame)
-    return frames
+    return frame_seq
 
 
-def fake_shape_dependents(div_seq, sid_seq=None, sid_dim=None):
+def fake_shape_dependents(
+    div_seq,
+    sid_seq=None,
+    sid_dim=None,
+    ipp_seq=None,
+    slice_dim=None,
+    flip_ipp_idx_corr=False,
+):
     """Make a fake dictionary of data that ``image_shape`` is dependent on.
 
+    If you are providing the ``ipp_seq`` argument, they should be generated using
+    a slice normal aligned with the z-axis (i.e. iop == (0, 1, 0, 1, 0, 0)).
+
     Parameters
     ----------
     div_seq : list of tuples
@@ -407,39 +418,86 @@ def fake_shape_dependents(div_seq, sid_seq=None, sid_dim=None):
         list of values to use for the `StackID` of each frame.
     sid_dim : int
         the index of the column in 'div_seq' to use as 'sid_seq'
+    ipp_seq : list of tuples
+        list of values to use for `ImagePositionPatient` for each frame
+    slice_dim : int
+        the index of the column in 'div_seq' corresponding to slices
+    flip_ipp_idx_corr : bool
+        generate ipp values so slice location is negatively correlated with slice index
     """
 
-    class DimIdxSeqElem:
+    class PrintBase:
+        def __repr__(self):
+            attr_strs = []
+            for attr in dir(self):
+                if attr[0].isupper():
+                    attr_strs.append(f'{attr}={getattr(self, attr)}')
+            return f"{self.__class__.__name__}({', '.join(attr_strs)})"
+
+    class DimIdxSeqElem(PrintBase):
         def __init__(self, dip=(0, 0), fgp=None):
             self.DimensionIndexPointer = dip
             if fgp is not None:
                 self.FunctionalGroupPointer = fgp
 
-    class FrmContSeqElem:
+    class FrmContSeqElem(PrintBase):
         def __init__(self, div, sid):
             self.DimensionIndexValues = div
             self.StackID = sid
 
-    class PerFrmFuncGrpSeqElem:
-        def __init__(self, div, sid):
+    class PlnPosSeqElem(PrintBase):
+        def __init__(self, ipp):
+            self.ImagePositionPatient = ipp
+
+    class PlnOrientSeqElem(PrintBase):
+        def __init__(self, iop):
+            self.ImageOrientationPatient = iop
+
+    class PerFrmFuncGrpSeqElem(PrintBase):
+        def __init__(self, div, sid, ipp, iop):
             self.FrameContentSequence = [FrmContSeqElem(div, sid)]
+            self.PlanePositionSequence = [PlnPosSeqElem(ipp)]
+            self.PlaneOrientationSequence = [PlnOrientSeqElem(iop)]
 
     # if no StackID values passed in then use the values at index 'sid_dim' in
     # the value for DimensionIndexValues for it
+    n_indices = len(div_seq[0])
     if sid_seq is None:
         if sid_dim is None:
             sid_dim = 0
         sid_seq = [div[sid_dim] for div in div_seq]
-    # create the DimensionIndexSequence
+    # Determine slice_dim and create per-slice ipp information
+    if slice_dim is None:
+        slice_dim = 1 if sid_dim == 0 else 0
     num_of_frames = len(div_seq)
-    dim_idx_seq = [DimIdxSeqElem()] * num_of_frames
+    frame_slc_indices = np.array(div_seq)[:, slice_dim]
+    uniq_slc_indices = np.unique(frame_slc_indices)
+    n_slices = len(uniq_slc_indices)
+    assert num_of_frames % n_slices == 0
+    iop_seq = [(0.0, 1.0, 0.0, 1.0, 0.0, 0.0) for _ in range(num_of_frames)]
+    if ipp_seq is None:
+        slc_locs = np.linspace(-1.0, 1.0, n_slices)
+        if flip_ipp_idx_corr:
+            slc_locs = slc_locs[::-1]
+        slc_idx_loc = {
+            div_idx: slc_locs[arr_idx] for arr_idx, div_idx in enumerate(np.sort(uniq_slc_indices))
+        }
+        ipp_seq = [(-1.0, -1.0, slc_idx_loc[idx]) for idx in frame_slc_indices]
+    else:
+        assert flip_ipp_idx_corr is False  # caller can flip it themselves
+        assert len(ipp_seq) == num_of_frames
+    # create the DimensionIndexSequence
+    dim_idx_seq = [DimIdxSeqElem()] * n_indices
     # add an entry for StackID into the DimensionIndexSequence
     if sid_dim is not None:
         sid_tag = pydicom.datadict.tag_for_keyword('StackID')
         fcs_tag = pydicom.datadict.tag_for_keyword('FrameContentSequence')
         dim_idx_seq[sid_dim] = DimIdxSeqElem(sid_tag, fcs_tag)
     # create the PerFrameFunctionalGroupsSequence
-    frames = [PerFrmFuncGrpSeqElem(div, sid) for div, sid in zip(div_seq, sid_seq)]
+    frames = [
+        PerFrmFuncGrpSeqElem(div, sid, ipp, iop)
+        for div, sid, ipp, iop in zip(div_seq, sid_seq, ipp_seq, iop_seq)
+    ]
     return {
         'NumberOfFrames': num_of_frames,
         'DimensionIndexSequence': dim_idx_seq,
@@ -480,7 +538,15 @@ def test_shape(self):
         # PerFrameFunctionalGroupsSequence does not match NumberOfFrames
         with pytest.raises(AssertionError):
             dw.image_shape
-        # check 3D shape when StackID index is 0
+        # check 2D shape with StackID index is 0
+        div_seq = ((1, 1),)
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
+        assert MFW(fake_mf).image_shape == (32, 64)
+        # Check 2D shape with extraneous extra indices
+        div_seq = ((1, 1, 2),)
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
+        assert MFW(fake_mf).image_shape == (32, 64)
+        # Check 3D shape when StackID index is 0
         div_seq = ((1, 1), (1, 2), (1, 3), (1, 4))
         fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
         assert MFW(fake_mf).image_shape == (32, 64, 4)
@@ -541,6 +607,18 @@ def test_shape(self):
         div_seq = ((1, 1, 1), (2, 1, 1), (1, 1, 2), (2, 1, 2), (1, 1, 3), (2, 1, 3))
         fake_mf.update(fake_shape_dependents(div_seq, sid_dim=1))
         assert MFW(fake_mf).image_shape == (32, 64, 2, 3)
+        # Test with combo indices, here with the last two needing to be combined into
+        # a single index corresponding to [(1, 1), (1, 1), (2, 1), (2, 1), (2, 2), (2, 2)]
+        div_seq = (
+            (1, 1, 1, 1),
+            (1, 2, 1, 1),
+            (1, 1, 2, 1),
+            (1, 2, 2, 1),
+            (1, 1, 2, 2),
+            (1, 2, 2, 2),
+        )
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
+        assert MFW(fake_mf).image_shape == (32, 64, 2, 3)
 
     def test_iop(self):
         # Test Image orient patient for multiframe
@@ -608,22 +686,30 @@ def test_image_position(self):
         with pytest.raises(didw.WrapperError):
             dw.image_position
         # Make a fake frame
-        fake_frame = fake_frames(
-            'PlanePositionSequence', 'ImagePositionPatient', [[-2.0, 3.0, 7]]
-        )[0]
-        fake_mf['SharedFunctionalGroupsSequence'] = [fake_frame]
+        iop = [0, 1, 0, 1, 0, 0]
+        frames = fake_frames('PlaneOrientationSequence', 'ImageOrientationPatient', [iop])
+        frames = fake_frames(
+            'PlanePositionSequence', 'ImagePositionPatient', [[-2.0, 3.0, 7]], frames
+        )
+        fake_mf['SharedFunctionalGroupsSequence'] = frames
         assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 7])
         fake_mf['SharedFunctionalGroupsSequence'] = [None]
         with pytest.raises(didw.WrapperError):
             MFW(fake_mf).image_position
-        fake_mf['PerFrameFunctionalGroupsSequence'] = [fake_frame]
+        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
         assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 7])
         # Check lists of Decimals work
-        fake_frame.PlanePositionSequence[0].ImagePositionPatient = [
+        frames[0].PlanePositionSequence[0].ImagePositionPatient = [
             Decimal(str(v)) for v in [-2, 3, 7]
         ]
         assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 7])
         assert MFW(fake_mf).image_position.dtype == float
+        # We should get minimum along slice normal with multiple frames
+        frames = fake_frames('PlaneOrientationSequence', 'ImageOrientationPatient', [iop] * 2)
+        ipps = [[-2.0, 3.0, 7], [-2.0, 3.0, 6]]
+        frames = fake_frames('PlanePositionSequence', 'ImagePositionPatient', ipps, frames)
+        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 6])
 
     @dicom_test
     @pytest.mark.xfail(reason='Not packaged in install', raises=FileNotFoundError)
@@ -644,7 +730,7 @@ def test_data_real(self):
         if endian_codes[data.dtype.byteorder] == '>':
             data = data.byteswap()
         dat_str = data.tobytes()
-        assert sha1(dat_str).hexdigest() == '149323269b0af92baa7508e19ca315240f77fa8c'
+        assert sha1(dat_str).hexdigest() == 'dc011bb49682fb78f3cebacf965cb65cc9daba7d'
 
     @dicom_test
     def test_slicethickness_fallback(self):
@@ -665,7 +751,7 @@ def test_data_derived_shape(self):
     def test_data_trace(self):
         # Test that a standalone trace volume is found and not dropped
         dw = didw.wrapper_from_file(DATA_FILE_SIEMENS_TRACE)
-        assert dw.image_shape == (72, 72, 39, 1)
+        assert dw.image_shape == (72, 72, 39)
 
     @dicom_test
     @needs_nibabel_data('nitest-dicom')