"""Metrics for tracks, where tracks are arrays of points"""
import numpy as np
from scipy.interpolate import splev, splprep
from dipy.testing.decorators import warning_for_keywords
[docs]
def winding(xyz):
"""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
----------
xyz : array-like shape (N,3)
Array representing x,y,z of N points in a track.
Returns
-------
a : scalar
Total turning angle in degrees.
"""
U, s, V = np.linalg.svd(xyz - np.mean(xyz, axis=0), 0)
proj = np.dot(U[:, 0:2], np.diag(s[0:2]))
turn = 0
for j in range(len(xyz) - 1):
v0 = proj[j]
v1 = proj[j + 1]
v = np.dot(v0, v1) / (np.linalg.norm(v0) * np.linalg.norm(v1))
v = np.clip(v, -1, 1)
tmp = np.arccos(v)
turn += tmp
return np.rad2deg(turn)
[docs]
@warning_for_keywords()
def length(xyz, *, along=False):
"""Euclidean length of track line
This will give length in mm if tracks are expressed in world coordinates.
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
along : bool, optional
If True, return array giving cumulative length along track,
otherwise (default) return scalar giving total length.
Returns
-------
L : scalar 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])
"""
xyz = np.asarray(xyz)
if xyz.shape[0] < 2:
if along:
return np.array([0])
return 0
dists = np.sqrt((np.diff(xyz, axis=0) ** 2).sum(axis=1))
if along:
return np.cumsum(dists)
return np.sum(dists)
[docs]
def bytes(xyz):
"""Size of track in bytes.
Parameters
----------
xyz : array-like shape (N,3)
Array representing x,y,z of N points in a track.
Returns
-------
b : int
Number of bytes.
"""
return xyz.nbytes
[docs]
def midpoint(xyz):
"""Midpoint of track
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
Returns
-------
mp : array 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. ])
"""
xyz = np.asarray(xyz)
n_pts = xyz.shape[0]
if n_pts == 0:
raise ValueError("xyz array cannot be empty")
if n_pts == 1:
return xyz.copy().squeeze()
cumlen = np.zeros(n_pts)
cumlen[1:] = length(xyz, along=True)
midlen = cumlen[-1] / 2.0
ind = np.where((cumlen - midlen) > 0)[0][0]
len0 = cumlen[ind - 1]
len1 = cumlen[ind]
Ds = midlen - len0
Lambda = Ds / (len1 - len0)
return Lambda * xyz[ind] + (1 - Lambda) * xyz[ind - 1]
[docs]
def center_of_mass(xyz):
"""Center of mass of streamline
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
Returns
-------
com : array 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.])
"""
xyz = np.asarray(xyz)
if xyz.size == 0:
raise ValueError("xyz array cannot be empty")
return np.mean(xyz, axis=0)
[docs]
@warning_for_keywords()
def magn(xyz, *, n=1):
"""magnitude of vector"""
mag = np.sum(xyz**2, axis=1) ** 0.5
imag = np.where(mag == 0)
mag[imag] = np.finfo(float).eps
if n > 1:
return np.tile(mag, (n, 1)).T
return mag.reshape(len(mag), 1)
[docs]
def frenet_serret(xyz):
r"""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
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
Returns
-------
T : array shape (N,3)
array representing the tangent of the curve xyz
N : array shape (N,3)
array representing the normal of the curve xyz
B : array shape (N,3)
array representing the binormal of the curve xyz
k : array shape (N,1)
array representing the curvature of the curve xyz
t : array shape (N,1)
array representing the torsion of the curve xyz
Examples
--------
Create a helix and calculate its tangent, normal, binormal, curvature
and torsion
>>> from dipy.tracking import metrics as tm
>>> import numpy as np
>>> theta = 2*np.pi*np.linspace(0,2,100)
>>> x=np.cos(theta)
>>> y=np.sin(theta)
>>> z=theta/(2*np.pi)
>>> xyz=np.vstack((x,y,z)).T
>>> T,N,B,k,t=tm.frenet_serret(xyz)
"""
xyz = np.asarray(xyz)
n_pts = xyz.shape[0]
if n_pts == 0:
raise ValueError("xyz array cannot be empty")
dxyz = np.gradient(xyz)[0]
ddxyz = np.gradient(dxyz)[0]
# Tangent
T = np.divide(dxyz, magn(dxyz, n=3))
# Derivative of Tangent
dT = np.gradient(T)[0]
# Normal
N = np.divide(dT, magn(dT, n=3))
# Binormal
B = np.cross(T, N)
# Curvature
k = magn(np.cross(dxyz, ddxyz), n=1) / (magn(dxyz, n=1) ** 3)
# Torsion
# (In matlab was t=dot(-B,N,2))
t = np.sum(-B * N, axis=1)
# return T,N,B,k,t,dxyz,ddxyz,dT
return T, N, B, k, t
[docs]
def mean_curvature(xyz):
"""Calculates the mean curvature of a curve
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a curve
Returns
-------
m : float
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
"""
xyz = np.asarray(xyz)
n_pts = xyz.shape[0]
if n_pts == 0:
raise ValueError("xyz array cannot be empty")
dxyz = np.gradient(xyz)[0]
ddxyz = np.gradient(dxyz)[0]
# Curvature
k = magn(np.cross(dxyz, ddxyz), n=1) / (magn(dxyz, n=1) ** 3)
return np.mean(k)
[docs]
def mean_orientation(xyz):
"""
Calculates the mean orientation of a curve
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a curve
Returns
-------
m : float
Mean orientation.
"""
xyz = np.asarray(xyz)
n_pts = xyz.shape[0]
if n_pts == 0:
raise ValueError("xyz array cannot be empty")
dxyz = np.gradient(xyz)[0]
return np.mean(dxyz, axis=0)
[docs]
def generate_combinations(items, n):
"""Combine sets of size n from items
Parameters
----------
items : sequence
n : int
Returns
-------
ic : iterator
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]
"""
if n == 0:
yield []
elif n == 2:
# if n=2 non_recursive
for i in range(len(items) - 1):
for j in range(i + 1, len(items)):
yield [i, j]
else:
# if n>2 uses recursion
for i in range(len(items)):
for cc in generate_combinations(items[i + 1 :], n - 1):
yield [items[i]] + cc
[docs]
@warning_for_keywords()
def longest_track_bundle(bundle, *, sort=False):
"""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
----------
bundle : sequence
of tracks as arrays, shape (N1,3) ... (Nm,3)
sort : bool, optional
If False (default) return longest track. If True, return length
sorted indices for tracks in bundle
Returns
-------
longest_or_indices : array
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) #doctest: +ELLIPSIS
array([0, 1]...)
"""
alllengths = [length(t) for t in bundle]
alllengths = np.array(alllengths)
if sort:
ilongest = alllengths.argsort()
return ilongest
else:
ilongest = alllengths.argmax()
return bundle[ilongest]
[docs]
def 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
Parameters
----------
xyz : array, shape (N,3)
representing x,y,z of the N points of the track
center : array, shape (3,)
center of the sphere
radius : float
radius of the sphere
Returns
-------
tf : {True, False}
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
"""
center = np.array(center)
# print center
lt = xyz.shape[0]
for i in range(lt - 1):
# first point
x1 = xyz[i]
# second point
x2 = xyz[i + 1]
# do the calculations as given in the Notes
x = x2 - x1
a = np.inner(x, x)
x1c = x1 - center
b = 2 * np.inner(x, x1c)
c = (
np.inner(center, center)
+ np.inner(x1, x1)
- 2 * np.inner(center, x1)
- radius**2
)
bb4ac = b * b - 4 * a * c
# print 'bb4ac',bb4ac
if abs(a) < np.finfo(float).eps or bb4ac < 0: # too small segment or
# no intersection
continue
if bb4ac == 0: # one intersection point p
mu = -b / 2 * a
p = x1 + mu * x
# check if point is inside the segment
# print 'p',p
if np.inner(p - x1, p - x1) <= a:
return True
if bb4ac > 0: # two intersection points p1 and p2
mu = (-b + np.sqrt(bb4ac)) / (2 * a)
p1 = x1 + mu * x
mu = (-b - np.sqrt(bb4ac)) / (2 * a)
p2 = x1 + mu * x
# check if points are inside the line segment
# print 'p1,p2',p1,p2
if np.inner(p1 - x1, p1 - x1) <= a or np.inner(p2 - x1, p2 - x1) <= a:
return True
return False
[docs]
def inside_sphere(xyz, center, radius):
r"""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
----------
xyz : array, shape (N,3)
representing x,y,z of the N points of the track
center : array, shape (3,)
center of the sphere
radius : float
radius of the sphere
Returns
-------
tf : {True,False}
Whether point is inside sphere.
Examples
--------
>>> from dipy.tracking.metrics import inside_sphere
>>> line=np.array(([0,0,0],[1,1,1],[2,2,2]))
>>> sph_cent=np.array([1,1,1])
>>> sph_radius = 1
>>> inside_sphere(line,sph_cent,sph_radius)
True
"""
return (np.sqrt(np.sum((xyz - center) ** 2, axis=1)) <= radius).any()
[docs]
def inside_sphere_points(xyz, center, radius):
r"""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
----------
xyz : array, shape (N,3)
representing x,y,z of the N points of the track
center : array, shape (3,)
center of the sphere
radius : float
radius of the sphere
Returns
-------
xyzn : array, shape(M,3)
array representing x,y,z of the M points inside the sphere
Examples
--------
>>> from dipy.tracking.metrics import inside_sphere_points
>>> line=np.array(([0,0,0],[1,1,1],[2,2,2]))
>>> sph_cent=np.array([1,1,1])
>>> sph_radius = 1
>>> inside_sphere_points(line,sph_cent,sph_radius)
array([[1, 1, 1]])
"""
return xyz[(np.sqrt(np.sum((xyz - center) ** 2, axis=1)) <= radius)]
[docs]
@warning_for_keywords()
def spline(xyz, *, s=3, k=2, nest=-1):
"""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
----------
xyz : array, shape (N,3)
array representing x,y,z of N points in 3d space
s : float, 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.
k : int, 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.
nest : None 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
-------
xyzn : array, shape (M,3)
array representing x,y,z of the M points inside the sphere
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
See Also
--------
scipy.interpolate.splprep
scipy.interpolate.splev
"""
# find the knot points
tckp, u = splprep([xyz[:, 0], xyz[:, 1], xyz[:, 2]], s=s, k=k, nest=nest)
# evaluate spline, including interpolated points
xnew, ynew, znew = splev(np.linspace(0, 1, 400), tckp)
return np.vstack((xnew, ynew, znew)).T
[docs]
def startpoint(xyz):
"""First point of the track
Parameters
----------
xyz : array, shape(N,3)
Track.
Returns
-------
sp : array, 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
"""
return xyz[0]
[docs]
def endpoint(xyz):
"""
Parameters
----------
xyz : array, shape(N,3)
Track.
Returns
-------
ep : array, 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
"""
return xyz[-1]
[docs]
def arbitrarypoint(xyz, distance):
"""Select an arbitrary point along distance on the track (curve)
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
distance : float
float representing distance travelled from the xyz[0] point of
the curve along the curve.
Returns
-------
ap : array 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)
"""
xyz = np.asarray(xyz)
n_pts = xyz.shape[0]
if n_pts == 0:
raise ValueError("xyz array cannot be empty")
if n_pts == 1:
return xyz.copy().squeeze()
cumlen = np.zeros(n_pts)
cumlen[1:] = length(xyz, along=True)
if cumlen[-1] < distance:
raise ValueError("Given distance is bigger than the length of the curve")
ind = np.where((cumlen - distance) > 0)[0][0]
len0 = cumlen[ind - 1]
len1 = cumlen[ind]
Ds = distance - len0
Lambda = Ds / (len1 - len0)
return Lambda * xyz[ind] + (1 - Lambda) * xyz[ind - 1]
def _extrap(xyz, cumlen, distance):
"""Helper function for extrapolate"""
ind = np.where((cumlen - distance) > 0)[0][0]
len0 = cumlen[ind - 1]
len1 = cumlen[ind]
Ds = distance - len0
Lambda = Ds / (len1 - len0)
return Lambda * xyz[ind] + (1 - Lambda) * xyz[ind - 1]
[docs]
def principal_components(xyz):
"""We use PCA to calculate the 3 principal directions for a track
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
Returns
-------
va : array_like
eigenvalues
ve : array_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
"""
C = np.cov(xyz.T)
va, ve = np.linalg.eig(C)
return va, ve
[docs]
def midpoint2point(xyz, p):
"""Calculate distance from midpoint of a curve to arbitrary point p
Parameters
----------
xyz : array-like shape (N,3)
array representing x,y,z of N points in a track
p : array shape (3,)
array representing an arbitrary point with x,y,z coordinates in
space.
Returns
-------
d : float
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]))
"""
mid = midpoint(xyz)
return np.sqrt(np.sum((xyz - mid) ** 2))