Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Determine number of shells #64

Closed
oesteban opened this issue Jan 17, 2020 · 2 comments · May be fixed by #150
Closed

Determine number of shells #64

oesteban opened this issue Jan 17, 2020 · 2 comments · May be fixed by #150
Milestone
0.6.0

Comments

@oesteban
Copy link
Member

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)

@edickie edickie mentioned this issue Feb 3, 2020
Open
@oesteban oesteban added this to the 0.3.0 milestone Mar 23, 2020
@oesteban oesteban mentioned this issue Oct 13, 2020
Closed
@oesteban oesteban modified the milestones: 0.3.0, 0.4.0 Oct 13, 2020
@oesteban
Copy link
Member Author

Just to note it down: what happens with cartesian samplings of the q-space? (e.g., DSI) - will this give a proper output @edickie?

/cc @josephmje for what this affects the tutorial

@oesteban oesteban modified the milestones: 0.4.0, 0.5.0 Dec 10, 2020
@josephmje josephmje mentioned this issue Feb 9, 2021
Open
5 tasks
@oesteban oesteban mentioned this issue May 10, 2023
Merged
@oesteban
Copy link
Member Author

Done in MRIQC - we should upstream somewhere to share with dMRIPrep. Closing this for the time being.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
0.6.0
Development

Successfully merging a pull request may close this issue.

ENH: Add number of volumes per *b*-value to QC report josephmje/dmriprep
1 participant
@oesteban