Mercurial > hg > orthanc-stone
comparison OrthancStone/UnitTestsSources/GeometryToolboxTests.cpp @ 1877:a2955abe4c2e
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 } |