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

tracks : sequence
of tracks as arrays, shape (N,3) .. (N,3) where N=points
d_thr : float
average euclidean distance threshold
Returns

C : dict
Clusters. 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)
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)
See Also

dipy.tracking.metrics.downsample
.
This is an implementation of the Linear Fascicle Evaluation (LiFE) algorithm
described in:
Pestilli, F., Yeatman, J, Rokem, A. Kay, K. and Wandell B.A. (2014). Validation
and statistical inference in living connectomes. Nature Methods 11:
10581063. doi:10.1038/nmeth.3098
Parameters
seed
array, float64 shape (3,)
Point where the tracking starts.
ref
cnp.npy_intp int
Index of peak to follow first.
qa
array, float64 shape (X, Y, Z, Np)
Anisotropy matrix, where Np is the number of maximum allowed peaks.
ind
array, float64 shape(x, y, z, Np)
Index of the track orientation.
odf_vertices
double array shape (N, 3)
Sampling directions on the sphere.
qa_thr
float
Threshold for QA, we want everything higher than this threshold.
ang_thr
float
Angle threshold, we only select fiber orientation within this range.
:17: (WARNING/2) Definition list ends without a blank line; unexpected unindent.
step_sz : double
total_weight : double
max_points : cnp.npy_intp
AnatomicallyConstrained Tractography (ACT) stopping criterion from [1]_. This implements the use of partial volume fraction (PVE) maps to determine when the tracking stops. The proposed ([1]_) method that cuts streamlines going through subcortical gray matter regions is not implemented here. The backtracking technique for streamlines reaching INVALIDPOINT is not implemented either. cdef: double interp_out_double[1] double[:] interp_out_view = interp_out_view double[:, :, :] include_map, exclude_map.
Various tools related to creating and working with streamlines.
This module provides tools for targeting streamlines using ROIs, for making
connectivity matrices from whole brain fiber tracking and some other tools that
allow streamlines to interact with image data.
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:
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
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.
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.
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.
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.
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.
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
———
tracks : sequence
of tracks as arrays, shape (N,3) .. (N,3) where N=points
d_thr : float
average euclidean distance threshold
Returns
——
C : dict
Clusters.
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)
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)
See Also
——–
dipy.tracking.metrics.downsample
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)
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
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
A collection of streamlines, each n by 3, with n being the number of
nodes in the fiber.
kernelKernel object
A diffusion kernel object created from EnhancementKernel.
min_fiberlengthint
Fibers with fewer points than minimum_length are excluded from FBC
computation.
max_windowsizeint
The maximal window size used to calculate the average LFBC region
num_threadsint, optional
Number of threads to be used for OpenMP parallelization. If None
(default) the value of OMP_NUM_THREADS environment variable is used
if it is set, otherwise all available threads are used. If < 0 the
maximal number of threads minus num_threads + 1 is used (enter 1
to use as many threads as possible). 0 raises an error.
[Meesters2016_HBM] S. Meesters, G. Sanguinetti, E. Garyfallidis,
J. Portegies, P. Ossenblok, R. Duits. (2016) Cleaning
output of tractography via fiber to bundle coherence,
a new open source implementation. Human Brain Mapping
conference 2016.
[Portegies2015b] J. Portegies, R. Fick, G. Sanguinetti, S. Meesters,
G.Girard, and R. Duits. (2015) Improving Fiber Alignment
in HARDI by Combining Contextual PDE flow with
Constrained Spherical Deconvolution. PLoS One.
The output contains:
1) a collection of streamlines, each n by 3, with n
being the number of nodes in the fiber that remain after filtering
2) the r,g,b values of the local fiber to bundle coherence (LFBC)
3) the relative fiber to bundle coherence (RFBC)
The n data points of dimension m to be indexed. This array is
not copied unless this is necessary to produce a contiguous
array of doubles, and so modifying this data will result in
bogus results. The data are also copied if the kdtree is built
with copy_data=True.
leafsizepositive int, optional
The number of points at which the algorithm switches over to
bruteforce. Default: 10.
compact_nodesbool, optional
If True, the kdtree is built to shrink the hyperrectangles to
the actual data range. This usually gives a more compact tree that
is robust against degenerated input data and gives faster queries
at the expense of longer build time. Default: True.
copy_databool, optional
If True the data is always copied to protect the kdtree against
data corruption. Default: False.
balanced_treebool, optional
If True, the median is used to split the hyperrectangles instead of
the midpoint. This usually gives a more compact tree and
faster queries at the expense of longer build time. Default: True.
boxsizearray_like or scalar, optional
Apply a md toroidal topology to the KDTree.. The topology is generated
by \(x_i + n_i L_i\) where \(n_i\) are integers and \(L_i\)
is the boxsize along ith dimension. The input data shall be wrapped
into \([0, L_i)\). A ValueError is raised if any of the data is
outside of this bound.
The algorithm used is described in Maneewongvatana and Mount 1999.
The general idea is that the kdtree is a binary tree, each of whose
nodes represents an axisaligned hyperrectangle. Each node specifies
an axis and splits the set of points based on whether their coordinate
along that axis is greater than or less than a particular value.
During construction, the axis and splitting point are chosen by the
“sliding midpoint” rule, which ensures that the cells do not all
become long and thin.
The tree can be queried for the r closest neighbors of any given point
(optionally returning only those within some maximum distance of the
point). It can also be queried, with a substantial gain in efficiency,
for the r approximate closest neighbors.
For large dimensions (20 is already large) do not expect this to run
significantly faster than brute force. Highdimensional nearestneighbor
queries are a substantial open problem in computer science.
The n data points of dimension m to be indexed. This array is
not copied unless this is necessary to produce a contiguous
array of doubles. The data are also copied if the kdtree is built
with copy_data=True.
leafsizepositive int
The number of points at which the algorithm switches over to
bruteforce.
mint
The dimension of a single datapoint.
nint
The number of data points.
maxesndarray, shape (m,)
The maximum value in each dimension of the n data points.
minsndarray, shape (m,)
The minimum value in each dimension of the n data points.
The other tree to draw points from, can be the same tree as self.
rfloat or onedimensional array of floats
The radius to produce a count for. Multiple radii are searched with
a single tree traversal.
If the count is noncumulative(cumulative=False), r defines
the edges of the bins, and must be nondecreasing.
pfloat, optional
1<=p<=infinity.
Which Minkowski pnorm to use.
Default 2.0.
A finite large p may cause a ValueError if overflow can occur.
weightstuple, array_like, or None, optional
If None, the paircounting is unweighted.
If given as a tuple, weights[0] is the weights of points in
self, and weights[1] is the weights of points in other;
either can be None to indicate the points are unweighted.
If given as an array_like, weights is the weights of points in
self and other. For this to make sense, self and
other must be the same tree. If self and other are two
different trees, a ValueError is raised.
Default: None
New in version 1.6.0.
cumulativebool, optional
Whether the returned counts are cumulative. When cumulative is set
to False the algorithm is optimized to work with a large number
of bins (>10) specified by r. When cumulative is set to
True, the algorithm is optimized to work with a small number of
r. Default: True
The number of pairs. For unweighted counts, the result is integer.
For weighted counts, the result is float.
If cumulative is False, result[i] contains the counts with
(infifi==0elser[i1])<R<=r[i]
Paircounting is the basic operation used to calculate the two point
correlation functions from a data set composed of position of objects.
Two point correlation function measures the clustering of objects and
is widely used in cosmology to quantify the large scale structure
in our Universe, but it may be useful for data analysis in other fields
where selfsimilar assembly of objects also occur.
The LandySzalay estimator for the two point correlation function of
D measures the clustering signal in D. [2]
For example, given the position of two sets of objects,
objects D (data) contains the clustering signal, and
where the brackets represents counting pairs between two data sets
in a finite bin around r (distance), corresponding to setting
cumulative=False, and f=float(len(D))/float(len(R)) is the
ratio between number of objects from data and random.
The algorithm implemented here is loosely based on the dualtree
algorithm described in [1]_. We switch between two different
paircumulation scheme depending on the setting of cumulative.
The computing time of the method we use when for
cumulative==False does not scale with the total number of bins.
The algorithm for cumulative==True scales linearly with the
number of bins, though it is slightly faster when only
1 or 2 bins are used. [5].
As an extension to the naive paircounting,
weighted paircounting counts the product of weights instead
of number of pairs.
Weighted paircounting is used to estimate marked correlation functions
([3], section 2.2),
or to properly calculate the average of data per distance bin
(e.g. [4], section 2.1 on redshift).
Either the number of nearest neighbors to return, or a list of the
kth nearest neighbors to return, starting from 1.
epsnonnegative float, optional
Return approximate nearest neighbors; the kth returned value
is guaranteed to be no further than (1+eps) times the
distance to the real kth nearest neighbor.
pfloat, 1<=p<=infinity, optional
Which Minkowski pnorm to use.
1 is the sumofabsolutevalues distance (“Manhattan” distance).
2 is the usual Euclidean distance.
infinity is the maximumcoordinatedifference distance.
A large, finite p may cause a ValueError if overflow can occur.
distance_upper_boundnonnegative float, optional
Return only neighbors within this distance. This is used to prune
tree searches, so if you are doing a series of nearestneighbor
queries, it may help to supply the distance to the nearest neighbor
of the most recent point.
workersint, optional
Number of workers to use for parallel processing. If 1 is given
all CPU threads are used. Default: 1.
The distances to the nearest neighbors.
If x has shape tuple+(self.m,), then d has shape
tuple+(k,).
When k == 1, the last dimension of the output is squeezed.
Missing neighbors are indicated with infinite distances.
Hits are sorted by distance (nearest first).
Changed in version 1.9.0: Previously if k=None, then d was an object array of
shape tuple, containing lists of distances. This behavior
has been removed, use query_ball_point instead.
iinteger or array of integers
The index of each neighbor in self.data.
i is the same shape as d.
Missing neighbors are indicated with self.n.
The radius of points to return, must broadcast to the length of x.
pfloat, optional
Which Minkowski pnorm to use. Should be in the range [1, inf].
A finite large p may cause a ValueError if overflow can occur.
epsnonnegative float, optional
Approximate search. Branches of the tree are not explored if their
nearest points are further than r/(1+eps), and branches are
added in bulk if their furthest points are nearer than
r*(1+eps).
workersint, optional
Number of jobs to schedule for parallel processing. If 1 is given
all processors are used. Default: 1.
New in version 1.6.0.
return_sortedbool, optional
Sorts returned indicies if True and does not sort them if False. If
None, does not sort single point queries, but does sort
multipoint queries which was the behavior before this option
was added.
New in version 1.6.0.
return_lengthbool, optional
Return the number of points inside the radius instead of a list
of the indices.
If x is a single point, returns a list of the indices of the
neighbors of x. If x is an array of points, returns an object
array of shape tuple containing lists of neighbors.
If you have many points whose neighbors you want to find, you may save
substantial amounts of time by putting them in a KDTree and using
query_ball_tree.
Which Minkowski norm to use. p has to meet the condition
1<=p<=infinity.
epsfloat, optional
Approximate search. Branches of the tree are not explored
if their nearest points are further than r/(1+eps), and
branches are added in bulk if their furthest points are nearer
than r*(1+eps). eps has to be nonnegative.
Which Minkowski norm to use. p has to meet the condition
1<=p<=infinity.
epsfloat, optional
Approximate search. Branches of the tree are not explored
if their nearest points are further than r/(1+eps), and
branches are added in bulk if their furthest points are nearer
than r*(1+eps). eps has to be nonnegative.
output_typestring, optional
Choose the output container, ‘set’ or ‘ndarray’. Default: ‘set’
Sparse matrix representing the results in “dictionary of keys”
format. If a dict is returned the keys are (i,j) tuples of indices.
If output_type is ‘ndarray’ a record array with fields ‘i’, ‘j’,
and ‘v’ is returned,
classdipy.tracking.fbcmeasures.interp1d(x, y, kind='linear', axis=1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)#
Bases: _Interpolator1D
Interpolate a 1D function.
x and y are arrays of values used to approximate some function f:
y=f(x). This class returns a function whose call method uses
interpolation to find the value of new points.
A ND array of real values. The length of y along the interpolation
axis must be equal to the length of x.
kindstr or int, optional
Specifies the kind of interpolation as a string or as an integer
specifying the order of the spline interpolator to use.
The string has to be one of ‘linear’, ‘nearest’, ‘nearestup’, ‘zero’,
‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’, or ‘next’. ‘zero’,
‘slinear’, ‘quadratic’ and ‘cubic’ refer to a spline interpolation of
zeroth, first, second or third order; ‘previous’ and ‘next’ simply
return the previous or next value of the point; ‘nearestup’ and
‘nearest’ differ when interpolating halfintegers (e.g. 0.5, 1.5)
in that ‘nearestup’ rounds up and ‘nearest’ rounds down. Default
is ‘linear’.
axisint, optional
Specifies the axis of y along which to interpolate.
Interpolation defaults to the last axis of y.
copybool, optional
If True, the class makes internal copies of x and y.
If False, references to x and y are used. The default is to copy.
bounds_errorbool, optional
If True, a ValueError is raised any time interpolation is attempted on
a value outside of the range of x (where extrapolation is
necessary). If False, out of bounds values are assigned fill_value.
By default, an error is raised unless fill_value="extrapolate".
fill_valuearraylike or (arraylike, array_like) or “extrapolate”, optional
if a ndarray (or float), this value will be used to fill in for
requested points outside of the data range. If not provided, then
the default is NaN. The arraylike must broadcast properly to the
dimensions of the noninterpolation axes.
If a twoelement tuple, then the first element is used as a
fill value for x_new<x[0] and the second element is used for
x_new>x[1]. Anything that is not a 2element tuple (e.g.,
list or ndarray, regardless of shape) is taken to be a single
arraylike argument meant to be used for both bounds as
below,above=fill_value,fill_value. Using a twoelement tuple
or ndarray requires bounds_error=False.
New in version 0.17.0.
If “extrapolate”, then points outside the data range will be
extrapolated.
New in version 0.17.0.
assume_sortedbool, optional
If False, values of x can be in any order and they are sorted first.
If True, x has to be an array of monotonically increasing values.
Calling interp1d with NaNs present in input values results in
undefined behaviour.
Input values x and y must be convertible to float values like
int or float.
If the values in x are not unique, the resulting behavior is
undefined and specific to the choice of kind, i.e., changing
kind will change the behavior for duplicates.
Determine the effective number of threads to be used for OpenMP calls
For num_threads=None,
 if the OMP_NUM_THREADS environment variable is set, return that
value
 otherwise, return the maximum number of cores retrieved by
openmp.opm_get_num_procs().
For num_threads>0, return this value.
For num_threads<0, return the maximal number of threads minus
num_threads+1. In particular num_threads=1 will use as
many threads as there are available cores on the machine.
track2trackarray (N,2) where N is len(indices) of int
it shows the correspondence in the following way:
the first column is the current index in tracks1
the second column is the corresponding index in tracks2
To find the corresponding tracks we use mam_distances with ‘avg’ option.
Then we calculate the argmin of all the calculated distances and return it
for every index. (See 3rd column of arr in the example given below.)
of int showing the correspondence in th following way
the first column is the current index of tracks1
the second column is the corresponding index in tracks2
To find the corresponding tracks we use mam_distances with ‘avg’ option.
Then we calculate the argmin of all the calculated distances and return it
for every index. (See 3rd column of arr in the example given below.)
The gradient table on which the signal is calculated.
evalsarraylike of 3 items, optional
The eigenvalues of the canonical tensor to use in calculating the
signal.
spheredipy.core.Sphere class instance, optional
The discrete sphere to use as an approximation for the continuous
sphere on which the signal is represented. If integer  we will use
an instance of one of the symmetric spheres cached in
dps.get_sphere. If a ‘dipy.core.Sphere’ class instance is
provided, we will use this object. Default: the dipy.data
symmetric sphere with 724 vertices
The mapping from voxel coordinates to streamline points.
The voxel_to_rasmm matrix, typically from a NIFTI file.
evalsarraylike (optional)
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, optional
Whether to approximate (and cache) the signal on a discrete
sphere. This may confer a significant speedup in setting up the
problem, but is not as accurate. If False, we use the exact
gradients along the streamlines to calculate the matrix, instead of
an approximation.
Set up the necessary components for the LiFE model: the matrix of
fibercontributions to the DWI signal, and the coordinates of voxels
for which the equations will be solved
The mapping from voxel coordinates to streamline points.
The voxel_to_rasmm matrix, typically from a NIFTI file.
evalsarraylike (3 items, optional)
The eigenvalues of the canonical tensor used as a response
function. Default:[0.001, 0, 0].
sphere: dipy.core.Sphere instance, optional
Whether to approximate (and cache) the signal on a discrete
sphere. This may confer a significant speedup in setting up the
problem, but is not as accurate. If False, we use the exact
gradients along the streamlines to calculate the matrix, instead of
an approximation. Defaults to use the 724vertex symmetric sphere
from dipy.data
The gradient is computed using central differences in the interior
and first differences at the boundaries. The returned gradient hence has
the same shape as the input array.
This is a simplified implementation of gradient that is part of numpy
1.8. In order to mitigate the effects of changes added to this
implementation in version 1.9 of numpy, we include this implementation
here.
Estimates of the radial/axial diffusivities may rely on
empirical measurements (for example, the AD in the Corpus Callosum), or
may be based on a biophysical model of some kind.
The first dict in the tuple answers the question: Given a voxel (from
the unique indices in this model), which fibers pass through it?
The second answers the question: Given a streamline, for each voxel that
this streamline passes through, which nodes of that streamline are in that
voxel?
Identifies endpoints and invalid points to inform tracking.
seedsarray (N, 3), optional
Points to seed the tracking. Seed points should be given in point
space of the track (see affine).
affinearray (4, 4), optional
Coordinate space for the streamline point with respect to voxel
indices of input data. This affine can contain scaling, rotational,
and translational components but should not contain any shearing.
An identity matrix can be used to generate streamlines in “voxel
coordinates” as long as isotropic voxels were used to acquire the
data.
step_sizefloat, optional
Step size used for tracking.
max_crossint or None, optional
The maximum number of direction to track from each seed in crossing
voxels. By default all initial directions are tracked.
maxlenint, optional
Maximum length of generated streamlines. Longer streamlines will be
discarted if return_all=False.
minlenint, optional
Minimum length of generated streamlines. Shorter streamlines will
be discarted if return_all=False.
fixedstepbool, optional
If true, a fixed stepsize is used, otherwise a variable step size
is used.
return_allbool, optional
If true, return all generated streamlines, otherwise only
streamlines reaching end points or exiting the image.
random_seedint, optional
The seed for the random seed generator (numpy.random.seed and
random.seed).
save_seedsbool, optional
If True, return seeds alongside streamlines
unidirectionalbool, optional
If true, the tracking is performed only in the forward direction.
The seed position will be the first point of all streamlines.
randomize_forward_directionbool, optional
If true, the forward direction is randomized (multiplied by 1
or 1). Otherwise, the provided forward direction is used.
Initial direction to follow from the seed position. If
max_cross is None, one streamline will be generated per peak
per voxel. If None, direction_getter.initial_direction is used.
direction_getterinstance of ProbabilisticDirectionGetter
Used to get directions for fiber tracking.
stopping_criterioninstance of AnatomicalStoppingCriterion
Identifies endpoints and invalid points to inform tracking.
seedsarray (N, 3)
Points to seed the tracking. Seed points should be given in point
space of the track (see affine).
affinearray (4, 4)
Coordinate space for the streamline point with respect to voxel
indices of input data. This affine can contain scaling, rotational,
and translational components but should not contain any shearing.
An identity matrix can be used to generate streamlines in “voxel
coordinates” as long as isotropic voxels were used to acquire the
data.
step_sizefloat
Step size used for tracking.
max_crossint or None, optional
The maximum number of direction to track from each seed in crossing
voxels. By default all initial directions are tracked.
maxlenint, optional
Maximum length of generated streamlines. Longer streamlines will be
discarted if return_all=False.
minlenint, optional
Minimum length of generated streamlines. Shorter streamlines will
be discarted if return_all=False.
pft_back_tracking_distfloat
Distance in mm to back track before starting the particle filtering
tractography. The total particle filtering tractography distance is
equal to back_tracking_dist + front_tracking_dist.
By default this is set to 2 mm.
pft_front_tracking_distfloat, optional
Distance in mm to run the particle filtering tractography after the
the back track distance. The total particle filtering tractography
distance is equal to back_tracking_dist + front_tracking_dist. By
default this is set to 1 mm.
pft_max_trialint, optional
Maximum number of trial for the particle filtering tractography
(Prevents infinite loops).
particle_countint, optional
Number of particles to use in the particle filter.
return_allbool, optional
If true, return all generated streamlines, otherwise only
streamlines reaching end points or exiting the image.
random_seedint, optional
The seed for the random seed generator (numpy.random.seed and
random.seed).
save_seedsbool, optional
If True, return seeds alongside streamlines
min_wm_pve_before_stoppingint, optional
Minimum white matter pve (1  stopping_criterion.include_map 
stopping_criterion.exclude_map) to reach before allowing the
tractography to stop.
unidirectionalbool, optional
If true, the tracking is performed only in the forward direction.
The seed position will be the first point of all streamlines.
randomize_forward_directionbool, optional
If true, the forward direction is randomized (multiplied by 1
or 1). Otherwise, the provided forward direction is used.
Initial direction to follow from the seed position. If
max_cross is None, one streamline will be generated per peak
per voxel. If None, direction_getter.initial_direction is used.
Middle point of line, such that, if L is the line length then
np is the point such that the length xyz[0] to mp and from
mp to xyz[1] is L/2. If the middle point is not a point in
xyz, then we take the interpolation between the two nearest
xyz points. If xyz is empty, return a ValueError
downsample for a specific number of points along the streamline
Uses the length of the curve. It works in a similar fashion to
midpoint and arbitrarypoint but it also reduces the number of segments
of a streamline.
dipy.tracking.metrics.downsample is deprecated, Please use dipy.tracking.streamline.set_number_of_points instead
deprecated from version: 1.2
Raises <class ‘dipy.utils.deprecator.ExpiredDeprecationError’> as of version: 1.4
FrenetSerret Space Curve Invariants
Calculates the 3 vector and 2 scalar invariants of a space curve
defined by vectors r = (x,y,z). If z is omitted (i.e. the array xyz has
shape (N,2)), then the curve is
only 2D (planar), but the equations are still valid.
Similar to
https://www.mathworks.com/matlabcentral/fileexchange/11169frenet
In the following equations the prime (\('\)) indicates differentiation
with respect to the parameter \(s\) of a parametrised curve \(\mathbf{r}(s)\).
 \(\mathbf{T}=\mathbf{r'}/\mathbf{r'}\qquad\) (Tangent vector)}
 \(\mathbf{N}=\mathbf{T'}/\mathbf{T'}\qquad\) (Normal vector)
 \(\mathbf{B}=\mathbf{T}\times\mathbf{N}\qquad\) (Binormal vector)
 \(\kappa=\mathbf{T'}\qquad\) (Curvature)
 \(\mathrm{\tau}=\mathbf{B'}\cdot\mathbf{N}\) (Torsion)
Parameters
———
xyz : arraylike shape (N,3)
array representing x,y,z of N points in a track
Returns
——
T : array shape (N,3)
array representing the tangent of the curve xyz
N : array shape (N,3)
array representing the normal of the curve xyz
B : array shape (N,3)
array representing the binormal of the curve xyz
k : array shape (N,1)
array representing the curvature of the curve xyz
t : array shape (N,1)
array representing the torsion of the curve xyz
Examples
——–
Create a helix and calculate its tangent, normal, binormal, curvature
and torsion
>>> from dipy.tracking import metrics as tm
>>> import numpy as np
>>> theta = 2*np.pi*np.linspace(0,2,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=theta/(2*np.pi)
>>> xyz=np.vstack((x,y,z)).T
>>> T,N,B,k,t=tm.frenet_serret(xyz)
dipy.tracking.metrics.inside_sphere(xyz, center, radius)#
If any point of the track is inside a sphere of a specified
center and radius return True otherwise False. Mathematically this
can be simply described by \(xc\le r\) where \(x\) a point \(c\) the
center of the sphere and \(r\) the radius of the sphere.
Parameters
———
xyz : array, shape (N,3)
representing x,y,z of the N points of the track
center : array, shape (3,)
center of the sphere
radius : float
radius of the sphere
Returns
——
tf : {True,False}
Whether point is inside sphere.
Examples
——–
>>> from dipy.tracking.metrics import inside_sphere
>>> line=np.array(([0,0,0],[1,1,1],[2,2,2]))
>>> sph_cent=np.array([1,1,1])
>>> sph_radius = 1
>>> inside_sphere(line,sph_cent,sph_radius)
True
dipy.tracking.metrics.inside_sphere_points(xyz, center, radius)#
If a track intersects with a sphere of a specified center and
radius return the points that are inside the sphere otherwise False.
Mathematically this can be simply described by \(xc \le r\) where \(x\)
a point \(c\) the center of the sphere and \(r\) the radius of the
sphere.
Parameters
———
xyz : array, shape (N,3)
representing x,y,z of the N points of the track
center : array, shape (3,)
center of the sphere
radius : float
radius of the sphere
Returns
——
xyzn : array, shape(M,3)
array representing x,y,z of the M points inside the sphere
Examples
——–
>>> from dipy.tracking.metrics import inside_sphere_points
>>> line=np.array(([0,0,0],[1,1,1],[2,2,2]))
>>> sph_cent=np.array([1,1,1])
>>> sph_radius = 1
>>> inside_sphere_points(line,sph_cent,sph_radius)
array([[1, 1, 1]])
The scipy.interpolate packages wraps the netlib FITPACK routines
(Dierckx) for calculating smoothing splines for various kinds of
data and geometries. Although the data is evenly spaced in this
example, it need not be so to use this routine.
A smoothing condition. The amount of smoothness is determined by
satisfying the conditions: sum((w * (y  g))**2,axis=0) <= s
where g(x) is the smoothed interpolation of (x,y). The user can
use s to control the tradeoff between closeness and smoothness of
fit. Larger satisfying the conditions: sum((w * (y 
g))**2,axis=0) <= s where g(x) is the smoothed interpolation of
(x,y). The user can use s to control the tradeoff between
closeness and smoothness of fit. Larger s means more smoothing
while smaller values of s indicate less smoothing. Recommended
values of s depend on the weights, w. If the weights represent
the inverse of the standarddeviation of y, then a: good s value
should be found in the range (msqrt(2*m),m+sqrt(2*m)) where m is
the number of datapoints in x, y, and w.
kint, optional
Degree of the spline. Cubic splines are recommended. Even
values of k should be avoided especially with a small svalue.
for the same set of data. If task=1 find the weighted least
square spline for a given set of knots, t.
nestNone or int, optional
An overestimate of the total number of knots of the spline to
help in determining the storage space. None results in value
m+2*k. 1 results in m+k+1. Always large enough is nest=m+k+1.
Default is 1.
Arbitrary point of line, such that, if the arbitrary point is not
a point in xyz, then we take the interpolation between the two
nearest xyz points. If xyz is empty, return a ValueError
AnatomicallyConstrained Tractography (ACT) stopping criterion from [1]_.
This implements the use of partial volume fraction (PVE) maps to
determine when the tracking stops. The proposed ([1]_) method that
cuts streamlines going through subcortical gray matter regions is
not implemented here. The backtracking technique for
streamlines reaching INVALIDPOINT is not implemented either.
cdef:
Abstract class that takes as input included and excluded tissue maps.
The ‘include_map’ defines when the streamline reached a ‘valid’ stopping
region (e.g. gray matter partial volume estimation (PVE) map) and the
‘exclude_map’ defines when the streamline reached an ‘invalid’ stopping
region (e.g. corticospinal fluid PVE map). The background of the anatomical
image should be added to the ‘include_map’ to keep streamlines exiting the
brain (e.g. through the brain stem).
Continuous map criterion (CMC) stopping criterion from [1]_.
This implements the use of partial volume fraction (PVE) maps to
determine when the tracking stops.
Select streamlines based on logical relations with several regions of
interest (ROIs). For example, select streamlines that pass near ROI1,
but only if they do not pass near ROI2.
The mapping from voxel coordinates to streamline points.
The voxel_to_rasmm matrix, typically from a NIFTI file.
roislist or ndarray
A list of 3D arrays, each with shape (x, y, z) corresponding to the
shape of the brain volume, or a 4D array with shape (n_rois, x, y,
z). Nonzeros in each volume are considered to be within the region
includearray or list
A list or 1D array of boolean values marking inclusion or exclusion
criteria. If a streamline is near any of the inclusion ROIs, it
should evaluate to True, unless it is also near any of the exclusion
ROIs.
modestring, optional
One of {“any”, “all”, “either_end”, “both_end”}, where a streamline is
associated with an ROI if:
“any” : any point is within tol from ROI. Default.
“all” : all points are within tol from ROI.
“either_end” : either of the endpoints is within tol from ROI
“both_end” : both end points are within tol from ROI.
tolfloat
Distance (in the units of the streamlines, usually mm). If any
coordinate in the streamline is within this distance from the center
of any voxel in the ROI, the filtering criterion is set to True for
this streamline, otherwise False. Defaults to the distance between
the center of each voxel and the corner of the voxel.
Computes the cluster confidence index (cci), which is an
estimation of the support a set of streamlines gives to
a particular pathway.
Ex: A single streamline with no others in the dataset
following a similar pathway has a low cci. A streamline
in a bundle of 100 streamlines that follow similar
pathways has a high cci.
See: Jordan et al. 2017
(Based on streamline MDF distance from Garyfallidis et al. 2012)
A sequence of streamlines of length N (# streamlines)
max_mdfint
The maximum MDF distance (mm) that will be considered a
“supporting” streamline and included in cci calculation
subsample: int
The number of points that are considered for each streamline
in the calculation. To save on calculation time, each
streamline is subsampled to subsampleN points.
power: int
The power to which the MDF distance for each streamline
will be raised to determine how much it contributes to
the cci. High values of power make the contribution value
degrade much faster. E.g., a streamline with 5mm MDF
similarity contributes 1/5 to the cci if power is 1, but
only contributes 1/5^2 = 1/25 if power is 2.
override: bool, False by default
override means that the cci calculation will still occur even
though there are short streamlines in the dataset that may alter
expected behaviour.
[Jordan17] Jordan K. Et al., Cluster Confidence Index: A StreamlineWise
Pathway Reproducibility Metric for DiffusionWeighted MRI Tractography,
Journal of Neuroimaging, vol 28, no 1, 2017.
[Garyfallidis12] Garyfallidis E. et al., QuickBundles a method for
tractography simplification, Frontiers in Neuroscience,
vol 6, no 175, 2012.
List or generator of 2d arrays of 3d coordinates. Each array contains
the xyz coordinates of a single streamline.
affinearray_like (4, 4)
The mapping from voxel coordinates to streamline points.
The voxel_to_rasmm matrix, typically from a NIFTI file.
roi1, roi2ndarray
Binary masks designating the location of the regions of interest, or
coordinate arrays (nby3 array with ROI coordinate in each row).
in_placebool
Whether to make the change inplace in the original list
(and return a reference to the list), or to make a copy of the list
and return this copy, with the relevant streamlines reoriented.
Default: False.
as_generatorbool
Whether to return a generator as output. Default: False
This provides the standard orientation according to which the
streamlines in the provided bundle should be reoriented.
n_points: int, optional
The number of samples to apply to each of the streamlines.
in_placebool
Whether to make the change inplace in the original input
(and return a reference), or to make a copy of the list
and return this copy, with the relevant streamlines reoriented.
Default: False.
as_generatorbool
Whether to return a generator as output. Default: False
Values are extracted from the image based on the 3D coordinates of the
nodes that comprise the points in the streamline, without any interpolation
into segments between the nodes. Using this function with streamlines that
have been resampled into a very small number of nodes will result in very
few values.
Compress streamlines by linearization as in [Presseau15].
The compression consists in merging consecutive segments that are
nearly collinear. The merging is achieved by removing the point the two
segments have in common.
The linearization process [Presseau15] ensures that every point being
removed are within a certain margin (in mm) of the resulting streamline.
Recommendations for setting this margin can be found in [Presseau15]
(in which they called it tolerance error).
The compression also ensures that two consecutive points won’t be too far
from each other (precisely less or equal than `max_segment_length`mm).
This is a tradeoff to speed up the linearization process [Rheault15]. A low
value will result in a faster linearization but low compression, whereas
a high value will result in a slower linearization but high compression.
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.
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.
Change the number of points of streamlines in order to obtain
nb_points1 segments of equal length. Points of streamlines will be
modified along the curve.
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.
A streamline can pass through a voxel even if one of the points of the
streamline does not lie in the voxel. For example a step from [0,0,0] to
[0,0,2] passes through [0,0,1]. Consider subsegmenting the streamlines when
the edges of the voxels are smaller than the steps of the streamlines.
The number of connection between each pair of regions in
label_volume.
mappingdefaultdict(list)
mapping[i,j] returns all the streamlines that connect region i
to region j. If symmetric is True mapping will only have one key
for each start end pair such that if i<j mapping will have key
(i,j) but not key (j,i).
Split the segments of the streamlines into small segments.
Replaces each segment of each of the streamlines with the smallest possible
number of equally sized smaller segments such that no segment is longer
than max_segment_length. Among other things, this can useful for getting
streamline counts on a grid that is smaller than the length of the
streamline segments.
A binary array specifying where to place the seeds for fiber tracking.
affinearray, (4, 4)
The mapping between voxel indices and the point space for seeds.
The voxel_to_rasmm matrix, typically from a NIFTI file.
A seed point at the center the voxel [i,j,k]
will be represented as [x,y,z] where
[x,y,z,1]==np.dot(affine,[i,j,k,1]).
densityint or array_like (3,)
Specifies the number of seeds to place along each dimension. A
density of 2 is the same as [2,2,2] and will result in a
total of 8 seeds per voxel.
Create randomly placed seeds for fiber tracking from a binary mask.
Seeds points are placed randomly distributed in voxels of mask
which are True.
If seed_count_per_voxel is True, this function is
similar to seeds_from_mask(), with the difference that instead of
evenly distributing the seeds, it randomly places the seeds within the
voxels specified by the mask.
A binary array specifying where to place the seeds for fiber tracking.
affinearray, (4, 4)
The mapping between voxel indices and the point space for seeds.
The voxel_to_rasmm matrix, typically from a NIFTI file.
A seed point at the center the voxel [i,j,k]
will be represented as [x,y,z] where
[x,y,z,1]==np.dot(affine,[i,j,k,1]).
seeds_countint
The number of seeds to generate. If seed_count_per_voxel is True,
specifies the number of seeds to place in each voxel. Otherwise,
specifies the total number of seeds to place in the mask.
seed_count_per_voxel: bool
If True, seeds_count is per voxel, else seeds_count is the total number
of seeds.
random_seedint
The seed for the random seed generator (numpy.random.Generator).
Filter streamlines based on whether or not they pass through a ROI,
using a linebased algorithm. Mostly used as a replacement of target
for compressed streamlines.
This function never returns singlepoint streamlines, whatever the
value of include.
ROI coordinates transformed to the streamline coordinate frame.
tolfloat
Distance (in the units of the streamlines, usually mm). If any
coordinate in the streamline is within this distance from the center
of any voxel in the ROI, this function returns True.
modestring
One of {“any”, “all”, “either_end”, “both_end”}, where return True
if:
“any” : any point is within tol from ROI.
“all” : all points are within tol from ROI.
“either_end” : either of the endpoints is within tol from ROI
“both_end” : both end points are within tol from ROI.
A sequence of streamlines. Each streamline should be a (N, 3) array,
where N is the length of the streamline.
affinearray (4, 4)
The mapping between voxel indices and the point space for seeds.
The voxel_to_rasmm matrix, typically from a NIFTI file.
region_of_interestndarray
A mask used as a target. Nonzero values are considered to be within
the target region.
tolfloat
Distance (in the units of the streamlines, usually mm). If any
coordinate in the streamline is within this distance from the center
of any voxel in the ROI, the filtering criterion is set to True for
this streamline, otherwise False. Defaults to the distance between
the center of each voxel and the corner of the voxel.
modestring, optional
One of {“any”, “all”, “either_end”, “both_end”}, where return True
if:
“any” : any point is within tol from ROI. Default.
“all” : all points are within tol from ROI.
“either_end” : either of the endpoints is within tol from ROI
“both_end” : both end points are within tol from ROI.
A list of 3D arrays, each with shape (x, y, z) corresponding to the
shape of the brain volume, or a 4D array with shape (n_rois, x, y,
z). Nonzeros in each volume are considered to be within the region.
includearray or list
A list or 1D array of boolean marking inclusion or exclusion
criteria.
The include_roi and exclude_roi can be used to perform the operation: “(A
or B or …) and not (X or Y or …)”, where A, B are inclusion regions
and X, Y are exclusion regions.
The mapping from voxel coordinates to streamline coordinates.
The streamline values are assumed to be in voxel coordinates.
IE [0,0,0] is the center of the first voxel and the voxel size
is [1,1,1].
mapping_as_streamlinesbool, optional, False by default
If True voxel indices map to lists of streamline objects. Otherwise
voxel indices map to lists of integers.
Counts of points in tracks that pass through voxels in volume
We find whether a point passed through a track by rounding the mm
point values to voxels. For a track that passes through a voxel more
than once, we only record counts and elements for the first point in
the line that enters the voxel.
sequence of T tracks. One track is an ndarray of shape (N, 3), where N
is the number of points in that track, and tracks[t][n] is the nth
point in the tth track. Points are of form x, y, z in voxel mm
coordinates. That is, if i,j,k is the possibly noninteger voxel
coordinate of the track point, and vox_sizes are 3 floats giving voxel
sizes of dimensions 0, 1, 2 respectively, then the voxel mm coordinates
x,y,z are simply i*vox_sizes[0],j*vox_sizes[1],k*vox_sizes[2]. This convention derives from trackviz. To pass in
tracks as voxel coordinates, just pass vox_sizes=(1,1,1) (see
below).
vol_dimssequence length 3
volume dimensions in voxels, x, y, z.
vox_sizesoptional, sequence length 3
voxel sizes in mm. Default is (1,1,1)
return_elements{True, False}, optional
If True, also return object array with one list per voxel giving
track indices and point indices passing through the voxel (see
below)
An array where entry tcs[x,y,z] is the number of tracks
that passed through voxel at voxel coordinate x, y, z
tesndarray dtype np.object, shape vol_dim
If return_elements is True, we also return an object array with
one object per voxel. The objects at each voxel are a list of
integers, where the integers are the indices of the track that
passed through the voxel.
Imagine you have a volume (voxel) space of dimension (10,20,30).
Imagine you had voxel coordinate tracks in vs. To just fill an array
with the counts of how many tracks pass through each voxel: