Optimized track distances, similarities and clustering algorithms using track distances
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 : | xyz : array(N,3)
alpha : float
|
---|---|
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
>>> from fos.tracking.distances import approx_polygon_track
>>> #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
Implementation of Lee et al Approximate Trajectory Partitioning Algorithm
This is base on the minimum description length principle
Parameters : | xyz : array(N,3)
alpha : float
|
---|---|
Returns : | characteristic_points : list of M array(3,) points |
Calculate distances between list of tracks A and list of tracks B
Parameters : | tracksA : sequence
tracksB : sequence
metric : str
|
---|---|
Returns : | DM : array, shape (len(tracksA), len(tracksB))
|
Extract divergence vectors and points of intersection between planes normal to the reference fiber and other tracks
Parameters : | tracks : sequence
ref : array, shape (N,3)
|
---|---|
Returns : | hits : sequence
|
Examples
>>> from dipy.tracking.distances import cut_plane
>>> 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')]
>>> cut_plane(bundlex,refx)
[array([[ 1. , 1.5 , 0. , 0.70710683, 0,]], dtype=float32),
array([[ 2. , 2.5 , 0. , 0.70710677, 0.]], dtype=float32)]
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.
Intersect Segment S(t) = sa +t(sb-sa), 0 <=t<= 1 against cylinder specified by p,q and r Look p.197 from Real Time Collision Detection C. Ericson
Examples
>>> from dipy.tracking.distances import intersect_segment_cylinder as isc
>>> # Define cylinder using a segment defined by
>>> p=np.array([0,0,0],dtype=float32)
>>> q=np.array([1,0,0],dtype=float32)
>>> r=0.5
>>> # Define segment
>>> sa=np.array([0.5,1 ,0],dtype=float32)
>>> sb=np.array([0.5,-1,0],dtype=float32)
>>> isc(sa,sb,p,q,r)
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 ———– C : graph with clusters
of indices 3tracks (tracks consisting of 3 points only)
Returns : | C : dict
|
---|
Generate a first pass clustering using 3 points (first, mid and last) on the tracks only.
Parameters : | tracks : sequence
indices : sequence
trh : float
|
---|---|
Returns : | C : dict, 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
>>> import numpy as np
>>> from dipy.viz import fvtk
>>> from dipy.tracking.distances as pf
>>> 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=pf.larch_fast_split(tracks,None,0.5)
Here is an example of how to visualize the clustering above
r=fvtk.ren() fvtk.add(r,fvtk.line(tracks,fvtk.red)) fvtk.show(r) for c in C:
- color=np.random.rand(3)
- for i in C[c][‘indices’]:
- fos.add(r,fvtk.line(tracks[i],color))
fvtk.show(r) for c in C:
fvtk.add(r,fos.line(C[c][‘rep3’]/C[c][‘N’],fos.white))
fvtk.show(r)
This function assumes that norm(end0-start0)>norm(end1-start1) i.e. that the first segment will be bigger than the second one.
Parameters : | start0 : float array(3,) end0 : float array(3,) start1 : float array(3,) end1 : float array(3,) |
---|---|
Returns : | angle_distance : float |
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
>>> import dipy.tracking.distances as pf
>>> pf.lee_angle_distance([0,0,0],[1,0,0],[3,4,5],[5,4,3])
2.0
This function assumes that norm(end0-start0)>norm(end1-start1) i.e. that the first segment will be bigger than the second one. Parameters ————
start0 : float array(3,) end0 : float array(3,) start1 : float array(3,) end1 : float 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))
Examples
>>> import dipy.core.performance as pf
>>> pf.lee_perpendicular_distance([0,0,0],[1,0,0],[3,4,5],[5,4,3])
5.9380966767403436
Efficient tractography clustering
Every track can needs to have the same number of points. Use dipy.tracking.metrics.downsample to restrict the number of points
Parameters : | tracks : sequence
d_thr : float, average euclidean distance threshold |
---|---|
Returns : | C : dict |
See also
Notes
The distance calculated between two tracks
t_1 t_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 fliped
Examples
>>> from dipy.tracking.distances import local_skeleton_clustering
>>> 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,3)
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 : | tracks : sequence
d_thr : float, average euclidean distance threshold |
---|---|
Returns : | C : dict |
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
>>> from dipy.viz import fvtk
>>> 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)
Min/Max/Mean Average Minimum Distance between tracks xyz1 and xyz2
Based on the metrics in Zhang, Correia, Laidlaw 2008 http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4479455 which in turn are based on those of Corouge et al. 2004
Parameters : | xyz1 : array, shape (N1,3), dtype float32 xyz2 : array, shape (N2,3), dtype float32
metrics : {‘avg’,’min’,’max’,’all’}
|
---|---|
Returns : | avg_mcd : float
min_mcd : float
max_mcd : float
|
Notes
Algorithmic description
Lets 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)
Find the minimum distance between two curves xyz1, xyz2
Parameters : | xyz1 : array, shape (N1,3), dtype float32 xyz2 : array, shape (N2,3), dtype float32
|
---|---|
Returns : | md : minimum distance |
Notes
Algorithmic description
Lets 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
Find the most similar track in a bundle using distances calculated from Zhang et. al 2008.
Parameters : | tracks : sequence
metric : str
|
---|---|
Returns : | si : int
s : array, shape (len(tracks),)
|
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
Euclidean (L2) norm of length 3 vector
Parameters : | vec : array-like shape (3,) |
---|---|
Returns : | norm : float
|
Return normalized 3D vector
Vector divided by Euclidean (L2) norm
Parameters : | vec : array-like shape (3,) |
---|---|
Returns : | vec_out : array shape (3,) |
Calculate the squared distance from a point c to a finite line segment ab.
Examples
>>> from dipy.tracking.distances import point_segment_sq_distance
>>> a=np.array([0,0,0],dtype=float32)
>>> b=np.array([1,0,0],dtype=float32)
>>> c=np.array([0,1,0],dtype=float32)
>>> point_segment_sq_distance(a,b,c)
1.0
>>> c=np.array([0,3,0],dtype=float32)
>>> pf.point_segment_sq_distance(a,b,c)
9.0
>>> c=np.array([-1,1,0],dtype=float32)
>>> pf.point_segment_sq_distance(a,b,c)
2.0
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
>>> from dipy.tracking.distances import point_track_sq_distance_check
>>> 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
Calculate the euclidean distance between two 3pt tracks both direct and flip distances are calculated but only the smallest is returned
Parameters : | a : array, shape (3,3) a three point track : b : array, shape (3,3) a three point track : |
---|---|
Returns : | dist :float : |
Examples
>>> import numpy as np
>>> from dipy.tracking.distances import track_dist_3pts
>>> 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]])
>>> track_dist_3pts(a,b)
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
>>> from dipy.tracking.distances import track_roi_intersection_check
>>> 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