Line data Source code
1 : /**
2 : * @file test_Geometry.cpp
3 : *
4 : * @brief Test for Geometry class
5 : *
6 : * @author Tobias Lasser - initial code
7 : */
8 :
9 : #include "doctest/doctest.h"
10 : #include "Geometry.h"
11 : #include "VolumeDescriptor.h"
12 :
13 : using namespace elsa;
14 : using namespace doctest;
15 :
16 : TEST_SUITE_BEGIN("core");
17 :
18 16 : TEST_CASE("Geometry: Testing 2D geometries")
19 : {
20 : using namespace geometry;
21 32 : GIVEN("some 2D setup")
22 : {
23 32 : IndexVector_t volCoeff(2);
24 16 : volCoeff << 5, 5;
25 32 : VolumeDescriptor ddVol(volCoeff);
26 :
27 32 : IndexVector_t detCoeff(2);
28 16 : detCoeff << 5, 1;
29 32 : VolumeDescriptor ddDet(detCoeff);
30 :
31 16 : real_t s2c = 10;
32 16 : real_t c2d = 4;
33 :
34 20 : WHEN("testing geometry without rotation/offsets")
35 : {
36 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Radian{0},
37 8 : VolumeData2D{Size2D{volCoeff}}, SinogramData2D{Size2D{detCoeff}});
38 :
39 5 : THEN("a copy is the same")
40 : {
41 2 : Geometry gCopy(g);
42 1 : REQUIRE_EQ(gCopy, g);
43 : }
44 :
45 5 : THEN("then P and Pinv are inverse")
46 : {
47 2 : RealMatrix_t id(2, 2);
48 1 : id.setIdentity();
49 :
50 2 : auto P = g.getProjectionMatrix();
51 2 : auto Pinv = g.getInverseProjectionMatrix();
52 1 : RealMatrix_t result = (P * Pinv) - id;
53 1 : REQUIRE_EQ(result.sum(), Approx(0));
54 : }
55 :
56 5 : THEN("then the rotation matrix is correct")
57 : {
58 1 : RealMatrix_t id(2, 2);
59 1 : id.setIdentity();
60 1 : REQUIRE_EQ((g.getRotationMatrix() - id).sum(), Approx(0));
61 : }
62 :
63 5 : THEN("the camera center is correct")
64 : {
65 2 : auto c = g.getCameraCenter();
66 1 : auto o = ddVol.getLocationOfOrigin();
67 1 : REQUIRE_EQ((c[0] - o[0]), Approx(0));
68 1 : REQUIRE_EQ((c[1] - o[1] + s2c), Approx(0));
69 : }
70 : }
71 :
72 20 : WHEN("testing geometry without rotation but with principal point offset")
73 : {
74 4 : real_t px = 2;
75 :
76 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Radian{0},
77 8 : VolumeData2D{Size2D{volCoeff}}, SinogramData2D{Size2D{detCoeff}},
78 16 : PrincipalPointOffset{px});
79 :
80 5 : THEN("a copy is the same")
81 : {
82 2 : Geometry gCopy(g);
83 1 : REQUIRE_EQ(gCopy, g);
84 : }
85 :
86 5 : THEN("then P and Pinv are inverse")
87 : {
88 2 : RealMatrix_t id(2, 2);
89 1 : id.setIdentity();
90 :
91 2 : auto P = g.getProjectionMatrix();
92 2 : auto Pinv = g.getInverseProjectionMatrix();
93 1 : RealMatrix_t result = (P * Pinv) - id;
94 1 : REQUIRE_EQ(result.sum(), Approx(0));
95 : }
96 :
97 5 : THEN("then the rotation matrix is correct")
98 : {
99 1 : RealMatrix_t id(2, 2);
100 1 : id.setIdentity();
101 1 : REQUIRE_EQ((g.getRotationMatrix() - id).sum(), Approx(0));
102 : }
103 :
104 5 : THEN("the camera center is correct")
105 : {
106 2 : auto c = g.getCameraCenter();
107 1 : auto o = ddVol.getLocationOfOrigin();
108 1 : REQUIRE_EQ((c[0] - o[0]), Approx(0));
109 1 : REQUIRE_EQ((c[1] - o[1] + s2c), Approx(0));
110 : }
111 : }
112 :
113 20 : WHEN("testing geometry with 90 degree rotation and no offsets")
114 : {
115 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{90},
116 8 : VolumeData2D{Size2D{volCoeff}}, SinogramData2D{Size2D{detCoeff}});
117 :
118 5 : THEN("a copy is the same")
119 : {
120 2 : Geometry gCopy(g);
121 1 : REQUIRE_EQ(gCopy, g);
122 : }
123 :
124 5 : THEN("then P and Pinv are inverse")
125 : {
126 2 : RealMatrix_t id(2, 2);
127 1 : id.setIdentity();
128 :
129 2 : auto P = g.getProjectionMatrix();
130 2 : auto Pinv = g.getInverseProjectionMatrix();
131 1 : RealMatrix_t result = (P * Pinv) - id;
132 1 : REQUIRE_EQ(result.sum(), Approx(0.0).epsilon(0.01));
133 : }
134 :
135 5 : THEN("then the rotation matrix is correct")
136 : {
137 1 : RealMatrix_t rot(2, 2);
138 1 : rot << 0, -1, 1, 0;
139 1 : REQUIRE_EQ((g.getRotationMatrix() - rot).sum(), Approx(0).epsilon(0.01));
140 : }
141 :
142 5 : THEN("the camera center is correct")
143 : {
144 2 : auto c = g.getCameraCenter();
145 1 : auto o = ddVol.getLocationOfOrigin();
146 :
147 1 : REQUIRE_EQ((c[0] - o[0] + s2c), Approx(0));
148 1 : REQUIRE_EQ((c[1] - o[1]), Approx(0).epsilon(0.01));
149 : }
150 : }
151 :
152 20 : WHEN("testing geometry with 45 degree rotation and offset center of rotation")
153 : {
154 4 : real_t angle = pi_t / 4; // 45 degrees
155 4 : real_t cx = -1;
156 4 : real_t cy = 2;
157 :
158 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
159 8 : Radian{angle}, VolumeData2D{Size2D{volCoeff}},
160 8 : SinogramData2D{Size2D{detCoeff}}, PrincipalPointOffset{0},
161 16 : RotationOffset2D{cx, cy});
162 :
163 5 : THEN("a copy is the same")
164 : {
165 2 : Geometry gCopy(g);
166 1 : REQUIRE_EQ(gCopy, g);
167 : }
168 :
169 5 : THEN("then P and Pinv are inverse")
170 : {
171 2 : RealMatrix_t id(2, 2);
172 1 : id.setIdentity();
173 :
174 2 : auto P = g.getProjectionMatrix();
175 2 : auto Pinv = g.getInverseProjectionMatrix();
176 1 : RealMatrix_t result = (P * Pinv) - id;
177 1 : REQUIRE_EQ(result.sum(), Approx(0.0).epsilon(0.01));
178 : }
179 :
180 5 : THEN("then the rotation matrix is correct")
181 : {
182 1 : RealMatrix_t rot(2, 2);
183 1 : rot << std::cos(angle), -std::sin(angle), std::sin(angle), std::cos(angle);
184 1 : REQUIRE_EQ((g.getRotationMatrix() - rot).sum(), Approx(0));
185 : }
186 :
187 5 : THEN("the camera center is correct")
188 : {
189 2 : auto c = g.getCameraCenter();
190 1 : auto o = ddVol.getLocationOfOrigin();
191 :
192 1 : real_t oldX = 0;
193 1 : real_t oldY = -s2c;
194 1 : real_t newX = std::cos(angle) * oldX + std::sin(angle) * oldY + o[0] + cx;
195 1 : real_t newY = -std::sin(angle) * oldX + std::cos(angle) * oldY + o[1] + cy;
196 :
197 1 : REQUIRE_EQ((c[0] - newX), Approx(0).epsilon(0.01));
198 1 : REQUIRE_EQ((c[1] - newY), Approx(0).epsilon(0.01));
199 : }
200 : }
201 : }
202 16 : }
203 :
204 20 : TEST_CASE("Geometry: Testing 3D geometries")
205 : {
206 : using namespace geometry;
207 40 : GIVEN("some 3D setup")
208 : {
209 40 : IndexVector_t volCoeff(3);
210 20 : volCoeff << 5, 5, 5;
211 40 : VolumeDescriptor ddVol(volCoeff);
212 :
213 40 : IndexVector_t detCoeff(3);
214 20 : detCoeff << 5, 5, 1;
215 40 : VolumeDescriptor ddDet(detCoeff);
216 :
217 20 : real_t s2c = 10;
218 20 : real_t c2d = 4;
219 :
220 24 : WHEN("testing geometry without rotation/offsets")
221 : {
222 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
223 8 : VolumeData3D{Size3D{volCoeff}}, SinogramData3D{Size3D{detCoeff}},
224 16 : RotationAngles3D{Radian{0}, Radian{0}, Radian{0}});
225 :
226 5 : THEN("a copy is the same")
227 : {
228 2 : Geometry gCopy(g);
229 1 : REQUIRE_EQ(gCopy, g);
230 : }
231 :
232 5 : THEN("then P and Pinv are inverse")
233 : {
234 2 : RealMatrix_t id(3, 3);
235 1 : id.setIdentity();
236 :
237 2 : auto P = g.getProjectionMatrix();
238 2 : auto Pinv = g.getInverseProjectionMatrix();
239 1 : RealMatrix_t result = (P * Pinv) - id;
240 1 : REQUIRE_EQ(result.sum(), Approx(0));
241 : }
242 :
243 5 : THEN("then the rotation matrix is correct")
244 : {
245 1 : RealMatrix_t id(3, 3);
246 1 : id.setIdentity();
247 1 : REQUIRE_EQ((g.getRotationMatrix() - id).sum(), Approx(0));
248 : }
249 :
250 5 : THEN("the camera center is correct")
251 : {
252 2 : auto c = g.getCameraCenter();
253 1 : auto o = ddVol.getLocationOfOrigin();
254 1 : REQUIRE_EQ((c[0] - o[0]), Approx(0));
255 1 : REQUIRE_EQ((c[1] - o[1]), Approx(0));
256 1 : REQUIRE_EQ((c[2] - o[2] + s2c), Approx(0));
257 : }
258 : }
259 :
260 24 : WHEN("testing geometry without rotation but with principal point offsets")
261 : {
262 4 : real_t px = -1;
263 4 : real_t py = 3;
264 :
265 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
266 8 : VolumeData3D{Size3D{volCoeff}}, SinogramData3D{Size3D{detCoeff}},
267 4 : RotationAngles3D{Radian{0}, Radian{0}, Radian{0}},
268 20 : PrincipalPointOffset2D{px, py});
269 :
270 5 : THEN("a copy is the same")
271 : {
272 2 : Geometry gCopy(g);
273 1 : REQUIRE_EQ(gCopy, g);
274 : }
275 :
276 5 : THEN("then P and Pinv are inverse")
277 : {
278 2 : RealMatrix_t id(3, 3);
279 1 : id.setIdentity();
280 :
281 2 : auto P = g.getProjectionMatrix();
282 2 : auto Pinv = g.getInverseProjectionMatrix();
283 1 : RealMatrix_t result = (P * Pinv) - id;
284 1 : REQUIRE_EQ(result.sum(), Approx(0));
285 : }
286 :
287 5 : THEN("then the rotation matrix is correct")
288 : {
289 1 : RealMatrix_t id(3, 3);
290 1 : id.setIdentity();
291 1 : REQUIRE_EQ((g.getRotationMatrix() - id).sum(), Approx(0));
292 : }
293 :
294 5 : THEN("the camera center is correct")
295 : {
296 2 : auto c = g.getCameraCenter();
297 1 : auto o = ddVol.getLocationOfOrigin();
298 1 : REQUIRE_EQ((c[0] - o[0]), Approx(0));
299 1 : REQUIRE_EQ((c[1] - o[1]), Approx(0));
300 1 : REQUIRE_EQ((c[2] - o[2] + s2c), Approx(0));
301 : }
302 : }
303 :
304 24 : WHEN("testing geometry with 90 degree rotation and no offsets")
305 : {
306 4 : real_t angle = pi_t / 2;
307 :
308 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
309 8 : VolumeData3D{Size3D{volCoeff}}, SinogramData3D{Size3D{detCoeff}},
310 16 : RotationAngles3D{Radian{angle}, Radian{0}, Radian{0}});
311 :
312 5 : THEN("a copy is the same")
313 : {
314 2 : Geometry gCopy(g);
315 1 : REQUIRE_EQ(gCopy, g);
316 : }
317 :
318 5 : THEN("then P and Pinv are inverse")
319 : {
320 2 : RealMatrix_t id(3, 3);
321 1 : id.setIdentity();
322 :
323 2 : auto P = g.getProjectionMatrix();
324 2 : auto Pinv = g.getInverseProjectionMatrix();
325 1 : RealMatrix_t result = (P * Pinv) - id;
326 1 : REQUIRE_EQ(result.sum(), Approx(0).epsilon(0.01));
327 : }
328 :
329 5 : THEN("then the rotation matrix is correct")
330 : {
331 1 : RealMatrix_t rot(3, 3);
332 2 : rot << std::cos(angle), 0, std::sin(angle), 0, 1, 0, -std::sin(angle), 0,
333 2 : std::cos(angle);
334 1 : REQUIRE_EQ((g.getRotationMatrix() - rot).sum(), Approx(0));
335 : }
336 :
337 5 : THEN("the camera center is correct")
338 : {
339 2 : auto c = g.getCameraCenter();
340 1 : auto o = ddVol.getLocationOfOrigin();
341 :
342 1 : REQUIRE_EQ((c[0] - o[0] - s2c), Approx(0));
343 1 : REQUIRE_EQ((c[1] - o[1]), Approx(0).epsilon(0.01));
344 1 : REQUIRE_EQ((c[2] - o[2]), Approx(0).epsilon(0.01));
345 : }
346 : }
347 :
348 24 : WHEN("testing geometry with 45/22.5 degree rotation and offset center of rotation")
349 : {
350 4 : real_t angle1 = pi_t / 4;
351 4 : real_t angle2 = pi_t / 2;
352 8 : RealVector_t offset(3);
353 4 : offset << 1, -2, -1;
354 :
355 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
356 8 : VolumeData3D{Size3D{volCoeff}}, SinogramData3D{Size3D{detCoeff}},
357 4 : RotationAngles3D{Radian{angle1}, Radian{angle2}, Radian{0}},
358 20 : PrincipalPointOffset2D{0, 0}, RotationOffset3D{offset});
359 :
360 5 : THEN("a copy is the same")
361 : {
362 2 : Geometry gCopy(g);
363 1 : REQUIRE_EQ(gCopy, g);
364 : }
365 :
366 5 : THEN("then P and Pinv are inverse")
367 : {
368 2 : RealMatrix_t id(3, 3);
369 1 : id.setIdentity();
370 :
371 2 : auto P = g.getProjectionMatrix();
372 2 : auto Pinv = g.getInverseProjectionMatrix();
373 1 : RealMatrix_t result = (P * Pinv) - id;
374 1 : REQUIRE_EQ(result.sum(), Approx(0).epsilon(0.01));
375 : }
376 :
377 5 : THEN("then the rotation matrix is correct")
378 : {
379 2 : RealMatrix_t rot1(3, 3);
380 2 : rot1 << std::cos(angle1), 0, std::sin(angle1), 0, 1, 0, -std::sin(angle1), 0,
381 2 : std::cos(angle1);
382 1 : RealMatrix_t rot2(3, 3);
383 2 : rot2 << std::cos(angle2), -std::sin(angle2), 0, std::sin(angle2), std::cos(angle2),
384 2 : 0, 0, 0, 1;
385 1 : REQUIRE_EQ((g.getRotationMatrix() - rot1 * rot2).sum(), Approx(0));
386 : }
387 :
388 5 : THEN("the camera center is correct")
389 : {
390 2 : auto c = g.getCameraCenter();
391 2 : auto o = ddVol.getLocationOfOrigin();
392 :
393 2 : RealVector_t src(3);
394 1 : src[0] = 0;
395 1 : src[1] = 0;
396 1 : src[2] = -s2c;
397 1 : RealVector_t rotSrc = g.getRotationMatrix().transpose() * src + o + offset;
398 :
399 1 : REQUIRE_EQ((rotSrc - c).sum(), Approx(0).epsilon(0.01));
400 : }
401 : }
402 :
403 24 : WHEN("testing geometry with 45/22.5/12.25 degree rotation as a rotation matrix")
404 : {
405 4 : real_t angle1 = pi_t / 4;
406 4 : real_t angle2 = pi_t / 2;
407 4 : real_t angle3 = pi_t / 8;
408 8 : RealMatrix_t rot1(3, 3);
409 8 : rot1 << std::cos(angle1), 0, std::sin(angle1), 0, 1, 0, -std::sin(angle1), 0,
410 8 : std::cos(angle1);
411 8 : RealMatrix_t rot2(3, 3);
412 8 : rot2 << std::cos(angle2), -std::sin(angle2), 0, std::sin(angle2), std::cos(angle2), 0,
413 8 : 0, 0, 1;
414 8 : RealMatrix_t rot3(3, 3);
415 8 : rot3 << std::cos(angle3), 0, std::sin(angle3), 0, 1, 0, -std::sin(angle3), 0,
416 8 : std::cos(angle3);
417 8 : RealMatrix_t R = rot1 * rot2 * rot3;
418 :
419 8 : Geometry g(s2c, c2d, ddVol, ddDet, R);
420 :
421 5 : THEN("a copy is the same")
422 : {
423 2 : Geometry gCopy(g);
424 1 : REQUIRE_EQ(gCopy, g);
425 : }
426 :
427 5 : THEN("then P and Pinv are inverse")
428 : {
429 2 : RealMatrix_t id(3, 3);
430 1 : id.setIdentity();
431 :
432 2 : auto P = g.getProjectionMatrix();
433 2 : auto Pinv = g.getInverseProjectionMatrix();
434 1 : RealMatrix_t result = (P * Pinv) - id;
435 1 : REQUIRE_EQ(result.sum(), Approx(0).epsilon(0.01));
436 : }
437 :
438 5 : THEN("then the rotation matrix is correct")
439 : {
440 1 : REQUIRE_EQ((g.getRotationMatrix() - R).sum(), Approx(0));
441 : }
442 :
443 5 : THEN("the camera center is correct")
444 : {
445 2 : auto c = g.getCameraCenter();
446 2 : auto o = ddVol.getLocationOfOrigin();
447 :
448 2 : RealVector_t src(3);
449 1 : src[0] = 0;
450 1 : src[1] = 0;
451 1 : src[2] = -s2c;
452 1 : RealVector_t rotSrc = g.getRotationMatrix().transpose() * src + o;
453 :
454 1 : REQUIRE_EQ((rotSrc - c).sum(), Approx(0).epsilon(0.01));
455 : }
456 : }
457 : }
458 20 : }
459 :
460 : TEST_SUITE_END();
|