Line data Source code
1 : #include "doctest/doctest.h"
2 :
3 : #include "SiddonsMethodCUDA.h"
4 : #include "SiddonsMethod.h"
5 : #include "Geometry.h"
6 : #include "testHelpers.h"
7 : #include "Logger.h"
8 : #include "VolumeDescriptor.h"
9 : #include "PlanarDetectorDescriptor.h"
10 : #include "testHelpers.h"
11 :
12 : #include <array>
13 :
14 : using namespace elsa;
15 : using namespace elsa::geometry;
16 : using namespace doctest;
17 :
18 : // TODO(dfrank): remove this and replace with checkApproxEq
19 : using doctest::Approx;
20 :
21 : /*
22 : * this function declaration can be used in conjunction with decltype to deduce the
23 : * template parameter of a templated class at compile time
24 : *
25 : * the template parameter must be a typename
26 : */
27 : template <template <typename> typename T, typename data_t>
28 : constexpr data_t return_data_t(const T<data_t>&);
29 :
30 116 : TEST_CASE_TEMPLATE("SiddonsMethodCUDA: Calls to functions of super class", TestType,
31 : SiddonsMethodCUDA<float>, SiddonsMethodCUDA<double>, SiddonsMethod<float>,
32 : SiddonsMethod<double>)
33 : {
34 : // Turn logger of
35 16 : Logger::setLevel(Logger::LogLevel::OFF);
36 :
37 : using data_t = decltype(return_data_t(std::declval<TestType>()));
38 32 : GIVEN("A projector")
39 : {
40 32 : IndexVector_t volumeDims(2), sinoDims(2);
41 16 : const index_t volSize = 10;
42 16 : const index_t detectorSize = 10;
43 16 : const index_t numImgs = 10;
44 16 : volumeDims << volSize, volSize;
45 16 : sinoDims << detectorSize, numImgs;
46 32 : VolumeDescriptor volumeDescriptor(volumeDims);
47 32 : DataContainer<data_t> volume(volumeDescriptor);
48 16 : volume = 1;
49 :
50 16 : auto stc = SourceToCenterOfRotation{20 * volSize};
51 16 : auto ctr = CenterOfRotationToDetector{volSize};
52 :
53 32 : std::vector<Geometry> geom;
54 176 : for (index_t i = 0; i < numImgs; i++) {
55 160 : real_t angle = static_cast<real_t>(i) * 2 * pi_t / 10;
56 160 : geom.emplace_back(stc, ctr, Radian{angle}, VolumeData2D{Size2D{volumeDims}},
57 : SinogramData2D{Size2D{sinoDims}});
58 : }
59 :
60 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
61 32 : DataContainer<data_t> sino(sinoDescriptor);
62 16 : sino = 0;
63 :
64 32 : TestType op(volumeDescriptor, sinoDescriptor);
65 :
66 32 : WHEN("Projector is cloned")
67 : {
68 32 : auto opClone = op.clone();
69 32 : auto sinoClone = DataContainer<data_t>(sinoDescriptor);
70 32 : auto volumeClone = DataContainer<data_t>(volumeDescriptor);
71 :
72 32 : THEN("Results do not change (may still be slightly different due to summation being "
73 : "performed in a different order)")
74 : {
75 16 : op.apply(volume, sino);
76 16 : opClone->apply(volume, sinoClone);
77 16 : REQUIRE_UNARY(isApprox<data_t>(sino, sinoClone, epsilon));
78 :
79 16 : op.applyAdjoint(sino, volume);
80 16 : opClone->applyAdjoint(sino, volumeClone);
81 16 : REQUIRE_UNARY(isApprox<data_t>(volume, volumeClone, epsilon));
82 : }
83 : }
84 : }
85 16 : }
86 :
87 388 : TEST_CASE_TEMPLATE("SiddonsMethodCUDA: 2D setup with a single ray", TestType,
88 : SiddonsMethodCUDA<float>, SiddonsMethodCUDA<double>, SiddonsMethod<float>,
89 : SiddonsMethod<double>)
90 : {
91 : // Turn logger of
92 288 : Logger::setLevel(Logger::LogLevel::OFF);
93 :
94 : using data_t = decltype(return_data_t(std::declval<TestType>()));
95 :
96 448 : GIVEN("A 5x5 volume and detector of size 1, with 1 angle")
97 : {
98 : // Domain setup
99 160 : const index_t volSize = 5;
100 :
101 320 : IndexVector_t volumeDims(2);
102 160 : volumeDims << volSize, volSize;
103 :
104 320 : VolumeDescriptor volumeDescriptor(volumeDims);
105 320 : DataContainer<data_t> volume(volumeDescriptor);
106 :
107 : // range setup
108 160 : const index_t detectorSize = 1;
109 160 : const index_t numImgs = 1;
110 :
111 320 : IndexVector_t sinoDims(2);
112 160 : sinoDims << detectorSize, numImgs;
113 :
114 : // Setup geometry
115 160 : auto stc = SourceToCenterOfRotation{20 * volSize};
116 160 : auto ctr = CenterOfRotationToDetector{volSize};
117 480 : auto volData = VolumeData2D{Size2D{volumeDims}};
118 480 : auto sinoData = SinogramData2D{Size2D{sinoDims}};
119 :
120 320 : std::vector<Geometry> geom;
121 :
122 192 : GIVEN("A basic geometry setup")
123 : {
124 32 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData));
125 :
126 64 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
127 64 : DataContainer<data_t> sino(sinoDescriptor);
128 :
129 64 : TestType op(volumeDescriptor, sinoDescriptor);
130 :
131 48 : WHEN("Sinogram conatainer is not zero initialized and we project through an empty "
132 : "volume")
133 : {
134 16 : volume = 0;
135 16 : sino = 1;
136 :
137 32 : THEN("Result is zero")
138 : {
139 16 : op.apply(volume, sino);
140 16 : DataContainer<data_t> zero(sinoDescriptor);
141 16 : zero = 0;
142 16 : REQUIRE_UNARY(isApprox(sino, zero, epsilon));
143 : }
144 : }
145 :
146 48 : WHEN("Volume container is not zero initialized and we backproject from an empty "
147 : "sinogram")
148 : {
149 16 : sino = 0;
150 16 : volume = 1;
151 :
152 32 : THEN("Result is zero")
153 : {
154 16 : op.applyAdjoint(sino, volume);
155 16 : DataContainer<data_t> zero(volumeDescriptor);
156 16 : zero = 0;
157 16 : REQUIRE_UNARY(isApprox(volume, zero, epsilon));
158 : }
159 : }
160 : }
161 :
162 224 : GIVEN("Scenario: Rays not intersecting the bounding box are present")
163 : {
164 80 : WHEN("Tracing along a y-axis-aligned ray with a negative x-coordinate of origin")
165 : {
166 16 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
167 : PrincipalPointOffset{}, RotationOffset2D{volSize, 0});
168 :
169 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
170 32 : DataContainer<data_t> sino(sinoDescriptor);
171 :
172 32 : TestType op(volumeDescriptor, sinoDescriptor);
173 :
174 32 : THEN("Result of forward projection is zero")
175 : {
176 32 : DataContainer<data_t> zero(sinoDescriptor);
177 16 : zero = 0;
178 :
179 16 : op.apply(volume, sino);
180 16 : REQUIRE_UNARY(isApprox(sino, zero));
181 :
182 32 : AND_THEN("Result of backprojection is zero")
183 : {
184 16 : DataContainer<data_t> zero(volumeDescriptor);
185 16 : zero = 0;
186 :
187 16 : op.applyAdjoint(sino, volume);
188 16 : REQUIRE_UNARY(isApprox(volume, zero));
189 : }
190 : }
191 : }
192 :
193 80 : WHEN("Tracing along a y-axis-aligned ray with a x-coordinate of origin beyond the "
194 : "bounding box")
195 : {
196 16 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
197 : PrincipalPointOffset{}, RotationOffset2D{-volSize, 0});
198 :
199 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
200 32 : DataContainer<data_t> sino(sinoDescriptor);
201 :
202 32 : TestType op(volumeDescriptor, sinoDescriptor);
203 :
204 32 : THEN("Result of forward projection is zero")
205 : {
206 32 : DataContainer<data_t> zero(sinoDescriptor);
207 16 : zero = 0;
208 :
209 16 : op.apply(volume, sino);
210 16 : REQUIRE_UNARY(isApprox(sino, zero));
211 :
212 32 : AND_THEN("Result of backprojection is zero")
213 : {
214 16 : DataContainer<data_t> zero(volumeDescriptor);
215 16 : zero = 0;
216 :
217 16 : op.applyAdjoint(sino, volume);
218 16 : REQUIRE_UNARY(isApprox(volume, zero));
219 : }
220 : }
221 : }
222 :
223 80 : WHEN("Tracing along a x-axis-aligned ray with a negative y-coordinate of origin")
224 : {
225 32 : geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData),
226 16 : std::move(sinoData), PrincipalPointOffset{},
227 : RotationOffset2D{0, volSize});
228 :
229 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
230 32 : DataContainer<data_t> sino(sinoDescriptor);
231 :
232 32 : TestType op(volumeDescriptor, sinoDescriptor);
233 :
234 32 : THEN("Result of forward projection is zero")
235 : {
236 32 : DataContainer<data_t> zero(sinoDescriptor);
237 16 : zero = 0;
238 :
239 16 : op.apply(volume, sino);
240 16 : REQUIRE_UNARY(isApprox(sino, zero));
241 :
242 32 : AND_THEN("Result of backprojection is zero")
243 : {
244 16 : DataContainer<data_t> zero(volumeDescriptor);
245 16 : zero = 0;
246 :
247 16 : op.applyAdjoint(sino, volume);
248 16 : REQUIRE_UNARY(isApprox(volume, zero));
249 : }
250 : }
251 : }
252 :
253 80 : WHEN("Tracing along a x-axis-aligned ray with a y-coordinate of origin beyond the "
254 : "bounding box")
255 : {
256 32 : geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData),
257 16 : std::move(sinoData), PrincipalPointOffset{},
258 : RotationOffset2D{0, -volSize});
259 :
260 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
261 32 : DataContainer<data_t> sino(sinoDescriptor);
262 :
263 32 : TestType op(volumeDescriptor, sinoDescriptor);
264 :
265 32 : THEN("Result of forward projection is zero")
266 : {
267 16 : op.apply(volume, sino);
268 32 : DataContainer<data_t> zero(sinoDescriptor);
269 16 : zero = 0;
270 :
271 16 : REQUIRE_UNARY(isApprox(sino, zero));
272 :
273 32 : AND_THEN("Result of backprojection is zero")
274 : {
275 16 : op.applyAdjoint(sino, volume);
276 :
277 16 : DataContainer<data_t> zero(volumeDescriptor);
278 16 : zero = 0;
279 :
280 16 : REQUIRE_UNARY(isApprox(volume, zero));
281 : }
282 : }
283 : }
284 : }
285 :
286 : // Expected results
287 960 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProjections[2];
288 160 : backProjections[0].resize(volSize * volSize);
289 160 : backProjections[1].resize(volSize * volSize);
290 :
291 160 : constexpr index_t numCases = 4;
292 160 : const std::array<real_t, numCases> angles = {0., 90., 180., 270.};
293 :
294 : // clang-format off
295 320 : backProjections[0] << 0, 0, 1, 0, 0,
296 160 : 0, 0, 1, 0, 0,
297 160 : 0, 0, 1, 0, 0,
298 160 : 0, 0, 1, 0, 0,
299 320 : 0, 0, 1, 0, 0;
300 :
301 320 : backProjections[1] << 0, 0, 0, 0, 0,
302 160 : 0, 0, 0, 0, 0,
303 160 : 1, 1, 1, 1, 1,
304 160 : 0, 0, 0, 0, 0,
305 320 : 0, 0, 0, 0, 0;
306 : // clang-format on
307 :
308 800 : for (index_t i = 0; i < numCases; i++) {
309 656 : WHEN("An axis-aligned ray with a fixed angle passes through the center of a pixel")
310 : {
311 32 : INFO("An axis-aligned ray with an angle of ", angles[asUnsigned(i)],
312 : " radians passes through the center of a pixel");
313 :
314 16 : geom.emplace_back(stc, ctr, Degree{angles[asUnsigned(i)]}, std::move(volData),
315 16 : std::move(sinoData));
316 :
317 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
318 32 : DataContainer<data_t> sino(sinoDescriptor);
319 :
320 32 : TestType op(volumeDescriptor, sinoDescriptor);
321 :
322 32 : THEN("The result of projecting through a pixel is exactly the pixel value")
323 : {
324 96 : for (index_t j = 0; j < volSize; j++) {
325 80 : volume = 0;
326 80 : if (i % 2 == 0)
327 80 : volume(volSize / 2, j) = 1;
328 : else
329 0 : volume(j, volSize / 2) = 1;
330 :
331 80 : op.apply(volume, sino);
332 80 : REQUIRE_EQ(sino[0], Approx(1));
333 : }
334 :
335 32 : AND_THEN("The backprojection sets the values of all hit pixels to the detector "
336 : "value")
337 : {
338 16 : op.applyAdjoint(sino, volume);
339 16 : REQUIRE_UNARY(
340 : isApprox(volume, DataContainer<data_t>(volumeDescriptor,
341 : backProjections[i % 2])));
342 : }
343 : }
344 : }
345 : }
346 :
347 176 : WHEN("A y-axis-aligned ray runs along the left voxel boundary")
348 : {
349 32 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, Radian{0},
350 16 : std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
351 : RotationOffset2D{-0.5, 0});
352 :
353 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
354 32 : DataContainer<data_t> sino(sinoDescriptor);
355 :
356 32 : TestType op(volumeDescriptor, sinoDescriptor);
357 :
358 32 : THEN("The result of projecting through a pixel is the value of the pixel with the "
359 : "higher index")
360 : {
361 96 : for (index_t j = 0; j < volSize; j++) {
362 80 : volume = 0;
363 80 : volume(volSize / 2, j) = 1;
364 :
365 80 : op.apply(volume, sino);
366 80 : REQUIRE_EQ(sino[0], Approx(1.0));
367 : }
368 :
369 32 : AND_THEN("The backprojection yields the exact adjoint")
370 : {
371 16 : sino[0] = 1;
372 16 : op.applyAdjoint(sino, volume);
373 16 : REQUIRE_UNARY(isApprox(
374 : volume, DataContainer<data_t>(volumeDescriptor, backProjections[0])));
375 : }
376 : }
377 : }
378 :
379 176 : WHEN("A y-axis-aligned ray runs along the right volume boundary")
380 : {
381 : // For Siddons's values in the range [0,boxMax) are considered, a ray running along
382 : // boxMax should be ignored
383 32 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, Radian{0},
384 16 : std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
385 : RotationOffset2D{volSize * 0.5, 0});
386 :
387 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
388 32 : DataContainer<data_t> sino(sinoDescriptor);
389 :
390 32 : TestType op(volumeDescriptor, sinoDescriptor);
391 :
392 32 : THEN("The result of projecting is zero")
393 : {
394 16 : volume = 0;
395 16 : op.apply(volume, sino);
396 16 : REQUIRE_EQ(sino[0], Approx(0));
397 :
398 32 : AND_THEN("The result of backprojection is also zero")
399 : {
400 16 : sino[0] = 1;
401 :
402 16 : op.applyAdjoint(sino, volume);
403 16 : DataContainer<data_t> zero(volumeDescriptor);
404 16 : zero = 0;
405 16 : REQUIRE_UNARY(isApprox(volume, zero));
406 : }
407 : }
408 : }
409 :
410 : // clang-format off
411 320 : backProjections[0] << 1, 0, 0, 0, 0,
412 160 : 1, 0, 0, 0, 0,
413 160 : 1, 0, 0, 0, 0,
414 160 : 1, 0, 0, 0, 0,
415 320 : 1, 0, 0, 0, 0;
416 : // clang-format on
417 :
418 176 : WHEN("A y-axis-aligned ray runs along the left volume boundary")
419 : {
420 32 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, Radian{0},
421 16 : std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
422 : RotationOffset2D{-volSize * 0.5, 0});
423 :
424 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
425 32 : DataContainer<data_t> sino(sinoDescriptor);
426 :
427 32 : TestType op(volumeDescriptor, sinoDescriptor);
428 :
429 32 : THEN("The result of projecting through a pixel is exactly the pixel's value")
430 : {
431 96 : for (index_t j = 0; j < volSize; j++) {
432 80 : volume = 0;
433 80 : volume(0, j) = 1;
434 :
435 80 : op.apply(volume, sino);
436 80 : REQUIRE_EQ(sino[0], Approx(1));
437 : }
438 :
439 32 : AND_THEN("The backprojection yields the exact adjoint")
440 : {
441 16 : sino[0] = 1;
442 16 : op.applyAdjoint(sino, volume);
443 16 : REQUIRE_UNARY(isApprox(
444 : volume, DataContainer<data_t>(volumeDescriptor, backProjections[0])));
445 : }
446 : }
447 : }
448 : }
449 416 : GIVEN("a 4x4 Volume")
450 : {
451 : // Domain setup
452 128 : const index_t volSize = 4;
453 :
454 256 : IndexVector_t volumeDims(2);
455 128 : volumeDims << volSize, volSize;
456 :
457 256 : VolumeDescriptor volumeDescriptor(volumeDims);
458 256 : DataContainer<data_t> volume(volumeDescriptor);
459 :
460 : // range setup
461 128 : const index_t detectorSize = 1;
462 128 : const index_t numImgs = 1;
463 :
464 256 : IndexVector_t sinoDims(2);
465 128 : sinoDims << detectorSize, numImgs;
466 :
467 : // Setup geometry
468 128 : auto stc = SourceToCenterOfRotation{20 * volSize};
469 128 : auto ctr = CenterOfRotationToDetector{volSize};
470 384 : auto volData = VolumeData2D{Size2D{volumeDims}};
471 384 : auto sinoData = SinogramData2D{Size2D{sinoDims}};
472 :
473 256 : std::vector<Geometry> geom;
474 :
475 128 : real_t sqrt3r = std::sqrt(static_cast<real_t>(3));
476 128 : data_t sqrt3d = std::sqrt(static_cast<data_t>(3));
477 :
478 192 : GIVEN("An angle of -30 degrees")
479 : {
480 64 : auto angle = Degree{-30};
481 80 : WHEN("A ray goes through center of volume")
482 : {
483 : // In this case the ray enters and exits the volume through the borders along the
484 : // main direction
485 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData));
486 :
487 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
488 32 : DataContainer<data_t> sino(sinoDescriptor);
489 :
490 32 : TestType op(volumeDescriptor, sinoDescriptor);
491 :
492 32 : THEN("Ray intersects the correct pixels")
493 : {
494 16 : volume = 1;
495 16 : volume(3, 0) = 0;
496 16 : volume(2, 0) = 0;
497 16 : volume(2, 1) = 0;
498 :
499 16 : volume(1, 2) = 0;
500 16 : volume(1, 3) = 0;
501 16 : volume(0, 3) = 0;
502 16 : volume(2, 2) = 0;
503 :
504 32 : DataContainer<data_t> zero(sinoDescriptor);
505 16 : zero = 0;
506 :
507 16 : op.apply(volume, sino);
508 16 : REQUIRE_EQ(zero[0], Approx(0).epsilon(epsilon));
509 :
510 32 : AND_THEN("The correct weighting is applied")
511 : {
512 16 : volume(3, 0) = 1;
513 16 : volume(2, 0) = 2;
514 16 : volume(2, 1) = 3;
515 :
516 16 : op.apply(volume, sino);
517 16 : REQUIRE_EQ(sino[0], Approx(2 * sqrt3d + 2));
518 :
519 : // on the other side of the center
520 16 : volume = 0;
521 16 : volume(1, 2) = 3;
522 16 : volume(1, 3) = 2;
523 16 : volume(0, 3) = 1;
524 :
525 16 : op.apply(volume, sino);
526 16 : REQUIRE_EQ(sino[0], Approx(2 * sqrt3d + 2));
527 :
528 16 : sino[0] = 1;
529 :
530 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
531 :
532 : // clang-format off
533 32 : expected << 0, 0, 2 - 2 / sqrt3d, 4 / sqrt3d - 2,
534 16 : 0, 0, 2 / sqrt3d, 0,
535 16 : 0, 2 / sqrt3d, 0, 0,
536 32 : 4 / sqrt3d - 2, 2 - 2 / sqrt3d, 0, 0;
537 : // clang-format on
538 :
539 16 : op.applyAdjoint(sino, volume);
540 16 : REQUIRE_UNARY(isApprox(
541 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
542 : }
543 : }
544 : }
545 :
546 80 : WHEN("A ray enters through the right border")
547 : {
548 : // In this case the ray exits through a border along the main ray direction, but
549 : // enters through a border not along the main direction First pixel should be
550 : // weighted differently
551 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
552 : PrincipalPointOffset{0}, RotationOffset2D{sqrt3r, 0});
553 :
554 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
555 32 : DataContainer<data_t> sino(sinoDescriptor);
556 :
557 32 : TestType op(volumeDescriptor, sinoDescriptor);
558 :
559 32 : THEN("Ray intersects the correct pixels")
560 : {
561 16 : volume = 1;
562 16 : volume(3, 1) = 0;
563 16 : volume(3, 2) = 0;
564 16 : volume(3, 3) = 0;
565 16 : volume(2, 3) = 0;
566 :
567 16 : op.apply(volume, sino);
568 32 : DataContainer<data_t> zero(sinoDescriptor);
569 16 : zero = 0;
570 16 : REQUIRE_UNARY(isApprox(sino, zero, epsilon));
571 :
572 32 : AND_THEN("The correct weighting is applied")
573 : {
574 16 : volume(3, 1) = 4;
575 16 : volume(3, 2) = 3;
576 16 : volume(3, 3) = 2;
577 16 : volume(2, 3) = 1;
578 :
579 16 : op.apply(volume, sino);
580 16 : REQUIRE_EQ(sino[0], Approx(14 - 4 * sqrt3d));
581 :
582 16 : sino[0] = 1;
583 :
584 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
585 :
586 : // clang-format off
587 32 : expected << 0, 0, 0, 0,
588 16 : 0, 0, 0, 4 - 2 * sqrt3d,
589 16 : 0, 0, 0, 2 / sqrt3d,
590 32 : 0, 0, 2 - 2 / sqrt3d, 4 / sqrt3d - 2;
591 : // clang-format on
592 :
593 16 : op.applyAdjoint(sino, volume);
594 16 : REQUIRE_UNARY(isApprox(
595 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
596 : }
597 : }
598 : }
599 :
600 80 : WHEN("A ray exits through the left border")
601 : {
602 : // In this case the ray enters through a border along the main ray direction, but
603 : // exits through a border not along the main direction
604 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
605 : PrincipalPointOffset{0}, RotationOffset2D{-sqrt3r, 0});
606 :
607 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
608 32 : DataContainer<data_t> sino(sinoDescriptor);
609 :
610 32 : TestType op(volumeDescriptor, sinoDescriptor);
611 :
612 32 : THEN("Ray intersects the correct pixels")
613 : {
614 16 : volume = 1;
615 16 : volume(0, 0) = 0;
616 16 : volume(1, 0) = 0;
617 16 : volume(0, 1) = 0;
618 16 : volume(0, 2) = 0;
619 :
620 16 : op.apply(volume, sino);
621 32 : DataContainer<data_t> zero(sinoDescriptor);
622 16 : zero = 0;
623 16 : REQUIRE_UNARY(isApprox(sino, zero, epsilon));
624 :
625 32 : AND_THEN("The correct weighting is applied")
626 : {
627 16 : volume(1, 0) = 1;
628 16 : volume(0, 0) = 2;
629 16 : volume(0, 1) = 3;
630 16 : volume(0, 2) = 4;
631 :
632 16 : op.apply(volume, sino);
633 16 : REQUIRE_EQ(sino[0], Approx(14 - 4 * sqrt3d));
634 :
635 16 : sino[0] = 1;
636 :
637 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
638 :
639 : // clang-format off
640 32 : expected << 4 / sqrt3d - 2, 2 - 2 / sqrt3d, 0, 0,
641 16 : 2 / sqrt3d, 0, 0, 0,
642 16 : 4 - 2 * sqrt3d, 0, 0, 0,
643 32 : 0, 0, 0, 0;
644 : // clang-format on
645 :
646 16 : op.applyAdjoint(sino, volume);
647 16 : REQUIRE_UNARY(isApprox(
648 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
649 : }
650 : }
651 : }
652 :
653 80 : WHEN("A ray only intersects a single pixel")
654 : {
655 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
656 16 : PrincipalPointOffset{0}, RotationOffset2D{-2 - sqrt3r / 2, 0});
657 :
658 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
659 32 : DataContainer<data_t> sino(sinoDescriptor);
660 :
661 32 : TestType op(volumeDescriptor, sinoDescriptor);
662 :
663 32 : THEN("Ray intersects the correct pixels")
664 : {
665 16 : volume = 1;
666 16 : volume(0, 0) = 0;
667 :
668 16 : op.apply(volume, sino);
669 32 : DataContainer<data_t> zero(sinoDescriptor);
670 16 : zero = 0;
671 16 : REQUIRE_UNARY(isApprox(sino, zero, epsilon));
672 :
673 32 : AND_THEN("The correct weighting is applied")
674 : {
675 16 : volume(0, 0) = 1;
676 :
677 16 : op.apply(volume, sino);
678 16 : REQUIRE_EQ(sino[0], Approx(1 / sqrt3d));
679 :
680 16 : sino[0] = 1;
681 :
682 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
683 : // clang-format off
684 32 : expected << 1 / sqrt3d, 0, 0, 0,
685 16 : 0, 0, 0, 0,
686 16 : 0, 0, 0, 0,
687 32 : 0, 0, 0, 0;
688 : // clang-format on
689 :
690 16 : op.applyAdjoint(sino, volume);
691 16 : REQUIRE_UNARY(isApprox(
692 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
693 : }
694 : }
695 : }
696 : }
697 :
698 192 : GIVEN("An angle of 120 degrees")
699 : {
700 64 : auto angle = Degree{-120};
701 :
702 80 : WHEN("A ray goes through center of volume")
703 : {
704 : // In this case the ray enters and exits the volume through the borders along the
705 : // main direction
706 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData));
707 :
708 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
709 32 : DataContainer<data_t> sino(sinoDescriptor);
710 :
711 32 : TestType op(volumeDescriptor, sinoDescriptor);
712 :
713 32 : THEN("Ray intersects the correct pixels")
714 : {
715 16 : volume = 1;
716 16 : volume(0, 0) = 0;
717 16 : volume(0, 1) = 0;
718 16 : volume(1, 1) = 0;
719 16 : volume(2, 2) = 0;
720 16 : volume(3, 2) = 0;
721 16 : volume(3, 3) = 0;
722 : // volume(1,2) hit due to numerical error
723 16 : volume(1, 2) = 0;
724 :
725 16 : op.apply(volume, sino);
726 32 : DataContainer<data_t> zero(sinoDescriptor);
727 16 : zero = 0;
728 16 : REQUIRE_EQ(zero[0], Approx(0).epsilon(epsilon));
729 :
730 32 : AND_THEN("The correct weighting is applied")
731 : {
732 16 : volume(0, 0) = 1;
733 16 : volume(0, 1) = 2;
734 16 : volume(1, 1) = 3;
735 :
736 16 : op.apply(volume, sino);
737 16 : REQUIRE_EQ(sino[0], Approx(2 * sqrt3d + 2));
738 :
739 : // on the other side of the center
740 16 : volume = 0;
741 16 : volume(2, 2) = 3;
742 16 : volume(3, 2) = 2;
743 16 : volume(3, 3) = 1;
744 :
745 16 : op.apply(volume, sino);
746 16 : REQUIRE_EQ(sino[0], Approx(2 * sqrt3d + 2));
747 :
748 16 : sino[0] = 1;
749 :
750 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
751 :
752 : // clang-format off
753 32 : expected << 4 / sqrt3d - 2, 0, 0, 0,
754 16 : 2 - 2 / sqrt3d, 2 / sqrt3d, 0, 0,
755 16 : 0, 0, 2 / sqrt3d, 2 - 2 / sqrt3d,
756 32 : 0, 0, 0, 4 / sqrt3d - 2;
757 : // clang-format on
758 :
759 16 : op.applyAdjoint(sino, volume);
760 16 : REQUIRE_UNARY(isApprox(
761 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
762 : }
763 : }
764 : }
765 :
766 80 : WHEN("A ray enters through the top border")
767 : {
768 : // In this case the ray exits through a border along the main ray direction, but
769 : // enters through a border not along the main direction
770 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
771 16 : PrincipalPointOffset{0}, RotationOffset2D{0, std::sqrt(3.f)});
772 :
773 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
774 32 : DataContainer<data_t> sino(sinoDescriptor);
775 :
776 32 : TestType op(volumeDescriptor, sinoDescriptor);
777 :
778 32 : THEN("Ray intersects the correct pixels")
779 : {
780 16 : volume = 1;
781 16 : volume(0, 2) = 0;
782 16 : volume(0, 3) = 0;
783 16 : volume(1, 3) = 0;
784 16 : volume(2, 3) = 0;
785 :
786 16 : op.apply(volume, sino);
787 32 : DataContainer<data_t> zero(sinoDescriptor);
788 16 : zero = 0;
789 16 : REQUIRE_EQ(zero[0], Approx(0).epsilon(epsilon));
790 :
791 32 : AND_THEN("The correct weighting is applied")
792 : {
793 16 : volume(0, 2) = 1;
794 16 : volume(0, 3) = 2;
795 16 : volume(1, 3) = 3;
796 16 : volume(2, 3) = 4;
797 :
798 16 : op.apply(volume, sino);
799 16 : REQUIRE_EQ(sino[0], Approx(14 - 4 * sqrt3d));
800 :
801 16 : sino[0] = 1;
802 :
803 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
804 :
805 : // clang-format off
806 32 : expected << 0, 0, 0, 0,
807 16 : 0, 0, 0, 0,
808 16 : 2 - 2 / sqrt3d, 0, 0, 0,
809 32 : 4 / sqrt3d - 2, 2 / sqrt3d, 4 - 2 * sqrt3d, 0;
810 : // clang-format on
811 :
812 16 : op.applyAdjoint(sino, volume);
813 16 : REQUIRE_UNARY(isApprox(
814 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
815 : }
816 : }
817 : }
818 :
819 80 : WHEN("A ray exits through the bottom border")
820 :
821 : {
822 : // In this case the ray enters through a border along the main ray direction, but
823 : // exits through a border not along the main direction
824 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
825 16 : PrincipalPointOffset{0}, RotationOffset2D{0, -std::sqrt(3.f)});
826 :
827 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
828 32 : DataContainer<data_t> sino(sinoDescriptor);
829 :
830 32 : TestType op(volumeDescriptor, sinoDescriptor);
831 :
832 32 : THEN("Ray intersects the correct pixels")
833 : {
834 16 : volume = 1;
835 16 : volume(1, 0) = 0;
836 16 : volume(2, 0) = 0;
837 16 : volume(3, 0) = 0;
838 16 : volume(3, 1) = 0;
839 :
840 16 : op.apply(volume, sino);
841 32 : DataContainer<data_t> zero(sinoDescriptor);
842 16 : zero = 0;
843 16 : REQUIRE_EQ(zero[0], Approx(0).epsilon(epsilon));
844 :
845 32 : AND_THEN("The correct weighting is applied")
846 : {
847 16 : volume(1, 0) = 4;
848 16 : volume(2, 0) = 3;
849 16 : volume(3, 0) = 2;
850 16 : volume(3, 1) = 1;
851 :
852 16 : op.apply(volume, sino);
853 16 : REQUIRE_EQ(sino[0], Approx(14 - 4 * sqrt3d));
854 :
855 16 : sino[0] = 1;
856 :
857 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
858 :
859 : // clang-format off
860 32 : expected << 0, 4 - 2 * sqrt3d, 2 / sqrt3d, 4 / sqrt3d - 2, 0, 0, 0,
861 32 : 2 - 2 / sqrt3d, 0, 0, 0, 0, 0, 0, 0, 0;
862 : // clang-format on
863 :
864 16 : op.applyAdjoint(sino, volume);
865 16 : REQUIRE_UNARY(isApprox(
866 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
867 : }
868 : }
869 : }
870 :
871 80 : WHEN("A ray only intersects a single pixel")
872 : {
873 : // This is a special case that is handled separately in both forward and
874 : // backprojection
875 16 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
876 : PrincipalPointOffset{0},
877 16 : RotationOffset2D{0, -2 - std::sqrt(3.f) / 2});
878 :
879 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
880 32 : DataContainer<data_t> sino(sinoDescriptor);
881 :
882 32 : TestType op(volumeDescriptor, sinoDescriptor);
883 :
884 32 : THEN("Ray intersects the correct pixels")
885 : {
886 16 : volume = 1;
887 16 : volume(3, 0) = 0;
888 :
889 16 : op.apply(volume, sino);
890 32 : DataContainer<data_t> zero(sinoDescriptor);
891 16 : zero = 0;
892 16 : REQUIRE_EQ(zero[0], Approx(0).epsilon(epsilon));
893 :
894 32 : AND_THEN("The correct weighting is applied")
895 : {
896 16 : volume(3, 0) = 1;
897 :
898 16 : op.apply(volume, sino);
899 16 : REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
900 :
901 16 : sino[0] = 1;
902 :
903 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> expected(volSize * volSize);
904 :
905 : // clang-format off
906 16 : expected << 0, 0, 0, 1 / sqrt3d, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
907 : // clang-format on
908 :
909 16 : op.applyAdjoint(sino, volume);
910 16 : REQUIRE_UNARY(isApprox(
911 : volume, DataContainer<data_t>(volumeDescriptor, expected), epsilon));
912 : }
913 : }
914 : }
915 : }
916 : }
917 288 : }
918 :
919 116 : TEST_CASE_TEMPLATE("SiddonsMethodCUDA: 2D setup with a multiple rays", TestType,
920 : SiddonsMethodCUDA<float>, SiddonsMethodCUDA<double>, SiddonsMethod<float>,
921 : SiddonsMethod<double>)
922 : {
923 : // Turn logger of
924 16 : Logger::setLevel(Logger::LogLevel::OFF);
925 :
926 : using data_t = decltype(return_data_t(std::declval<TestType>()));
927 :
928 32 : GIVEN("Given a 5x5 volume")
929 : {
930 : // Domain setup
931 16 : const index_t volSize = 5;
932 :
933 32 : IndexVector_t volumeDims(2);
934 16 : volumeDims << volSize, volSize;
935 :
936 32 : VolumeDescriptor volumeDescriptor(volumeDims);
937 32 : DataContainer<data_t> volume(volumeDescriptor);
938 :
939 : // range setup
940 16 : const index_t detectorSize = 1;
941 16 : const index_t numImgs = 4;
942 :
943 32 : IndexVector_t sinoDims(2);
944 16 : sinoDims << detectorSize, numImgs;
945 :
946 : // Setup geometry
947 16 : auto stc = SourceToCenterOfRotation{20 * volSize};
948 16 : auto ctr = CenterOfRotationToDetector{volSize};
949 48 : auto volData = VolumeData2D{Size2D{volumeDims}};
950 48 : auto sinoData = SinogramData2D{Size2D{sinoDims}};
951 :
952 32 : std::vector<Geometry> geom;
953 :
954 32 : WHEN("Both x- and y-axis-aligned rays are present")
955 : {
956 16 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{volumeDims}},
957 : SinogramData2D{Size2D{sinoDims}});
958 16 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{volumeDims}},
959 : SinogramData2D{Size2D{sinoDims}});
960 16 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{volumeDims}},
961 : SinogramData2D{Size2D{sinoDims}});
962 16 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{volumeDims}},
963 : SinogramData2D{Size2D{sinoDims}});
964 :
965 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
966 32 : DataContainer<data_t> sino(sinoDescriptor);
967 :
968 32 : TestType op(volumeDescriptor, sinoDescriptor);
969 :
970 32 : THEN("Values are accumulated correctly along each ray's path")
971 : {
972 16 : volume = 0;
973 :
974 : // set only values along the rays' path to one to make sure interpolation is
975 : // dones correctly
976 96 : for (index_t i = 0; i < volSize; i++) {
977 80 : volume(i, volSize / 2) = 1;
978 80 : volume(volSize / 2, i) = 1;
979 : }
980 :
981 16 : op.apply(volume, sino);
982 80 : for (index_t i = 0; i < numImgs; i++)
983 64 : REQUIRE_EQ(sino[i], Approx(5.0));
984 :
985 32 : AND_THEN("Backprojection yields the exact adjoint")
986 : {
987 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> cmp(volSize * volSize);
988 :
989 : // clang-format off
990 32 : cmp << 0, 0, 10, 0, 0,
991 16 : 0, 0, 10, 0, 0,
992 16 : 10, 10, 20, 10, 10,
993 16 : 0, 0, 10, 0, 0,
994 32 : 0, 0, 10, 0, 0;
995 : // clang-format on
996 :
997 16 : op.applyAdjoint(sino, volume);
998 16 : REQUIRE_UNARY(isApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
999 : }
1000 : }
1001 : }
1002 : }
1003 16 : }
1004 :
1005 292 : TEST_CASE_TEMPLATE("SiddonsMethodCUDA: 3D setup with a single ray", TestType,
1006 : SiddonsMethodCUDA<float>, SiddonsMethodCUDA<double>, SiddonsMethod<float>,
1007 : SiddonsMethod<double>)
1008 : {
1009 : // Turn logger of
1010 192 : Logger::setLevel(Logger::LogLevel::OFF);
1011 :
1012 : using data_t = decltype(return_data_t(std::declval<TestType>()));
1013 :
1014 352 : GIVEN("Given a 3x3x3 volume")
1015 : {
1016 : // Domain setup
1017 160 : const index_t volSize = 3;
1018 :
1019 320 : IndexVector_t volumeDims(3);
1020 160 : volumeDims << volSize, volSize, volSize;
1021 :
1022 320 : VolumeDescriptor volumeDescriptor(volumeDims);
1023 320 : DataContainer<data_t> volume(volumeDescriptor);
1024 :
1025 : // range setup
1026 160 : const index_t detectorSize = 1;
1027 160 : const index_t numImgs = 1;
1028 :
1029 320 : IndexVector_t sinoDims(3);
1030 160 : sinoDims << detectorSize, detectorSize, numImgs;
1031 :
1032 : // Setup geometry
1033 160 : auto stc = SourceToCenterOfRotation{20 * volSize};
1034 160 : auto ctr = CenterOfRotationToDetector{volSize};
1035 480 : auto volData = VolumeData3D{Size3D{volumeDims}};
1036 480 : auto sinoData = SinogramData3D{Size3D{sinoDims}};
1037 :
1038 320 : std::vector<Geometry> geom;
1039 :
1040 192 : GIVEN("A basic 3D setup")
1041 : {
1042 32 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1043 : RotationAngles3D{Gamma{0}});
1044 :
1045 64 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1046 64 : DataContainer<data_t> sino(sinoDescriptor);
1047 :
1048 64 : TestType op(volumeDescriptor, sinoDescriptor);
1049 :
1050 48 : WHEN("Sinogram conatainer is not zero initialized and we project through an empty "
1051 : "volume")
1052 : {
1053 16 : volume = 0;
1054 16 : sino = 1;
1055 :
1056 32 : THEN("Result is zero")
1057 : {
1058 16 : DataContainer<data_t> zero(sinoDescriptor);
1059 16 : zero = 0;
1060 :
1061 16 : op.apply(volume, sino);
1062 16 : REQUIRE_UNARY(isApprox(sino, zero, epsilon));
1063 : }
1064 : }
1065 :
1066 48 : WHEN("Volume container is not zero initialized and we backproject from an empty "
1067 : "sinogram")
1068 : {
1069 16 : sino = 0;
1070 16 : volume = 1;
1071 :
1072 32 : THEN("Result is zero")
1073 : {
1074 16 : DataContainer<data_t> zero(volumeDescriptor);
1075 16 : zero = 0;
1076 :
1077 16 : op.applyAdjoint(sino, volume);
1078 16 : REQUIRE_UNARY(isApprox(volume, zero, epsilon));
1079 : }
1080 : }
1081 : }
1082 :
1083 224 : GIVEN("A rays along different axes")
1084 : {
1085 64 : const index_t numCases = 6;
1086 :
1087 : using RealArray = std::array<real_t, numCases>;
1088 : using StrArray = std::array<std::string, numCases>;
1089 :
1090 64 : RealArray beta = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
1091 64 : RealArray gamma = {0.0, pi_t, pi_t / 2, 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
1092 128 : StrArray al = {"z", "-z", "x", "-x", "y", "-y"};
1093 :
1094 448 : Eigen::Matrix<data_t, volSize * volSize * volSize, 1> backProj[numCases];
1095 :
1096 : // clang-format off
1097 128 : backProj[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0,
1098 64 : 0, 0, 0, 0, 1, 0, 0, 0, 0,
1099 128 : 0, 0, 0, 0, 1, 0, 0, 0, 0;
1100 :
1101 128 : backProj[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0,
1102 64 : 0, 0, 0, 1, 1, 1, 0, 0, 0,
1103 128 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1104 :
1105 128 : backProj[2] << 0, 0, 0, 0, 0, 0, 0, 0, 0,
1106 64 : 0, 1, 0, 0, 1, 0, 0, 1, 0,
1107 128 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1108 : // clang-format on
1109 :
1110 448 : for (index_t i = 0; i < numCases; i++) {
1111 400 : WHEN("An axis-aligned ray passes through the center of a pixel")
1112 : {
1113 32 : INFO("A ", al[asUnsigned(i)],
1114 : "-axis-aligned ray passes through the center of a pixel");
1115 :
1116 32 : geom.emplace_back(
1117 16 : stc, ctr, std::move(volData), std::move(sinoData),
1118 : RotationAngles3D{Gamma{gamma[asUnsigned(i)]}, Beta{beta[asUnsigned(i)]}});
1119 :
1120 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1121 32 : DataContainer<data_t> sino(sinoDescriptor);
1122 :
1123 32 : TestType op(volumeDescriptor, sinoDescriptor);
1124 :
1125 32 : THEN("The result of projecting through a voxel is exactly the voxel value")
1126 : {
1127 64 : for (index_t j = 0; j < volSize; j++) {
1128 48 : volume = 0;
1129 48 : if (i / 2 == 0)
1130 48 : volume(volSize / 2, volSize / 2, j) = 1;
1131 0 : else if (i / 2 == 1)
1132 0 : volume(j, volSize / 2, volSize / 2) = 1;
1133 0 : else if (i / 2 == 2)
1134 0 : volume(volSize / 2, j, volSize / 2) = 1;
1135 :
1136 48 : op.apply(volume, sino);
1137 48 : REQUIRE_EQ(sino[0], Approx(1));
1138 : }
1139 :
1140 32 : AND_THEN(
1141 : "The backprojection sets the values of all hit voxels to the detector "
1142 : "value")
1143 : {
1144 16 : op.applyAdjoint(sino, volume);
1145 16 : REQUIRE_UNARY(isApprox(
1146 : volume, DataContainer<data_t>(volumeDescriptor, backProj[i / 2])));
1147 : }
1148 : }
1149 : }
1150 : }
1151 :
1152 : RealArray offsetx;
1153 : RealArray offsety;
1154 :
1155 64 : offsetx[0] = volSize / 2.0;
1156 64 : offsetx[1] = 0.0;
1157 64 : offsetx[2] = (volSize / 2.0);
1158 64 : offsetx[3] = -(volSize / 2.0);
1159 64 : offsetx[4] = 0.0;
1160 64 : offsetx[5] = -(volSize / 2.0);
1161 :
1162 64 : offsety[0] = 0.0;
1163 64 : offsety[1] = volSize / 2.0;
1164 64 : offsety[2] = (volSize / 2.0);
1165 64 : offsety[3] = 0.0;
1166 64 : offsety[4] = -(volSize / 2.0);
1167 64 : offsety[5] = -(volSize / 2.0);
1168 :
1169 : // clang-format off
1170 128 : backProj[0] << 0, 0, 0, 1, 0, 0, 0, 0, 0,
1171 64 : 0, 0, 0, 1, 0, 0, 0, 0, 0,
1172 128 : 0, 0, 0, 1, 0, 0, 0, 0, 0;
1173 :
1174 128 : backProj[1] << 0, 1, 0, 0, 0, 0, 0, 0, 0,
1175 64 : 0, 1, 0, 0, 0, 0, 0, 0, 0,
1176 128 : 0, 1, 0, 0, 0, 0, 0, 0, 0;
1177 :
1178 128 : backProj[2] << 1, 0, 0, 0, 0, 0, 0, 0, 0,
1179 64 : 1, 0, 0, 0, 0, 0, 0, 0, 0,
1180 128 : 1, 0, 0, 0, 0, 0, 0, 0, 0;
1181 : // clang-format on
1182 :
1183 64 : al[0] = "left border";
1184 64 : al[1] = "bottom border";
1185 64 : al[2] = "bottom left border";
1186 64 : al[3] = "right border";
1187 64 : al[4] = "top border";
1188 64 : al[5] = "top right edge";
1189 :
1190 256 : for (unsigned int i = 0; i < numCases / 2; i++) {
1191 208 : WHEN("A z-axis-aligned ray runs along the corners and edges of the volume")
1192 : {
1193 32 : INFO("A z-axis-aligned ray runs along the ", al[i], " of the volume");
1194 :
1195 : // x-ray source must be very far from the volume center to make testing of the
1196 : // op backprojection simpler
1197 32 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr,
1198 16 : std::move(volData), std::move(sinoData),
1199 : RotationAngles3D{Gamma{0}}, PrincipalPointOffset2D{0, 0},
1200 16 : RotationOffset3D{-offsetx[i], -offsety[i], 0});
1201 :
1202 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1203 32 : DataContainer<data_t> sino(sinoDescriptor);
1204 :
1205 32 : TestType op(volumeDescriptor, sinoDescriptor);
1206 :
1207 32 : THEN("The result of projecting through a voxel is exactly the voxel's value")
1208 : {
1209 64 : for (index_t j = 0; j < volSize; j++) {
1210 48 : volume = 0;
1211 48 : switch (i) {
1212 48 : case 0:
1213 48 : volume(0, volSize / 2, j) = 1;
1214 48 : break;
1215 0 : case 1:
1216 0 : volume(volSize / 2, 0, j) = 1;
1217 0 : break;
1218 0 : case 2:
1219 0 : volume(0, 0, j) = 1;
1220 0 : break;
1221 0 : default:
1222 0 : break;
1223 : }
1224 :
1225 48 : op.apply(volume, sino);
1226 48 : REQUIRE_EQ(sino[0], Approx(1));
1227 : }
1228 :
1229 32 : AND_THEN("The backprojection yields the exact adjoints")
1230 : {
1231 16 : sino[0] = 1;
1232 16 : op.applyAdjoint(sino, volume);
1233 :
1234 16 : REQUIRE_UNARY(isApprox(
1235 : volume, DataContainer<data_t>(volumeDescriptor, backProj[i])));
1236 : }
1237 : }
1238 : }
1239 : }
1240 :
1241 256 : for (unsigned i = numCases / 2; i < numCases; i++) {
1242 208 : WHEN("A z-axis-aligned ray runs along the corners and edges of the volume")
1243 : {
1244 32 : INFO("A z-axis-aligned ray runs along the ", al[i], " of the volume");
1245 : // x-ray source must be very far from the volume center to make testing of the
1246 : // op backprojection simpler
1247 32 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr,
1248 16 : std::move(volData), std::move(sinoData),
1249 : RotationAngles3D{Gamma{0}}, PrincipalPointOffset2D{0, 0},
1250 16 : RotationOffset3D{-offsetx[i], -offsety[i], 0});
1251 :
1252 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1253 32 : DataContainer<data_t> sino(sinoDescriptor);
1254 :
1255 32 : TestType op(volumeDescriptor, sinoDescriptor);
1256 :
1257 32 : THEN("The result of projecting is zero")
1258 : {
1259 16 : volume = 1;
1260 :
1261 16 : op.apply(volume, sino);
1262 16 : REQUIRE_EQ(sino[0], Approx(0));
1263 :
1264 32 : AND_THEN("The result of backprojection is also zero")
1265 : {
1266 16 : sino[0] = 1;
1267 16 : op.applyAdjoint(sino, volume);
1268 :
1269 16 : DataContainer<data_t> zero(volumeDescriptor);
1270 16 : zero = 0;
1271 16 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
1272 : }
1273 : }
1274 : }
1275 : }
1276 : }
1277 :
1278 224 : GIVEN("An angle of 30 degrees")
1279 : {
1280 64 : data_t sqrt3d = std::sqrt(static_cast<data_t>(3));
1281 :
1282 64 : Eigen::Matrix<data_t, volSize * volSize * volSize, 1> backProj;
1283 :
1284 64 : auto gamma = Gamma{Degree{30}};
1285 :
1286 80 : WHEN("A ray goes through the center of the volume")
1287 : {
1288 : // In this case the ray enters and exits the volume along the main direction
1289 16 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1290 : RotationAngles3D{gamma});
1291 :
1292 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1293 32 : DataContainer<data_t> sino(sinoDescriptor);
1294 :
1295 32 : TestType op(volumeDescriptor, sinoDescriptor);
1296 :
1297 32 : THEN("The ray intersects the correct voxels")
1298 : {
1299 16 : volume = 1;
1300 16 : volume(1, 1, 1) = 0;
1301 16 : volume(2, 1, 0) = 0;
1302 16 : volume(1, 1, 0) = 0;
1303 16 : volume(0, 1, 2) = 0;
1304 16 : volume(1, 1, 2) = 0;
1305 :
1306 16 : op.apply(volume, sino);
1307 16 : REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
1308 :
1309 32 : AND_THEN("The correct weighting is applied")
1310 : {
1311 16 : volume(1, 1, 1) = 1;
1312 16 : volume(2, 1, 0) = 3;
1313 16 : volume(1, 1, 2) = 2;
1314 :
1315 16 : op.apply(volume, sino);
1316 16 : REQUIRE_EQ(sino[0], Approx(3 * sqrt3d - 1).epsilon(epsilon));
1317 :
1318 16 : sino[0] = 1;
1319 : // clang-format off
1320 32 : backProj << 0, 0, 0, 0, 1 - 1 / sqrt3d, sqrt3d - 1, 0, 0, 0,
1321 16 : 0, 0, 0, 0, 2 / sqrt3d, 0, 0, 0, 0,
1322 32 : 0, 0, 0, sqrt3d - 1, 1 - 1 / sqrt3d, 0, 0, 0, 0;
1323 : // clang-format on
1324 :
1325 16 : op.applyAdjoint(sino, volume);
1326 16 : REQUIRE_UNARY(isApprox(
1327 : volume, DataContainer<data_t>(volumeDescriptor, backProj), epsilon));
1328 : }
1329 : }
1330 : }
1331 :
1332 80 : WHEN("A ray with an angle of 30 degrees enters through the right border")
1333 : {
1334 : // In this case the ray enters through a border orthogonal to a non-main
1335 : // direction
1336 16 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1337 : RotationAngles3D{Gamma{pi_t / 6}}, PrincipalPointOffset2D{0, 0},
1338 : RotationOffset3D{1, 0, 0});
1339 :
1340 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1341 32 : DataContainer<data_t> sino(sinoDescriptor);
1342 :
1343 32 : TestType op(volumeDescriptor, sinoDescriptor);
1344 :
1345 32 : THEN("The ray intersects the correct voxels")
1346 : {
1347 16 : volume = 1;
1348 16 : volume(2, 1, 1) = 0;
1349 16 : volume(2, 1, 0) = 0;
1350 16 : volume(2, 1, 2) = 0;
1351 16 : volume(1, 1, 2) = 0;
1352 :
1353 16 : op.apply(volume, sino);
1354 16 : REQUIRE_EQ(sino[0], Approx(0).epsilon(epsilon));
1355 :
1356 32 : AND_THEN("The correct weighting is applied")
1357 : {
1358 16 : volume(2, 1, 0) = 4;
1359 16 : volume(1, 1, 2) = 3;
1360 16 : volume(2, 1, 1) = 1;
1361 :
1362 16 : op.apply(volume, sino);
1363 16 : REQUIRE_EQ(sino[0], Approx(1 - 2 / sqrt3d + 3 * sqrt3d));
1364 :
1365 16 : sino[0] = 1;
1366 :
1367 : // clang-format off
1368 32 : backProj << 0, 0, 0, 0, 0, 1 - 1 / sqrt3d, 0, 0, 0,
1369 16 : 0, 0, 0, 0, 0, 2 / sqrt3d, 0, 0, 0,
1370 32 : 0, 0, 0, 0, sqrt3d - 1, 1 - 1 / sqrt3d, 0, 0, 0;
1371 : // clang-format on
1372 :
1373 16 : op.applyAdjoint(sino, volume);
1374 16 : REQUIRE_UNARY(isApprox(
1375 : volume, DataContainer<data_t>(volumeDescriptor, backProj), epsilon));
1376 : }
1377 : }
1378 : }
1379 :
1380 80 : WHEN("A ray with an angle of 30 degrees exits through the left border")
1381 : {
1382 : // In this case the ray exit through a border orthogonal to a non-main direction
1383 16 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1384 : RotationAngles3D{Gamma{pi_t / 6}}, PrincipalPointOffset2D{0, 0},
1385 : RotationOffset3D{-1, 0, 0});
1386 :
1387 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1388 32 : DataContainer<data_t> sino(sinoDescriptor);
1389 :
1390 32 : TestType op(volumeDescriptor, sinoDescriptor);
1391 :
1392 32 : THEN("The ray intersects the correct voxels")
1393 : {
1394 16 : volume = 1;
1395 16 : volume(0, 1, 0) = 0;
1396 16 : volume(1, 1, 0) = 0;
1397 16 : volume(0, 1, 1) = 0;
1398 16 : volume(0, 1, 2) = 0;
1399 :
1400 16 : op.apply(volume, sino);
1401 16 : REQUIRE_EQ(sino[0], Approx(0).epsilon(epsilon));
1402 :
1403 32 : AND_THEN("The correct weighting is applied")
1404 : {
1405 16 : volume(0, 1, 2) = 4;
1406 16 : volume(1, 1, 0) = 3;
1407 16 : volume(0, 1, 1) = 1;
1408 :
1409 16 : op.apply(volume, sino);
1410 16 : REQUIRE_EQ(sino[0], Approx(3 * sqrt3d + 1 - 2 / sqrt3d).epsilon(epsilon));
1411 :
1412 16 : sino[0] = 1;
1413 :
1414 : // clang-format off
1415 32 : backProj << 0, 0, 0, 1 - 1 / sqrt3d, sqrt3d - 1, 0, 0, 0, 0,
1416 16 : 0, 0, 0, 2 / sqrt3d, 0, 0, 0, 0, 0,
1417 32 : 0, 0, 0, 1 - 1 / sqrt3d, 0, 0, 0, 0, 0;
1418 : // clang-format on
1419 :
1420 16 : op.applyAdjoint(sino, volume);
1421 16 : REQUIRE_UNARY(isApprox(
1422 : volume, DataContainer<data_t>(volumeDescriptor, backProj), epsilon));
1423 : }
1424 : }
1425 : }
1426 :
1427 80 : WHEN("A ray with an angle of 30 degrees only intersects a single voxel")
1428 : {
1429 : // special case - no interior voxels, entry and exit voxels are the same
1430 16 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1431 : RotationAngles3D{Gamma{pi_t / 6}}, PrincipalPointOffset2D{0, 0},
1432 : RotationOffset3D{-2, 0, 0});
1433 :
1434 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1435 32 : DataContainer<data_t> sino(sinoDescriptor);
1436 :
1437 32 : TestType op(volumeDescriptor, sinoDescriptor);
1438 :
1439 32 : THEN("The ray intersects the correct voxels")
1440 : {
1441 16 : volume = 1;
1442 16 : volume(0, 1, 0) = 0;
1443 :
1444 16 : op.apply(volume, sino);
1445 16 : REQUIRE_EQ(sino[0], Approx(0).epsilon(epsilon));
1446 :
1447 32 : AND_THEN("The correct weighting is applied")
1448 : {
1449 16 : volume(0, 1, 0) = 1;
1450 :
1451 16 : op.apply(volume, sino);
1452 16 : REQUIRE_EQ(sino[0], Approx(sqrt3d - 1).epsilon(epsilon));
1453 :
1454 16 : sino[0] = 1;
1455 :
1456 : // clang-format off
1457 32 : backProj << 0, 0, 0, sqrt3d - 1, 0, 0, 0, 0, 0,
1458 16 : 0, 0, 0, 0, 0, 0, 0, 0, 0,
1459 32 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1460 : // clang-format on
1461 :
1462 16 : op.applyAdjoint(sino, volume);
1463 16 : REQUIRE_UNARY(isApprox(
1464 : volume, DataContainer<data_t>(volumeDescriptor, backProj), epsilon));
1465 : }
1466 : }
1467 : }
1468 : }
1469 : }
1470 :
1471 224 : GIVEN("Given a 5x5x5 volume")
1472 : {
1473 : // Domain setup
1474 32 : const index_t volSize = 5;
1475 :
1476 64 : IndexVector_t volumeDims(3);
1477 32 : volumeDims << volSize, volSize, volSize;
1478 :
1479 64 : VolumeDescriptor volumeDescriptor(volumeDims);
1480 64 : DataContainer<data_t> volume(volumeDescriptor);
1481 :
1482 : // range setup
1483 32 : const index_t detectorSize = 1;
1484 32 : const index_t numImgs = 1;
1485 :
1486 64 : IndexVector_t sinoDims(3);
1487 32 : sinoDims << detectorSize, detectorSize, numImgs;
1488 :
1489 : // Setup geometry
1490 32 : auto stc = SourceToCenterOfRotation{20 * volSize};
1491 32 : auto ctr = CenterOfRotationToDetector{volSize};
1492 96 : auto volData = VolumeData3D{Size3D{volumeDims}};
1493 96 : auto sinoData = SinogramData3D{Size3D{sinoDims}};
1494 :
1495 64 : std::vector<Geometry> geom;
1496 :
1497 64 : GIVEN("Tracing rays along axis")
1498 : {
1499 32 : const index_t numCases = 9;
1500 : using Array = std::array<real_t, numCases>;
1501 : using StrArray = std::array<std::string, numCases>;
1502 :
1503 32 : Array alpha = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1504 32 : Array beta = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, pi_t / 2, pi_t / 2, pi_t / 2};
1505 32 : Array gamma = {0.0, 0.0, 0.0, pi_t / 2, pi_t / 2,
1506 : pi_t / 2, pi_t / 2, pi_t / 2, pi_t / 2};
1507 :
1508 32 : Array offsetx = {volSize, 0.0, volSize, 0.0, 0.0, 0.0, volSize, 0.0, volSize};
1509 32 : Array offsety = {0.0, volSize, volSize, volSize, 0.0, volSize, 0.0, 0.0, 0.0};
1510 32 : Array offsetz = {0.0, 0.0, 0.0, 0.0, volSize, volSize, 0.0, volSize, volSize};
1511 :
1512 64 : StrArray neg = {"x", "y", "x and y", "y", "z", "y and z", "x", "z", "x and z"};
1513 64 : StrArray ali = {"z", "z", "z", "x", "x", "x", "y", "y", "y"};
1514 :
1515 320 : for (unsigned i = 0; i < numCases; i++) {
1516 304 : WHEN("Tracing along a fixed axis-aligned ray with negative coordinate of origin")
1517 : {
1518 32 : INFO("Tracing along a ", ali[i], "-axis-aligned ray with negative ", neg[i],
1519 : "-coodinate of origin");
1520 32 : geom.emplace_back(
1521 16 : stc, ctr, std::move(volData), std::move(sinoData),
1522 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}, Alpha{alpha[i]}},
1523 : PrincipalPointOffset2D{0, 0},
1524 16 : RotationOffset3D{-offsetx[i], -offsety[i], -offsetz[i]});
1525 :
1526 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1527 32 : DataContainer<data_t> sino(sinoDescriptor);
1528 :
1529 32 : TestType op(volumeDescriptor, sinoDescriptor);
1530 :
1531 32 : THEN("Result of forward projection is zero")
1532 : {
1533 16 : op.apply(volume, sino);
1534 32 : DataContainer<data_t> zero(sinoDescriptor);
1535 16 : zero = 0;
1536 16 : REQUIRE_UNARY(isApprox(sino, zero, epsilon));
1537 :
1538 32 : AND_THEN("Result of backprojection is zero")
1539 : {
1540 16 : op.applyAdjoint(sino, volume);
1541 16 : DataContainer<data_t> zero(volumeDescriptor);
1542 16 : zero = 0;
1543 16 : REQUIRE_UNARY(isApprox(volume, zero, epsilon));
1544 : }
1545 : }
1546 : }
1547 : }
1548 : }
1549 : }
1550 192 : }
1551 :
1552 116 : TEST_CASE_TEMPLATE("SiddonsMethodCUDA: 3D setup with a multiple rays", TestType,
1553 : SiddonsMethodCUDA<float>, SiddonsMethodCUDA<double>, SiddonsMethod<float>,
1554 : SiddonsMethod<double>)
1555 : {
1556 : // Turn logger of
1557 16 : Logger::setLevel(Logger::LogLevel::OFF);
1558 :
1559 : using data_t = decltype(return_data_t(std::declval<TestType>()));
1560 :
1561 32 : GIVEN("A 3D setting with multiple projection angles")
1562 : {
1563 : // Domain setup
1564 16 : const index_t volSize = 3;
1565 :
1566 32 : IndexVector_t volumeDims(3);
1567 16 : volumeDims << volSize, volSize, volSize;
1568 :
1569 32 : VolumeDescriptor volumeDescriptor(volumeDims);
1570 32 : DataContainer<data_t> volume(volumeDescriptor);
1571 :
1572 : // range setup
1573 16 : const index_t detectorSize = 1;
1574 16 : const index_t numImgs = 6;
1575 :
1576 32 : IndexVector_t sinoDims(3);
1577 16 : sinoDims << detectorSize, detectorSize, numImgs;
1578 :
1579 : // Setup geometry
1580 16 : auto stc = SourceToCenterOfRotation{20 * volSize};
1581 16 : auto ctr = CenterOfRotationToDetector{volSize};
1582 48 : auto volData = VolumeData3D{Size3D{volumeDims}};
1583 48 : auto sinoData = SinogramData3D{Size3D{sinoDims}};
1584 :
1585 32 : std::vector<Geometry> geom;
1586 :
1587 32 : WHEN("x-, y and z-axis-aligned rays are present")
1588 : {
1589 16 : real_t beta[numImgs] = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
1590 16 : real_t gamma[numImgs] = {0.0, pi_t, pi_t / 2, 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
1591 :
1592 112 : for (index_t i = 0; i < numImgs; i++)
1593 96 : geom.emplace_back(stc, ctr, VolumeData3D{Size3D{volumeDims}},
1594 : SinogramData3D{Size3D{sinoDims}},
1595 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
1596 :
1597 32 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1598 32 : DataContainer<data_t> sino(sinoDescriptor);
1599 :
1600 32 : TestType op(volumeDescriptor, sinoDescriptor);
1601 :
1602 32 : THEN("Values are accumulated correctly along each ray's path")
1603 : {
1604 16 : volume = 0;
1605 :
1606 : // set only values along the rays' path to one to make sure interpolation is dones
1607 : // correctly
1608 64 : for (index_t i = 0; i < volSize; i++) {
1609 48 : volume(i, volSize / 2, volSize / 2) = 1;
1610 48 : volume(volSize / 2, i, volSize / 2) = 1;
1611 48 : volume(volSize / 2, volSize / 2, i) = 1;
1612 : }
1613 :
1614 16 : op.apply(volume, sino);
1615 112 : for (index_t i = 0; i < numImgs; i++)
1616 96 : REQUIRE_EQ(sino[i], Approx(3.0));
1617 :
1618 32 : AND_THEN("Backprojection yields the exact adjoint")
1619 : {
1620 16 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> cmp(volSize * volSize * volSize);
1621 :
1622 : // clang-format off
1623 32 : cmp << 0, 0, 0, 0, 6, 0, 0, 0, 0,
1624 16 : 0, 6, 0, 6, 18, 6, 0, 6, 0,
1625 32 : 0, 0, 0, 0, 6, 0, 0, 0, 0;
1626 : // clang-format on
1627 :
1628 16 : op.applyAdjoint(sino, volume);
1629 16 : REQUIRE_UNARY(isApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
1630 : }
1631 : }
1632 : }
1633 : }
1634 16 : }
|
