Merge pull request #1485 from MPIB/enh/filter_frames_broken_stackid

ENH: Multi frame filter to handle stack with multiple orientations

Author: effigies

.gitmodules

@@ -22,3 +22,6 @@
 [submodule "nibabel-data/dcm_qa_xa30"]
 	path = nibabel-data/dcm_qa_xa30
 	url = https://github.com/neurolabusc/dcm_qa_xa30.git
+[submodule "nibabel-data/dcm_qa_xa60"]
+	path = nibabel-data/dcm_qa_xa60
+	url = https://github.com/neurolabusc/dcm_qa_xa60.git

nibabel-data/dcm_qa_xa60

@@ -597,7 +597,46 @@ def keep(self, frame) -> bool:
         return frame.MRDiffusionSequence[0].DiffusionDirectionality != 'ISOTROPIC'
 
 
-DEFAULT_FRAME_FILTERS = (FilterMultiStack(), FilterDwiIso())
+class FilterMultiOrient(FrameFilter):
+    """Filter out all but one orientation when a stack contains multiple orientations.
+
+    DICOM permits a stack to contain a group of frames with different
+    orientations, but nibabel can only process one orientation at a time.
+    This filter retains only the frames that belong to the first or a specific
+    group index.
+    """
+
+    def __init__(self, keep_group=None):
+        self._keep_group = keep_group if keep_group is not None else 0
+        self._implicit = keep_group is None
+
+    def _frame_iop(self, frame):
+        try:
+            iop = frame.PlaneOrientationSequence[0].ImageOrientationPatient
+            return tuple(round(float(v), 4) for v in iop)
+        except AttributeError:
+            return None
+
+    def applies(self, dcm_wrp) -> bool:
+        iops = dict.fromkeys(self._frame_iop(f) for f in dcm_wrp.frames)
+        if None in iops or len(iops) <= 1:
+            return False
+        if self._keep_group >= len(iops) or self._keep_group < 0:
+            raise WrapperError(f'MultiOrientation group index must be in [0,{len(iops) - 1}]')
+        self._selected = list(iops.keys())[self._keep_group]
+        if self._implicit:
+            warnings.warn(
+                'Multiple frame orientations found in a single stack; '
+                'nibabel can only process one orientation at a time. '
+                'Retaining only frames matching the first orientation.'
+            )
+        return True
+
+    def keep(self, frame) -> bool:
+        return self._frame_iop(frame) == self._selected
+
+
+DEFAULT_FRAME_FILTERS = (FilterMultiStack(), FilterDwiIso(), FilterMultiOrient())
 
 
 class MultiframeWrapper(Wrapper):

nibabel/nicom/dicomwrappers.py

@@ -12,6 +12,7 @@
 import pytest
 from numpy.testing import assert_array_almost_equal, assert_array_equal
 
+from ...openers import ImageOpener
 from ...tests.nibabel_data import get_nibabel_data, needs_nibabel_data
 from ...volumeutils import endian_codes
 from .. import dicomreaders as didr
@@ -39,6 +40,11 @@
     'dcm_qa_xa30',
     'In/20_DWI_dir80_AP/0001_1.3.12.2.1107.5.2.43.67093.2022071112140611403312307.dcm',
 )
+DATA_FILE_XA60_MULTI_ORIENT = pjoin(
+    get_nibabel_data(),
+    'dcm_qa_xa60',
+    'In/XA60/DICOM/24100413/39280000/75739321',
+)
 
 # This affine from our converted image was shown to match our image spatially
 # with an image from SPM DICOM conversion. We checked the matching with SPM
@@ -873,6 +879,57 @@ def test_data_trace(self):
         dw = didw.wrapper_from_file(DATA_FILE_SIEMENS_TRACE)
         assert dw.image_shape == (72, 72, 39)
 
+    @dicom_test
+    @needs_nibabel_data('dcm_qa_xa60')
+    def test_data_multi_orient_same_stack(self):
+        # Test that a file with 3 orientations sharing the same StackID returns
+        # only one stack/orientation
+        dw = didw.wrapper_from_file(DATA_FILE_XA60_MULTI_ORIENT)
+        assert dw.image_shape == (512, 512, 7)
+
+    @dicom_test
+    @needs_nibabel_data('dcm_qa_xa60')
+    def test_data_multi_orient_extract_all(self):
+        with ImageOpener(DATA_FILE_XA60_MULTI_ORIENT) as fobj:
+            dcm_data = pydicom.dcmread(fobj)
+        for keep_group, nslices in enumerate([7, 3, 3]):
+            dw = didw.wrapper_from_data(
+                dcm_data, frame_filters=[didw.FilterMultiOrient(keep_group=keep_group)]
+            )
+            assert dw.image_shape == (512, 512, nslices)
+
+    @dicom_test
+    def test_filter_multi_orient_unit(self):
+        iop_a = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]
+        iop_b = [0.0, 1.0, 0.0, 1.0, 0.0, 0.0]
+
+        class FakeFrame(pydicom.Dataset):
+            def __init__(self, iop=None):
+                super().__init__()
+                if iop is not None:
+                    elem = pydicom.Dataset()
+                    elem.ImageOrientationPatient = iop
+                    self.PlaneOrientationSequence = [elem]
+
+        class FakeWrp:
+            def __init__(self, frames):
+                self.frames = frames
+
+        filt = didw.FilterMultiOrient()
+        # Single orientation -> does not apply
+        assert not filt.applies(FakeWrp([FakeFrame(iop_a), FakeFrame(iop_a)]))
+        # Frame missing PlaneOrientationSequence -> does not apply
+        assert not filt.applies(FakeWrp([FakeFrame(), FakeFrame()]))
+        # Mixed orientations -> applies, keeps first group
+        wrp = FakeWrp([FakeFrame(iop_a), FakeFrame(iop_b)])
+        with pytest.warns(UserWarning, match='Multiple frame orientations'):
+            assert filt.applies(wrp)
+        assert filt.keep(wrp.frames[0])
+        assert not filt.keep(wrp.frames[1])
+        # Out-of-range group index
+        with pytest.raises(didw.WrapperError):
+            didw.FilterMultiOrient(keep_group=5).applies(wrp)
+
     @dicom_test
     @needs_nibabel_data('nitest-dicom')
     def test_data_unreadable_private_headers(self):