Mercurial > hg > orthanc-stone

comparison OrthancStone/UnitTestsSources/UnitTestsMain.cpp @ 1512:244ad1e4e76a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 }