|
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) |