Mercurial > hg > orthanc-stone

comparison OrthancStone/UnitTestsSources/GeometryToolboxTests.cpp @ 1877:a2955abe4c2e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
skeleton for the RenderingPlugin
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 12 Jan 2022 08:23:38 +0100
parents UnitTestsSources/GeometryToolboxTests.cpp@7053b8a0aaec
children cbf54cd28d59
comparison
equal deleted inserted replaced
1876:b1f510e601d2 1877:a2955abe4c2e
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-2022 Osimis S.A., Belgium
6 * Copyright (C) 2021-2022 Sebastien Jodogne, ICTEAM UCLouvain, Belgium
7 *
8 * This program is free software: you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public License
10 * as published by the Free Software Foundation, either version 3 of
11 * the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful, but
14 * WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with this program. If not, see
20 * <http://www.gnu.org/licenses/>.
21 **/
22
23
24 #include <gtest/gtest.h>
25
26 #include "../Sources/Toolbox/DicomInstanceParameters.h"
27 #include "../Sources/Toolbox/FiniteProjectiveCamera.h"
28 #include "../Sources/Toolbox/GenericToolbox.h"
29 #include "../Sources/Toolbox/GeometryToolbox.h"
30 #include "../Sources/Toolbox/SlicesSorter.h"
31
32 #include <Logging.h>
33 #include <OrthancException.h>
34
35 #include <boost/math/special_functions/round.hpp>
36
37
38 TEST(GeometryToolbox, Interpolation)
39 {
40 using namespace OrthancStone::GeometryToolbox;
41
42 // https://en.wikipedia.org/wiki/Bilinear_interpolation#Application_in_image_processing
43 ASSERT_DOUBLE_EQ(146.1f, ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95));
44
45 ASSERT_DOUBLE_EQ(91, ComputeBilinearInterpolationUnitSquare(0, 0, 91, 210, 162, 95));
46 ASSERT_DOUBLE_EQ(210, ComputeBilinearInterpolationUnitSquare(1, 0, 91, 210, 162, 95));
47 ASSERT_DOUBLE_EQ(162, ComputeBilinearInterpolationUnitSquare(0, 1, 91, 210, 162, 95));
48 ASSERT_DOUBLE_EQ(95, ComputeBilinearInterpolationUnitSquare(1, 1, 91, 210, 162, 95));
49
50 ASSERT_DOUBLE_EQ(123.35f, ComputeTrilinearInterpolationUnitSquare
51 (0.5f, 0.2f, 0.7f,
52 91, 210, 162, 95,
53 51, 190, 80, 92));
54
55 ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95),
56 ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 0,
57 91, 210, 162, 95,
58 51, 190, 80, 92));
59
60 ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 51, 190, 80, 92),
61 ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 1,
62 91, 210, 162, 95,
63 51, 190, 80, 92));
64 }
65
66
67 static bool CompareMatrix(const OrthancStone::Matrix& a,
68 const OrthancStone::Matrix& b,
69 double threshold = 0.00000001)
70 {
71 if (a.size1() != b.size1() ||
72 a.size2() != b.size2())
73 {
74 return false;
75 }
76
77 for (size_t i = 0; i < a.size1(); i++)
78 {
79 for (size_t j = 0; j < a.size2(); j++)
80 {
81 if (fabs(a(i, j) - b(i, j)) > threshold)
82 {
83 LOG(ERROR) << "Too large difference in component ("
84 << i << "," << j << "): " << a(i,j) << " != " << b(i,j);
85 return false;
86 }
87 }
88 }
89
90 return true;
91 }
92
93
94 static bool CompareVector(const OrthancStone::Vector& a,
95 const OrthancStone::Vector& b,
96 double threshold = 0.00000001)
97 {
98 if (a.size() != b.size())
99 {
100 return false;
101 }
102
103 for (size_t i = 0; i < a.size(); i++)
104 {
105 if (fabs(a(i) - b(i)) > threshold)
106 {
107 LOG(ERROR) << "Too large difference in component "
108 << i << ": " << a(i) << " != " << b(i);
109 return false;
110 }
111 }
112
113 return true;
114 }
115
116
117
118 TEST(FiniteProjectiveCamera, Decomposition1)
119 {
120 // Example 6.2 of "Multiple View Geometry in Computer Vision - 2nd
121 // edition" (page 163)
122 const double p[12] = {
123 3.53553e+2, 3.39645e+2, 2.77744e+2, -1.44946e+6,
124 -1.03528e+2, 2.33212e+1, 4.59607e+2, -6.32525e+5,
125 7.07107e-1, -3.53553e-1, 6.12372e-1, -9.18559e+2
126 };
127
128 OrthancStone::FiniteProjectiveCamera camera(p);
129 ASSERT_EQ(3u, camera.GetMatrix().size1());
130 ASSERT_EQ(4u, camera.GetMatrix().size2());
131 ASSERT_EQ(3u, camera.GetIntrinsicParameters().size1());
132 ASSERT_EQ(3u, camera.GetIntrinsicParameters().size2());
133 ASSERT_EQ(3u, camera.GetRotation().size1());
134 ASSERT_EQ(3u, camera.GetRotation().size2());
135 ASSERT_EQ(3u, camera.GetCenter().size());
136
137 ASSERT_NEAR(1000.0, camera.GetCenter()[0], 0.01);
138 ASSERT_NEAR(2000.0, camera.GetCenter()[1], 0.01);
139 ASSERT_NEAR(1500.0, camera.GetCenter()[2], 0.01);
140
141 ASSERT_NEAR(468.2, camera.GetIntrinsicParameters() (0, 0), 0.1);
142 ASSERT_NEAR(91.2, camera.GetIntrinsicParameters() (0, 1), 0.1);
143 ASSERT_NEAR(300.0, camera.GetIntrinsicParameters() (0, 2), 0.1);
144 ASSERT_NEAR(427.2, camera.GetIntrinsicParameters() (1, 1), 0.1);
145 ASSERT_NEAR(200.0, camera.GetIntrinsicParameters() (1, 2), 0.1);
146 ASSERT_NEAR(1.0, camera.GetIntrinsicParameters() (2, 2), 0.1);
147
148 ASSERT_NEAR(0, camera.GetIntrinsicParameters() (1, 0), 0.0000001);
149 ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 0), 0.0000001);
150 ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 1), 0.0000001);
151
152 ASSERT_NEAR(0.41380, camera.GetRotation() (0, 0), 0.00001);
153 ASSERT_NEAR(0.90915, camera.GetRotation() (0, 1), 0.00001);
154 ASSERT_NEAR(0.04708, camera.GetRotation() (0, 2), 0.00001);
155 ASSERT_NEAR(-0.57338, camera.GetRotation() (1, 0), 0.00001);
156 ASSERT_NEAR(0.22011, camera.GetRotation() (1, 1), 0.00001);
157 ASSERT_NEAR(0.78917, camera.GetRotation() (1, 2), 0.00001);
158 ASSERT_NEAR(0.70711, camera.GetRotation() (2, 0), 0.00001);
159 ASSERT_NEAR(-0.35355, camera.GetRotation() (2, 1), 0.00001);
160 ASSERT_NEAR(0.61237, camera.GetRotation() (2, 2), 0.00001);
161
162 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation()));
163
164 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(),
165 camera.GetRotation(),
166 camera.GetCenter());
167
168 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix()));
169 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters()));
170 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation()));
171 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter()));
172 }
173
174
175 TEST(FiniteProjectiveCamera, Decomposition2)
176 {
177 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 };
178 const double k[] = { -1528.494743, 0.000000, 256.000000, 0.000000, 1528.494743, 256.000000, 0.000000, 0.000000, 1.000000 };
179 const double r[] = { -0.858893, -0.330733, 0.391047, -0.158171, 0.897503, 0.411668, -0.487118, 0.291726, -0.823172 };
180 const double c[] = { 243.558936, -145.863085, 411.585964 };
181
182 OrthancStone::FiniteProjectiveCamera camera(p);
183 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation()));
184
185 OrthancStone::FiniteProjectiveCamera camera2(k, r, c);
186 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix(), 1));
187 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters(), 0.001));
188 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation(), 0.000001));
189 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter(), 0.0001));
190 }
191
192
193 TEST(FiniteProjectiveCamera, Decomposition3)
194 {
195 const double p[] = { 10, 0, 0, 0,
196 0, 20, 0, 0,
197 0, 0, 30, 0 };
198
199 OrthancStone::FiniteProjectiveCamera camera(p);
200 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation()));
201 ASSERT_DOUBLE_EQ(10, camera.GetIntrinsicParameters() (0, 0));
202 ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1));
203 ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2));
204 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0));
205 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1));
206 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2));
207 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0));
208 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1));
209 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2));
210 }
211
212
213 TEST(FiniteProjectiveCamera, Decomposition4)
214 {
215 const double p[] = { 1, 0, 0, 10,
216 0, 1, 0, 20,
217 0, 0, 1, 30 };
218
219 OrthancStone::FiniteProjectiveCamera camera(p);
220 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation()));
221 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (0, 0));
222 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (1, 1));
223 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (2, 2));
224 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0));
225 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1));
226 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2));
227 ASSERT_DOUBLE_EQ(-10, camera.GetCenter() (0));
228 ASSERT_DOUBLE_EQ(-20, camera.GetCenter() (1));
229 ASSERT_DOUBLE_EQ(-30, camera.GetCenter() (2));
230 }
231
232
233 TEST(FiniteProjectiveCamera, Decomposition5)
234 {
235 const double p[] = { 0, 0, 10, 0,
236 0, 20, 0, 0,
237 30, 0, 0, 0 };
238
239 OrthancStone::FiniteProjectiveCamera camera(p);
240 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation()));
241 ASSERT_DOUBLE_EQ(-10, camera.GetIntrinsicParameters() (0, 0));
242 ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1));
243 ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2));
244 ASSERT_DOUBLE_EQ(-1, camera.GetRotation() (0, 2));
245 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1));
246 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 0));
247 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0));
248 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1));
249 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2));
250
251 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(),
252 camera.GetRotation(),
253 camera.GetCenter());
254 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix()));
255 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters()));
256 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation()));
257 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter()));
258 }
259
260
261 static double GetCosAngle(const OrthancStone::Vector& a,
262 const OrthancStone::Vector& b)
263 {
264 // Returns the cosine of the angle between two vectors
265 // https://en.wikipedia.org/wiki/Dot_product#Geometric_definition
266 return boost::numeric::ublas::inner_prod(a, b) /
267 (boost::numeric::ublas::norm_2(a) * boost::numeric::ublas::norm_2(b));
268 }
269
270
271 TEST(FiniteProjectiveCamera, Ray)
272 {
273 const double pp[] = { -1499.650894, 2954.618773, -259.737419, 637891.819097,
274 -2951.517707, -1501.019129, -285.785281, 637891.819097,
275 0.008528, 0.003067, -0.999959, 2491.764918 };
276
277 const OrthancStone::FiniteProjectiveCamera camera(pp);
278
279 ASSERT_NEAR(-21.2492, camera.GetCenter() (0), 0.0001);
280 ASSERT_NEAR(-7.64234, camera.GetCenter() (1), 0.00001);
281 ASSERT_NEAR(2491.66, camera.GetCenter() (2), 0.01);
282
283 // Image plane that led to these parameters, with principal point at
284 // (256,256). The image has dimensions 512x512.
285 OrthancStone::Vector o =
286 OrthancStone::LinearAlgebra::CreateVector(7.009620, 2.521030, -821.942000);
287 OrthancStone::Vector ax =
288 OrthancStone::LinearAlgebra::CreateVector(-0.453219, 0.891399, -0.001131);
289 OrthancStone::Vector ay =
290 OrthancStone::LinearAlgebra::CreateVector(-0.891359, -0.453210, -0.008992);
291
292 OrthancStone::CoordinateSystem3D imagePlane(o, ax, ay);
293
294 // Back-projection of the principal point
295 {
296 OrthancStone::Vector ray = camera.GetRayDirection(256, 256);
297
298 // The principal axis vector is orthogonal to the image plane
299 // (i.e. parallel to the plane normal), in the opposite direction
300 // ("-1" corresponds to "cos(pi)").
301 ASSERT_NEAR(-1, GetCosAngle(ray, imagePlane.GetNormal()), 0.0000001);
302
303 // Forward projection of principal axis, resulting in the principal point
304 double x, y;
305 camera.ApplyFinite(x, y, camera.GetCenter() - ray);
306
307 ASSERT_NEAR(256, x, 0.00001);
308 ASSERT_NEAR(256, y, 0.00001);
309 }
310
311 // Back-projection of the 4 corners of the image
312 std::vector<double> cx, cy;
313 cx.push_back(0);
314 cy.push_back(0);
315 cx.push_back(512);
316 cy.push_back(0);
317 cx.push_back(512);
318 cy.push_back(512);
319 cx.push_back(0);
320 cy.push_back(512);
321
322 bool first = true;
323 double angle;
324
325 for (size_t i = 0; i < cx.size(); i++)
326 {
327 OrthancStone::Vector ray = camera.GetRayDirection(cx[i], cy[i]);
328
329 // Check that the angle wrt. principal axis is the same for all
330 // the 4 corners
331 double a = GetCosAngle(ray, imagePlane.GetNormal());
332 if (first)
333 {
334 first = false;
335 angle = a;
336 }
337 else
338 {
339 ASSERT_NEAR(angle, a, 0.000001);
340 }
341
342 // Forward projection of the ray, going back to the original point
343 double x, y;
344 camera.ApplyFinite(x, y, camera.GetCenter() - ray);
345
346 ASSERT_NEAR(cx[i], x, 0.00001);
347 ASSERT_NEAR(cy[i], y, 0.00001);
348
349 // Alternative construction, by computing the intersection of the
350 // ray with the image plane
351 OrthancStone::Vector p;
352 ASSERT_TRUE(imagePlane.IntersectLine(p, camera.GetCenter(), -ray));
353 imagePlane.ProjectPoint(x, y, p);
354 ASSERT_NEAR(cx[i], x + 256, 0.01);
355 ASSERT_NEAR(cy[i], y + 256, 0.01);
356 }
357 }
358
359
360 TEST(Matrix, Inverse1)
361 {
362 OrthancStone::Matrix a, b;
363
364 a.resize(0, 0);
365 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
366 ASSERT_EQ(0u, b.size1());
367 ASSERT_EQ(0u, b.size2());
368
369 a.resize(2, 3);
370 ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException);
371
372 a.resize(1, 1);
373 a(0, 0) = 45.0;
374
375 ASSERT_DOUBLE_EQ(45, OrthancStone::LinearAlgebra::ComputeDeterminant(a));
376 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
377 ASSERT_EQ(1u, b.size1());
378 ASSERT_EQ(1u, b.size2());
379 ASSERT_DOUBLE_EQ(1.0 / 45.0, b(0, 0));
380
381 a(0, 0) = 0;
382 ASSERT_DOUBLE_EQ(0, OrthancStone::LinearAlgebra::ComputeDeterminant(a));
383 ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException);
384 }
385
386
387 TEST(Matrix, Inverse2)
388 {
389 OrthancStone::Matrix a, b;
390 a.resize(2, 2);
391 a(0, 0) = 4;
392 a(0, 1) = 3;
393 a(1, 0) = 3;
394 a(1, 1) = 2;
395
396 ASSERT_DOUBLE_EQ(-1, OrthancStone::LinearAlgebra::ComputeDeterminant(a));
397 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
398 ASSERT_EQ(2u, b.size1());
399 ASSERT_EQ(2u, b.size2());
400
401 ASSERT_DOUBLE_EQ(-2, b(0, 0));
402 ASSERT_DOUBLE_EQ(3, b(0, 1));
403 ASSERT_DOUBLE_EQ(3, b(1, 0));
404 ASSERT_DOUBLE_EQ(-4, b(1, 1));
405
406 a(0, 0) = 1;
407 a(0, 1) = 2;
408 a(1, 0) = 3;
409 a(1, 1) = 4;
410
411 ASSERT_DOUBLE_EQ(-2, OrthancStone::LinearAlgebra::ComputeDeterminant(a));
412 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
413
414 ASSERT_DOUBLE_EQ(-2, b(0, 0));
415 ASSERT_DOUBLE_EQ(1, b(0, 1));
416 ASSERT_DOUBLE_EQ(1.5, b(1, 0));
417 ASSERT_DOUBLE_EQ(-0.5, b(1, 1));
418 }
419
420
421 TEST(Matrix, Inverse3)
422 {
423 OrthancStone::Matrix a, b;
424 a.resize(3, 3);
425 a(0, 0) = 7;
426 a(0, 1) = 2;
427 a(0, 2) = 1;
428 a(1, 0) = 0;
429 a(1, 1) = 3;
430 a(1, 2) = -1;
431 a(2, 0) = -3;
432 a(2, 1) = 4;
433 a(2, 2) = -2;
434
435 ASSERT_DOUBLE_EQ(1, OrthancStone::LinearAlgebra::ComputeDeterminant(a));
436 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
437 ASSERT_EQ(3u, b.size1());
438 ASSERT_EQ(3u, b.size2());
439
440 ASSERT_DOUBLE_EQ(-2, b(0, 0));
441 ASSERT_DOUBLE_EQ(8, b(0, 1));
442 ASSERT_DOUBLE_EQ(-5, b(0, 2));
443 ASSERT_DOUBLE_EQ(3, b(1, 0));
444 ASSERT_DOUBLE_EQ(-11, b(1, 1));
445 ASSERT_DOUBLE_EQ(7, b(1, 2));
446 ASSERT_DOUBLE_EQ(9, b(2, 0));
447 ASSERT_DOUBLE_EQ(-34, b(2, 1));
448 ASSERT_DOUBLE_EQ(21, b(2, 2));
449
450
451 a(0, 0) = 1;
452 a(0, 1) = 2;
453 a(0, 2) = 2;
454 a(1, 0) = 1;
455 a(1, 1) = 0;
456 a(1, 2) = 1;
457 a(2, 0) = 1;
458 a(2, 1) = 2;
459 a(2, 2) = 1;
460
461 ASSERT_DOUBLE_EQ(2, OrthancStone::LinearAlgebra::ComputeDeterminant(a));
462 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
463 ASSERT_EQ(3u, b.size1());
464 ASSERT_EQ(3u, b.size2());
465
466 ASSERT_DOUBLE_EQ(-1, b(0, 0));
467 ASSERT_DOUBLE_EQ(1, b(0, 1));
468 ASSERT_DOUBLE_EQ(1, b(0, 2));
469 ASSERT_DOUBLE_EQ(0, b(1, 0));
470 ASSERT_DOUBLE_EQ(-0.5, b(1, 1));
471 ASSERT_DOUBLE_EQ(0.5, b(1, 2));
472 ASSERT_DOUBLE_EQ(1, b(2, 0));
473 ASSERT_DOUBLE_EQ(0, b(2, 1));
474 ASSERT_DOUBLE_EQ(-1, b(2, 2));
475 }
476
477
478 TEST(Matrix, Inverse4)
479 {
480 OrthancStone::Matrix a, b;
481 a.resize(4, 4);
482 a(0, 0) = 2;
483 a(0, 1) = 1;
484 a(0, 2) = 2;
485 a(0, 3) = -3;
486 a(1, 0) = -2;
487 a(1, 1) = 2;
488 a(1, 2) = -1;
489 a(1, 3) = -1;
490 a(2, 0) = 2;
491 a(2, 1) = 2;
492 a(2, 2) = -3;
493 a(2, 3) = -1;
494 a(3, 0) = 3;
495 a(3, 1) = -2;
496 a(3, 2) = -3;
497 a(3, 3) = -1;
498
499 OrthancStone::LinearAlgebra::InvertMatrix(b, a);
500 ASSERT_EQ(4u, b.size1());
501 ASSERT_EQ(4u, b.size2());
502
503 b *= 134.0; // This is the determinant
504
505 ASSERT_DOUBLE_EQ(8, b(0, 0));
506 ASSERT_DOUBLE_EQ(-44, b(0, 1));
507 ASSERT_DOUBLE_EQ(30, b(0, 2));
508 ASSERT_DOUBLE_EQ(-10, b(0, 3));
509 ASSERT_DOUBLE_EQ(2, b(1, 0));
510 ASSERT_DOUBLE_EQ(-11, b(1, 1));
511 ASSERT_DOUBLE_EQ(41, b(1, 2));
512 ASSERT_DOUBLE_EQ(-36, b(1, 3));
513 ASSERT_DOUBLE_EQ(16, b(2, 0));
514 ASSERT_DOUBLE_EQ(-21, b(2, 1));
515 ASSERT_DOUBLE_EQ(-7, b(2, 2));
516 ASSERT_DOUBLE_EQ(-20, b(2, 3));
517 ASSERT_DOUBLE_EQ(-28, b(3, 0));
518 ASSERT_DOUBLE_EQ(-47, b(3, 1));
519 ASSERT_DOUBLE_EQ(29, b(3, 2));
520 ASSERT_DOUBLE_EQ(-32, b(3, 3));
521 }
522
523
524 TEST(FiniteProjectiveCamera, Calibration)
525 {
526 unsigned int volumeWidth = 512;
527 unsigned int volumeHeight = 512;
528 unsigned int volumeDepth = 110;
529
530 OrthancStone::Vector camera = OrthancStone::LinearAlgebra::CreateVector
531 (-1000, -5000, -static_cast<double>(volumeDepth) * 32);
532
533 OrthancStone::Vector principalPoint = OrthancStone::LinearAlgebra::CreateVector
534 (volumeWidth/2, volumeHeight/2, volumeDepth * 2);
535
536 OrthancStone::FiniteProjectiveCamera c(camera, principalPoint, 0, 512, 512, 1, 1);
537
538 double swapv[9] = { 1, 0, 0,
539 0, -1, 512,
540 0, 0, 1 };
541 OrthancStone::Matrix swap;
542 OrthancStone::LinearAlgebra::FillMatrix(swap, 3, 3, swapv);
543
544 OrthancStone::Matrix p = OrthancStone::LinearAlgebra::Product(swap, c.GetMatrix());
545 p /= p(2,3);
546
547 ASSERT_NEAR( 1.04437, p(0,0), 0.00001);
548 ASSERT_NEAR(-0.0703111, p(0,1), 0.00000001);
549 ASSERT_NEAR(-0.179283, p(0,2), 0.000001);
550 ASSERT_NEAR( 61.7431, p(0,3), 0.0001);
551 ASSERT_NEAR( 0.11127, p(1,0), 0.000001);
552 ASSERT_NEAR(-0.595541, p(1,1), 0.000001);
553 ASSERT_NEAR( 0.872211, p(1,2), 0.000001);
554 ASSERT_NEAR( 203.748, p(1,3), 0.001);
555 ASSERT_NEAR( 3.08593e-05, p(2,0), 0.0000000001);
556 ASSERT_NEAR( 0.000129138, p(2,1), 0.000000001);
557 ASSERT_NEAR( 9.18901e-05, p(2,2), 0.0000000001);
558 ASSERT_NEAR( 1, p(2,3), 0.0000001);
559 }
560
561
562 static bool IsEqualRotationVector(OrthancStone::Vector a,
563 OrthancStone::Vector b)
564 {
565 if (a.size() != b.size() ||
566 a.size() != 3)
567 {
568 return false;
569 }
570 else
571 {
572 OrthancStone::LinearAlgebra::NormalizeVector(a);
573 OrthancStone::LinearAlgebra::NormalizeVector(b);
574 return OrthancStone::LinearAlgebra::IsCloseToZero(boost::numeric::ublas::norm_2(a - b));
575 }
576 }
577
578
579 TEST(GeometryToolbox, AlignVectorsWithRotation)
580 {
581 OrthancStone::Vector a, b;
582 OrthancStone::Matrix r;
583
584 OrthancStone::LinearAlgebra::AssignVector(a, -200, 200, -846.63);
585 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1);
586
587 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b);
588 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r));
589 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b));
590
591 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, b, a);
592 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r));
593 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, b), a));
594
595 OrthancStone::LinearAlgebra::AssignVector(a, 1, 0, 0);
596 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1);
597 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b);
598 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r));
599 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b));
600
601 OrthancStone::LinearAlgebra::AssignVector(a, 0, 1, 0);
602 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1);
603 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b);
604 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r));
605 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b));
606
607 OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 1);
608 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1);
609 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b);
610 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r));
611 ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b));
612
613 OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 0);
614 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1);
615 ASSERT_THROW(OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b), Orthanc::OrthancException);
616
617 // TODO: Deal with opposite vectors
618
619 /*
620 OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, -1);
621 OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1);
622 OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b);
623 OrthancStone::LinearAlgebra::Print(r);
624 OrthancStone::LinearAlgebra::Print(boost::numeric::ublas::prod(r, a));
625 */
626 }
627
628
629 static bool IsEqualVectorL1(OrthancStone::Vector a,
630 OrthancStone::Vector b)
631 {
632 if (a.size() != b.size())
633 {
634 return false;
635 }
636 else
637 {
638 for (size_t i = 0; i < a.size(); i++)
639 {
640 if (!OrthancStone::LinearAlgebra::IsNear(a[i], b[i], 0.0001))
641 {
642 return false;
643 }
644 }
645
646 return true;
647 }
648 }
649
650
651 TEST(VolumeImageGeometry, Basic)
652 {
653 using namespace OrthancStone;
654
655 VolumeImageGeometry g;
656 g.SetSizeInVoxels(10, 20, 30);
657 g.SetVoxelDimensions(1, 2, 3);
658
659 Vector p = g.GetCoordinates(0, 0, 0);
660 ASSERT_EQ(3u, p.size());
661 ASSERT_DOUBLE_EQ(-1.0 / 2.0, p[0]);
662 ASSERT_DOUBLE_EQ(-2.0 / 2.0, p[1]);
663 ASSERT_DOUBLE_EQ(-3.0 / 2.0, p[2]);
664
665 p = g.GetCoordinates(1, 1, 1);
666 ASSERT_DOUBLE_EQ(-1.0 / 2.0 + 10.0 * 1.0, p[0]);
667 ASSERT_DOUBLE_EQ(-2.0 / 2.0 + 20.0 * 2.0, p[1]);
668 ASSERT_DOUBLE_EQ(-3.0 / 2.0 + 30.0 * 3.0, p[2]);
669
670 VolumeProjection proj;
671 bool isOpposite;
672 ASSERT_TRUE(g.DetectProjection(proj, isOpposite, g.GetAxialGeometry().GetNormal()));
673 ASSERT_EQ(VolumeProjection_Axial, proj);
674 ASSERT_FALSE(isOpposite);
675 ASSERT_TRUE(g.DetectProjection(proj, isOpposite, g.GetCoronalGeometry().GetNormal()));
676 ASSERT_EQ(VolumeProjection_Coronal, proj);
677 ASSERT_FALSE(isOpposite);
678 ASSERT_TRUE(g.DetectProjection(proj, isOpposite, g.GetSagittalGeometry().GetNormal()));
679 ASSERT_EQ(VolumeProjection_Sagittal, proj);
680 ASSERT_FALSE(isOpposite);
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(GenericToolbox, NormalizeUuid)
816 {
817 std::string ref("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb");
818
819 {
820 std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb");
821 OrthancStone::GenericToolbox::NormalizeUuid(u);
822 ASSERT_EQ(ref, u);
823 }
824
825 // space left
826 {
827 std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb");
828 OrthancStone::GenericToolbox::NormalizeUuid(u);
829 ASSERT_EQ(ref, u);
830 }
831
832 // space right
833 {
834 std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb ");
835 OrthancStone::GenericToolbox::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::GenericToolbox::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::GenericToolbox::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::GenericToolbox::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::GenericToolbox::NormalizeUuid(u);
864 ASSERT_EQ(ref, u);
865 }
866
867 // no
868 {
869 std::string u(" 44ca5051-14ef-4d2f-8bd7- db20bfb61fbb");
870 OrthancStone::GenericToolbox::NormalizeUuid(u);
871 ASSERT_NE(ref, u);
872 }
873
874 // no
875 {
876 std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fb");
877 OrthancStone::GenericToolbox::NormalizeUuid(u);
878 ASSERT_NE(ref, u);
879 }
880 }
881
882
883 TEST(CoordinateSystem3D, Basic)
884 {
885 using namespace OrthancStone;
886
887 {
888 CoordinateSystem3D c;
889 ASSERT_FALSE(c.IsValid());
890 ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0);
891 ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0);
892 ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1);
893
894 ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 0)));
895 ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(5, 0, 0)));
896 ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 5, 0)));
897 ASSERT_DOUBLE_EQ(5, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 5)));
898 }
899
900 {
901 CoordinateSystem3D c("nope1", "nope2");
902 ASSERT_FALSE(c.IsValid());
903 ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0);
904 ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0);
905 ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1);
906 }
907
908 {
909 // https://www.vedantu.com/maths/perpendicular-distance-of-a-point-from-a-plane
910 CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, 4, -4, -6);
911 ASSERT_DOUBLE_EQ(3, c.ComputeDistance(LinearAlgebra::CreateVector(0, 3, 6)));
912 }
913
914 {
915 // https://mathinsight.org/distance_point_plane_examples
916 CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, -2, 5, 8);
917 ASSERT_DOUBLE_EQ(39.0 / sqrt(33.0), c.ComputeDistance(LinearAlgebra::CreateVector(4, -4, 3)));
918 }
919
920 {
921 // https://www.ck12.org/calculus/distance-between-a-point-and-a-plane/lesson/Distance-Between-a-Point-and-a-Plane-MAT-ALY/
922 const Vector a = LinearAlgebra::CreateVector(3, 6, 9);
923 const Vector b = LinearAlgebra::CreateVector(9, 6, 3);
924 const Vector c = LinearAlgebra::CreateVector(6, -9, 9);
925 CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c);
926 ASSERT_DOUBLE_EQ(0, d.ComputeDistance(a));
927 ASSERT_DOUBLE_EQ(0, d.ComputeDistance(b));
928 ASSERT_DOUBLE_EQ(0, d.ComputeDistance(c));
929 }
930
931 {
932 // https://tutorial.math.lamar.edu/classes/calcii/eqnsofplanes.aspx
933 const Vector a = LinearAlgebra::CreateVector(1, -2, 0);
934 const Vector b = LinearAlgebra::CreateVector(3, 1, 4);
935 const Vector c = LinearAlgebra::CreateVector(0, -1, 2);
936 CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c);
937 double r = d.GetNormal() [0] / 2.0;
938 ASSERT_DOUBLE_EQ(-8 * r, d.GetNormal() [1]);
939 ASSERT_DOUBLE_EQ(5 * r, d.GetNormal() [2]);
940 }
941 }
942
943
944 TEST(SlicesSorter, HFP)
945 {
946 // 2021-04-27-repro-bug-HFP-HFS-cartman
947
948 {
949 // This is the last instance in the CT series ("InstanceNumber" is 368):
950 // CT1.2.752.243.1.1.20210202150623868.3730.61448.dcm
951 const OrthancStone::CoordinateSystem3D system("300\\302.5\\323.11", "-1\\0\\0\\0\\-1\\0");
952
953 /**
954 * The first instance in the series ("InstanceNumber" is 1) is
955 * CT1.2.752.243.1.1.20210202150623381.2000.76318.dcm, and its
956 * "ImagePositionPatient" is "300\\302.5\\690.11". It cannot be
957 * taken as the origin of the volume, otherwise the Z axis is
958 * shifted by the depth of the volume.
959 **/
960
961 const double spacingXY = 1.171875;
962 const double width = 512;
963 const double height = 512;
964 const double depth = 368; // Number of instances in the series
965
966 OrthancStone::VolumeImageGeometry geometry;
967 geometry.SetAxialGeometry(system);
968 geometry.SetSizeInVoxels(width, height, depth);
969 geometry.SetVoxelDimensions(spacingXY, spacingXY, 1 /* pixel spacing Z */);
970
971 OrthancStone::Vector p, q;
972
973 p = OrthancStone::LinearAlgebra::CreateVector(0.5 / double(width),
974 0.5 / double(height),
975 0.5 / double(depth), 1);
976 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
977 ASSERT_FLOAT_EQ(300, q[0]);
978 ASSERT_FLOAT_EQ(302.5, q[1]);
979 ASSERT_FLOAT_EQ(323.11, q[2]);
980 ASSERT_FLOAT_EQ(1, q[3]);
981
982 p = OrthancStone::LinearAlgebra::CreateVector((width - 0.5) / double(width),
983 (height - 0.5) / double(height),
984 (depth - 0.5) / double(depth), 1);
985 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
986
987 ASSERT_FLOAT_EQ(300.0 - (width - 1.0) * spacingXY, q[0]); // "X" is swapped
988 ASSERT_FLOAT_EQ(302.5 - (height - 1.0) * spacingXY, q[1]); // "Y" is swapped
989 ASSERT_FLOAT_EQ(690.11, q[2]);
990 ASSERT_FLOAT_EQ(1, q[3]);
991 }
992
993 {
994 // DOSE instance: RD1.2.752.243.1.1.20210202150624529.3790.85357_DoseTPS.dcm
995 OrthancStone::CoordinateSystem3D system("-217.0492\\-161.4141\\376.61", "1\\0\\0\\0\\1\\0");
996 const double spacingXY = 3;
997 const double width = 146; // Columns
998 const double height = 84; // Row
999 const double depth = 86; // Number of frames, same as the length of "GridFrameOffsetVector"
1000
1001 OrthancStone::VolumeImageGeometry geometry;
1002 geometry.SetAxialGeometry(system);
1003 geometry.SetSizeInVoxels(width, height, depth);
1004 geometry.SetVoxelDimensions(spacingXY, spacingXY, 3 /* pixel spacing Z, cf. "GridFrameOffsetVector" */);
1005
1006 OrthancStone::Vector p, q;
1007
1008 p = OrthancStone::LinearAlgebra::CreateVector(0.5 / double(width), 0.5 / double(height), 0.5 / double(depth), 1);
1009 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
1010 ASSERT_FLOAT_EQ(-217.0492, q[0]);
1011 ASSERT_FLOAT_EQ(-161.4141, q[1]);
1012 ASSERT_FLOAT_EQ(376.61, q[2]);
1013 ASSERT_FLOAT_EQ(1, q[3]);
1014
1015 p = OrthancStone::LinearAlgebra::CreateVector((width - 0.5) / double(width),
1016 (height - 0.5) / double(height),
1017 (depth - 0.5) / double(depth), 1);
1018 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
1019
1020 ASSERT_FLOAT_EQ(-217.0492 + (width - 1.0) * spacingXY, q[0]);
1021 ASSERT_FLOAT_EQ(-161.4141 + (height - 1.0) * spacingXY, q[1]);
1022 ASSERT_FLOAT_EQ(376.61 + 255.0 /* last item in "GridFrameOffsetVector" */, q[2]);
1023 ASSERT_FLOAT_EQ(1, q[3]);
1024 }
1025
1026 for (unsigned int upward = 0; upward < 2; upward++)
1027 {
1028 OrthancStone::SlicesSorter slices;
1029
1030 for (unsigned int i = 0; i < 368; i++)
1031 {
1032 unsigned int z = (upward ? 323 + i : 690 - i);
1033 OrthancStone::CoordinateSystem3D p("300\\302.5\\" + boost::lexical_cast<std::string>(z) + ".11",
1034 "-1\\0\\0\\0\\-1\\0");
1035 slices.AddSlice(p, new Orthanc::SingleValueObject<unsigned int>(z));
1036 }
1037
1038 slices.Sort();
1039
1040 double spacingZ;
1041 ASSERT_TRUE(slices.ComputeSpacingBetweenSlices(spacingZ));
1042 ASSERT_FLOAT_EQ(1, spacingZ);
1043 ASSERT_TRUE(slices.AreAllSlicesDistinct());
1044
1045 ASSERT_FLOAT_EQ(300.0, slices.GetSliceGeometry(0).GetOrigin() [0]);
1046 ASSERT_FLOAT_EQ(302.5, slices.GetSliceGeometry(0).GetOrigin() [1]);
1047 ASSERT_FLOAT_EQ(323.11, slices.GetSliceGeometry(0).GetOrigin() [2]);
1048
1049 ASSERT_FLOAT_EQ(300.0, slices.GetSliceGeometry(367).GetOrigin() [0]);
1050 ASSERT_FLOAT_EQ(302.5, slices.GetSliceGeometry(367).GetOrigin() [1]);
1051 ASSERT_FLOAT_EQ(690.11, slices.GetSliceGeometry(367).GetOrigin() [2]);
1052
1053 ASSERT_EQ(323u, dynamic_cast<const Orthanc::SingleValueObject<unsigned int>&>(slices.GetSlicePayload(0)).GetValue());
1054 ASSERT_EQ(690u, dynamic_cast<const Orthanc::SingleValueObject<unsigned int>&>(slices.GetSlicePayload(367)).GetValue());
1055
1056 const unsigned int width = 512;
1057 const unsigned int height = 512;
1058 const unsigned int depth = 368;
1059 const double spacingXY = 1.171875;
1060
1061 OrthancStone::VolumeImageGeometry geometry;
1062 geometry.SetSizeInVoxels(width, height, depth);
1063 geometry.SetVoxelDimensions(spacingXY, spacingXY, 1);
1064 geometry.SetAxialGeometry(slices.GetSliceGeometry(0));
1065
1066 OrthancStone::Vector q = geometry.GetCoordinates(0, 0, 0);
1067 ASSERT_EQ(3u, q.size());
1068 ASSERT_FLOAT_EQ(300 + spacingXY / 2.0, q[0]);
1069 ASSERT_FLOAT_EQ(302.5 + spacingXY / 2.0, q[1]);
1070 ASSERT_FLOAT_EQ(323.11 - 0.5, q[2]);
1071
1072 q = geometry.GetCoordinates(1, 1, 1);
1073 ASSERT_EQ(3u, q.size());
1074 ASSERT_FLOAT_EQ(300 - double(width) * spacingXY + spacingXY / 2.0, q[0]);
1075 ASSERT_FLOAT_EQ(302.5 - double(height) * spacingXY + spacingXY / 2.0, q[1]);
1076 ASSERT_FLOAT_EQ(323.11 + 368.0 - 0.5, q[2]);
1077
1078 OrthancStone::VolumeProjection projection;
1079 unsigned int slice;
1080 ASSERT_TRUE(geometry.DetectSlice(projection, slice, OrthancStone::CoordinateSystem3D("300\\302.5\\690.11", "-1\\0\\0\\0\\-1\\0")));
1081 ASSERT_EQ(OrthancStone::VolumeProjection_Axial, projection);
1082 ASSERT_EQ(367u, slice);
1083 ASSERT_TRUE(geometry.DetectSlice(projection, slice, OrthancStone::CoordinateSystem3D("300\\302.5\\323.11", "-1\\0\\0\\0\\-1\\0")));
1084 ASSERT_EQ(OrthancStone::VolumeProjection_Axial, projection);
1085 ASSERT_EQ(0u, slice);
1086
1087 // DOSE instance: RD1.2.752.243.1.1.20210202150624529.3790.85357_DoseTPS.dcm
1088 ASSERT_TRUE(geometry.DetectSlice(projection, slice,
1089 OrthancStone::CoordinateSystem3D("-217.0492\\-161.4141\\376.61", "1\\0\\0\\0\\1\\0")));
1090 ASSERT_EQ(OrthancStone::VolumeProjection_Axial, projection);
1091 ASSERT_EQ(376u - 323u, slice);
1092 }
1093 }