Merge pull request #439 from embaker/bf-stackid-assumption

MRG: Remove assumption about StackID index

In determining the shape of a multi-frame DICOM, the first value in the DimensionIndexValues was taken to be the StackID. The DICOM standard does not appear to guarantee this to always be true. Replaced the assumption by using the DimensionIndexSequence to determine the position of the StackID in the DimensionIndexValues.

Author: matthew-brett

nibabel/nicom/dicomwrappers.py

@@ -21,7 +21,7 @@
 from .dwiparams import B2q, nearest_pos_semi_def, q2bg
 from ..openers import ImageOpener
 from ..onetime import setattr_on_read as one_time
-
+from ..pydicom_compat import pydicom
 
 class WrapperError(Exception):
     pass
@@ -461,7 +461,29 @@ def __init__(self, dcm_data):
 
     @one_time
     def image_shape(self):
-        """The array shape as it will be returned by ``get_data()``"""
+        """The array shape as it will be returned by ``get_data()``
+
+        The shape is determined by the *Rows* DICOM attribute, *Columns*
+        DICOM attribute, and the set of frame indices given by the
+        *FrameContentSequence[0].DimensionIndexValues* DICOM attribute of each
+        element in the *PerFrameFunctionalGroupsSequence*.  The first two
+        axes of the returned shape correspond to the rows, and columns
+        respectively. The remaining axes correspond to those of the frame
+        indices with order preserved.
+
+        What each axis in the frame indices refers to is given by the
+        corresponding entry in the *DimensionIndexSequence* DICOM attribute.
+        **WARNING**: Any axis refering to the *StackID* DICOM attribute will
+        have been removed from the frame indices in determining the shape. This
+        is because only a file containing a single stack is currently allowed by
+        this wrapper.
+
+        References
+        ----------
+
+        * C.7.6.16 Multi-Frame Functional Groups Module: http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_C.7.6.16
+        * C.7.6.17 Multi-Frame Dimension Module: http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_C.7.6.17
+        """
         rows, cols = self.get('Rows'), self.get('Columns')
         if None in (rows, cols):
             raise WrapperError("Rows and/or Columns are empty.")
@@ -471,13 +493,25 @@ def image_shape(self):
         frame_indices = np.array(
             [frame.FrameContentSequence[0].DimensionIndexValues
              for frame in self.frames])
-        n_dim = frame_indices.shape[1] + 1
-        # Check there is only one multiframe stack index
-        if np.any(np.diff(frame_indices[:, 0])):
-            raise WrapperError("File contains more than one StackID. Cannot "
-                               "handle multi-stack files")
+        # Check that there is only one multiframe stack index
+        stack_ids = set(frame.FrameContentSequence[0].StackID
+                        for frame in self.frames)
+        if len(stack_ids) > 1:
+            raise WrapperError("File contains more than one StackID. "
+                               "Cannot handle multi-stack files")
+        # Determine if one of the dimension indices refers to the stack id
+        dim_seq = [dim.DimensionIndexPointer
+                   for dim in self.get('DimensionIndexSequence')]
+        stackid_tag = pydicom.datadict.tag_for_name('StackID')
+        # remove the stack id axis if present
+        if stackid_tag in dim_seq:
+            stackid_dim_idx = dim_seq.index(stackid_tag)
+            frame_indices = np.delete(frame_indices, stackid_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
         # Store frame indices
-        self._frame_indices = frame_indices[:, 1:]
+        self._frame_indices = frame_indices
         if n_dim < 4:  # 3D volume
             return rows, cols, n_frames
         # More than 3 dimensions

nibabel/nicom/tests/test_dicomwrappers.py

@@ -5,7 +5,7 @@
 import gzip
 from hashlib import sha1
 from decimal import Decimal
-from copy import copy
+from copy import copy
 
 import numpy as np
 
@@ -379,8 +379,54 @@ class Fake(object):
         setattr(fake_frame, seq_name, [fake_element])
         frames.append(fake_frame)
     return frames
-
-
+
+
+def fake_shape_dependents(div_seq, sid_seq=None, sid_dim=None):
+    """ Make a fake dictionary of data that ``image_shape`` is dependent on.
+
+    Parameters
+    ----------
+    div_seq : list of tuples
+        list of values to use for the `DimensionIndexValues` of each frame.
+    sid_seq : list of int
+        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'
+    """
+    class DimIdxSeqElem(object):
+        def __init__(self, dip=(0, 0), fgp=None):
+            self.DimensionIndexPointer = dip
+            if fgp is not None:
+                self.FunctionalGroupPointer = fgp
+    class FrmContSeqElem(object):
+        def __init__(self, div, sid):
+            self.DimensionIndexValues = div
+            self.StackID = sid
+    class PerFrmFuncGrpSeqElem(object):
+        def __init__(self, div, sid):
+            self.FrameContentSequence = [FrmContSeqElem(div, sid)]
+    # if no StackID values passed in then use the values at index 'sid_dim' in
+    # the value for DimensionIndexValues for it
+    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
+    num_of_frames = len(div_seq)
+    dim_idx_seq = [DimIdxSeqElem()] * num_of_frames
+    # add an entry for StackID into the DimensionIndexSequence
+    if sid_dim is not None:
+        sid_tag = pydicom.datadict.tag_for_name('StackID')
+        fcs_tag = pydicom.datadict.tag_for_name('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)]
+    return {'NumberOfFrames' : num_of_frames,
+            'DimensionIndexSequence' : dim_idx_seq,
+            'PerFrameFunctionalGroupsSequence' : frames}
+
+
 class TestMultiFrameWrapper(TestCase):
     # Test MultiframeWrapper
     MINIMAL_MF = {
@@ -406,47 +452,66 @@ def test_shape(self):
         assert_raises(AssertionError, getattr, dw, 'image_shape')
         fake_mf['NumberOfFrames'] = 4
         # PerFrameFunctionalGroupsSequence does not match NumberOfFrames
-        assert_raises(AssertionError, getattr, dw, 'image_shape')
-        # Make some fake frame data for 3D
-
-        def my_fake_frames(div_seq):
-            return fake_frames('FrameContentSequence',
-                               'DimensionIndexValues',
-                               div_seq)
-        div_seq = ((1, 1), (1, 2), (1, 3), (1, 4))
-        frames = my_fake_frames(div_seq)
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        assert_raises(AssertionError, getattr, dw, 'image_shape')
+        # 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_equal(MFW(fake_mf).image_shape, (32, 64, 4))
-        # Check stack number matching
+        # Check stack number matching when StackID index is 0
         div_seq = ((1, 1), (1, 2), (1, 3), (2, 4))
-        frames = my_fake_frames(div_seq)
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
         assert_raises(didw.WrapperError, getattr, MFW(fake_mf), 'image_shape')
-        # Make some fake frame data for 4D
-        fake_mf['NumberOfFrames'] = 6
+        # Make some fake frame data for 4D when StackID index is 0
         div_seq = ((1, 1, 1), (1, 2, 1), (1, 1, 2), (1, 2, 2),
-                   (1, 1, 3), (1, 2, 3))
-        frames = my_fake_frames(div_seq)
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+                (1, 1, 3), (1, 2, 3))
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
         assert_equal(MFW(fake_mf).image_shape, (32, 64, 2, 3))
-        # Check stack number matching for 4D
+        # Check stack number matching for 4D when StackID index is 0
         div_seq = ((1, 1, 1), (1, 2, 1), (1, 1, 2), (1, 2, 2),
-                   (1, 1, 3), (2, 2, 3))
-        frames = my_fake_frames(div_seq)
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+                (1, 1, 3), (2, 2, 3))
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
         assert_raises(didw.WrapperError, getattr, MFW(fake_mf), 'image_shape')
-        # Check indices can be non-contiguous
+        # Check indices can be non-contiguous when StackID index is 0
         div_seq = ((1, 1, 1), (1, 2, 1), (1, 1, 3), (1, 2, 3))
-        frames = my_fake_frames(div_seq)
-        fake_mf['NumberOfFrames'] = 4
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
         assert_equal(MFW(fake_mf).image_shape, (32, 64, 2, 2))
-        # Check indices can include zero
+        # Check indices can include zero when StackID index is 0
         div_seq = ((1, 1, 0), (1, 2, 0), (1, 1, 3), (1, 2, 3))
-        frames = my_fake_frames(div_seq)
-        fake_mf['NumberOfFrames'] = 4
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0))
         assert_equal(MFW(fake_mf).image_shape, (32, 64, 2, 2))
+        # check 3D shape when there is no StackID index
+        div_seq = ((1,), (2,), (3,), (4,))
+        sid_seq = (1, 1, 1, 1)
+        fake_mf.update(fake_shape_dependents(div_seq, sid_seq=sid_seq))
+        assert_equal(MFW(fake_mf).image_shape, (32, 64, 4))
+        # check 3D stack number matching when there is no StackID index
+        div_seq = ((1,), (2,), (3,), (4,))
+        sid_seq = (1, 1, 1, 2)
+        fake_mf.update(fake_shape_dependents(div_seq, sid_seq=sid_seq))
+        assert_raises(didw.WrapperError, getattr, MFW(fake_mf), 'image_shape')
+        # check 4D shape when there is no StackID index
+        div_seq = ((1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3))
+        sid_seq = (1, 1, 1, 1, 1, 1)
+        fake_mf.update(fake_shape_dependents(div_seq, sid_seq=sid_seq))
+        assert_equal(MFW(fake_mf).image_shape, (32, 64, 2, 3))
+        # check 4D stack number matching when there is no StackID index
+        div_seq = ((1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3))
+        sid_seq = (1, 1, 1, 1, 1, 2)
+        fake_mf.update(fake_shape_dependents(div_seq, sid_seq=sid_seq))
+        assert_raises(didw.WrapperError, getattr, MFW(fake_mf), 'image_shape')
+        # check 3D shape when StackID index is 1
+        div_seq = ((1, 1), (2, 1), (3, 1), (4, 1))
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=1))
+        assert_equal(MFW(fake_mf).image_shape, (32, 64, 4))
+        # Check stack number matching when StackID index is 1
+        div_seq = ((1, 1), (2, 1), (3, 2), (4, 1))
+        fake_mf.update(fake_shape_dependents(div_seq, sid_dim=1))
+        assert_raises(didw.WrapperError, getattr, MFW(fake_mf), 'image_shape')
+        # Make some fake frame data for 4D when StackID index is 1
+        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_equal(MFW(fake_mf).image_shape, (32, 64, 2, 3))
 
     def test_iop(self):
         # Test Image orient patient for multiframe
@@ -562,12 +627,9 @@ def test_data_fake(self):
         assert_raises(didw.WrapperError, dw.get_data)
         # Make shape and indices
         fake_mf['Rows'] = 2
-        fake_mf['Columns'] = 3
-        fake_mf['NumberOfFrames'] = 4
-        frames = fake_frames('FrameContentSequence',
-                             'DimensionIndexValues',
-                             ((1, 1), (1, 2), (1, 3), (1, 4)))
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        fake_mf['Columns'] = 3
+        dim_idxs = ((1, 1), (1, 2), (1, 3), (1, 4))
+        fake_mf.update(fake_shape_dependents(dim_idxs, sid_dim=0))
         assert_equal(MFW(fake_mf).image_shape, (2, 3, 4))
         # Still fails - no data
         assert_raises(didw.WrapperError, dw.get_data)
@@ -582,11 +644,9 @@ def test_data_fake(self):
         fake_mf['RescaleSlope'] = 2.0
         fake_mf['RescaleIntercept'] = -1
         assert_array_equal(MFW(fake_mf).get_data(), data * 2.0 - 1)
-        # Check slice sorting
-        frames = fake_frames('FrameContentSequence',
-                             'DimensionIndexValues',
-                             ((1, 4), (1, 2), (1, 3), (1, 1)))
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
+        # Check slice sorting
+        dim_idxs = ((1, 4), (1, 2), (1, 3), (1, 1))
+        fake_mf.update(fake_shape_dependents(dim_idxs, sid_dim=0))
         sorted_data = data[..., [3, 1, 2, 0]]
         fake_mf['pixel_array'] = np.rollaxis(sorted_data, 2)
         assert_array_equal(MFW(fake_mf).get_data(), data * 2.0 - 1)
@@ -608,11 +668,7 @@ def test_data_fake(self):
             [1, 2, 1, 2],
             [1, 3, 1, 2],
             [1, 1, 1, 2]]
-        frames = fake_frames('FrameContentSequence',
-                             'DimensionIndexValues',
-                             dim_idxs)
-        fake_mf['PerFrameFunctionalGroupsSequence'] = frames
-        fake_mf['NumberOfFrames'] = len(frames)
+        fake_mf.update(fake_shape_dependents(dim_idxs, sid_dim=0))
         shape = (2, 3, 4, 2, 2)
         data = np.arange(np.prod(shape)).reshape(shape)
         sorted_data = data.reshape(shape[:2] + (-1,), order='F')