Merge pull request #1943 from frheault/increase_affine_consistency_pa…

…rt_3

Increase affine consistency in dipy.tracking.life and dipy.stats.analysis

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.

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,15 +178,15 @@ def test_afq_profile():
     bundle._data = bundle._data + affine[:3, 3]
     profile = afq_profile(data,
                           bundle,
-                          affine=affine,
+                          affine,
                           n_points=10,
                           weights=None)
 
     npt.assert_equal(profile, np.ones(10))
 
     # 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__':

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,8 +262,7 @@ def streamline_signal(self, streamline):
         return sig_out
 
 
-def voxel2streamline(streamline, transformed=False, affine=None,
-                     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
@@ -272,15 +272,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 +289,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 +312,7 @@ class FiberModel(ReconstModel):
         B.A. (2014). Validation and statistical inference in living
         connectomes. Nature Methods.
     """
+
     def __init__(self, gtab):
         """
         Parameters
@@ -343,8 +333,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,17 +352,15 @@ 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)
         vox_coords = unique_rows(np.round(all_coords).astype(np.intp))
         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 +405,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 +436,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 +445,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 +481,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):

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()

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``

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.
 """
@@ -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)