Line data Source code
1 : /**
2 : * @file test_Intersection.cpp
3 : *
4 : * @brief Test for Intersection class
5 : *
6 : * @author David Frank - initial code
7 : * @author Maximilian Hornung - modularization
8 : * @author Tobias Lasser - consistency changes
9 : */
10 :
11 : #include "doctest/doctest.h"
12 :
13 : #include "Intersection.h"
14 :
15 : using namespace elsa;
16 : using namespace doctest;
17 :
18 : #include "PrettyPrint/Stl.h"
19 : #include "PrettyPrint/Eigen.h"
20 :
21 : TEST_CASE("Intersection: Intersect corners of pixels")
22 1 : {
23 1 : size_t dim = 2;
24 :
25 1 : IndexVector_t voxel(dim);
26 1 : voxel << 1, 0;
27 1 : BoundingBox aabb(voxel);
28 :
29 : // top left corner
30 1 : RealVector_t ro(dim);
31 1 : ro << -3, -3;
32 1 : RealVector_t rd(dim);
33 1 : rd << 1.0, 1.0;
34 1 : rd.normalize();
35 1 : RealRay_t r(ro, rd);
36 :
37 1 : REQUIRE_UNARY(intersectRay(aabb, r));
38 :
39 : // top right corner
40 1 : ro << 1, 2;
41 1 : rd << 1.0, -1.0;
42 1 : rd.normalize();
43 1 : r = RealRay_t(ro, rd);
44 :
45 1 : REQUIRE_UNARY_FALSE(intersectRay(aabb, r));
46 :
47 : // bottom left corner
48 1 : ro << 3, -2;
49 1 : rd << -1.0, 1.0;
50 1 : rd.normalize();
51 1 : r = RealRay_t(ro, rd);
52 :
53 1 : REQUIRE_UNARY(intersectRay(aabb, r));
54 :
55 : // bottom right corner
56 1 : ro << 3, 1;
57 1 : rd << -1.0, -1.0;
58 1 : rd.normalize();
59 1 : r = RealRay_t(ro, rd);
60 :
61 1 : REQUIRE_UNARY_FALSE(intersectRay(aabb, r));
62 1 : }
63 :
64 : TEST_CASE("Intersection: Intersect edges of voxels")
65 9 : {
66 9 : GIVEN("A ray which intersects the edge of a voxel")
67 9 : {
68 9 : size_t dim = 2;
69 9 : RealVector_t ro(dim);
70 9 : ro << 132, 30;
71 9 : RealVector_t rd(dim);
72 9 : rd << -1.0, 0;
73 9 : RealRay_t r(ro, rd);
74 :
75 9 : IndexVector_t voxel(dim);
76 :
77 : // horizontal check
78 9 : voxel << 30, 31;
79 9 : THEN("the ray intersects")
80 9 : {
81 1 : BoundingBox aabb(voxel);
82 1 : REQUIRE_UNARY(intersectRay(aabb, r));
83 1 : }
84 :
85 9 : voxel << 40, 30;
86 9 : THEN("the ray intersects")
87 9 : {
88 1 : BoundingBox aabb(voxel);
89 1 : REQUIRE_UNARY(intersectRay(aabb, r));
90 1 : }
91 :
92 9 : voxel << 40, 29;
93 9 : THEN("the ray does not intersect")
94 9 : {
95 1 : BoundingBox aabb(voxel);
96 1 : REQUIRE_UNARY_FALSE(intersectRay(aabb, r));
97 1 : }
98 :
99 : // vertical check
100 9 : ro << 30, -35;
101 9 : rd << 0, 1.0;
102 9 : r = RealRay_t(ro, rd);
103 :
104 9 : voxel << 31, 30;
105 9 : THEN("the ray intersects")
106 9 : {
107 1 : BoundingBox aabb(voxel);
108 1 : REQUIRE_UNARY(intersectRay(aabb, r));
109 1 : }
110 :
111 9 : voxel << 30, 40;
112 9 : THEN("the ray intersects")
113 9 : {
114 1 : BoundingBox aabb(voxel);
115 1 : REQUIRE_UNARY(intersectRay(aabb, r));
116 1 : }
117 :
118 9 : voxel << 29, 40;
119 9 : THEN("the ray does not intersect")
120 9 : {
121 1 : BoundingBox aabb(voxel);
122 1 : REQUIRE_UNARY_FALSE(intersectRay(aabb, r));
123 1 : }
124 :
125 9 : rd << 0.0, 1.0;
126 9 : ro << 1.5, -10;
127 9 : r = RealRay_t(ro, rd);
128 :
129 9 : voxel << 0, 1;
130 9 : THEN("the ray does not intersect")
131 9 : {
132 1 : BoundingBox aabb(voxel);
133 1 : REQUIRE_UNARY_FALSE(intersectRay(aabb, r));
134 1 : }
135 :
136 9 : voxel << 1, 1;
137 9 : THEN("the ray does not intersect")
138 9 : {
139 1 : BoundingBox aabb(voxel);
140 1 : REQUIRE_UNARY_FALSE(intersectRay(aabb, r));
141 1 : }
142 :
143 9 : voxel << 2, 1;
144 9 : THEN("the ray intersects")
145 9 : {
146 1 : BoundingBox aabb(voxel);
147 1 : REQUIRE_UNARY(intersectRay(aabb, r));
148 1 : }
149 9 : }
150 9 : }
151 :
152 : TEST_SUITE_BEGIN("Intersection::voxelCenter");
153 :
154 : Eigen::IOFormat vecfmt(8, 0, ", ", ", ", "", "", "[", "]");
155 :
156 : TEST_CASE("Intersection: Ray with voxel center")
157 6 : {
158 6 : IndexVector_t upperbound({{5, 5}});
159 6 : BoundingBox aabb(upperbound);
160 :
161 6 : GIVEN("Ray entering the bottom of the volume")
162 6 : {
163 : // Check multiple different rays all coming from the bottom
164 3 : auto steps = 30;
165 96 : for (int i = 0; i <= steps; ++i) {
166 93 : auto inc = (1.f / static_cast<float>(steps)) * static_cast<float>(i);
167 93 : auto dirinc = static_cast<float>(i) / (2 * static_cast<float>(steps));
168 :
169 93 : RealVector_t origin({{0.75f + inc, 0}});
170 93 : RealVector_t dir({{1, 0.5f + dirinc}});
171 93 : dir.normalize();
172 :
173 93 : CAPTURE(inc);
174 93 : CAPTURE(origin.format(vecfmt));
175 93 : CAPTURE(dir.format(vecfmt));
176 :
177 93 : RealRay_t ray(origin, dir);
178 :
179 93 : auto hit = intersectXPlanes(aabb, ray);
180 :
181 93 : CAPTURE(hit);
182 93 : REQUIRE(hit);
183 :
184 93 : THEN("entries x-coord lies in the middle of a voxel")
185 93 : {
186 1 : RealVector_t enter = ray.pointAt(hit->_tmin);
187 1 : CAPTURE(enter.format(vecfmt));
188 :
189 1 : real_t trash = 0;
190 1 : CHECK_EQ(0.5, doctest::Approx(std::modf(enter[0], &trash)));
191 :
192 1 : CHECK((enter.array() < upperbound.template cast<real_t>().array()).all());
193 1 : CHECK((enter.array() > 0).all());
194 1 : }
195 :
196 93 : THEN("x-coord lies in the middle of a voxel")
197 93 : {
198 1 : RealVector_t leave = ray.pointAt(hit->_tmax);
199 1 : CAPTURE(leave.format(vecfmt));
200 :
201 1 : real_t trash = 0;
202 1 : CHECK_EQ(0.5, doctest::Approx(std::modf(leave[0], &trash)));
203 :
204 1 : CHECK((leave.array() < upperbound.template cast<real_t>().array()).all());
205 1 : CHECK((leave.array() > 0).all());
206 1 : }
207 93 : }
208 3 : }
209 :
210 6 : GIVEN("Ray entering the top of the volume")
211 6 : {
212 : // Check multiple different rays all coming from the top
213 3 : auto steps = 30;
214 96 : for (int i = 0; i <= steps; ++i) {
215 93 : auto inc = (1.f / static_cast<float>(steps)) * static_cast<float>(i);
216 93 : auto dirinc = static_cast<float>(i) / (2 * static_cast<float>(steps));
217 :
218 93 : RealVector_t origin({{0.75f + inc, 6}});
219 93 : RealVector_t dir({{1, -0.5f - dirinc}});
220 93 : dir.normalize();
221 :
222 93 : CAPTURE(inc);
223 93 : CAPTURE(origin.format(vecfmt));
224 93 : CAPTURE(dir.format(vecfmt));
225 :
226 93 : RealRay_t ray(origin, dir);
227 :
228 93 : auto hit = intersectXPlanes(aabb, ray);
229 :
230 93 : CAPTURE(hit);
231 93 : REQUIRE(hit);
232 :
233 93 : THEN("entries x-coord lies in the middle of a voxel")
234 93 : {
235 1 : RealVector_t enter = ray.pointAt(hit->_tmin);
236 1 : CAPTURE(enter.format(vecfmt));
237 :
238 1 : real_t trash = 0;
239 1 : CHECK_EQ(0.5, doctest::Approx(std::modf(enter[0], &trash)));
240 :
241 1 : CHECK((enter.array() < upperbound.template cast<real_t>().array()).all());
242 1 : CHECK((enter.array() > 0).all());
243 1 : }
244 :
245 93 : THEN("x-coord lies in the middle of a voxel")
246 93 : {
247 1 : RealVector_t leave = ray.pointAt(hit->_tmax);
248 1 : CAPTURE(leave.format(vecfmt));
249 :
250 1 : real_t trash = 0;
251 1 : CHECK_EQ(0.5, doctest::Approx(std::modf(leave[0], &trash)));
252 :
253 1 : CHECK((leave.array() < upperbound.template cast<real_t>().array()).all());
254 1 : CHECK((leave.array() > 0).all());
255 1 : }
256 93 : }
257 3 : }
258 6 : }
259 :
260 : TEST_CASE("Intersection: Quick bug tests")
261 2 : {
262 2 : const IndexVector_t size({{3, 3}});
263 :
264 2 : BoundingBox aabb(size);
265 2 : aabb.min() = RealVector_t({{-1.5f, -1.5f}});
266 2 : aabb.max() = RealVector_t({{1.5f, 1.5f}});
267 :
268 2 : const RealVector_t ro({{-2, -3}});
269 2 : RealVector_t rd({{1, 1}});
270 2 : rd.normalize();
271 2 : const RealRay_t ray(ro, rd);
272 :
273 2 : auto hit = intersectXPlanes(aabb, ray);
274 :
275 2 : CAPTURE(hit);
276 2 : REQUIRE(hit);
277 :
278 2 : auto dist_to_integer = [&](auto f) {
279 2 : real_t aabbmin = static_cast<int>(std::round(f));
280 2 : auto frac = std::abs(f - aabbmin);
281 2 : return frac;
282 2 : };
283 :
284 2 : THEN("entries x-coord lies in the middle of a voxel")
285 2 : {
286 1 : RealVector_t entry = ray.pointAt(hit->_tmin);
287 1 : CAPTURE(entry.format(vecfmt));
288 :
289 : // real_t trash = 0;
290 1 : CHECK_EQ(0, doctest::Approx(dist_to_integer(entry[0])));
291 :
292 1 : CHECK((entry.array() < 1.5).all());
293 1 : CHECK((entry.array() > -1.5).all());
294 1 : }
295 :
296 2 : THEN("x-coord lies in the middle of a voxel")
297 2 : {
298 1 : RealVector_t leave = ray.pointAt(hit->_tmax);
299 1 : CAPTURE(leave.format(vecfmt));
300 :
301 : // real_t trash = 0;
302 : // CHECK_EQ(0, doctest::Approx(std::modf(leave[0], &trash)));
303 1 : CHECK_EQ(0, doctest::Approx(dist_to_integer(leave[0])));
304 :
305 1 : CHECK((leave.array() < 1.5).all());
306 1 : CHECK((leave.array() > -1.5).all());
307 1 : }
308 2 : }
309 :
310 : TEST_CASE("Intersection: Quick bug tests")
311 2 : {
312 2 : const IndexVector_t size({{3, 3}});
313 :
314 2 : BoundingBox aabb(size);
315 2 : aabb.min() = RealVector_t({{-1.5f, -1.5f}});
316 2 : aabb.max() = RealVector_t({{1.5f, 1.5f}});
317 :
318 2 : const RealVector_t ro({{-2, -3}});
319 2 : RealVector_t rd({{1, 1}});
320 2 : rd.normalize();
321 2 : const RealRay_t ray(ro, rd);
322 :
323 2 : auto hit = intersectXPlanes(aabb, ray);
324 :
325 2 : CAPTURE(hit);
326 2 : REQUIRE(hit);
327 :
328 2 : auto dist_to_integer = [&](auto f) {
329 2 : real_t aabbmin = static_cast<int>(std::round(f));
330 2 : auto frac = std::abs(f - aabbmin);
331 2 : return frac;
332 2 : };
333 :
334 2 : THEN("entries x-coord lies in the middle of a voxel")
335 2 : {
336 1 : RealVector_t entry = ray.pointAt(hit->_tmin);
337 1 : CAPTURE(entry.format(vecfmt));
338 :
339 : // real_t trash = 0;
340 1 : CHECK_EQ(0, doctest::Approx(dist_to_integer(entry[0])));
341 :
342 1 : CHECK((entry.array() < 1.5).all());
343 1 : CHECK((entry.array() > -1.5).all());
344 1 : }
345 :
346 2 : THEN("x-coord lies in the middle of a voxel")
347 2 : {
348 1 : RealVector_t leave = ray.pointAt(hit->_tmax);
349 1 : CAPTURE(leave.format(vecfmt));
350 :
351 : // real_t trash = 0;
352 : // CHECK_EQ(0, doctest::Approx(std::modf(leave[0], &trash)));
353 1 : CHECK_EQ(0, doctest::Approx(dist_to_integer(leave[0])));
354 :
355 1 : CHECK((leave.array() < 1.5).all());
356 1 : CHECK((leave.array() > -1.5).all());
357 1 : }
358 2 : }
359 :
360 : // Redefine GIVEN such that it's nicely usable inside an loop
361 : #undef GIVEN
362 6250000 : #define GIVEN(...) DOCTEST_SUBCASE((std::string(" Given: ") + std::string(__VA_ARGS__)).c_str())
363 :
364 : TEST_CASE("Intersection: 3D rays with bounding box")
365 2500 : {
366 2500 : const IndexVector_t upperbound({{5, 5, 5}});
367 2500 : const BoundingBox aabb(upperbound);
368 :
369 65000 : for (int z = 0; z < 25; ++z) {
370 1625000 : for (int y = 0; y < 25; ++y) {
371 1562500 : const auto zinc = static_cast<real_t>(z) / 5.f;
372 1562500 : const auto yinc = static_cast<real_t>(y) / 5.f;
373 1562500 : GIVEN("ray starting at (6, " + std::to_string(yinc) + ", " + std::to_string(zinc) + ")")
374 1562500 : {
375 1250 : const RealVector_t origin({{-1, yinc, zinc}});
376 1250 : RealVector_t dir({{1, 0, 0}});
377 :
378 1250 : CAPTURE(yinc);
379 1250 : CAPTURE(zinc);
380 1250 : CAPTURE(origin.format(vecfmt));
381 1250 : CAPTURE(dir.format(vecfmt));
382 :
383 1250 : RealRay_t ray(origin, dir);
384 :
385 1250 : auto hit = intersectXPlanes(aabb, ray);
386 :
387 1250 : CAPTURE(hit);
388 1250 : REQUIRE(hit);
389 :
390 1250 : THEN("entries x-coord lies in the middle of a voxel")
391 1250 : {
392 625 : RealVector_t enter = ray.pointAt(hit->_tmin);
393 625 : CAPTURE(enter.format(vecfmt));
394 :
395 625 : real_t trash = 0;
396 625 : CHECK_EQ(0.5, doctest::Approx(std::modf(enter[0], &trash)));
397 :
398 625 : CHECK((enter.array() <= upperbound.template cast<real_t>().array()).all());
399 625 : CHECK((enter.array() >= 0).all());
400 625 : }
401 :
402 1250 : THEN("x-coord lies in the middle of a voxel")
403 1250 : {
404 625 : RealVector_t leave = ray.pointAt(hit->_tmax);
405 625 : CAPTURE(leave.format(vecfmt));
406 :
407 625 : real_t trash = 0;
408 625 : CHECK_EQ(0.5, doctest::Approx(std::modf(leave[0], &trash)));
409 :
410 625 : CHECK((leave.array() <= upperbound.template cast<real_t>().array()).all());
411 625 : CHECK((leave.array() >= 0).all());
412 625 : }
413 1250 : }
414 1562500 : }
415 62500 : }
416 :
417 65000 : for (int z = 0; z < 25; ++z) {
418 1625000 : for (int y = 0; y < 25; ++y) {
419 1562500 : const auto zinc = static_cast<real_t>(z) / 5.f;
420 1562500 : const auto yinc = static_cast<real_t>(y) / 5.f;
421 1562500 : GIVEN("ray starting at (6, " + std::to_string(yinc) + ", " + std::to_string(zinc) + ")")
422 1562500 : {
423 1250 : const RealVector_t origin({{6, yinc, zinc}});
424 1250 : RealVector_t dir({{-1, 0, 0}});
425 :
426 1250 : CAPTURE(yinc);
427 1250 : CAPTURE(zinc);
428 1250 : CAPTURE(origin.format(vecfmt));
429 1250 : CAPTURE(dir.format(vecfmt));
430 :
431 1250 : RealRay_t ray(origin, dir);
432 :
433 1250 : auto hit = intersectXPlanes(aabb, ray);
434 :
435 1250 : CAPTURE(hit);
436 1250 : REQUIRE(hit);
437 :
438 1250 : THEN("entries x-coord lies in the middle of a voxel")
439 1250 : {
440 625 : RealVector_t enter = ray.pointAt(hit->_tmin);
441 625 : CAPTURE(enter.format(vecfmt));
442 :
443 625 : real_t trash = 0;
444 625 : CHECK_EQ(0.5, doctest::Approx(std::modf(enter[0], &trash)));
445 :
446 625 : CHECK((enter.array() <= upperbound.template cast<real_t>().array()).all());
447 625 : CHECK((enter.array() >= 0).all());
448 625 : }
449 :
450 1250 : THEN("x-coord lies in the middle of a voxel")
451 1250 : {
452 625 : RealVector_t leave = ray.pointAt(hit->_tmax);
453 625 : CAPTURE(leave.format(vecfmt));
454 :
455 625 : real_t trash = 0;
456 625 : CHECK_EQ(0.5, doctest::Approx(std::modf(leave[0], &trash)));
457 :
458 625 : CHECK((leave.array() <= upperbound.template cast<real_t>().array()).all());
459 625 : CHECK((leave.array() >= 0).all());
460 625 : }
461 1250 : }
462 1562500 : }
463 62500 : }
464 2500 : }
465 :
466 : TEST_CASE("Intersection: 3D rays with bounding box")
467 2500 : {
468 2500 : const IndexVector_t upperbound({{5, 5, 5}});
469 2500 : BoundingBox aabb(upperbound);
470 2500 : aabb.min() = RealVector_t({{-2.5f, -2.5f, -2.5f}});
471 2500 : aabb.max() = RealVector_t({{2.5f, 2.5f, 2.5f}});
472 :
473 65000 : for (int z = 0; z < 25; ++z) {
474 1625000 : for (int y = 0; y < 25; ++y) {
475 1562500 : const auto zinc = static_cast<real_t>(z) / 5.f - 2.5f;
476 1562500 : const auto yinc = static_cast<real_t>(y) / 5.f - 2.5f;
477 :
478 1562500 : GIVEN("ray starting at (6, " + std::to_string(yinc) + ", " + std::to_string(zinc) + ")")
479 1562500 : {
480 1250 : const RealVector_t origin({{-5, yinc, zinc}});
481 1250 : RealVector_t dir({{1, 0, 0}});
482 :
483 1250 : CAPTURE(yinc);
484 1250 : CAPTURE(zinc);
485 1250 : CAPTURE(origin.format(vecfmt));
486 1250 : CAPTURE(dir.format(vecfmt));
487 :
488 1250 : RealRay_t ray(origin, dir);
489 :
490 1250 : auto hit = intersectXPlanes(aabb, ray);
491 :
492 1250 : CAPTURE(hit);
493 1250 : REQUIRE(hit);
494 :
495 1250 : THEN("entries x-coord lies in the middle of a voxel")
496 1250 : {
497 625 : RealVector_t enter = ray.pointAt(hit->_tmin);
498 625 : CAPTURE(enter.format(vecfmt));
499 :
500 625 : CHECK_EQ(-2, doctest::Approx(enter[0]));
501 :
502 625 : CHECK((enter.array() <= 2.5).all());
503 625 : CHECK((enter.array() >= -2.5).all());
504 625 : }
505 :
506 1250 : THEN("x-coord lies in the middle of a voxel")
507 1250 : {
508 625 : RealVector_t leave = ray.pointAt(hit->_tmax);
509 625 : CAPTURE(leave.format(vecfmt));
510 :
511 625 : CHECK_EQ(2, doctest::Approx(leave[0]));
512 :
513 625 : CHECK((leave.array() <= 2.5).all());
514 625 : CHECK((leave.array() >= -2.5).all());
515 625 : }
516 1250 : }
517 :
518 1562500 : GIVEN("ray starting at (6, " + std::to_string(yinc) + ", " + std::to_string(zinc) + ")")
519 1562500 : {
520 1250 : const RealVector_t origin({{6, yinc, zinc}});
521 1250 : RealVector_t dir({{-1, 0, 0}});
522 :
523 1250 : CAPTURE(yinc);
524 1250 : CAPTURE(zinc);
525 1250 : CAPTURE(origin.format(vecfmt));
526 1250 : CAPTURE(dir.format(vecfmt));
527 :
528 1250 : RealRay_t ray(origin, dir);
529 :
530 1250 : auto hit = intersectXPlanes(aabb, ray);
531 :
532 1250 : CAPTURE(hit);
533 1250 : REQUIRE(hit);
534 :
535 1250 : THEN("entries x-coord lies in the middle of a voxel")
536 1250 : {
537 625 : RealVector_t enter = ray.pointAt(hit->_tmin);
538 625 : CAPTURE(enter.format(vecfmt));
539 :
540 625 : CHECK_EQ(2, doctest::Approx(enter[0]));
541 :
542 625 : CHECK((enter.array() <= 2.5).all());
543 625 : CHECK((enter.array() >= -2.5).all());
544 625 : }
545 :
546 1250 : THEN("x-coord lies in the middle of a voxel")
547 1250 : {
548 625 : RealVector_t leave = ray.pointAt(hit->_tmax);
549 625 : CAPTURE(leave.format(vecfmt));
550 :
551 625 : CHECK_EQ(-2, doctest::Approx(leave[0]));
552 :
553 625 : CHECK((leave.array() <= 2.5).all());
554 625 : CHECK((leave.array() >= -2.5).all());
555 625 : }
556 1250 : }
557 1562500 : }
558 62500 : }
559 2500 : }
560 :
561 : TEST_SUITE_END();
|
