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
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))
|
Calculate distances between list of tracks A and list of tracks B
All tracks need to have the same number of points
Parameters: | tracksA : sequence
tracksB : sequence
|
---|---|
Returns: | DM : array, shape (len(tracksA), len(tracksB))
|
See also
dipy.metrics.downsample
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
|
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
>>> print(res[0])
[[ 1. 1.5 0. 0.70710683 0. ]]
>>> print(res[1])
[[ 2. 2.5 0. 0.70710677 0. ]]
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)
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
sqd_trh: float :
|
---|---|
Returns: | C : dict
|
Generate a first pass clustering using 3 points on the tracks only.
Parameters: | tracks : sequence
indices : None or sequence, optional
trh : float, optional
|
---|---|
Returns: | C : dict
|
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 fvtk
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)
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: | 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
>>> lee_angle_distance([0,0,0],[1,0,0],[3,4,5],[5,4,3])
2.0
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: | 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
>>> d = lee_perpendicular_distance([0,0,0],[1,0,0],[3,4,5],[5,4,3])
>>> print('%.6f' % d)
5.787888
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
|
---|---|
Returns: | C : dict
|
See also
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 fvtk module:
from dipy.viz import fvtk
r=fvtk.ren()
for c in C:
color=np.random.rand(3)
for i in C[c]['indices']:
fvtk.add(r,fvtk.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(tracks, d_thr=0.5)
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
|
---|---|
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
>>> 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)
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
>>> 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
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
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
>>> 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
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