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