Bunch of multiframe fixes

Changelog

@@ -25,6 +25,33 @@ Eric Larson (EL), Demian Wassermann, Stephan Gerhard and Ross Markello (RM).
 
 References like "pr/298" refer to github pull request numbers.
 
+Upcoming release (To be determined)
+===================================
+
+New features
+------------
+
+Enhancements
+------------
+ * Ability to read data from many multiframe DICOM files that previously generated errors
+
+Bug fixes
+---------
+ * Fixed multiframe DICOM issue where data could be flipped along slice dimension relative to the
+   affine
+ * Fixed multiframe DICOM issue where ``image_position`` and the translation component in the
+   ``affine`` could be incorrect
+
+Documentation
+-------------
+
+Maintenance
+-----------
+
+API changes and deprecations
+----------------------------
+
+
 5.2.1 (Monday 26 February 2024)
 ===============================
 

nibabel/nicom/dicomwrappers.py

@@ -467,6 +467,25 @@
             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 = [f.PlanePositionSequence[0].ImagePositionPatient for f in self.frames]
+        except AttributeError:
+            try:
+                frame_ipps = [self.shared.PlanePositionSequence[0].ImagePositionPatient]
+            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 @@
         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,23 +575,85 @@
                 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)
-        # account for the 2 additional dimensions (row and column) not included
-        # in the indices
-        n_dim = frame_indices.shape[1] + 2
+            dim_seq.pop(derived_dim_idx)
+        # Determine the shape and which indices to use
+        shape = [rows, cols]
+        curr_parts = n_frames
+        frames_per_part = 1
+        del_indices = {}
+        stackpos_tag = pydicom.datadict.tag_for_keyword('InStackPositionNumber')
+        slice_dim_idx = dim_seq.index(stackpos_tag)
+        for row_idx, row in enumerate(frame_indices.T):
+            unique = np.unique(row)
+            count = len(unique)
+            if curr_parts == 1 or (count == 1 and row_idx != slice_dim_idx):
+                del_indices[row_idx] = count
+                continue
+            # Replace slice indices with order determined from slice positions along normal
+            if row_idx == slice_dim_idx:
+                if len(shape) > 2:
+                    raise WrapperError('Non-singular index precedes the slice index')
+                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")
+            elif count == n_frames:
+                if shape[-1] == 'remaining':
+                    raise WrapperError('At most one index have ambiguous size')
+                shape.append('remaining')
+                continue
+            new_parts, leftover = divmod(curr_parts, count)
+            expected = new_parts * frames_per_part
+            if leftover != 0 or any(np.count_nonzero(row == val) != expected for val in unique):
+                if row_idx == slice_dim_idx:
+                    raise WrapperError('Missing slices from multiframe')
+                del_indices[row_idx] = count
+                continue
+            if shape[-1] == 'remaining':
+                shape[-1] = new_parts
+                frames_per_part *= shape[-1]
+                new_parts = 1
+            frames_per_part *= count
+            shape.append(count)
+            curr_parts = new_parts
+        if shape[-1] == 'remaining':
+            if curr_parts > 1:
+                shape[-1] = curr_parts
+                curr_parts = 1
+            else:
+                del_indices[len(shape)] = 1
+                shape = shape[:-1]
+        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
-        if n_dim < 4:  # 3D volume
-            return rows, cols, n_frames
-        # More than 3 dimensions
-        ns_unique = [len(np.unique(row)) for row in self._frame_indices.T]
-        shape = (rows, cols) + tuple(ns_unique)
-        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
@@ -610,18 +693,11 @@
         # 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):
@@ -640,10 +716,11 @@
             raise WrapperError('No valid information for image shape')
         data = self.get_pixel_array()
         # Roll frames axis to last
-        data = data.transpose((1, 2, 0))
-        # Sort frames with first index changing fastest, last slowest
-        sorted_indices = np.lexsort(self._frame_indices.T)
-        data = data[..., sorted_indices]
+        if len(data.shape) > 2:
+            data = data.transpose((1, 2, 0))
+            # Sort frames with first index changing fastest, last slowest
+            sorted_indices = np.lexsort(self._frame_indices.T)
+            data = data[..., sorted_indices]
         data = data.reshape(shape, order='F')
         return self._scale_data(data)