ENH: PAR/REC data ordering and b-vector normalization

bin/parrec2nii

@@ -8,6 +8,8 @@ import numpy as np
 import sys
 import os
 import gzip
+import json
+from copy import deepcopy
 import nibabel
 import nibabel.parrec as pr
 from nibabel.parrec import one_line
@@ -47,6 +49,12 @@ def get_opt_parser():
                help=one_line(
                    """Output bvals/bvecs files in addition to NIFTI
                    image.""")))
+    p.add_option(
+        Option("--norm-bvs", action="store_true", dest="norm_bvs",
+               default=False,
+               help=one_line(
+                   """Normalize bvectors to 1.0, rescaling b-values
+                   appropriately.""")))
     p.add_option(
         Option("-d", "--dwell-time", action="store_true", default=False,
                dest="dwell_time",
@@ -63,6 +71,17 @@ 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", "--dimension-info", action="store_true", dest="dim_info",
+               default=False,
+               help=one_line(
+                   """Export .PAR dimension labels corresponding to the fourth
+                   dimension of the data.  The dimension info will be stored in
+                   JSON format with a filename matching the input .PAR, but
+                   with the extension .ordering.json.  Only labels that vary
+                   along the 4th dimension are exported.  (e.g. for a single
+                   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(
@@ -118,6 +137,16 @@ def error(msg, exit_code):
     sys.exit(exit_code)
 
 
+def convert_arrays_to_lists(input_dict):
+    # convert dictionary of numpy arrays to a dictionary of lists
+    output_dict = deepcopy(input_dict)
+    for key, value in output_dict.items():
+        if not isinstance(value, np.ndarray):
+            raise ValueError("all dictionary entries must be numpy arrays")
+        output_dict[key] = value.tolist()
+    return output_dict
+
+
 def proc_file(infile, opts):
     # figure out the output filename, and see if it exists
     basefilename = splitext_addext(os.path.basename(infile))[0]
@@ -207,7 +236,7 @@ def proc_file(infile, opts):
                   offset=offset)
     # write out bvals/bvecs if requested
     if opts.bvs:
-        bvals, bvecs = pr_hdr.get_bvals_bvecs()
+        bvals, bvecs = pr_hdr.get_bvals_bvecs(normalize_bvecs=opts.norm_bvs)
         if (bvals, bvecs) == (None, None):
             verbose('No DTI volumes detected, bvals and bvecs not written')
         else:
@@ -223,6 +252,14 @@ def proc_file(infile, opts):
                         fid.write('%s ' % val)
                     fid.write('\n')
 
+    # export data labels varying along the 4th dimensions
+    if opts.dim_info:
+        labels = pr_img.header.get_dimension_labels(collapse_slices=True)
+        if len(labels) > 0:
+            labels = convert_arrays_to_lists(labels)
+            with open(basefilename + '.ordering.json', 'w') as fid:
+                json.dump(labels, fid, sort_keys=True, indent=4)
+
     # write out dwell time if requested
     if opts.dwell_time:
         try:

nibabel/parrec.py

@@ -216,6 +216,18 @@ class PARRECError(Exception):
     pass
 
 
+def _unique_rows(a):
+    """
+    find unique rows of a 2D array.  based on discussion at:
+    http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array
+    """
+    if a.ndim != 2:
+        raise ValueError("expected 2D input")
+    b = np.ascontiguousarray(a).view(
+        np.dtype((np.void, a.dtype.itemsize * a.shape[1])))
+    return np.unique(b).view(a.dtype).reshape(-1, a.shape[1])
+
+
 def _split_header(fobj):
     """ Split header into `version`, `gen_dict`, `image_lines` """
     version = None
@@ -411,7 +423,7 @@ def parse_PAR_header(fobj):
     -------
     general_info : dict
         Contains all "General Information" from the header file
-    image_info : ndarray
+    image_defs : ndarray
         Structured array with fields giving all "Image information" in the
         header
     """
@@ -489,7 +501,7 @@ def __init__(self, file_like, header, mmap=True, scaling='dv'):
             True gives the same behavior as ``mmap='c'``.  If `file_like`
             cannot be memory-mapped, ignore `mmap` value and read array from
             file.
-        scaling : {'fp', 'dv'}, optional, keyword only
+        scaling : {'fp', 'dv', None}, optional, keyword only
             Type of scaling to use - see header ``get_data_scaling`` method.
         """
         if mmap not in (True, False, 'c', 'r'):
@@ -499,8 +511,11 @@ def __init__(self, file_like, header, mmap=True, scaling='dv'):
         self._shape = header.get_data_shape()
         self._dtype = header.get_data_dtype()
         self._slice_indices = header.get_sorted_slice_indices()
-        self._mmap=mmap
-        self._slice_scaling = header.get_data_scaling(scaling)
+        self._mmap = mmap
+        if scaling is None:
+            self._slice_scaling = None
+        else:
+            self._slice_scaling = header.get_data_scaling(scaling)
         self._rec_shape = header.get_rec_shape()
 
     @property
@@ -606,21 +621,27 @@ def get_echo_train_length(self):
         """Echo train length of the recording"""
         return self.general_info['epi_factor']
 
-    def get_q_vectors(self):
+    def get_q_vectors(self, normalize_bvecs=False):
         """Get Q vectors from the data
 
+        Parameters
+        ----------
+        normalize_bvecs : bool, optional
+            whether to scale the b-values by the norm of the b_vectors and then
+            renormalize any non-zero b_vectors to 1.0.
+
         Returns
         -------
         q_vectors : None or array
             Array of q vectors (bvals * bvecs), or None if not a diffusion
             acquisition.
         """
-        bvals, bvecs = self.get_bvals_bvecs()
+        bvals, bvecs = self.get_bvals_bvecs(normalize_bvecs=normalize_bvecs)
         if bvals is None and bvecs is None:
             return None
         return bvecs * bvals[:, np.newaxis]
 
-    def get_bvals_bvecs(self):
+    def get_bvals_bvecs(self, normalize_bvecs=False):
         """Get bvals and bvecs from data
 
         Returns
@@ -631,6 +652,9 @@ def get_bvals_bvecs(self):
         b_vectors : None or array
             Array of b vectors, shape (n_directions, 3), or None if not a
             diffusion acquisition.
+        normalize_bvecs : bool, optional
+            whether to scale the b-values by the norm of the b_vectors and then
+            renormalize any non-zero b_vectors to 1.0.
         """
         if self.general_info['diffusion'] == 0:
             return None, None
@@ -646,6 +670,20 @@ def get_bvals_bvecs(self):
         # All 3 values of bvecs should be same within volume
         assert not np.any(np.diff(bvecs, axis=0))
         bvecs = bvecs[0]
+        if normalize_bvecs:
+            # scale bvals by the norms of the bvecs then normalize any nonzero
+            # bvecs
+            bvec_norms = np.sqrt(np.sum(bvecs*bvecs, axis=1))
+            non_unity_indices = np.where(np.abs(bvec_norms - 1.0) > 1e-2)[0]
+            if len(non_unity_indices) > 0:
+                warnings.warn('Not all bvecs were normalized to 1.0. '
+                              'bvals will be scaled by the bvec norms')
+                bvals[non_unity_indices] *= bvec_norms[non_unity_indices]
+                # normalize any non-zero b-vectors to 1
+                non_zeros = np.where(bvec_norms[non_unity_indices] > 1e-2)[0]
+                if len(non_zeros) > 0:
+                    bvecs[non_unity_indices[non_zeros]] /= \
+                        bvec_norms[non_unity_indices[non_zeros]][:, np.newaxis]
         # rotate bvecs to match stored image orientation
         permute_to_psl = ACQ_TO_PSL[self.get_slice_orientation()]
         bvecs = apply_affine(np.linalg.inv(permute_to_psl), bvecs)
@@ -878,7 +916,7 @@ def get_data_scaling(self, method="dv"):
             slope = 1.0 / scale_slope
             intercept = rescale_intercept / (rescale_slope * scale_slope)
         else:
-            raise ValueError("Unknown scling method '%s'." % method)
+            raise ValueError("Unknown scaling method '%s'." % method)
         reorder = self.get_sorted_slice_indices()
         slope = slope[reorder]
         intercept = intercept[reorder]
@@ -919,6 +957,88 @@ def get_sorted_slice_indices(self):
         n_used = np.prod(self.get_data_shape()[2:])
         return np.lexsort(keys)[:n_used]
 
+    def get_dimension_labels(self, collapse_slices=True):
+        """ 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.
+
+        Parameters
+        ----------
+        collapse_slices : bool, optional
+            If True, only return volume indices (corresponding to the first
+            slice).  If False, return indices corresponding to individual
+            slices as well (with shape matching self.shape[2:]).
+
+        Returns
+        -------
+        sort_info : dict
+            Each key corresponds to a dynamically varying sequence dimension.
+            The ordering of each value corresponds to that returned by
+            ``self.get_sorted_slice_indices``.
+        """
+
+        # define which keys to store sorting info for
+        dynamic_keys = ['slice number',
+                        'cardiac phase number',
+                        'echo number',
+                        'label type',
+                        'image_type_mr',
+                        'dynamic scan number',
+                        'slice orientation',
+                        'image_display_orientation',  # ????
+                        'image angulation',
+                        'scanning sequence',
+                        'gradient orientation number',
+                        'diffusion b value number']
+
+        sorted_indices = self.get_sorted_slice_indices()
+        image_defs = self.image_defs
+
+        if collapse_slices:
+            dynamic_keys.remove('slice 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]))
+            elif ndim == 2:
+                # for 2D cases, e.g. 'image angulation'
+                num_unique = len(_unique_rows(image_defs[key]))
+            else:
+                raise ValueError("unexpected image_defs shape > 2D")
+            if num_unique > 1:
+                non_unique_keys.append(key)
+
+        if collapse_slices:
+            if 'slice orientation' in non_unique_keys:
+                raise ValueError("for non-unique slice orientation, need "
+                                 "collapse_slice=False")
+            if 'image angulation' in non_unique_keys:
+                raise ValueError("for non-unique image angulation, need "
+                                 "collapse_slice=False")
+            if 'image_display_orientation' in non_unique_keys:  # ???
+                raise ValueError("for non-unique display orientation, need "
+                                 "collapse_slice=False")
+            sl1_indices = image_defs['slice number'][sorted_indices] == 1
+
+        sort_info = {}
+        for key in non_unique_keys:
+            if collapse_slices:
+                sort_info[key] = image_defs[key][sorted_indices][sl1_indices]
+            else:
+                value = image_defs[key][sorted_indices]
+                sort_info[key] = value.reshape(self.get_data_shape()[2:],
+                                               order='F')
+        return sort_info
+
 
 class PARRECImage(SpatialImage):
     """PAR/REC image"""
@@ -948,7 +1068,7 @@ def from_file_map(klass, file_map, mmap=True, permit_truncated=False,
         permit_truncated : {False, True}, optional, keyword-only
             If False, raise an error for an image where the header shows signs
             that fewer slices / volumes were recorded than were expected.
-        scaling : {'dv', 'fp'}, optional, keyword-only
+        scaling : {'dv', 'fp', None}, optional, keyword-only
             Scaling method to apply to data (see
             :meth:`PARRECHeader.get_data_scaling`).
         """
@@ -982,7 +1102,7 @@ def from_filename(klass, filename, mmap=True, permit_truncated=False,
         permit_truncated : {False, True}, optional, keyword-only
             If False, raise an error for an image where the header shows signs
             that fewer slices / volumes were recorded than were expected.
-        scaling : {'dv', 'fp'}, optional, keyword-only
+        scaling : {'dv', 'fp', None}, optional, keyword-only
             Scaling method to apply to data (see
             :meth:`PARRECHeader.get_data_scaling`).
         """

nibabel/tests/test_parrec.py

@@ -170,6 +170,24 @@ def test_header_scaling():
     assert_false(np.all(fp_scaling == dv_scaling))
 
 
+def test_header_dimension_labels():
+    hdr = PARRECHeader(HDR_INFO, HDR_DEFS)
+    # check volume labels
+    vol_labels = hdr.get_dimension_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))
+    # check case with individual slice labels
+    slice_vol_labels = hdr.get_dimension_labels(collapse_slices=False)
+    # verify that both expected keys are present
+    assert_true('slice number' in slice_vol_labels)
+    assert_true('dynamic scan number' in slice_vol_labels)
+    # verify shape of labels matches final dimensions of data
+    assert_equal(slice_vol_labels['slice number'].shape,
+                 hdr.get_data_shape()[2:])
+
+
 def test_orientation():
     hdr = PARRECHeader(HDR_INFO, HDR_DEFS)
     assert_array_equal(HDR_DEFS['slice orientation'], 1)
@@ -320,6 +338,16 @@ def test_diffusion_parameters():
     assert_almost_equal(bvecs, DTI_PAR_BVECS[:, [2, 0, 1]])
     # Check q vectors
     assert_almost_equal(dti_hdr.get_q_vectors(), bvals[:, None] * bvecs)
+    # test bval/bvec normalization
+    bvals_norm, bvecs_norm = dti_hdr.get_bvals_bvecs(normalize_bvecs=True)
+    # bvecs were already normalized, so expect their values to remain unchanged
+    assert_almost_equal(bvecs_norm, DTI_PAR_BVECS[:, [2, 0, 1]])
+    bnorm = np.sqrt(np.sum(bvecs_norm*bvecs_norm, axis=1))
+    # all bvals where bvecs = [0, 0, 0] got set to 0?
+    assert_equal(np.abs(bvals_norm[bnorm == 0]).sum(), 0)
+    # all other bvals should not have changed
+    assert_almost_equal(bvals_norm[bnorm > 0],
+                        np.asarray(DTI_PAR_BVALS)[bnorm > 0])
 
 
 def test_null_diffusion_params():
@@ -520,6 +548,9 @@ def get_rec_shape(self):
         n_slices = np.prod(self._shape[2:])
         return self._shape[:2] + (n_slices,)
 
+    def sorted_labels(self, sorted_indices, collapse_slices):
+        return np.arange(self._shape[-1])
+
 
 def test_parrec_proxy():
     # Test PAR / REC proxy class, including mmap flags