Line data Source code
1 : /**
2 : * @file test_TraverseAABB.cpp
3 : *
4 : * @brief Test for TraverseAABB class
5 : *
6 : * @author David Frank - initial code
7 : * @author Maximilian Hornung - modularization, fixes
8 : * @author Tobias Lasser - minor fixes
9 : */
10 :
11 : #include "doctest/doctest.h"
12 :
13 : #include "TraverseAABBBranchless.h"
14 : #include "Intersection.h"
15 :
16 : using namespace elsa;
17 : using namespace doctest;
18 :
19 : bool intersect(const RealVector_t& voxel, const RealRay_t& r)
20 3255 : {
21 : // pre-check parallel rays
22 8997 : for (index_t i = 0; i < r.dim(); ++i) {
23 6254 : real_t tmp = std::abs(r.origin()(i) - voxel(i));
24 :
25 6254 : if (std::abs(r.direction()(i)) < 0.0000001 && tmp >= 0.0 && tmp < 1.0)
26 512 : return true;
27 6254 : }
28 :
29 : // check if ray intersects pixel
30 3255 : IndexVector_t ones(voxel.size());
31 2743 : ones.setOnes();
32 2743 : BoundingBox bb(ones);
33 2743 : bb.min() += voxel;
34 2743 : bb.max() += voxel;
35 :
36 2743 : return intersectRay(bb, r).operator bool();
37 3255 : }
38 :
39 : TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
40 28 : {
41 : // setup
42 28 : const size_t dim = 2;
43 28 : index_t x = 3;
44 28 : index_t y = 3;
45 28 : IndexVector_t volumeDims(dim);
46 28 : volumeDims << x, y;
47 28 : IndexVector_t productOfCoefficientsPerDimension(2);
48 28 : productOfCoefficientsPerDimension << 1, 3;
49 :
50 28 : RealVector_t spacing(dim);
51 :
52 28 : RealVector_t ro(dim);
53 28 : RealVector_t rd(dim);
54 :
55 28 : GIVEN("A 3x3 aabb with standard spacing")
56 28 : {
57 28 : BoundingBox aabb(volumeDims);
58 :
59 : //================================
60 : // intersection from straight rays from the bottom
61 : //================================
62 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (0.5, "
63 28 : "-0.5) and direction (0, 1)")
64 28 : {
65 1 : ro << 0.5, -0.5;
66 1 : rd << 0.0, 1.0;
67 1 : RealRay_t r(ro, rd);
68 :
69 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
70 1 : CHECK_UNARY(traverse.isInBoundingBox());
71 1 : THEN("The ray intersects the aabb at the bottom left pixel")
72 1 : {
73 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
74 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
75 1 : }
76 1 : }
77 :
78 28 : WHEN("A ray with origin = (1.0, -0.5) and direction (0, 1), hits the boundary between 2 "
79 28 : "voxels")
80 28 : {
81 1 : ro << 1.0, -0.5;
82 1 : rd << 0.0, 1.0;
83 1 : RealRay_t r(ro, rd);
84 :
85 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
86 1 : CHECK_UNARY(traverse.isInBoundingBox());
87 1 : THEN("The ray intersects the aabb at the bottom left pixel")
88 1 : {
89 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 1);
90 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
91 1 : }
92 1 : }
93 :
94 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (1.5, "
95 28 : "-0.5) and direction (0, 1)")
96 28 : {
97 1 : ro << 1.5, -0.5;
98 1 : rd << 0.0, 1.0;
99 1 : RealRay_t r(ro, rd);
100 :
101 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
102 1 : CHECK_UNARY(traverse.isInBoundingBox());
103 1 : THEN("The ray intersects the aabb at the bottom left pixel")
104 1 : {
105 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 1);
106 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
107 1 : }
108 1 : }
109 :
110 28 : WHEN("A ray with origin = (2.0, -0.5) and direction (0, 1), hits the boundary between 2 "
111 28 : "voxels")
112 28 : {
113 1 : ro << 2.0, -0.5;
114 1 : rd << 0.0, 1.0;
115 1 : RealRay_t r(ro, rd);
116 :
117 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
118 1 : CHECK_UNARY(traverse.isInBoundingBox());
119 1 : THEN("The ray intersects the aabb at the bottom left pixel")
120 1 : {
121 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
122 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
123 1 : }
124 1 : }
125 :
126 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (2.5, "
127 28 : "-0.5) and direction (0, 1)")
128 28 : {
129 1 : ro << 2.5, -0.5;
130 1 : rd << 0.0, 1.0;
131 1 : RealRay_t r(ro, rd);
132 :
133 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
134 1 : CHECK_UNARY(traverse.isInBoundingBox());
135 1 : THEN("The ray intersects the aabb at the bottom left pixel")
136 1 : {
137 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
138 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
139 1 : }
140 1 : }
141 :
142 : //================================
143 : // intersection from straight rays from the left
144 : //================================
145 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (-0.5, "
146 28 : "0.5) and direction (1, 0)")
147 28 : {
148 1 : ro << -0.5, 0.5;
149 1 : rd << 1.0, 0.0;
150 1 : RealRay_t r(ro, rd);
151 :
152 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
153 1 : CHECK_UNARY(traverse.isInBoundingBox());
154 1 : THEN("The ray intersects the aabb at the bottom left pixel")
155 1 : {
156 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
157 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
158 1 : }
159 1 : }
160 :
161 28 : WHEN("A ray with origin = (1.0, -0.5) and direction (0, 1), hits the boundary between 2 "
162 28 : "voxels")
163 28 : {
164 1 : ro << -0.5, 1.0;
165 1 : rd << 1.0, 0.0;
166 1 : RealRay_t r(ro, rd);
167 :
168 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
169 1 : CHECK_UNARY(traverse.isInBoundingBox());
170 1 : THEN("The ray intersects the aabb at the bottom left pixel")
171 1 : {
172 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
173 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 1);
174 1 : }
175 1 : }
176 :
177 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (-0.5, "
178 28 : "1.5) and direction (1, 0)")
179 28 : {
180 1 : ro << -0.5, 1.5;
181 1 : rd << 1.0, 0.0;
182 1 : RealRay_t r(ro, rd);
183 :
184 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
185 1 : CHECK_UNARY(traverse.isInBoundingBox());
186 1 : THEN("The ray intersects the aabb at the bottom left pixel")
187 1 : {
188 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
189 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 1);
190 1 : }
191 1 : }
192 :
193 28 : WHEN("A ray with origin = (2.0, -0.5) and direction (0, 1), hits the boundary between 2 "
194 28 : "voxels")
195 28 : {
196 1 : ro << -0.5, 2.0;
197 1 : rd << 1.0, 0.0;
198 1 : RealRay_t r(ro, rd);
199 :
200 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
201 1 : CHECK_UNARY(traverse.isInBoundingBox());
202 1 : THEN("The ray intersects the aabb at the bottom left pixel")
203 1 : {
204 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
205 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
206 1 : }
207 1 : }
208 :
209 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (-0.5, "
210 28 : "2.5) and direction (1, 0)")
211 28 : {
212 1 : ro << -0.5, 2.5;
213 1 : rd << 1.0, 0.0;
214 1 : RealRay_t r(ro, rd);
215 :
216 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
217 1 : CHECK_UNARY(traverse.isInBoundingBox());
218 1 : THEN("The ray intersects the aabb at the bottom left pixel")
219 1 : {
220 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
221 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
222 1 : }
223 1 : }
224 :
225 : //================================
226 : // intersection from straight rays from the right
227 : //================================
228 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (3.5, "
229 28 : "0.5) and direction (-1, 0)")
230 28 : {
231 1 : ro << 3.5, 0.5;
232 1 : rd << -1.0, 0.0;
233 1 : RealRay_t r(ro, rd);
234 :
235 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
236 1 : CHECK_UNARY(traverse.isInBoundingBox());
237 1 : THEN("The ray intersects the aabb at the bottom left pixel")
238 1 : {
239 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
240 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
241 1 : }
242 1 : }
243 :
244 28 : WHEN("A ray with origin = (3.5, 1.0) and direction (-1, 0), hits the boundary between 2 "
245 28 : "voxels")
246 28 : {
247 1 : ro << 3.5, 1.0;
248 1 : rd << -1.0, 0.0;
249 1 : RealRay_t r(ro, rd);
250 :
251 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
252 1 : CHECK_UNARY(traverse.isInBoundingBox());
253 1 : THEN("The ray intersects the aabb at the bottom left pixel")
254 1 : {
255 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
256 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 1);
257 1 : }
258 1 : }
259 :
260 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (3.5, "
261 28 : "1.5) and direction (-1, 0)")
262 28 : {
263 1 : ro << 3.5, 1.5;
264 1 : rd << -1.0, 0.0;
265 1 : RealRay_t r(ro, rd);
266 :
267 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
268 1 : CHECK_UNARY(traverse.isInBoundingBox());
269 1 : THEN("The ray intersects the aabb at the bottom left pixel")
270 1 : {
271 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
272 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 1);
273 1 : }
274 1 : }
275 :
276 28 : WHEN("A ray with origin = (3.5, 2.0) and direction (-1, 0), hits the boundary between 2 "
277 28 : "voxels")
278 28 : {
279 1 : ro << 3.5, 2.0;
280 1 : rd << -1.0, 0.0;
281 1 : RealRay_t r(ro, rd);
282 :
283 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
284 1 : CHECK_UNARY(traverse.isInBoundingBox());
285 1 : THEN("The ray intersects the aabb at the bottom left pixel")
286 1 : {
287 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
288 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
289 1 : }
290 1 : }
291 :
292 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (3.5, "
293 28 : "2.5) and direction (-1, 0)")
294 28 : {
295 1 : ro << 3.5, 2.5;
296 1 : rd << -1.0, 0.0;
297 1 : RealRay_t r(ro, rd);
298 :
299 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
300 1 : CHECK_UNARY(traverse.isInBoundingBox());
301 1 : THEN("The ray intersects the aabb at the bottom left pixel")
302 1 : {
303 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
304 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
305 1 : }
306 1 : }
307 :
308 : //================================
309 : // intersection from straight rays from the top
310 : //================================
311 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (0.5, "
312 28 : "3.5) and direction (0, -1)")
313 28 : {
314 1 : ro << 0.5, 3.5;
315 1 : rd << 0.0, -1.0;
316 1 : RealRay_t r(ro, rd);
317 :
318 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
319 1 : CHECK_UNARY(traverse.isInBoundingBox());
320 1 : THEN("The ray intersects the aabb at the bottom left pixel")
321 1 : {
322 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
323 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
324 1 : }
325 1 : }
326 :
327 28 : WHEN("A ray with origin = (1.0, 3.5) and direction (-1, 0), hits the boundary between 2 "
328 28 : "voxels")
329 28 : {
330 1 : ro << 1.0, 3.5;
331 1 : rd << 0.0, -1.0;
332 1 : RealRay_t r(ro, rd);
333 :
334 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
335 1 : CHECK_UNARY(traverse.isInBoundingBox());
336 1 : THEN("The ray intersects the aabb at the bottom left pixel")
337 1 : {
338 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 1);
339 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
340 1 : }
341 1 : }
342 :
343 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (1.5, "
344 28 : "3.5) and direction (0, -1)")
345 28 : {
346 1 : ro << 1.5, 3.5;
347 1 : rd << 0.0, -1.0;
348 1 : RealRay_t r(ro, rd);
349 :
350 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
351 1 : CHECK_UNARY(traverse.isInBoundingBox());
352 1 : THEN("The ray intersects the aabb at the bottom left pixel")
353 1 : {
354 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 1);
355 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
356 1 : }
357 1 : }
358 :
359 28 : WHEN("A ray with origin = (2.0, 3.5) and direction (-1, 0), hits the boundary between 2 "
360 28 : "voxels")
361 28 : {
362 1 : ro << 2.0, 3.5;
363 1 : rd << 0.0, -1.0;
364 1 : RealRay_t r(ro, rd);
365 :
366 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
367 1 : CHECK_UNARY(traverse.isInBoundingBox());
368 1 : THEN("The ray intersects the aabb at the bottom left pixel")
369 1 : {
370 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
371 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
372 1 : }
373 1 : }
374 :
375 28 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (2.5, "
376 28 : "3.5) and direction (0, -1)")
377 28 : {
378 1 : ro << 2.5, 3.5;
379 1 : rd << 0.0, -1.0;
380 1 : RealRay_t r(ro, rd);
381 :
382 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
383 1 : CHECK_UNARY(traverse.isInBoundingBox());
384 1 : THEN("The ray intersects the aabb at the bottom left pixel")
385 1 : {
386 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
387 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
388 1 : }
389 1 : }
390 :
391 : //
392 : // Some edge cases
393 : //
394 28 : WHEN("A ray with origin = (-0.5, 0.0) and direction (1, 0) hits the left edge of aabb")
395 28 : {
396 1 : ro << -0.5, 0.0;
397 1 : rd << 1.0, 0.0;
398 1 : RealRay_t r(ro, rd);
399 :
400 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
401 1 : CHECK_UNARY(traverse.isInBoundingBox());
402 1 : THEN("The ray intersects the aabb at the bottom left pixel")
403 1 : {
404 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
405 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
406 1 : }
407 1 : }
408 :
409 28 : WHEN("A ray with origin = (3.5, 0.0) and direction (-1, 0) hits the left edge of aabb")
410 28 : {
411 1 : ro << 3.5, 0.0;
412 1 : rd << -1.0, 0.0;
413 1 : RealRay_t r(ro, rd);
414 :
415 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
416 1 : CHECK_UNARY(traverse.isInBoundingBox());
417 1 : THEN("The ray intersects the aabb at the top left pixel")
418 1 : {
419 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 2);
420 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
421 1 : }
422 1 : }
423 :
424 28 : WHEN("A ray with origin = (-0.5, 3.0) and direction (1, 0) hits the top edge of aabb")
425 28 : {
426 1 : ro << -0.5, 3.0;
427 1 : rd << 1.0, 0.0;
428 1 : RealRay_t r(ro, rd);
429 :
430 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
431 :
432 1 : THEN("the the aabb is not hit")
433 1 : {
434 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
435 1 : }
436 1 : }
437 :
438 28 : WHEN("A ray with origin = (3.5, 3.0) and direction (-1, 0) hits the top edge of aabb")
439 28 : {
440 1 : ro << 3.5, 3.0;
441 1 : rd << -1.0, 0.0;
442 1 : RealRay_t r(ro, rd);
443 :
444 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
445 :
446 1 : THEN("the the aabb is not hit")
447 1 : {
448 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
449 1 : }
450 1 : }
451 :
452 28 : WHEN("A ray with origin = (0.0, -0.5) and direction (0, 1) hits the bottom edge of aabb)")
453 28 : {
454 1 : ro << 0.0, -0.5;
455 1 : rd << 0.0, 1.0;
456 1 : RealRay_t r(ro, rd);
457 :
458 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
459 1 : CHECK_UNARY(traverse.isInBoundingBox());
460 1 : THEN("The ray intersects the aabb at the bottom left pixel")
461 1 : {
462 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
463 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
464 1 : }
465 1 : }
466 :
467 28 : WHEN("A ray with origin = (0.0, 3.5) and direction (0, -1) hits the top edge of aabb)")
468 28 : {
469 1 : ro << 0.0, 3.5;
470 1 : rd << 0.0, -1.0;
471 1 : RealRay_t r(ro, rd);
472 :
473 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
474 1 : CHECK_UNARY(traverse.isInBoundingBox());
475 1 : THEN("The ray intersects the aabb at the top left pixel")
476 1 : {
477 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
478 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 2);
479 1 : }
480 1 : }
481 :
482 28 : WHEN("A ray with origin = (3.0, -0.5) and direction (0, 1) hits the right edge of aabb)")
483 28 : {
484 1 : ro << 3.0, -0.5;
485 1 : rd << 0.0, 1.0;
486 1 : RealRay_t r(ro, rd);
487 :
488 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
489 :
490 1 : THEN("the the aabb is not hit")
491 1 : {
492 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
493 1 : }
494 1 : }
495 :
496 28 : WHEN("A ray with origin = (3.0, 3.5) and direction (0, -1) hits the top edge of aabb)")
497 28 : {
498 1 : ro << 3.0, 3.5;
499 1 : rd << 0.0, -1.0;
500 1 : RealRay_t r(ro, rd);
501 :
502 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
503 :
504 1 : THEN("the the aabb is not hit")
505 1 : {
506 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
507 1 : }
508 1 : }
509 28 : }
510 28 : }
511 :
512 : TEST_CASE("TraverseAABB: Construction of a 3D traversal object")
513 1 : {
514 : // setup
515 1 : const size_t dim = 3;
516 1 : IndexVector_t volumeDims(dim);
517 1 : volumeDims << 3, 3, 3;
518 1 : IndexVector_t productOfCoefficientsPerDimension(3);
519 1 : productOfCoefficientsPerDimension << 1, 3, 9;
520 :
521 1 : RealVector_t spacing(dim);
522 :
523 1 : RealVector_t ro(dim);
524 1 : RealVector_t rd(dim);
525 :
526 1 : GIVEN("a 3x3x3 aabb with standard spacing")
527 1 : {
528 1 : BoundingBox aabb(volumeDims);
529 :
530 1 : WHEN("A traversal algorithms is initialised with the aabb and a ray with origin = (0.5, "
531 1 : "-0.5, 0.5) and direction = (0, 1, 0)")
532 1 : {
533 1 : ro << 0.5, -0.5, 0.5;
534 1 : rd << 0.0, 1.0, 0.0;
535 1 : RealRay_t r(ro, rd);
536 :
537 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
538 1 : CHECK_UNARY(traverse.isInBoundingBox());
539 :
540 1 : THEN("The ray intersects the aabb at the voxel (0, 0, 0)")
541 1 : {
542 1 : REQUIRE_UNARY(traverse.isInBoundingBox());
543 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
544 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
545 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(2), 0);
546 1 : }
547 1 : }
548 1 : }
549 1 : }
550 :
551 : TEST_CASE("TraverseAABB: Traverse a minimal 3D volume of size 1x1x1")
552 1 : {
553 : // setup
554 1 : const size_t dim = 3;
555 1 : index_t x = 1;
556 1 : index_t y = 1;
557 1 : index_t z = 1;
558 1 : IndexVector_t volumeDims(dim);
559 1 : volumeDims << x, y, z;
560 1 : IndexVector_t productOfCoefficientsPerDimension(3);
561 1 : productOfCoefficientsPerDimension << 1, 1, 1;
562 :
563 1 : RealVector_t spacing(dim);
564 1 : RealVector_t ro(dim);
565 1 : RealVector_t rd(dim);
566 :
567 1 : GIVEN("A 1x1x1 volume with uniform scaling")
568 1 : {
569 1 : BoundingBox aabb(volumeDims);
570 1 : spacing << 1.0, 1.0, 1.0;
571 :
572 1 : WHEN("The volume is traversed with a ray with origin = (-0.5, 0.5, 0.5) and a direction = "
573 1 : "(0, 1, 0)")
574 1 : {
575 1 : ro << 0.5, -0.5, 0.5;
576 1 : rd << 0.0, 1.0, 0.0;
577 :
578 1 : RealRay_t r(ro, rd);
579 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
580 1 : CHECK_UNARY(traverse.isInBoundingBox());
581 :
582 1 : traverse.updateTraverse();
583 :
584 1 : THEN("The algorithms left the volume and the voxel it left the box is (0, 1, 0)")
585 1 : {
586 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
587 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
588 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 1);
589 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(2), 0);
590 1 : }
591 1 : }
592 1 : }
593 1 : }
594 :
595 : TEST_CASE("TraverseAABB: Traverse a 2D volume and only check that the endpoint is correct")
596 1 : {
597 : // setup
598 1 : const size_t dim = 2;
599 1 : index_t x = 10;
600 1 : index_t y = 10;
601 1 : IndexVector_t volumeDims(dim);
602 1 : volumeDims << x, y;
603 1 : IndexVector_t productOfCoefficientsPerDimension(2);
604 1 : productOfCoefficientsPerDimension << 1, 10;
605 :
606 1 : RealVector_t spacing(dim);
607 1 : RealVector_t ro(dim);
608 1 : RealVector_t rd(dim);
609 :
610 1 : GIVEN("A 10x10 volume with uniform scaling")
611 1 : {
612 1 : BoundingBox aabb(volumeDims);
613 :
614 1 : WHEN("The volume is traversed with a ray with origin = (-0.5, 0.5, 0.5) and a direction = "
615 1 : "(0, 1, 0)")
616 1 : {
617 1 : ro << -1, 4.5;
618 1 : rd << 1.0, 0;
619 :
620 1 : RealRay_t r(ro, rd);
621 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
622 1 : CHECK_UNARY(traverse.isInBoundingBox());
623 :
624 11 : while (traverse.isInBoundingBox())
625 10 : traverse.updateTraverse();
626 :
627 1 : THEN("The endpoint should be (10,4)")
628 1 : {
629 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
630 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 10);
631 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 4);
632 1 : }
633 1 : }
634 1 : }
635 1 : }
636 :
637 : TEST_CASE("TraverseAABB: Traverse a 3D Volume diagonally")
638 2 : {
639 : // TODO: run through all 4 diagonals
640 : // TODO: make a non cube volume and run through all 4 diagonals
641 : // TODO: make non uniform scaling and run through all 4 diagonals
642 :
643 2 : const size_t dim = 3;
644 2 : IndexVector_t volumeDims(dim);
645 2 : volumeDims << 10, 10, 10;
646 2 : IndexVector_t productOfCoefficientsPerDimension(3);
647 2 : productOfCoefficientsPerDimension << 1, 10, 100;
648 :
649 2 : RealVector_t spacing(dim);
650 2 : RealVector_t ro(dim);
651 2 : RealVector_t rd(dim);
652 :
653 2 : GIVEN("A 10x10 volume with uniform scaling")
654 2 : {
655 2 : BoundingBox aabb(volumeDims);
656 2 : WHEN("Start at (-1, -1, -1) (so bottom left front) and run to (10, 10, 10) (so top right "
657 2 : "back)")
658 2 : {
659 2 : ro << -1.0, -1.0, -1.0;
660 2 : rd << 1.0, 1.0, 1.0;
661 2 : rd.normalize();
662 :
663 2 : RealRay_t r(ro, rd);
664 :
665 2 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
666 2 : CHECK_UNARY(traverse.isInBoundingBox());
667 :
668 2 : THEN("You entered at (0, 0, 0)")
669 2 : {
670 1 : CHECK_EQ(traverse.getCurrentVoxel()(0), 0);
671 1 : CHECK_EQ(traverse.getCurrentVoxel()(1), 0);
672 1 : CHECK_EQ(traverse.getCurrentVoxel()(2), 0);
673 1 : }
674 22 : while (traverse.isInBoundingBox())
675 20 : traverse.updateTraverse();
676 :
677 2 : THEN("You leave the volume at (10, 10, 10)")
678 2 : {
679 1 : REQUIRE_UNARY_FALSE(traverse.isInBoundingBox());
680 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 10);
681 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 10);
682 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(2), 10);
683 1 : }
684 2 : }
685 2 : }
686 2 : }
687 :
688 : TEST_CASE("TraverseAABB: Check that the first step into the 2D Volume is correct")
689 4 : {
690 4 : const size_t dim = 2;
691 4 : index_t x = 5;
692 4 : index_t y = 5;
693 4 : IndexVector_t volumeDims(dim);
694 4 : volumeDims << x, y;
695 4 : IndexVector_t productOfCoefficientsPerDimension(2);
696 4 : productOfCoefficientsPerDimension << 1, 5;
697 :
698 4 : RealVector_t ro(dim);
699 4 : RealVector_t rd(dim);
700 :
701 4 : GIVEN("A 5x5 volume with uniform scaling")
702 4 : {
703 4 : BoundingBox aabb(volumeDims);
704 :
705 4 : WHEN("The ray direction has the biggest value on the y axis")
706 4 : {
707 2 : ro << 0, 0;
708 2 : rd << 0.5f, 0.7f;
709 2 : rd.normalize();
710 :
711 2 : RealRay_t r(ro, rd);
712 2 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
713 2 : CHECK_UNARY(traverse.isInBoundingBox());
714 :
715 2 : THEN("The traversal is initially at (0, 0)")
716 2 : {
717 1 : auto voxel = traverse.getCurrentVoxel();
718 1 : CHECK_EQ(voxel(0), 0);
719 1 : CHECK_EQ(voxel(1), 0);
720 1 : }
721 :
722 2 : traverse.updateTraverse();
723 :
724 2 : THEN("The first step is in y direction")
725 2 : {
726 1 : REQUIRE_UNARY(traverse.isInBoundingBox());
727 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 0);
728 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 1);
729 1 : }
730 2 : }
731 :
732 4 : WHEN("The ray direction has the biggest value on the x axis")
733 4 : {
734 2 : ro << 0, 0;
735 2 : rd << 0.7f, 0.5f;
736 2 : rd.normalize();
737 :
738 2 : RealRay_t r(ro, rd);
739 2 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
740 2 : CHECK_UNARY(traverse.isInBoundingBox());
741 :
742 2 : THEN("The traversal is initially at (0, 0)")
743 2 : {
744 1 : auto voxel = traverse.getCurrentVoxel();
745 1 : CHECK_EQ(voxel(0), 0);
746 1 : CHECK_EQ(voxel(1), 0);
747 1 : }
748 :
749 2 : traverse.updateTraverse();
750 :
751 2 : THEN("The first step is in y direction")
752 2 : {
753 1 : REQUIRE_UNARY(traverse.isInBoundingBox());
754 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(0), 1);
755 1 : REQUIRE_EQ(traverse.getCurrentVoxel()(1), 0);
756 1 : }
757 2 : }
758 4 : }
759 4 : }
760 :
761 : TEST_CASE("TraverseAABB: Traverse_Volume_2D_EachPointIsTested")
762 1 : {
763 : // setup
764 1 : const size_t dim = 2;
765 1 : index_t x = 128;
766 1 : index_t y = 128;
767 1 : IndexVector_t volumeDims(dim);
768 1 : volumeDims << x, y;
769 1 : BoundingBox aabb(volumeDims);
770 1 : IndexVector_t productOfCoefficientsPerDimension(2);
771 1 : productOfCoefficientsPerDimension << 1, 128;
772 :
773 1 : RealVector_t ro(dim);
774 1 : ro << -168.274f, -143.397f;
775 1 : RealVector_t rd(dim);
776 1 : rd << 0.761124909f, 0.648605406f;
777 1 : rd.normalize();
778 1 : RealRay_t r(ro, rd);
779 :
780 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
781 1 : CHECK_UNARY(traverse.isInBoundingBox());
782 :
783 1 : size_t iter = 0;
784 238 : while (traverse.isInBoundingBox()) {
785 237 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
786 237 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
787 :
788 237 : REQUIRE_UNARY(intersect(voxel, r));
789 237 : traverse.updateTraverse();
790 237 : iter++;
791 237 : }
792 1 : }
793 :
794 : TEST_CASE("TraverseAABB: Traversal through 2D volume should be equal to a ray voxel intersection "
795 : "for every voxel along the way")
796 8 : {
797 : // TODO make this a stronger test, for first some "easy" direction (parallel ones)
798 : // TODO Then make some harder ones
799 : // setup
800 8 : const size_t dim = 2;
801 8 : index_t x = 128;
802 8 : index_t y = 128;
803 8 : IndexVector_t volumeDims(dim);
804 8 : volumeDims << x, y;
805 8 : BoundingBox aabb(volumeDims);
806 8 : IndexVector_t productOfCoefficientsPerDimension(2);
807 8 : productOfCoefficientsPerDimension << 1, 128;
808 :
809 8 : RealVector_t ro(dim);
810 8 : RealVector_t rd(dim);
811 :
812 8 : GIVEN("a point at the bottom left of the volume and a ray with leading dimension x")
813 8 : {
814 1 : ro << -168.274f, -143.397f;
815 :
816 1 : rd << 0.761124909f, 0.648605406f;
817 1 : rd.normalize();
818 :
819 1 : RealRay_t r(ro, rd);
820 :
821 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
822 1 : {
823 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
824 1 : CHECK_UNARY(traverse.isInBoundingBox());
825 :
826 1 : size_t iter = 0;
827 238 : while (traverse.isInBoundingBox()) {
828 237 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
829 237 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
830 :
831 237 : REQUIRE_UNARY(intersect(voxel, r));
832 237 : traverse.updateTraverse();
833 237 : iter++;
834 237 : }
835 1 : }
836 1 : }
837 :
838 8 : GIVEN("a point at the bottom left of the volume and a ray with leading dimension y")
839 8 : {
840 1 : ro << 0, 0;
841 :
842 1 : rd << 0.648605406f, 0.761124909f;
843 1 : rd.normalize();
844 :
845 1 : RealRay_t r(ro, rd);
846 :
847 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
848 1 : {
849 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
850 1 : CHECK_UNARY(traverse.isInBoundingBox());
851 :
852 1 : size_t iter = 0;
853 238 : while (traverse.isInBoundingBox()) {
854 237 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
855 237 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
856 :
857 237 : REQUIRE_UNARY(intersect(voxel, r));
858 237 : traverse.updateTraverse();
859 237 : iter++;
860 237 : }
861 1 : }
862 1 : }
863 :
864 8 : GIVEN("a ray going through the border of voxel column 0 and 1")
865 8 : {
866 1 : ro << 1, -0.5f;
867 :
868 1 : rd << 0, 1;
869 1 : rd.normalize();
870 :
871 1 : RealRay_t r(ro, rd);
872 :
873 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
874 1 : {
875 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
876 1 : CHECK_UNARY(traverse.isInBoundingBox());
877 :
878 1 : size_t iter = 0;
879 129 : while (traverse.isInBoundingBox()) {
880 128 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
881 128 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
882 :
883 128 : REQUIRE_UNARY(intersect(voxel, r));
884 128 : traverse.updateTraverse();
885 128 : iter++;
886 128 : }
887 1 : }
888 1 : }
889 :
890 8 : GIVEN("a ray going through the border of voxel row 0 and 1")
891 8 : {
892 1 : ro << -0.5f, 1;
893 :
894 1 : rd << 1, 0;
895 1 : rd.normalize();
896 :
897 1 : RealRay_t r(ro, rd);
898 :
899 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
900 1 : {
901 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
902 1 : CHECK_UNARY(traverse.isInBoundingBox());
903 :
904 1 : size_t iter = 0;
905 129 : while (traverse.isInBoundingBox()) {
906 128 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
907 128 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
908 :
909 128 : REQUIRE_UNARY(intersect(voxel, r));
910 128 : traverse.updateTraverse();
911 128 : iter++;
912 128 : }
913 1 : }
914 1 : }
915 :
916 8 : GIVEN("A ray going diagonally through the volume")
917 8 : {
918 1 : ro << -0.5f, -0.5f;
919 :
920 1 : rd << 1, 1;
921 1 : rd.normalize();
922 :
923 1 : RealRay_t r(ro, rd);
924 :
925 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
926 1 : {
927 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
928 1 : CHECK_UNARY(traverse.isInBoundingBox());
929 :
930 1 : size_t iter = 0;
931 129 : while (traverse.isInBoundingBox()) {
932 128 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
933 128 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
934 :
935 128 : REQUIRE_UNARY(intersect(voxel, r));
936 128 : traverse.updateTraverse();
937 128 : iter++;
938 128 : }
939 1 : }
940 1 : }
941 :
942 8 : GIVEN("A ray going diagonally through the volume")
943 8 : {
944 1 : ro << -0.5f, 32;
945 :
946 1 : rd << 1, 1;
947 1 : rd.normalize();
948 :
949 1 : RealRay_t r(ro, rd);
950 :
951 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
952 1 : {
953 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
954 1 : CHECK_UNARY(traverse.isInBoundingBox());
955 :
956 1 : size_t iter = 0;
957 192 : while (traverse.isInBoundingBox()) {
958 191 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
959 191 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
960 :
961 191 : REQUIRE_UNARY(intersect(voxel, r));
962 191 : traverse.updateTraverse();
963 191 : iter++;
964 191 : }
965 1 : }
966 1 : }
967 :
968 8 : GIVEN("A ray going diagonally through the volume")
969 8 : {
970 1 : ro << -0.5f, -0.5f;
971 :
972 1 : rd << 0.699428f, 0.472203f;
973 1 : rd.normalize();
974 :
975 1 : RealRay_t r(ro, rd);
976 :
977 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
978 1 : {
979 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
980 1 : CHECK_UNARY(traverse.isInBoundingBox());
981 :
982 1 : size_t iter = 0;
983 215 : while (traverse.isInBoundingBox()) {
984 214 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
985 214 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
986 :
987 214 : REQUIRE_UNARY(intersect(voxel, r));
988 214 : traverse.updateTraverse();
989 214 : iter++;
990 214 : }
991 1 : }
992 1 : }
993 :
994 8 : GIVEN("A ray going diagonally through the volume")
995 8 : {
996 1 : volumeDims << 64, 64;
997 1 : BoundingBox aabb(volumeDims);
998 :
999 1 : ro << 32.0002823f, 3232;
1000 :
1001 1 : rd << -0.000723466626f, -0.999999762f;
1002 1 : rd.normalize();
1003 :
1004 1 : RealRay_t r(ro, rd);
1005 :
1006 1 : THEN("Then all points the traversal visits are also hit by the intersection algorithm")
1007 1 : {
1008 1 : TraverseAABBBranchless<dim> traverse(aabb, r, productOfCoefficientsPerDimension);
1009 1 : traverse.isInBoundingBox();
1010 :
1011 1 : size_t iter = 0;
1012 65 : while (traverse.isInBoundingBox()) {
1013 64 : RealVector_t voxel = traverse.getCurrentVoxel().template cast<real_t>();
1014 64 : INFO("Current Voxel: (" << voxel(0) << ", " << voxel(1) << ") in iter: " << iter);
1015 :
1016 64 : REQUIRE_UNARY(intersect(voxel, r));
1017 64 : traverse.updateTraverse();
1018 64 : iter++;
1019 64 : }
1020 1 : }
1021 1 : }
1022 8 : }
|
