From a737431ab6d1249cb1b810b55bee55b7809f4cb7 Mon Sep 17 00:00:00 2001 From: dPys Date: Wed, 16 Oct 2019 14:35:28 -0500 Subject: [PATCH 01/12] ENH: Vector-checking utilities Update dmriprep/utils/vectors.py Co-Authored-By: Ariel Rokem Update dmriprep/utils/vectors.py Co-Authored-By: Ariel Rokem Update dmriprep/utils/vectors.py Co-Authored-By: Ariel Rokem add tests, corrections to check_vecs & rescale_bvals Update vectors.py Change dir variable name to dvol Update dmriprep/utils/vectors.py to use Pathlib universally for all file naming Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/tests/test_vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban parent 4e0233f452832d4d33ddad7badb90f5996375835 author Derek Pisner <16432683+dPys@users.noreply.github.com> 1572104610 -0500 committer dPys 1572240020 -0500 Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban # This is a combination of 5 commits.tree 38e99c741b7d5fea0f6fe4d21948a9482d37a21b parent 4e0233f452832d4d33ddad7badb90f5996375835 author Derek Pisner <16432683+dPys@users.noreply.github.com> 1572104610 -0500 committer dPys 1572240020 -0500 Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban update interface to utilize various class methods Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban parent 4e0233f452832d4d33ddad7badb90f5996375835 author Derek Pisner <16432683+dPys@users.noreply.github.com> 1572104610 -0500 committer dPys 1572240020 -0500 Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban # This is a combination of 5 commits.tree 38e99c741b7d5fea0f6fe4d21948a9482d37a21b parent 4e0233f452832d4d33ddad7badb90f5996375835 author Derek Pisner <16432683+dPys@users.noreply.github.com> 1572104610 -0500 committer dPys 1572240020 -0500 Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/utils/vectors.py Co-Authored-By: Oscar Esteban Update dmriprep/interfaces/vectors.py Co-Authored-By: Oscar Esteban update interface to utilize various class methods add reorientation to interface, add signal clustering to dwi/bval check add reorientation to interface, add signal clustering to dwi/bval check Remove image reorientation, add vector inversion routines that force to RAS when saving gtab to .tsv and image is not in RAS, remove image_gradient_consistency_check function in order to create a separate PR for it, update tests Remove image reorientation, add vector inversion routines that force to RAS when saving gtab to .tsv and image is not in RAS, remove image_gradient_consistency_check function in order to create a separate PR for it, update tests --- dmriprep/interfaces/vectors.py | 91 ++++++ dmriprep/utils/tests/__init__.py | 0 dmriprep/utils/tests/test_vectors.py | 192 +++++++++++++ dmriprep/utils/vectors.py | 402 +++++++++++++++++++++++++++ 4 files changed, 685 insertions(+) create mode 100644 dmriprep/interfaces/vectors.py create mode 100644 dmriprep/utils/tests/__init__.py create mode 100644 dmriprep/utils/tests/test_vectors.py create mode 100644 dmriprep/utils/vectors.py diff --git a/dmriprep/interfaces/vectors.py b/dmriprep/interfaces/vectors.py new file mode 100644 index 00000000..05227b4c --- /dev/null +++ b/dmriprep/interfaces/vectors.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +from nipype.interfaces.base import SimpleInterface, BaseInterfaceInputSpec, TraitedSpec, File, traits, Directory, isdefined + + +class CleanVecsInputSpec(BaseInterfaceInputSpec): + """Input interface wrapper for CleanVecs""" + basedir = Directory(exists=True, mandatory=False) + dwi_file = File(exists=True, mandatory=True) + fbval = File(exists=False, mandatory=False) + fbvec = File(exists=False, mandatory=False) + rasb_tsv_file = File(exists=False, mandatory=False) + rescale = traits.Bool(True, usedefault=True) + save_fsl_style = traits.Bool(True, usedefault=True) + B0_NORM_EPSILON = traits.Float(50) + bval_scaling = traits.Bool(True, usedefault=True) + + +class CleanVecsOutputSpec(TraitedSpec): + """Output interface wrapper for CleanVecs""" + from dipy.core.gradients import GradientTable + gtab = traits.Instance(GradientTable) + is_hemi = traits.Bool() + pole = traits.Array() + b0_ixs = traits.Array() + slm = traits.Str() + fbval_out_file = traits.Any() + fbvec_out_file = traits.Any() + rasb_tsv_out_file = File(exists=True, mandatory=True) + + +class CleanVecs(SimpleInterface): + """Interface wrapper for CleanVecs""" + input_spec = CleanVecsInputSpec + output_spec = CleanVecsOutputSpec + + def _run_interface(self, runtime): + from dmriprep.utils.vectors import VectorTools + if isdefined(self.inputs.fbval) and isdefined(self.inputs.fbvec): + fbval = self.inputs.fbval + fbvec = self.inputs.fbvec + else: + fbval = None + fbvec = None + + if isdefined(self.inputs.rasb_tsv_file): + rasb_tsv_file = self.inputs.rasb_tsv_file + else: + rasb_tsv_file = None + + if isdefined(self.inputs.basedir): + basedir = self.inputs.basedir + else: + basedir = runtime.cwd + + vt = VectorTools(basedir, self.inputs.dwi_file, fbval, fbvec, rasb_tsv_file, self.inputs.B0_NORM_EPSILON, + self.inputs.bval_scaling) + + # Read in vectors + if rasb_tsv_file is not None: + vt.read_rasb_tsv() + elif fbval is not None and fbvec is not None: + vt.read_bvals_bvecs() + else: + raise OSError('VectorTools methods require either fbval/fbvec input files or an input rasb .tsv file.') + + # Rescale vectors + if self.inputs.rescale is True: + vt.rescale_vecs() + + # Build gradient table + gtab, b0_ixs = vt.build_gradient_table() + + # Check vector integrity + is_hemi, pole, slm = vt.checkvecs() + + if self.inputs.save_fsl_style is True: + fbval_out_file, fbvec_out_file = vt.save_vecs_fsl() + else: + fbval_out_file = None + fbvec_out_file = None + + self._results['gtab'] = gtab + self._results['is_hemi'] = is_hemi + self._results['pole'] = pole + self._results['slm'] = slm + self._results['b0_ixs'] = b0_ixs + self._results['fbval_out_file'] = fbval_out_file + self._results['fbvec_out_file'] = fbvec_out_file + self._results['rasb_tsv_out_file'] = vt.write_rasb_tsv() + return runtime diff --git a/dmriprep/utils/tests/__init__.py b/dmriprep/utils/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py new file mode 100644 index 00000000..97968697 --- /dev/null +++ b/dmriprep/utils/tests/test_vectors.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +import os +import pytest +import numpy as np +from pathlib import Path +from dmriprep.utils import vectors + + +def test_make_and_save_gradient_table(): + """Test for gradient table generation and saving from bval/bvec.""" + from dipy.data import fetch_sherbrooke_3shell + from dipy.io import read_bvals_bvecs + from os.path import expanduser, join + from dipy.core.gradients import GradientTable + + fetch_sherbrooke_3shell() + home = expanduser('~') + basedir = join(home, '.dipy', 'sherbrooke_3shell') + fdwi = join(basedir, 'HARDI193.nii.gz') + fbval = join(basedir, 'HARDI193.bval') + fbvec = join(basedir, 'HARDI193.bvec') + + # load bvecs/bvals. + bvals, bvecs = read_bvals_bvecs(fbval, fbvec) + + b0_thr = 50 + # Create vector directory if it does not already exist. + vector_dir = str(Path(basedir)/'vectors') + os.makedirs(vector_dir, exist_ok=True) + gtab_tsv_file = str(Path(vector_dir) / 'gradients.tsv') + + [gtab, b0_ixs] = vectors.make_gradient_table(bvals, bvecs, b0_thr) + vectors.write_gradients_tsv(bvecs, bvals, gtab_tsv_file) + + [bvals_read, bvecs_read] = vectors.read_gradients_tsv(gtab_tsv_file) + + assert gtab is not None and gtab.__class__ == GradientTable + assert gtab.b0_threshold == 50 + assert np.sum(gtab.b0s_mask) == len(b0_ixs) + assert os.path.isfile(gtab_tsv_file) is True + assert np.isclose(bvecs_read.all(), bvecs.all()) + assert np.isclose(bvals_read.all(), bvals.all()) + os.remove(gtab_tsv_file) + + +def test_rescale_vectors(): + """Test vector rescaling.""" + from dipy.data import fetch_sherbrooke_3shell + from os.path import expanduser, join + from dipy.io import read_bvals_bvecs + + fetch_sherbrooke_3shell() + home = expanduser('~') + basedir = join(home, '.dipy', 'sherbrooke_3shell') + fbval = join(basedir, 'HARDI193.bval') + fbvec = join(basedir, 'HARDI193.bvec') + + bval_scaling = True + b0_thr = 50 + + # load bvecs/bvals. + bvals, bvecs = read_bvals_bvecs(fbval, fbvec) + + # Get b order of magnitude + bmag = int(np.log10(np.max(bvals))) - 1 + + [bval_normed, bvec_normed] = vectors.rescale_vectors(bvals, bvecs, bval_scaling, bmag, b0_thr) + assert np.max(bvec_normed) < 1.0 and np.min(bvec_normed) > -1.0 + assert len(bvecs) == len(bvec_normed) + assert np.sum(np.count_nonzero(bvecs, axis=1).astype('bool')) == np.sum(np.count_nonzero(bvec_normed, + axis=1).astype('bool')) + assert len(bvals) == len(bval_normed) + assert np.sum(np.count_nonzero(bvals, axis=0).astype('bool')) == np.sum(np.count_nonzero(bval_normed, + axis=0).astype('bool')) + + +def test_vector_checkers(): + """Test for vector corruption.""" + from dipy.data import fetch_sherbrooke_3shell + from os.path import expanduser, join + from dipy.io import read_bvals_bvecs + import itertools + + fetch_sherbrooke_3shell() + home = expanduser('~') + basedir = join(home, '.dipy', 'sherbrooke_3shell') + fbval = join(basedir, 'HARDI193.bval') + fbvec = join(basedir, 'HARDI193.bvec') + + b0_thr = 50 + + # load bvecs/bvals. + bvals, bvecs = read_bvals_bvecs(fbval, fbvec) + + # Perform various corruption checks using synthetic corrupted bval-bvec. + bval_short = bvals[:-1] + bval_long = np.hstack([bvals, bvals[-1]]) + bval_no_b0 = bvals.copy() + bval_no_b0[0] = 51 + bval_odd_b0 = bvals.copy() + bval_odd_b0[np.where(abs(bval_odd_b0) == 0)] = 0.00000001 + bvec_short = bvecs[:-1] + bvec_long = np.vstack([bvecs, 0.5*bvecs[-1]]) + bvec_no_b0 = bvecs.copy() + bvec_no_b0[np.where(np.all(abs(bvec_no_b0) == 0, axis=1) == True)] = [1, 1, 1] + bvec_odd_b0 = bvecs.copy() + bvec_odd_b0[np.where(np.all(abs(bvec_odd_b0) == 0, axis=1) == True)] = [10, 10, 10] + + bval_lists = [bvals, bval_short, bval_long, bval_no_b0, bval_odd_b0] + bvec_lists = [bvecs, bvec_short, bvec_long, bvec_no_b0, bvec_odd_b0] + for bval, bvec in list(itertools.product(bval_lists, bvec_lists)): + [bvecs_checked, bvals_checked] = vectors.check_corruption(bvec, bval, b0_thr) + if (bval is bvals) and (bvec is bvecs): + assert bvecs_checked is not None + assert bvals_checked is not None + else: + with pytest.raises(Exception) as e_info: + assert bvecs_checked is None + assert bvals_checked is None + + +def test_hemisphericity(): + """Test vector hemisphere coverage and second-level model setting based on that coverage.""" + from dipy.data import fetch_sherbrooke_3shell + from dipy.io import read_bvals_bvecs + from os.path import expanduser, join + + fetch_sherbrooke_3shell() + home = expanduser('~') + basedir = join(home, '.dipy', 'sherbrooke_3shell') + fbval = join(basedir, 'HARDI193.bval') + fbvec = join(basedir, 'HARDI193.bvec') + + b0_thr = 50 + + # load bvecs/bvals. + bvals, bvecs = read_bvals_bvecs(fbval, fbvec) + [gtab, _] = vectors.make_gradient_table(bvals, bvecs, b0_thr) + + # Check hemisphere coverage. + [is_hemi, pole] = vectors.is_hemispherical(np.round(gtab.bvecs, 8)[~(np.round(gtab.bvecs, 8) == 0).all(1)]) + + # Use sampling scheme to set second-level model for eddy correction (this step would conceivably apply to any + # non-FSL form of eddy correction as well. + slm = vectors.set_slm(gtab, is_hemi) + assert is_hemi is False + assert isinstance(pole, (np.ndarray, np.generic)) is True + assert slm == 'none' + + +def test_clean_vecs(): + """Test CleanVecs interface functionality.""" + from dipy.data import fetch_sherbrooke_3shell + from os.path import expanduser, join + from dipy.core.gradients import GradientTable + from dmriprep.interfaces.vectors import CleanVecs + import itertools + + fetch_sherbrooke_3shell() + home = expanduser('~') + basedir = join(home, '.dipy', 'sherbrooke_3shell') + fdwi = join(basedir, 'HARDI193.nii.gz') + fbval = join(basedir, 'HARDI193.bval') + fbvec = join(basedir, 'HARDI193.bvec') + + # Smoke test parameter combos + for bool_list in list(itertools.product([True, False], repeat=3)): + cv = CleanVecs() + cv.inputs.B0_NORM_EPSILON = 50 + cv.inputs.bval_scaling = bool_list[0] + cv.inputs.rescale = bool_list[1] + cv.inputs.save_fsl_style = bool_list[2] + cv.inputs.dwi_file = fdwi + cv.inputs.fbval = fbval + cv.inputs.fbvec = fbvec + cv.inputs.basedir = basedir + + res = cv.run() + + assert res.outputs.gtab is not None and res.outputs.gtab.__class__ == GradientTable + assert res.outputs.is_hemi is False + assert isinstance(res.outputs.pole, (np.ndarray, np.generic)) is True + assert np.sum(res.outputs.gtab.b0s_mask) == len(res.outputs.b0_ixs) + assert os.path.isfile(res.outputs.rasb_tsv_out_file) is True + if cv.inputs.save_fsl_style is True: + assert os.path.isfile(res.outputs.fbval_out_file) is True + assert os.path.isfile(res.outputs.fbvec_out_file) is True + os.remove(res.outputs.fbval_out_file) + os.remove(res.outputs.fbvec_out_file) + assert res.outputs.slm == 'none' + os.remove(res.outputs.rasb_tsv_out_file) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py new file mode 100644 index 00000000..30fd4cc8 --- /dev/null +++ b/dmriprep/utils/vectors.py @@ -0,0 +1,402 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +import os +import nibabel as nib +import numpy as np + + +def make_gradient_table(bvals, bvecs, B0_NORM_EPSILON): + """ + Create gradient table from bval/bvec. + + Parameters + ---------- + bvals : 1d array + Raw b-values float array. + bvecs : m x n 2d array + Raw b-vectors array. + B0_NORM_EPSILON : float + Gradient threshold below which volumes and vectors are considered B0's. + + Returns + ------- + gtab : Obj + DiPy object storing diffusion gradient information. + b0_ixs : 1d array + Array containing the indices of the b0s. + """ + from dipy.core.gradients import gradient_table + from dmriprep.utils.vectors import check_corruption + + # Perform initial corruption check, correcting for obvious format/character encoding anomalies. + bvecs, bvals = check_corruption(bvecs, bvals, B0_NORM_EPSILON) + + # Creating the gradient table + gtab = gradient_table(bvals, bvecs, b0_threshold=B0_NORM_EPSILON) + + b0_ixs = np.where(gtab.b0s_mask == True)[0] + + return gtab, b0_ixs + + +def read_gradients_tsv(rasb_tsv_file): + """ + Read bvals/bvecs from a RAS-B tsv file. + + Parameters + ---------- + rasb_tsv_file : pathlike or str + Path to .tsv file containing RAS-B gradient table. + + Returns + ------- + bvals : 1d array + Raw b-values float array. + bvecs : m x n 2d array + Raw b-vectors array. + """ + RASB_table = np.genfromtxt(rasb_tsv_file, delimiter="\t") + bvecs = RASB_table[1:, 0:3] + bvals = RASB_table[1:, 3] + return bvals, bvecs + + +def write_gradients_tsv(bvecs, bvals, out_file): + """ + Save bvals/bveccs as RAS-B tsv. + + Parameters + ---------- + bvals : 1d array + Raw b-values float array. + bvecs : m x n 2d array + Raw b-vectors array. + out_file : pathlike or str + Path to .tsv file containing RAS-B gradient table. + """ + RASB_table = np.column_stack([bvecs, bvals]) + with open(out_file, 'wb') as f: + f.write(b'R\tA\tS\tB\n') + np.savetxt(f, RASB_table.astype('float'), fmt='%.8f', delimiter="\t") + f.close() + return + + +def save_corrected_bval_bvec(vector_dir, bvec_corrected, bval_corrected, bval_scaling): + """ + Normalize b-vectors to be of unit length for the non-zero b-values. + If the b-vector row reflects a B0, the vector is untouched. + + Parameters + ---------- + vector_dir : str + Path the output vector directory to save the corrected bval/bvec files. + bvec_corrected : m x n 2d array + Corrected b-vectors array. + bval_corrected : 1d array + Corrected b-values float array. + bval_scaling : bool + If True, then normalizes b-val by the square of the vector amplitude. + + Returns + ------- + fbval_corrected : pathlike or str + File path of the corrected b-values. + fbvec_corrected : pathlike or str + File path of the corrected b-vectors. + """ + from pathlib import Path + + fbvec_corrected = str(Path(vector_dir)/'bvec_resc_corr.bvec') + np.savetxt(fbvec_corrected, bvec_corrected.T.astype('float'), fmt='%.8f') + + if bval_scaling is True: + fbval_corrected = str(Path(vector_dir) / 'bval_resc_corr.bval') + else: + fbval_corrected = str(Path(vector_dir) / 'bval_corr.bval') + np.savetxt(fbval_corrected, bval_corrected.T.astype('int'), fmt='%d', newline=' ') + return fbval_corrected, fbvec_corrected + + +def rescale_bvec(bvec, B0_NORM_EPSILON): + """ + Normalize b-vectors to be of unit length for the non-zero b-values. + If the b-vector row reflects a B0, the vector is untouched. + + Parameters + ---------- + bvec : m x n 2d array + Raw b-vectors array. + B0_NORM_EPSILON : float + Gradient threshold below which volumes and vectors are considered B0's. + + Returns + ------- + bvec : m x n 2d array + Unit-normed b-vectors array. + """ + axis = int(bvec.shape[0] == 3) + b0s = np.linalg.norm(bvec) < B0_NORM_EPSILON + + # Rescale b-vecs, skipping b0's, on the appropriate axis to unit-norm length. + bvec[~b0s] = bvec[~b0s] / np.linalg.norm(bvec[~b0s], axis=axis) + return bvec + + +def rescale_bval(bval, bvec, B0_NORM_EPSILON): + """ + Normalize b-val by the square of the vector amplitude. + If the b-vector row reflects a B0, the corresponding b-val entry is untouched. + + Parameters + ---------- + bval : 1d array + Raw b-values float array. + bvec : m x n 2d array + Raw b-vectors array. + B0_NORM_EPSILON : float + Gradient threshold below which volumes and vectors are considered B0's. + + Returns + ------- + bval : 1d int array + Vector amplitude square normed b-values array. + """ + axis = int(bvec.shape[0] == 3) + b0s_loc = np.where(bval < B0_NORM_EPSILON) + + # Iterate through b-vec rows to perform amplitude correction to each corresponding b-val entry. + bval = bval.copy() + for dvol in range(len(bvec)): + if dvol not in b0s_loc: + bval[dvol] = bval[dvol] * np.linalg.norm(bvec[dvol], axis=axis)**2 + return np.round(bval, 0) + + +def is_hemispherical(vecs): + """ + Determine whether all points on a unit sphere lie in the same hemisphere. + + Parameters + ---------- + vecs : numpy.ndarray + 2D numpy array with shape (N, 3) where N is the number of points. + All points must lie on the unit sphere. + + Returns + ------- + is_hemi : bool + If True, one can find a hemisphere that contains all the points. + If False, then the points do not lie in any hemisphere + pole : numpy.ndarray + If `is_hemi == True`, then pole is the "central" pole of the + input vectors. Otherwise, pole is the zero vector. + + References + ---------- + https://rstudio-pubs-static.s3.amazonaws.com/27121_a22e51b47c544980bad594d5e0bb2d04.html + """ + import itertools + if vecs.shape[1] != 3: + raise ValueError("Input vectors must be 3D vectors") + if not np.allclose(1, np.linalg.norm(vecs, axis=1)): + raise ValueError("Input vectors must be unit vectors") + + # Generate all pairwise cross products. + v0, v1 = zip(*[p for p in itertools.permutations(vecs, 2)]) + cross_prods = np.cross(v0, v1) + + # Normalize them. + cross_prods /= np.linalg.norm(cross_prods, axis=1)[:, np.newaxis] + + # `cross_prods` now contains all candidate vertex points for "the polygon" in the reference. "The polygon" is a + # subset. Find which points belong to the polygon using a dot product test with each of the original vectors. + angles = np.arccos(np.dot(cross_prods, vecs.transpose())) + + # And test whether it is orthogonal or less. + dot_prod_test = angles <= np.pi / 2.0 + + # If there is at least one point that is orthogonal or less to each input vector, then the points lie on some + # hemisphere. + is_hemi = len(vecs) in np.sum(dot_prod_test.astype(int), axis=1) + + if is_hemi: + vertices = cross_prods[np.sum(dot_prod_test.astype(int), axis=1) == len(vecs)] + pole = np.mean(vertices, axis=0) + pole /= np.linalg.norm(pole) + else: + pole = np.array([0.0, 0.0, 0.0]) + return is_hemi, pole + + +def set_slm(gtab, is_hemi): + """ + Evaluate eddy compatibility and set second level model that specifies the mathematical form for how the + diffusion gradients cause eddy currents. If the data has few directions and/or is has not been sampled on the whole + sphere it can be advantageous to specify --slm=linear. + + Parameters + ---------- + gtab : Obj + DiPy object storing diffusion gradient information. + is_hemi : bool + If True, one can find a hemisphere that contains all the points. + If False, then the points do not lie in any hemisphere. + + Returns + ------- + slm : str + Second-level model to use for eddy correction. + """ + # Warn user if gradient table is not conducive to eddy correction. + if ((np.sum(np.invert(gtab.b0s_mask)) < 10) and (max(gtab.bvals) < 1500)) or ((np.sum(np.invert(gtab.b0s_mask)) < 30) and (max(gtab.bvals) < 5000)): + raise UserWarning('Too few directions in this data. Use of eddy is not recommended.') + + if is_hemi is True: + slm = 'linear' + print("B-vectors for this data are hemispherical. To ensure adequate eddy current correction, eddy will be run " + "using the --slm=linear option.") + elif np.sum(np.invert(gtab.b0s_mask)) < 30: + slm = 'linear' + print('Fewer than 30 diffusion-weighted directions detected. To ensure adequate eddy current correction, eddy ' + 'will be run using the --slm=linear option') + else: + slm = 'none' + return slm + + +def check_corruption(bvecs, bvals, B0_NORM_EPSILON): + """ + Performs a variety of basic check to ensure bval/bvec files are not corrupted. + + Parameters + ---------- + bvecs : m x n 2d array + Raw b-vectors array. + bvals : 1d array + Raw b-values float array. + B0_NORM_EPSILON : float + Gradient threshold below which volumes and vectors are considered B0's. + + Returns + ------- + bvecs : m x n 2d array + Raw b-vectors array. + bvals : 1d array + Raw b-values float array. + """ + # Correct any b0's in bvecs misstated as 10's. + bvecs[np.where(np.any(abs(bvecs) >= 10, axis=1) == True)] = [1, 0, 0] + + # Check for bval-bvec discrepancy. + if len(bvecs[np.where(np.logical_and(bvals > B0_NORM_EPSILON, np.where(np.all(abs(bvecs)) == 0)))]) > 0: + raise ValueError("Corrupted bval/bvecs detected. Check to ensure volumes with a diffusion weighting are not " + "listed as B0s in the bvecs.") + + # Ensure that B0's in the bvecs are represented are zero floats. + bvecs[np.where(bvals < B0_NORM_EPSILON)] = float(0) + + return bvecs, bvals + + +def rescale_vectors(bvals, bvecs, bval_scaling, bmag, B0_NORM_EPSILON): + """ + Performs a variety of rescaling routines on the b vectors and values. + + Parameters + ---------- + bvals : 1d array + Raw b-values float array. + bvecs : m x n 2d array + Raw b-vectors array. + bval_scaling : bool + If True, then normalizes b-val by the square of the vector amplitude. + bmag : int + The order of magnitude to round the b-values. + B0_NORM_EPSILON : float + Gradient threshold below which volumes and vectors are considered B0's. + + Returns + ------- + bval_mag_norm : 1d array + Rescaled b-values float array. + bvec_rescaled : m x n 2d array + Rescaled b-vectors array. + """ + from dipy.core.gradients import round_bvals + from dmriprep.utils.vectors import rescale_bvec, rescale_bval + # Rescale bvals by vector amplitudes. + if bval_scaling is True: + bvals = rescale_bval(bvals.astype('float'), bvecs.astype('float'), B0_NORM_EPSILON) + + # Round bvals by bmag. + bval_mag_norm = round_bvals(bvals, bmag=bmag) + + # Rescale bvecs. + bvec_rescaled = rescale_bvec(bvecs, B0_NORM_EPSILON) + return bval_mag_norm, bvec_rescaled + + +class VectorTools(object): + """ + Class of vector-handling methods for io and integrity checking. + + Parameters + ---------- + basedir : str + Path the output base directory for derivative dmriprep data. + dwi_file : str + Optionally provide a file path to the diffusion-weighted image series to which the + bvecs/bvals should correspond. + fbval : str + File path of the b-values. + fbvec : str + File path of the b-vectors. + rasb_tsv_file : pathlike or str + Path to .tsv file containing RAS-B gradient table. + B0_NORM_EPSILON : float + Gradient threshold below which volumes and vectors are considered B0's. Default is 50. + bval_scaling : bool + If True, then normalizes b-val by the square of the vector amplitude. Default is True. + """ + from pathlib import Path + from dipy.io import read_bvals_bvecs + from dmriprep.utils.vectors import make_gradient_table, is_hemispherical, check_corruption, save_corrected_bval_bvec, make_gradients_tsv, set_slm + + # Create vector directory if it does not already exist. + vector_dir = str(Path(sesdir)/'vectors') + os.makedirs(vector_dir, exist_ok=True) + + # load bvecs/bvals. + bvals, bvecs = read_bvals_bvecs(fbval, fbvec) + + # Get b order of magnitude + bmag = int(np.log10(np.max(bvals))) - 1 + + # Perform initial corruption check, correcting for obvious anomalies. + [bvecs, bvals] = check_corruption(bvecs, bvals, B0_NORM_EPSILON) + + # Perform b vector/value rescaling. + [bval_mag_norm, bvec_rescaled] = rescale_vectors(bvals, bvecs, bval_scaling, bmag, B0_NORM_EPSILON) + + # Save corrected bval/bvec to file in FSL style to ensure-backwards compatibility. + [fbval_rescaled, fbvec_rescaled] = save_corrected_bval_bvec(vector_dir, bvec_rescaled, bval_mag_norm, bval_scaling) + + # Make gradient table. + [gtab, b0_ixs] = make_gradient_table(fbval_rescaled, fbvec_rescaled, B0_NORM_EPSILON) + + # Save gradient table to tsv. + gtab_tsv_file = str(Path(vector_dir)/'gradients.tsv') + make_gradients_tsv(gtab, gtab_tsv_file) + + # Check consistency of dwi image and gradient vectors. + image_gradient_consistency_check(dwi_file, b0_ixs, gtab, len(np.unique(gtab.bvals)), bmag) + + # Check hemisphere coverage. + [is_hemi, pole] = is_hemispherical(np.round(gtab.bvecs, 8)[~(np.round(gtab.bvecs, 8) == 0).all(1)]) + + # Use sampling scheme to set second-level model for eddy correction (this step would conceivably apply to any + # non-FSL form of eddy correction as well. + slm = set_slm(gtab, is_hemi) + + return gtab, is_hemi, pole, b0_ixs, gtab_tsv_file, slm From 4cd6731518d4b0a5bd2942045481ad1967616229 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 1 Nov 2019 00:30:55 -0700 Subject: [PATCH 02/12] Revision of the branch I have revised the utilities and the nipype interface, trying to simplify the code (without missing functionality). Some assumptions that can be made on inputs: * RASB tables will be BIDS-valid, and hence normalized, with absolute b-vals, etc. * bvec+bval will be in FSL format (as mandated by BIDS), and hence in image coordinates. In general, I've tried to remove repetition of code sections, added doctests that will serve for documentation, minimized dependencies, checked code style. I haven't gone through the tests, but they will need a deep revision. I would recommend using pytest fixtures to minimize the lines of code and automate some clerical tasks (e.g., setting up data, changing to a temporal folder, etc.). --- dmriprep/interfaces/vectors.py | 184 +++++----- dmriprep/utils/vectors.py | 626 +++++++++++++++------------------ 2 files changed, 382 insertions(+), 428 deletions(-) diff --git a/dmriprep/interfaces/vectors.py b/dmriprep/interfaces/vectors.py index 05227b4c..74c594a5 100644 --- a/dmriprep/interfaces/vectors.py +++ b/dmriprep/interfaces/vectors.py @@ -1,91 +1,101 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -from nipype.interfaces.base import SimpleInterface, BaseInterfaceInputSpec, TraitedSpec, File, traits, Directory, isdefined - - -class CleanVecsInputSpec(BaseInterfaceInputSpec): - """Input interface wrapper for CleanVecs""" - basedir = Directory(exists=True, mandatory=False) +"""Handling the gradient table.""" +from pathlib import Path +import numpy as np +from nipype.utils.filemanip import fname_presuffix +from nipype.interfaces.base import ( + SimpleInterface, BaseInterfaceInputSpec, TraitedSpec, + File, traits, isdefined +) +from ..utils.vectors import DiffusionGradientTable, B0_THRESHOLD, BVEC_NORM_EPSILON + + +class _CheckGradientTableInputSpec(BaseInterfaceInputSpec): dwi_file = File(exists=True, mandatory=True) - fbval = File(exists=False, mandatory=False) - fbvec = File(exists=False, mandatory=False) - rasb_tsv_file = File(exists=False, mandatory=False) - rescale = traits.Bool(True, usedefault=True) - save_fsl_style = traits.Bool(True, usedefault=True) - B0_NORM_EPSILON = traits.Float(50) - bval_scaling = traits.Bool(True, usedefault=True) - - -class CleanVecsOutputSpec(TraitedSpec): - """Output interface wrapper for CleanVecs""" - from dipy.core.gradients import GradientTable - gtab = traits.Instance(GradientTable) - is_hemi = traits.Bool() - pole = traits.Array() - b0_ixs = traits.Array() - slm = traits.Str() - fbval_out_file = traits.Any() - fbvec_out_file = traits.Any() - rasb_tsv_out_file = File(exists=True, mandatory=True) - - -class CleanVecs(SimpleInterface): - """Interface wrapper for CleanVecs""" - input_spec = CleanVecsInputSpec - output_spec = CleanVecsOutputSpec + in_bvec = File(exists=True, xor=['in_rasb']) + in_bval = File(exists=True, xor=['in_rasb']) + in_rasb = File(exists=True, xor=['in_bval', 'in_bvec']) + b0_threshold = traits.Float(B0_THRESHOLD, usedefault=True) + bvec_norm_epsilon = traits.Float(BVEC_NORM_EPSILON, usedefault=True) + b_scale = traits.Bool(True, usedefault=True) + + +class _CheckGradientTableOutputSpec(TraitedSpec): + out_rasb = File(exists=True) + out_bval = File(exists=True) + out_bvec = File(exists=True) + full_sphere = traits.Bool() + pole = traits.Tuple(traits.Float, traits.Float, traits.Float) + b0_ixs = traits.List(traits.Int) + + +class CheckGradientTable(SimpleInterface): + """ + Ensure the correctness of the gradient table. + + Example + ------- + + >>> os.chdir(tmpdir) + >>> check = CheckGradientTable( + ... dwi_file=str(data_dir / 'dwi.nii.gz'), + ... in_rasb=str(data_dir / 'dwi.tsv')).run() + >>> check.outputs.pole + (0.0, 0.0, 0.0) + >>> check.outputs.full_sphere + True + + >>> check = CheckGradientTable( + ... dwi_file=str(data_dir / 'dwi.nii.gz'), + ... in_bvec=str(data_dir / 'bvec'), + ... in_bval=str(data_dir / 'bval')).run() + >>> check.outputs.pole + (0.0, 0.0, 0.0) + >>> check.outputs.full_sphere + True + >>> newrasb = np.loadtxt(check.outputs.out_rasb, skiprows=1) + >>> oldrasb = np.loadtxt(str(data_dir / 'dwi.tsv'), skiprows=1) + >>> np.allclose(newrasb, oldrasb, rtol=1.e-3) + True + + """ + + input_spec = _CheckGradientTableInputSpec + output_spec = _CheckGradientTableOutputSpec def _run_interface(self, runtime): - from dmriprep.utils.vectors import VectorTools - if isdefined(self.inputs.fbval) and isdefined(self.inputs.fbvec): - fbval = self.inputs.fbval - fbvec = self.inputs.fbvec - else: - fbval = None - fbvec = None - - if isdefined(self.inputs.rasb_tsv_file): - rasb_tsv_file = self.inputs.rasb_tsv_file - else: - rasb_tsv_file = None - - if isdefined(self.inputs.basedir): - basedir = self.inputs.basedir - else: - basedir = runtime.cwd - - vt = VectorTools(basedir, self.inputs.dwi_file, fbval, fbvec, rasb_tsv_file, self.inputs.B0_NORM_EPSILON, - self.inputs.bval_scaling) - - # Read in vectors - if rasb_tsv_file is not None: - vt.read_rasb_tsv() - elif fbval is not None and fbvec is not None: - vt.read_bvals_bvecs() - else: - raise OSError('VectorTools methods require either fbval/fbvec input files or an input rasb .tsv file.') - - # Rescale vectors - if self.inputs.rescale is True: - vt.rescale_vecs() - - # Build gradient table - gtab, b0_ixs = vt.build_gradient_table() - - # Check vector integrity - is_hemi, pole, slm = vt.checkvecs() - - if self.inputs.save_fsl_style is True: - fbval_out_file, fbvec_out_file = vt.save_vecs_fsl() - else: - fbval_out_file = None - fbvec_out_file = None - - self._results['gtab'] = gtab - self._results['is_hemi'] = is_hemi - self._results['pole'] = pole - self._results['slm'] = slm - self._results['b0_ixs'] = b0_ixs - self._results['fbval_out_file'] = fbval_out_file - self._results['fbvec_out_file'] = fbvec_out_file - self._results['rasb_tsv_out_file'] = vt.write_rasb_tsv() + rasb_file = _undefined(self.inputs, 'in_rasb') + + table = DiffusionGradientTable( + self.inputs.dwi_file, + bvecs=_undefined(self.inputs, 'in_bvec'), + bvals=_undefined(self.inputs, 'in_bval'), + rasb_file=rasb_file, + b_scale=self.inputs.b_scale, + bvec_norm_epsilon=self.inputs.bvec_norm_epsilon, + b0_threshold=self.inputs.b0_threshold, + ) + pole = table.pole + self._results['pole'] = tuple(pole) + self._results['full_sphere'] = np.all(pole == 0.0) + self._results['b0_ixs'] = np.where(table.b0mask)[0].tolist() + + if rasb_file is None: + rasb_file = fname_presuffix( + self.inputs.dwi_file, use_ext=False, suffix='.tsv', + newpath=runtime.cwd) + table.to_filename(rasb_file) + self._results['out_rasb'] = rasb_file + + cwd = Path(runtime.cwd).absolute() + + table.to_filename(bvecs=cwd / 'bvec', bvals=cwd / 'bval') + self._results['out_bval'] = str(cwd / 'bval') + self._results['out_bvec'] = str(cwd / 'bvec') return runtime + + +def _undefined(objekt, name, default=None): + value = getattr(objekt, name) + if not isdefined(value): + return default + return value diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 30fd4cc8..ba720ecd 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -1,402 +1,346 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -import os -import nibabel as nib +"""Utilities to operate on diffusion gradients.""" +from pathlib import Path +from itertools import permutations +import nibabel as nb import numpy as np - - -def make_gradient_table(bvals, bvecs, B0_NORM_EPSILON): +from dipy.core.gradients import round_bvals + +B0_THRESHOLD = 50 +BVEC_NORM_EPSILON = 0.1 + + +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 + self.generate_vecval() + elif dwi_file and bvecs and bvals: + self.bvecs = bvecs + self.bvals = bvals + self.generate_rasb() + + @property + def affine(self): + return self._affine + + @property + def gradients(self): + """Get gradient table (rasb).""" + return self._gradients + + @property + def bvecs(self): + return self._bvecs + + @property + def bvals(self): + return self._bvals + + @property + def normalized(self): + 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) + self._bvecs = value + + @bvals.setter + def bvals(self, value): + if isinstance(value, (str, Path)): + value = np.loadtxt(str(value)).flatten() + self._bvals = np.array(value) + + def b0mask(self): + """Get a mask of low-b frames.""" + return np.squeeze(self.gradients[..., -1] > self._b0_thres) + + def normalize(self): + 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): + if not self.bvecs: + 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, path=None, bvecs=None, bvals=None): + if path: + np.savetxt(str(path), self.gradients, + delimiter='\t', header='\t'.join('RASB'), + fmt=['%.8f'] * 3 + ['%g']) + if bvecs: + np.savetxt(str(bvecs), self.bvecs.T, fmt='%.6f') + if bvals: + np.savetxt(str(bvals), self.bvals, fmt='%.6f') + + +def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, + bvec_norm_epsilon=BVEC_NORM_EPSILON, b_scale=True): """ - Create gradient table from bval/bvec. + Normalize b-vectors and b-values. + + The resulting b-vectors will be of unit length for the non-zero b-values. + The resultinb b-values will be normalized by the square of the + corresponding vector amplitude. Parameters ---------- + bvec : m x n 2d array + Raw b-vectors array. bvals : 1d array Raw b-values float array. - bvecs : m x n 2d array - Raw b-vectors array. - B0_NORM_EPSILON : float + b0_threshold : float Gradient threshold below which volumes and vectors are considered B0's. Returns ------- - gtab : Obj - DiPy object storing diffusion gradient information. - b0_ixs : 1d array - Array containing the indices of the b0s. - """ - from dipy.core.gradients import gradient_table - from dmriprep.utils.vectors import check_corruption - - # Perform initial corruption check, correcting for obvious format/character encoding anomalies. - bvecs, bvals = check_corruption(bvecs, bvals, B0_NORM_EPSILON) - - # Creating the gradient table - gtab = gradient_table(bvals, bvecs, b0_threshold=B0_NORM_EPSILON) - - b0_ixs = np.where(gtab.b0s_mask == True)[0] - - return gtab, b0_ixs - - -def read_gradients_tsv(rasb_tsv_file): - """ - Read bvals/bvecs from a RAS-B tsv file. - - Parameters - ---------- - rasb_tsv_file : pathlike or str - Path to .tsv file containing RAS-B gradient table. - - Returns - ------- - bvals : 1d array - Raw b-values float array. - bvecs : m x n 2d array - Raw b-vectors array. - """ - RASB_table = np.genfromtxt(rasb_tsv_file, delimiter="\t") - bvecs = RASB_table[1:, 0:3] - bvals = RASB_table[1:, 3] - return bvals, bvecs - - -def write_gradients_tsv(bvecs, bvals, out_file): - """ - Save bvals/bveccs as RAS-B tsv. - - Parameters - ---------- - bvals : 1d array - Raw b-values float array. bvecs : m x n 2d array - Raw b-vectors array. - out_file : pathlike or str - Path to .tsv file containing RAS-B gradient table. - """ - RASB_table = np.column_stack([bvecs, bvals]) - with open(out_file, 'wb') as f: - f.write(b'R\tA\tS\tB\n') - np.savetxt(f, RASB_table.astype('float'), fmt='%.8f', delimiter="\t") - f.close() - return + Unit-normed b-vectors array. + bvals : 1d int array + Vector amplitude square normed b-values array. + Examples + -------- + >>> bvecs = np.vstack((np.zeros(3), 2.0 * np.eye(3), -0.8 * np.eye(3), np.ones(3))) + >>> bvals = np.array([1000] * bvecs.shape[0]) + >>> normalize_gradients(bvecs, bvals, 50) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ValueError: + + >>> bvals[0] = 0.0 + >>> norm_vecs, norm_vals = normalize_gradients(bvecs, bvals) + >>> np.all(norm_vecs[0] == 0) + True + + >>> norm_vecs[1, ...].tolist() + [1.0, 0.0, 0.0] + + >>> norm_vals[0] + 0 + >>> norm_vals[1] + 4000 + >>> norm_vals[-2] + 600 + >>> norm_vals[-1] + 3000 -def save_corrected_bval_bvec(vector_dir, bvec_corrected, bval_corrected, bval_scaling): """ - Normalize b-vectors to be of unit length for the non-zero b-values. - If the b-vector row reflects a B0, the vector is untouched. + bvals = np.array(bvals, dtype='float32') + bvecs = np.array(bvecs, dtype='float32') - Parameters - ---------- - vector_dir : str - Path the output vector directory to save the corrected bval/bvec files. - bvec_corrected : m x n 2d array - Corrected b-vectors array. - bval_corrected : 1d array - Corrected b-values float array. - bval_scaling : bool - If True, then normalizes b-val by the square of the vector amplitude. - - Returns - ------- - fbval_corrected : pathlike or str - File path of the corrected b-values. - fbvec_corrected : pathlike or str - File path of the corrected b-vectors. - """ - from pathlib import Path - - fbvec_corrected = str(Path(vector_dir)/'bvec_resc_corr.bvec') - np.savetxt(fbvec_corrected, bvec_corrected.T.astype('float'), fmt='%.8f') - - if bval_scaling is True: - fbval_corrected = str(Path(vector_dir) / 'bval_resc_corr.bval') - else: - fbval_corrected = str(Path(vector_dir) / 'bval_corr.bval') - np.savetxt(fbval_corrected, bval_corrected.T.astype('int'), fmt='%d', newline=' ') - return fbval_corrected, fbvec_corrected + b0s = bvals < b0_threshold + b0_vecs = np.linalg.norm(bvecs, axis=1) < bvec_norm_epsilon + # Check for bval-bvec discrepancy. + if not np.all(b0s == b0_vecs): + raise ValueError( + 'Inconsistent bvals and bvecs (%d, %d low-b, respectively).' % + (b0s.sum(), b0_vecs.sum())) -def rescale_bvec(bvec, B0_NORM_EPSILON): - """ - Normalize b-vectors to be of unit length for the non-zero b-values. - If the b-vector row reflects a B0, the vector is untouched. + # Rescale b-vals if requested + if b_scale: + bvals[~b0s] *= np.linalg.norm(bvecs[~b0s], axis=1) ** 2 - Parameters - ---------- - bvec : m x n 2d array - Raw b-vectors array. - B0_NORM_EPSILON : float - Gradient threshold below which volumes and vectors are considered B0's. + # Ensure b0s have (0, 0, 0) vectors + bvecs[b0s, :3] = np.zeros(3) - Returns - ------- - bvec : m x n 2d array - Unit-normed b-vectors array. - """ - axis = int(bvec.shape[0] == 3) - b0s = np.linalg.norm(bvec) < B0_NORM_EPSILON + # Round bvals + bvals = round_bvals(bvals) # Rescale b-vecs, skipping b0's, on the appropriate axis to unit-norm length. - bvec[~b0s] = bvec[~b0s] / np.linalg.norm(bvec[~b0s], axis=axis) - return bvec + bvecs[~b0s] /= np.linalg.norm(bvecs[~b0s], axis=1)[..., np.newaxis] + return bvecs, bvals.astype('uint16') -def rescale_bval(bval, bvec, B0_NORM_EPSILON): +def calculate_pole(bvecs, bvec_norm_epsilon=BVEC_NORM_EPSILON): """ - Normalize b-val by the square of the vector amplitude. - If the b-vector row reflects a B0, the corresponding b-val entry is untouched. + Check whether the b-vecs cover a hemisphere, and if so, calculate the pole. Parameters ---------- - bval : 1d array - Raw b-values float array. - bvec : m x n 2d array - Raw b-vectors array. - B0_NORM_EPSILON : float - Gradient threshold below which volumes and vectors are considered B0's. - - Returns - ------- - bval : 1d int array - Vector amplitude square normed b-values array. - """ - axis = int(bvec.shape[0] == 3) - b0s_loc = np.where(bval < B0_NORM_EPSILON) - - # Iterate through b-vec rows to perform amplitude correction to each corresponding b-val entry. - bval = bval.copy() - for dvol in range(len(bvec)): - if dvol not in b0s_loc: - bval[dvol] = bval[dvol] * np.linalg.norm(bvec[dvol], axis=axis)**2 - return np.round(bval, 0) - - -def is_hemispherical(vecs): - """ - Determine whether all points on a unit sphere lie in the same hemisphere. - - Parameters - ---------- - vecs : numpy.ndarray + bvecs : numpy.ndarray 2D numpy array with shape (N, 3) where N is the number of points. All points must lie on the unit sphere. Returns ------- - is_hemi : bool - If True, one can find a hemisphere that contains all the points. - If False, then the points do not lie in any hemisphere pole : numpy.ndarray - If `is_hemi == True`, then pole is the "central" pole of the - input vectors. Otherwise, pole is the zero vector. + A zero-vector if ``bvecs`` covers the full sphere, and the unit vector + locating the hemisphere pole othewise. + + Examples + -------- + >>> bvecs = [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), + ... (0.0, 1.0, 0.0), (0.0, -1.0, 0.0), + ... (0.0, 0.0, 1.0)] # Just half a sphere + >>> calculate_pole(bvecs).tolist() + [0.0, 0.0, 1.0] + + >>> bvecs += [(0.0, 0.0, -1.0)] # Make it a full-sphere + >>> calculate_pole(bvecs).tolist() + [0.0, 0.0, 0.0] References ---------- https://rstudio-pubs-static.s3.amazonaws.com/27121_a22e51b47c544980bad594d5e0bb2d04.html + """ - import itertools - if vecs.shape[1] != 3: - raise ValueError("Input vectors must be 3D vectors") - if not np.allclose(1, np.linalg.norm(vecs, axis=1)): - raise ValueError("Input vectors must be unit vectors") + bvecs = np.array(bvecs, dtype='float32') # Normalize inputs + b0s = np.linalg.norm(bvecs, axis=1) < bvec_norm_epsilon + bvecs = bvecs[~b0s] # Generate all pairwise cross products. - v0, v1 = zip(*[p for p in itertools.permutations(vecs, 2)]) - cross_prods = np.cross(v0, v1) + pairs = np.array(list(permutations(bvecs, 2))) + pairs = np.swapaxes(pairs, 0, 1) + cross_prods = np.cross(pairs[0, ...], pairs[1, ...]) # Normalize them. - cross_prods /= np.linalg.norm(cross_prods, axis=1)[:, np.newaxis] - - # `cross_prods` now contains all candidate vertex points for "the polygon" in the reference. "The polygon" is a - # subset. Find which points belong to the polygon using a dot product test with each of the original vectors. - angles = np.arccos(np.dot(cross_prods, vecs.transpose())) + cross_norms = np.linalg.norm(cross_prods, axis=1) + cross_zeros = cross_norms < 1.0e-4 + cross_prods = cross_prods[~cross_zeros] + cross_prods /= cross_norms[~cross_zeros, np.newaxis] + + # `cross_prods` now contains all candidate vertex points for "the polygon" + # in the reference. + # "The polygon" is a subset. + # Find which points belong to the polygon using a dot product test + # with each of the original vectors. + angles = np.arccos(cross_prods.dot(bvecs.T)) # And test whether it is orthogonal or less. dot_prod_test = angles <= np.pi / 2.0 + ntests = dot_prod_test.sum(axis=1) == len(bvecs) - # If there is at least one point that is orthogonal or less to each input vector, then the points lie on some - # hemisphere. - is_hemi = len(vecs) in np.sum(dot_prod_test.astype(int), axis=1) + # If there is at least one point that is orthogonal or less to each + # input vector, then the points lie on some hemisphere. + is_hemi = np.any(ntests) + pole = np.zeros(3) if is_hemi: - vertices = cross_prods[np.sum(dot_prod_test.astype(int), axis=1) == len(vecs)] + vertices = cross_prods[ntests] pole = np.mean(vertices, axis=0) pole /= np.linalg.norm(pole) - else: - pole = np.array([0.0, 0.0, 0.0]) - return is_hemi, pole - - -def set_slm(gtab, is_hemi): - """ - Evaluate eddy compatibility and set second level model that specifies the mathematical form for how the - diffusion gradients cause eddy currents. If the data has few directions and/or is has not been sampled on the whole - sphere it can be advantageous to specify --slm=linear. + return pole - Parameters - ---------- - gtab : Obj - DiPy object storing diffusion gradient information. - is_hemi : bool - If True, one can find a hemisphere that contains all the points. - If False, then the points do not lie in any hemisphere. - Returns - ------- - slm : str - Second-level model to use for eddy correction. +def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): """ - # Warn user if gradient table is not conducive to eddy correction. - if ((np.sum(np.invert(gtab.b0s_mask)) < 10) and (max(gtab.bvals) < 1500)) or ((np.sum(np.invert(gtab.b0s_mask)) < 30) and (max(gtab.bvals) < 5000)): - raise UserWarning('Too few directions in this data. Use of eddy is not recommended.') - - if is_hemi is True: - slm = 'linear' - print("B-vectors for this data are hemispherical. To ensure adequate eddy current correction, eddy will be run " - "using the --slm=linear option.") - elif np.sum(np.invert(gtab.b0s_mask)) < 30: - slm = 'linear' - print('Fewer than 30 diffusion-weighted directions detected. To ensure adequate eddy current correction, eddy ' - 'will be run using the --slm=linear option') - else: - slm = 'none' - return slm - - -def check_corruption(bvecs, bvals, B0_NORM_EPSILON): - """ - Performs a variety of basic check to ensure bval/bvec files are not corrupted. + Convert b-vectors given in image coordinates to RAS+. - Parameters - ---------- - bvecs : m x n 2d array - Raw b-vectors array. - bvals : 1d array - Raw b-values float array. - B0_NORM_EPSILON : float - Gradient threshold below which volumes and vectors are considered B0's. + Examples + -------- + >>> bvecs2ras(2.0 * np.eye(3), [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)]).tolist() + [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0]] - Returns - ------- - bvecs : m x n 2d array - Raw b-vectors array. - bvals : 1d array - Raw b-values float array. - """ - # Correct any b0's in bvecs misstated as 10's. - bvecs[np.where(np.any(abs(bvecs) >= 10, axis=1) == True)] = [1, 0, 0] - - # Check for bval-bvec discrepancy. - if len(bvecs[np.where(np.logical_and(bvals > B0_NORM_EPSILON, np.where(np.all(abs(bvecs)) == 0)))]) > 0: - raise ValueError("Corrupted bval/bvecs detected. Check to ensure volumes with a diffusion weighting are not " - "listed as B0s in the bvecs.") + >>> affine = np.eye(4) + >>> affine[0, 0] *= -1.0 # Make it LAS + >>> bvecs2ras(affine, [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)]).tolist() + [[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]] - # Ensure that B0's in the bvecs are represented are zero floats. - bvecs[np.where(bvals < B0_NORM_EPSILON)] = float(0) + >>> affine = np.eye(3) + >>> affine[:2, :2] *= -1.0 # Make it LPS + >>> bvecs2ras(affine, [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)]).tolist() + [[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]] - return bvecs, bvals + >>> affine[:2, :2] = [[0.0, 1.0], [1.0, 0.0]] # Make it ARS + >>> bvecs2ras(affine, [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)]).tolist() + [[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]] + >>> bvecs2ras(affine, [(0.0, 0.0, 0.0)]).tolist() + [[0.0, 0.0, 0.0]] -def rescale_vectors(bvals, bvecs, bval_scaling, bmag, B0_NORM_EPSILON): """ - Performs a variety of rescaling routines on the b vectors and values. - - Parameters - ---------- - bvals : 1d array - Raw b-values float array. - bvecs : m x n 2d array - Raw b-vectors array. - bval_scaling : bool - If True, then normalizes b-val by the square of the vector amplitude. - bmag : int - The order of magnitude to round the b-values. - B0_NORM_EPSILON : float - Gradient threshold below which volumes and vectors are considered B0's. - - Returns - ------- - bval_mag_norm : 1d array - Rescaled b-values float array. - bvec_rescaled : m x n 2d array - Rescaled b-vectors array. - """ - from dipy.core.gradients import round_bvals - from dmriprep.utils.vectors import rescale_bvec, rescale_bval - # Rescale bvals by vector amplitudes. - if bval_scaling is True: - bvals = rescale_bval(bvals.astype('float'), bvecs.astype('float'), B0_NORM_EPSILON) - - # Round bvals by bmag. - bval_mag_norm = round_bvals(bvals, bmag=bmag) - - # Rescale bvecs. - bvec_rescaled = rescale_bvec(bvecs, B0_NORM_EPSILON) - return bval_mag_norm, bvec_rescaled - - -class VectorTools(object): - """ - Class of vector-handling methods for io and integrity checking. - - Parameters - ---------- - basedir : str - Path the output base directory for derivative dmriprep data. - dwi_file : str - Optionally provide a file path to the diffusion-weighted image series to which the - bvecs/bvals should correspond. - fbval : str - File path of the b-values. - fbvec : str - File path of the b-vectors. - rasb_tsv_file : pathlike or str - Path to .tsv file containing RAS-B gradient table. - B0_NORM_EPSILON : float - Gradient threshold below which volumes and vectors are considered B0's. Default is 50. - bval_scaling : bool - If True, then normalizes b-val by the square of the vector amplitude. Default is True. - """ - from pathlib import Path - from dipy.io import read_bvals_bvecs - from dmriprep.utils.vectors import make_gradient_table, is_hemispherical, check_corruption, save_corrected_bval_bvec, make_gradients_tsv, set_slm - - # Create vector directory if it does not already exist. - vector_dir = str(Path(sesdir)/'vectors') - os.makedirs(vector_dir, exist_ok=True) - - # load bvecs/bvals. - bvals, bvecs = read_bvals_bvecs(fbval, fbvec) - - # Get b order of magnitude - bmag = int(np.log10(np.max(bvals))) - 1 - - # Perform initial corruption check, correcting for obvious anomalies. - [bvecs, bvals] = check_corruption(bvecs, bvals, B0_NORM_EPSILON) - - # Perform b vector/value rescaling. - [bval_mag_norm, bvec_rescaled] = rescale_vectors(bvals, bvecs, bval_scaling, bmag, B0_NORM_EPSILON) - - # Save corrected bval/bvec to file in FSL style to ensure-backwards compatibility. - [fbval_rescaled, fbvec_rescaled] = save_corrected_bval_bvec(vector_dir, bvec_rescaled, bval_mag_norm, bval_scaling) - - # Make gradient table. - [gtab, b0_ixs] = make_gradient_table(fbval_rescaled, fbvec_rescaled, B0_NORM_EPSILON) - - # Save gradient table to tsv. - gtab_tsv_file = str(Path(vector_dir)/'gradients.tsv') - make_gradients_tsv(gtab, gtab_tsv_file) - - # Check consistency of dwi image and gradient vectors. - image_gradient_consistency_check(dwi_file, b0_ixs, gtab, len(np.unique(gtab.bvals)), bmag) - - # Check hemisphere coverage. - [is_hemi, pole] = is_hemispherical(np.round(gtab.bvecs, 8)[~(np.round(gtab.bvecs, 8) == 0).all(1)]) - - # Use sampling scheme to set second-level model for eddy correction (this step would conceivably apply to any - # non-FSL form of eddy correction as well. - slm = set_slm(gtab, is_hemi) - - return gtab, is_hemi, pole, b0_ixs, gtab_tsv_file, slm + if affine.shape == (4, 4): + affine = affine[:3, :3] + + bvecs = np.array(bvecs, dtype='float32') # Normalize inputs + rotated_bvecs = affine[np.newaxis, ...].dot(bvecs.T)[0].T + norms_bvecs = np.linalg.norm(rotated_bvecs, axis=1) + b0s = norms_bvecs < bvec_norm_epsilon + rotated_bvecs[~b0s] /= norms_bvecs[~b0s, np.newaxis] + rotated_bvecs[b0s] = np.zeros(3) + return rotated_bvecs From 3be002e93ff9723aa65c3f1a6ed182dd5a42045f Mon Sep 17 00:00:00 2001 From: dPys Date: Sun, 3 Nov 2019 13:51:45 -0600 Subject: [PATCH 03/12] correct condition on bvecs2ras --- dmriprep/utils/vectors.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index ba720ecd..39dcda0d 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -34,7 +34,6 @@ def __init__(self, dwi_file=None, bvecs=None, bvals=None, rasb_file=None, 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 @@ -169,7 +168,7 @@ def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, Parameters ---------- - bvec : m x n 2d array + bvecs : m x n 2d array Raw b-vectors array. bvals : 1d array Raw b-values float array. @@ -339,8 +338,9 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): bvecs = np.array(bvecs, dtype='float32') # Normalize inputs rotated_bvecs = affine[np.newaxis, ...].dot(bvecs.T)[0].T - norms_bvecs = np.linalg.norm(rotated_bvecs, axis=1) - b0s = norms_bvecs < bvec_norm_epsilon - rotated_bvecs[~b0s] /= norms_bvecs[~b0s, np.newaxis] - rotated_bvecs[b0s] = np.zeros(3) + if norm is True: + norms_bvecs = np.linalg.norm(rotated_bvecs, axis=1) + b0s = norms_bvecs < bvec_norm_epsilon + rotated_bvecs[~b0s] /= norms_bvecs[~b0s, np.newaxis] + rotated_bvecs[b0s] = np.zeros(3) return rotated_bvecs From 374b4eefffc7327d3e290bdafa8e336fe8c4797a Mon Sep 17 00:00:00 2001 From: oesteban Date: Sun, 3 Nov 2019 12:04:34 -0800 Subject: [PATCH 04/12] fix: add dipy dependency (it has to be >= 1.0.0 to use round_bvals) --- setup.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.cfg b/setup.cfg index f46dcfbc..af763bc3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -20,6 +20,7 @@ classifiers = [options] python_requires = >=3.5 install_requires = + dipy >= 1.0.0 indexed_gzip >=0.8.8 nibabel >=2.2.1 nilearn !=0.5.0, !=0.5.1 From 711d357d65d66b29a450272373812628a55877d6 Mon Sep 17 00:00:00 2001 From: oesteban Date: Sun, 3 Nov 2019 22:54:56 -0800 Subject: [PATCH 05/12] enh: revise tests, add checks to bval/bvecs setter --- dmriprep/conftest.py | 24 +++ dmriprep/utils/tests/test_vectors.py | 220 ++++++--------------------- dmriprep/utils/vectors.py | 9 +- 3 files changed, 76 insertions(+), 177 deletions(-) diff --git a/dmriprep/conftest.py b/dmriprep/conftest.py index fdb523d5..cb59a51d 100644 --- a/dmriprep/conftest.py +++ b/dmriprep/conftest.py @@ -19,3 +19,27 @@ def doctest_autoimport(doctest_namespace): doctest_namespace['tmpdir'] = tmpdir.name yield tmpdir.cleanup() + + +@pytest.fixture(scope="session") +def dipy_test_data(tmpdir_factory): + """Create a temporal directory shared across tests to pull data in.""" + from dipy.data.fetcher import _make_fetcher, UW_RW_URL + + datadir = Path(tmpdir_factory.mktemp("data")) + _make_fetcher( + "fetch_sherbrooke_3shell", + str(datadir), + UW_RW_URL + "1773/38475/", + ['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'], + ['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'], + ['0b735e8f16695a37bfbd66aab136eb66', + 'e9b9bb56252503ea49d31fb30a0ac637', + '0c83f7e8b917cd677ad58a078658ebb7'], + doc="Download a 3shell HARDI dataset with 192 gradient direction")() + + return { + 'dwi_file': str(datadir / "HARDI193.nii.gz"), + 'bvecs': str(datadir / "HARDI193.bvec"), + 'bvals': str(datadir / "HARDI193.bval"), + } diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index 97968697..b9655d6a 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -1,192 +1,62 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -import os +"""Test vector utilities.""" import pytest import numpy as np -from pathlib import Path -from dmriprep.utils import vectors +from dmriprep.utils import vectors as v -def test_make_and_save_gradient_table(): - """Test for gradient table generation and saving from bval/bvec.""" - from dipy.data import fetch_sherbrooke_3shell - from dipy.io import read_bvals_bvecs - from os.path import expanduser, join - from dipy.core.gradients import GradientTable +def test_corruption(dipy_test_data): + """Check whether b-value rescaling is operational.""" - fetch_sherbrooke_3shell() - home = expanduser('~') - basedir = join(home, '.dipy', 'sherbrooke_3shell') - fdwi = join(basedir, 'HARDI193.nii.gz') - fbval = join(basedir, 'HARDI193.bval') - fbvec = join(basedir, 'HARDI193.bvec') + bvecs = np.loadtxt(dipy_test_data['bvecs']).T + bvals = np.loadtxt(dipy_test_data['bvals']) - # load bvecs/bvals. - bvals, bvecs = read_bvals_bvecs(fbval, fbvec) - - b0_thr = 50 - # Create vector directory if it does not already exist. - vector_dir = str(Path(basedir)/'vectors') - os.makedirs(vector_dir, exist_ok=True) - gtab_tsv_file = str(Path(vector_dir) / 'gradients.tsv') - - [gtab, b0_ixs] = vectors.make_gradient_table(bvals, bvecs, b0_thr) - vectors.write_gradients_tsv(bvecs, bvals, gtab_tsv_file) - - [bvals_read, bvecs_read] = vectors.read_gradients_tsv(gtab_tsv_file) - - assert gtab is not None and gtab.__class__ == GradientTable - assert gtab.b0_threshold == 50 - assert np.sum(gtab.b0s_mask) == len(b0_ixs) - assert os.path.isfile(gtab_tsv_file) is True - assert np.isclose(bvecs_read.all(), bvecs.all()) - assert np.isclose(bvals_read.all(), bvals.all()) - os.remove(gtab_tsv_file) - - -def test_rescale_vectors(): - """Test vector rescaling.""" - from dipy.data import fetch_sherbrooke_3shell - from os.path import expanduser, join - from dipy.io import read_bvals_bvecs - - fetch_sherbrooke_3shell() - home = expanduser('~') - basedir = join(home, '.dipy', 'sherbrooke_3shell') - fbval = join(basedir, 'HARDI193.bval') - fbvec = join(basedir, 'HARDI193.bvec') - - bval_scaling = True - b0_thr = 50 - - # load bvecs/bvals. - bvals, bvecs = read_bvals_bvecs(fbval, fbvec) - - # Get b order of magnitude - bmag = int(np.log10(np.max(bvals))) - 1 - - [bval_normed, bvec_normed] = vectors.rescale_vectors(bvals, bvecs, bval_scaling, bmag, b0_thr) - assert np.max(bvec_normed) < 1.0 and np.min(bvec_normed) > -1.0 - assert len(bvecs) == len(bvec_normed) - assert np.sum(np.count_nonzero(bvecs, axis=1).astype('bool')) == np.sum(np.count_nonzero(bvec_normed, - axis=1).astype('bool')) - assert len(bvals) == len(bval_normed) - assert np.sum(np.count_nonzero(bvals, axis=0).astype('bool')) == np.sum(np.count_nonzero(bval_normed, - axis=0).astype('bool')) - - -def test_vector_checkers(): - """Test for vector corruption.""" - from dipy.data import fetch_sherbrooke_3shell - from os.path import expanduser, join - from dipy.io import read_bvals_bvecs - import itertools - - fetch_sherbrooke_3shell() - home = expanduser('~') - basedir = join(home, '.dipy', 'sherbrooke_3shell') - fbval = join(basedir, 'HARDI193.bval') - fbvec = join(basedir, 'HARDI193.bvec') - - b0_thr = 50 + # Perform various corruption checks using synthetic corrupted bval-bvec. + dgt = v.DiffusionGradientTable() + dgt.bvecs = bvecs + with pytest.raises(ValueError): + dgt.bvals = bvals[:-1] - # load bvecs/bvals. - bvals, bvecs = read_bvals_bvecs(fbval, fbvec) + dgt = v.DiffusionGradientTable() + dgt.bvals = bvals + with pytest.raises(ValueError): + dgt.bvecs = bvecs[:-1] - # Perform various corruption checks using synthetic corrupted bval-bvec. - bval_short = bvals[:-1] - bval_long = np.hstack([bvals, bvals[-1]]) + # Missing b0 bval_no_b0 = bvals.copy() bval_no_b0[0] = 51 - bval_odd_b0 = bvals.copy() - bval_odd_b0[np.where(abs(bval_odd_b0) == 0)] = 0.00000001 - bvec_short = bvecs[:-1] - bvec_long = np.vstack([bvecs, 0.5*bvecs[-1]]) + with pytest.raises(ValueError): + dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'], + bvals=bval_no_b0, bvecs=bvecs) bvec_no_b0 = bvecs.copy() - bvec_no_b0[np.where(np.all(abs(bvec_no_b0) == 0, axis=1) == True)] = [1, 1, 1] - bvec_odd_b0 = bvecs.copy() - bvec_odd_b0[np.where(np.all(abs(bvec_odd_b0) == 0, axis=1) == True)] = [10, 10, 10] - - bval_lists = [bvals, bval_short, bval_long, bval_no_b0, bval_odd_b0] - bvec_lists = [bvecs, bvec_short, bvec_long, bvec_no_b0, bvec_odd_b0] - for bval, bvec in list(itertools.product(bval_lists, bvec_lists)): - [bvecs_checked, bvals_checked] = vectors.check_corruption(bvec, bval, b0_thr) - if (bval is bvals) and (bvec is bvecs): - assert bvecs_checked is not None - assert bvals_checked is not None - else: - with pytest.raises(Exception) as e_info: - assert bvecs_checked is None - assert bvals_checked is None - - -def test_hemisphericity(): - """Test vector hemisphere coverage and second-level model setting based on that coverage.""" - from dipy.data import fetch_sherbrooke_3shell - from dipy.io import read_bvals_bvecs - from os.path import expanduser, join - - fetch_sherbrooke_3shell() - home = expanduser('~') - basedir = join(home, '.dipy', 'sherbrooke_3shell') - fbval = join(basedir, 'HARDI193.bval') - fbvec = join(basedir, 'HARDI193.bvec') - - b0_thr = 50 - - # load bvecs/bvals. - bvals, bvecs = read_bvals_bvecs(fbval, fbvec) - [gtab, _] = vectors.make_gradient_table(bvals, bvecs, b0_thr) - - # Check hemisphere coverage. - [is_hemi, pole] = vectors.is_hemispherical(np.round(gtab.bvecs, 8)[~(np.round(gtab.bvecs, 8) == 0).all(1)]) - - # Use sampling scheme to set second-level model for eddy correction (this step would conceivably apply to any - # non-FSL form of eddy correction as well. - slm = vectors.set_slm(gtab, is_hemi) - assert is_hemi is False - assert isinstance(pole, (np.ndarray, np.generic)) is True - assert slm == 'none' + bvec_no_b0[0] = np.array([1.0, 0.0, 0.0]) + with pytest.raises(ValueError): + dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'], + bvals=bvals, bvecs=bvec_no_b0) + # Corrupt b0 b-val + bval_odd_b0 = bvals.copy() + bval_odd_b0[bval_odd_b0 == 0] = 1e-8 + dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'], + bvals=bval_odd_b0, bvecs=bvecs) + assert dgt.bvals[0] == 0 -def test_clean_vecs(): - """Test CleanVecs interface functionality.""" - from dipy.data import fetch_sherbrooke_3shell - from os.path import expanduser, join - from dipy.core.gradients import GradientTable - from dmriprep.interfaces.vectors import CleanVecs - import itertools + # Corrupt b0 b-vec + bvec_odd_b0 = bvecs.copy() + b0mask = np.all(bvec_odd_b0 == 0, axis=1) + bvec_odd_b0[b0mask] = [10, 10, 10] + dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'], + bvals=bvals, bvecs=bvec_odd_b0) + assert np.all(dgt.bvecs[b0mask] == [0., 0., 0.]) - fetch_sherbrooke_3shell() - home = expanduser('~') - basedir = join(home, '.dipy', 'sherbrooke_3shell') - fdwi = join(basedir, 'HARDI193.nii.gz') - fbval = join(basedir, 'HARDI193.bval') - fbvec = join(basedir, 'HARDI193.bvec') + # Test normalization + bvecs_factor = 2.0 * bvecs + dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'], + bvals=bvals, bvecs=bvecs_factor) + assert -1.0 <= np.max(np.abs(dgt.gradients[..., :-1])) <= 1.0 - # Smoke test parameter combos - for bool_list in list(itertools.product([True, False], repeat=3)): - cv = CleanVecs() - cv.inputs.B0_NORM_EPSILON = 50 - cv.inputs.bval_scaling = bool_list[0] - cv.inputs.rescale = bool_list[1] - cv.inputs.save_fsl_style = bool_list[2] - cv.inputs.dwi_file = fdwi - cv.inputs.fbval = fbval - cv.inputs.fbvec = fbvec - cv.inputs.basedir = basedir - res = cv.run() +def test_hemisphericity(dipy_test_data): + """Test vector hemisphere coverage.""" - assert res.outputs.gtab is not None and res.outputs.gtab.__class__ == GradientTable - assert res.outputs.is_hemi is False - assert isinstance(res.outputs.pole, (np.ndarray, np.generic)) is True - assert np.sum(res.outputs.gtab.b0s_mask) == len(res.outputs.b0_ixs) - assert os.path.isfile(res.outputs.rasb_tsv_out_file) is True - if cv.inputs.save_fsl_style is True: - assert os.path.isfile(res.outputs.fbval_out_file) is True - assert os.path.isfile(res.outputs.fbvec_out_file) is True - os.remove(res.outputs.fbval_out_file) - os.remove(res.outputs.fbvec_out_file) - assert res.outputs.slm == 'none' - os.remove(res.outputs.rasb_tsv_out_file) + dgt = v.DiffusionGradientTable(**dipy_test_data) + assert np.all(dgt.pole == [0., 0., 0.]) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 39dcda0d..79bb24c1 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -48,8 +48,9 @@ def __init__(self, dwi_file=None, bvecs=None, bvals=None, rasb_file=None, self.affine = dwi_file if rasb_file is not None: self.gradients = rasb_file - self.generate_vecval() - elif dwi_file and bvecs and bvals: + if self.affine is not None: + self.generate_vecval() + elif dwi_file and bvecs is not None and bvals is not None: self.bvecs = bvecs self.bvals = bvals self.generate_rasb() @@ -100,12 +101,16 @@ def bvecs(self, value): # 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) def b0mask(self): From a5df719607d3cb2a4d2989388f130fc196cb828a Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 4 Nov 2019 08:37:36 -0800 Subject: [PATCH 06/12] fix: tests were failing for an obscure numpy problem --- .travis.yml | 1 + dmriprep/conftest.py | 47 +++++++++++++++------------- dmriprep/utils/tests/test_vectors.py | 6 ++-- 3 files changed, 29 insertions(+), 25 deletions(-) diff --git a/.travis.yml b/.travis.yml index 8dd653dd..647dc893 100644 --- a/.travis.yml +++ b/.travis.yml @@ -6,6 +6,7 @@ language: python cache: directories: - $HOME/.cache/pip + - $HOME/.cache/data python: - 3.6 diff --git a/dmriprep/conftest.py b/dmriprep/conftest.py index cb59a51d..ea35c7ba 100644 --- a/dmriprep/conftest.py +++ b/dmriprep/conftest.py @@ -5,6 +5,28 @@ import numpy as np import nibabel as nb import pytest +from dipy.data.fetcher import _make_fetcher, UW_RW_URL + +_dipy_datadir_root = os.getenv('DMRIPREP_TESTS_DATA') or Path.home() +dipy_datadir = Path(_dipy_datadir_root) / '.cache' / 'data' +dipy_datadir.mkdir(parents=True, exist_ok=True) + +_make_fetcher( + "fetch_sherbrooke_3shell", + str(dipy_datadir), + UW_RW_URL + "1773/38475/", + ['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'], + ['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'], + ['0b735e8f16695a37bfbd66aab136eb66', + 'e9b9bb56252503ea49d31fb30a0ac637', + '0c83f7e8b917cd677ad58a078658ebb7'], + doc="Download a 3shell HARDI dataset with 192 gradient direction")() + +_sherbrooke_data = { + 'dwi_file': dipy_datadir / "HARDI193.nii.gz", + 'bvecs': np.loadtxt(dipy_datadir / "HARDI193.bvec").T, + 'bvals': np.loadtxt(dipy_datadir / "HARDI193.bval"), +} @pytest.fixture(autouse=True) @@ -15,31 +37,14 @@ def doctest_autoimport(doctest_namespace): doctest_namespace['os'] = os doctest_namespace['Path'] = Path doctest_namespace['data_dir'] = Path(__file__).parent / 'data' / 'tests' + doctest_namespace['dipy_datadir'] = dipy_datadir tmpdir = tempfile.TemporaryDirectory() doctest_namespace['tmpdir'] = tmpdir.name yield tmpdir.cleanup() -@pytest.fixture(scope="session") -def dipy_test_data(tmpdir_factory): +@pytest.fixture() +def dipy_test_data(scope='session'): """Create a temporal directory shared across tests to pull data in.""" - from dipy.data.fetcher import _make_fetcher, UW_RW_URL - - datadir = Path(tmpdir_factory.mktemp("data")) - _make_fetcher( - "fetch_sherbrooke_3shell", - str(datadir), - UW_RW_URL + "1773/38475/", - ['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'], - ['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'], - ['0b735e8f16695a37bfbd66aab136eb66', - 'e9b9bb56252503ea49d31fb30a0ac637', - '0c83f7e8b917cd677ad58a078658ebb7'], - doc="Download a 3shell HARDI dataset with 192 gradient direction")() - - return { - 'dwi_file': str(datadir / "HARDI193.nii.gz"), - 'bvecs': str(datadir / "HARDI193.bvec"), - 'bvals': str(datadir / "HARDI193.bval"), - } + return _sherbrooke_data diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index b9655d6a..72573415 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -6,9 +6,8 @@ def test_corruption(dipy_test_data): """Check whether b-value rescaling is operational.""" - - bvecs = np.loadtxt(dipy_test_data['bvecs']).T - bvals = np.loadtxt(dipy_test_data['bvals']) + bvals = dipy_test_data['bvals'] + bvecs = dipy_test_data['bvecs'] # Perform various corruption checks using synthetic corrupted bval-bvec. dgt = v.DiffusionGradientTable() @@ -57,6 +56,5 @@ def test_corruption(dipy_test_data): def test_hemisphericity(dipy_test_data): """Test vector hemisphere coverage.""" - dgt = v.DiffusionGradientTable(**dipy_test_data) assert np.all(dgt.pole == [0., 0., 0.]) From 8596a2153a7b1d2a47d95a480d549b3785707c7e Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 4 Nov 2019 08:50:13 -0800 Subject: [PATCH 07/12] sty: add missing docstrings [skip ci] --- dmriprep/utils/vectors.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 79bb24c1..acd5b99e 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -34,6 +34,7 @@ def __init__(self, dwi_file=None, bvecs=None, bvals=None, rasb_file=None, 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 @@ -57,6 +58,7 @@ def __init__(self, dwi_file=None, bvecs=None, bvals=None, rasb_file=None, @property def affine(self): + """Get the affine for RAS+/image-coordinates conversions.""" return self._affine @property @@ -66,14 +68,17 @@ def gradients(self): @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 @@ -118,6 +123,7 @@ def b0mask(self): return np.squeeze(self.gradients[..., -1] > self._b0_thres) def normalize(self): + """Normalize (l2-norm) b-vectors.""" if self._normalized: return @@ -136,7 +142,8 @@ def generate_rasb(self): self.gradients = np.hstack((_ras, self.bvals[..., np.newaxis])) def generate_vecval(self): - if not self.bvecs: + """Compose a bvec/bval pair in image coordinates.""" + if self.bvecs is not None: self._bvecs = bvecs2ras(np.linalg.inv(self.affine), self.gradients[..., :-1]) self._bvals = self.gradients[..., -1].flatten() @@ -152,6 +159,7 @@ def pole(self): return calculate_pole(self.gradients[..., :-1], bvec_norm_epsilon=self._bvec_norm_epsilon) def to_filename(self, path=None, bvecs=None, bvals=None): + """Write files (RASB, bvecs/bvals) to a given path.""" if path: np.savetxt(str(path), self.gradients, delimiter='\t', header='\t'.join('RASB'), From 96a25c75856edd22b2c487b057d3fdf8b9973fe7 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 4 Nov 2019 12:38:06 -0800 Subject: [PATCH 08/12] fix: bad condition check in generate_vecval, improved to_filename API addresses https://github.com/nipreps/dmriprep/pull/26#discussion_r342240815 --- dmriprep/interfaces/vectors.py | 12 +++++------- dmriprep/utils/vectors.py | 17 +++++++++-------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/dmriprep/interfaces/vectors.py b/dmriprep/interfaces/vectors.py index 74c594a5..6d9da11b 100644 --- a/dmriprep/interfaces/vectors.py +++ b/dmriprep/interfaces/vectors.py @@ -79,18 +79,16 @@ def _run_interface(self, runtime): self._results['full_sphere'] = np.all(pole == 0.0) self._results['b0_ixs'] = np.where(table.b0mask)[0].tolist() + cwd = Path(runtime.cwd).absolute() if rasb_file is None: rasb_file = fname_presuffix( self.inputs.dwi_file, use_ext=False, suffix='.tsv', - newpath=runtime.cwd) + newpath=str(cwd)) table.to_filename(rasb_file) self._results['out_rasb'] = rasb_file - - cwd = Path(runtime.cwd).absolute() - - table.to_filename(bvecs=cwd / 'bvec', bvals=cwd / 'bval') - self._results['out_bval'] = str(cwd / 'bval') - self._results['out_bvec'] = str(cwd / 'bvec') + table.to_filename('%s/dwi' % cwd, filetype='fsl') + self._results['out_bval'] = str(cwd / 'dwi.bval') + self._results['out_bvec'] = str(cwd / 'dwi.bvec') return runtime diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index acd5b99e..6194ac61 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -143,7 +143,7 @@ def generate_rasb(self): def generate_vecval(self): """Compose a bvec/bval pair in image coordinates.""" - if self.bvecs is not None: + if self.bvecs is None or self.bvals is None: self._bvecs = bvecs2ras(np.linalg.inv(self.affine), self.gradients[..., :-1]) self._bvals = self.gradients[..., -1].flatten() @@ -158,16 +158,17 @@ def pole(self): self.generate_rasb() return calculate_pole(self.gradients[..., :-1], bvec_norm_epsilon=self._bvec_norm_epsilon) - def to_filename(self, path=None, bvecs=None, bvals=None): + def to_filename(self, filename, filetype='rasb'): """Write files (RASB, bvecs/bvals) to a given path.""" - if path: - np.savetxt(str(path), self.gradients, + if filetype == 'rasb': + np.savetxt(str(filename), self.gradients, delimiter='\t', header='\t'.join('RASB'), fmt=['%.8f'] * 3 + ['%g']) - if bvecs: - np.savetxt(str(bvecs), self.bvecs.T, fmt='%.6f') - if bvals: - np.savetxt(str(bvals), self.bvals, fmt='%.6f') + elif filetype == 'fsl': + 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) def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, From e4ec7e99e98cb85749a3ea7c0b20c5a866098219 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 4 Nov 2019 12:40:20 -0800 Subject: [PATCH 09/12] fix: reversed comparison addresses https://github.com/nipreps/dmriprep/pull/26#discussion_r342250064 --- dmriprep/utils/vectors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 6194ac61..0ff16583 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -120,7 +120,7 @@ def bvals(self, value): def b0mask(self): """Get a mask of low-b frames.""" - return np.squeeze(self.gradients[..., -1] > self._b0_thres) + return np.squeeze(self.gradients[..., -1] < self._b0_thres) def normalize(self): """Normalize (l2-norm) b-vectors.""" From 166e723526bd5e94536f2662519d6d71773f6c41 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 4 Nov 2019 12:41:06 -0800 Subject: [PATCH 10/12] Update dmriprep/utils/vectors.py Co-Authored-By: Michael Joseph --- dmriprep/utils/vectors.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index acd5b99e..e21307a6 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -118,6 +118,7 @@ def bvals(self, value): 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) From b5e34e242a0b76d915fcd385e16120a21db786f7 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 4 Nov 2019 13:27:20 -0800 Subject: [PATCH 11/12] enh: increase code coverage with more test cases --- dmriprep/utils/tests/test_vectors.py | 36 ++++++++++++++++++++++++---- dmriprep/utils/vectors.py | 20 ++++++++++++++-- 2 files changed, 49 insertions(+), 7 deletions(-) diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index 72573415..c2674cc4 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -2,13 +2,29 @@ import pytest import numpy as np from dmriprep.utils import vectors as v +from collections import namedtuple -def test_corruption(dipy_test_data): +def test_corruption(tmpdir, dipy_test_data, monkeypatch): """Check whether b-value rescaling is operational.""" + tmpdir.chdir() + bvals = dipy_test_data['bvals'] bvecs = dipy_test_data['bvecs'] + # Test vector hemisphere coverage + dgt = v.DiffusionGradientTable(**dipy_test_data) + assert np.all(dgt.pole == [0., 0., 0.]) + + dgt.to_filename('dwi.tsv') + dgt = v.DiffusionGradientTable(rasb_file='dwi.tsv') + assert dgt.normalized is False + with pytest.raises(TypeError): + dgt.to_filename('dwi', filetype='fsl') # You can do this iff the affine is set. + + aff = namedtuple('Affine', ['affine'])(dgt.affine) # check accessing obj.affine + dgt = v.DiffusionGradientTable(dwi_file=aff) + # Perform various corruption checks using synthetic corrupted bval-bvec. dgt = v.DiffusionGradientTable() dgt.bvecs = bvecs @@ -52,9 +68,19 @@ def test_corruption(dipy_test_data): dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'], bvals=bvals, bvecs=bvecs_factor) assert -1.0 <= np.max(np.abs(dgt.gradients[..., :-1])) <= 1.0 + assert dgt.normalized is True + def mock_func(*args, **kwargs): + return 'called!' -def test_hemisphericity(dipy_test_data): - """Test vector hemisphere coverage.""" - dgt = v.DiffusionGradientTable(**dipy_test_data) - assert np.all(dgt.pole == [0., 0., 0.]) + with monkeypatch.context() as m: + m.setattr(v, 'normalize_gradients', mock_func) + assert dgt.normalize() is None # Test nothing is executed. + + with monkeypatch.context() as m: + m.setattr(v, 'bvecs2ras', mock_func) + assert dgt.generate_vecval() is None # Test nothing is executed. + + # Miscellaneous tests + with pytest.raises(ValueError): + dgt.to_filename('path', filetype='mrtrix') diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 638f96e4..e21555b6 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -145,6 +145,10 @@ def generate_rasb(self): 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() @@ -161,11 +165,13 @@ def pole(self): def to_filename(self, filename, filetype='rasb'): """Write files (RASB, bvecs/bvals) to a given path.""" - if filetype == 'rasb': + 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 == 'fsl': + 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: @@ -222,6 +228,12 @@ def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, >>> norm_vals[-1] 3000 + >>> norm_vecs, norm_vals = normalize_gradients(bvecs, bvals, b_scale=False) + >>> norm_vals[0] + 0 + >>> np.all(norm_vals[1:] == 1000) + True + """ bvals = np.array(bvals, dtype='float32') bvecs = np.array(bvecs, dtype='float32') @@ -347,6 +359,10 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): >>> bvecs2ras(affine, [(0.0, 0.0, 0.0)]).tolist() [[0.0, 0.0, 0.0]] + >>> bvecs2ras(2.0 * np.eye(3), [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)], + ... norm=False).tolist() + [[2.0, 0.0, 0.0], [-2.0, 0.0, 0.0]] + """ if affine.shape == (4, 4): affine = affine[:3, :3] From d2bc103cbd197f5f34ef468496f0d7feb448b334 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 4 Nov 2019 14:27:19 -0800 Subject: [PATCH 12/12] fix: push code coverage to 100% of utils.vectors --- dmriprep/utils/tests/test_vectors.py | 11 ++++++++--- dmriprep/utils/vectors.py | 2 +- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index c2674cc4..46a57917 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -12,8 +12,10 @@ def test_corruption(tmpdir, dipy_test_data, monkeypatch): bvals = dipy_test_data['bvals'] bvecs = dipy_test_data['bvecs'] - # Test vector hemisphere coverage dgt = v.DiffusionGradientTable(**dipy_test_data) + affine = dgt.affine.copy() + + # Test vector hemisphere coverage assert np.all(dgt.pole == [0., 0., 0.]) dgt.to_filename('dwi.tsv') @@ -22,8 +24,11 @@ def test_corruption(tmpdir, dipy_test_data, monkeypatch): with pytest.raises(TypeError): dgt.to_filename('dwi', filetype='fsl') # You can do this iff the affine is set. - aff = namedtuple('Affine', ['affine'])(dgt.affine) # check accessing obj.affine - dgt = v.DiffusionGradientTable(dwi_file=aff) + # check accessing obj.affine + dgt = v.DiffusionGradientTable(dwi_file=namedtuple('Affine', ['affine'])(affine)) + assert np.all(dgt.affine == affine) + dgt = v.DiffusionGradientTable(dwi_file=affine) + assert np.all(dgt.affine == affine) # Perform various corruption checks using synthetic corrupted bval-bvec. dgt = v.DiffusionGradientTable() diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index e21555b6..d9284ef3 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -51,7 +51,7 @@ def __init__(self, dwi_file=None, bvecs=None, bvals=None, rasb_file=None, self.gradients = rasb_file if self.affine is not None: self.generate_vecval() - elif dwi_file and bvecs is not None and bvals is not None: + 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()