Merge branch 'malu-ms'

dipy/tracking/tests/test_utils.py

@@ -0,0 +1,83 @@
+import numpy as np
+from dipy.io.bvectxt import orientation_from_string
+from dipy.tracking.utils import reorder_voxels_affine, move_streamlines
+from numpy.testing import assert_array_almost_equal, assert_array_equal
+
+def make_streamlines():
+    streamlines = [ np.array([[0, 0, 0],
+                              [1, 1, 1],
+                              [2, 2, 2],
+                              [5, 10, 12]], 'float'),
+                    np.array([[1, 2, 3],
+                              [3, 2, 0],
+                              [5, 20, 33],
+                              [40, 80, 120]], 'float') ]
+    return streamlines
+
+def test_move_streamlines():
+    streamlines = make_streamlines()
+    affine = np.eye(4)
+    new_streamlines = move_streamlines(streamlines, affine)
+    for i, test_sl in enumerate(new_streamlines):
+        assert_array_equal(test_sl, streamlines[i])
+
+    affine[:3,3] += (4,5,6)
+    new_streamlines = move_streamlines(streamlines, affine)
+    for i, test_sl in enumerate(new_streamlines):
+        assert_array_equal(test_sl, streamlines[i]+(4, 5, 6))
+
+    affine = np.eye(4)
+    affine = affine[[2,1,0,3]]
+    new_streamlines = move_streamlines(streamlines, affine)
+    for i, test_sl in enumerate(new_streamlines):
+        assert_array_equal(test_sl, streamlines[i][:, [2, 1, 0]])
+
+def test_voxel_ornt():
+    sh = (40,40,40)
+    sz = (1, 2, 3)
+    I4 = np.eye(4)
+
+    ras = orientation_from_string('ras')
+    sra = orientation_from_string('sra')
+    lpi = orientation_from_string('lpi')
+    srp = orientation_from_string('srp')
+
+    affine = reorder_voxels_affine(ras, ras, sh, sz)
+    assert_array_equal(affine, I4)
+    affine = reorder_voxels_affine(sra, sra, sh, sz)
+    assert_array_equal(affine, I4)
+    affine = reorder_voxels_affine(lpi, lpi, sh, sz)
+    assert_array_equal(affine, I4)
+    affine = reorder_voxels_affine(srp, srp, sh, sz)
+    assert_array_equal(affine, I4)
+
+    streamlines = make_streamlines()
+    box = np.array(sh)*sz
+
+    sra_affine = reorder_voxels_affine(ras, sra, sh, sz)
+    toras_affine = reorder_voxels_affine(sra, ras, sh, sz)
+    assert_array_equal(np.dot(toras_affine, sra_affine), I4)
+    expected_sl = (sl[:, [2, 0, 1]] for sl in streamlines)
+    test_sl = move_streamlines(streamlines, sra_affine)
+    for ii in xrange(len(streamlines)):
+        assert_array_equal(test_sl.next(), expected_sl.next())
+
+    lpi_affine = reorder_voxels_affine(ras, lpi, sh, sz)
+    toras_affine = reorder_voxels_affine(lpi, ras, sh, sz)
+    assert_array_equal(np.dot(toras_affine, lpi_affine), I4)
+    expected_sl = (box - sl for sl in streamlines)
+    test_sl = move_streamlines(streamlines, lpi_affine)
+    for ii in xrange(len(streamlines)):
+        assert_array_equal(test_sl.next(), expected_sl.next())
+
+    srp_affine = reorder_voxels_affine(ras, srp, sh, sz)
+    toras_affine = reorder_voxels_affine(srp, ras, (40,40,40), (3,1,2))
+    assert_array_equal(np.dot(toras_affine, srp_affine), I4)
+    expected_sl = [sl.copy() for sl in streamlines]
+    for sl in expected_sl:
+        sl[:, 1] = box[1] - sl[:, 1]
+    expected_sl = (sl[:, [2, 0, 1]] for sl in expected_sl)
+    test_sl = move_streamlines(streamlines, srp_affine)
+    for ii in xrange(len(streamlines)):
+        assert_array_equal(test_sl.next(), expected_sl.next())
+

dipy/tracking/utils.py

@@ -1,7 +1,9 @@
 from __future__ import division
 from numpy import asarray, array, atleast_3d, ceil, concatenate, empty, \
-        mgrid, sqrt, zeros
+        eye, mgrid, sqrt, zeros, linalg, diag, dot
 from collections import defaultdict
+from dipy.io.bvectxt import ornt_mapping
+
 
 def streamline_counts(streamlines, vol_dims, voxel_size):
     """Counts the number of unique streamlines that pass though each voxel
@@ -81,7 +83,7 @@ def subsegment(streamlines, max_segment_length):
 
     Returns:
     --------
-    new_streamlines : generator
+    output_streamlines : generator
         A set of streamlines.
 
     Notes:
@@ -114,29 +116,29 @@ def subsegment(streamlines, max_segment_length):
         length = sqrt((diff*diff).sum(-1))
         num_segments = ceil(length/max_segment_length).astype('int')
 
-        new_sl = empty((num_segments.sum()+1, 3), 'float')
-        new_sl[0] = sl[0]
+        output_sl = empty((num_segments.sum()+1, 3), 'float')
+        output_sl[0] = sl[0]
 
         count = 1
         for ii in xrange(len(num_segments)):
             ns = num_segments[ii]
             if ns == 1:
-                new_sl[count] = sl[ii+1]
+                output_sl[count] = sl[ii+1]
                 count += 1
             elif ns > 1:
                 small_d = diff[ii]/ns
                 point = sl[ii]
                 for jj in xrange(ns):
                     point = point + small_d
-                    new_sl[count] = point
+                    output_sl[count] = point
                     count += 1
             elif ns == 0:
                 pass
                 #repeated point
             else:
                 #this should never happen because ns should be a posative int
                 assert(ns >= 0)
-        yield new_sl
+        yield output_sl
 
 def seeds_from_mask(mask, density, voxel_size=(1,1,1)):
     """Takes a binary mask and returns seeds in voxels != 0
@@ -151,7 +153,8 @@ def seeds_from_mask(mask, density, voxel_size=(1,1,1)):
     >>> mask = zeros((3,3,3), 'bool')
     >>> mask[0,0,0] = 1
     >>> seeds_from_mask(mask, [1,1,1], [1,1,1])
-    array([[ 0.5,  0.5,  0.5]])
+    array([[ 0.5,  0.5,  0.5]])counts_mask = counts > 0
+
     >>> seeds_from_mask(mask, [1,2,3], [1,1,1])
     array([[ 0.5       ,  0.25      ,  0.16666667],
            [ 0.5       ,  0.25      ,  0.5       ],
@@ -211,7 +214,8 @@ def target(streamlines, target_mask, voxel_size):
     IndexError
         When the points of the streamlines lie outside of the target_mask
 
-    See Also:
+    See Also:from numpy import linalg as LA
+
     ---------
     streamline_counts
 
@@ -235,7 +239,8 @@ def target(streamlines, target_mask, voxel_size):
 def merge_streamlines(backward, forward):
     """Merges two sets of streamlines seeded at the same points
 
-    Because the first point of each streamline pair should be the same, only
+   
+ Because the first point of each streamline pair should be the same, only
     one is kept
 
     Parameters:
@@ -272,3 +277,81 @@ def merge_streamlines(backward, forward):
     F = iter(forward)
     while True:
         yield concatenate((B.next()[:0:-1], F.next()))
+
+
+
+def move_streamlines(streamlines, affine):
+    """Applies a linear transformation, given by affine, to streamlines
+
+    Parameters:
+    -----------
+    streamlines : sequence
+        A set of streamlines to be transformed.
+    affine : array (4, 4)
+        A linear tranformation to be applied to the streamlines. The last row
+        of affine should be [0, 0, 0, 1].
+
+    Returns:
+    --------
+    streamlines : generator
+        A sequence of transformed streamlines
+    """
+    for sl in streamlines:
+        yield dot(sl, affine[:3,:3].T) + affine[:3,3]
+
+def reorder_voxels_affine(input_ornt, output_ornt, shape, voxel_size):
+    """Calculates a linear tranformation equivelent to chaning voxel order
+
+    Calculates a linear tranformation A such that [a, b, c, 1] = A[x, y, z, 1].
+    where [x, y, z] is a point in the coordinate system defined by input_ornt
+    and [a, b, c] is the same point in the coordinate system defined by
+    output_ornt.
+
+    Parameters:
+    -----------
+    input_ornt : array (n, 2)
+        A description of the orientation of a point in n-space. See
+        nibabel.orientation or dipy.io.bvectxt for more information.
+    output_ornt : array (n, 2)
+        A description of the orientation of a point in n-space.
+    shape
+        shape of the image in the input orientation. map = ornt_mapping(input_ornt, output_ornt)
+
+    voxel_size
+        voxel_size of the image in the input orientation.
+
+    Returns:
+    --------
+    A : array (n+1, n+1)
+        Affine matrix of the transformation between input_ornt and output_ornt.
+
+    See Also:
+    ---------
+    nibabel.orientation
+    dipy.io.bvectxt.orientation_to_string
+    dipy.io.bvectxt.orientation_from_string
+    """
+    map = ornt_mapping(input_ornt, output_ornt)
+    if input_ornt.shape != output_ornt.shape:
+        raise ValueError("input_ornt and output_ornt must have the same shape")
+    affine = eye(len(input_ornt)+1)
+    affine[:3] = affine[map[:, 0]]
+    corner = asarray(voxel_size) * shape
+    affine[:3, 3] = (map[:, 1] < 0) * corner[map[:, 0]]
+    #multiply the rows of affine to get right sign
+    affine[:3, :3] *= map[:, 1:]
+    return affine
+
+def affine_from_fsl_mat_file(mat_affine, input_voxsz, output_voxsz):
+    """It takes the affine matrix from flirt (FSLdot) and the voxel size of the
+       input and output images and it returns the adjusted affine matrix for trackvis
+    """
+    input_voxsz = asarray(input_voxsz)
+    output_voxsz = asarray(output_voxsz)
+    shift = eye(4)
+    shift[:3,3] = -input_voxsz/2
+
+    affine = dot(mat_affine, shift)
+    affine[:3,3] += output_voxsz/2
+
+    return affine