tracking#

Tracking objects

Module: tracking._utils#

This is a helper module for dipy.tracking.utils.

Module: tracking.direction_getter#

Module: tracking.distances#

Optimized track distances, similarities and distanch clustering algorithms

add_3vecs(vec1, vec2)

approx_polygon_track(xyz[, alpha])

Fast and simple trajectory approximation algorithm by Eleftherios and Ian

approximate_mdl_trajectory(xyz[, alpha])

Implementation of Lee et al Approximate Trajectory Partitioning Algorithm

bundles_distances_mam(tracksA, tracksB[, metric])

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

bundles_distances_mdf(tracksA, tracksB)

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

cut_plane(tracks, ref)

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

inner_3vecs(vec1, vec2)

intersect_segment_cylinder(sa, sb, p, q, r)

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

larch_3merge(C[, thr])

Reassign tracks to existing clusters by merging clusters that their representative tracks are not very distant i.e. less than sqd_thr.

larch_3split(tracks[, indices, thr])

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

lee_angle_distance(start0, end0, start1, end1)

Calculates angle distance metric for the distance between two line segments

lee_perpendicular_distance(start0, end0, ...)

Calculates perpendicular distance metric for the distance between two line segments

local_skeleton_clustering(tracks[, d_thr])

Efficient tractography clustering

local_skeleton_clustering_3pts(tracks[, d_thr])

Does a first pass clustering

mam_distances(xyz1, xyz2[, metric])

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

minimum_closest_distance(xyz1, xyz2)

Find the minimum distance between two curves xyz1, xyz2

most_similar_track_mam(tracks[, metric])

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

mul_3vec(a, vec)

mul_3vecs(vec1, vec2)

norm_3vec(vec)

Euclidean (L2) norm of length 3 vector

normalized_3vec(vec)

Return normalized 3D vector

point_segment_sq_distance(a, b, c)

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

point_track_sq_distance_check(track, point, ...)

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

sub_3vecs(vec1, vec2)

track_dist_3pts(tracka, trackb)

Calculate the euclidean distance between two 3pt tracks

track_roi_intersection_check(track, roi, ...)

Check if a track is intersecting a region of interest

warn()

Issue a warning, or maybe ignore it or raise an exception.

Module: tracking.fbcmeasures#

FBCMeasures

KDTree(data[, leafsize, compact_nodes, ...])

kd-tree for quick nearest-neighbor lookup.

interp1d(x, y[, kind, axis, copy, ...])

Interpolate a 1-D function.

compute_rfbc(streamlines_length, ...[, ...])

Compute the relative fiber to bundle coherence (RFBC)

determine_num_threads(num_threads)

Determine the effective number of threads to be used for OpenMP calls

min_moving_average(a, n)

Return the lowest cumulative sum for the score of a streamline segment

Module: tracking.learning#

Learning algorithms for tractography

detect_corresponding_tracks(indices, ...)

Detect corresponding tracks from list tracks1 to list tracks2 where tracks1 & tracks2 are lists of tracks

detect_corresponding_tracks_plus(indices, ...)

Detect corresponding tracks from 1 to 2 where tracks1 & tracks2 are sequences of tracks

Module: tracking.life#

This is an implementation of the Linear Fascicle Evaluation (LiFE) algorithm described in [1].

References#

LifeSignalMaker(gtab, *[, evals, sphere])

A class for generating signals from streamlines in an efficient and speedy manner.

FiberModel(gtab)

A class for representing and solving predictive models based on tractography solutions.

FiberFit(fiber_model, life_matrix, ...)

A fit of the LiFE model to diffusion data

gradient(f)

Return the gradient of an N-dimensional array.

streamline_gradients(streamline)

Calculate the gradients of the streamline along the spatial dimension

grad_tensor(grad, evals)

Calculate the 3 by 3 tensor for a given spatial gradient, given a canonical tensor shape (also as a 3 by 3), pointing at [1,0,0]

streamline_tensors(streamline, *[, evals])

The tensors generated by this fiber.

streamline_signal(streamline, gtab, *[, evals])

The signal from a single streamline estimate along each of its nodes.

voxel2streamline(streamline, affine, *[, ...])

Maps voxels to streamlines and streamlines to voxels, for setting up the LiFE equations matrix

Module: tracking.local_tracking#

LocalTracking(direction_getter, ...[, ...])

ParticleFilteringTracking(direction_getter, ...)

Module: tracking.localtrack#

local_tracker(dg, sc, seed_pos, first_step, ...)

Tracks one direction from a seed.

pft_tracker(dg, sc, seed_pos, first_step, ...)

Tracks one direction from a seed using the particle filtering algorithm.

random(/)

Module: tracking.mesh#

random_coordinates_from_surface(...[, ...])

Generate random triangles_indices and trilinear_coord

seeds_from_surface_coordinates(triangles, ...)

Compute points from triangles_indices and trilinear_coord

triangles_area(triangles, vts)

Compute the local area of each triangle

vertices_to_triangles_values(triangles, ...)

Change from values per vertex to values per triangle

Module: tracking.metrics#

Metrics for tracks, where tracks are arrays of points

winding(xyz)

Total turning angle projected.

length(xyz, *[, along])

Euclidean length of track line

bytes(xyz)

Size of track in bytes.

midpoint(xyz)

Midpoint of track

center_of_mass(xyz)

Center of mass of streamline

magn(xyz, *[, n])

magnitude of vector

frenet_serret(xyz)

Frenet-Serret Space Curve Invariants

mean_curvature(xyz)

Calculates the mean curvature of a curve

mean_orientation(xyz)

Calculates the mean orientation of a curve

generate_combinations(items, n)

Combine sets of size n from items

longest_track_bundle(bundle, *[, sort])

Return longest track or length sorted track indices in bundle

intersect_sphere(xyz, center, radius)

If any segment of the track is intersecting with a sphere of specific center and radius return True otherwise False

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.

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.

spline(xyz, *[, s, k, nest])

Generate B-splines as documented in https://scipy-cookbook.readthedocs.io/items/Interpolation.html

startpoint(xyz)

First point of the track

endpoint(xyz)

arbitrarypoint(xyz, distance)

Select an arbitrary point along distance on the track (curve)

principal_components(xyz)

We use PCA to calculate the 3 principal directions for a track

midpoint2point(xyz, p)

Calculate distance from midpoint of a curve to arbitrary point p

Module: tracking.propspeed#

Track propagation performance functions

eudx_both_directions(seed, ref, qa, ind, ...)

ndarray_offset(indices, strides, lenind, ...)

Find offset in an N-dimensional ndarray using strides

Module: tracking.stopping_criterion#

ActStoppingCriterion

Anatomically-Constrained Tractography (ACT) stopping criterion.

AnatomicalStoppingCriterion

Abstract class that takes as input included and excluded tissue maps.

BinaryStoppingCriterion

cdef:

CmcStoppingCriterion

Continuous map criterion (CMC) stopping criterion.

StoppingCriterion

StreamlineStatus(value[, names, module, ...])

ThresholdStoppingCriterion

Module: tracking.streamline#

unlist_streamlines(streamlines)

Return the streamlines not as a list but as an array and an offset

relist_streamlines(points, offsets)

Given a representation of a set of streamlines as a large array and an offsets array return the streamlines as a list of shorter arrays.

center_streamlines(streamlines)

Move streamlines to the origin

deform_streamlines(streamlines, ...)

Apply deformation field to streamlines

transform_streamlines(streamlines, mat, *[, ...])

Apply affine transformation to streamlines

select_random_set_of_streamlines(...[, rng])

Select a random set of streamlines

select_by_rois(streamlines, affine, rois, ...)

Select streamlines based on logical relations with several regions of interest (ROIs).

cluster_confidence(streamlines, *[, ...])

Computes the cluster confidence index (cci), which is an estimation of the support a set of streamlines gives to a particular pathway.

orient_by_rois(streamlines, affine, roi1, ...)

Orient a set of streamlines according to a pair of ROIs

orient_by_streamline(streamlines, standard, *)

Orient a bundle of streamlines to a standard streamline.

values_from_volume(data, streamlines, affine)

Extract values of a scalar/vector along each streamline from a volume.

nbytes(streamlines)

Module: tracking.streamlinespeed#

Streamlines

alias of ArraySequence

as_native_array(arr)

Return arr as native byteordered array

compress_streamlines(streamlines[, ...])

Compress streamlines by linearization.

length(streamlines)

Euclidean length of streamlines

set_number_of_points(streamlines[, nb_points])

Change the number of points of streamlines

Module: tracking.utils#

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.

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.]

density_map(streamlines, affine, vol_dims)

Count the number of unique streamlines that pass through each voxel.

connectivity_matrix(streamlines, affine, ...)

Count the streamlines that start and end at each label pair.

ndbincount(x, *[, weights, shape])

Like bincount, but for nd-indices.

reduce_labels(label_volume)

Reduce an array of labels to the integers from 0 to n with smallest possible n.

subsegment(streamlines, max_segment_length)

Split the segments of the streamlines into small segments.

seeds_from_mask(mask, affine, *[, density])

Create seeds for fiber tracking from a binary mask.

random_seeds_from_mask(mask, affine, *[, ...])

Create randomly placed seeds for fiber tracking from a binary mask.

target(streamlines, affine, target_mask, *)

Filter streamlines based on whether or not they pass through an ROI.

target_line_based(streamlines, affine, ...)

Filter streamlines based on whether or not they pass through a ROI, using a line-based algorithm.

streamline_near_roi(streamline, roi_coords, ...)

Is a streamline near an ROI.

near_roi(streamlines, affine, ...[, tol, mode])

Provide filtering criteria for a set of streamlines based on whether they fall within a tolerance distance from an ROI.

length(streamlines)

Calculate the lengths of many streamlines in a bundle.

unique_rows(in_array, *[, dtype])

Find the unique rows in an array.

transform_tracking_output(tracking_output, ...)

Apply a linear transformation, given by affine, to streamlines.

reduce_rois(rois, include)

Reduce multiple ROIs to one inclusion and one exclusion ROI.

path_length(streamlines, affine, aoi, *[, ...])

Compute the shortest path, along any streamline, between aoi and each voxel.

max_angle_from_curvature(...)

Get the maximum deviation angle from the minimum radius curvature.

min_radius_curvature_from_angle(max_angle, ...)

Get minimum radius of curvature from a deviation angle.

seeds_directions_pairs(positions, peaks, *)

Pair each seed to the corresponding peaks.

Module: tracking.vox2track#

This module contains the parts of dipy.tracking.utils that need to be implemented in cython.

streamline_mapping(streamlines[, affine, ...])

Creates a mapping from voxel indices to streamlines.

track_counts(tracks, vol_dims[, vox_sizes, ...])

Counts of points in tracks that pass through voxels in volume

DirectionGetter#

class dipy.tracking.direction_getter.DirectionGetter#

Bases: object

Methods

generate_streamline

get_direction

initial_direction

generate_streamline(seed, direction, voxel_size, step_size, stopping_criterion, streamline, stream_status, fixedstep)#
get_direction(point, direction)#
initial_direction(point)#

add_3vecs#

dipy.tracking.distances.add_3vecs(vec1, vec2)#

approx_polygon_track#

dipy.tracking.distances.approx_polygon_track(xyz, alpha=0.392)#

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(xyz, alpha=1.0)#

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(tracksA, tracksB, metric='avg')#

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(tracksA, tracksB)#

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(tracks, ref)#

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(vec1, vec2)#

intersect_segment_cylinder#

dipy.tracking.distances.intersect_segment_cylinder(sa, sb, p, q, r)#

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(C, thr=10.0)#

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(tracks, indices=None, thr=10.0)#

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(start0, end0, start1, end1)#

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(start0, end0, start1, end1)#

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(tracks, d_thr=10)#

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(tracks, d_thr=10)#

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(xyz1, xyz2, metric='all')#

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

Based on the metrics in Zhang, Correia, Laidlaw 2008 https://ieeexplore.ieee.org/document/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(xyz1, xyz2)#

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(tracks, metric='avg')#

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(a, vec)#

mul_3vecs#

dipy.tracking.distances.mul_3vecs(vec1, vec2)#

norm_3vec#

dipy.tracking.distances.norm_3vec(vec)#

Euclidean (L2) norm of length 3 vector

Parameters:
vecarray-like shape (3,)
Returns:
normfloat

Euclidean norm

normalized_3vec#

dipy.tracking.distances.normalized_3vec(vec)#

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(a, b, c)#

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(track, point, sq_dist_thr)#

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(vec1, vec2)#

track_dist_3pts#

dipy.tracking.distances.track_dist_3pts(tracka, trackb)#

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(track, roi, sq_dist_thr)#

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

warn#

warn(/, message, category=None, stacklevel=1, source=None, *,
warn     skip_file_prefixes=<unrepresentable>)

Issue a warning, or maybe ignore it or raise an exception.

message

Text of the warning message.

category

The Warning category subclass. Defaults to UserWarning.

stacklevel

How far up the call stack to make this warning appear. A value of 2 for example attributes the warning to the caller of the code calling warn().

source

If supplied, the destroyed object which emitted a ResourceWarning

skip_file_prefixes

An optional tuple of module filename prefixes indicating frames to skip during stacklevel computations for stack frame attribution.

FBCMeasures#

class dipy.tracking.fbcmeasures.FBCMeasures#

Bases: object

Methods

get_points_rfbc_thresholded(threshold[, ...])

Set a threshold on the RFBC to remove spurious fibers.

get_points_rfbc_thresholded(threshold, emphasis=0.5, verbose=False)#

Set a threshold on the RFBC to remove spurious fibers.

Parameters:
thresholdfloat

The threshold to set on the RFBC, should be within 0 and 1.

emphasisfloat

Enhances the coloring of the fibers by LFBC. Increasing emphasis will stress spurious fibers by logarithmic weighting.

verboseboolean

Prints info about the found RFBC for the set of fibers such as median, mean, min and max values.

Returns:
outputtuple with 3 elements

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)

KDTree#

class dipy.tracking.fbcmeasures.KDTree(data, leafsize=10, compact_nodes=True, copy_data=False, balanced_tree=True, boxsize=None)[source]#

Bases: cKDTree

kd-tree for quick nearest-neighbor lookup.

This class provides an index into a set of k-dimensional points which can be used to rapidly look up the nearest neighbors of any point.

Parameters:
dataarray_like, shape (n,m)

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 kd-tree is built with copy_data=True.

leafsizepositive int, optional

The number of points at which the algorithm switches over to brute-force. Default: 10.

compact_nodesbool, optional

If True, the kd-tree 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 kd-tree 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 m-d 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 i-th 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.

Attributes:
datandarray, shape (n,m)

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 kd-tree is built with copy_data=True.

leafsizepositive int

The number of points at which the algorithm switches over to brute-force.

mint

The dimension of a single data-point.

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.

sizeint

The number of nodes in the tree.

Methods

count_neighbors(other, r[, p, weights, ...])

Count how many nearby pairs can be formed.

query(x[, k, eps, p, distance_upper_bound, ...])

Query the kd-tree for nearest neighbors.

query_ball_point(x, r[, p, eps, workers, ...])

Find all points within distance r of point(s) x.

query_ball_tree(other, r[, p, eps])

Find all pairs of points between self and other whose distance is at most r.

query_pairs(r[, p, eps, output_type])

Find all pairs of points in self whose distance is at most r.

sparse_distance_matrix(other, max_distance)

Compute a sparse distance matrix.

innernode

leafnode

node

Notes

The algorithm used is described in Maneewongvatana and Mount 1999. The general idea is that the kd-tree is a binary tree, each of whose nodes represents an axis-aligned 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. High-dimensional nearest-neighbor queries are a substantial open problem in computer science.

count_neighbors(other, r, p=2.0, weights=None, cumulative=True)[source]#

Count how many nearby pairs can be formed.

Count the number of pairs (x1,x2) can be formed, with x1 drawn from self and x2 drawn from other, and where distance(x1, x2, p) <= r.

Data points on self and other are optionally weighted by the weights argument. (See below)

This is adapted from the “two-point correlation” algorithm described by Gray and Moore [1]. See notes for further discussion.

Parameters:
otherKDTree

The other tree to draw points from, can be the same tree as self.

rfloat or one-dimensional array of floats

The radius to produce a count for. Multiple radii are searched with a single tree traversal. If the count is non-cumulative(cumulative=False), r defines the edges of the bins, and must be non-decreasing.

pfloat, optional

1<=p<=infinity. Which Minkowski p-norm 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 pair-counting 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

Added 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

Added in version 1.6.0.

Returns:
resultscalar or 1-D array

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 (-inf if i == 0 else r[i-1]) < R <= r[i]

Notes

Pair-counting 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 self-similar assembly of objects also occur.

The Landy-Szalay 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

  • objects R (random) that contains no signal,

\[\xi(r) = \frac{<D, D> - 2 f <D, R> + f^2<R, R>}{f^2<R, R>},\]

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 dual-tree algorithm described in [1]. We switch between two different pair-cumulation 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 pair-counting, weighted pair-counting counts the product of weights instead of number of pairs. Weighted pair-counting 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).

[1] (1,2)

Gray and Moore, “N-body problems in statistical learning”, Mining the sky, 2000, https://arxiv.org/abs/astro-ph/0012333

[2]

Landy and Szalay, “Bias and variance of angular correlation functions”, The Astrophysical Journal, 1993, http://adsabs.harvard.edu/abs/1993ApJ…412…64L

[3]

Sheth, Connolly and Skibba, “Marked correlations in galaxy formation models”, Arxiv e-print, 2005, https://arxiv.org/abs/astro-ph/0511773

[4]

Hawkins, et al., “The 2dF Galaxy Redshift Survey: correlation functions, peculiar velocities and the matter density of the Universe”, Monthly Notices of the Royal Astronomical Society, 2002, http://adsabs.harvard.edu/abs/2003MNRAS.346…78H

Examples

You can count neighbors number between two kd-trees within a distance:

>>> import numpy as np
>>> from scipy.spatial import KDTree
>>> rng = np.random.default_rng()
>>> points1 = rng.random((5, 2))
>>> points2 = rng.random((5, 2))
>>> kd_tree1 = KDTree(points1)
>>> kd_tree2 = KDTree(points2)
>>> kd_tree1.count_neighbors(kd_tree2, 0.2)
1

This number is same as the total pair number calculated by query_ball_tree:

>>> indexes = kd_tree1.query_ball_tree(kd_tree2, r=0.2)
>>> sum([len(i) for i in indexes])
1
class innernode(ckdtreenode)[source]#

Bases: node

Attributes:
children
split
split_dim
property children#
property split#
property split_dim#
class leafnode(ckdtree_node=None)[source]#

Bases: node

Attributes:
children
idx
property children#
property idx#
class node(ckdtree_node=None)[source]#

Bases: object

query(x, k=1, eps=0, p=2, distance_upper_bound=inf, workers=1)[source]#

Query the kd-tree for nearest neighbors.

Parameters:
xarray_like, last dimension self.m

An array of points to query.

kint or Sequence[int], optional

Either the number of nearest neighbors to return, or a list of the k-th 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 p-norm to use. 1 is the sum-of-absolute-values distance (“Manhattan” distance). 2 is the usual Euclidean distance. infinity is the maximum-coordinate-difference 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 nearest-neighbor 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.

Added in version 1.6.0.

Returns:
dfloat or array of floats

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.

Examples

>>> import numpy as np
>>> from scipy.spatial import KDTree
>>> x, y = np.mgrid[0:5, 2:8]
>>> tree = KDTree(np.c_[x.ravel(), y.ravel()])

To query the nearest neighbours and return squeezed result, use

>>> dd, ii = tree.query([[0, 0], [2.2, 2.9]], k=1)
>>> print(dd, ii, sep='\n')
[2.         0.2236068]
[ 0 13]

To query the nearest neighbours and return unsqueezed result, use

>>> dd, ii = tree.query([[0, 0], [2.2, 2.9]], k=[1])
>>> print(dd, ii, sep='\n')
[[2.        ]
 [0.2236068]]
[[ 0]
 [13]]

To query the second nearest neighbours and return unsqueezed result, use

>>> dd, ii = tree.query([[0, 0], [2.2, 2.9]], k=[2])
>>> print(dd, ii, sep='\n')
[[2.23606798]
 [0.80622577]]
[[ 6]
 [19]]

To query the first and second nearest neighbours, use

>>> dd, ii = tree.query([[0, 0], [2.2, 2.9]], k=2)
>>> print(dd, ii, sep='\n')
[[2.         2.23606798]
 [0.2236068  0.80622577]]
[[ 0  6]
 [13 19]]

or, be more specific

>>> dd, ii = tree.query([[0, 0], [2.2, 2.9]], k=[1, 2])
>>> print(dd, ii, sep='\n')
[[2.         2.23606798]
 [0.2236068  0.80622577]]
[[ 0  6]
 [13 19]]
query_ball_point(x, r, p=2.0, eps=0, workers=1, return_sorted=None, return_length=False)[source]#

Find all points within distance r of point(s) x.

Parameters:
xarray_like, shape tuple + (self.m,)

The point or points to search for neighbors of.

rarray_like, float

The radius of points to return, must broadcast to the length of x.

pfloat, optional

Which Minkowski p-norm 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.

Added in version 1.6.0.

return_sortedbool, optional

Sorts returned indices if True and does not sort them if False. If None, does not sort single point queries, but does sort multi-point queries which was the behavior before this option was added.

Added in version 1.6.0.

return_lengthbool, optional

Return the number of points inside the radius instead of a list of the indices.

Added in version 1.6.0.

Returns:
resultslist or array of lists

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.

Notes

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.

Examples

>>> import numpy as np
>>> from scipy import spatial
>>> x, y = np.mgrid[0:5, 0:5]
>>> points = np.c_[x.ravel(), y.ravel()]
>>> tree = spatial.KDTree(points)
>>> sorted(tree.query_ball_point([2, 0], 1))
[5, 10, 11, 15]

Query multiple points and plot the results:

>>> import matplotlib.pyplot as plt
>>> points = np.asarray(points)
>>> plt.plot(points[:,0], points[:,1], '.')
>>> for results in tree.query_ball_point(([2, 0], [3, 3]), 1):
...     nearby_points = points[results]
...     plt.plot(nearby_points[:,0], nearby_points[:,1], 'o')
>>> plt.margins(0.1, 0.1)
>>> plt.show()
query_ball_tree(other, r, p=2.0, eps=0)[source]#

Find all pairs of points between self and other whose distance is at most r.

Parameters:
otherKDTree instance

The tree containing points to search against.

rfloat

The maximum distance, has to be positive.

pfloat, optional

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 non-negative.

Returns:
resultslist of lists

For each element self.data[i] of this tree, results[i] is a list of the indices of its neighbors in other.data.

Examples

You can search all pairs of points between two kd-trees within a distance:

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from scipy.spatial import KDTree
>>> rng = np.random.default_rng()
>>> points1 = rng.random((15, 2))
>>> points2 = rng.random((15, 2))
>>> plt.figure(figsize=(6, 6))
>>> plt.plot(points1[:, 0], points1[:, 1], "xk", markersize=14)
>>> plt.plot(points2[:, 0], points2[:, 1], "og", markersize=14)
>>> kd_tree1 = KDTree(points1)
>>> kd_tree2 = KDTree(points2)
>>> indexes = kd_tree1.query_ball_tree(kd_tree2, r=0.2)
>>> for i in range(len(indexes)):
...     for j in indexes[i]:
...         plt.plot([points1[i, 0], points2[j, 0]],
...             [points1[i, 1], points2[j, 1]], "-r")
>>> plt.show()
query_pairs(r, p=2.0, eps=0, output_type='set')[source]#

Find all pairs of points in self whose distance is at most r.

Parameters:
rpositive float

The maximum distance.

pfloat, optional

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 non-negative.

output_typestring, optional

Choose the output container, ‘set’ or ‘ndarray’. Default: ‘set’

Added in version 1.6.0.

Returns:
resultsset or ndarray

Set of pairs (i,j), with i < j, for which the corresponding positions are close. If output_type is ‘ndarray’, an ndarry is returned instead of a set.

Examples

You can search all pairs of points in a kd-tree within a distance:

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from scipy.spatial import KDTree
>>> rng = np.random.default_rng()
>>> points = rng.random((20, 2))
>>> plt.figure(figsize=(6, 6))
>>> plt.plot(points[:, 0], points[:, 1], "xk", markersize=14)
>>> kd_tree = KDTree(points)
>>> pairs = kd_tree.query_pairs(r=0.2)
>>> for (i, j) in pairs:
...     plt.plot([points[i, 0], points[j, 0]],
...             [points[i, 1], points[j, 1]], "-r")
>>> plt.show()
sparse_distance_matrix(other, max_distance, p=2.0, output_type='dok_matrix')[source]#

Compute a sparse distance matrix.

Computes a distance matrix between two KDTrees, leaving as zero any distance greater than max_distance.

Parameters:
otherKDTree
max_distancepositive float
pfloat, 1<=p<=infinity

Which Minkowski p-norm to use. A finite large p may cause a ValueError if overflow can occur.

output_typestring, optional

Which container to use for output data. Options: ‘dok_matrix’, ‘coo_matrix’, ‘dict’, or ‘ndarray’. Default: ‘dok_matrix’.

Added in version 1.6.0.

Returns:
resultdok_matrix, coo_matrix, dict or ndarray

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,

Examples

You can compute a sparse distance matrix between two kd-trees:

>>> import numpy as np
>>> from scipy.spatial import KDTree
>>> rng = np.random.default_rng()
>>> points1 = rng.random((5, 2))
>>> points2 = rng.random((5, 2))
>>> kd_tree1 = KDTree(points1)
>>> kd_tree2 = KDTree(points2)
>>> sdm = kd_tree1.sparse_distance_matrix(kd_tree2, 0.3)
>>> sdm.toarray()
array([[0.        , 0.        , 0.12295571, 0.        , 0.        ],
   [0.        , 0.        , 0.        , 0.        , 0.        ],
   [0.28942611, 0.        , 0.        , 0.2333084 , 0.        ],
   [0.        , 0.        , 0.        , 0.        , 0.        ],
   [0.24617575, 0.29571802, 0.26836782, 0.        , 0.        ]])

You can check distances above the max_distance are zeros:

>>> from scipy.spatial import distance_matrix
>>> distance_matrix(points1, points2)
array([[0.56906522, 0.39923701, 0.12295571, 0.8658745 , 0.79428925],
   [0.37327919, 0.7225693 , 0.87665969, 0.32580855, 0.75679479],
   [0.28942611, 0.30088013, 0.6395831 , 0.2333084 , 0.33630734],
   [0.31994999, 0.72658602, 0.71124834, 0.55396483, 0.90785663],
   [0.24617575, 0.29571802, 0.26836782, 0.57714465, 0.6473269 ]])
property tree#

interp1d#

class dipy.tracking.fbcmeasures.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)[source]#

Bases: _Interpolator1D

Interpolate a 1-D 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.

Parameters:
x(npoints, ) array_like

A 1-D array of real values.

y(…, npoints, …) array_like

A N-D array of real values. The length of y along the interpolation axis must be equal to the length of x. Use the axis parameter to select correct axis. Unlike other interpolators, the default interpolation axis is the last axis of y.

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’, ‘nearest-up’, ‘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; ‘nearest-up’ and ‘nearest’ differ when interpolating half-integers (e.g. 0.5, 1.5) in that ‘nearest-up’ rounds up and ‘nearest’ rounds down. Default is ‘linear’.

axisint, optional

Axis in the y array corresponding to the x-coordinate values. Unlike other interpolators, defaults to axis=-1.

copybool, optional

If True, the class makes internal copies of x and y. If False, references to x and y are used if possible. 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_valuearray-like or (array-like, 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 array-like must broadcast properly to the dimensions of the non-interpolation axes.

  • If a two-element 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 2-element tuple (e.g., list or ndarray, regardless of shape) is taken to be a single array-like argument meant to be used for both bounds as below, above = fill_value, fill_value. Using a two-element tuple or ndarray requires bounds_error=False.

    Added in version 0.17.0.

  • If “extrapolate”, then points outside the data range will be extrapolated.

    Added 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.

Attributes:
fill_value

The fill value.

Methods

__call__(x)

Evaluate the interpolant

See also

splrep, splev

Spline interpolation/smoothing based on FITPACK.

UnivariateSpline

An object-oriented wrapper of the FITPACK routines.

interp2d

2-D interpolation

Notes

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.

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate
>>> x = np.arange(0, 10)
>>> y = np.exp(-x/3.0)
>>> f = interpolate.interp1d(x, y)
>>> xnew = np.arange(0, 9, 0.1)
>>> ynew = f(xnew)   # use interpolation function returned by `interp1d`
>>> plt.plot(x, y, 'o', xnew, ynew, '-')
>>> plt.show()
dtype#
property fill_value#

The fill value.

compute_rfbc#

dipy.tracking.fbcmeasures.compute_rfbc(streamlines_length, streamline_scores, max_windowsize=7)#

Compute the relative fiber to bundle coherence (RFBC)

Parameters:
streamlines_length1D int array

Contains the length of each streamline

streamlines_scores2D double array

Contains the local fiber to bundle coherence (LFBC) for each streamline element.

max_windowsizeint

The maximal window size used to calculate the average LFBC region

Returns:
output: normalized lowest average LFBC region along the fiber

determine_num_threads#

dipy.tracking.fbcmeasures.determine_num_threads(num_threads)#

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.

  • For num_threads = 0 a ValueError is raised.

Parameters:
num_threadsint or None

Desired number of threads to be used.

min_moving_average#

dipy.tracking.fbcmeasures.min_moving_average(a, n)#

Return the lowest cumulative sum for the score of a streamline segment

Parameters:
aarray

Input array

nint

Length of the segment

Returns:
output: normalized lowest average LFBC region along the fiber

detect_corresponding_tracks#

dipy.tracking.learning.detect_corresponding_tracks(indices, tracks1, tracks2)[source]#

Detect corresponding tracks from list tracks1 to list tracks2 where tracks1 & tracks2 are lists of tracks

Parameters:
indicessequence

of indices of tracks1 that are to be detected in tracks2

tracks1sequence

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

tracks2sequence

of tracks as arrays, shape (M1,3) .. (Mm,3)

Returns:
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

Notes

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

Examples

>>> import numpy as np
>>> import dipy.tracking.learning as tl
>>> A = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]])
>>> B = np.array([[1, 0, 0], [2, 0, 0], [3, 0, 0]])
>>> C = np.array([[0, 0, -1], [0, 0, -2], [0, 0, -3]])
>>> bundle1 = [A, B, C]
>>> bundle2 = [B, A]
>>> indices = [0, 1]
>>> arr = tl.detect_corresponding_tracks(indices, bundle1, bundle2)

detect_corresponding_tracks_plus#

dipy.tracking.learning.detect_corresponding_tracks_plus(indices, tracks1, indices2, tracks2)[source]#

Detect corresponding tracks from 1 to 2 where tracks1 & tracks2 are sequences of tracks

Parameters:
indicessequence

of indices of tracks1 that are to be detected in tracks2

tracks1sequence

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

indices2sequence

of indices of tracks2 in the initial brain

tracks2sequence

of tracks as arrays, shape (M1,3) .. (Mm,3)

Returns:
track2trackarray (N,2) where N is len(indices)

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

See also

distances.mam_distances

Notes

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

Examples

>>> import numpy as np
>>> import dipy.tracking.learning as tl
>>> A = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]])
>>> B = np.array([[1, 0, 0], [2, 0, 0], [3, 0, 0]])
>>> C = np.array([[0, 0, -1], [0, 0, -2], [0, 0, -3]])
>>> bundle1 = [A, B, C]
>>> bundle2 = [B, A]
>>> indices = [0, 1]
>>> indices2 = indices
>>> arr = tl.detect_corresponding_tracks_plus(indices, bundle1, indices2, bundle2)

LifeSignalMaker#

class dipy.tracking.life.LifeSignalMaker(gtab, *, evals=(0.001, 0, 0), sphere=None)[source]#

Bases: object

A class for generating signals from streamlines in an efficient and speedy manner.

Methods

streamline_signal(streamline)

Approximate the signal for a given streamline

calc_signal

calc_signal(xyz)[source]#
streamline_signal(streamline)[source]#

Approximate the signal for a given streamline

FiberModel#

class dipy.tracking.life.FiberModel(gtab)[source]#

Bases: ReconstModel

A class for representing and solving predictive models based on tractography solutions.

Methods

fit(data, streamline, affine, *[, evals, sphere])

Fit the LiFE FiberModel for data and a set of streamlines associated with this data

setup(streamline, affine, *[, evals, sphere])

Set up the necessary components for the LiFE model: the matrix of fiber-contributions to the DWI signal, and the coordinates of voxels for which the equations will be solved

Notes

This is an implementation of the LiFE model described in [1].

References

fit(data, streamline, affine, *, evals=(0.001, 0, 0), sphere=None)[source]#

Fit the LiFE FiberModel for data and a set of streamlines associated with this data

Parameters:
data4D array

Diffusion-weighted data

streamlinelist

A bunch of streamlines

affinearray_like (4, 4)

The mapping from voxel coordinates to streamline points. The voxel_to_rasmm matrix, typically from a NIFTI file.

evalsarray-like, 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 speed-up 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.

Returns:
FiberFit class instance
setup(streamline, affine, *, evals=(0.001, 0, 0), sphere=None)[source]#

Set up the necessary components for the LiFE model: the matrix of fiber-contributions to the DWI signal, and the coordinates of voxels for which the equations will be solved

Parameters:
streamlinelist

Streamlines, each is an array of shape (n, 3)

affinearray_like (4, 4)

The mapping from voxel coordinates to streamline points. The voxel_to_rasmm matrix, typically from a NIFTI file.

evalsarray-like (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 speed-up 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 724-vertex symmetric sphere from dipy.data

FiberFit#

class dipy.tracking.life.FiberFit(fiber_model, life_matrix, vox_coords, to_fit, beta, weighted_signal, b0_signal, relative_signal, mean_sig, vox_data, streamline, affine, evals)[source]#

Bases: ReconstFit

A fit of the LiFE model to diffusion data

Methods

predict(*[, gtab, S0])

Predict the signal

predict(*, gtab=None, S0=None)[source]#

Predict the signal

Parameters:
gtabGradientTable

Default: use self.gtab

S0float or array

The non-diffusion-weighted signal in the voxels for which a prediction is made. Default: use self.b0_signal

Returns:
predictionndarray of shape (voxels, bvecs)

An array with a prediction of the signal in each voxel/direction

gradient#

dipy.tracking.life.gradient(f)[source]#

Return the gradient of an N-dimensional array.

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.

Parameters:
farray_like

An N-dimensional array containing samples of a scalar function.

Returns:
gradientndarray

N arrays of the same shape as f giving the derivative of f with respect to each dimension.

Notes

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.

Examples

>>> x = np.array([1, 2, 4, 7, 11, 16], dtype=float)
>>> gradient(x)
array([ 1. ,  1.5,  2.5,  3.5,  4.5,  5. ])
>>> gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
[array([[ 2.,  2., -1.],
       [ 2.,  2., -1.]]), array([[ 1. ,  2.5,  4. ],
       [ 1. ,  1. ,  1. ]])]

streamline_gradients#

dipy.tracking.life.streamline_gradients(streamline)[source]#

Calculate the gradients of the streamline along the spatial dimension

Parameters:
streamlinearray-like of shape (n, 3)

The 3d coordinates of a single streamline

Returns:
Array of shape (3, n): Spatial gradients along the length of the
streamline.

grad_tensor#

dipy.tracking.life.grad_tensor(grad, evals)[source]#

Calculate the 3 by 3 tensor for a given spatial gradient, given a canonical tensor shape (also as a 3 by 3), pointing at [1,0,0]

Parameters:
grad1d array of shape (3,)

The spatial gradient (e.g between two nodes of a streamline).

evals: 1d array of shape (3,)

The eigenvalues of a canonical tensor to be used as a response function.

streamline_tensors#

dipy.tracking.life.streamline_tensors(streamline, *, evals=(0.001, 0, 0))[source]#

The tensors generated by this fiber.

Parameters:
streamlinearray-like of shape (n, 3)

The 3d coordinates of a single streamline

evalsiterable with three entries, optional

The estimated eigenvalues of a single fiber tensor.

Returns:
An n_nodes by 3 by 3 array with the tensor for each node in the fiber.

Notes

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.

streamline_signal#

dipy.tracking.life.streamline_signal(streamline, gtab, *, evals=(0.001, 0, 0))[source]#

The signal from a single streamline estimate along each of its nodes.

Parameters:
streamlinea single streamline

Streamline data.

gtabGradientTable class instance

Gradient table.

evalsarray-like of length 3, optional

The eigenvalues of the canonical tensor used as an estimate of the signal generated by each node of the streamline.

voxel2streamline#

dipy.tracking.life.voxel2streamline(streamline, affine, *, unique_idx=None)[source]#

Maps voxels to streamlines and streamlines to voxels, for setting up the LiFE equations matrix

Parameters:
streamlinelist

A collection of streamlines, each n by 3, with n being the number of nodes in the fiber.

affinearray_like (4, 4)

The mapping from voxel coordinates to streamline points. The voxel_to_rasmm matrix, typically from a NIFTI file.

unique_idxarray, optional.

The unique indices in the streamlines

Returns:
v2f, v2fntuple of dicts
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?

LocalTracking#

class dipy.tracking.local_tracking.LocalTracking(direction_getter, stopping_criterion, seeds, affine, step_size, *, max_cross=None, maxlen=500, minlen=2, fixedstep=True, return_all=True, random_seed=None, save_seeds=False, unidirectional=False, randomize_forward_direction=False, initial_directions=None)[source]#

Bases: object

ParticleFilteringTracking#

class dipy.tracking.local_tracking.ParticleFilteringTracking(direction_getter, stopping_criterion, seeds, affine, step_size, *, max_cross=None, maxlen=500, minlen=2, pft_back_tracking_dist=2, pft_front_tracking_dist=1, pft_max_trial=20, particle_count=15, return_all=True, random_seed=None, save_seeds=False, min_wm_pve_before_stopping=0, unidirectional=False, randomize_forward_direction=False, initial_directions=None)[source]#

Bases: LocalTracking

local_tracker#

dipy.tracking.localtrack.local_tracker(dg, sc, seed_pos, first_step, voxel_size, streamline, step_size, fixedstep)#

Tracks one direction from a seed.

This function is the main workhorse of the LocalTracking class defined in dipy.tracking.local_tracking.

Parameters:
dgDirectionGetter

Used to choosing tracking directions.

scStoppingCriterion

Used to check the streamline status (e.g. endpoint) along path.

seed_posarray, float, 1d, (3,)

First point of the (partial) streamline.

first_steparray, float, 1d, (3,)

Initial seeding direction. Used as prev_dir for selecting the step direction from the seed point.

voxel_sizearray, float, 1d, (3,)

Size of voxels in the data set.

streamlinearray, float, 2d, (N, 3)

Output of tracking will be put into this array. The length of this array, N, will set the maximum allowable length of the streamline.

step_sizefloat

Size of tracking steps in mm if fixed_step.

fixedstepint

If greater than 0, a fixed step_size is used, otherwise a variable step size is used.

Returns:
endint

Length of the tracked streamline

stream_statusStreamlineStatus

Ending state of the streamlines as determined by the StoppingCriterion.

pft_tracker#

dipy.tracking.localtrack.pft_tracker(dg, sc, seed_pos, first_step, voxel_size, streamline, directions, step_size, pft_max_nbr_back_steps, pft_max_nbr_front_steps, pft_max_trials, particle_count, particle_paths, particle_dirs, particle_weights, particle_steps, particle_stream_statuses, min_wm_pve_before_stopping)#

Tracks one direction from a seed using the particle filtering algorithm.

This function is the main workhorse of the ParticleFilteringTracking class defined in dipy.tracking.local_tracking.

Parameters:
dgDirectionGetter

Used to choosing tracking directions.

scAnatomicalStoppingCriterion

Used to check the streamline status (e.g. endpoint) along path.

seed_posarray, float, 1d, (3,)

First point of the (partial) streamline.

first_steparray, float, 1d, (3,)

Initial seeding direction. Used as prev_dir for selecting the step direction from the seed point.

voxel_sizearray, float, 1d, (3,)

Size of voxels in the data set.

streamlinearray, float, 2d, (N, 3)

Output of tracking will be put into this array. The length of this array, N, will set the maximum allowable length of the streamline.

directionsarray, float, 2d, (N, 3)

Output of tracking directions will be put into this array. The length of this array, N, will set the maximum allowable length of the streamline.

step_sizefloat

Size of tracking steps in mm if fixed_step.

pft_max_nbr_back_stepsint

Number of tracking steps to back track before starting the particle filtering tractography.

pft_max_nbr_front_stepsint

Number of additional tracking steps to track.

pft_max_trialsint

Maximum number of trials for the particle filtering tractography (Prevents infinite loops).

particle_countint

Number of particles to use in the particle filter.

particle_pathsarray, float, 4d, (2, particle_count, pft_max_steps, 3)

Temporary array for paths followed by all particles.

particle_dirsarray, float, 4d, (2, particle_count, pft_max_steps, 3)

Temporary array for directions followed by particles.

particle_weightsarray, float, 1d (particle_count)

Temporary array for the weights of particles.

particle_stepsarray, float, (2, particle_count)

Temporary array for the number of steps of particles.

particle_stream_statusesarray, float, (2, particle_count)

Temporary array for the stream status of particles.

min_wm_pve_before_stoppingint, optional

Minimum white matter pve (1 - sc.include_map - sc.exclude_map) to reach before allowing the tractography to stop.

Returns:
endint

Length of the tracked streamline

stream_statusStreamlineStatus

Ending state of the streamlines as determined by the StoppingCriterion.

random#

dipy.tracking.localtrack.random(/)#

random_coordinates_from_surface#

dipy.tracking.mesh.random_coordinates_from_surface(nb_triangles, nb_seed, *, triangles_mask=None, triangles_weight=None, rand_gen=None)[source]#

Generate random triangles_indices and trilinear_coord

Triangles_indices probability are weighted by triangles_weight,

for each triangles inside the given triangles_mask

Parameters:
nb_trianglesint (n)

The amount of triangles in the mesh

nb_seedint

The number of random indices and coordinates generated.

triangles_mask[n] numpy array

Specifies which triangles should be chosen (or not)

triangles_weight[n] numpy array

Specifies the weight/probability of choosing each triangle

rand_genint

The seed for the random seed generator (numpy.random.seed).

Returns:
triangles_idx: [s] array

Randomly chosen triangles_indices

trilin_coord: [s,3] array

Randomly chosen trilinear_coordinates

See also

seeds_from_surface_coordinates, random_seeds_from_mask

seeds_from_surface_coordinates#

dipy.tracking.mesh.seeds_from_surface_coordinates(triangles, vts_values, triangles_idx, trilinear_coord)[source]#

Compute points from triangles_indices and trilinear_coord

Parameters:
triangles[n, 3] -> m array

A list of triangles from a mesh

vts_values[m, .] array

List of values to interpolates from coordinates along vertices, (vertices, vertices_normal, vertices_colors …)

triangles_idx[s] array

Specifies which triangles should be chosen (or not)

trilinear_coord[s, 3] array

Specifies the weight/probability of choosing each triangle

Returns:
pts[s, …] array

Interpolated values of vertices with triangles_idx and trilinear_coord

triangles_area#

dipy.tracking.mesh.triangles_area(triangles, vts)[source]#

Compute the local area of each triangle

Parameters:
triangles[n, 3] -> m array

A list of triangles from a mesh

vts[m, .] array

List of vertices

Returns:
triangles_area[m] array

List of area for each triangle in the mesh

vertices_to_triangles_values#

dipy.tracking.mesh.vertices_to_triangles_values(triangles, vts_values)[source]#

Change from values per vertex to values per triangle

Parameters:
triangles[n, 3] -> m array

A list of triangles from a mesh

vts_values[m, .] array

List of values to interpolates from coordinates along vertices, (vertices, vertices_normal, vertices_colors …)

Returns:
triangles_values[m] array

List of values for each triangle in the mesh

winding#

dipy.tracking.metrics.winding(xyz)[source]#

Total turning angle projected.

Project space curve to best fitting plane. Calculate the cumulative signed angle between each line segment and the previous one.

Parameters:
xyzarray-like shape (N,3)

Array representing x,y,z of N points in a track.

Returns:
ascalar

Total turning angle in degrees.

length#

dipy.tracking.metrics.length(xyz, *, along=False)[source]#

Euclidean length of track line

This will give length in mm if tracks are expressed in world coordinates.

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

alongbool, optional

If True, return array giving cumulative length along track, otherwise (default) return scalar giving total length.

Returns:
Lscalar or array shape (N-1,)

scalar in case of along == False, giving total length, array if along == True, giving cumulative lengths.

Examples

>>> from dipy.tracking.metrics import length
>>> xyz = np.array([[1,1,1],[2,3,4],[0,0,0]])
>>> expected_lens = np.sqrt([1+2**2+3**2, 2**2+3**2+4**2])
>>> length(xyz) == expected_lens.sum()
True
>>> len_along = length(xyz, along=True)
>>> np.allclose(len_along, expected_lens.cumsum())
True
>>> length([])
0
>>> length([[1, 2, 3]])
0
>>> length([], along=True)
array([0])

bytes#

dipy.tracking.metrics.bytes(xyz)[source]#

Size of track in bytes.

Parameters:
xyzarray-like shape (N,3)

Array representing x,y,z of N points in a track.

Returns:
bint

Number of bytes.

midpoint#

dipy.tracking.metrics.midpoint(xyz)[source]#

Midpoint of track

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

Returns:
mparray shape (3,)

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

Examples

>>> from dipy.tracking.metrics import midpoint
>>> midpoint([])
Traceback (most recent call last):
   ...
ValueError: xyz array cannot be empty
>>> midpoint([[1, 2, 3]])
array([1, 2, 3])
>>> xyz = np.array([[1,1,1],[2,3,4]])
>>> midpoint(xyz)
array([ 1.5,  2. ,  2.5])
>>> xyz = np.array([[0,0,0],[1,1,1],[2,2,2]])
>>> midpoint(xyz)
array([ 1.,  1.,  1.])
>>> xyz = np.array([[0,0,0],[1,0,0],[3,0,0]])
>>> midpoint(xyz)
array([ 1.5,  0. ,  0. ])
>>> xyz = np.array([[0,9,7],[1,9,7],[3,9,7]])
>>> midpoint(xyz)
array([ 1.5,  9. ,  7. ])

center_of_mass#

dipy.tracking.metrics.center_of_mass(xyz)[source]#

Center of mass of streamline

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

Returns:
comarray shape (3,)

center of mass of streamline

Examples

>>> from dipy.tracking.metrics import center_of_mass
>>> center_of_mass([])
Traceback (most recent call last):
   ...
ValueError: xyz array cannot be empty
>>> center_of_mass([[1,1,1]])
array([ 1.,  1.,  1.])
>>> xyz = np.array([[0,0,0],[1,1,1],[2,2,2]])
>>> center_of_mass(xyz)
array([ 1.,  1.,  1.])

magn#

dipy.tracking.metrics.magn(xyz, *, n=1)[source]#

magnitude of vector

frenet_serret#

dipy.tracking.metrics.frenet_serret(xyz)[source]#

Frenet-Serret 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/11169-frenet

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:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

Returns:
Tarray shape (N,3)

array representing the tangent of the curve xyz

Narray shape (N,3)

array representing the normal of the curve xyz

Barray shape (N,3)

array representing the binormal of the curve xyz

karray shape (N,1)

array representing the curvature of the curve xyz

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

mean_curvature#

dipy.tracking.metrics.mean_curvature(xyz)[source]#

Calculates the mean curvature of a curve

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a curve

Returns:
mfloat

Mean curvature.

Examples

Create a straight line and a semi-circle and print their mean curvatures

>>> from dipy.tracking import metrics as tm
>>> import numpy as np
>>> x=np.linspace(0,1,100)
>>> y=0*x
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> m=tm.mean_curvature(xyz) #mean curvature straight line
>>> theta=np.pi*np.linspace(0,1,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> _= tm.mean_curvature(xyz) #mean curvature for semi-circle

mean_orientation#

dipy.tracking.metrics.mean_orientation(xyz)[source]#

Calculates the mean orientation of a curve

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a curve

Returns:
mfloat

Mean orientation.

generate_combinations#

dipy.tracking.metrics.generate_combinations(items, n)[source]#

Combine sets of size n from items

Parameters:
itemssequence
nint
Returns:
iciterator

Examples

>>> from dipy.tracking.metrics import generate_combinations
>>> ic=generate_combinations(range(3),2)
>>> for i in ic: print(i)
[0, 1]
[0, 2]
[1, 2]

longest_track_bundle#

dipy.tracking.metrics.longest_track_bundle(bundle, *, sort=False)[source]#

Return longest track or length sorted track indices in bundle

If sort == True, return the indices of the sorted tracks in the bundle, otherwise return the longest track.

Parameters:
bundlesequence

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

sortbool, optional

If False (default) return longest track. If True, return length sorted indices for tracks in bundle

Returns:
longest_or_indicesarray

longest track - shape (N,3) - (if sort is False), or indices of length sorted tracks (if sort is True)

Examples

>>> from dipy.tracking.metrics import longest_track_bundle
>>> import numpy as np
>>> bundle = [np.array([[0,0,0],[2,2,2]]),np.array([[0,0,0],[4,4,4]])]
>>> longest_track_bundle(bundle)
array([[0, 0, 0],
       [4, 4, 4]])
>>> longest_track_bundle(bundle, sort=True) 
array([0, 1]...)

intersect_sphere#

dipy.tracking.metrics.intersect_sphere(xyz, center, radius)[source]#

If any segment of the track is intersecting with a sphere of specific center and radius return True otherwise False

Parameters:
xyzarray, shape (N,3)

representing x,y,z of the N points of the track

centerarray, shape (3,)

center of the sphere

radiusfloat

radius of the sphere

Returns:
tf{True, False}

True if track xyz intersects sphere

>>> from dipy.tracking.metrics import intersect_sphere
    ..
>>> line=np.array(([0,0,0],[1,1,1],[2,2,2]))
    ..
>>> sph_cent=np.array([1,1,1])
    ..
>>> sph_radius = 1
    ..
>>> intersect_sphere(line,sph_cent,sph_radius)
    ..
True

Notes

The ray to sphere intersection method used here is similar with https://paulbourke.net/geometry/circlesphere/ https://paulbourke.net/geometry/circlesphere/source.cpp we just applied it for every segment neglecting the intersections where the intersecting points are not inside the segment

inside_sphere#

dipy.tracking.metrics.inside_sphere(xyz, center, radius)[source]#

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 \(|x-c|\le r\) where \(x\) a point \(c\) the center of the sphere and \(r\) the radius of the sphere.

Parameters:
xyzarray, shape (N,3)

representing x,y,z of the N points of the track

centerarray, shape (3,)

center of the sphere

radiusfloat

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

inside_sphere_points#

dipy.tracking.metrics.inside_sphere_points(xyz, center, radius)[source]#

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 \(|x-c| \le r\) where \(x\) a point \(c\) the center of the sphere and \(r\) the radius of the sphere.

Parameters:
xyzarray, shape (N,3)

representing x,y,z of the N points of the track

centerarray, shape (3,)

center of the sphere

radiusfloat

radius of the sphere

Returns:
xyznarray, 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]])

spline#

dipy.tracking.metrics.spline(xyz, *, s=3, k=2, nest=-1)[source]#

Generate B-splines as documented in https://scipy-cookbook.readthedocs.io/items/Interpolation.html

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.

Parameters:
xyzarray, shape (N,3)

array representing x,y,z of N points in 3d space

sfloat, optional

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 standard-deviation of y, then a: good s value should be found in the range (m-sqrt(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 s-value. 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 over-estimate 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.

Returns:
xyznarray, shape (M,3)

array representing x,y,z of the M points inside the sphere

See also

scipy.interpolate.splprep
scipy.interpolate.splev

Examples

>>> import numpy as np
>>> t=np.linspace(0,1.75*2*np.pi,100)# make ascending spiral in 3-space
>>> x = np.sin(t)
>>> y = np.cos(t)
>>> z = t
>>> x+= np.random.normal(scale=0.1, size=x.shape) # add noise
>>> y+= np.random.normal(scale=0.1, size=y.shape)
>>> z+= np.random.normal(scale=0.1, size=z.shape)
>>> xyz=np.vstack((x,y,z)).T
>>> xyzn=spline(xyz,s=3,k=2,nest=-1)
>>> len(xyzn) > len(xyz)
True

startpoint#

dipy.tracking.metrics.startpoint(xyz)[source]#

First point of the track

Parameters:
xyzarray, shape(N,3)

Track.

Returns:
sparray, shape(3,)

First track point.

Examples

>>> from dipy.tracking.metrics import startpoint
>>> import numpy as np
>>> theta=np.pi*np.linspace(0,1,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> sp=startpoint(xyz)
>>> sp.any()==xyz[0].any()
True

endpoint#

dipy.tracking.metrics.endpoint(xyz)[source]#
Parameters:
xyzarray, shape(N,3)

Track.

Returns:
eparray, shape(3,)

First track point.

Examples

>>> from dipy.tracking.metrics import endpoint
>>> import numpy as np
>>> theta=np.pi*np.linspace(0,1,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> ep=endpoint(xyz)
>>> ep.any()==xyz[-1].any()
True

arbitrarypoint#

dipy.tracking.metrics.arbitrarypoint(xyz, distance)[source]#

Select an arbitrary point along distance on the track (curve)

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

distancefloat

float representing distance travelled from the xyz[0] point of the curve along the curve.

Returns:
aparray shape (3,)

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

Examples

>>> import numpy as np
>>> from dipy.tracking.metrics import arbitrarypoint, length
>>> theta=np.pi*np.linspace(0,1,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> ap=arbitrarypoint(xyz,length(xyz)/3)

principal_components#

dipy.tracking.metrics.principal_components(xyz)[source]#

We use PCA to calculate the 3 principal directions for a track

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

Returns:
vaarray_like

eigenvalues

vearray_like

eigenvectors

Examples

>>> import numpy as np
>>> from dipy.tracking.metrics import principal_components
>>> theta=np.pi*np.linspace(0,1,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> va, ve = principal_components(xyz)
>>> np.allclose(va, [0.51010101, 0.09883545, 0])
True

midpoint2point#

dipy.tracking.metrics.midpoint2point(xyz, p)[source]#

Calculate distance from midpoint of a curve to arbitrary point p

Parameters:
xyzarray-like shape (N,3)

array representing x,y,z of N points in a track

parray shape (3,)

array representing an arbitrary point with x,y,z coordinates in space.

Returns:
dfloat

a float number representing Euclidean distance

Examples

>>> import numpy as np
>>> from dipy.tracking.metrics import midpoint2point, midpoint
>>> theta=np.pi*np.linspace(0,1,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=0*x
>>> xyz=np.vstack((x,y,z)).T
>>> dist=midpoint2point(xyz,np.array([0,0,0]))

eudx_both_directions#

dipy.tracking.propspeed.eudx_both_directions(seed, ref, qa, ind, odf_vertices, qa_thr, ang_thr, step_sz, total_weight, max_points)#
Parameters:
seedarray, float64 shape (3,)

Point where the tracking starts.

refcnp.npy_intp int

Index of peak to follow first.

qaarray, float64 shape (X, Y, Z, Np)

Anisotropy matrix, where Np is the number of maximum allowed peaks.

indarray, float64 shape(x, y, z, Np)

Index of the track orientation.

odf_verticesdouble array shape (N, 3)

Sampling directions on the sphere.

qa_thrfloat

Threshold for QA, we want everything higher than this threshold.

ang_thrfloat

Angle threshold, we only select fiber orientation within this range.

step_szdouble
total_weightdouble
max_pointscnp.npy_intp
Returns:
trackarray, shape (N,3)

ndarray_offset#

dipy.tracking.propspeed.ndarray_offset(indices, strides, lenind, typesize)#

Find offset in an N-dimensional ndarray using strides

Parameters:
indicesarray, npy_intp shape (N,)

Indices of the array which we want to find the offset.

stridesarray, shape (N,)

Strides of array.

lenindint

len of the indices array.

typesizeint

Number of bytes for data type e.g. if 8 for double, 4 for int32

Returns:
offsetinteger

Index position in flattened array

Examples

>>> import numpy as np
>>> from dipy.tracking.propspeed import ndarray_offset
>>> I=np.array([1,1])
>>> A=np.array([[1,0,0],[0,2,0],[0,0,3]])
>>> S=np.array(A.strides)
>>> ndarray_offset(I,S,2,A.dtype.itemsize)
4
>>> A.ravel()[4]==A[1,1]
True

ActStoppingCriterion#

class dipy.tracking.stopping_criterion.ActStoppingCriterion#

Bases: AnatomicalStoppingCriterion

Anatomically-Constrained Tractography (ACT) stopping criterion.

See [2] for further details about the method.

This implements the use of partial volume fraction (PVE) maps to determine when the tracking stops. The proposed method [2] that cuts streamlines going through subcortical gray matter regions is not implemented here. The backtracking technique for streamlines reaching INVALIDPOINT is not implemented either.

Methods

from_pve(wm_map, gm_map, csf_map, **kw)

AnatomicalStoppingCriterion from partial volume fraction (PVE) maps.

check_point

get_exclude

get_include

References

AnatomicalStoppingCriterion#

class dipy.tracking.stopping_criterion.AnatomicalStoppingCriterion#

Bases: StoppingCriterion

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

cdef:

double interp_out_double[1] double[:] interp_out_view = interp_out_view double[:, :, :] include_map, exclude_map

Methods

from_pve(wm_map, gm_map, csf_map, **kw)

AnatomicalStoppingCriterion from partial volume fraction (PVE) maps.

check_point

get_exclude

get_include

classmethod from_pve(wm_map, gm_map, csf_map, **kw)#

AnatomicalStoppingCriterion from partial volume fraction (PVE) maps.

Parameters:
wm_maparray

The partial volume fraction of white matter at each voxel.

gm_maparray

The partial volume fraction of gray matter at each voxel.

csf_maparray

The partial volume fraction of corticospinal fluid at each voxel.

get_exclude(point)#
get_include(point)#

BinaryStoppingCriterion#

class dipy.tracking.stopping_criterion.BinaryStoppingCriterion#

Bases: StoppingCriterion

cdef:

unsigned char[:, :, :] mask

Methods

check_point

CmcStoppingCriterion#

class dipy.tracking.stopping_criterion.CmcStoppingCriterion#

Bases: AnatomicalStoppingCriterion

Continuous map criterion (CMC) stopping criterion.

This implements the use of partial volume fraction (PVE) maps to determine when the tracking stops [3].

cdef:

double interp_out_double[1] double[:] interp_out_view = interp_out_view double[:, :, :] include_map, exclude_map double step_size double average_voxel_size double correction_factor

Methods

from_pve(wm_map, gm_map, csf_map, **kw)

AnatomicalStoppingCriterion from partial volume fraction (PVE) maps.

check_point

get_exclude

get_include

References

StoppingCriterion#

class dipy.tracking.stopping_criterion.StoppingCriterion#

Bases: object

Methods

check_point

check_point(point)#

StreamlineStatus#

class dipy.tracking.stopping_criterion.StreamlineStatus(value, names=<not given>, *values, module=None, qualname=None, type=None, start=1, boundary=None)#

Bases: IntFlag

Attributes:
denominator

the denominator of a rational number in lowest terms

imag

the imaginary part of a complex number

numerator

the numerator of a rational number in lowest terms

real

the real part of a complex number

Methods

as_integer_ratio(/)

Return a pair of integers, whose ratio is equal to the original int.

bit_count(/)

Number of ones in the binary representation of the absolute value of self.

bit_length(/)

Number of bits necessary to represent self in binary.

conjugate

Returns self, the complex conjugate of any int.

from_bytes(/, bytes[, byteorder, signed])

Return the integer represented by the given array of bytes.

is_integer(/)

Returns True.

to_bytes(/[, length, byteorder, signed])

Return an array of bytes representing an integer.

ENDPOINT = 2#
INVALIDPOINT = 0#
OUTSIDEIMAGE = -1#
PYERROR = -2#
TRACKPOINT = 1#

ThresholdStoppingCriterion#

class dipy.tracking.stopping_criterion.ThresholdStoppingCriterion#

Bases: StoppingCriterion

Methods

check_point

unlist_streamlines#

dipy.tracking.streamline.unlist_streamlines(streamlines)[source]#

Return the streamlines not as a list but as an array and an offset

Parameters:
streamlines: sequence

Streamlines.

Returns:
pointsarray

Streamlines’ points.

offsetsarray

Streamlines’ offsets.

relist_streamlines#

dipy.tracking.streamline.relist_streamlines(points, offsets)[source]#

Given a representation of a set of streamlines as a large array and an offsets array return the streamlines as a list of shorter arrays.

Parameters:
pointsarray

Streamlines’ points.

offsetsarray

Streamlines’ offsets.

Returns:
streamlines: sequence

Streamlines.

center_streamlines#

dipy.tracking.streamline.center_streamlines(streamlines)[source]#

Move streamlines to the origin

Parameters:
streamlineslist

List of 2D ndarrays of shape[-1]==3

Returns:
new_streamlineslist

List of 2D ndarrays of shape[-1]==3

inv_shiftndarray

Translation in x,y,z to go back in the initial position

deform_streamlines#

dipy.tracking.streamline.deform_streamlines(streamlines, deform_field, stream_to_current_grid, current_grid_to_world, stream_to_ref_grid, ref_grid_to_world)[source]#

Apply deformation field to streamlines

Parameters:
streamlineslist

List of 2D ndarrays of shape[-1]==3

deform_field4D numpy array

x,y,z displacements stored in volume, shape[-1]==3

stream_to_current_gridarray, (4, 4)

transform matrix voxmm space to original grid space

current_grid_to_worldarray (4, 4)

transform matrix original grid space to world coordinates

stream_to_ref_gridarray (4, 4)

transform matrix voxmm space to new grid space

ref_grid_to_worldarray(4, 4)

transform matrix new grid space to world coordinates

Returns:
new_streamlineslist

List of the transformed 2D ndarrays of shape[-1]==3

transform_streamlines#

dipy.tracking.streamline.transform_streamlines(streamlines, mat, *, in_place=False)[source]#

Apply affine transformation to streamlines

Parameters:
streamlinesStreamlines

Streamlines object

matarray, (4, 4)

transformation matrix

in_placebool

If True then change data in place. Be careful changes input streamlines.

Returns:
new_streamlinesStreamlines

Sequence transformed 2D ndarrays of shape[-1]==3

select_random_set_of_streamlines#

dipy.tracking.streamline.select_random_set_of_streamlines(streamlines, select, *, rng=None)[source]#

Select a random set of streamlines

Parameters:
streamlinesStreamlines

Object of 2D ndarrays of shape[-1]==3

selectint

Number of streamlines to select. If there are less streamlines than select then select=len(streamlines).

rngRandomState

Default None.

Returns:
selected_streamlineslist

Notes

The same streamline will not be selected twice.

select_by_rois#

dipy.tracking.streamline.select_by_rois(streamlines, affine, rois, include, *, mode=None, tol=None)[source]#

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.

Parameters:
streamlineslist

A list of candidate streamlines for selection

affinearray_like (4, 4)

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). Non-zeros 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 end-points 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.

Returns:
generator

Generates the streamlines to be included based on these criteria.

Notes

The only operation currently possible is “(A or B or …) and not (X or Y or …)”, where A, B are inclusion regions and X, Y are exclusion regions.

Examples

>>> streamlines = [np.array([[0, 0., 0.9],
...                          [1.9, 0., 0.]]),
...                np.array([[0., 0., 0],
...                          [0, 1., 1.],
...                          [0, 2., 2.]]),
...                np.array([[2, 2, 2],
...                          [3, 3, 3]])]
>>> mask1 = np.zeros((4, 4, 4), dtype=bool)
>>> mask2 = np.zeros_like(mask1)
>>> mask1[0, 0, 0] = True
>>> mask2[1, 0, 0] = True
>>> selection = select_by_rois(streamlines, np.eye(4), [mask1, mask2],
...                            [True, True],
...                            tol=1)
>>> list(selection) # The result is a generator
[array([[ 0. ,  0. ,  0.9],
       [ 1.9,  0. ,  0. ]]), array([[ 0.,  0.,  0.],
       [ 0.,  1.,  1.],
       [ 0.,  2.,  2.]])]
>>> selection = select_by_rois(streamlines, np.eye(4), [mask1, mask2],
...                            [True, False],
...                            tol=0.87)
>>> list(selection)
[array([[ 0.,  0.,  0.],
       [ 0.,  1.,  1.],
       [ 0.,  2.,  2.]])]
>>> selection = select_by_rois(streamlines, np.eye(4), [mask1, mask2],
...                            [True, True],
...                            mode="both_end",
...                            tol=1.0)
>>> list(selection)
[array([[ 0. ,  0. ,  0.9],
       [ 1.9,  0. ,  0. ]])]
>>> mask2[0, 2, 2] = True
>>> selection = select_by_rois(streamlines, np.eye(4), [mask1, mask2],
...                            [True, True],
...                            mode="both_end",
...                            tol=1.0)
>>> list(selection)
[array([[ 0. ,  0. ,  0.9],
       [ 1.9,  0. ,  0. ]]), array([[ 0.,  0.,  0.],
       [ 0.,  1.,  1.],
       [ 0.,  2.,  2.]])]

cluster_confidence#

dipy.tracking.streamline.cluster_confidence(streamlines, *, max_mdf=5, subsample=12, power=1, override=False)[source]#

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 [4] (based on the streamline MDF distance from Garyfallidis et al.[5]).

Parameters:
streamlineslist of 2D (N, 3) arrays

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.

Returns:
Returns an array of CCI scores

References

orient_by_rois#

dipy.tracking.streamline.orient_by_rois(streamlines, affine, roi1, roi2, *, in_place=False, as_generator=False)[source]#

Orient a set of streamlines according to a pair of ROIs

Parameters:
streamlineslist or generator

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 (n-by-3 array with ROI coordinate in each row).

in_placebool

Whether to make the change in-place 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.

as_generatorbool

Whether to return a generator as output.

Returns:
streamlineslist or generator

The same 3D arrays as a list or generator, but reoriented with respect to the ROIs

Examples

>>> streamlines = [np.array([[0, 0., 0],
...                          [1, 0., 0.],
...                          [2, 0., 0.]]),
...                np.array([[2, 0., 0.],
...                          [1, 0., 0],
...                          [0, 0,  0.]])]
>>> roi1 = np.zeros((4, 4, 4), dtype=bool)
>>> roi2 = np.zeros_like(roi1)
>>> roi1[0, 0, 0] = True
>>> roi2[1, 0, 0] = True
>>> orient_by_rois(streamlines, np.eye(4), roi1, roi2)
[array([[ 0.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 2.,  0.,  0.]]), array([[ 0.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 2.,  0.,  0.]])]

orient_by_streamline#

dipy.tracking.streamline.orient_by_streamline(streamlines, standard, *, n_points=12, in_place=False, as_generator=False)[source]#

Orient a bundle of streamlines to a standard streamline.

Parameters:
streamlinesStreamlines, list

The input streamlines to orient.

standardStreamlines, list, or ndarrray

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 in-place 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.

as_generatorbool

Whether to return a generator as output.

Returns:
Streamlineswith each individual array oriented to be as similar as

possible to the standard.

values_from_volume#

dipy.tracking.streamline.values_from_volume(data, streamlines, affine)[source]#

Extract values of a scalar/vector along each streamline from a volume.

Parameters:
data3D or 4D array

Scalar (for 3D) and vector (for 4D) values to be extracted. For 4D data, interpolation will be done on the 3 spatial dimensions in each volume.

streamlinesndarray or list

If array, of shape (n_streamlines, n_nodes, 3) If list, len(n_streamlines) with (n_nodes, 3) array in each element of the list.

affinearray_like (4, 4)

The mapping from voxel coordinates to streamline points. The voxel_to_rasmm matrix, typically from a NIFTI file.

Returns:
array or list (depending on the input)values interpolate to each

coordinate along the length of each streamline.

Notes

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.

nbytes#

dipy.tracking.streamline.nbytes(streamlines)[source]#

Streamlines#

dipy.tracking.streamlinespeed.Streamlines#

alias of ArraySequence

as_native_array#

dipy.tracking.streamlinespeed.as_native_array(arr)[source]#

Return arr as native byteordered array

If arr is already native byte ordered, return unchanged. If it is opposite endian, then make a native byte ordered copy and return that

Parameters:
arrndarray
Returns:
native_arrndarray

If arr was native order, this is just arr. Otherwise it’s a new array such that np.all(native_arr == arr), with native byte ordering.

compress_streamlines#

dipy.tracking.streamlinespeed.compress_streamlines(streamlines, tol_error=0.01, max_segment_length=10)#

Compress streamlines by linearization.

The compression [6] 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 [6] 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 [6] (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 :footcite:p:`Rheault2015. 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.

Parameters:
streamlinesone or a list of array-like of shape (N,3)

Array representing x,y,z of N points in a streamline.

tol_errorfloat, optional

Tolerance error in mm. A rule of thumb is to set it to 0.01mm for deterministic streamlines and 0.1mm for probabilitic streamlines.

max_segment_lengthfloat, optional

Maximum length in mm of any given segment produced by the compression. The default is 10mm. (In [6] they used a value of np.inf).

Returns:
compressed_streamlinesone or a list of array-like

Results of the linearization process.

Notes

Be aware that compressed streamlines have variable step sizes. One needs to be careful when computing streamlines-based metrics [7].

References

Examples

>>> from dipy.tracking.streamlinespeed import compress_streamlines
>>> import numpy as np
>>> # One streamline: a wiggling line
>>> rng = np.random.RandomState(42)
>>> streamline = np.linspace(0, 10, 100*3).reshape((100, 3))
>>> streamline += 0.2 * rng.rand(100, 3)
>>> c_streamline = compress_streamlines(streamline, tol_error=0.2)
>>> len(streamline)
100
>>> len(c_streamline)
10
>>> # Multiple streamlines
>>> streamlines = [streamline, streamline[::2]]
>>> c_streamlines = compress_streamlines(streamlines, tol_error=0.2)
>>> [len(s) for s in streamlines]
[100, 50]
>>> [len(s) for s in c_streamlines]
[10, 7]

length#

dipy.tracking.streamlinespeed.length(streamlines)#

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

set_number_of_points#

dipy.tracking.streamlinespeed.set_number_of_points(streamlines, nb_points=3)#
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]

density_map#

dipy.tracking.utils.density_map(streamlines, affine, vol_dims)[source]#

Count the number of unique streamlines that pass through each voxel.

Parameters:
streamlinesiterable

A sequence of streamlines.

affinearray_like (4, 4)

The mapping from voxel coordinates to streamline points. The voxel_to_rasmm matrix, typically from a NIFTI file.

vol_dims3 ints

The shape of the volume to be returned containing the streamlines counts

Returns:
image_volumendarray, shape=vol_dims

The number of streamline points in each voxel of volume.

Raises:
IndexError

When the points of the streamlines lie outside of the return volume.

Notes

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.

connectivity_matrix#

dipy.tracking.utils.connectivity_matrix(streamlines, affine, label_volume, *, inclusive=False, symmetric=True, return_mapping=False, mapping_as_streamlines=False)[source]#

Count the streamlines that start and end at each label pair.

Parameters:
streamlinessequence

A sequence of streamlines.

affinearray_like (4, 4)

The mapping from voxel coordinates to streamline coordinates. The voxel_to_rasmm matrix, typically from a NIFTI file.

label_volumendarray

An image volume with an integer data type, where the intensities in the volume map to anatomical structures.

inclusive: bool

Whether to analyze the entire streamline, as opposed to just the endpoints. False by default.

symmetricbool, True by default

Symmetric means we don’t distinguish between start and end points. If symmetric is True, matrix[i, j] == matrix[j, i].

return_mappingbool, False by default

If True, a mapping is returned which maps matrix indices to streamlines.

mapping_as_streamlinesbool, False by default

If True voxel indices map to lists of streamline objects. Otherwise voxel indices map to lists of integers.

Returns:
matrixndarray

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

ndbincount#

dipy.tracking.utils.ndbincount(x, *, weights=None, shape=None)[source]#

Like bincount, but for nd-indices.

Parameters:
xarray_like (N, M)

M indices to a an Nd-array

weightsarray_like (M,), optional

Weights associated with indices

shapeoptional

the shape of the output

reduce_labels#

dipy.tracking.utils.reduce_labels(label_volume)[source]#

Reduce an array of labels to the integers from 0 to n with smallest possible n.

Examples

>>> labels = np.array([[1, 3, 9],
...                    [1, 3, 8],
...                    [1, 3, 7]])
>>> new_labels, lookup = reduce_labels(labels)
>>> lookup
array([1, 3, 7, 8, 9])
>>> new_labels 
array([[0, 1, 4],
       [0, 1, 3],
       [0, 1, 2]]...)
>>> (lookup[new_labels] == labels).all()
True

subsegment#

dipy.tracking.utils.subsegment(streamlines, max_segment_length)[source]#

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.

Parameters:
streamlinessequence of ndarrays

The streamlines to be subsegmented.

max_segment_lengthfloat

The longest allowable segment length.

Returns:
output_streamlinesgenerator

A set of streamlines.

Notes

Segments of 0 length are removed. If unchanged

Examples

>>> streamlines = [np.array([[0,0,0],[2,0,0],[5,0,0]])]
>>> list(subsegment(streamlines, 3.))
[array([[ 0.,  0.,  0.],
       [ 2.,  0.,  0.],
       [ 5.,  0.,  0.]])]
>>> list(subsegment(streamlines, 1))
[array([[ 0.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 2.,  0.,  0.],
       [ 3.,  0.,  0.],
       [ 4.,  0.,  0.],
       [ 5.,  0.,  0.]])]
>>> list(subsegment(streamlines, 1.6))
[array([[ 0. ,  0. ,  0. ],
       [ 1. ,  0. ,  0. ],
       [ 2. ,  0. ,  0. ],
       [ 3.5,  0. ,  0. ],
       [ 5. ,  0. ,  0. ]])]

seeds_from_mask#

dipy.tracking.utils.seeds_from_mask(mask, affine, *, density=(1, 1, 1))[source]#

Create seeds for fiber tracking from a binary mask.

Seeds points are placed evenly distributed in all voxels of mask which are True.

Parameters:
maskbinary 3d array_like

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.

Raises:
ValueError

When mask is not a three-dimensional array

Examples

>>> mask = np.zeros((3,3,3), 'bool')
>>> mask[0,0,0] = 1
>>> seeds_from_mask(mask, np.eye(4), density=[1,1,1])
array([[ 0.,  0.,  0.]])

random_seeds_from_mask#

dipy.tracking.utils.random_seeds_from_mask(mask, affine, *, seeds_count=1, seed_count_per_voxel=True, random_seed=None)[source]#

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.

Parameters:
maskbinary 3d array_like

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

Raises:
ValueError

When mask is not a three-dimensional array

See also

seeds_from_mask

Examples

>>> mask = np.zeros((3,3,3), 'bool')
>>> mask[0,0,0] = 1
>>> random_seeds_from_mask(mask, np.eye(4), seeds_count=1,
... seed_count_per_voxel=True, random_seed=1)
array([[-0.23838787, -0.20150886,  0.31422574]])
>>> random_seeds_from_mask(mask, np.eye(4), seeds_count=6,
... seed_count_per_voxel=True, random_seed=1)
array([[-0.23838787, -0.20150886,  0.31422574],
       [-0.41435083, -0.26318949,  0.30127447],
       [ 0.44305611,  0.01132755,  0.47624371],
       [ 0.30500292,  0.30794079,  0.01532556],
       [ 0.03816435, -0.15672913, -0.13093276],
       [ 0.12509547,  0.3972138 ,  0.27568569]])
>>> mask[0,1,2] = 1
>>> random_seeds_from_mask(mask, np.eye(4),
... seeds_count=2, seed_count_per_voxel=True, random_seed=1)
array([[ 0.30500292,  1.30794079,  2.01532556],
       [-0.23838787, -0.20150886,  0.31422574],
       [ 0.3702492 ,  0.78681721,  2.10314815],
       [-0.41435083, -0.26318949,  0.30127447]])

target#

dipy.tracking.utils.target(streamlines, affine, target_mask, *, include=True)[source]#

Filter streamlines based on whether or not they pass through an ROI.

Parameters:
streamlinesiterable

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.

target_maskarray-like

A mask used as a target. Non-zero values are considered to be within the target region.

includebool, default True

If True, streamlines passing through target_mask are kept. If False, the streamlines not passing through target_mask are kept.

Returns:
streamlinesgenerator

A sequence of streamlines that pass through target_mask.

Raises:
ValueError

When the points of the streamlines lie outside of the target_mask.

See also

density_map

target_line_based#

dipy.tracking.utils.target_line_based(streamlines, affine, target_mask, *, include=True)[source]#

Filter streamlines based on whether or not they pass through a ROI, using a line-based algorithm. Mostly used as a replacement of target for compressed streamlines.

This function never returns single-point streamlines, whatever the value of include.

See [8] and [7] for further details about the method.

Parameters:
streamlinesiterable

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.

target_maskarray-like

A mask used as a target. Non-zero values are considered to be within the target region.

includebool, default True

If True, streamlines passing through target_mask are kept. If False, the streamlines not passing through target_mask are kept.

Returns:
streamlinesgenerator

A sequence of streamlines that pass through target_mask.

See also

dipy.tracking.utils.density_map
dipy.tracking.streamline.compress_streamlines

References

streamline_near_roi#

dipy.tracking.utils.streamline_near_roi(streamline, roi_coords, tol, *, mode='any')[source]#

Is a streamline near an ROI.

Implements the inner loops of the near_roi() function.

Parameters:
streamlinearray, shape (N, 3)

A single streamline

roi_coordsarray, shape (M, 3)

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 end-points is within tol from ROI

“both_end” : both end points are within tol from ROI.

Returns:
outboolean

near_roi#

dipy.tracking.utils.near_roi(streamlines, affine, region_of_interest, *, tol=None, mode='any')[source]#

Provide filtering criteria for a set of streamlines based on whether they fall within a tolerance distance from an ROI.

Parameters:
streamlineslist or generator

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. Non-zero 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 end-points is within tol from ROI

“both_end” : both end points are within tol from ROI.

Returns:
1D array of boolean dtype, shape (len(streamlines), )
This contains True for indices corresponding to each streamline
that passes within a tolerance distance from the target ROI, False
otherwise.

length#

dipy.tracking.utils.length(streamlines)[source]#

Calculate the lengths of many streamlines in a bundle.

Parameters:
streamlineslist

Each item in the list is an array with 3D coordinates of a streamline.

Returns:
Iterator object which then computes the length of each
streamline in the bundle, upon iteration.

unique_rows#

dipy.tracking.utils.unique_rows(in_array, *, dtype='f4')[source]#

Find the unique rows in an array.

Parameters:
in_array: ndarray

The array for which the unique rows should be found

dtype: str, optional

This determines the intermediate representation used for the values. Should at least preserve the values of the input array.

Returns:
u_return: ndarray

Array with the unique rows of the original array.

transform_tracking_output#

dipy.tracking.utils.transform_tracking_output(tracking_output, affine, *, save_seeds=False)[source]#

Apply a linear transformation, given by affine, to streamlines.

Parameters:
tracking_outputStreamlines generator

Either streamlines (list, ArraySequence) or a tuple with streamlines and seeds together

affinearray (4, 4)

The mapping between voxel indices and the point space for seeds. The voxel_to_rasmm matrix, typically from a NIFTI file.

save_seedsbool, optional

If set, seeds associated to streamlines will be also moved and returned

Returns:
streamlinesgenerator

A generator for the sequence of transformed streamlines. If save_seeds is True, also return a generator for the transformed seeds.

reduce_rois#

dipy.tracking.utils.reduce_rois(rois, include)[source]#

Reduce multiple ROIs to one inclusion and one exclusion ROI.

Parameters:
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). Non-zeros 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.

Returns:
include_roiboolean 3D array

An array marking the inclusion mask.

exclude_roiboolean 3D array

An array marking the exclusion mask

Notes

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.

path_length#

dipy.tracking.utils.path_length(streamlines, affine, aoi, *, fill_value=-1)[source]#

Compute the shortest path, along any streamline, between aoi and each voxel.

Parameters:
streamlinesseq of (N, 3) arrays

A sequence of streamlines, path length is given in mm along the curve of the streamline.

aoiarray, 3d

A mask (binary array) of voxels from which to start computing distance.

affinearray (4, 4)

The mapping between voxel indices and the point space for seeds. The voxel_to_rasmm matrix, typically from a NIFTI file.

fill_valuefloat

The value of voxel in the path length map that are not connected to the aoi.

Returns:
plmarray

Same shape as aoi. The minimum distance between every point and aoi along the path of a streamline.

max_angle_from_curvature#

dipy.tracking.utils.max_angle_from_curvature(min_radius_curvature, step_size)[source]#

Get the maximum deviation angle from the minimum radius curvature.

See [9] for further details about the method.

Parameters:
min_radius_curvature: float

Minimum radius of curvature in mm.

step_size: float

The tracking step size in mm.

Returns:
max_angle: float

The maximum deviation angle in radian, given the radius curvature and the step size.

References

min_radius_curvature_from_angle#

dipy.tracking.utils.min_radius_curvature_from_angle(max_angle, step_size)[source]#

Get minimum radius of curvature from a deviation angle.

See [9] for further details about the method.

Parameters:
max_angle: float

The maximum deviation angle in radian. theta should be between [0 - pi/2] otherwise default will be pi/2.

step_size: float

The tracking step size in mm.

Returns:
min_radius_curvature: float

Minimum radius of curvature in mm, given the maximum deviation angle theta and the step size.

References

seeds_directions_pairs#

dipy.tracking.utils.seeds_directions_pairs(positions, peaks, *, max_cross=-1)[source]#

Pair each seed to the corresponding peaks. If multiple peaks are available the seed is repeated for each.

Parameters:
positionsarray, (N, 3)

Coordinates of the N positions.

peaksarray (N, M, 3)

Peaks at each position

max_crossint, optional

The maximum number of direction to track from each seed in crossing voxels. By default all voxel peaks are used.

Returns:
seedsarray (K, 3)
directionsarray (K, 3)

streamline_mapping#

dipy.tracking.vox2track.streamline_mapping(streamlines, affine=None, mapping_as_streamlines=False)#

Creates a mapping from voxel indices to streamlines.

Returns a dictionary where each key is a 3d voxel index and the associated value is a list of the streamlines that pass through that voxel.

Parameters:
streamlinessequence

A sequence of streamlines.

affinearray_like (4, 4), optional

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.

Returns:
mappingdefaultdict(list)

A mapping from voxel indices to the streamlines that pass through that voxel.

Examples

>>> streamlines = [np.array([[0., 0., 0.],
...                          [1., 1., 1.],
...                          [2., 3., 4.]]),
...                np.array([[0., 0., 0.],
...                          [1., 2., 3.]])]
>>> mapping = streamline_mapping(streamlines, affine=np.eye(4))
>>> mapping[0, 0, 0]
[0, 1]
>>> mapping[1, 1, 1]
[0]
>>> mapping[1, 2, 3]
[1]
>>> mapping.get((3, 2, 1), 'no streamlines')
'no streamlines'
>>> mapping = streamline_mapping(streamlines, affine=np.eye(4),
...                              mapping_as_streamlines=True)
>>> mapping[1, 2, 3][0] is streamlines[1]
True

track_counts#

dipy.tracking.vox2track.track_counts(tracks, vol_dims, vox_sizes=(1, 1, 1), return_elements=True)#

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.

Parameters:
trackssequence

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 n-th point in the t-th track. Points are of form x, y, z in voxel mm coordinates. That is, if i, j, k is the possibly non-integer 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)

Returns:
tcsndarray shape vol_dim

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.

Examples

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:

>>> vox_track0 = np.array([[0, 0, 0], [1.1, 2.2, 3.3], [2.2, 4.4, 6.6]])
>>> vox_track1 = np.array([[0, 0, 0], [0, 0, 1], [0, 0, 2]])
>>> vs = (vox_track0, vox_track1)
>>> vox_dim = (10, 20, 30) # original voxel array size
>>> tcs = track_counts(vs, vox_dim, (1, 1, 1), False)
>>> tcs.shape
(10, 20, 30)
>>> tcs[0, 0, 0:4]
array([2, 1, 1, 0])
>>> tcs[1, 2, 3], tcs[2, 4, 7]
(1, 1)

You can also use the routine to count into larger-than-voxel boxes. To do this, increase the voxel size and decrease the vox_dim accordingly:

>>> tcs=track_counts(vs, (10/2., 20/2., 30/2.), (2,2,2), False)
>>> tcs.shape
(5, 10, 15)
>>> tcs[1,1,2], tcs[1,2,3]
(1, 1)