tracking

Important Notes

Dipy uses affine matrices to represent the relationship between streamline points, which are defined as points in a continuous 3d space, and image voxels, which are typically arranged in a discrete 3d grid. Dipy uses a convention similar to nifti files to interpret these affine matrices. This convention is that the point at the center of voxel [i, j, k] is represented by the point [x, y, z] where [x, y, z, 1] = affine * [i, j, k, 1]. Also when the phrase “voxel coordinates” is used, it is understood to be the same as affine = eye(4).

As an example, let’s take a 2d image where the affine is:

[[1., 0., 0.],
 [0., 2., 0.],
 [0., 0., 1.]]

The pixels of an image with this affine would look something like:

A------------
|   |   |   |
| C |   |   |
|   |   |   |
----B--------
|   |   |   |
|   |   |   |
|   |   |   |
-------------
|   |   |   |
|   |   |   |
|   |   |   |
------------D

And the letters A-D represent the following points in “real world coordinates”:

A = [-.5, -1.]
B = [ .5,  1.]
C = [ 0.,  0.]
D = [ 2.5,  5.]

Streamlines

dipy.tracking.Streamlines

alias of ArraySequence

bench

dipy.tracking.bench(label='fast', verbose=1, extra_argv=None)

Run benchmarks for module using nose.

Parameters:

label{‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the benchmarks to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

• ‘fast’ - the default - which corresponds to the nosetests -A option of ‘not slow’.

• ‘full’ - fast (as above) and slow benchmarks as in the ‘no -A’ option to nosetests - this is the same as ‘’.

• None or ‘’ - run all tests.

• attribute_identifier - string passed directly to nosetests as ‘-A’.

verboseint, optional

Verbosity value for benchmark outputs, in the range 1-10. Default is 1.

extra_argvlist, optional

List with any extra arguments to pass to nosetests.

Returns:

successbool

Returns True if running the benchmarks works, False if an error occurred.

Notes

Benchmarks are like tests, but have names starting with “bench” instead of “test”, and can be found under the “benchmarks” sub-directory of the module.

Each NumPy module exposes bench in its namespace to run all benchmarks for it.

Examples

>>> success = np.lib.bench() 
Running benchmarks for numpy.lib
...
using 562341 items:
unique:
0.11
unique1d:
0.11
ratio: 1.0
nUnique: 56230 == 56230
...
OK

>>> success 
True

test

dipy.tracking.test(label='fast', verbose=1, extra_argv=None, doctests=False, coverage=False, raise_warnings=None, timer=False)

Run tests for module using nose.

Parameters:

label{‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the tests to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

• ‘fast’ - the default - which corresponds to the nosetests -A option of ‘not slow’.

• ‘full’ - fast (as above) and slow tests as in the ‘no -A’ option to nosetests - this is the same as ‘’.

• None or ‘’ - run all tests.

• attribute_identifier - string passed directly to nosetests as ‘-A’.

verboseint, optional

Verbosity value for test outputs, in the range 1-10. Default is 1.

extra_argvlist, optional

List with any extra arguments to pass to nosetests.

doctestsbool, optional

If True, run doctests in module. Default is False.

coveragebool, optional

If True, report coverage of NumPy code. Default is False. (This requires the coverage module).

raise_warningsNone, str or sequence of warnings, optional

This specifies which warnings to configure as ‘raise’ instead of being shown once during the test execution. Valid strings are:

• “develop” : equals (Warning,)

• “release” : equals (), do not raise on any warnings.

timerbool or int, optional

Timing of individual tests with nose-timer (which needs to be installed). If True, time tests and report on all of them. If an integer (say N), report timing results for N slowest tests.

Returns:

resultobject

Returns the result of running the tests as a nose.result.TextTestResult object.

Notes

Each NumPy module exposes test in its namespace to run all tests for it. For example, to run all tests for numpy.lib:

>>> np.lib.test() 

Examples

>>> result = np.lib.test() 
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s

OK

>>> result.errors 
[]
>>> result.knownfail 
[]

Streamlines

dipy.tracking.benchmarks.bench_streamline.Streamlines

alias of ArraySequence

assert_array_almost_equal

dipy.tracking.benchmarks.bench_streamline.assert_array_almost_equal(x, y, decimal=6, err_msg='', verbose=True)

Raises an AssertionError if two objects are not equal up to desired precision.

Note

It is recommended to use one of assert_allclose, assert_array_almost_equal_nulp or assert_array_max_ulp instead of this function for more consistent floating point comparisons.

The test verifies identical shapes and that the elements of actual and desired satisfy.

abs(desired-actual) < 1.5 * 10**(-decimal)

That is a looser test than originally documented, but agrees with what the actual implementation did up to rounding vagaries. An exception is raised at shape mismatch or conflicting values. In contrast to the standard usage in numpy, NaNs are compared like numbers, no assertion is raised if both objects have NaNs in the same positions.

Parameters:

xarray_like

The actual object to check.

yarray_like

The desired, expected object.

decimalint, optional

Desired precision, default is 6.

err_msgstr, optional

The error message to be printed in case of failure.

verbosebool, optional

If True, the conflicting values are appended to the error message.

Raises:

AssertionError

If actual and desired are not equal up to specified precision.

See also

assert_allclose

Compare two array_like objects for equality with desired relative and/or absolute precision.

assert_array_almost_equal_nulp, assert_array_max_ulp, assert_equal

Examples

the first assert does not raise an exception

>>> np.testing.assert_array_almost_equal([1.0,2.333,np.nan],
...                                      [1.0,2.333,np.nan])

>>> np.testing.assert_array_almost_equal([1.0,2.33333,np.nan],
...                                      [1.0,2.33339,np.nan], decimal=5)
Traceback (most recent call last):
    ...
AssertionError:
Arrays are not almost equal to 5 decimals

Mismatched elements: 1 / 3 (33.3%)
Max absolute difference: 6.e-05
Max relative difference: 2.57136612e-05
 x: array([1.     , 2.33333,     nan])
 y: array([1.     , 2.33339,     nan])

>>> np.testing.assert_array_almost_equal([1.0,2.33333,np.nan],
...                                      [1.0,2.33333, 5], decimal=5)
Traceback (most recent call last):
    ...
AssertionError:
Arrays are not almost equal to 5 decimals

x and y nan location mismatch:
 x: array([1.     , 2.33333,     nan])
 y: array([1.     , 2.33333, 5.     ])

assert_array_equal

dipy.tracking.benchmarks.bench_streamline.assert_array_equal(x, y, err_msg='', verbose=True, *, strict=False)

Raises an AssertionError if two array_like objects are not equal.

Given two array_like objects, check that the shape is equal and all elements of these objects are equal (but see the Notes for the special handling of a scalar). An exception is raised at shape mismatch or conflicting values. In contrast to the standard usage in numpy, NaNs are compared like numbers, no assertion is raised if both objects have NaNs in the same positions.

The usual caution for verifying equality with floating point numbers is advised.

Parameters:

xarray_like

The actual object to check.

yarray_like

The desired, expected object.

err_msgstr, optional

The error message to be printed in case of failure.

verbosebool, optional

If True, the conflicting values are appended to the error message.

strictbool, optional

If True, raise an AssertionError when either the shape or the data type of the array_like objects does not match. The special handling for scalars mentioned in the Notes section is disabled.

Raises:

AssertionError

If actual and desired objects are not equal.

See also

assert_allclose

Compare two array_like objects for equality with desired relative and/or absolute precision.

assert_array_almost_equal_nulp, assert_array_max_ulp, assert_equal

Notes

When one of x and y is a scalar and the other is array_like, the function checks that each element of the array_like object is equal to the scalar. This behaviour can be disabled with the strict parameter.

Examples

The first assert does not raise an exception:

>>> np.testing.assert_array_equal([1.0,2.33333,np.nan],
...                               [np.exp(0),2.33333, np.nan])

Assert fails with numerical imprecision with floats:

>>> np.testing.assert_array_equal([1.0,np.pi,np.nan],
...                               [1, np.sqrt(np.pi)**2, np.nan])
Traceback (most recent call last):
    ...
AssertionError:
Arrays are not equal

Mismatched elements: 1 / 3 (33.3%)
Max absolute difference: 4.4408921e-16
Max relative difference: 1.41357986e-16
 x: array([1.      , 3.141593,      nan])
 y: array([1.      , 3.141593,      nan])

Use assert_allclose or one of the nulp (number of floating point values) functions for these cases instead:

>>> np.testing.assert_allclose([1.0,np.pi,np.nan],
...                            [1, np.sqrt(np.pi)**2, np.nan],
...                            rtol=1e-10, atol=0)

As mentioned in the Notes section, assert_array_equal has special handling for scalars. Here the test checks that each value in x is 3:

>>> x = np.full((2, 5), fill_value=3)
>>> np.testing.assert_array_equal(x, 3)

Use strict to raise an AssertionError when comparing a scalar with an array:

>>> np.testing.assert_array_equal(x, 3, strict=True)
Traceback (most recent call last):
    ...
AssertionError:
Arrays are not equal

(shapes (2, 5), () mismatch)
 x: array([[3, 3, 3, 3, 3],
       [3, 3, 3, 3, 3]])
 y: array(3)

The strict parameter also ensures that the array data types match:

>>> x = np.array([2, 2, 2])
>>> y = np.array([2., 2., 2.], dtype=np.float32)
>>> np.testing.assert_array_equal(x, y, strict=True)
Traceback (most recent call last):
    ...
AssertionError:
Arrays are not equal

(dtypes int64, float32 mismatch)
 x: array([2, 2, 2])
 y: array([2., 2., 2.], dtype=float32)

bench_compress_streamlines

dipy.tracking.benchmarks.bench_streamline.bench_compress_streamlines()

bench_length

dipy.tracking.benchmarks.bench_streamline.bench_length()

bench_set_number_of_points

dipy.tracking.benchmarks.bench_streamline.bench_set_number_of_points()

generate_streamlines

dipy.tracking.benchmarks.bench_streamline.generate_streamlines(nb_streamlines, min_nb_points, max_nb_points, rng)

get_fnames

dipy.tracking.benchmarks.bench_streamline.get_fnames(name='small_64D')

Provide full paths to example or test datasets.

Parameters:

namestr

the filename/s of which dataset to return, one of:

• ‘small_64D’ small region of interest nifti,bvecs,bvals 64 directions

• ‘small_101D’ small region of interest nifti, bvecs, bvals 101 directions

• ‘aniso_vox’ volume with anisotropic voxel size as Nifti

• ‘fornix’ 300 tracks in Trackvis format (from Pittsburgh Brain Competition)

• ‘gqi_vectors’ the scanner wave vectors needed for a GQI acquisitions of 101 directions tested on Siemens 3T Trio

• ‘small_25’ small ROI (10x8x2) DTI data (b value 2000, 25 directions)

• ‘test_piesno’ slice of N=8, K=14 diffusion data

• ‘reg_c’ small 2D image used for validating registration

• ‘reg_o’ small 2D image used for validation registration

• ‘cb_2’ two vectorized cingulum bundles

Returns:

fnamestuple

filenames for dataset

Examples

>>> import numpy as np
>>> from dipy.io.image import load_nifti
>>> from dipy.data import get_fnames
>>> fimg, fbvals, fbvecs = get_fnames('small_101D')
>>> bvals=np.loadtxt(fbvals)
>>> bvecs=np.loadtxt(fbvecs).T
>>> data, affine = load_nifti(fimg)
>>> data.shape == (6, 10, 10, 102)
True
>>> bvals.shape == (102,)
True
>>> bvecs.shape == (102, 3)
True

length

dipy.tracking.benchmarks.bench_streamline.length()

Euclidean length of streamlines

Length is in mm only if streamlines are expressed in world coordinates.

Parameters:

streamlinesndarray or a list or dipy.tracking.Streamlines

If ndarray, must have shape (N,3) where N is the number of points of the streamline. If list, each item must be ndarray shape (Ni,3) where Ni is the number of points of streamline i. If dipy.tracking.Streamlines, its common_shape must be 3.

Returns:

lengthsscalar or ndarray shape (N,)

If there is only one streamline, a scalar representing the length of the streamline. If there are several streamlines, ndarray containing the length of every streamline.

Examples

>>> from dipy.tracking.streamline import length
>>> import numpy as np
>>> streamline = np.array([[1, 1, 1], [2, 3, 4], [0, 0, 0]])
>>> expected_length = np.sqrt([1+2**2+3**2, 2**2+3**2+4**2]).sum()
>>> length(streamline) == expected_length
True
>>> streamlines = [streamline, np.vstack([streamline, streamline[::-1]])]
>>> expected_lengths = [expected_length, 2*expected_length]
>>> lengths = [length(streamlines[0]), length(streamlines[1])]
>>> np.allclose(lengths, expected_lengths)
True
>>> length([])
0.0
>>> length(np.array([[1, 2, 3]]))
0.0

length_python

dipy.tracking.benchmarks.bench_streamline.length_python(xyz, along=False)

load_tractogram

dipy.tracking.benchmarks.bench_streamline.load_tractogram(filename, reference, to_space=Space.RASMM, to_origin=Origin.NIFTI, bbox_valid_check=True, trk_header_check=True)

Load the stateful tractogram from any format (trk/tck/vtk/vtp/fib/dpy)

Parameters:

filenamestring

Filename with valid extension

referenceNifti or Trk filename, Nifti1Image or TrkFile, Nifti1Header or

trk.header (dict), or ‘same’ if the input is a trk file. Reference that provides the spatial attribute. Typically a nifti-related object from the native diffusion used for streamlines generation

to_spaceEnum (dipy.io.stateful_tractogram.Space)

Space to which the streamlines will be transformed after loading

to_originEnum (dipy.io.stateful_tractogram.Origin)

Origin to which the streamlines will be transformed after loading

NIFTI standard, default (center of the voxel) TRACKVIS standard (corner of the voxel)

bbox_valid_checkbool

Verification for negative voxel coordinates or values above the volume dimensions. Default is True, to enforce valid file.

trk_header_checkbool

Verification that the reference has the same header as the spatial attributes as the input tractogram when a Trk is loaded

Returns:

outputStatefulTractogram

The tractogram to load (must have been saved properly)

measure

dipy.tracking.benchmarks.bench_streamline.measure(code_str, times=1, label=None)

Return elapsed time for executing code in the namespace of the caller.

The supplied code string is compiled with the Python builtin compile. The precision of the timing is 10 milli-seconds. If the code will execute fast on this timescale, it can be executed many times to get reasonable timing accuracy.

Parameters:

code_strstr

The code to be timed.

timesint, optional

The number of times the code is executed. Default is 1. The code is only compiled once.

labelstr, optional

A label to identify code_str with. This is passed into compile as the second argument (for run-time error messages).

Returns:

elapsedfloat

Total elapsed time in seconds for executing code_str times times.

Examples

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  
Time for a single execution :  0.005 s

set_number_of_points

dipy.tracking.benchmarks.bench_streamline.set_number_of_points()

Change the number of points of streamlines

(either by downsampling or upsampling)

Change the number of points of streamlines in order to obtain nb_points-1 segments of equal length. Points of streamlines will be modified along the curve.

Parameters:

streamlinesndarray or a list or dipy.tracking.Streamlines

If ndarray, must have shape (N,3) where N is the number of points of the streamline. If list, each item must be ndarray shape (Ni,3) where Ni is the number of points of streamline i. If dipy.tracking.Streamlines, its common_shape must be 3.

nb_pointsint

integer representing number of points wanted along the curve.

Returns:

new_streamlinesndarray or a list or dipy.tracking.Streamlines

Results of the downsampling or upsampling process.

Examples

>>> from dipy.tracking.streamline import set_number_of_points
>>> import numpy as np

One streamline, a semi-circle:

>>> theta = np.pi*np.linspace(0, 1, 100)
>>> x = np.cos(theta)
>>> y = np.sin(theta)
>>> z = 0 * x
>>> streamline = np.vstack((x, y, z)).T
>>> modified_streamline = set_number_of_points(streamline, 3)
>>> len(modified_streamline)
3

Multiple streamlines:

>>> streamlines = [streamline, streamline[::2]]
>>> new_streamlines = set_number_of_points(streamlines, 10)
>>> [len(s) for s in streamlines]
[100, 50]
>>> [len(s) for s in new_streamlines]
[10, 10]

set_number_of_points_python

dipy.tracking.benchmarks.bench_streamline.set_number_of_points_python(xyz, n_pols=3)

setup

dipy.tracking.benchmarks.bench_streamline.setup()

DirectionGetter

class dipy.tracking.direction_getter.DirectionGetter

Bases: object

Methods

generate_streamline
get_direction
initial_direction

__init__(*args, **kwargs)

generate_streamline()

get_direction()

initial_direction()

add_3vecs

dipy.tracking.distances.add_3vecs()

approx_polygon_track

dipy.tracking.distances.approx_polygon_track()

Fast and simple trajectory approximation algorithm by Eleftherios and Ian

It will reduce the number of points of the track by keeping intact the start and endpoints of the track and trying to remove as many points as possible without distorting much the shape of the track

Parameters:

xyzarray(N,3)

initial trajectory

alphafloat

smoothing parameter (<0.392 smoother, >0.392 rougher) if the trajectory was a smooth circle then with alpha =0.393 ~=pi/8. the circle would be approximated with an decahexagon if alpha = 0.7853 ~=pi/4. with an octagon.

Returns:

characteristic_points: list of M array(3,) points

Notes

Assuming that a good approximation for a circle is an octagon then that means that the points of the octagon will have angle alpha = 2*pi/8 = pi/4 . We calculate the angle between every two neighbour segments of a trajectory and if the angle is higher than pi/4 we choose that point as a characteristic point otherwise we move at the next point.

Examples

Approximating a helix:

>>> t=np.linspace(0,1.75*2*np.pi,100)
>>> x = np.sin(t)
>>> y = np.cos(t)
>>> z = t
>>> xyz=np.vstack((x,y,z)).T
>>> xyza = approx_polygon_track(xyz)
>>> len(xyza) < len(xyz)
True

approximate_mdl_trajectory

dipy.tracking.distances.approximate_mdl_trajectory()

Implementation of Lee et al Approximate Trajectory Partitioning Algorithm

This is base on the minimum description length principle

Parameters:

xyzarray(N,3)

initial trajectory

alphafloat

smoothing parameter (>1 smoother, <1 rougher)

Returns:

characteristic_pointslist of M array(3,) points

bundles_distances_mam

dipy.tracking.distances.bundles_distances_mam()

Calculate distances between list of tracks A and list of tracks B

Parameters:

tracksAsequence

of tracks as arrays, shape (N1,3) .. (Nm,3)

tracksBsequence

of tracks as arrays, shape (N1,3) .. (Nm,3)

metricstr

‘avg’, ‘min’, ‘max’

Returns:

DMarray, shape (len(tracksA), len(tracksB))

distances between tracksA and tracksB according to metric

See also

dipy.tracking.streamline.set_number_of_points

bundles_distances_mdf

dipy.tracking.distances.bundles_distances_mdf()

Calculate distances between list of tracks A and list of tracks B

All tracks need to have the same number of points

Parameters:

tracksAsequence

of tracks as arrays, [(N,3) .. (N,3)]

tracksBsequence

of tracks as arrays, [(N,3) .. (N,3)]

Returns:

DMarray, shape (len(tracksA), len(tracksB))

distances between tracksA and tracksB according to metric

See also

dipy.tracking.streamline.set_number_of_points

cut_plane

dipy.tracking.distances.cut_plane()

Extract divergence vectors and points of intersection between planes normal to the reference fiber and other tracks

Parameters:

trackssequence

of tracks as arrays, shape (N1,3) .. (Nm,3)

refarray, shape (N,3)

reference track

Returns:

hitssequence

list of points and rcds (radial coefficient of divergence)

Notes

The orthogonality relationship np.inner(hits[p][q][0:3]-ref[p+1],ref[p+2]-ref[r][p+1]) will hold throughout for every point q in the hits plane at point (p+1) on the reference track.

Examples

>>> refx = np.array([[0,0,0],[1,0,0],[2,0,0],[3,0,0]],dtype='float32')
>>> bundlex = [np.array([[0.5,1,0],[1.5,2,0],[2.5,3,0]],dtype='float32')]
>>> res = cut_plane(bundlex,refx)
>>> len(res)
2
>>> np.allclose(res[0], np.array([[1., 1.5, 0., 0.70710683, 0. ]]))
True
>>> np.allclose(res[1], np.array([[2., 2.5, 0., 0.70710677, 0. ]]))
True

inner_3vecs

dipy.tracking.distances.inner_3vecs()

intersect_segment_cylinder

dipy.tracking.distances.intersect_segment_cylinder()

Intersect Segment S(t) = sa +t(sb-sa), 0 <=t<= 1 against cylinder specified by p,q and r

See p.197 from Real Time Collision Detection by C. Ericson

Examples

Define cylinder using a segment defined by

>>> p=np.array([0,0,0],dtype=np.float32)
>>> q=np.array([1,0,0],dtype=np.float32)
>>> r=0.5

Define segment

>>> sa=np.array([0.5,1 ,0],dtype=np.float32)
>>> sb=np.array([0.5,-1,0],dtype=np.float32)

Intersection

>>> intersect_segment_cylinder(sa, sb, p, q, r)
(1.0, 0.25, 0.75)

larch_3merge

dipy.tracking.distances.larch_3merge()

Reassign tracks to existing clusters by merging clusters that their representative tracks are not very distant i.e. less than sqd_thr. Using tracks consisting of 3 points (first, mid and last). This is necessary after running larch_fast_split after multiple split in different levels (squared thresholds) as some of them have created independent clusters.

Parameters:

Cgraph with clusters

of indices 3tracks (tracks consisting of 3 points only)

sqd_trh: float

squared euclidean distance threshold

Returns:

Cdict

a tree graph containing the clusters

larch_3split

dipy.tracking.distances.larch_3split()

Generate a first pass clustering using 3 points on the tracks only.

Parameters:

trackssequence

of tracks as arrays, shape (N1,3) .. (Nm,3), where 3 points are (first, mid and last)

indicesNone or sequence, optional

Sequence of integer indices of tracks

trhfloat, optional

squared euclidean distance threshold

Returns:

Cdict

A tree graph containing the clusters.

Notes

If a 3 point track (3track) is far away from all clusters then add a new cluster and assign this 3track as the rep(representative) track for the new cluster. Otherwise the rep 3track of each cluster is the average track of the cluster.

Examples

>>> tracks=[np.array([[0,0,0],[1,0,0,],[2,0,0]],dtype=np.float32),
...         np.array([[3,0,0],[3.5,1,0],[4,2,0]],dtype=np.float32),
...         np.array([[3.2,0,0],[3.7,1,0],[4.4,2,0]],dtype=np.float32),
...         np.array([[3.4,0,0],[3.9,1,0],[4.6,2,0]],dtype=np.float32),
...         np.array([[0,0.2,0],[1,0.2,0],[2,0.2,0]],dtype=np.float32),
...         np.array([[2,0.2,0],[1,0.2,0],[0,0.2,0]],dtype=np.float32),
...         np.array([[0,0,0],[0,1,0],[0,2,0]],dtype=np.float32),
...         np.array([[0.2,0,0],[0.2,1,0],[0.2,2,0]],dtype=np.float32),
...         np.array([[-0.2,0,0],[-0.2,1,0],[-0.2,2,0]],dtype=np.float32)]
>>> C = larch_3split(tracks, None, 0.5)

Here is an example of how to visualize the clustering above:

from dipy.viz import window, actor
scene = window.Scene()
scene.add(actor.line(tracks,window.colors.red))
window.show(scene)
for c in C:
    color=np.random.rand(3)
    for i in C[c]['indices']:
        scene.add(actor.line(tracks[i],color))
window.show(scene)
for c in C:
    scene.add(actor.line(C[c]['rep3']/C[c]['N'],
                         window.colors.white))
window.show(scene)

lee_angle_distance

dipy.tracking.distances.lee_angle_distance()

Calculates angle distance metric for the distance between two line segments

Based on Lee , Han & Whang SIGMOD07.

This function assumes that norm(end0-start0)>norm(end1-start1) i.e. that the first segment will be bigger than the second one.

Parameters:

start0float array(3,)

end0float array(3,)

start1float array(3,)

end1float array(3,)

Returns:

angle_distancefloat

Notes

l_0 = np.inner(end0-start0,end0-start0) l_1 = np.inner(end1-start1,end1-start1)

cos_theta_squared = np.inner(end0-start0,end1-start1)**2/ (l_0*l_1) return np.sqrt((1-cos_theta_squared)*l_1)

Examples

>>> lee_angle_distance([0,0,0],[1,0,0],[3,4,5],[5,4,3])
2.0

lee_perpendicular_distance

dipy.tracking.distances.lee_perpendicular_distance()

Calculates perpendicular distance metric for the distance between two line segments

Based on Lee , Han & Whang SIGMOD07.

This function assumes that norm(end0-start0)>norm(end1-start1) i.e. that the first segment will be bigger than the second one.

Parameters:

start0float array(3,)

end0float array(3,)

start1float array(3,)

end1float array(3,)

Returns:

perpendicular_distance: float

Notes

l0 = np.inner(end0-start0,end0-start0) l1 = np.inner(end1-start1,end1-start1)

k0=end0-start0

u1 = np.inner(start1-start0,k0)/l0 u2 = np.inner(end1-start0,k0)/l0

ps = start0+u1*k0 pe = start0+u2*k0

lperp1 = np.sqrt(np.inner(ps-start1,ps-start1)) lperp2 = np.sqrt(np.inner(pe-end1,pe-end1))

if lperp1+lperp2 > 0.:

return (lperp1**2+lperp2**2)/(lperp1+lperp2)

else:

return 0.

Examples

>>> d = lee_perpendicular_distance([0,0,0],[1,0,0],[3,4,5],[5,4,3])
>>> print('%.6f' % d)
5.787888

local_skeleton_clustering

dipy.tracking.distances.local_skeleton_clustering()

Efficient tractography clustering

Every track can needs to have the same number of points. Use dipy.tracking.metrics.downsample to restrict the number of points

Parameters:

trackssequence

of tracks as arrays, shape (N,3) .. (N,3) where N=points

d_thrfloat

average euclidean distance threshold

Returns:

Cdict

Clusters.

See also

dipy.tracking.metrics.downsample

Notes

The distance calculated between two tracks:

t_1       t_2

0*   a    *0
  \       |
   \      |
   1*     |
    |  b  *1
    |      \
    2*      \
        c    *2

is equal to \((a+b+c)/3\) where \(a\) the euclidean distance between t_1[0] and t_2[0], \(b\) between t_1[1] and t_2[1] and \(c\) between t_1[2] and t_2[2]. Also the same with t2 flipped (so t_1[0] compared to t_2[2] etc).

Visualization:

It is possible to visualize the clustering C from the example above using the dipy.viz module:

from dipy.viz import window, actor
scene = window.Scene()
for c in C:
    color=np.random.rand(3)
    for i in C[c]['indices']:
        scene.add(actor.line(tracks[i],color))
window.show(scene)

Examples

>>> tracks=[np.array([[0,0,0],[1,0,0,],[2,0,0]]),
...         np.array([[3,0,0],[3.5,1,0],[4,2,0]]),
...         np.array([[3.2,0,0],[3.7,1,0],[4.4,2,0]]),
...         np.array([[3.4,0,0],[3.9,1,0],[4.6,2,0]]),
...         np.array([[0,0.2,0],[1,0.2,0],[2,0.2,0]]),
...         np.array([[2,0.2,0],[1,0.2,0],[0,0.2,0]]),
...         np.array([[0,0,0],[0,1,0],[0,2,0]])]
>>> C = local_skeleton_clustering(tracks, d_thr=0.5)

local_skeleton_clustering_3pts

dipy.tracking.distances.local_skeleton_clustering_3pts()

Does a first pass clustering

Every track can only have 3 pts neither less or more. Use dipy.tracking.metrics.downsample to restrict the number of points

Parameters:

trackssequence

of tracks as arrays, shape (N,3) .. (N,3) where N=3

d_thrfloat

Average euclidean distance threshold

Returns:

Cdict

Clusters.

Notes

It is possible to visualize the clustering C from the example above using the fvtk module:

r=fvtk.ren()
for c in C:
    color=np.random.rand(3)
    for i in C[c]['indices']:
        fvtk.add(r,fos.line(tracks[i],color))
fvtk.show(r)

Examples

>>> tracks=[np.array([[0,0,0],[1,0,0,],[2,0,0]]),
...         np.array([[3,0,0],[3.5,1,0],[4,2,0]]),
...         np.array([[3.2,0,0],[3.7,1,0],[4.4,2,0]]),
...         np.array([[3.4,0,0],[3.9,1,0],[4.6,2,0]]),
...         np.array([[0,0.2,0],[1,0.2,0],[2,0.2,0]]),
...         np.array([[2,0.2,0],[1,0.2,0],[0,0.2,0]]),
...         np.array([[0,0,0],[0,1,0],[0,2,0]])]
>>> C=local_skeleton_clustering_3pts(tracks,d_thr=0.5)

mam_distances

dipy.tracking.distances.mam_distances()

Min/Max/Mean Average Minimum Distance between tracks xyz1 and xyz2

Based on the metrics in Zhang, Correia, Laidlaw 2008 http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4479455 which in turn are based on those of Corouge et al. 2004

Parameters:

xyz1array, shape (N1,3), dtype float32

xyz2array, shape (N2,3), dtype float32

arrays representing x,y,z of the N1 and N2 points of two tracks

metrics{‘avg’,’min’,’max’,’all’}

Metric to calculate. {‘avg’,’min’,’max’} return a scalar. ‘all’ returns a tuple

Returns:

avg_mcdfloat

average_mean_closest_distance

min_mcdfloat

minimum_mean_closest_distance

max_mcdfloat

maximum_mean_closest_distance

Notes

Algorithmic description

Let’s say we have curves A and B.

For every point in A calculate the minimum distance from every point in B stored in minAB

For every point in B calculate the minimum distance from every point in A stored in minBA

find average of minAB stored as avg_minAB find average of minBA stored as avg_minBA

if metric is ‘avg’ then return (avg_minAB + avg_minBA)/2.0 if metric is ‘min’ then return min(avg_minAB,avg_minBA) if metric is ‘max’ then return max(avg_minAB,avg_minBA)

minimum_closest_distance

dipy.tracking.distances.minimum_closest_distance()

Find the minimum distance between two curves xyz1, xyz2

Parameters:

xyz1array, shape (N1,3), dtype float32

xyz2array, shape (N2,3), dtype float32

arrays representing x,y,z of the N1 and N2 points of two tracks

Returns:

mdminimum distance

Notes

Algorithmic description

Let’s say we have curves A and B

for every point in A calculate the minimum distance from every point in B stored in minAB for every point in B calculate the minimum distance from every point in A stored in minBA find min of minAB stored in min_minAB find min of minBA stored in min_minBA

Then return (min_minAB + min_minBA)/2.0

most_similar_track_mam

dipy.tracking.distances.most_similar_track_mam()

Find the most similar track in a bundle using distances calculated from Zhang et. al 2008.

Parameters:

trackssequence

of tracks as arrays, shape (N1,3) .. (Nm,3)

metricstr

‘avg’, ‘min’, ‘max’

Returns:

siint

index of the most similar track in tracks. This can be used as a reference track for a bundle.

sarray, shape (len(tracks),)

similarities between tracks[si] and the rest of the tracks in the bundle

Notes

A vague description of this function is given below:

for (i,j) in tracks_combinations_of_2:

calculate the mean_closest_distance from i to j (mcd_i) calculate the mean_closest_distance from j to i (mcd_j)

if ‘avg’:

s holds the average similarities

if ‘min’:

s holds the minimum similarities

if ‘max’:

s holds the maximum similarities

si holds the index of the track with min {avg,min,max} average metric

mul_3vec

dipy.tracking.distances.mul_3vec()

mul_3vecs

dipy.tracking.distances.mul_3vecs()

norm_3vec

dipy.tracking.distances.norm_3vec()

Euclidean (L2) norm of length 3 vector

Parameters:

vecarray-like shape (3,)

Returns:

normfloat

Euclidean norm

normalized_3vec

dipy.tracking.distances.normalized_3vec()

Return normalized 3D vector

Vector divided by Euclidean (L2) norm

Parameters:

vecarray-like shape (3,)

Returns:

vec_outarray shape (3,)

point_segment_sq_distance

dipy.tracking.distances.point_segment_sq_distance()

Calculate the squared distance from a point c to a finite line segment ab.

Examples

>>> a=np.array([0,0,0], dtype=np.float32)
>>> b=np.array([1,0,0], dtype=np.float32)
>>> c=np.array([0,1,0], dtype=np.float32)
>>> point_segment_sq_distance(a, b, c)
1.0
>>> c = np.array([0,3,0], dtype=np.float32)
>>> point_segment_sq_distance(a,b,c)
9.0
>>> c = np.array([-1,1,0], dtype=np.float32)
>>> point_segment_sq_distance(a, b, c)
2.0

point_track_sq_distance_check

dipy.tracking.distances.point_track_sq_distance_check()

Check if square distance of track from point is smaller than threshold

Parameters:

track: array,float32, shape (N,3)

point: array,float32, shape (3,)

sq_dist_thr: double, threshold

Returns:

bool: True, if sq_distance <= sq_dist_thr, otherwise False.

Examples

>>> t=np.random.rand(10,3).astype(np.float32)
>>> p=np.array([0.5,0.5,0.5],dtype=np.float32)
>>> point_track_sq_distance_check(t,p,2**2)
True
>>> t=np.array([[0,0,0],[1,1,1],[2,2,2]],dtype='f4')
>>> p=np.array([-1,-1.,-1],dtype='f4')
>>> point_track_sq_distance_check(t,p,.2**2)
False
>>> point_track_sq_distance_check(t,p,2**2)
True

sub_3vecs

dipy.tracking.distances.sub_3vecs()

track_dist_3pts

dipy.tracking.distances.track_dist_3pts()

Calculate the euclidean distance between two 3pt tracks

Both direct and flip distances are calculated but only the smallest is returned

Parameters:

aarray, shape (3,3)

a three point track

barray, shape (3,3)

a three point track

Returns:

dist :float

Examples

>>> a = np.array([[0,0,0],[1,0,0,],[2,0,0]])
>>> b = np.array([[3,0,0],[3.5,1,0],[4,2,0]])
>>> c = track_dist_3pts(a, b)
>>> print('%.6f' % c)
2.721573

track_roi_intersection_check

dipy.tracking.distances.track_roi_intersection_check()

Check if a track is intersecting a region of interest

Parameters:

track: array,float32, shape (N,3)

roi: array,float32, shape (M,3)

sq_dist_thr: double, threshold, check squared euclidean distance from every roi point

Returns:

bool: True, if sq_distance <= sq_dist_thr, otherwise False.

Examples

>>> roi=np.array([[0,0,0],[1,0,0],[2,0,0]],dtype='f4')
>>> t=np.array([[0,0,0],[1,1,1],[2,2,2]],dtype='f4')
>>> track_roi_intersection_check(t,roi,1)
True
>>> track_roi_intersection_check(t,np.array([[10,0,0]],dtype='f4'),1)
False