From 5f51f9d602a13f293330d54adcf42f929e86ce0d Mon Sep 17 00:00:00 2001 From: dPys Date: Mon, 13 Jan 2020 12:43:25 -0800 Subject: [PATCH 1/6] [ENH] Add API to vectors for reorienting bvecs from rasb .tsv and affine list --- .idea/dmriprep-1.iml | 13 +++ .idea/libraries/R_User_Library.xml | 6 ++ .idea/misc.xml | 7 ++ .idea/modules.xml | 8 ++ .idea/vcs.xml | 6 ++ .idea/workspace.xml | 160 +++++++++++++++++++++++++++++ dmriprep/utils/vectors.py | 29 ++++++ 7 files changed, 229 insertions(+) create mode 100644 .idea/dmriprep-1.iml create mode 100644 .idea/libraries/R_User_Library.xml create mode 100644 .idea/misc.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/vcs.xml create mode 100644 .idea/workspace.xml diff --git a/.idea/dmriprep-1.iml b/.idea/dmriprep-1.iml new file mode 100644 index 00000000..86abb4ca --- /dev/null +++ b/.idea/dmriprep-1.iml @@ -0,0 +1,13 @@ + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/libraries/R_User_Library.xml b/.idea/libraries/R_User_Library.xml new file mode 100644 index 00000000..71f5ff74 --- /dev/null +++ b/.idea/libraries/R_User_Library.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 00000000..1d8a5d64 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 00000000..f0a271af --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 00000000..94a25f7f --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/.idea/workspace.xml b/.idea/workspace.xml new file mode 100644 index 00000000..4f8793bd --- /dev/null +++ b/.idea/workspace.xml @@ -0,0 +1,160 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - 1578948143379 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file From f5ca8c1c5b85fa4212a5ba8e0ee542f84e5f94e5 Mon Sep 17 00:00:00 2001 From: dPys Date: Tue, 11 Feb 2020 16:13:01 -0600 Subject: [PATCH 3/6] Add signal prediction vector handling functions --- dmriprep/utils/vectors.py | 66 +++++++++++++++++++++++++++------------ 1 file changed, 46 insertions(+), 20 deletions(-) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index b63afad6..4063015e 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -377,30 +377,56 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): return rotated_bvecs -def reorient_bvecs_from_ras_b(ras_b, affines): +def _nonoverlapping_qspace_samples( + prediction_bval, prediction_bvec, all_bvals, all_bvecs, cutoff +): + """Ensure that none of the training samples are too close to the sample to predict. + Parameters + """ + min_bval = min(min(all_bvals), prediction_bval) + all_qvals = np.sqrt(all_bvals - min_bval) + prediction_qval = np.sqrt(prediction_bval - min_bval) + + # Convert q values to percent of maximum qval + max_qval = max(max(all_qvals), prediction_qval) + all_qvals_scaled = all_qvals / max_qval * 100 + scaled_qvecs = all_bvecs * all_qvals_scaled[:, np.newaxis] + scaled_prediction_qvec = prediction_bvec * (prediction_qval / max_qval * 100) + + # Calculate the distance between the sampled qvecs and the prediction qvec + ok_samples = ( + np.linalg.norm(scaled_qvecs - scaled_prediction_qvec, axis=1) > cutoff + ) * (np.linalg.norm(scaled_qvecs + scaled_prediction_qvec, axis=1) > cutoff) + + return ok_samples + + +def _rasb_to_bvec_list(in_rasb): """ - Reorient the vectors from a rasb .tsv file. - When correcting for motion, rotation of the diffusion-weighted volumes - might cause systematic bias in rotationally invariant measures, such as FA - and MD, and also cause characteristic biases in tractography, unless the - gradient directions are appropriately reoriented to compensate for this - effect [Leemans2009]_. + Create a list of b-vectors from a rasb gradient table. Parameters ---------- - 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. - - affines : list or ndarray of shape (n, 4, 4) or (n, 3, 3) - Each entry in this list or array contain either an affine - transformation (4,4) or a rotation matrix (3, 3). - In both cases, the transformations encode the rotation that was applied - to the image corresponding to one of the non-zero gradient directions. + in_rasb : str or os.pathlike + File path to a RAS-B gradient table. """ - from dipy.core.gradients import gradient_table_from_bvals_bvecs, reorient_bvecs + import numpy as np - ras_b_mat = np.genfromtxt(ras_b, delimiter='\t') - gt = gradient_table_from_bvals_bvecs(ras_b_mat[:,3], ras_b_mat[:,0:3], b0_threshold=50) + ras_b_mat = np.genfromtxt(in_rasb, delimiter="\t") + bvec = [vec for vec in ras_b_mat[:, 0:3] if not np.isclose(all(vec), 0)] + return list(bvec) + + +def _rasb_to_bval_floats(in_rasb): + """ + Create a list of b-values from a rasb gradient table. + + Parameters + ---------- + in_rasb : str or os.pathlike + File path to a RAS-B gradient table. + """ + import numpy as np - return reorient_bvecs(gt, affines) \ No newline at end of file + ras_b_mat = np.genfromtxt(in_rasb, delimiter="\t") + return [float(bval) for bval in ras_b_mat[:, 3] if bval > 0] \ No newline at end of file From 66f8b4857322e3ab684a022daf394ecaee04e49a Mon Sep 17 00:00:00 2001 From: dPys Date: Tue, 24 Mar 2020 13:01:44 -0500 Subject: [PATCH 4/6] Add doctest to non_overlapping_qspace_samples, made general-purpose --- dmriprep/utils/vectors.py | 71 +++++++++++++++++++++++++++++++++------ 1 file changed, 61 insertions(+), 10 deletions(-) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 4063015e..19ff10ff 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -377,26 +377,70 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): return rotated_bvecs -def _nonoverlapping_qspace_samples( - prediction_bval, prediction_bvec, all_bvals, all_bvecs, cutoff -): +def nonoverlapping_qspace_samples(sample_bval, sample_bvec, all_bvals, + all_bvecs, cutoff=2): """Ensure that none of the training samples are too close to the sample to predict. Parameters + + Parameters + ---------- + sample_bval : int + A single b-value sampled along the sphere. + sample_bvec : int + A single b-vector sampled along the sphere. + Should correspond to `sample_bval`. + all_bvals : ndarray + A 1D vector of all b-values from the diffusion series. + all_bvecs: ndarray + A 3 x n vector of all vectors from the diffusion series, + where n is the total number of samples. + cutoff : float + A minimal allowable q-space distance between points on + the sphere. + + Returns + ------- + ok_samples : boolean ndarray + True for q-vectors whose spatial distribution along + the sphere is non-overlapping, else False. + + Examples + -------- + >>> bvec1 = np.array([1, 0, 0]) + >>> bvec2 = np.array([1, 0, 0]) + >>> bvec3 = np.array([0, 1, 0]) + >>> bval1 = 1000 + >>> bval2 = 1000 + >>> bval3 = 1000 + >>> all_bvals = np.array([0, bval2, bval3]) + >>> all_bvecs = np.array([np.zeros(3), bvec2, bvec3]) + >>> # Case 1: overlapping + >>> nonoverlapping_qspace_samples(bval1, bvec1, all_bvals, all_bvecs, cutoff=2) + array([ True, False, True]) + >>> all_bvals = np.array([0, bval1, bval2]) + >>> all_bvecs = np.array([np.zeros(3), bvec1, bvec2]) + >>> # Case 2: non-overlapping + >>> nonoverlapping_qspace_samples(bval3, bvec3, all_bvals, all_bvecs, cutoff=2) + array([ True, True, True]) """ - min_bval = min(min(all_bvals), prediction_bval) + min_bval = min(min(all_bvals), sample_bval) + max_bval = max(max(all_bvals), sample_bval) + if min_bval == max_bval: + raise ValueError('All b-values are identical') + all_qvals = np.sqrt(all_bvals - min_bval) - prediction_qval = np.sqrt(prediction_bval - min_bval) + sample_qval = np.sqrt(sample_bval - min_bval) # Convert q values to percent of maximum qval - max_qval = max(max(all_qvals), prediction_qval) + max_qval = max(max(all_qvals), sample_qval) all_qvals_scaled = all_qvals / max_qval * 100 scaled_qvecs = all_bvecs * all_qvals_scaled[:, np.newaxis] - scaled_prediction_qvec = prediction_bvec * (prediction_qval / max_qval * 100) + scaled_sample_qvec = sample_bvec * (sample_qval / max_qval * 100) - # Calculate the distance between the sampled qvecs and the prediction qvec + # Calculate the distance between all qvecs and the sample qvec ok_samples = ( - np.linalg.norm(scaled_qvecs - scaled_prediction_qvec, axis=1) > cutoff - ) * (np.linalg.norm(scaled_qvecs + scaled_prediction_qvec, axis=1) > cutoff) + np.linalg.norm(scaled_qvecs - scaled_sample_qvec, axis=1) > cutoff + ) * (np.linalg.norm(scaled_qvecs + scaled_sample_qvec, axis=1) > cutoff) return ok_samples @@ -409,6 +453,9 @@ def _rasb_to_bvec_list(in_rasb): ---------- in_rasb : str or os.pathlike File path to a RAS-B gradient table. + Returns + ------- + List of b-vectors as floats. """ import numpy as np @@ -425,6 +472,10 @@ def _rasb_to_bval_floats(in_rasb): ---------- in_rasb : str or os.pathlike File path to a RAS-B gradient table. + + Returns + ------- + List of b-values as floats. """ import numpy as np From 7e1b11e829d6bdc5382c9d166f45fb113c153bd4 Mon Sep 17 00:00:00 2001 From: Derek Pisner <16432683+dPys@users.noreply.github.com> Date: Tue, 24 Mar 2020 13:09:51 -0500 Subject: [PATCH 5/6] Update dmriprep/utils/vectors.py Co-Authored-By: Ariel Rokem --- dmriprep/utils/vectors.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 19ff10ff..84e6fa47 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -379,8 +379,9 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): def nonoverlapping_qspace_samples(sample_bval, sample_bvec, all_bvals, all_bvecs, cutoff=2): - """Ensure that none of the training samples are too close to the sample to predict. - Parameters + """ + Checks the q-space overlap (within some distance) between a sample + and a collection of q-space points. Parameters ---------- @@ -453,6 +454,7 @@ def _rasb_to_bvec_list(in_rasb): ---------- in_rasb : str or os.pathlike File path to a RAS-B gradient table. + Returns ------- List of b-vectors as floats. @@ -480,4 +482,4 @@ def _rasb_to_bval_floats(in_rasb): import numpy as np ras_b_mat = np.genfromtxt(in_rasb, delimiter="\t") - return [float(bval) for bval in ras_b_mat[:, 3] if bval > 0] \ No newline at end of file + return [float(bval) for bval in ras_b_mat[:, 3] if bval > 0] From 8b365f8f67c5020fa442d6b777d39e75e2bb3697 Mon Sep 17 00:00:00 2001 From: dPys Date: Tue, 24 Mar 2020 15:09:26 -0500 Subject: [PATCH 6/6] Change exception error to user warning in normalize gradients --- dmriprep/utils/registration.py | 169 +++++++++++++++++++++++++++++++++ dmriprep/utils/vectors.py | 2 +- 2 files changed, 170 insertions(+), 1 deletion(-) create mode 100644 dmriprep/utils/registration.py diff --git a/dmriprep/utils/registration.py b/dmriprep/utils/registration.py new file mode 100644 index 00000000..5d1dcd0f --- /dev/null +++ b/dmriprep/utils/registration.py @@ -0,0 +1,169 @@ +""" +Linear affine registration tools for motion correction. +""" +import numpy as np +import nibabel as nb +from dipy.align.metrics import CCMetric, EMMetric, SSDMetric +from dipy.align.imaffine import ( + transform_centers_of_mass, + AffineMap, + MutualInformationMetric, + AffineRegistration, +) +from dipy.align.transforms import ( + TranslationTransform3D, + RigidTransform3D, + AffineTransform3D, +) +from nipype.utils.filemanip import fname_presuffix + +syn_metric_dict = {"CC": CCMetric, "EM": EMMetric, "SSD": SSDMetric} + +__all__ = [ + "c_of_mass", + "translation", + "rigid", + "affine", + "affine_registration", +] + + +def apply_affine(moving, static, transform_affine, invert=False): + """Apply an affine to transform an image from one space to another. + + Parameters + ---------- + moving : array + The image to be resampled + + static : array + + Returns + ------- + warped_img : the moving array warped into the static array's space. + + """ + affine_map = AffineMap( + transform_affine, static.shape, static.affine, moving.shape, moving.affine + ) + if invert is True: + warped_arr = affine_map.transform_inverse(np.asarray(moving.dataobj)) + else: + warped_arr = affine_map.transform(np.asarray(moving.dataobj)) + + return nb.Nifti1Image(warped_arr, static.affine) + + +def average_affines(transforms): + affine_list = [np.load(aff) for aff in transforms] + average_affine_file = fname_presuffix( + transforms[0], use_ext=False, suffix="_average.npy" + ) + np.save(average_affine_file, np.mean(affine_list, axis=0)) + return average_affine_file + + +# Affine registration pipeline: +affine_metric_dict = {"MI": MutualInformationMetric, "CC": CCMetric} + + +def c_of_mass( + moving, static, static_affine, moving_affine, reg, starting_affine, params0=None +): + transform = transform_centers_of_mass(static, static_affine, moving, moving_affine) + transformed = transform.transform(moving) + return transformed, transform.affine + + +def translation( + moving, static, static_affine, moving_affine, reg, starting_affine, params0=None +): + transform = TranslationTransform3D() + translation = reg.optimize( + static, + moving, + transform, + params0, + static_affine, + moving_affine, + starting_affine=starting_affine, + ) + + return translation.transform(moving), translation.affine + + +def rigid( + moving, static, static_affine, moving_affine, reg, starting_affine, params0=None +): + transform = RigidTransform3D() + rigid = reg.optimize( + static, + moving, + transform, + params0, + static_affine, + moving_affine, + starting_affine=starting_affine, + ) + return rigid.transform(moving), rigid.affine + + +def affine( + moving, static, static_affine, moving_affine, reg, starting_affine, params0=None +): + transform = AffineTransform3D() + affine = reg.optimize( + static, + moving, + transform, + params0, + static_affine, + moving_affine, + starting_affine=starting_affine, + ) + + return affine.transform(moving), affine.affine + + +def affine_registration( + moving, + static, + nbins, + sampling_prop, + metric, + pipeline, + level_iters, + sigmas, + factors, + params0, +): + """ + Find the affine transformation between two 3D images. + + Parameters + ---------- + + """ + # Define the Affine registration object we'll use with the chosen metric: + use_metric = affine_metric_dict[metric](nbins, sampling_prop) + affreg = AffineRegistration( + metric=use_metric, level_iters=level_iters, sigmas=sigmas, factors=factors + ) + + if not params0: + starting_affine = np.eye(4) + else: + starting_affine = params0 + + # Go through the selected transformation: + for func in pipeline: + transformed, starting_affine = func( + np.asarray(moving.dataobj), + np.asarray(static.dataobj), + static.affine, + moving.affine, + affreg, + starting_affine, + params0, + ) + return nb.Nifti1Image(np.array(transformed), static.affine), starting_affine diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 84e6fa47..e409d187 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -243,7 +243,7 @@ def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, # Check for bval-bvec discrepancy. if not np.all(b0s == b0_vecs): - raise ValueError( + raise UserWarning( 'Inconsistent bvals and bvecs (%d, %d low-b, respectively).' % (b0s.sum(), b0_vecs.sum()))