-
Notifications
You must be signed in to change notification settings - Fork 24
Closed
Milestone
Description
I know @edickie and @mattcieslak were working on this.
I would suggest adding this as a cached method (i.e., only called once and then cached) of the DiffusionGradientTable object
dmriprep/dmriprep/utils/vectors.py
Lines 12 to 178 in 0a609b4
| class DiffusionGradientTable: | |
| """Data structure for DWI gradients.""" | |
| __slots__ = ['_affine', '_gradients', '_b_scale', '_bvecs', '_bvals', '_normalized', | |
| '_b0_thres', '_bvec_norm_epsilon'] | |
| def __init__(self, dwi_file=None, bvecs=None, bvals=None, rasb_file=None, | |
| b_scale=True, b0_threshold=B0_THRESHOLD, bvec_norm_epsilon=BVEC_NORM_EPSILON): | |
| """ | |
| Create a new table of diffusion gradients. | |
| Parameters | |
| ---------- | |
| dwi_file : str or os.pathlike or nibabel.spatialimage | |
| File path to the diffusion-weighted image series to which the | |
| bvecs/bvals correspond. | |
| bvals : str or os.pathlike or numpy.ndarray | |
| File path of the b-values. | |
| bvecs : str or os.pathlike or numpy.ndarray | |
| File path of the b-vectors. | |
| rasb_file : str or os.pathlike | |
| File path to a RAS-B gradient table. If rasb_file is provided, | |
| then bvecs and bvals will be dismissed. | |
| b_scale : bool | |
| Whether b-values should be normalized. | |
| """ | |
| self._b_scale = b_scale | |
| self._b0_thres = b0_threshold | |
| self._bvec_norm_epsilon = bvec_norm_epsilon | |
| self._gradients = None | |
| self._bvals = None | |
| self._bvecs = None | |
| self._affine = None | |
| self._normalized = False | |
| if dwi_file is not None: | |
| self.affine = dwi_file | |
| if rasb_file is not None: | |
| self.gradients = rasb_file | |
| if self.affine is not None: | |
| self.generate_vecval() | |
| elif dwi_file is not None and bvecs is not None and bvals is not None: | |
| self.bvecs = bvecs | |
| self.bvals = bvals | |
| self.generate_rasb() | |
| @property | |
| def affine(self): | |
| """Get the affine for RAS+/image-coordinates conversions.""" | |
| return self._affine | |
| @property | |
| def gradients(self): | |
| """Get gradient table (rasb).""" | |
| return self._gradients | |
| @property | |
| def bvecs(self): | |
| """Get the N x 3 list of bvecs.""" | |
| return self._bvecs | |
| @property | |
| def bvals(self): | |
| """Get the N b-values.""" | |
| return self._bvals | |
| @property | |
| def normalized(self): | |
| """Return whether b-vecs have been normalized.""" | |
| return self._normalized | |
| @affine.setter | |
| def affine(self, value): | |
| if isinstance(value, (str, Path)): | |
| dwi_file = nb.load(str(value)) | |
| self._affine = dwi_file.affine.copy() | |
| return | |
| if hasattr(value, 'affine'): | |
| self._affine = value.affine | |
| self._affine = np.array(value) | |
| @gradients.setter | |
| def gradients(self, value): | |
| if isinstance(value, (str, Path)): | |
| value = np.loadtxt(value, skiprows=1) | |
| self._gradients = value | |
| @bvecs.setter | |
| def bvecs(self, value): | |
| if isinstance(value, (str, Path)): | |
| value = np.loadtxt(str(value)).T | |
| else: | |
| value = np.array(value, dtype='float32') | |
| # Correct any b0's in bvecs misstated as 10's. | |
| value[np.any(abs(value) >= 10, axis=1)] = np.zeros(3) | |
| if self.bvals is not None and value.shape[0] != self.bvals.shape[0]: | |
| raise ValueError('The number of b-vectors and b-values do not match') | |
| self._bvecs = value | |
| @bvals.setter | |
| def bvals(self, value): | |
| if isinstance(value, (str, Path)): | |
| value = np.loadtxt(str(value)).flatten() | |
| if self.bvecs is not None and value.shape[0] != self.bvecs.shape[0]: | |
| raise ValueError('The number of b-vectors and b-values do not match') | |
| self._bvals = np.array(value) | |
| @property | |
| def b0mask(self): | |
| """Get a mask of low-b frames.""" | |
| return np.squeeze(self.gradients[..., -1] < self._b0_thres) | |
| def normalize(self): | |
| """Normalize (l2-norm) b-vectors.""" | |
| if self._normalized: | |
| return | |
| self._bvecs, self._bvals = normalize_gradients( | |
| self.bvecs, self.bvals, | |
| b0_threshold=self._b0_thres, | |
| bvec_norm_epsilon=self._bvec_norm_epsilon, | |
| b_scale=self._b_scale) | |
| self._normalized = True | |
| def generate_rasb(self): | |
| """Compose RAS+B gradient table.""" | |
| if self.gradients is None: | |
| self.normalize() | |
| _ras = bvecs2ras(self.affine, self.bvecs) | |
| self.gradients = np.hstack((_ras, self.bvals[..., np.newaxis])) | |
| def generate_vecval(self): | |
| """Compose a bvec/bval pair in image coordinates.""" | |
| if self.bvecs is None or self.bvals is None: | |
| if self.affine is None: | |
| raise TypeError( | |
| "Cannot generate b-vectors & b-values in image coordinates. " | |
| "Please set the corresponding DWI image's affine matrix.") | |
| self._bvecs = bvecs2ras(np.linalg.inv(self.affine), self.gradients[..., :-1]) | |
| self._bvals = self.gradients[..., -1].flatten() | |
| @property | |
| def pole(self): | |
| """ | |
| Check whether the b-vectors cover a full or just half a shell. | |
| If pole is all-zeros then the b-vectors cover a full sphere. | |
| """ | |
| self.generate_rasb() | |
| return calculate_pole(self.gradients[..., :-1], bvec_norm_epsilon=self._bvec_norm_epsilon) | |
| def to_filename(self, filename, filetype='rasb'): | |
| """Write files (RASB, bvecs/bvals) to a given path.""" | |
| if filetype.lower() == 'rasb': | |
| self.generate_rasb() | |
| np.savetxt(str(filename), self.gradients, | |
| delimiter='\t', header='\t'.join('RASB'), | |
| fmt=['%.8f'] * 3 + ['%g']) | |
| elif filetype.lower() == 'fsl': | |
| self.generate_vecval() | |
| np.savetxt('%s.bvec' % filename, self.bvecs.T, fmt='%.6f') | |
| np.savetxt('%s.bval' % filename, self.bvals, fmt='%.6f') | |
| else: | |
| raise ValueError('Unknown filetype "%s"' % filetype) |
Metadata
Metadata
Assignees
Labels
No labels