|
3 | 3 |
|
4 | 4 |
|
5 | 5 | @pytest.fixture |
6 | | -def ellipsoid_points(): |
| 6 | +def generate_pointclouds(): |
7 | 7 | """ |
8 | | - Generates points on the surface of an ellipsoid. |
| 8 | + Generates synthetic point clouds for testing. |
9 | 9 | """ |
10 | | - a, b, c = 5, 3, 2 # Semi-axes of the ellipsoid |
11 | | - u = np.linspace(0, 2 * np.pi, 30) # Longitude-like angles |
12 | | - v = np.linspace(0, np.pi, 30) # Latitude-like angles |
13 | | - u, v = np.meshgrid(u, v) |
14 | | - x = a * np.cos(u) * np.sin(v) |
15 | | - y = b * np.sin(u) * np.sin(v) |
16 | | - z = c * np.cos(v) |
17 | | - return np.column_stack((x.ravel(), y.ravel(), z.ravel())) |
18 | | - |
19 | | - |
20 | | -@pytest.fixture |
21 | | -def droplet_points(): |
22 | | - """ |
23 | | - Create a droplet point cloud. |
24 | | - """ |
25 | | - from napari_stress import sample_data |
26 | | - |
27 | | - pointcloud = sample_data.get_droplet_point_cloud()[0][0][:, 1:] |
28 | | - |
29 | | - return pointcloud |
30 | | - |
31 | | - |
32 | | -def test_spherical_harmonics_fit(droplet_points): |
| 10 | + def ellipsoid_points(a=5, b=3, c=2): |
| 11 | + u = np.linspace(0, 2 * np.pi, 30) # Longitude-like angles |
| 12 | + v = np.linspace(0, np.pi, 30) # Latitude-like angles |
| 13 | + u, v = np.meshgrid(u, v) |
| 14 | + x = a * np.cos(u) * np.sin(v) |
| 15 | + y = b * np.sin(u) * np.sin(v) |
| 16 | + z = c * np.cos(v) |
| 17 | + return np.column_stack((x.ravel(), y.ravel(), z.ravel())) |
| 18 | + |
| 19 | + def droplet_points(): |
| 20 | + from napari_stress import sample_data |
| 21 | + return sample_data.get_droplet_point_cloud()[0][0][:, 1:] |
| 22 | + |
| 23 | + def random_ellipsoid_points(a0=10, a1=20, a2=30, x0=(0, 0, 0)): |
| 24 | + import vedo |
| 25 | + ellipsoid = vedo.Ellipsoid( |
| 26 | + pos=x0, |
| 27 | + axis1=(a0, 0, 0), |
| 28 | + axis2=(0, a1, 0), |
| 29 | + axis3=(0, 0, a2), |
| 30 | + ) |
| 31 | + points = ellipsoid.vertices |
| 32 | + # Random rotations |
| 33 | + angle_z = np.radians(np.random.randint(0, 360)) |
| 34 | + angle_x = np.radians(np.random.randint(0, 180)) |
| 35 | + R_z = np.array([[np.cos(angle_z), -np.sin(angle_z), 0], |
| 36 | + [np.sin(angle_z), np.cos(angle_z), 0], |
| 37 | + [0, 0, 1]]) |
| 38 | + R_x = np.array([[1, 0, 0], |
| 39 | + [0, np.cos(angle_x), -np.sin(angle_x)], |
| 40 | + [0, np.sin(angle_x), np.cos(angle_x)]]) |
| 41 | + return points @ R_z @ R_x |
| 42 | + |
| 43 | + return ellipsoid_points, droplet_points, random_ellipsoid_points |
| 44 | + |
| 45 | + |
| 46 | +def test_spherical_harmonics(generate_pointclouds): |
33 | 47 | """ |
34 | | - Test fitting spherical harmonics to ellipsoid points. |
| 48 | + Test spherical harmonics fitting and expansion. |
35 | 49 | """ |
36 | 50 | from napari_stress import approximation |
37 | 51 |
|
38 | | - expander = approximation.SphericalHarmonicsExpander( |
39 | | - max_degree=5, expansion_type="cartesian" |
40 | | - ) |
41 | | - expander.fit(droplet_points) |
42 | | - |
43 | | - assert expander.coefficients_ is not None |
44 | | - assert expander.coefficients_.shape == (3, (5 + 1), (5 + 1)) |
45 | | - |
46 | | - |
47 | | -def test_spherical_harmonics_expand(droplet_points): |
48 | | - """ |
49 | | - Test expanding points using spherical harmonics. |
50 | | - """ |
51 | | - from napari_stress import approximation |
| 52 | + _, droplet_points, _ = generate_pointclouds |
| 53 | + points = droplet_points() |
| 54 | + degree = 20 |
52 | 55 |
|
| 56 | + # Test fitting |
53 | 57 | expander = approximation.SphericalHarmonicsExpander( |
54 | | - max_degree=10, expansion_type="cartesian" |
| 58 | + max_degree=degree, expansion_type="cartesian" |
55 | 59 | ) |
56 | | - expander.fit(droplet_points) |
57 | | - expanded_points = expander.expand(droplet_points) |
| 60 | + expander.fit(points) |
| 61 | + assert expander.coefficients_ is not None |
| 62 | + assert expander.coefficients_.shape == (3, (degree + 1), (degree + 1)) |
58 | 63 |
|
59 | | - # Assert the expanded points are close to the original points |
60 | | - np.testing.assert_allclose(expanded_points, droplet_points, rtol=0.01) |
| 64 | + # Test expansion |
| 65 | + expanded_points = expander.expand(points) |
| 66 | + np.testing.assert_allclose(expanded_points, points, rtol=0.01) |
61 | 67 |
|
| 68 | + # Test radial expansion |
62 | 69 | expander_radial = approximation.SphericalHarmonicsExpander( |
63 | 70 | max_degree=10, expansion_type="radial" |
64 | 71 | ) |
65 | | - expander_radial.fit(droplet_points) |
66 | | - expanded_points_radial = expander_radial.expand(droplet_points) |
67 | | - |
68 | | - # Assert the expanded points are close to the original points |
69 | | - np.testing.assert_allclose( |
70 | | - expanded_points_radial, droplet_points, rtol=0.01 |
71 | | - ) |
72 | | - |
73 | | - |
74 | | -def generate_ellipsoidal_pointclouds( |
75 | | - a0: float = 10, a1: float = 20, a2: float = 30, x0: tuple = (0, 0, 0) |
76 | | -): |
77 | | - import vedo |
78 | | - |
79 | | - ellipsoid = vedo.Ellipsoid( |
80 | | - pos=x0, |
81 | | - axis1=(a0, 0, 0), |
82 | | - axis2=(0, a1, 0), |
83 | | - axis3=(0, 0, a2), |
84 | | - ) |
85 | | - |
86 | | - points = ellipsoid.vertices |
| 72 | + expander_radial.fit(points) |
| 73 | + expanded_points_radial = expander_radial.expand(points) |
| 74 | + np.testing.assert_allclose(expanded_points_radial, points, rtol=0.01) |
87 | 75 |
|
88 | | - # rotate all points around the z-axis by random angle between |
89 | | - # 0 and 360 degrees using only numpy |
90 | | - angle = np.radians(np.random.randint(0, 360)) |
91 | | - c, s = np.cos(angle), np.sin(angle) |
92 | | - R = np.array(((c, -s, 0), (s, c, 0), (0, 0, 1))) |
93 | | - points = points @ R |
94 | 76 |
|
95 | | - # rotate also around x-axis by random angle between |
96 | | - # 0 and 180 degrees using only numpy |
97 | | - angle = np.radians(np.random.randint(0, 180)) |
98 | | - c, s = np.cos(angle), np.sin(angle) |
99 | | - R = np.array(((1, 0, 0), (0, c, -s), (0, s, c))) |
100 | | - points = points @ R |
101 | | - |
102 | | - return points |
103 | | - |
104 | | - |
105 | | -def test_lsq_ellipsoid0(n_tests=10): |
| 77 | +def test_ellipsoid_fitting(generate_pointclouds): |
| 78 | + """ |
| 79 | + Test ellipsoid fitting and properties. |
| 80 | + """ |
106 | 81 | from napari_stress import approximation |
107 | 82 |
|
108 | | - a0 = 10 |
109 | | - a1 = 20 |
110 | | - a2 = 30 |
| 83 | + _, _, random_ellipsoid_points = generate_pointclouds |
| 84 | + a0, a1, a2 = 10, 20, 30 |
111 | 85 | x0 = np.asarray([0, 0, 0]) |
112 | 86 |
|
113 | | - max_mean_curvatures = [] |
114 | | - min_mean_curvatures = [] |
115 | | - for i in range(n_tests): |
116 | | - points = generate_ellipsoidal_pointclouds(a0=a0, a1=a1, a2=a2, x0=x0) |
117 | | - |
118 | | - expander = approximation.EllipsoidExpander() |
119 | | - expander.fit(points) |
120 | | - |
121 | | - expanded_points = expander.expand(points) |
122 | | - center = expander.center_ |
123 | | - axes_lengths = np.sort(expander.axes_) |
124 | | - |
125 | | - max_mean_curvatures.append( |
126 | | - expander.properties["maximum_mean_curvature"] |
127 | | - ) |
128 | | - min_mean_curvatures.append( |
129 | | - expander.properties["minimum_mean_curvature"] |
130 | | - ) |
131 | | - |
132 | | - # The mean curvature of an ellipsoid is given by the formula |
133 | | - # H = a / (2 * c^2) + a / (2 * b^2) for the maximum mean curvature |
134 | | - # and H = c / (2 * b^2) + c / (2 * a^2) for the minimum mean curvature |
135 | | - # where a, b, c are the semi-axes of the ellipsoid |
136 | | - axes_sorted = np.sort(expander.axes_) |
137 | | - a = axes_sorted[2] |
138 | | - b = axes_sorted[1] |
139 | | - c = axes_sorted[0] |
140 | | - h_max_theory = a / (2 * c**2) + a / (2 * b**2) |
141 | | - h_min_theory = c / (2 * b**2) + c / (2 * a**2) |
142 | | - assert expander.properties["maximum_mean_curvature"] == h_max_theory |
143 | | - assert expander.properties["minimum_mean_curvature"] == h_min_theory |
144 | | - assert ( |
145 | | - expander.properties["maximum_mean_curvature"] |
146 | | - > expander.properties["minimum_mean_curvature"] |
147 | | - ) |
148 | | - |
149 | | - assert np.allclose(center, x0) |
150 | | - assert np.allclose(a2, axes_lengths[2]) |
151 | | - assert np.allclose(a1, axes_lengths[1]) |
152 | | - assert np.allclose(a0, axes_lengths[0]) |
153 | | - |
154 | | - assert np.allclose( |
155 | | - expander.properties["residuals"].mean(), 0, atol=0.01 |
156 | | - ) |
157 | | - assert len(expanded_points) == len(points) |
158 | | - |
159 | | - # all elements in either list should be the same |
160 | | - assert np.allclose(max_mean_curvatures, max_mean_curvatures[0]) |
161 | | - assert np.allclose(min_mean_curvatures, min_mean_curvatures[0]) |
162 | | - |
163 | | - |
164 | | -def test_lsq_ellipsoid1(): |
165 | | - "Test whether pproperties in class are correctly set." |
166 | | - from napari_stress import approximation |
| 87 | + points = random_ellipsoid_points(a0, a1, a2, x0) |
167 | 88 |
|
| 89 | + # Test fitting |
168 | 90 | expander = approximation.EllipsoidExpander() |
| 91 | + expander.fit(points) |
| 92 | + expanded_points = expander.expand(points) |
| 93 | + center = expander.center_ |
| 94 | + axes_lengths = np.sort(expander.axes_) |
| 95 | + |
| 96 | + # Validate properties |
| 97 | + assert np.allclose(center, x0) |
| 98 | + assert np.allclose(a2, axes_lengths[2]) |
| 99 | + assert np.allclose(a1, axes_lengths[1]) |
| 100 | + assert np.allclose(a0, axes_lengths[0]) |
| 101 | + assert np.allclose(expander.properties["residuals"].mean(), 0, atol=0.01) |
| 102 | + assert len(expanded_points) == len(points) |
| 103 | + |
| 104 | + # Validate mean curvature |
| 105 | + max_curvature = expander.properties["maximum_mean_curvature"] |
| 106 | + min_curvature = expander.properties["minimum_mean_curvature"] |
| 107 | + axes_sorted = np.sort(expander.axes_) |
| 108 | + a, b, c = axes_sorted[2], axes_sorted[1], axes_sorted[0] |
| 109 | + h_max_theory = a / (2 * c**2) + a / (2 * b**2) |
| 110 | + h_min_theory = c / (2 * b**2) + c / (2 * a**2) |
| 111 | + assert np.allclose(max_curvature, h_max_theory) |
| 112 | + assert np.allclose(min_curvature, h_min_theory) |
| 113 | + assert max_curvature > min_curvature |
| 114 | + |
| 115 | + |
| 116 | +def test_curvature_and_normals(generate_pointclouds, make_napari_viewer): |
| 117 | + """ |
| 118 | + Test curvature and normals on ellipsoid. |
| 119 | + """ |
| 120 | + from napari_stress import approximation, measurements, types |
169 | 121 |
|
170 | | - random_coefficients = np.random.random((3, 2, 3)) |
171 | | - expander.coefficients_ = random_coefficients |
172 | | - |
173 | | - assert np.array_equal(expander.center_, random_coefficients[0, 0]) |
174 | | - assert np.array_equal( |
175 | | - expander.axes_, np.linalg.norm(random_coefficients[:, 1], axis=1) |
176 | | - ) |
177 | | - |
178 | | - |
179 | | -def test_lsq_ellipsoid(): |
180 | | - from napari_stress import approximation, get_droplet_point_cloud |
181 | | - |
182 | | - pointcloud = get_droplet_point_cloud()[0][0][:, 1:] |
183 | | - ellipsoid = approximation.least_squares_ellipsoid(pointcloud) |
184 | | - |
185 | | - fitted_points = approximation.expand_points_on_ellipse( |
186 | | - ellipsoid, pointcloud |
187 | | - ) |
188 | | - |
189 | | - assert fitted_points is not None |
190 | | - |
191 | | - |
192 | | -def test_ellipse_normals(): |
193 | | - from napari_stress import approximation, get_droplet_point_cloud |
194 | | - |
195 | | - pointcloud = get_droplet_point_cloud()[0][0][:, 1:] |
196 | | - |
197 | | - normals = approximation.normals_on_ellipsoid(pointcloud) |
198 | | - |
199 | | - assert normals is not None |
200 | | - |
| 122 | + _, droplet_points, _ = generate_pointclouds |
| 123 | + points = droplet_points() |
201 | 124 |
|
202 | | -def test_curvature_on_ellipsoid(make_napari_viewer): |
203 | | - from napari_stress import ( |
204 | | - approximation, |
205 | | - get_droplet_point_cloud, |
206 | | - measurements, |
207 | | - types, |
208 | | - ) |
| 125 | + # Fit ellipsoid |
| 126 | + ellipsoid = approximation.least_squares_ellipsoid(points) |
| 127 | + fitted_points = approximation.expand_points_on_ellipse(ellipsoid, points) |
209 | 128 |
|
210 | | - pointcloud = get_droplet_point_cloud()[0][0][:, 1:] |
211 | | - ellipsoid_stress = approximation.least_squares_ellipsoid(pointcloud) |
212 | | - fitted_points_stress = approximation.expand_points_on_ellipse( |
213 | | - ellipsoid_stress, pointcloud |
214 | | - ) |
| 129 | + # Test curvature |
215 | 130 | data, features, metadata = measurements.curvature_on_ellipsoid( |
216 | | - ellipsoid_stress, fitted_points_stress |
| 131 | + ellipsoid, fitted_points |
217 | 132 | ) |
218 | | - |
219 | 133 | assert data is not None |
220 | 134 | assert features is not None |
221 | 135 | assert metadata is not None |
222 | 136 |
|
223 | | - viewer = make_napari_viewer() |
224 | | - viewer.add_points(fitted_points_stress) |
225 | | - viewer.add_vectors(ellipsoid_stress) |
226 | | - results = measurements.curvature_on_ellipsoid( |
227 | | - ellipsoid_stress, fitted_points_stress |
228 | | - ) |
| 137 | + # Test normals |
| 138 | + normals = approximation.normals_on_ellipsoid(points) |
| 139 | + assert normals is not None |
229 | 140 |
|
| 141 | + # Validate metadata keys |
| 142 | + viewer = make_napari_viewer() |
| 143 | + viewer.add_points(fitted_points) |
| 144 | + viewer.add_vectors(ellipsoid) |
| 145 | + results = measurements.curvature_on_ellipsoid(ellipsoid, fitted_points) |
230 | 146 | assert types._METADATAKEY_H_E123_ELLIPSOID in results[1]["metadata"].keys() |
231 | 147 | assert types._METADATAKEY_H0_ELLIPSOID in results[1]["metadata"].keys() |
232 | 148 | assert types._METADATAKEY_MEAN_CURVATURE in results[1]["features"].keys() |
0 commit comments