From 256d0d93437a82ddcff28b4977c4c93b96a23590 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 17 Nov 2014 19:32:18 -0500 Subject: [PATCH 01/18] docstring/typo fix and allow scaling=None --- nibabel/parrec.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/nibabel/parrec.py b/nibabel/parrec.py index d64887ac25..ead81dea16 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -411,7 +411,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 +489,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 +499,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 @@ -878,7 +881,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] @@ -948,7 +951,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 +985,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`). """ From 85ae8dc75a41e148e023c5eb4420d2358a5c9c89 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 17 Nov 2014 19:35:24 -0500 Subject: [PATCH 02/18] add option to normalize b-vectors and b-values --- bin/parrec2nii | 8 +++++++- nibabel/parrec.py | 19 ++++++++++++++++++- 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/bin/parrec2nii b/bin/parrec2nii index bf5832792d..e10744a3cc 100755 --- a/bin/parrec2nii +++ b/bin/parrec2nii @@ -47,6 +47,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", @@ -207,7 +213,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: diff --git a/nibabel/parrec.py b/nibabel/parrec.py index ead81dea16..875d419baf 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -623,7 +623,7 @@ def get_q_vectors(self): 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 @@ -634,6 +634,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 @@ -649,6 +652,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.linalg.norm(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]][:, None] # 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) From f5ef5fd73dcdfad5ecafc5a55408c24f225d7abc Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 17 Nov 2014 19:39:16 -0500 Subject: [PATCH 03/18] add sorted_labels function to PARRECHeader and corresponding info to PARRECArrayProxy --- bin/parrec2nii | 13 ++++++++ nibabel/parrec.py | 84 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+) diff --git a/bin/parrec2nii b/bin/parrec2nii index e10744a3cc..8e1c2755e7 100755 --- a/bin/parrec2nii +++ b/bin/parrec2nii @@ -8,6 +8,7 @@ import numpy as np import sys import os import gzip +import json import nibabel import nibabel.parrec as pr from nibabel.parrec import one_line @@ -229,6 +230,18 @@ def proc_file(infile, opts): fid.write('%s ' % val) fid.write('\n') + opts.dim_info = True # TODO: remove hard-coded value + if opts.dim_info and pr_img.dataobj._dim_4_labels is not None: + if len(list(pr_img.dataobj._dim_4_labels.keys())) > 0: + with open(basefilename + '.ordering.json', 'w') as fid: + json.dump(pr_img.dataobj._dim_4_labels, fid, sort_keys=True, + indent=4) + elif opts.dim_info and pr_img.dataobj._dim_3_4_labels is not None: + if len(list(pr_img.dataobj._dim_3_4_labels.keys())) > 0: + with open(basefilename + '.ordering.json', 'w') as fid: + json.dump(pr_img.dataobj._dim_3_4_labels, fid, sort_keys=True, + indent=4) + # write out dwell time if requested if opts.dwell_time: try: diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 875d419baf..1e6418100b 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -506,6 +506,21 @@ def __init__(self, file_like, header, mmap=True, scaling='dv'): self._slice_scaling = header.get_data_scaling(scaling) self._rec_shape = header.get_rec_shape() + try: + # store additional slice & volume label information + self._dim_4_labels = header.sorted_labels( + sorted_indices=self._slice_indices, + collapse_slices=True) + self._dim_3_4_labels = None + except: + # if non-unique slice orientations or image angulations, need to + # also keep labels for the 3rd (slice) dimension + self._dim_4_labels = None + self._dim_3_4_labels = header.sorted_labels( + sorted_indices=self._slice_indices, + collapse_slices=False) + + @property def shape(self): return self._shape @@ -939,6 +954,75 @@ def get_sorted_slice_indices(self): n_used = np.prod(self.get_data_shape()[2:]) return np.lexsort(keys)[:n_used] + def sorted_labels(self, sorted_indices=None, collapse_slices=True): + """ return a dictionary of lists of the varying values in + ``self.image_defs``. This is useful for custom data sorting. only + keys that have more than 1 unique value across the dataset will exist + in the returned `sort_info` dictionary. + + Parameters + ---------- + sorted_indices : array, optional + sorted slice indices as returned by + ``self.get_sorted_slice_indices``. + collapse_slices : bool, optional + if True, only return indices corresponding to the first slice + + Returns + ------- + sort_info : dict + Each key corresponds to a dynamically varying sequence dimension. + The ordering of values corresponds to that in sorted_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'] + + if sorted_indices is None: + sorted_indices = self.get_sorted_slice_indices() + + if collapse_slices: + dynamic_keys.remove('slice number') + + non_unique_keys = [] + for key in dynamic_keys: + if len(np.unique(self.image_defs[key])) > 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 image_display_orientation' in non_unique_keys: # ??? + raise ValueError("for non-unique display orientation, need " + "collapse_slice=False") + sl1_indices = self.image_defs['slice number'][sorted_indices] == 1 + + sort_info = {} + for key in non_unique_keys: + if collapse_slices: + sort_info[key] = list( + self.image_defs[key][sorted_indices][sl1_indices]) + else: + sort_info[key] = list( + self.image_defs[key][sorted_indices]) + + return sort_info + class PARRECImage(SpatialImage): """PAR/REC image""" From b223a89bb66594222b63f8ffb780ec6cd9ce3af5 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 18 Nov 2014 10:35:24 -0500 Subject: [PATCH 04/18] add sorted_labels property to FakeHeader in test_parrec.py --- nibabel/tests/test_parrec.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index 3a11660766..ccd323e815 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -520,6 +520,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 From 994fc5889ce85e7118597a4f0ff57ecdeb8a7bc5 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 18 Nov 2014 12:01:16 -0500 Subject: [PATCH 05/18] bugfix to replace np.unique() with _unique_rows() for 2D inputs needed for sorting vector properties. proper 2D ndarray to list of lists for JSON export in parrec2nii --- bin/parrec2nii | 35 +++++++++++++++++++++++++++-------- nibabel/parrec.py | 31 +++++++++++++++++++++++++------ 2 files changed, 52 insertions(+), 14 deletions(-) diff --git a/bin/parrec2nii b/bin/parrec2nii index 8e1c2755e7..4391340fa8 100755 --- a/bin/parrec2nii +++ b/bin/parrec2nii @@ -9,6 +9,7 @@ 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 @@ -125,6 +126,24 @@ def error(msg, exit_code): sys.exit(exit_code) +def sort_info_to_lists(sort_info): + sort_info = deepcopy(sort_info) + # convert from elements from numpy array to lists for easy JSON export + for key in sort_info: + ndim = sort_info[key].ndim + if ndim == 1: + sort_info[key] = list(sort_info[key]) + elif ndim == 2: + k = sort_info[key] + k_list = [] + for row in range(k.shape[0]): + k_list.append(list(k[row])) + sort_info[key] = k_list + else: + raise ValueError("only 1D and 2D arrays supported") + return sort_info + + def proc_file(infile, opts): # figure out the output filename, and see if it exists basefilename = splitext_addext(os.path.basename(infile))[0] @@ -232,15 +251,15 @@ def proc_file(infile, opts): opts.dim_info = True # TODO: remove hard-coded value if opts.dim_info and pr_img.dataobj._dim_4_labels is not None: - if len(list(pr_img.dataobj._dim_4_labels.keys())) > 0: - with open(basefilename + '.ordering.json', 'w') as fid: - json.dump(pr_img.dataobj._dim_4_labels, fid, sort_keys=True, - indent=4) + labels = pr_img.dataobj._dim_4_labels elif opts.dim_info and pr_img.dataobj._dim_3_4_labels is not None: - if len(list(pr_img.dataobj._dim_3_4_labels.keys())) > 0: - with open(basefilename + '.ordering.json', 'w') as fid: - json.dump(pr_img.dataobj._dim_3_4_labels, fid, sort_keys=True, - indent=4) + labels = pr_img.dataobj._dim_3_4_labels + else: + labels = None + if labels is not None and len(labels) > 0: + sort_info = sort_info_to_lists(labels) + with open(basefilename + '.ordering.json', 'w') as fid: + json.dump(sort_info, fid, sort_keys=True, indent=4) # write out dwell time if requested if opts.dwell_time: diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 1e6418100b..e40f26c318 100644 --- a/nibabel/parrec.py +++ b/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 @@ -996,8 +1008,17 @@ def sorted_labels(self, sorted_indices=None, collapse_slices=True): dynamic_keys.remove('slice number') non_unique_keys = [] + for key in dynamic_keys: - if len(np.unique(self.image_defs[key])) > 1: + ndim = self.image_defs[key].ndim + if ndim == 1: + num_unique = len(np.unique(self.image_defs[key])) + elif ndim == 2: + # for 2D cases, e.g. 'image angulation' + num_unique = len(_unique_rows(self.image_defs[key])) + else: + raise ValueError("unexpected image_defs shape > 2D") + if num_unique > 1: non_unique_keys.append(key) if collapse_slices: @@ -1014,13 +1035,11 @@ def sorted_labels(self, sorted_indices=None, collapse_slices=True): sort_info = {} for key in non_unique_keys: + if collapse_slices: - sort_info[key] = list( - self.image_defs[key][sorted_indices][sl1_indices]) + sort_info[key] = self.image_defs[key][sorted_indices][sl1_indices] else: - sort_info[key] = list( - self.image_defs[key][sorted_indices]) - + sort_info[key] = self.image_defs[key][sorted_indices] return sort_info From 23acdedae86e7215c413c7e3bdea854881f78c4b Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 10:31:55 -0500 Subject: [PATCH 06/18] clarify/streamline code based on feedback --- bin/parrec2nii | 56 +++++++++++++++++++------------------- nibabel/parrec.py | 68 +++++++++++++++++++++-------------------------- 2 files changed, 59 insertions(+), 65 deletions(-) diff --git a/bin/parrec2nii b/bin/parrec2nii index 4391340fa8..14497ba139 100755 --- a/bin/parrec2nii +++ b/bin/parrec2nii @@ -71,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=True, + 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( @@ -126,22 +137,14 @@ def error(msg, exit_code): sys.exit(exit_code) -def sort_info_to_lists(sort_info): - sort_info = deepcopy(sort_info) - # convert from elements from numpy array to lists for easy JSON export - for key in sort_info: - ndim = sort_info[key].ndim - if ndim == 1: - sort_info[key] = list(sort_info[key]) - elif ndim == 2: - k = sort_info[key] - k_list = [] - for row in range(k.shape[0]): - k_list.append(list(k[row])) - sort_info[key] = k_list - else: - raise ValueError("only 1D and 2D arrays supported") - return sort_info +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): @@ -249,17 +252,16 @@ def proc_file(infile, opts): fid.write('%s ' % val) fid.write('\n') - opts.dim_info = True # TODO: remove hard-coded value - if opts.dim_info and pr_img.dataobj._dim_4_labels is not None: - labels = pr_img.dataobj._dim_4_labels - elif opts.dim_info and pr_img.dataobj._dim_3_4_labels is not None: - labels = pr_img.dataobj._dim_3_4_labels - else: - labels = None - if labels is not None and len(labels) > 0: - sort_info = sort_info_to_lists(labels) - with open(basefilename + '.ordering.json', 'w') as fid: - json.dump(sort_info, fid, sort_keys=True, indent=4) + # export data labels varying along the 4th dimensions + if opts.dim_info: + if pr_img.dataobj._dim_4_labels is not None: + labels = pr_img.header.get_dimension_labels(collapse_slices=True) + else: + labels = None + if labels is not None and 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: diff --git a/nibabel/parrec.py b/nibabel/parrec.py index e40f26c318..95f53c67ec 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -224,7 +224,7 @@ def _unique_rows(a): if a.ndim != 2: raise ValueError("expected 2D input") b = np.ascontiguousarray(a).view( - np.dtype((np.void, a.dtype.itemsize * a.shape[1]))) + np.dtype((np.void, a.dtype.itemsize * a.shape[1]))) return np.unique(b).view(a.dtype).reshape(-1, a.shape[1]) @@ -518,21 +518,6 @@ def __init__(self, file_like, header, mmap=True, scaling='dv'): self._slice_scaling = header.get_data_scaling(scaling) self._rec_shape = header.get_rec_shape() - try: - # store additional slice & volume label information - self._dim_4_labels = header.sorted_labels( - sorted_indices=self._slice_indices, - collapse_slices=True) - self._dim_3_4_labels = None - except: - # if non-unique slice orientations or image angulations, need to - # also keep labels for the 3rd (slice) dimension - self._dim_4_labels = None - self._dim_3_4_labels = header.sorted_labels( - sorted_indices=self._slice_indices, - collapse_slices=False) - - @property def shape(self): return self._shape @@ -692,7 +677,7 @@ def get_bvals_bvecs(self, normalize_bvecs=False): 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]][:, None] + 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) @@ -966,25 +951,27 @@ def get_sorted_slice_indices(self): n_used = np.prod(self.get_data_shape()[2:]) return np.lexsort(keys)[:n_used] - def sorted_labels(self, sorted_indices=None, collapse_slices=True): - """ return a dictionary of lists of the varying values in - ``self.image_defs``. This is useful for custom data sorting. only - keys that have more than 1 unique value across the dataset will exist - in the returned `sort_info` dictionary. + 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 ---------- - sorted_indices : array, optional - sorted slice indices as returned by - ``self.get_sorted_slice_indices``. collapse_slices : bool, optional - if True, only return indices corresponding to the first slice + 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 values corresponds to that in sorted_indices. + The ordering of each value corresponds to that returned by + ``self.get_sorted_slice_indices``. """ # define which keys to store sorting info for @@ -1001,21 +988,26 @@ def sorted_labels(self, sorted_indices=None, collapse_slices=True): 'gradient orientation number', 'diffusion b value number'] - if sorted_indices is None: - sorted_indices = self.get_sorted_slice_indices() + sorted_indices = self.get_sorted_slice_indices() + image_defs = self.image_defs if collapse_slices: dynamic_keys.remove('slice number') - non_unique_keys = [] + # remove any dynamic keys that may not be implemented in older .PAR + # versions + for key in dynamic_keys: + if key not in image_defs: + dynamic_keys.remove(key) + non_unique_keys = [] for key in dynamic_keys: - ndim = self.image_defs[key].ndim + ndim = image_defs[key].ndim if ndim == 1: - num_unique = len(np.unique(self.image_defs[key])) + num_unique = len(np.unique(image_defs[key])) elif ndim == 2: # for 2D cases, e.g. 'image angulation' - num_unique = len(_unique_rows(self.image_defs[key])) + num_unique = len(_unique_rows(image_defs[key])) else: raise ValueError("unexpected image_defs shape > 2D") if num_unique > 1: @@ -1028,18 +1020,18 @@ def sorted_labels(self, sorted_indices=None, collapse_slices=True): if 'image angulation' in non_unique_keys: raise ValueError("for non-unique image angulation, need " "collapse_slice=False") - if 'image image_display_orientation' in non_unique_keys: # ??? + if 'image_display_orientation' in non_unique_keys: # ??? raise ValueError("for non-unique display orientation, need " "collapse_slice=False") - sl1_indices = self.image_defs['slice number'][sorted_indices] == 1 + sl1_indices = image_defs['slice number'][sorted_indices] == 1 sort_info = {} for key in non_unique_keys: - if collapse_slices: - sort_info[key] = self.image_defs[key][sorted_indices][sl1_indices] + sort_info[key] = image_defs[key][sorted_indices][sl1_indices] else: - sort_info[key] = self.image_defs[key][sorted_indices] + value = image_defs[key][sorted_indices] + sort_info[key] = value.reshape(self.shape[2:]) return sort_info From 68743ed40ba219f5fb66a57ae7930a79fe8d2c4b Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 10:53:55 -0500 Subject: [PATCH 07/18] bugfix to newly refactored code --- nibabel/parrec.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 95f53c67ec..60d0a78873 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -994,10 +994,9 @@ def get_dimension_labels(self, collapse_slices=True): if collapse_slices: dynamic_keys.remove('slice number') - # remove any dynamic keys that may not be implemented in older .PAR - # versions + # remove dynamic keys that may not be present in older .PAR versions for key in dynamic_keys: - if key not in image_defs: + if key not in image_defs.dtype.fields: dynamic_keys.remove(key) non_unique_keys = [] @@ -1031,7 +1030,8 @@ def get_dimension_labels(self, collapse_slices=True): sort_info[key] = image_defs[key][sorted_indices][sl1_indices] else: value = image_defs[key][sorted_indices] - sort_info[key] = value.reshape(self.shape[2:]) + sort_info[key] = value.reshape(self.get_data_shape()[2:], + order='F') return sort_info From 28f19a3f9a920e3925fd894d42efdcab08b9de58 Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 11:02:11 -0500 Subject: [PATCH 08/18] bugfix to newly refactored code --- bin/parrec2nii | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/bin/parrec2nii b/bin/parrec2nii index 14497ba139..1f5fa374f7 100755 --- a/bin/parrec2nii +++ b/bin/parrec2nii @@ -73,7 +73,7 @@ def get_opt_parser(): because it is not encoded in the PAR/REC format."""))) p.add_option( Option("-i", "--dimension-info", action="store_true", dest="dim_info", - default=True, + default=False, help=one_line( """Export .PAR dimension labels corresponding to the fourth dimension of the data. The dimension info will be stored in @@ -254,11 +254,8 @@ def proc_file(infile, opts): # export data labels varying along the 4th dimensions if opts.dim_info: - if pr_img.dataobj._dim_4_labels is not None: - labels = pr_img.header.get_dimension_labels(collapse_slices=True) - else: - labels = None - if labels is not None and len(labels) > 0: + 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) From bb64b087f2914a11e00e8fb08391ac71a9ed84f0 Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 11:18:14 -0500 Subject: [PATCH 09/18] parrec: test_header_dimension_labels() added --- nibabel/tests/test_parrec.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index ccd323e815..80322ed4e1 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -170,6 +170,22 @@ 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_equal(vol_labels['dynamic scan number'], [1, 2, 3]) + # check case with individual slice labels + slice_vol_labels = hdr.get_dimension_labels(collapse_slices=False) + # verify the 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) From 816f55838a0080311c513b261b83841b323b47f7 Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 11:21:23 -0500 Subject: [PATCH 10/18] slight revision to test --- nibabel/tests/test_parrec.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index 80322ed4e1..101740dfdb 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -176,9 +176,11 @@ def test_header_dimension_labels(): vol_labels = hdr.get_dimension_labels() assert_equal(list(vol_labels.keys()), ['dynamic scan number']) assert_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 the expected keys are present + # 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 From 9933117ba68515778d957a408bade1c638c8a755 Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 11:41:03 -0500 Subject: [PATCH 11/18] expand test_diffusion_parameters to test normalization case --- nibabel/tests/test_parrec.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index 101740dfdb..7bed4c346d 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -338,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.linalg.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) + # the other b-vectors were already normalized and should not have changed + assert_almost_equal(bvals_norm[bnorm > 0], + np.asarray(DTI_PAR_BVALS)[bnorm > 0]) def test_null_diffusion_params(): From c8d946e4d6400ad342959202fbd2e303bbeca112 Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 11:42:58 -0500 Subject: [PATCH 12/18] fix comment test_diffusion_parameters --- nibabel/tests/test_parrec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index 7bed4c346d..83eededbad 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -345,7 +345,7 @@ def test_diffusion_parameters(): bnorm = np.linalg.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) - # the other b-vectors were already normalized and should not have changed + # all other bvals should not have changed assert_almost_equal(bvals_norm[bnorm > 0], np.asarray(DTI_PAR_BVALS)[bnorm > 0]) From 7b04b221678da8553183fe487e4043d5ab50712c Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 11:51:45 -0500 Subject: [PATCH 13/18] added normalize_bvecs option to get_q_vectors --- nibabel/parrec.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 60d0a78873..e3c4e487fa 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -621,16 +621,22 @@ 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] From ef45763ecd4eb42bedc31b5c42338a2f0663b5b3 Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 19 Nov 2014 14:25:52 -0500 Subject: [PATCH 14/18] test bugfix: assert_equal to asssert_array_equal in new test --- nibabel/tests/test_parrec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index 83eededbad..cbb9432d1a 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -175,7 +175,7 @@ def test_header_dimension_labels(): # check volume labels vol_labels = hdr.get_dimension_labels() assert_equal(list(vol_labels.keys()), ['dynamic scan number']) - assert_equal(vol_labels['dynamic scan number'], [1, 2, 3]) + 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 From 4cb0f0b225d08441274987a021da58de8cd07d71 Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 20 Nov 2014 08:34:27 -0500 Subject: [PATCH 15/18] trim dynamic_keys via list comprehension --- nibabel/parrec.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/nibabel/parrec.py b/nibabel/parrec.py index e3c4e487fa..cf56f62987 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -1001,9 +1001,8 @@ def get_dimension_labels(self, collapse_slices=True): dynamic_keys.remove('slice number') # remove dynamic keys that may not be present in older .PAR versions - for key in dynamic_keys: - if key not in image_defs.dtype.fields: - dynamic_keys.remove(key) + dynamic_keys = [d for d in dynamic_keys if d in + image_defs.dtype.fields] non_unique_keys = [] for key in dynamic_keys: From 0cbec64fbde35bd3b56cde2a0394169a21a97b5e Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 20 Nov 2014 09:59:09 -0500 Subject: [PATCH 16/18] fix numpy1.5 compatibility of tests --- nibabel/tests/test_parrec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index cbb9432d1a..3ff4c9d273 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -342,7 +342,7 @@ def test_diffusion_parameters(): 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.linalg.norm(bvecs_norm, axis=1) + bnorm = np.sqrt(np.sum(bvecs_norm**2, 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 From 466ec30970115cda25303ad56121c8bd4841e703 Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 20 Nov 2014 10:09:30 -0500 Subject: [PATCH 17/18] x**2 -> x*x --- nibabel/tests/test_parrec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/tests/test_parrec.py b/nibabel/tests/test_parrec.py index 3ff4c9d273..a99c115d13 100644 --- a/nibabel/tests/test_parrec.py +++ b/nibabel/tests/test_parrec.py @@ -342,7 +342,7 @@ def test_diffusion_parameters(): 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**2, axis=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 From a23d976cdabbd55eb718b6a5d688480fa6465165 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 21 Nov 2014 11:21:02 -0500 Subject: [PATCH 18/18] remove another instance of np.linalg.norm for numpy 1.5 --- nibabel/parrec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/parrec.py b/nibabel/parrec.py index cf56f62987..2b81133550 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -673,7 +673,7 @@ def get_bvals_bvecs(self, normalize_bvecs=False): if normalize_bvecs: # scale bvals by the norms of the bvecs then normalize any nonzero # bvecs - bvec_norms = np.linalg.norm(bvecs, axis=1) + 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. '