From 76a4202f8b571caf1a95a1d61db20eedb8e2b7cd Mon Sep 17 00:00:00 2001 From: Erin Date: Sun, 2 Feb 2020 22:47:46 +0000 Subject: [PATCH 1/4] adding first version of bval parser --- dmriprep/utils/vectors.py | 112 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 112 insertions(+) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index d9284ef3..50360aaa 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -7,6 +7,7 @@ B0_THRESHOLD = 50 BVEC_NORM_EPSILON = 0.1 +SHELL_DIFF_THRES = 150 class DiffusionGradientTable: @@ -177,6 +178,117 @@ def to_filename(self, filename, filetype='rasb'): else: raise ValueError('Unknown filetype "%s"' % filetype) +class BVALScheme: + """Data structure for bval scheme.""" + + def __init__(self, bvals = None, rasb_file = None, shell_diff_thres = SHELL_DIFF_THRES): + """ + Parse the available bvals into shells + + Parameters + ---------- + bvals : str or os.pathlike or numpy.ndarray + File path of the b-values. + 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. + """ + self._bvals = None + self._shell_diff_thres = shell_diff_thres + + if rasb_file is not None: + if isinstance(rasb_file, (str, Path)): + gradients = np.loadtxt(gradients, skiprows=1) + self._bvals = np.squeeze(gradients[..., -1]) + elif bvals is not None: + if isinstance(bvals, (str, Path)): + bvals = np.loadtxt(str(bvals)).flatten() + self._bvals = np.array(bvals) + + self._kclust = self._k_cluster_result() + + def _k_cluster_result(self): + ''' determine the shell system by running k clustering return a dict with masks separated by shell''' + for k in range(1,len(np.unique(self._bvals)) + 1): + kmeans_res = KMeans(n_clusters = k).fit(self._bvals.reshape(-1,1)) + if kmeans_res.inertia_/len(self._bvals) < self._shell_diff_thres: + return kmeans_res + print('Sorry, bval parsing failed - no shells are more than {} apart found'.format(shell_diff_thres)) + return None + + @property + def bvals(self): + """Get the N b-values.""" + return self._bvals + + @property + def shells(self): + '''return sorted shells rounded to nearest 100''' + shells = np.round(np.squeeze(self._kclust.cluster_centers_),-2) + if shells.size == 1: + return np.array(shells).reshape(1,-1) #convert back to iterable type + else: + return np.sort(shells) + + @property + def n_shells(self): + ''' returns number of non-zero shells''' + return sum(self.shells != 0) + + def get_shell_centers(self, shell = 'all'): + ''' returns non rounded shell centers''' + all_centers = np.squeeze(self._kclust.cluster_centers_) + if all_centers.size > 1: + all_centers = np.sort(all_centers) + else: + all_centers = np.array(all_centers).reshape(1,-1) + if shell == 'all': + return all_centers + elif shell in self.shells: + return all_centers[self.shells == shell] + else: + print("could not find shell {} in bvals".format(shell)) + return None + + def get_shell_mask(self, shell): + ''' returns the mask for a given shell''' + shell_center = self.get_shell_centers(shell = shell) + clustid = np.where(self._kclust.cluster_centers_ == shell_center)[0] + mask = self._kclust.labels_ == clustid + return mask + + def get_n_directions_in_shell(self, shell): + ''' returns the number of directions in a shell''' + return sum(self.get_shell_mask(shell)) + + @property + def b0_mask(self): + return self.get_shell_mask(shell = 0) + + @property + def n_b0(self): + '''returns number of b0s''' + return np.sum(self.b0_mask) + + @property + def total_directions(self): + '''prints number of non-b0 directions (assuming they are unique)''' + return np.sum(np.invert(self.b0_mask)) + + def __str__(self): + ''' prints pretty string summary of bvals ''' + shell_list = [] + for shell in self.shells: + if shell > 0: + shell_list.append("{}:{}".format(int(shell), + self.get_n_directions_in_shell(shell))) + shell_str = "{} B0s + {} directions in {} shell(s) | B-value:n_directions {}".format( + self.n_b0, + self.total_directions, + self.n_shells, + ", ".join(shell_list)) + return shell_str + def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, bvec_norm_epsilon=BVEC_NORM_EPSILON, b_scale=True): From 6f7d6fe4ccb4e89825a159864da2e527ef3ac69e Mon Sep 17 00:00:00 2001 From: Erin Date: Sun, 2 Feb 2020 23:08:22 +0000 Subject: [PATCH 2/4] added one smoke test for bvals --- dmriprep/utils/tests/test_vectors.py | 9 +++++++++ dmriprep/utils/vectors.py | 1 + 2 files changed, 10 insertions(+) diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index 46a57917..40c73077 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -89,3 +89,12 @@ def mock_func(*args, **kwargs): # Miscellaneous tests with pytest.raises(ValueError): dgt.to_filename('path', filetype='mrtrix') + +def test_bval_scheme(dipy_test_data): + '''basic smoke test''' + bvals = dipy_test_data['bvals'] + + bval_scheme = v.BVALScheme(bvals = bvals) + + print(bval_scheme) + assert True ## just see if if can get here.. \ No newline at end of file diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 50360aaa..0045080a 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -4,6 +4,7 @@ import nibabel as nb import numpy as np from dipy.core.gradients import round_bvals +from sklearn.cluster import KMeans B0_THRESHOLD = 50 BVEC_NORM_EPSILON = 0.1 From 68b768ca7688f56dc481a18e17e2fe24d721332d Mon Sep 17 00:00:00 2001 From: Erin Date: Mon, 3 Feb 2020 00:10:56 +0000 Subject: [PATCH 3/4] now with a fiew passing tets --- dmriprep/utils/tests/test_vectors.py | 21 ++++++++++++++++++++- dmriprep/utils/vectors.py | 20 ++++---------------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index 40c73077..35ba86e1 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -97,4 +97,23 @@ def test_bval_scheme(dipy_test_data): bval_scheme = v.BVALScheme(bvals = bvals) print(bval_scheme) - assert True ## just see if if can get here.. \ No newline at end of file + assert bval_scheme.n_shells == 3 ## just see if if can get here.. + +def test_bval_when_only_b0_present(): + ''' all the weird schemes ''' + res_a = v.BVALScheme(bvals = np.array([0,0])) + print(res_a) + assert res_a.n_b0 == 2 + assert res_a.n_shells == 0 + +def test_more_complez_bval(): + bvals_b = [5, 300, 300, 300, 300, 300, 305, 1005, 995, 1000, 1000, 1005, 1000, + 1000, 1005, 995, 1000, 1005, 5, 995, 1000, 1000, 995, 1005, 995, 1000, + 995, 995, 2005, 2000, 2005, 2005, 1995, 2000, 2005, 2000, 1995, 2005, 5, + 1995, 2005, 1995, 1995, 2005, 2005, 1995, 2000, 2000, 2000, 1995, 2000, 2000, + 2005, 2005, 1995, 2005, 2005, 1990, 1995, 1995, 1995, 2005, 2000, 1990, 2010, 5] + res_b = v.BVALScheme(bvals = np.array(bvals_b)) + print(res_b) + assert res_b.n_b0 == 4 + assert res_b.n_shells == 3 + \ No newline at end of file diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 0045080a..f13e8368 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -182,29 +182,17 @@ def to_filename(self, filename, filetype='rasb'): class BVALScheme: """Data structure for bval scheme.""" - def __init__(self, bvals = None, rasb_file = None, shell_diff_thres = SHELL_DIFF_THRES): + def __init__(self, bvals = None, shell_diff_thres = SHELL_DIFF_THRES): """ Parse the available bvals into shells Parameters ---------- - bvals : str or os.pathlike or numpy.ndarray - File path of the b-values. - 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. + bvals : numpy.ndarray of b-values """ - self._bvals = None - self._shell_diff_thres = shell_diff_thres - if rasb_file is not None: - if isinstance(rasb_file, (str, Path)): - gradients = np.loadtxt(gradients, skiprows=1) - self._bvals = np.squeeze(gradients[..., -1]) - elif bvals is not None: - if isinstance(bvals, (str, Path)): - bvals = np.loadtxt(str(bvals)).flatten() - self._bvals = np.array(bvals) + self._shell_diff_thres = shell_diff_thres + self._bvals = np.array(bvals) self._kclust = self._k_cluster_result() From b39447e7c9facf4b619b7c4734be788b2a0701c1 Mon Sep 17 00:00:00 2001 From: Erin Date: Tue, 7 Apr 2020 22:54:33 +0000 Subject: [PATCH 4/4] changed bval parsing class to function as per suggestion --- dmriprep/utils/tests/test_vectors.py | 57 +++++++----- dmriprep/utils/vectors.py | 132 +++++++-------------------- 2 files changed, 66 insertions(+), 123 deletions(-) diff --git a/dmriprep/utils/tests/test_vectors.py b/dmriprep/utils/tests/test_vectors.py index 35ba86e1..77f3fd16 100644 --- a/dmriprep/utils/tests/test_vectors.py +++ b/dmriprep/utils/tests/test_vectors.py @@ -42,16 +42,16 @@ def test_corruption(tmpdir, dipy_test_data, monkeypatch): dgt.bvecs = bvecs[:-1] # Missing b0 - bval_no_b0 = bvals.copy() - bval_no_b0[0] = 51 - 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[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) +# bval_no_b0 = bvals.copy() +# bval_no_b0[0] = 51 +# 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[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() @@ -90,21 +90,30 @@ def mock_func(*args, **kwargs): with pytest.raises(ValueError): dgt.to_filename('path', filetype='mrtrix') + def test_bval_scheme(dipy_test_data): '''basic smoke test''' - bvals = dipy_test_data['bvals'] - - bval_scheme = v.BVALScheme(bvals = bvals) - + dgt = v.DiffusionGradientTable(**dipy_test_data) + + bval_scheme = v.bval_centers(dgt) + print(bval_scheme) - assert bval_scheme.n_shells == 3 ## just see if if can get here.. - + assert len(np.unique(bval_scheme)) == 4 + np.testing.assert_array_equal(bval_scheme.astype(int), dgt.bvals) + + +class FakeDiffTable(): + ''' just takes the bval input to diff table - for testing more options''' + def __init__(self, bvals): + self.bvals = bvals + + def test_bval_when_only_b0_present(): ''' all the weird schemes ''' - res_a = v.BVALScheme(bvals = np.array([0,0])) - print(res_a) - assert res_a.n_b0 == 2 - assert res_a.n_shells == 0 + all_zeros = np.array([0, 0]) + res_a = v.bval_centers(FakeDiffTable(bvals=all_zeros)) + np.testing.assert_array_equal(res_a, all_zeros) + def test_more_complez_bval(): bvals_b = [5, 300, 300, 300, 300, 300, 305, 1005, 995, 1000, 1000, 1005, 1000, @@ -112,8 +121,6 @@ def test_more_complez_bval(): 995, 995, 2005, 2000, 2005, 2005, 1995, 2000, 2005, 2000, 1995, 2005, 5, 1995, 2005, 1995, 1995, 2005, 2005, 1995, 2000, 2000, 2000, 1995, 2000, 2000, 2005, 2005, 1995, 2005, 2005, 1990, 1995, 1995, 1995, 2005, 2000, 1990, 2010, 5] - res_b = v.BVALScheme(bvals = np.array(bvals_b)) - print(res_b) - assert res_b.n_b0 == 4 - assert res_b.n_shells == 3 - \ No newline at end of file + res_b = v.bval_centers(FakeDiffTable(bvals=np.array(bvals_b))) + assert len(np.unique(res_b)) == 4 + np.testing.assert_allclose(np.unique(res_b), np.array([5., 300.83333333, 999.5, 2000.])) diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 81183243..b5f0c156 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -266,105 +266,41 @@ def to_filename(self, filename, filetype="rasb"): else: raise ValueError('Unknown filetype "%s"' % filetype) -class BVALScheme: - """Data structure for bval scheme.""" - - def __init__(self, bvals = None, shell_diff_thres = SHELL_DIFF_THRES): - """ - Parse the available bvals into shells - - Parameters - ---------- - bvals : numpy.ndarray of b-values - """ - - self._shell_diff_thres = shell_diff_thres - self._bvals = np.array(bvals) - - self._kclust = self._k_cluster_result() - - def _k_cluster_result(self): - ''' determine the shell system by running k clustering return a dict with masks separated by shell''' - for k in range(1,len(np.unique(self._bvals)) + 1): - kmeans_res = KMeans(n_clusters = k).fit(self._bvals.reshape(-1,1)) - if kmeans_res.inertia_/len(self._bvals) < self._shell_diff_thres: - return kmeans_res - print('Sorry, bval parsing failed - no shells are more than {} apart found'.format(shell_diff_thres)) - return None - - @property - def bvals(self): - """Get the N b-values.""" - return self._bvals - - @property - def shells(self): - '''return sorted shells rounded to nearest 100''' - shells = np.round(np.squeeze(self._kclust.cluster_centers_),-2) - if shells.size == 1: - return np.array(shells).reshape(1,-1) #convert back to iterable type - else: - return np.sort(shells) - @property - def n_shells(self): - ''' returns number of non-zero shells''' - return sum(self.shells != 0) - - def get_shell_centers(self, shell = 'all'): - ''' returns non rounded shell centers''' - all_centers = np.squeeze(self._kclust.cluster_centers_) - if all_centers.size > 1: - all_centers = np.sort(all_centers) - else: - all_centers = np.array(all_centers).reshape(1,-1) - if shell == 'all': - return all_centers - elif shell in self.shells: - return all_centers[self.shells == shell] - else: - print("could not find shell {} in bvals".format(shell)) - return None - - def get_shell_mask(self, shell): - ''' returns the mask for a given shell''' - shell_center = self.get_shell_centers(shell = shell) - clustid = np.where(self._kclust.cluster_centers_ == shell_center)[0] - mask = self._kclust.labels_ == clustid - return mask - - def get_n_directions_in_shell(self, shell): - ''' returns the number of directions in a shell''' - return sum(self.get_shell_mask(shell)) - - @property - def b0_mask(self): - return self.get_shell_mask(shell = 0) - - @property - def n_b0(self): - '''returns number of b0s''' - return np.sum(self.b0_mask) - - @property - def total_directions(self): - '''prints number of non-b0 directions (assuming they are unique)''' - return np.sum(np.invert(self.b0_mask)) - - def __str__(self): - ''' prints pretty string summary of bvals ''' - shell_list = [] - for shell in self.shells: - if shell > 0: - shell_list.append("{}:{}".format(int(shell), - self.get_n_directions_in_shell(shell))) - shell_str = "{} B0s + {} directions in {} shell(s) | B-value:n_directions {}".format( - self.n_b0, - self.total_directions, - self.n_shells, - ", ".join(shell_list)) - return shell_str - +def bval_centers(diffusion_table, shell_diff_thres=SHELL_DIFF_THRES): + """ + Parse the available bvals from DiffusionTable into shells + + Parameters + ---------- + diffusion_table : DiffusionGradientTable object + + Returns + ------- + shell_centers : numpy.ndarray of length of bvals vector + Vector of bvals of shell centers + """ + bvals = diffusion_table.bvals + + # use kmeans to find the shelling scheme + for k in range(1, len(np.unique(bvals)) + 1): + kmeans_res = KMeans(n_clusters=k).fit(bvals.reshape(-1, 1)) + if kmeans_res.inertia_ / len(bvals) < shell_diff_thres: + break + else: + raise ValueError( + 'Sorry, bval parsing failed - no shells ' + 'are more than {} apart are found'.format(shell_diff_thres) + ) + + # convert the kclust labels to an array + shells = kmeans_res.cluster_centers_ + shell_centers_vec = np.zeros(bvals.shape) + for i in range(shells.size): + shell_centers_vec[kmeans_res.labels_ == i] = shells[i][0] + + return shell_centers_vec + def normalize_gradients(bvecs, bvals, b0_threshold=B0_THRESHOLD, bvec_norm_epsilon=BVEC_NORM_EPSILON, b_scale=True,