Merge pull request #427 from grlee77/parrec_labeling

MRG: add a function to retrieve .PAR data labels for the 4th dimension

This was pulled out of the previously submitted #278 to make reviewing easier.  This complements the functionality in the recently merged #409.

The primary enhancement is a new method in the PARREC header called `get_volume_labels` that lists which .PAR fields the data stored along the 4th dimension correspond to.  (i.e. all data >4D is still collapsed into 4 dimensions exactly as before, but labels can be retrieved to make the data ordering unambiguous).

The script interface to parrec can write these labels out as a CSV file with one row per volume.

Author: matthew-brett

bin/parrec2nii

@@ -8,6 +8,7 @@ import numpy as np
 import numpy.linalg as npl
 import sys
 import os
+import csv
 import nibabel
 import nibabel.parrec as pr
 from nibabel.parrec import one_line
@@ -66,6 +67,18 @@ def get_opt_parser():
                    """The magnetic field strength of the recording, only needed
                    for --dwell-time. The field strength must be supplied
                    because it is not encoded in the PAR/REC format.""")))
+    p.add_option(
+        Option("-i", "--volume-info", action="store_true", dest="vol_info",
+               default=False,
+               help=one_line(
+                   """Export .PAR volume labels corresponding to the fourth
+                   dimension of the data.  The dimension info will be stored in
+                   CSV format with the first row containing dimension labels
+                   and the subsequent rows (one per volume), the corresponding
+                   indices.  Only labels that vary along the 4th dimension are
+                   exported (e.g. for a single volume structural scan there
+                   are no dynamic labels and no output file will be created).
+                   """)))
     p.add_option(
         Option("--origin", action="store", dest="origin", default="scanner",
                help=one_line(
@@ -262,6 +275,17 @@ def proc_file(infile, opts):
                         fid.write('%s ' % val)
                     fid.write('\n')
 
+    # export data labels varying along the 4th dimensions if requested
+    if opts.vol_info:
+        labels = pr_img.header.get_volume_labels()
+        if len(labels) > 0:
+            vol_keys = list(labels.keys())
+            with open(basefilename + '.ordering.csv', 'w') as csvfile:
+                csvwriter = csv.writer(csvfile, delimiter=',')
+                csvwriter.writerow(vol_keys)
+                for vals in zip(*[labels[k] for k in vol_keys]):
+                    csvwriter.writerow(vals)
+
     # write out dwell time if requested
     if opts.dwell_time:
         try:

nibabel/parrec.py

@@ -95,6 +95,7 @@
 import re
 from io import StringIO
 from locale import getpreferredencoding
+from collections import OrderedDict
 
 from .keywordonly import kw_only_meth
 from .spatialimages import SpatialHeader, SpatialImage
@@ -1123,6 +1124,57 @@ def get_sorted_slice_indices(self):
         n_used = np.prod(self.get_data_shape()[2:])
         return sort_order[:n_used]
 
+    def get_volume_labels(self):
+        """ Dynamic labels corresponding to the final data dimension(s).
+
+        This is useful for custom data sorting.  A subset of the info in
+        ``self.image_defs`` is returned in an order that matches the final
+        data dimension(s).  Only labels that have more than one unique value
+        across the dataset will be returned.
+
+        Returns
+        -------
+        sort_info : dict
+            Each key corresponds to volume labels for a dynamically varying
+            sequence dimension.  The ordering of the labels matches the volume
+            ordering determined via ``self.get_sorted_slice_indices``.
+        """
+        sorted_indices = self.get_sorted_slice_indices()
+        image_defs = self.image_defs
+
+        # define which keys which might vary across image volumes
+        dynamic_keys = ['cardiac phase number',
+                        'echo number',
+                        'label type',
+                        'image_type_mr',
+                        'dynamic scan number',
+                        'scanning sequence',
+                        'gradient orientation number',
+                        'diffusion b value number']
+
+        # remove dynamic keys that may not be present in older .PAR versions
+        dynamic_keys = [d for d in dynamic_keys if d in
+                        image_defs.dtype.fields]
+
+        non_unique_keys = []
+        for key in dynamic_keys:
+            ndim = image_defs[key].ndim
+            if ndim == 1:
+                num_unique = len(np.unique(image_defs[key]))
+            else:
+                raise ValueError("unexpected image_defs shape > 1D")
+            if num_unique > 1:
+                non_unique_keys.append(key)
+
+        # each key in dynamic keys will be identical across slices, so use
+        # the value at slice 1.
+        sl1_indices = image_defs['slice number'][sorted_indices] == 1
+
+        sort_info = OrderedDict()
+        for key in non_unique_keys:
+            sort_info[key] = image_defs[key][sorted_indices][sl1_indices]
+        return sort_info
+
 
 class PARRECImage(SpatialImage):
     """PAR/REC image"""

nibabel/tests/test_parrec.py

@@ -191,6 +191,16 @@ def test_header_scaling():
     assert_false(np.all(fp_scaling == dv_scaling))
 
 
+def test_header_volume_labels():
+    hdr = PARRECHeader(HDR_INFO, HDR_DEFS)
+    # check volume labels
+    vol_labels = hdr.get_volume_labels()
+    assert_equal(list(vol_labels.keys()), ['dynamic scan number'])
+    assert_array_equal(vol_labels['dynamic scan number'], [1, 2, 3])
+    # check that output is ndarray rather than list
+    assert_true(isinstance(vol_labels['dynamic scan number'], np.ndarray))
+
+
 def test_orientation():
     hdr = PARRECHeader(HDR_INFO, HDR_DEFS)
     assert_array_equal(HDR_DEFS['slice orientation'], 1)
@@ -285,6 +295,11 @@ def test_sorting_dual_echo_T1():
     # second half (volume 2) should all correspond to echo 2
     assert_equal(np.all(sorted_echos[n_half:] == 2), True)
 
+    # check volume labels
+    vol_labels = t1_hdr.get_volume_labels()
+    assert_equal(list(vol_labels.keys()), ['echo number'])
+    assert_array_equal(vol_labels['echo number'], [1, 2])
+
 
 def test_sorting_multiple_echos_and_contrasts():
     # This .PAR file has 3 echos and 4 image types (real, imaginary, magnitude,
@@ -327,6 +342,13 @@ def test_sorting_multiple_echos_and_contrasts():
     assert_equal(np.all(sorted_types[ntotal//2:3*ntotal//4] == 2), True)
     assert_equal(np.all(sorted_types[3*ntotal//4:ntotal] == 3), True)
 
+    # check volume labels
+    vol_labels = t1_hdr.get_volume_labels()
+    assert_equal(list(vol_labels.keys()), ['echo number', 'image_type_mr'])
+    assert_array_equal(vol_labels['echo number'], [1, 2, 3]*4)
+    assert_array_equal(vol_labels['image_type_mr'],
+                       [0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3])
+
 
 def test_sorting_multiecho_ASL():
     # For this .PAR file has 3 keys corresponding to volumes:
@@ -368,6 +390,15 @@ def test_sorting_multiecho_ASL():
     # check that slices vary fastest
     assert_array_equal(sorted_slices[:nslices], np.arange(1, nslices+1))
 
+    # check volume labels
+    vol_labels = asl_hdr.get_volume_labels()
+    assert_equal(list(vol_labels.keys()),
+                 ['echo number', 'label type', 'dynamic scan number'])
+    assert_array_equal(vol_labels['dynamic scan number'], [1]*6 + [2]*6)
+    assert_array_equal(vol_labels['label type'], [1]*3 + [2]*3 + [1]*3 + [2]*3)
+    assert_array_equal(vol_labels['echo number'], [1, 2, 3]*4)
+
+
 def test_vol_number():
     # Test algorithm for calculating volume number
     assert_array_equal(vol_numbers([1, 3, 0]), [0, 0, 0])

nibabel/tests/test_scripts.py

@@ -9,7 +9,7 @@
 import os
 from os.path import (dirname, join as pjoin, abspath, splitext, basename,
                      exists)
-import re
+import csv
 from glob import glob
 
 import numpy as np
@@ -18,8 +18,7 @@
 from ..loadsave import load
 from ..orientations import flip_axis, aff2axcodes, inv_ornt_aff
 
-from nose.tools import (assert_true, assert_false, assert_not_equal,
-                        assert_equal)
+from nose.tools import assert_true, assert_false, assert_equal
 from nose import SkipTest
 
 from numpy.testing import assert_almost_equal
@@ -333,3 +332,17 @@ def test_parrec2nii_with_data():
         data_sorted = img.get_data().copy()
         assert_almost_equal(data[..., np.argsort(DTI_PAR_BVALS)], data_sorted)
         del img
+
+        # Writes .ordering.csv if requested
+        run_command(['parrec2nii', '--overwrite', '--volume-info', dti_par])
+        assert_true(exists('DTI.ordering.csv'))
+        with open('DTI.ordering.csv', 'r') as csvfile:
+            csvreader = csv.reader(csvfile, delimiter=',')
+            csv_keys = csvreader.__next__()  # header row
+            nlines = 0  # count number of non-header rows
+            for line in csvreader:
+                nlines += 1
+
+        assert_equal(sorted(csv_keys), ['diffusion b value number',
+                                        'gradient orientation number'])
+        assert_equal(nlines, 8)  # 8 volumes present in DTI.PAR