ENH: Multiframe support

doc/conf.py

@@ -14,7 +14,10 @@
 import sys, os
 
 #Mock unavailable packages for ReadTheDocs
-import mock
+try:
+    from unitttest import mock
+except ImportError:
+    import mock
 
 MOCK_MODULES = ['numpy', 
                 'nibabel',

src/dcmstack/convert.py

@@ -0,0 +1,74 @@
+"""Provide some common conversions for DICOM elements"""
+import re
+from datetime import date, datetime, timedelta, timezone
+
+
+def tm_to_sec(time_str):
+    '''Convert a DICOM time value (value representation of 'TM') to the number
+    of seconds past midnight.
+
+    Parameters
+    ----------
+    time_str : str
+        The DICOM time value string
+
+    Returns
+    -------
+    A floating point representing the number of seconds past midnight
+    '''
+    #Allow ACR/NEMA style format by removing any colon chars
+    time_str = time_str.replace(':', '')
+    #Only the hours portion is required
+    result = int(time_str[:2]) * 3600
+    str_len = len(time_str)
+    if str_len > 2:
+        result += int(time_str[2:4]) * 60
+    if str_len > 4:
+        result += float(time_str[4:])
+    return float(result)
+
+
+def da_to_date(date_str):
+    '''Convert a DICOM date value (value representation of 'DA') to a python `date`
+
+    Parameters
+    ----------
+    date_str : str
+        The DICOM date value string
+    '''
+    #Allow ACR/NEMA style format by removing any dot chars
+    date_str = date_str.replace('.', '')
+    return date(int(date_str[:4]), int(date_str[4:6]), int(date_str[6:]))
+
+
+def dt_to_datetime(dt_str):
+    '''Convert a DICOM datetime value (value representation of 'DT') to a python `datetime`
+
+    Parameters
+    ----------
+    dt_str : str
+        The DICOM datetime value string
+    '''
+    dt_str, suffix = re.match("([0-9\.]+)([\+\-][0-9]{4})?", dt_str).groups()
+    year = int(dt_str[:4])
+    month = day = 1
+    hour = minute = seconds = microsecs = 0
+    if len(dt_str) > 4:
+        month = int(dt_str[4:6])
+    if len(dt_str) > 6:
+        day = int(dt_str[6:8])
+    if len(dt_str) > 8:
+        hour = int(dt_str[8:10])
+    if len(dt_str) > 10:
+        minute = int(dt_str[10:12])
+    if len(dt_str) > 12:
+        seconds = int(dt_str[12:14])
+    if len(dt_str) > 14:
+        microsecs = int(float(dt_str[14:]) * 1e6)
+    tzinfo = None
+    if suffix is not None:
+        td = timedelta(hours=int(suffix[1:3]), minutes=int(suffix[3:5]))
+        if suffix[0] == '-':
+            td *= -1
+        tzinfo = timezone(td)
+    return datetime(year, month, day, hour, minute, seconds, microsecs, tzinfo=tzinfo)

src/dcmstack/dcmmeta.py

@@ -3,6 +3,7 @@
 """
 from __future__ import print_function
 
+import itertools
 import sys, re, json, warnings
 from copy import deepcopy
 try:
@@ -1229,6 +1230,45 @@ def patch_dcm_ds_is(dcm):
                     elem.value = [int(val) for val in elem.value]
 
 
+def gen_simplified_sequences(meta_dict):
+    """Get rid of useless nesting of meta data from multiframe DICOM"""
+    for k, v in meta_dict.items():
+        if isinstance(v, list):
+            if len(v) == 0:
+                continue
+            if len(v) == 1:
+                for sub_key, sub_val in v[0].items():
+                    yield sub_key, sub_val
+                continue
+        yield k, v
+
+
+def _get_fg_const_and_varying(fg_seqs):
+    varying = {}
+    for idx, fg_seq in enumerate(fg_seqs):
+        for k, v in gen_simplified_sequences(fg_seq):
+            if k not in varying:
+                if idx == 0:
+                    varying[k] = [v]
+                else:
+                    varying[k] = [None] * idx
+                    varying[k].append(v)
+            else:
+                n_vals = len(varying[k])
+                if n_vals != idx:
+                    varying[k] += [None] * (idx - n_vals)
+                varying[k].append(v)
+    const = {}
+    for k, vals in varying.items():
+        if len(vals) != len(fg_seqs):
+            vals += [None] * (len(fg_seqs) - len(vals))
+        if all(x == vals[0] for x in vals):
+            const[k] = vals[0]
+    for k in const:
+        del varying[k]
+    return const, varying
+
+
 class NiftiWrapper(object):
     '''Wraps a Nifti1Image object containing a DcmMeta header extension.
     Provides access to the meta data and the ability to split or merge the
@@ -1528,42 +1568,148 @@ def from_dicom_wrapper(klass, dcm_wrp, meta_dict=None):
         meta_dict : dict
             An optional dictionary of meta data extracted from `dcm_data`. See
             the `extract` module for generating this dict.
-
         '''
-        data = dcm_wrp.get_data()
-
-        #The Nifti patient space flips the x and y directions
+        shape = dcm_wrp.image_shape
+        n_dims = len(shape)
+        if n_dims > 4:
+            raise ValueError("5D+ multiframe not supported")
+        # Figure out any data rescaling
+        scale_factors = dcm_wrp.scale_factors
+        if dcm_wrp.is_multiframe and scale_factors.shape[1] > 1:
+            slope, inter = 1, 0
+            data = dcm_wrp.get_data()
+        else:
+            slope, inter = scale_factors[0, :]
+            data = dcm_wrp.get_unscaled_data()
+        # The Nifti patient space flips the x and y directions
         affine = np.dot(np.diag([-1., -1., 1., 1.]), dcm_wrp.affine)
-
-        #Make 2D data 3D
-        if len(data.shape) == 2:
+        n_vols = 1
+        if n_dims == 2:
             data = data.reshape(data.shape + (1,))
-
-        #Create the nifti image and set header data
+            slices_per_vol = 1
+        elif n_dims >= 3:
+            data = np.squeeze(data)
+            slices_per_vol = shape[2]
+        if n_dims == 4:
+            n_vols = shape[3]
+        # Create the nifti image and set header data
         nii_img = nb.nifti1.Nifti1Image(data, affine)
         hdr = nii_img.header
+        if (slope, inter) != (1.0, 0.0):
+            hdr.set_slope_inter(slope, inter)
         hdr.set_xyzt_units('mm', 'sec')
-        dim_info = {'freq' : None,
-                    'phase' : None,
-                    'slice' : 2
-                   }
-        phase_dir = dcm_wrp.get('InPlanePhaseEncodingDirection')
+        # Determine phase encoding direction, set dimension info in header
+        phase_info = None
+        if dcm_wrp.is_multiframe:
+            phase_info = dcm_wrp.shared.get("MRFOVGeometrySequence")
+            if phase_info is None and "MRFOVGeometrySequence" in dcm_wrp.frames[0]:
+                phase_info = [f.get("MRFOVGeometrySequence")[0] for f in dcm_wrp.frames]
+        if phase_info is None:
+            phase_info = [dcm_wrp]
+        phase_dirs = set(d.get('InPlanePhaseEncodingDirection') for d in phase_info)
+        if len(phase_dirs) > 1:
+            phase_dir = None
+        else:
+            phase_dir = phase_dirs.pop()
+        dim_info = {'freq' : None, 'phase' : None, 'slice' : 2}
         if phase_dir:
             if phase_dir == 'ROW':
-                dim_info['phase'] = 1
-                dim_info['freq'] = 0
+                dim_info['phase'], dim_info['freq'] = 1, 0
             else:
-                dim_info['phase'] = 0
-                dim_info['freq'] = 1
+                dim_info['phase'], dim_info['freq'] = 0, 1
         hdr.set_dim_info(**dim_info)
-
-        #Embed the meta data extension
+        # Create result and embed any provided meta data
         result = klass(nii_img, make_empty=True)
-
         result.meta_ext.reorient_transform = np.eye(4)
         if meta_dict:
-            result.meta_ext.get_class_dict(('global', 'const')).update(meta_dict)
-
+            if not dcm_wrp.is_multiframe:
+                global_meta = result.meta_ext.get_class_dict(('global', 'const'))
+                global_meta.update(meta_dict)
+                if dcm_wrp.is_mosaic:
+                    # For mosaic images we move a few elems that provide per-slice meta
+                    slice_meta = result.meta_ext.get_class_dict(('global', 'slices'))
+                    for key in (
+                        "SIEMENS_MR_HEADER.MosaicRefAcqTimes", 
+                        "CsaImage.MosaicRefAcqTimes",
+                        "SourceImageSequence",
+                    ):
+                        vals = global_meta.get(key)
+                        if vals is not None:
+                            slice_meta[key] = vals
+                            del global_meta[key]
+            else:
+                # Unpack and sort meta data from multiframe file that varies
+                global_meta = meta_dict.copy()
+                del global_meta["SharedFunctionalGroupsSequence"]
+                del global_meta["PerFrameFunctionalGroupsSequence"]
+                assert len(meta_dict["SharedFunctionalGroupsSequence"]) == 1
+                for k, v in gen_simplified_sequences(
+                    meta_dict["SharedFunctionalGroupsSequence"][0]
+                ):
+                    if k in global_meta:
+                        if global_meta[k] == v:
+                            continue
+                        k = f"Shared.{k}"
+                    global_meta[k] = v
+                fg_seqs = meta_dict.get("PerFrameFunctionalGroupsSequence")
+                sorted_indices = dcm_wrp.frame_order
+                fg_keys = set()
+                slice_keys = set()
+                vol_results = []
+                for vol_idx in range(n_vols):
+                    start = vol_idx * slices_per_vol
+                    end = start + slices_per_vol
+                    vol_res = _get_fg_const_and_varying(
+                        [fg_seqs[x] for x in sorted_indices[start:end]]
+                    )
+                    vol_results.append(vol_res)
+                    for res in vol_res:
+                        for k in res:
+                            fg_keys.add(k)
+                    for k in vol_res[1]:
+                        slice_keys.add(k)
+                global_slices = result.meta_ext.get_class_dict(('global', 'slices'))
+                if n_vols > 1:
+                    for vconst, vslices in vol_results:
+                        for k in fg_keys:
+                            if k not in vconst and k not in vslices:
+                                vconst[k] = None
+                        for k in slice_keys:
+                            if k in vconst:
+                                vslices[k] = [vconst[k]] * slices_per_vol
+                                del vconst[k]
+                    time_samples = result.meta_ext.get_class_dict(('time', 'samples'))
+                    time_slices = result.meta_ext.get_class_dict(('time', 'slices'))
+                    for k, first_val in vol_results[0][0].items():
+                        dest_key = k if k not in global_meta else f"PerFrame.{k}"
+                        vals = [vconst[k] for vconst, _ in vol_results]
+                        if all(v == first_val for v in vals):
+                            if k in global_meta and global_meta[k] == first_val:
+                                continue
+                            global_meta[dest_key] = first_val
+                        else:
+                            time_samples[dest_key] = vals
+                    for k, first_val in vol_results[0][1].items():
+                        dest_key = k if k not in global_meta else f"PerFrame.{k}"
+                        vals = [vslices[k] for _, vslices in vol_results]
+                        if all(v == first_val for v in vals):
+                            time_slices[dest_key] = first_val
+                        else:
+                            global_slices[dest_key] = list(itertools.chain(*vals))
+                else:
+                    vconst, vslices = vol_results[0]
+                    for k, val in vconst.items():
+                        if k in global_meta:
+                            if global_meta[k] != val:
+                                global_meta[f"PerFrame.{k}"] = val
+                        else:
+                            global_meta[k] = val
+                    for k, vals in vslices.items():
+                        if k in global_meta:
+                            global_slices[f"PerFrame.{k}"] = vals
+                        else:
+                            global_slices[k] = vals
+                result.meta_ext.get_class_dict(('global', 'const')).update(global_meta)
         return result
 
     @classmethod