From 19634fa885f1322173266feea824547bff1991ce Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 31 Jul 2019 08:41:37 -0400 Subject: [PATCH 1/4] Force affine where required in afq and life --- dipy/stats/analysis.py | 13 +++---- dipy/stats/tests/test_analysis.py | 18 +++++----- dipy/tracking/life.py | 57 ++++++++++++------------------ dipy/tracking/tests/test_life.py | 10 +++--- doc/examples/afq_tract_profiles.py | 4 +-- 5 files changed, 42 insertions(+), 60 deletions(-) diff --git a/dipy/stats/analysis.py b/dipy/stats/analysis.py index 8b7faee411..c140b49ee8 100644 --- a/dipy/stats/analysis.py +++ b/dipy/stats/analysis.py @@ -345,7 +345,7 @@ def gaussian_weights(bundle, n_points=100, return_mahalnobis=False, return w / np.sum(w, 0) -def afq_profile(data, bundle, affine=None, n_points=100, +def afq_profile(data, bundle, affine, n_points=100, orient_by=None, weights=None, **weights_kwarg): """ Calculates a summarized profile of data for a bundle or tract @@ -362,23 +362,18 @@ def afq_profile(data, bundle, affine=None, n_points=100, The collection of streamlines (possibly already resampled into an array for each to have the same length) with which we are resampling. See Note below about orienting the streamlines. - - affine: 4-by-4 array, optional. - A transformation associated with the streamlines in the bundle. - Default: identity. - + affine : array_like (4, 4) + The mapping from voxel coordinates to streamline points. + The voxel_to_rasmm matrix, typically from a NIFTI file. n_points: int, optional The number of points to sample along the bundle. Default: 100. - orient_by: streamline, optional. A streamline to use as a standard to orient all of the streamlines in the bundle according to. - weights : 1D array or 2D array or callable (optional) Weight each streamline (1D) or each node (2D) when calculating the tract-profiles. Must sum to 1 across streamlines (in each node if relevant). If callable, this is a function that calculates weights. - weights_kwarg : key-word arguments Additional key-word arguments to pass to the weight-calculating function. Only to be used if weights is a callable. diff --git a/dipy/stats/tests/test_analysis.py b/dipy/stats/tests/test_analysis.py index 13e06de3ad..43702b9edd 100644 --- a/dipy/stats/tests/test_analysis.py +++ b/dipy/stats/tests/test_analysis.py @@ -141,34 +141,34 @@ def test_afq_profile(): [1, 0., 0], [2, 0, 0.]]])) - profile = afq_profile(data, bundle) + profile = afq_profile(data, bundle, np.eye(4)) npt.assert_equal(profile, np.ones(100)) - profile = afq_profile(data, bundle, affine=None, n_points=10, + profile = afq_profile(data, bundle, np.eye(4), n_points=10, weights=None) npt.assert_equal(profile, np.ones(10)) - profile = afq_profile(data, bundle, affine=None, + profile = afq_profile(data, bundle, np.eye(4), weights=gaussian_weights, stat=np.median) npt.assert_equal(profile, np.ones(100)) - profile = afq_profile(data, bundle, affine=None, orient_by=bundle[0], + profile = afq_profile(data, bundle, np.eye(4), orient_by=bundle[0], weights=gaussian_weights, stat=np.median) npt.assert_equal(profile, np.ones(100)) - profile = afq_profile(data, bundle, affine=None, n_points=10, + profile = afq_profile(data, bundle, np.eye(4), n_points=10, weights=None) npt.assert_equal(profile, np.ones(10)) - profile = afq_profile(data, bundle, affine=None, n_points=10, + profile = afq_profile(data, bundle, np.eye(4), n_points=10, weights=np.ones((2, 10)) * 0.5) npt.assert_equal(profile, np.ones(10)) # Disallow setting weights that don't sum to 1 across fibers/nodes: npt.assert_raises(ValueError, afq_profile, - data, bundle, affine=None, + data, bundle, np.eye(4), n_points=10, weights=np.ones((2, 10)) * 0.6) # Test using an affine: @@ -178,7 +178,7 @@ def test_afq_profile(): bundle._data = bundle._data + affine[:3, 3] profile = afq_profile(data, bundle, - affine=affine, + affine, n_points=10, weights=None) @@ -186,7 +186,7 @@ def test_afq_profile(): # Test for error-handling: empty_bundle = Streamlines([]) - npt.assert_raises(ValueError, afq_profile, data, empty_bundle) + npt.assert_raises(ValueError, afq_profile, data, empty_bundle, np.eye(4)) if __name__ == '__main__': diff --git a/dipy/tracking/life.py b/dipy/tracking/life.py index 8aa29a9dc5..a0d790575d 100644 --- a/dipy/tracking/life.py +++ b/dipy/tracking/life.py @@ -206,6 +206,7 @@ class LifeSignalMaker(object): A class for generating signals from streamlines in an efficient and speedy manner. """ + def __init__(self, gtab, evals=[0.001, 0, 0], sphere=None): """ Initialize a signal maker @@ -261,7 +262,7 @@ def streamline_signal(self, streamline): return sig_out -def voxel2streamline(streamline, transformed=False, affine=None, +def voxel2streamline(streamline, affine, transformed=False, unique_idx=None): """ Maps voxels to streamlines and streamlines to voxels, for setting up @@ -272,15 +273,9 @@ def voxel2streamline(streamline, transformed=False, affine=None, streamline : list A collection of streamlines, each n by 3, with n being the number of nodes in the fiber. - - affine : 4 by 4 array (optional) - Defines the spatial transformation from streamline to data. - Default: np.eye(4) - - transformed : bool (optional) - Whether the streamlines have been already transformed (in which case - they don't need to be transformed in here). - + affine : array_like (4, 4) + The mapping from voxel coordinates to streamline points. + The voxel_to_rasmm matrix, typically from a NIFTI file. unique_idx : array (optional). The unique indices in the streamlines @@ -295,12 +290,7 @@ def voxel2streamline(streamline, transformed=False, affine=None, this streamline passes through, which nodes of that streamline are in that voxel? """ - if transformed: - transformed_streamline = streamline - else: - if affine is None: - affine = np.eye(4) - transformed_streamline = transform_streamlines(streamline, affine) + transformed_streamline = transform_streamlines(streamline, affine) if unique_idx is None: all_coords = np.concatenate(transformed_streamline) @@ -323,6 +313,7 @@ class FiberModel(ReconstModel): B.A. (2014). Validation and statistical inference in living connectomes. Nature Methods. """ + def __init__(self, gtab): """ Parameters @@ -343,8 +334,9 @@ def setup(self, streamline, affine, evals=[0.001, 0, 0], sphere=None): ---------- streamline : list Streamlines, each is an array of shape (n, 3) - affine : 4 by 4 array - Mapping from the streamline coordinates to the data + affine : array_like (4, 4) + The mapping from voxel coordinates to streamline points. + The voxel_to_rasmm matrix, typically from a NIFTI file. evals : list (3 items, optional) The eigenvalues of the canonical tensor used as a response function. Default:[0.001, 0, 0]. @@ -361,8 +353,6 @@ def setup(self, streamline, affine, evals=[0.001, 0, 0], sphere=None): evals=evals, sphere=sphere) - if affine is None: - affine = np.eye(4) streamline = transform_streamlines(streamline, affine) # Assign some local variables, for shorthand: all_coords = np.concatenate(streamline) @@ -370,8 +360,8 @@ def setup(self, streamline, affine, evals=[0.001, 0, 0], sphere=None): del all_coords # We only consider the diffusion-weighted signals: n_bvecs = self.gtab.bvals[~self.gtab.b0s_mask].shape[0] - v2f, v2fn = voxel2streamline(streamline, transformed=True, - affine=affine, unique_idx=vox_coords) + v2f, v2fn = voxel2streamline(streamline, np.eye(4), + unique_idx=vox_coords) # How many fibers in each voxel (this will determine how many # components are in the matrix): n_unique_f = len(np.hstack(v2f.values())) @@ -416,7 +406,7 @@ def setup(self, streamline, affine, evals=[0.001, 0, 0], sphere=None): # Allocate the sparse matrix, using the more memory-efficient 'csr' # format: life_matrix = sps.csr_matrix((f_matrix_sig, - [f_matrix_row, f_matrix_col])) + [f_matrix_row, f_matrix_col])) return life_matrix, vox_coords @@ -447,7 +437,7 @@ def _signals(self, data, vox_coords): return (to_fit, weighted_signal, b0_signal, relative_signal, mean_sig, vox_data) - def fit(self, data, streamline, affine=None, evals=[0.001, 0, 0], + def fit(self, data, streamline, affine, evals=[0.001, 0, 0], sphere=None): """ Fit the LiFE FiberModel for data and a set of streamlines associated @@ -456,19 +446,15 @@ def fit(self, data, streamline, affine=None, evals=[0.001, 0, 0], Parameters ---------- data : 4D array - Diffusion-weighted data - + Diffusion-weighted data streamline : list - A bunch of streamlines - - affine: 4 by 4 array (optional) - The affine to go from the streamline coordinates to the data - coordinates. Defaults to use `np.eye(4)` - + A bunch of streamlines + affine : array_like (4, 4) + The mapping from voxel coordinates to streamline points. + The voxel_to_rasmm matrix, typically from a NIFTI file. evals : list (optional) - The eigenvalues of the tensor response function used in constructing - the model signal. Default: [0.001, 0, 0] - + The eigenvalues of the tensor response function used in constructing + the model signal. Default: [0.001, 0, 0] sphere: `dipy.core.Sphere` instance, or False Whether to approximate (and cache) the signal on a discrete sphere. This may confer a significant speed-up in setting up the @@ -496,6 +482,7 @@ class FiberFit(ReconstFit): """ A fit of the LiFE model to diffusion data """ + def __init__(self, fiber_model, life_matrix, vox_coords, to_fit, beta, weighted_signal, b0_signal, relative_signal, mean_sig, vox_data, streamline, affine, evals): diff --git a/dipy/tracking/tests/test_life.py b/dipy/tracking/tests/test_life.py index abfa8571c3..9f4305990d 100644 --- a/dipy/tracking/tests/test_life.py +++ b/dipy/tracking/tests/test_life.py @@ -77,7 +77,7 @@ def test_voxel2streamline(): streamline = [[[1.1, 2.4, 2.9], [4, 5, 3], [5, 6, 3], [6, 7, 3]], [[1, 2, 3], [4, 5, 3], [5, 6, 3]]] affine = np.eye(4) - v2f, v2fn = life.voxel2streamline(streamline, False, affine) + v2f, v2fn = life.voxel2streamline(streamline, affine) npt.assert_equal(v2f, {0: [0, 1], 1: [0, 1], 2: [0, 1], 3: [0]}) npt.assert_equal(v2fn, {0: {0: [0], 1: [1], 2: [2], 3: [3]}, 1: {0: [0], 1: [1], 2: [2]}}) @@ -87,7 +87,7 @@ def test_voxel2streamline(): [0, 0, 0, 1]]) xform_sl = life.transform_streamlines(streamline, np.linalg.inv(affine)) - v2f, v2fn = life.voxel2streamline(xform_sl, False, affine) + v2f, v2fn = life.voxel2streamline(xform_sl, affine) npt.assert_equal(v2f, {0: [0, 1], 1: [0, 1], 2: [0, 1], 3: [0]}) npt.assert_equal(v2fn, {0: {0: [0], 1: [1], 2: [2], 3: [3]}, 1: {0: [0], 1: [1], 2: [2]}}) @@ -127,7 +127,7 @@ def test_FiberFit(): streamline = [[[1, 2, 3], [4, 5, 3], [5, 6, 3], [6, 7, 3]], [[1, 2, 3], [4, 5, 3], [5, 6, 3]]] - fiber_matrix, vox_coords = FM.setup(streamline, None, evals) + fiber_matrix, vox_coords = FM.setup(streamline, np.eye(4), evals) w = np.array([0.5, 0.5]) sig = opt.spdot(fiber_matrix, w) + 1.0 # Add some isotropic stuff @@ -140,7 +140,7 @@ def test_FiberFit(): # Grab some realistic S0 values: this_data = np.concatenate([data[..., gtab.b0s_mask], this_data], -1) - fit = FM.fit(this_data, streamline) + fit = FM.fit(this_data, streamline, np.eye(4)) npt.assert_almost_equal(fit.predict()[1], fit.data[1], decimal=-1) @@ -165,7 +165,7 @@ def test_fit_data(): tensor_streamlines_vox = sft.streamlines life_model = life.FiberModel(gtab) - life_fit = life_model.fit(data, tensor_streamlines_vox) + life_fit = life_model.fit(data, tensor_streamlines_vox, np.eye(4)) model_error = life_fit.predict() - life_fit.data model_rmse = np.sqrt(np.mean(model_error ** 2, -1)) matlab_rmse, matlab_weights = dpd.matlab_life_results() diff --git a/doc/examples/afq_tract_profiles.py b/doc/examples/afq_tract_profiles.py index 32854ab6dd..e2324a5013 100644 --- a/doc/examples/afq_tract_profiles.py +++ b/doc/examples/afq_tract_profiles.py @@ -135,10 +135,10 @@ And then use the weights to calculate the tract profiles for each bundle """ -profile_cst_l = dsa.afq_profile(fa, oriented_cst_l, affine=img.affine, +profile_cst_l = dsa.afq_profile(fa, oriented_cst_l, img.affine, weights=w_cst_l) -profile_af_l = dsa.afq_profile(fa, oriented_af_l, affine=img.affine, +profile_af_l = dsa.afq_profile(fa, oriented_af_l, img.affine, weights=w_af_l) fig, (ax1, ax2) = plt.subplots(1, 2) From 20a692e79616e2d363d1e0bf7b5a1e855b748b32 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 31 Jul 2019 09:14:22 -0400 Subject: [PATCH 2/4] Flake8 and testing --- dipy/tracking/life.py | 4 ++-- doc/examples/afq_tract_profiles.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/dipy/tracking/life.py b/dipy/tracking/life.py index a0d790575d..1ce48dd202 100644 --- a/dipy/tracking/life.py +++ b/dipy/tracking/life.py @@ -453,8 +453,8 @@ def fit(self, data, streamline, affine, evals=[0.001, 0, 0], The mapping from voxel coordinates to streamline points. The voxel_to_rasmm matrix, typically from a NIFTI file. evals : list (optional) - The eigenvalues of the tensor response function used in constructing - the model signal. Default: [0.001, 0, 0] + The eigenvalues of the tensor response function used in + constructing the model signal. Default: [0.001, 0, 0] sphere: `dipy.core.Sphere` instance, or False Whether to approximate (and cache) the signal on a discrete sphere. This may confer a significant speed-up in setting up the diff --git a/doc/examples/afq_tract_profiles.py b/doc/examples/afq_tract_profiles.py index e2324a5013..d5380e713a 100644 --- a/doc/examples/afq_tract_profiles.py +++ b/doc/examples/afq_tract_profiles.py @@ -112,7 +112,7 @@ """ Read volumetric data from an image corresponding to this subject. -For the purpose of this, we've extracted only the FA within the bundles in +For the purpose of this, we've extracted only the FA within the bundles in question, but in real use, this is where you would add the FA map of your subject. """ From 6fb6f1fbc050188e6e2d7a0c31029908f620143c Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 31 Jul 2019 09:21:46 -0400 Subject: [PATCH 3/4] Updated the api_change file --- doc/api_changes.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/api_changes.rst b/doc/api_changes.rst index acc5a3d849..3dcb2c62d5 100644 --- a/doc/api_changes.rst +++ b/doc/api_changes.rst @@ -32,6 +32,14 @@ The API of ``dipy.io.streamlines.load_tractogram`` and When loading trk, tck, vtk, fib, or dpy) a reference nifti file is needed to guarantee proper spatial transformation handling. +**Spatial transformation handling** + +Functions from ``dipy.tracking.life`` were modified to enforce the +affine parameter and uniformize docstring. ``voxel2streamline``, +``setup`` and ``fit`` from class ``FiberModel`` were all modified. + +``afq_profile`` from ``dipy.stats.analysis`` was modified in a similar way. + **Simulation** - ``dipy.sims.voxel.SingleTensor`` has been replaced by ``dipy.sims.voxel.single_tensor`` From 4d03a2fa4d368576127fe92b8ba7e14388f5dbd4 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 31 Jul 2019 11:50:00 -0400 Subject: [PATCH 4/4] Remove forgotten parameters in life --- dipy/tracking/life.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dipy/tracking/life.py b/dipy/tracking/life.py index 1ce48dd202..bc00a5fedc 100644 --- a/dipy/tracking/life.py +++ b/dipy/tracking/life.py @@ -262,8 +262,7 @@ def streamline_signal(self, streamline): return sig_out -def voxel2streamline(streamline, affine, transformed=False, - unique_idx=None): +def voxel2streamline(streamline, affine, unique_idx=None): """ Maps voxels to streamlines and streamlines to voxels, for setting up the LiFE equations matrix