Mercurial > hg > orthanc-stone
comparison UnitTestsSources/GeometryToolboxTests.cpp @ 1679:5b8b88e5bfd6
successfully running unit tests in WebAssembly
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 24 Nov 2020 12:59:10 +0100 |
parents | |
children | 9ac2a65d4172 |
comparison
equal
deleted
inserted
replaced
1678:1393e3393a0b | 1679:5b8b88e5bfd6 |
---|---|
1 /** | |
2 * Stone of Orthanc | |
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics | |
4 * Department, University Hospital of Liege, Belgium | |
5 * Copyright (C) 2017-2020 Osimis S.A., Belgium | |
6 * | |
7 * This program is free software: you can redistribute it and/or | |
8 * modify it under the terms of the GNU Affero General Public License | |
9 * as published by the Free Software Foundation, either version 3 of | |
10 * the License, or (at your option) any later version. | |
11 * | |
12 * This program is distributed in the hope that it will be useful, but | |
13 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 * Affero General Public License for more details. | |
16 * | |
17 * You should have received a copy of the GNU Affero General Public License | |
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. | |
19 **/ | |
20 | |
21 | |
22 #include <gtest/gtest.h> | |
23 | |
24 #include "../OrthancStone/Sources/Toolbox/FiniteProjectiveCamera.h" | |
25 #include "../OrthancStone/Sources/Toolbox/GenericToolbox.h" | |
26 #include "../OrthancStone/Sources/Toolbox/GeometryToolbox.h" | |
27 | |
28 #include <Logging.h> | |
29 #include <OrthancException.h> | |
30 | |
31 #include <boost/math/special_functions/round.hpp> | |
32 | |
33 | |
34 TEST(GeometryToolbox, Interpolation) | |
35 { | |
36 using namespace OrthancStone::GeometryToolbox; | |
37 | |
38 // https://en.wikipedia.org/wiki/Bilinear_interpolation#Application_in_image_processing | |
39 ASSERT_DOUBLE_EQ(146.1f, ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95)); | |
40 | |
41 ASSERT_DOUBLE_EQ(91, ComputeBilinearInterpolationUnitSquare(0, 0, 91, 210, 162, 95)); | |
42 ASSERT_DOUBLE_EQ(210, ComputeBilinearInterpolationUnitSquare(1, 0, 91, 210, 162, 95)); | |
43 ASSERT_DOUBLE_EQ(162, ComputeBilinearInterpolationUnitSquare(0, 1, 91, 210, 162, 95)); | |
44 ASSERT_DOUBLE_EQ(95, ComputeBilinearInterpolationUnitSquare(1, 1, 91, 210, 162, 95)); | |
45 | |
46 ASSERT_DOUBLE_EQ(123.35f, ComputeTrilinearInterpolationUnitSquare | |
47 (0.5f, 0.2f, 0.7f, | |
48 91, 210, 162, 95, | |
49 51, 190, 80, 92)); | |
50 | |
51 ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95), | |
52 ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 0, | |
53 91, 210, 162, 95, | |
54 51, 190, 80, 92)); | |
55 | |
56 ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 51, 190, 80, 92), | |
57 ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 1, | |
58 91, 210, 162, 95, | |
59 51, 190, 80, 92)); | |
60 } | |
61 | |
62 | |
63 static bool CompareMatrix(const OrthancStone::Matrix& a, | |
64 const OrthancStone::Matrix& b, | |
65 double threshold = 0.00000001) | |
66 { | |
67 if (a.size1() != b.size1() || | |
68 a.size2() != b.size2()) | |
69 { | |
70 return false; | |
71 } | |
72 | |
73 for (size_t i = 0; i < a.size1(); i++) | |
74 { | |
75 for (size_t j = 0; j < a.size2(); j++) | |
76 { | |
77 if (fabs(a(i, j) - b(i, j)) > threshold) | |
78 { | |
79 LOG(ERROR) << "Too large difference in component (" | |
80 << i << "," << j << "): " << a(i,j) << " != " << b(i,j); | |
81 return false; | |
82 } | |
83 } | |
84 } | |
85 | |
86 return true; | |
87 } | |
88 | |
89 | |
90 static bool CompareVector(const OrthancStone::Vector& a, | |
91 const OrthancStone::Vector& b, | |
92 double threshold = 0.00000001) | |
93 { | |
94 if (a.size() != b.size()) | |
95 { | |
96 return false; | |
97 } | |
98 | |
99 for (size_t i = 0; i < a.size(); i++) | |
100 { | |
101 if (fabs(a(i) - b(i)) > threshold) | |
102 { | |
103 LOG(ERROR) << "Too large difference in component " | |
104 << i << ": " << a(i) << " != " << b(i); | |
105 return false; | |
106 } | |
107 } | |
108 | |
109 return true; | |
110 } | |
111 | |
112 | |
113 | |
114 TEST(FiniteProjectiveCamera, Decomposition1) | |
115 { | |
116 // Example 6.2 of "Multiple View Geometry in Computer Vision - 2nd | |
117 // edition" (page 163) | |
118 const double p[12] = { | |
119 3.53553e+2, 3.39645e+2, 2.77744e+2, -1.44946e+6, | |
120 -1.03528e+2, 2.33212e+1, 4.59607e+2, -6.32525e+5, | |
121 7.07107e-1, -3.53553e-1, 6.12372e-1, -9.18559e+2 | |
122 }; | |
123 | |
124 OrthancStone::FiniteProjectiveCamera camera(p); | |
125 ASSERT_EQ(3u, camera.GetMatrix().size1()); | |
126 ASSERT_EQ(4u, camera.GetMatrix().size2()); | |
127 ASSERT_EQ(3u, camera.GetIntrinsicParameters().size1()); | |
128 ASSERT_EQ(3u, camera.GetIntrinsicParameters().size2()); | |
129 ASSERT_EQ(3u, camera.GetRotation().size1()); | |
130 ASSERT_EQ(3u, camera.GetRotation().size2()); | |
131 ASSERT_EQ(3u, camera.GetCenter().size()); | |
132 | |
133 ASSERT_NEAR(1000.0, camera.GetCenter()[0], 0.01); | |
134 ASSERT_NEAR(2000.0, camera.GetCenter()[1], 0.01); | |
135 ASSERT_NEAR(1500.0, camera.GetCenter()[2], 0.01); | |
136 | |
137 ASSERT_NEAR(468.2, camera.GetIntrinsicParameters() (0, 0), 0.1); | |
138 ASSERT_NEAR(91.2, camera.GetIntrinsicParameters() (0, 1), 0.1); | |
139 ASSERT_NEAR(300.0, camera.GetIntrinsicParameters() (0, 2), 0.1); | |
140 ASSERT_NEAR(427.2, camera.GetIntrinsicParameters() (1, 1), 0.1); | |
141 ASSERT_NEAR(200.0, camera.GetIntrinsicParameters() (1, 2), 0.1); | |
142 ASSERT_NEAR(1.0, camera.GetIntrinsicParameters() (2, 2), 0.1); | |
143 | |
144 ASSERT_NEAR(0, camera.GetIntrinsicParameters() (1, 0), 0.0000001); | |
145 ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 0), 0.0000001); | |
146 ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 1), 0.0000001); | |
147 | |
148 ASSERT_NEAR(0.41380, camera.GetRotation() (0, 0), 0.00001); | |
149 ASSERT_NEAR(0.90915, camera.GetRotation() (0, 1), 0.00001); | |
150 ASSERT_NEAR(0.04708, camera.GetRotation() (0, 2), 0.00001); | |
151 ASSERT_NEAR(-0.57338, camera.GetRotation() (1, 0), 0.00001); | |
152 ASSERT_NEAR(0.22011, camera.GetRotation() (1, 1), 0.00001); | |
153 ASSERT_NEAR(0.78917, camera.GetRotation() (1, 2), 0.00001); | |
154 ASSERT_NEAR(0.70711, camera.GetRotation() (2, 0), 0.00001); | |
155 ASSERT_NEAR(-0.35355, camera.GetRotation() (2, 1), 0.00001); | |
156 ASSERT_NEAR(0.61237, camera.GetRotation() (2, 2), 0.00001); | |
157 | |
158 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
159 | |
160 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), | |
161 camera.GetRotation(), | |
162 camera.GetCenter()); | |
163 | |
164 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); | |
165 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); | |
166 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); | |
167 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); | |
168 } | |
169 | |
170 | |
171 TEST(FiniteProjectiveCamera, Decomposition2) | |
172 { | |
173 const double p[] = { 1188.111986, 580.205341, -808.445330, 128000.000000, -366.466264, 1446.510501, 418.499736, 128000.000000, -0.487118, 0.291726, -0.823172, 500.000000 }; | |
174 const double k[] = { -1528.494743, 0.000000, 256.000000, 0.000000, 1528.494743, 256.000000, 0.000000, 0.000000, 1.000000 }; | |
175 const double r[] = { -0.858893, -0.330733, 0.391047, -0.158171, 0.897503, 0.411668, -0.487118, 0.291726, -0.823172 }; | |
176 const double c[] = { 243.558936, -145.863085, 411.585964 }; | |
177 | |
178 OrthancStone::FiniteProjectiveCamera camera(p); | |
179 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
180 | |
181 OrthancStone::FiniteProjectiveCamera camera2(k, r, c); | |
182 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix(), 1)); | |
183 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters(), 0.001)); | |
184 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation(), 0.000001)); | |
185 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter(), 0.0001)); | |
186 } | |
187 | |
188 | |
189 TEST(FiniteProjectiveCamera, Decomposition3) | |
190 { | |
191 const double p[] = { 10, 0, 0, 0, | |
192 0, 20, 0, 0, | |
193 0, 0, 30, 0 }; | |
194 | |
195 OrthancStone::FiniteProjectiveCamera camera(p); | |
196 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
197 ASSERT_DOUBLE_EQ(10, camera.GetIntrinsicParameters() (0, 0)); | |
198 ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); | |
199 ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); | |
200 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); | |
201 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); | |
202 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); | |
203 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); | |
204 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); | |
205 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); | |
206 } | |
207 | |
208 | |
209 TEST(FiniteProjectiveCamera, Decomposition4) | |
210 { | |
211 const double p[] = { 1, 0, 0, 10, | |
212 0, 1, 0, 20, | |
213 0, 0, 1, 30 }; | |
214 | |
215 OrthancStone::FiniteProjectiveCamera camera(p); | |
216 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
217 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (0, 0)); | |
218 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (1, 1)); | |
219 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (2, 2)); | |
220 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); | |
221 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); | |
222 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); | |
223 ASSERT_DOUBLE_EQ(-10, camera.GetCenter() (0)); | |
224 ASSERT_DOUBLE_EQ(-20, camera.GetCenter() (1)); | |
225 ASSERT_DOUBLE_EQ(-30, camera.GetCenter() (2)); | |
226 } | |
227 | |
228 | |
229 TEST(FiniteProjectiveCamera, Decomposition5) | |
230 { | |
231 const double p[] = { 0, 0, 10, 0, | |
232 0, 20, 0, 0, | |
233 30, 0, 0, 0 }; | |
234 | |
235 OrthancStone::FiniteProjectiveCamera camera(p); | |
236 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
237 ASSERT_DOUBLE_EQ(-10, camera.GetIntrinsicParameters() (0, 0)); | |
238 ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); | |
239 ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); | |
240 ASSERT_DOUBLE_EQ(-1, camera.GetRotation() (0, 2)); | |
241 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); | |
242 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 0)); | |
243 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); | |
244 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); | |
245 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); | |
246 | |
247 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), | |
248 camera.GetRotation(), | |
249 camera.GetCenter()); | |
250 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); | |
251 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); | |
252 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); | |
253 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); | |
254 } | |
255 | |
256 | |
257 static double GetCosAngle(const OrthancStone::Vector& a, | |
258 const OrthancStone::Vector& b) | |
259 { | |
260 // Returns the cosine of the angle between two vectors | |
261 // https://en.wikipedia.org/wiki/Dot_product#Geometric_definition | |
262 return boost::numeric::ublas::inner_prod(a, b) / | |
263 (boost::numeric::ublas::norm_2(a) * boost::numeric::ublas::norm_2(b)); | |
264 } | |
265 | |
266 | |
267 TEST(FiniteProjectiveCamera, Ray) | |
268 { | |
269 const double pp[] = { -1499.650894, 2954.618773, -259.737419, 637891.819097, | |
270 -2951.517707, -1501.019129, -285.785281, 637891.819097, | |
271 0.008528, 0.003067, -0.999959, 2491.764918 }; | |
272 | |
273 const OrthancStone::FiniteProjectiveCamera camera(pp); | |
274 | |
275 ASSERT_NEAR(-21.2492, camera.GetCenter() (0), 0.0001); | |
276 ASSERT_NEAR(-7.64234, camera.GetCenter() (1), 0.00001); | |
277 ASSERT_NEAR(2491.66, camera.GetCenter() (2), 0.01); | |
278 | |
279 // Image plane that led to these parameters, with principal point at | |
280 // (256,256). The image has dimensions 512x512. | |
281 OrthancStone::Vector o = | |
282 OrthancStone::LinearAlgebra::CreateVector(7.009620, 2.521030, -821.942000); | |
283 OrthancStone::Vector ax = | |
284 OrthancStone::LinearAlgebra::CreateVector(-0.453219, 0.891399, -0.001131); | |
285 OrthancStone::Vector ay = | |
286 OrthancStone::LinearAlgebra::CreateVector(-0.891359, -0.453210, -0.008992); | |
287 | |
288 OrthancStone::CoordinateSystem3D imagePlane(o, ax, ay); | |
289 | |
290 // Back-projection of the principal point | |
291 { | |
292 OrthancStone::Vector ray = camera.GetRayDirection(256, 256); | |
293 | |
294 // The principal axis vector is orthogonal to the image plane | |
295 // (i.e. parallel to the plane normal), in the opposite direction | |
296 // ("-1" corresponds to "cos(pi)"). | |
297 ASSERT_NEAR(-1, GetCosAngle(ray, imagePlane.GetNormal()), 0.0000001); | |
298 | |
299 // Forward projection of principal axis, resulting in the principal point | |
300 double x, y; | |
301 camera.ApplyFinite(x, y, camera.GetCenter() - ray); | |
302 | |
303 ASSERT_NEAR(256, x, 0.00001); | |
304 ASSERT_NEAR(256, y, 0.00001); | |
305 } | |
306 | |
307 // Back-projection of the 4 corners of the image | |
308 std::vector<double> cx, cy; | |
309 cx.push_back(0); | |
310 cy.push_back(0); | |
311 cx.push_back(512); | |
312 cy.push_back(0); | |
313 cx.push_back(512); | |
314 cy.push_back(512); | |
315 cx.push_back(0); | |
316 cy.push_back(512); | |
317 | |
318 bool first = true; | |
319 double angle; | |
320 | |
321 for (size_t i = 0; i < cx.size(); i++) | |
322 { | |
323 OrthancStone::Vector ray = camera.GetRayDirection(cx[i], cy[i]); | |
324 | |
325 // Check that the angle wrt. principal axis is the same for all | |
326 // the 4 corners | |
327 double a = GetCosAngle(ray, imagePlane.GetNormal()); | |
328 if (first) | |
329 { | |
330 first = false; | |
331 angle = a; | |
332 } | |
333 else | |
334 { | |
335 ASSERT_NEAR(angle, a, 0.000001); | |
336 } | |
337 | |
338 // Forward projection of the ray, going back to the original point | |
339 double x, y; | |
340 camera.ApplyFinite(x, y, camera.GetCenter() - ray); | |
341 | |
342 ASSERT_NEAR(cx[i], x, 0.00001); | |
343 ASSERT_NEAR(cy[i], y, 0.00001); | |
344 | |
345 // Alternative construction, by computing the intersection of the | |
346 // ray with the image plane | |
347 OrthancStone::Vector p; | |
348 ASSERT_TRUE(imagePlane.IntersectLine(p, camera.GetCenter(), -ray)); | |
349 imagePlane.ProjectPoint(x, y, p); | |
350 ASSERT_NEAR(cx[i], x + 256, 0.01); | |
351 ASSERT_NEAR(cy[i], y + 256, 0.01); | |
352 } | |
353 } | |
354 | |
355 | |
356 TEST(Matrix, Inverse1) | |
357 { | |
358 OrthancStone::Matrix a, b; | |
359 | |
360 a.resize(0, 0); | |
361 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
362 ASSERT_EQ(0u, b.size1()); | |
363 ASSERT_EQ(0u, b.size2()); | |
364 | |
365 a.resize(2, 3); | |
366 ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException); | |
367 | |
368 a.resize(1, 1); | |
369 a(0, 0) = 45.0; | |
370 | |
371 ASSERT_DOUBLE_EQ(45, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); | |
372 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
373 ASSERT_EQ(1u, b.size1()); | |
374 ASSERT_EQ(1u, b.size2()); | |
375 ASSERT_DOUBLE_EQ(1.0 / 45.0, b(0, 0)); | |
376 | |
377 a(0, 0) = 0; | |
378 ASSERT_DOUBLE_EQ(0, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); | |
379 ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException); | |
380 } | |
381 | |
382 | |
383 TEST(Matrix, Inverse2) | |
384 { | |
385 OrthancStone::Matrix a, b; | |
386 a.resize(2, 2); | |
387 a(0, 0) = 4; | |
388 a(0, 1) = 3; | |
389 a(1, 0) = 3; | |
390 a(1, 1) = 2; | |
391 | |
392 ASSERT_DOUBLE_EQ(-1, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); | |
393 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
394 ASSERT_EQ(2u, b.size1()); | |
395 ASSERT_EQ(2u, b.size2()); | |
396 | |
397 ASSERT_DOUBLE_EQ(-2, b(0, 0)); | |
398 ASSERT_DOUBLE_EQ(3, b(0, 1)); | |
399 ASSERT_DOUBLE_EQ(3, b(1, 0)); | |
400 ASSERT_DOUBLE_EQ(-4, b(1, 1)); | |
401 | |
402 a(0, 0) = 1; | |
403 a(0, 1) = 2; | |
404 a(1, 0) = 3; | |
405 a(1, 1) = 4; | |
406 | |
407 ASSERT_DOUBLE_EQ(-2, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); | |
408 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
409 | |
410 ASSERT_DOUBLE_EQ(-2, b(0, 0)); | |
411 ASSERT_DOUBLE_EQ(1, b(0, 1)); | |
412 ASSERT_DOUBLE_EQ(1.5, b(1, 0)); | |
413 ASSERT_DOUBLE_EQ(-0.5, b(1, 1)); | |
414 } | |
415 | |
416 | |
417 TEST(Matrix, Inverse3) | |
418 { | |
419 OrthancStone::Matrix a, b; | |
420 a.resize(3, 3); | |
421 a(0, 0) = 7; | |
422 a(0, 1) = 2; | |
423 a(0, 2) = 1; | |
424 a(1, 0) = 0; | |
425 a(1, 1) = 3; | |
426 a(1, 2) = -1; | |
427 a(2, 0) = -3; | |
428 a(2, 1) = 4; | |
429 a(2, 2) = -2; | |
430 | |
431 ASSERT_DOUBLE_EQ(1, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); | |
432 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
433 ASSERT_EQ(3u, b.size1()); | |
434 ASSERT_EQ(3u, b.size2()); | |
435 | |
436 ASSERT_DOUBLE_EQ(-2, b(0, 0)); | |
437 ASSERT_DOUBLE_EQ(8, b(0, 1)); | |
438 ASSERT_DOUBLE_EQ(-5, b(0, 2)); | |
439 ASSERT_DOUBLE_EQ(3, b(1, 0)); | |
440 ASSERT_DOUBLE_EQ(-11, b(1, 1)); | |
441 ASSERT_DOUBLE_EQ(7, b(1, 2)); | |
442 ASSERT_DOUBLE_EQ(9, b(2, 0)); | |
443 ASSERT_DOUBLE_EQ(-34, b(2, 1)); | |
444 ASSERT_DOUBLE_EQ(21, b(2, 2)); | |
445 | |
446 | |
447 a(0, 0) = 1; | |
448 a(0, 1) = 2; | |
449 a(0, 2) = 2; | |
450 a(1, 0) = 1; | |
451 a(1, 1) = 0; | |
452 a(1, 2) = 1; | |
453 a(2, 0) = 1; | |
454 a(2, 1) = 2; | |
455 a(2, 2) = 1; | |
456 | |
457 ASSERT_DOUBLE_EQ(2, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); | |
458 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
459 ASSERT_EQ(3u, b.size1()); | |
460 ASSERT_EQ(3u, b.size2()); | |
461 | |
462 ASSERT_DOUBLE_EQ(-1, b(0, 0)); | |
463 ASSERT_DOUBLE_EQ(1, b(0, 1)); | |
464 ASSERT_DOUBLE_EQ(1, b(0, 2)); | |
465 ASSERT_DOUBLE_EQ(0, b(1, 0)); | |
466 ASSERT_DOUBLE_EQ(-0.5, b(1, 1)); | |
467 ASSERT_DOUBLE_EQ(0.5, b(1, 2)); | |
468 ASSERT_DOUBLE_EQ(1, b(2, 0)); | |
469 ASSERT_DOUBLE_EQ(0, b(2, 1)); | |
470 ASSERT_DOUBLE_EQ(-1, b(2, 2)); | |
471 } | |
472 | |
473 | |
474 TEST(Matrix, Inverse4) | |
475 { | |
476 OrthancStone::Matrix a, b; | |
477 a.resize(4, 4); | |
478 a(0, 0) = 2; | |
479 a(0, 1) = 1; | |
480 a(0, 2) = 2; | |
481 a(0, 3) = -3; | |
482 a(1, 0) = -2; | |
483 a(1, 1) = 2; | |
484 a(1, 2) = -1; | |
485 a(1, 3) = -1; | |
486 a(2, 0) = 2; | |
487 a(2, 1) = 2; | |
488 a(2, 2) = -3; | |
489 a(2, 3) = -1; | |
490 a(3, 0) = 3; | |
491 a(3, 1) = -2; | |
492 a(3, 2) = -3; | |
493 a(3, 3) = -1; | |
494 | |
495 OrthancStone::LinearAlgebra::InvertMatrix(b, a); | |
496 ASSERT_EQ(4u, b.size1()); | |
497 ASSERT_EQ(4u, b.size2()); | |
498 | |
499 b *= 134.0; // This is the determinant | |
500 | |
501 ASSERT_DOUBLE_EQ(8, b(0, 0)); | |
502 ASSERT_DOUBLE_EQ(-44, b(0, 1)); | |
503 ASSERT_DOUBLE_EQ(30, b(0, 2)); | |
504 ASSERT_DOUBLE_EQ(-10, b(0, 3)); | |
505 ASSERT_DOUBLE_EQ(2, b(1, 0)); | |
506 ASSERT_DOUBLE_EQ(-11, b(1, 1)); | |
507 ASSERT_DOUBLE_EQ(41, b(1, 2)); | |
508 ASSERT_DOUBLE_EQ(-36, b(1, 3)); | |
509 ASSERT_DOUBLE_EQ(16, b(2, 0)); | |
510 ASSERT_DOUBLE_EQ(-21, b(2, 1)); | |
511 ASSERT_DOUBLE_EQ(-7, b(2, 2)); | |
512 ASSERT_DOUBLE_EQ(-20, b(2, 3)); | |
513 ASSERT_DOUBLE_EQ(-28, b(3, 0)); | |
514 ASSERT_DOUBLE_EQ(-47, b(3, 1)); | |
515 ASSERT_DOUBLE_EQ(29, b(3, 2)); | |
516 ASSERT_DOUBLE_EQ(-32, b(3, 3)); | |
517 } | |
518 | |
519 | |
520 TEST(FiniteProjectiveCamera, Calibration) | |
521 { | |
522 unsigned int volumeWidth = 512; | |
523 unsigned int volumeHeight = 512; | |
524 unsigned int volumeDepth = 110; | |
525 | |
526 OrthancStone::Vector camera = OrthancStone::LinearAlgebra::CreateVector | |
527 (-1000, -5000, -static_cast<double>(volumeDepth) * 32); | |
528 | |
529 OrthancStone::Vector principalPoint = OrthancStone::LinearAlgebra::CreateVector | |
530 (volumeWidth/2, volumeHeight/2, volumeDepth * 2); | |
531 | |
532 OrthancStone::FiniteProjectiveCamera c(camera, principalPoint, 0, 512, 512, 1, 1); | |
533 | |
534 double swapv[9] = { 1, 0, 0, | |
535 0, -1, 512, | |
536 0, 0, 1 }; | |
537 OrthancStone::Matrix swap; | |
538 OrthancStone::LinearAlgebra::FillMatrix(swap, 3, 3, swapv); | |
539 | |
540 OrthancStone::Matrix p = OrthancStone::LinearAlgebra::Product(swap, c.GetMatrix()); | |
541 p /= p(2,3); | |
542 | |
543 ASSERT_NEAR( 1.04437, p(0,0), 0.00001); | |
544 ASSERT_NEAR(-0.0703111, p(0,1), 0.00000001); | |
545 ASSERT_NEAR(-0.179283, p(0,2), 0.000001); | |
546 ASSERT_NEAR( 61.7431, p(0,3), 0.0001); | |
547 ASSERT_NEAR( 0.11127, p(1,0), 0.000001); | |
548 ASSERT_NEAR(-0.595541, p(1,1), 0.000001); | |
549 ASSERT_NEAR( 0.872211, p(1,2), 0.000001); | |
550 ASSERT_NEAR( 203.748, p(1,3), 0.001); | |
551 ASSERT_NEAR( 3.08593e-05, p(2,0), 0.0000000001); | |
552 ASSERT_NEAR( 0.000129138, p(2,1), 0.000000001); | |
553 ASSERT_NEAR( 9.18901e-05, p(2,2), 0.0000000001); | |
554 ASSERT_NEAR( 1, p(2,3), 0.0000001); | |
555 } | |
556 | |
557 | |
558 static bool IsEqualRotationVector(OrthancStone::Vector a, | |
559 OrthancStone::Vector b) | |
560 { | |
561 if (a.size() != b.size() || | |
562 a.size() != 3) | |
563 { | |
564 return false; | |
565 } | |
566 else | |
567 { | |
568 OrthancStone::LinearAlgebra::NormalizeVector(a); | |
569 OrthancStone::LinearAlgebra::NormalizeVector(b); | |
570 return OrthancStone::LinearAlgebra::IsCloseToZero(boost::numeric::ublas::norm_2(a - b)); | |
571 } | |
572 } | |
573 | |
574 | |
575 TEST(GeometryToolbox, AlignVectorsWithRotation) | |
576 { | |
577 OrthancStone::Vector a, b; | |
578 OrthancStone::Matrix r; | |
579 | |
580 OrthancStone::LinearAlgebra::AssignVector(a, -200, 200, -846.63); | |
581 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); | |
582 | |
583 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); | |
584 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); | |
585 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); | |
586 | |
587 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, b, a); | |
588 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); | |
589 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, b), a)); | |
590 | |
591 OrthancStone::LinearAlgebra::AssignVector(a, 1, 0, 0); | |
592 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); | |
593 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); | |
594 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); | |
595 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); | |
596 | |
597 OrthancStone::LinearAlgebra::AssignVector(a, 0, 1, 0); | |
598 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); | |
599 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); | |
600 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); | |
601 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); | |
602 | |
603 OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 1); | |
604 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); | |
605 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); | |
606 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); | |
607 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); | |
608 | |
609 OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 0); | |
610 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); | |
611 ASSERT_THROW(OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b), Orthanc::OrthancException); | |
612 | |
613 // TODO: Deal with opposite vectors | |
614 | |
615 /* | |
616 OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, -1); | |
617 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); | |
618 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); | |
619 OrthancStone::LinearAlgebra::Print(r); | |
620 OrthancStone::LinearAlgebra::Print(boost::numeric::ublas::prod(r, a)); | |
621 */ | |
622 } | |
623 | |
624 | |
625 static bool IsEqualVectorL1(OrthancStone::Vector a, | |
626 OrthancStone::Vector b) | |
627 { | |
628 if (a.size() != b.size()) | |
629 { | |
630 return false; | |
631 } | |
632 else | |
633 { | |
634 for (size_t i = 0; i < a.size(); i++) | |
635 { | |
636 if (!OrthancStone::LinearAlgebra::IsNear(a[i], b[i], 0.0001)) | |
637 { | |
638 return false; | |
639 } | |
640 } | |
641 | |
642 return true; | |
643 } | |
644 } | |
645 | |
646 | |
647 TEST(VolumeImageGeometry, Basic) | |
648 { | |
649 using namespace OrthancStone; | |
650 | |
651 VolumeImageGeometry g; | |
652 g.SetSizeInVoxels(10, 20, 30); | |
653 g.SetVoxelDimensions(1, 2, 3); | |
654 | |
655 Vector p = g.GetCoordinates(0, 0, 0); | |
656 ASSERT_EQ(3u, p.size()); | |
657 ASSERT_DOUBLE_EQ(-1.0 / 2.0, p[0]); | |
658 ASSERT_DOUBLE_EQ(-2.0 / 2.0, p[1]); | |
659 ASSERT_DOUBLE_EQ(-3.0 / 2.0, p[2]); | |
660 | |
661 p = g.GetCoordinates(1, 1, 1); | |
662 ASSERT_DOUBLE_EQ(-1.0 / 2.0 + 10.0 * 1.0, p[0]); | |
663 ASSERT_DOUBLE_EQ(-2.0 / 2.0 + 20.0 * 2.0, p[1]); | |
664 ASSERT_DOUBLE_EQ(-3.0 / 2.0 + 30.0 * 3.0, p[2]); | |
665 | |
666 VolumeProjection proj; | |
667 ASSERT_TRUE(g.DetectProjection(proj, g.GetAxialGeometry().GetNormal())); | |
668 ASSERT_EQ(VolumeProjection_Axial, proj); | |
669 ASSERT_TRUE(g.DetectProjection(proj, g.GetCoronalGeometry().GetNormal())); | |
670 ASSERT_EQ(VolumeProjection_Coronal, proj); | |
671 ASSERT_TRUE(g.DetectProjection(proj, g.GetSagittalGeometry().GetNormal())); | |
672 ASSERT_EQ(VolumeProjection_Sagittal, proj); | |
673 | |
674 ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Axial)); | |
675 ASSERT_EQ(20u, g.GetProjectionHeight(VolumeProjection_Axial)); | |
676 ASSERT_EQ(30u, g.GetProjectionDepth(VolumeProjection_Axial)); | |
677 ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Coronal)); | |
678 ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Coronal)); | |
679 ASSERT_EQ(20u, g.GetProjectionDepth(VolumeProjection_Coronal)); | |
680 ASSERT_EQ(20u, g.GetProjectionWidth(VolumeProjection_Sagittal)); | |
681 ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Sagittal)); | |
682 ASSERT_EQ(10u, g.GetProjectionDepth(VolumeProjection_Sagittal)); | |
683 | |
684 p = g.GetVoxelDimensions(VolumeProjection_Axial); | |
685 ASSERT_EQ(3u, p.size()); | |
686 ASSERT_DOUBLE_EQ(1, p[0]); | |
687 ASSERT_DOUBLE_EQ(2, p[1]); | |
688 ASSERT_DOUBLE_EQ(3, p[2]); | |
689 p = g.GetVoxelDimensions(VolumeProjection_Coronal); | |
690 ASSERT_EQ(3u, p.size()); | |
691 ASSERT_DOUBLE_EQ(1, p[0]); | |
692 ASSERT_DOUBLE_EQ(3, p[1]); | |
693 ASSERT_DOUBLE_EQ(2, p[2]); | |
694 p = g.GetVoxelDimensions(VolumeProjection_Sagittal); | |
695 ASSERT_EQ(3u, p.size()); | |
696 ASSERT_DOUBLE_EQ(2, p[0]); | |
697 ASSERT_DOUBLE_EQ(3, p[1]); | |
698 ASSERT_DOUBLE_EQ(1, p[2]); | |
699 | |
700 // Loop over all the voxels of the volume | |
701 for (unsigned int z = 0; z < g.GetDepth(); z++) | |
702 { | |
703 const float zz = (0.5f + static_cast<float>(z)) / static_cast<float>(g.GetDepth()); // Z-center of the voxel | |
704 | |
705 for (unsigned int y = 0; y < g.GetHeight(); y++) | |
706 { | |
707 const float yy = (0.5f + static_cast<float>(y)) / static_cast<float>(g.GetHeight()); // Y-center of the voxel | |
708 | |
709 for (unsigned int x = 0; x < g.GetWidth(); x++) | |
710 { | |
711 const float xx = (0.5f + static_cast<float>(x)) / static_cast<float>(g.GetWidth()); // X-center of the voxel | |
712 | |
713 const float sx = 1.0f; | |
714 const float sy = 2.0f; | |
715 const float sz = 3.0f; | |
716 | |
717 Vector p = g.GetCoordinates(xx, yy, zz); | |
718 | |
719 Vector q = (g.GetAxialGeometry().MapSliceToWorldCoordinates( | |
720 static_cast<double>(x) * sx, | |
721 static_cast<double>(y) * sy) + | |
722 z * sz * g.GetAxialGeometry().GetNormal()); | |
723 ASSERT_TRUE(IsEqualVectorL1(p, q)); | |
724 | |
725 q = (g.GetCoronalGeometry().MapSliceToWorldCoordinates( | |
726 static_cast<double>(x) * sx, | |
727 static_cast<double>(g.GetDepth() - 1 - z) * sz) + | |
728 y * sy * g.GetCoronalGeometry().GetNormal()); | |
729 ASSERT_TRUE(IsEqualVectorL1(p, q)); | |
730 | |
731 /** | |
732 * WARNING: In sagittal geometry, the normal points to | |
733 * REDUCING X-axis in the 3D world. This is necessary to keep | |
734 * the right-hand coordinate system. Hence the "-". | |
735 **/ | |
736 q = (g.GetSagittalGeometry().MapSliceToWorldCoordinates( | |
737 static_cast<double>(y) * sy, | |
738 static_cast<double>(g.GetDepth() - 1 - z) * sz) + | |
739 x * sx * (-g.GetSagittalGeometry().GetNormal())); | |
740 ASSERT_TRUE(IsEqualVectorL1(p, q)); | |
741 } | |
742 } | |
743 } | |
744 | |
745 ASSERT_EQ(0, (int) VolumeProjection_Axial); | |
746 ASSERT_EQ(1, (int) VolumeProjection_Coronal); | |
747 ASSERT_EQ(2, (int) VolumeProjection_Sagittal); | |
748 | |
749 for (int p = 0; p < 3; p++) | |
750 { | |
751 VolumeProjection projection = (VolumeProjection) p; | |
752 const CoordinateSystem3D& s = g.GetProjectionGeometry(projection); | |
753 | |
754 ASSERT_THROW(g.GetProjectionSlice(projection, g.GetProjectionDepth(projection)), Orthanc::OrthancException); | |
755 | |
756 for (unsigned int i = 0; i < g.GetProjectionDepth(projection); i++) | |
757 { | |
758 CoordinateSystem3D plane = g.GetProjectionSlice(projection, i); | |
759 | |
760 if (projection == VolumeProjection_Sagittal) | |
761 { | |
762 ASSERT_TRUE(IsEqualVectorL1(plane.GetOrigin(), s.GetOrigin() + static_cast<double>(i) * | |
763 (-s.GetNormal()) * g.GetVoxelDimensions(projection)[2])); | |
764 } | |
765 else | |
766 { | |
767 ASSERT_TRUE(IsEqualVectorL1(plane.GetOrigin(), s.GetOrigin() + static_cast<double>(i) * | |
768 s.GetNormal() * g.GetVoxelDimensions(projection)[2])); | |
769 } | |
770 | |
771 ASSERT_TRUE(IsEqualVectorL1(plane.GetAxisX(), s.GetAxisX())); | |
772 ASSERT_TRUE(IsEqualVectorL1(plane.GetAxisY(), s.GetAxisY())); | |
773 | |
774 unsigned int slice; | |
775 VolumeProjection q; | |
776 ASSERT_TRUE(g.DetectSlice(q, slice, plane)); | |
777 ASSERT_EQ(projection, q); | |
778 ASSERT_EQ(i, slice); | |
779 } | |
780 } | |
781 } | |
782 | |
783 | |
784 TEST(LinearAlgebra, ParseVectorLocale) | |
785 { | |
786 OrthancStone::Vector v; | |
787 | |
788 ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "1.2")); | |
789 ASSERT_EQ(1u, v.size()); | |
790 ASSERT_DOUBLE_EQ(1.2, v[0]); | |
791 | |
792 ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "-1.2e+2")); | |
793 ASSERT_EQ(1u, v.size()); | |
794 ASSERT_DOUBLE_EQ(-120.0, v[0]); | |
795 | |
796 ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "-1e-2\\2")); | |
797 ASSERT_EQ(2u, v.size()); | |
798 ASSERT_DOUBLE_EQ(-0.01, v[0]); | |
799 ASSERT_DOUBLE_EQ(2.0, v[1]); | |
800 | |
801 ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "1.3671875\\1.3671875")); | |
802 ASSERT_EQ(2u, v.size()); | |
803 ASSERT_DOUBLE_EQ(1.3671875, v[0]); | |
804 ASSERT_DOUBLE_EQ(1.3671875, v[1]); | |
805 } | |
806 | |
807 TEST(GenericToolbox, NormalizeUuid) | |
808 { | |
809 std::string ref("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); | |
810 | |
811 { | |
812 std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); | |
813 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
814 ASSERT_EQ(ref, u); | |
815 } | |
816 | |
817 // space left | |
818 { | |
819 std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); | |
820 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
821 ASSERT_EQ(ref, u); | |
822 } | |
823 | |
824 // space right | |
825 { | |
826 std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb "); | |
827 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
828 ASSERT_EQ(ref, u); | |
829 } | |
830 | |
831 // space l & r | |
832 { | |
833 std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb "); | |
834 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
835 ASSERT_EQ(ref, u); | |
836 } | |
837 | |
838 // space left + case | |
839 { | |
840 std::string u(" 44CA5051-14ef-4d2f-8bd7-dB20bfb61fbb"); | |
841 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
842 ASSERT_EQ(ref, u); | |
843 } | |
844 | |
845 // space right + case | |
846 { | |
847 std::string u("44ca5051-14EF-4D2f-8bd7-db20bfb61fbB "); | |
848 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
849 ASSERT_EQ(ref, u); | |
850 } | |
851 | |
852 // space l & r + case | |
853 { | |
854 std::string u(" 44cA5051-14Ef-4d2f-8bD7-db20bfb61fbb "); | |
855 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
856 ASSERT_EQ(ref, u); | |
857 } | |
858 | |
859 // no | |
860 { | |
861 std::string u(" 44ca5051-14ef-4d2f-8bd7- db20bfb61fbb"); | |
862 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
863 ASSERT_NE(ref, u); | |
864 } | |
865 | |
866 // no | |
867 { | |
868 std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fb"); | |
869 OrthancStone::GenericToolbox::NormalizeUuid(u); | |
870 ASSERT_NE(ref, u); | |
871 } | |
872 } | |
873 | |
874 | |
875 TEST(CoordinateSystem3D, Basic) | |
876 { | |
877 using namespace OrthancStone; | |
878 | |
879 { | |
880 CoordinateSystem3D c; | |
881 ASSERT_FALSE(c.IsValid()); | |
882 ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0); | |
883 ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0); | |
884 ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1); | |
885 | |
886 ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 0))); | |
887 ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(5, 0, 0))); | |
888 ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 5, 0))); | |
889 ASSERT_DOUBLE_EQ(5, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 5))); | |
890 } | |
891 | |
892 { | |
893 CoordinateSystem3D c("nope1", "nope2"); | |
894 ASSERT_FALSE(c.IsValid()); | |
895 ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0); | |
896 ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0); | |
897 ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1); | |
898 } | |
899 | |
900 { | |
901 // https://www.vedantu.com/maths/perpendicular-distance-of-a-point-from-a-plane | |
902 CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, 4, -4, -6); | |
903 ASSERT_DOUBLE_EQ(3, c.ComputeDistance(LinearAlgebra::CreateVector(0, 3, 6))); | |
904 } | |
905 | |
906 { | |
907 // https://mathinsight.org/distance_point_plane_examples | |
908 CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, -2, 5, 8); | |
909 ASSERT_DOUBLE_EQ(39.0 / sqrt(33.0), c.ComputeDistance(LinearAlgebra::CreateVector(4, -4, 3))); | |
910 } | |
911 | |
912 { | |
913 // https://www.ck12.org/calculus/distance-between-a-point-and-a-plane/lesson/Distance-Between-a-Point-and-a-Plane-MAT-ALY/ | |
914 const Vector a = LinearAlgebra::CreateVector(3, 6, 9); | |
915 const Vector b = LinearAlgebra::CreateVector(9, 6, 3); | |
916 const Vector c = LinearAlgebra::CreateVector(6, -9, 9); | |
917 CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c); | |
918 ASSERT_DOUBLE_EQ(0, d.ComputeDistance(a)); | |
919 ASSERT_DOUBLE_EQ(0, d.ComputeDistance(b)); | |
920 ASSERT_DOUBLE_EQ(0, d.ComputeDistance(c)); | |
921 } | |
922 | |
923 { | |
924 // https://tutorial.math.lamar.edu/classes/calcii/eqnsofplanes.aspx | |
925 const Vector a = LinearAlgebra::CreateVector(1, -2, 0); | |
926 const Vector b = LinearAlgebra::CreateVector(3, 1, 4); | |
927 const Vector c = LinearAlgebra::CreateVector(0, -1, 2); | |
928 CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c); | |
929 double r = d.GetNormal() [0] / 2.0; | |
930 ASSERT_DOUBLE_EQ(-8 * r, d.GetNormal() [1]); | |
931 ASSERT_DOUBLE_EQ(5 * r, d.GetNormal() [2]); | |
932 } | |
933 } |