/usr/lib/python2.7/dist-packages/dipy/segment/tests/test_metric.py is in python-dipy 0.10.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 | import numpy as np
import dipy.segment.metric as dipymetric
import itertools
from nose.tools import assert_true, assert_false, assert_equal, assert_almost_equal
from numpy.testing import assert_array_equal, assert_raises, run_module_suite
def norm(x, ord=None, axis=None):
if axis is not None:
return np.apply_along_axis(np.linalg.norm, axis, x.astype(np.float64), ord)
return np.linalg.norm(x.astype(np.float64), ord=ord)
dtype = "float32"
# Create wiggling streamline
nb_points = 18
rng = np.random.RandomState(42)
x = np.linspace(0, 10, nb_points)
y = rng.rand(nb_points)
z = np.sin(np.linspace(0, np.pi, nb_points)) # Bending
s = np.array([x, y, z], dtype=dtype).T
# Create trivial streamlines
s1 = np.array([np.arange(10, dtype=dtype)]*3).T # 10x3
s2 = np.arange(3*10, dtype=dtype).reshape((-1, 3))[::-1] # 10x3
s3 = np.array([np.arange(5, dtype=dtype)]*4) # 5x4
s4 = np.array([np.arange(5, dtype=dtype)]*3) # 5x3
streamlines = [s, s1, s2, s3, s4]
def test_metric_minimum_average_direct_flip():
feature = dipymetric.IdentityFeature()
class MinimumAverageDirectFlipMetric(dipymetric.Metric):
def __init__(self, feature):
super(MinimumAverageDirectFlipMetric, self).__init__(feature=feature)
@property
def is_order_invariant(self):
return True # Ordering is handled in the distance computation
def are_compatible(self, shape1, shape2):
return shape1[0] == shape2[0]
def dist(self, v1, v2):
average_euclidean = lambda x, y: np.mean(norm(x-y, axis=1))
dist_direct = average_euclidean(v1, v2)
dist_flipped = average_euclidean(v1, v2[::-1])
return min(dist_direct, dist_flipped)
for metric in [MinimumAverageDirectFlipMetric(feature),
dipymetric.MinimumAverageDirectFlipMetric(feature)]:
# Test special cases of the MDF distance.
assert_equal(metric.dist(s, s), 0.)
assert_equal(metric.dist(s, s[::-1]), 0.)
# Translation
offset = np.array([0.8, 1.3, 5], dtype=dtype)
assert_almost_equal(metric.dist(s, s+offset), norm(offset), 5)
# Scaling
M_scaling = np.diag([1.2, 2.8, 3]).astype(dtype)
s_mean = np.mean(s, axis=0)
s_zero_mean = s - s_mean
s_scaled = np.dot(M_scaling, s_zero_mean.T).T + s_mean
d = np.mean(norm((np.diag(M_scaling)-1)*s_zero_mean, axis=1))
assert_almost_equal(metric.dist(s, s_scaled), d, 5)
# Rotation
from dipy.core.geometry import rodrigues_axis_rotation
rot_axis = np.array([1, 2, 3], dtype=dtype)
M_rotation = rodrigues_axis_rotation(rot_axis, 60.).astype(dtype)
s_mean = np.mean(s, axis=0)
s_zero_mean = s - s_mean
s_rotated = np.dot(M_rotation, s_zero_mean.T).T + s_mean
opposite = norm(np.cross(rot_axis, s_zero_mean), axis=1) / norm(rot_axis)
distances = np.sqrt(2*opposite**2 * (1 - np.cos(60.*np.pi/180.))).astype(dtype)
d = np.mean(distances)
assert_almost_equal(metric.dist(s, s_rotated), d, 5)
for s1, s2 in itertools.product(*[streamlines]*2): # All possible pairs
# Extract features since metric doesn't work directly on streamlines
f1 = metric.feature.extract(s1)
f2 = metric.feature.extract(s2)
# Test method are_compatible
same_nb_points = f1.shape[0] == f2.shape[0]
assert_equal(metric.are_compatible(f1.shape, f2.shape), same_nb_points)
# Test method dist if features are compatible
if metric.are_compatible(f1.shape, f2.shape):
distance = metric.dist(f1, f2)
if np.all(f1 == f2):
assert_equal(distance, 0.)
assert_almost_equal(distance, dipymetric.dist(metric, s1, s2))
assert_almost_equal(distance, dipymetric.mdf(s1, s2))
assert_true(distance >= 0.)
# This metric type is order invariant
assert_true(metric.is_order_invariant)
for s1, s2 in itertools.product(*[streamlines]*2): # All possible pairs
f1 = metric.feature.extract(s1)
f2 = metric.feature.extract(s2)
if not metric.are_compatible(f1.shape, f2.shape):
continue
f1_flip = metric.feature.extract(s1[::-1])
f2_flip = metric.feature.extract(s2[::-1])
distance = metric.dist(f1, f2)
assert_almost_equal(metric.dist(f1_flip, f2_flip), distance)
if not np.all(f1_flip == f2_flip):
assert_true(np.allclose(metric.dist(f1, f2_flip), distance))
assert_true(np.allclose(metric.dist(f1_flip, f2), distance))
def test_metric_cosine():
feature = dipymetric.VectorOfEndpointsFeature()
class CosineMetric(dipymetric.Metric):
def __init__(self, feature):
super(CosineMetric, self).__init__(feature=feature)
def are_compatible(self, shape1, shape2):
# Cosine metric works on vectors.
return shape1 == shape2 and shape1[0] == 1
def dist(self, v1, v2):
# Check if we have null vectors
if norm(v1) == 0:
return 0. if norm(v2) == 0 else 1.
v1_normed = v1.astype(np.float64) / norm(v1.astype(np.float64))
v2_normed = v2.astype(np.float64) / norm(v2.astype(np.float64))
cos_theta = np.dot(v1_normed, v2_normed.T)
# Make sure it's in [-1, 1], i.e. within domain of arccosine
cos_theta = np.minimum(cos_theta, 1.)
cos_theta = np.maximum(cos_theta, -1.)
return np.arccos(cos_theta) / np.pi # Normalized cosine distance
for metric in [CosineMetric(feature), dipymetric.CosineMetric(feature)]:
# Test special cases of the cosine distance.
v0 = np.array([[0, 0, 0]], dtype=np.float32)
v1 = np.array([[1, 2, 3]], dtype=np.float32)
v2 = np.array([[1, -1./2, 0]], dtype=np.float32)
v3 = np.array([[-1, -2, -3]], dtype=np.float32)
assert_equal(metric.dist(v0, v0), 0.) # dot-dot
assert_equal(metric.dist(v0, v1), 1.) # dot-line
assert_equal(metric.dist(v1, v1), 0.) # collinear
assert_equal(metric.dist(v1, v2), 0.5) # orthogonal
assert_equal(metric.dist(v1, v3), 1.) # opposite
for s1, s2 in itertools.product(*[streamlines]*2): # All possible pairs
# Extract features since metric doesn't work directly on streamlines
f1 = metric.feature.extract(s1)
f2 = metric.feature.extract(s2)
# Test method are_compatible
are_vectors = f1.shape[0] == 1 and f2.shape[0] == 1
same_dimension = f1.shape[1] == f2.shape[1]
assert_equal(metric.are_compatible(f1.shape, f2.shape),
are_vectors and same_dimension)
# Test method dist if features are compatible
if metric.are_compatible(f1.shape, f2.shape):
distance = metric.dist(f1, f2)
if np.all(f1 == f2):
assert_almost_equal(distance, 0.)
assert_almost_equal(distance, dipymetric.dist(metric, s1, s2))
assert_true(distance >= 0.)
assert_true(distance <= 1.)
# This metric type is not order invariant
assert_false(metric.is_order_invariant)
for s1, s2 in itertools.product(*[streamlines]*2): # All possible pairs
f1 = metric.feature.extract(s1)
f2 = metric.feature.extract(s2)
if not metric.are_compatible(f1.shape, f2.shape):
continue
f1_flip = metric.feature.extract(s1[::-1])
f2_flip = metric.feature.extract(s2[::-1])
distance = metric.dist(f1, f2)
assert_almost_equal(metric.dist(f1_flip, f2_flip), distance)
if not np.all(f1_flip == f2_flip):
assert_false(metric.dist(f1, f2_flip) == distance)
assert_false(metric.dist(f1_flip, f2) == distance)
def test_subclassing_metric():
class EmptyMetric(dipymetric.Metric):
pass
metric = EmptyMetric()
assert_raises(NotImplementedError, metric.are_compatible, None, None)
assert_raises(NotImplementedError, metric.dist, None, None)
def test_distance_matrix():
metric = dipymetric.SumPointwiseEuclideanMetric()
for dtype in [np.int32, np.int64, np.float32, np.float64]:
# Compute distances of all tuples spawn by the Cartesian product
# of `data` with itself.
data = (np.random.rand(4, 10, 3)*10).astype(dtype)
D = dipymetric.distance_matrix(metric, data)
assert_equal(D.shape, (len(data), len(data)))
assert_array_equal(np.diag(D), np.zeros(len(data)))
if metric.is_order_invariant:
# Distance matrix should be symmetric
assert_array_equal(D, D.T)
for i in range(len(data)):
for j in range(len(data)):
assert_equal(D[i, j], dipymetric.dist(metric, data[i], data[j]))
# Compute distances of all tuples spawn by the Cartesian product
# of `data` with `data2`.
data2 = (np.random.rand(3, 10, 3)*10).astype(dtype)
D = dipymetric.distance_matrix(metric, data, data2)
assert_equal(D.shape, (len(data), len(data2)))
for i in range(len(data)):
for j in range(len(data2)):
assert_equal(D[i, j], dipymetric.dist(metric, data[i], data2[j]))
if __name__ == '__main__':
run_module_suite()
|