Line data Source code
1 : #include "doctest/doctest.h"
2 :
3 : #include "JosephsMethodCUDA.h"
4 : #include "Geometry.h"
5 : #include "VolumeDescriptor.h"
6 : #include "Logger.h"
7 : #include "PlanarDetectorDescriptor.h"
8 : #include "testHelpers.h"
9 :
10 : #include <array>
11 :
12 : using namespace elsa;
13 : using namespace elsa::geometry;
14 : using namespace doctest;
15 :
16 : // TODO(dfrank): remove this and replace with checkApproxEq
17 : using doctest::Approx;
18 :
19 : /*
20 : * this function declaration can be used in conjunction with decltype to deduce the
21 : * template parameter of a templated class at compile time
22 : *
23 : * the template parameter must be a typename
24 : */
25 : template <template <typename> typename T, typename data_t>
26 : constexpr data_t return_data_t(const T<data_t>&);
27 :
28 34 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: Calls to functions of super class", TestType,
29 : JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
30 : {
31 : // Turn logger of
32 4 : Logger::setLevel(Logger::LogLevel::OFF);
33 :
34 : using data_t = decltype(return_data_t(std::declval<TestType>()));
35 8 : GIVEN("A projector")
36 : {
37 8 : IndexVector_t volumeDims(2), sinoDims(2);
38 4 : const index_t volSize = 10;
39 4 : const index_t detectorSize = 10;
40 4 : const index_t numImgs = 10;
41 4 : volumeDims << volSize, volSize;
42 4 : sinoDims << detectorSize, numImgs;
43 8 : VolumeDescriptor volumeDescriptor(volumeDims);
44 8 : DataContainer<data_t> volume(volumeDescriptor);
45 4 : volume = 0;
46 :
47 4 : auto stc = SourceToCenterOfRotation{20 * volSize};
48 4 : auto ctr = CenterOfRotationToDetector{volSize};
49 :
50 8 : std::vector<Geometry> geom;
51 44 : for (index_t i = 0; i < numImgs; i++) {
52 40 : real_t angle = static_cast<real_t>(i) * 2 * pi_t / 10;
53 40 : geom.emplace_back(stc, ctr, Radian{angle}, VolumeData2D{Size2D{volumeDims}},
54 : SinogramData2D{Size2D{sinoDims}});
55 : }
56 :
57 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
58 8 : DataContainer<data_t> sino(sinoDescriptor);
59 4 : sino = 0;
60 :
61 8 : TestType fast(volumeDescriptor, sinoDescriptor);
62 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
63 :
64 8 : WHEN("Projector is cloned")
65 : {
66 8 : auto fastClone = fast.clone();
67 8 : auto slowClone = slow.clone();
68 8 : auto sinoClone = DataContainer<data_t>(sinoDescriptor);
69 4 : sinoClone = 0;
70 8 : auto volumeClone = DataContainer<data_t>(volumeDescriptor);
71 4 : volumeClone = 0;
72 :
73 8 : THEN("Results do not change (may still be slightly different due to summation being "
74 : "performed in a different order)")
75 : {
76 4 : fast.apply(volume, sino);
77 4 : fastClone->apply(volume, sinoClone);
78 4 : REQUIRE_UNARY(isCwiseApprox(sino, sinoClone));
79 :
80 4 : slowClone->apply(volume, sinoClone);
81 4 : REQUIRE_UNARY(isCwiseApprox(sino, sinoClone));
82 :
83 4 : fast.applyAdjoint(sino, volume);
84 4 : fastClone->applyAdjoint(sino, volumeClone);
85 4 : REQUIRE_UNARY(isCwiseApprox(volume, volumeClone));
86 :
87 4 : slow.applyAdjoint(sino, volume);
88 4 : slowClone->applyAdjoint(sino, volumeClone);
89 4 : REQUIRE_UNARY(isCwiseApprox(volume, volumeClone));
90 : }
91 : }
92 : }
93 4 : }
94 :
95 110 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 2D setup with a single ray", TestType,
96 : JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
97 : {
98 : // Turn logger of
99 80 : Logger::setLevel(Logger::LogLevel::OFF);
100 :
101 : using data_t = decltype(return_data_t(std::declval<TestType>()));
102 :
103 128 : GIVEN("A 5x5 volume and detector of size 1 with 1 angle")
104 : {
105 : // Domain setup
106 48 : const index_t volSize = 5;
107 :
108 96 : IndexVector_t volumeDims(2);
109 48 : volumeDims << volSize, volSize;
110 :
111 96 : VolumeDescriptor volumeDescriptor(volumeDims);
112 96 : DataContainer<data_t> volume(volumeDescriptor);
113 :
114 : // range setup
115 48 : const index_t detectorSize = 1;
116 48 : const index_t numImgs = 1;
117 :
118 96 : IndexVector_t sinoDims(2);
119 48 : sinoDims << detectorSize, numImgs;
120 :
121 : // Setup geometry
122 48 : auto stc = SourceToCenterOfRotation{20 * volSize};
123 48 : auto ctr = CenterOfRotationToDetector{volSize};
124 144 : auto volData = VolumeData2D{Size2D{volumeDims}};
125 144 : auto sinoData = SinogramData2D{Size2D{sinoDims}};
126 :
127 96 : std::vector<Geometry> geom;
128 :
129 56 : GIVEN("A basic geometry setup")
130 : {
131 8 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData));
132 :
133 16 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
134 16 : DataContainer<data_t> sino(sinoDescriptor);
135 :
136 16 : TestType slow(volumeDescriptor, sinoDescriptor, false);
137 16 : TestType fast(volumeDescriptor, sinoDescriptor);
138 :
139 12 : WHEN("Sinogram conatainer is not zero initialized and we project through an empty "
140 : "volume")
141 : {
142 4 : volume = 0;
143 4 : sino = 1;
144 :
145 8 : THEN("Result is zero")
146 : {
147 4 : fast.apply(volume, sino);
148 4 : DataContainer<data_t> zero(sinoDescriptor);
149 4 : zero = 0;
150 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
151 :
152 4 : sino = 1;
153 4 : slow.apply(volume, sino);
154 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
155 : }
156 : }
157 :
158 12 : WHEN("Volume container is not zero initialized and we backproject from an empty "
159 : "sinogram")
160 : {
161 4 : sino = 0;
162 4 : volume = 1;
163 :
164 8 : THEN("Result is zero")
165 : {
166 4 : fast.applyAdjoint(sino, volume);
167 4 : DataContainer<data_t> zero(volumeDescriptor);
168 4 : zero = 0;
169 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
170 :
171 4 : volume = 1;
172 4 : slow.applyAdjoint(sino, volume);
173 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
174 : }
175 : }
176 : }
177 :
178 64 : GIVEN("Rays not intersecting the bounding box are present")
179 : {
180 20 : WHEN("Tracing along a y-axis-aligned ray with a negative x-coordinate of origin")
181 : {
182 4 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
183 : PrincipalPointOffset{}, RotationOffset2D{volSize, 0});
184 :
185 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
186 8 : DataContainer<data_t> sino(sinoDescriptor);
187 :
188 8 : TestType fast(volumeDescriptor, sinoDescriptor);
189 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
190 :
191 8 : THEN("Result of forward projection is zero")
192 : {
193 4 : fast.apply(volume, sino);
194 8 : DataContainer<data_t> zero(sinoDescriptor);
195 4 : zero = 0;
196 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
197 :
198 4 : slow.apply(volume, sino);
199 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
200 :
201 8 : AND_THEN("Result of backprojection is zero")
202 : {
203 4 : fast.applyAdjoint(sino, volume);
204 4 : DataContainer<data_t> zero(volumeDescriptor);
205 4 : zero = 0;
206 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
207 :
208 4 : slow.applyAdjoint(sino, volume);
209 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
210 : }
211 : }
212 : }
213 :
214 20 : WHEN("Tracing along a y-axis-aligned ray with a x-coordinate of origin beyond the "
215 : "bounding "
216 : "box")
217 : {
218 4 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
219 : PrincipalPointOffset{}, RotationOffset2D{-volSize, 0});
220 :
221 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
222 8 : DataContainer<data_t> sino(sinoDescriptor);
223 :
224 8 : TestType fast(volumeDescriptor, sinoDescriptor);
225 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
226 :
227 8 : THEN("Result of forward projection is zero")
228 : {
229 8 : DataContainer<data_t> zero(sinoDescriptor);
230 4 : zero = 0;
231 :
232 4 : fast.apply(volume, sino);
233 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
234 :
235 4 : slow.apply(volume, sino);
236 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
237 :
238 8 : AND_THEN("Result of backprojection is zero")
239 : {
240 4 : DataContainer<data_t> zero(volumeDescriptor);
241 4 : zero = 0;
242 :
243 4 : fast.applyAdjoint(sino, volume);
244 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
245 :
246 4 : slow.applyAdjoint(sino, volume);
247 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
248 : }
249 : }
250 : }
251 :
252 20 : WHEN("Tracing along a x-axis-aligned ray with a negative y-coordinate of origin")
253 : {
254 8 : geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData),
255 4 : std::move(sinoData), PrincipalPointOffset{},
256 : RotationOffset2D{0, volSize});
257 :
258 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
259 8 : DataContainer<data_t> sino(sinoDescriptor);
260 :
261 8 : TestType fast(volumeDescriptor, sinoDescriptor);
262 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
263 :
264 8 : THEN("Result of forward projection is zero")
265 : {
266 8 : DataContainer<data_t> zero(sinoDescriptor);
267 4 : zero = 0;
268 :
269 4 : fast.apply(volume, sino);
270 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
271 :
272 4 : slow.apply(volume, sino);
273 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
274 :
275 8 : AND_THEN("Result of backprojection is zero")
276 : {
277 4 : DataContainer<data_t> zero(volumeDescriptor);
278 4 : zero = 0;
279 :
280 4 : fast.applyAdjoint(sino, volume);
281 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
282 :
283 4 : slow.applyAdjoint(sino, volume);
284 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
285 : }
286 : }
287 : }
288 :
289 20 : WHEN("Tracing along a x-axis-aligned ray with a y-coordinate of origin beyond the "
290 : "bounding "
291 : "box")
292 : {
293 8 : geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData),
294 4 : std::move(sinoData), PrincipalPointOffset{},
295 : RotationOffset2D{0, -volSize});
296 :
297 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
298 8 : DataContainer<data_t> sino(sinoDescriptor);
299 :
300 8 : TestType fast(volumeDescriptor, sinoDescriptor);
301 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
302 :
303 8 : THEN("Result of forward projection is zero")
304 : {
305 8 : DataContainer<data_t> zero(sinoDescriptor);
306 4 : zero = 0;
307 :
308 4 : fast.apply(volume, sino);
309 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
310 :
311 4 : slow.apply(volume, sino);
312 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
313 :
314 8 : AND_THEN("Result of backprojection is zero")
315 : {
316 4 : DataContainer<data_t> zero(volumeDescriptor);
317 4 : zero = 0;
318 :
319 4 : fast.applyAdjoint(sino, volume);
320 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
321 :
322 4 : slow.applyAdjoint(sino, volume);
323 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
324 : }
325 : }
326 : }
327 : }
328 :
329 : // Expected results
330 288 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProjections[2];
331 48 : backProjections[0].resize(volSize * volSize);
332 48 : backProjections[1].resize(volSize * volSize);
333 :
334 48 : constexpr index_t numCases = 4;
335 48 : const std::array<real_t, numCases> angles = {0., 90., 180., 270.};
336 :
337 : // Testing for axis-aligned rays
338 56 : GIVEN("Given rays along the axes")
339 : {
340 :
341 16 : backProjections[0] << 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,
342 16 : 1, 0, 0;
343 16 : backProjections[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
344 16 : 0, 0, 0;
345 :
346 40 : for (std::size_t i = 0; i < numCases; i++) {
347 36 : WHEN("An axis-aligned ray with a fixed angle passes through the center of a pixel")
348 : {
349 8 : INFO("An axis-aligned ray with an angle of ", angles[i],
350 : " radians passes through the center of a pixel");
351 4 : geom.emplace_back(stc, ctr, Degree{angles[i]}, std::move(volData),
352 4 : std::move(sinoData));
353 :
354 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
355 8 : DataContainer<data_t> sino(sinoDescriptor);
356 :
357 8 : TestType fast(volumeDescriptor, sinoDescriptor);
358 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
359 :
360 8 : THEN("The result of projecting through a pixel is exactly the pixel value")
361 : {
362 24 : for (index_t j = 0; j < volSize; j++) {
363 20 : volume = 0;
364 20 : if (i % 2 == 0)
365 20 : volume(volSize / 2, j) = 1;
366 : else
367 0 : volume(j, volSize / 2) = 1;
368 :
369 20 : fast.apply(volume, sino);
370 : // Using doubles significantly increases interpolation accuracy
371 : // For example: when using floats, points very near the border (less
372 : // than ~1/500th of a pixel away from the border) are rounded to
373 : // actually lie on the border. This then yields more accurate results
374 : // when using floats in some of the axis-aligned test cases, despite the
375 : // lower interpolation accuracy.
376 : //
377 : // => different requirements for floats and doubles, looser
378 : // requirements for doubles
379 20 : REQUIRE_EQ(sino[0], Approx(1.0));
380 :
381 20 : slow.apply(volume, sino);
382 20 : REQUIRE_EQ(sino[0], Approx(1.0));
383 : }
384 :
385 8 : AND_THEN(
386 : "The backprojection sets the values of all hit pixels to the detector "
387 : "value")
388 : {
389 4 : fast.applyAdjoint(sino, volume);
390 4 : REQUIRE_UNARY(isCwiseApprox(
391 : volume,
392 : DataContainer<data_t>(volumeDescriptor, backProjections[i % 2])));
393 :
394 4 : slow.applyAdjoint(sino, volume);
395 4 : REQUIRE_UNARY(isCwiseApprox(
396 : volume,
397 : DataContainer<data_t>(volumeDescriptor, backProjections[i % 2])));
398 : }
399 : }
400 : }
401 : }
402 : }
403 :
404 : // Move the source far far back, so rays become nearly parallel
405 48 : stc = SourceToCenterOfRotation{volSize * 2000};
406 :
407 56 : GIVEN("Off center rays along axes")
408 : {
409 8 : std::array<real_t, numCases> offsetx = {-0.25, 0.0, -0.25, 0.0};
410 8 : std::array<real_t, numCases> offsety = {0.0, -0.25, 0.0, -0.25};
411 :
412 16 : backProjections[0] << 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0, 0,
413 16 : 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0;
414 :
415 16 : backProjections[1] << 0, 0, 0, 0, 0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.75, 0.75, 0.75,
416 16 : 0.75, 0.75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
417 :
418 40 : for (std::size_t i = 0; i < numCases; i++) {
419 36 : WHEN("An axis-aligned ray with a fixed angle, which does not pass through the "
420 : "center of a pixel")
421 : {
422 8 : INFO("An axis-aligned ray with an angle of ", angles[i],
423 : " radians not passing through the center of a pixel");
424 8 : geom.emplace_back(stc, ctr, Degree{angles[i]}, std::move(volData),
425 4 : std::move(sinoData), PrincipalPointOffset{0},
426 : RotationOffset2D{offsetx[i], offsety[i]});
427 :
428 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
429 8 : DataContainer<data_t> sino(sinoDescriptor);
430 :
431 8 : TestType fast(volumeDescriptor, sinoDescriptor);
432 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
433 :
434 8 : THEN("The result of projecting through a pixel is the interpolated value "
435 : "between "
436 : "the two pixels closest to the ray")
437 : {
438 24 : for (index_t j = 0; j < volSize; j++) {
439 20 : volume = 0;
440 20 : if (i % 2 == 0)
441 20 : volume(volSize / 2, j) = 1;
442 : else
443 0 : volume(j, volSize / 2) = 1;
444 :
445 20 : fast.apply(volume, sino);
446 20 : REQUIRE_EQ(sino[0], Approx(0.75));
447 :
448 20 : slow.apply(volume, sino);
449 20 : REQUIRE_EQ(sino[0], Approx(0.75));
450 : }
451 :
452 8 : AND_THEN("The slow backprojection yields the exact adjoint, the fast "
453 : "backprojection "
454 : "also yields the exact adjoint for a very distant x-ray source")
455 : {
456 4 : sino[0] = 1;
457 4 : slow.applyAdjoint(sino, volume);
458 4 : REQUIRE_UNARY(isCwiseApprox(
459 : volume,
460 : DataContainer<data_t>(volumeDescriptor, backProjections[i % 2])));
461 :
462 4 : fast.applyAdjoint(sino, volume);
463 :
464 : // Using doubles significantly increases interpolation accuracy
465 : // For example: when using floats, points very near the border (less
466 : // than
467 : // ~1/500th of a pixel away from the border) are rounded to actually lie
468 : // on the border. This then yields more accurate results when using
469 : // floats in some of the axis-aligned test cases, despite the lower
470 : // interpolation accuracy.
471 : //
472 : // => different requirements for floats and doubles, looser
473 : // requirements
474 : // for doubles
475 : if constexpr (std::is_same_v<data_t, float>)
476 2 : REQUIRE_UNARY(isCwiseApprox(
477 : volume, DataContainer<data_t>(volumeDescriptor,
478 : backProjections[i % 2])));
479 : else
480 2 : REQUIRE_UNARY(isApprox(
481 : volume,
482 : DataContainer<data_t>(volumeDescriptor, backProjections[i % 2]),
483 : static_cast<real_t>(0.001)));
484 : }
485 : }
486 : }
487 : }
488 : }
489 :
490 52 : GIVEN("Rays going running allong right bolume boundary")
491 : {
492 8 : backProjections[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
493 8 : 0, 0, 1;
494 :
495 8 : WHEN("A y-axis-aligned ray runs along the right volume boundary")
496 : {
497 4 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
498 : PrincipalPointOffset{0}, RotationOffset2D{volSize * 0.5, 0});
499 :
500 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
501 8 : DataContainer<data_t> sino(sinoDescriptor);
502 :
503 8 : TestType fast(volumeDescriptor, sinoDescriptor);
504 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
505 :
506 8 : THEN("The result of projecting through a pixel is exactly the pixel's value (we "
507 : "mirror "
508 : "values at the border for the purpose of interpolation)")
509 : {
510 24 : for (index_t j = 0; j < volSize; j++) {
511 20 : volume = 0;
512 20 : volume(volSize - 1, j) = 1;
513 :
514 20 : fast.apply(volume, sino);
515 20 : REQUIRE_EQ(sino[0], Approx(1));
516 :
517 20 : slow.apply(volume, sino);
518 20 : REQUIRE_EQ(sino[0], Approx(1));
519 : }
520 :
521 8 : AND_THEN(
522 : "The slow backprojection yields the exact adjoint, the fast backprojection "
523 : "also yields the exact adjoint for a very distant x-ray source")
524 : {
525 4 : sino[0] = 1;
526 4 : slow.applyAdjoint(sino, volume);
527 4 : REQUIRE_UNARY(isCwiseApprox(
528 : volume, DataContainer<data_t>(volumeDescriptor, backProjections[0])));
529 :
530 4 : fast.applyAdjoint(sino, volume);
531 : // Using doubles significantly increases interpolation accuracy
532 : // For example: when using floats, points very near the border (less than
533 : // ~1/500th of a pixel away from the border) are rounded to actually lie on
534 : // the border. This then yields more accurate results when using floats in
535 : // some of the axis-aligned test cases, despite the lower interpolation
536 : // accuracy.
537 : //
538 : // => different requirements for floats and doubles, looser requirements
539 : // for doubles
540 : if constexpr (std::is_same_v<data_t, float>)
541 2 : REQUIRE_UNARY(isCwiseApprox(
542 : volume, DataContainer<data_t>(volumeDescriptor,
543 : (backProjections[0] / 2).eval())));
544 : else
545 2 : REQUIRE_UNARY(
546 : isApprox(volume,
547 : DataContainer<data_t>(volumeDescriptor,
548 : (backProjections[0] / 2).eval()),
549 : static_cast<real_t>(0.001)));
550 : }
551 : }
552 : }
553 : }
554 :
555 52 : GIVEN("Rays going running along left bolume boundary")
556 : {
557 8 : backProjections[0] << 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0,
558 8 : 0, 0, 0;
559 :
560 8 : WHEN("A y-axis-aligned ray runs along the left volume boundary")
561 : {
562 4 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
563 : PrincipalPointOffset{0}, RotationOffset2D{-volSize * 0.5, 0});
564 :
565 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
566 8 : DataContainer<data_t> sino(sinoDescriptor);
567 :
568 8 : TestType fast(volumeDescriptor, sinoDescriptor);
569 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
570 :
571 : // TestType fast(volumeDescriptor, sinoDescriptor, geom);
572 : // TestType slow(volumeDescriptor, sinoDescriptor, geom, false);
573 8 : THEN("The result of projecting through a pixel is exactly the pixel's value (we "
574 : "mirror "
575 : "values at the border for the purpose of interpolation)")
576 : {
577 24 : for (index_t j = 0; j < volSize; j++) {
578 20 : volume = 0;
579 20 : volume(0, j) = 1;
580 :
581 20 : fast.apply(volume, sino);
582 20 : REQUIRE_EQ(sino[0], Approx(1));
583 :
584 20 : slow.apply(volume, sino);
585 20 : REQUIRE_EQ(sino[0], Approx(1));
586 : }
587 :
588 8 : AND_THEN(
589 : "The slow backprojection yields the exact adjoint, the fast backprojection "
590 : "also yields the exact adjoint for a very distant x-ray source")
591 : {
592 4 : sino[0] = 1;
593 4 : slow.applyAdjoint(sino, volume);
594 4 : REQUIRE_UNARY(isCwiseApprox(
595 : volume, DataContainer<data_t>(volumeDescriptor, backProjections[0])));
596 :
597 4 : fast.applyAdjoint(sino, volume);
598 : /** Using doubles significantly increases interpolation accuracy
599 : * For example: when using floats, points very near the border (less than
600 : * ~1/500th of a pixel away from the border) are rounded to actually lie on
601 : * the border. This then yields more accurate results when using floats in
602 : * some of the axis-aligned test cases, despite the lower interpolation
603 : * accuracy.
604 : *
605 : * => different requirements for floats and doubles, looser requirements
606 : * for doubles
607 : */
608 : if constexpr (std::is_same_v<data_t, float>)
609 2 : REQUIRE_UNARY(isCwiseApprox(
610 : volume, DataContainer<data_t>(volumeDescriptor,
611 : (backProjections[0] / 2).eval())));
612 : else
613 2 : REQUIRE_UNARY(
614 : isApprox(volume,
615 : DataContainer<data_t>(volumeDescriptor,
616 : (backProjections[0] / 2).eval()),
617 : static_cast<real_t>(0.001)));
618 : }
619 : }
620 : }
621 : }
622 : }
623 :
624 112 : GIVEN("A 4x4 volume and detector of size 1 with 1 angle at -30 degree")
625 : {
626 : // Domain setup
627 32 : const index_t volSize = 4;
628 :
629 64 : IndexVector_t volumeDims(2);
630 32 : volumeDims << volSize, volSize;
631 :
632 64 : VolumeDescriptor volumeDescriptor(volumeDims);
633 64 : DataContainer<data_t> volume(volumeDescriptor);
634 :
635 : // range setup
636 32 : const index_t detectorSize = 1;
637 32 : const index_t numImgs = 1;
638 :
639 64 : IndexVector_t sinoDims(2);
640 32 : sinoDims << detectorSize, numImgs;
641 :
642 : // Setup geometry
643 32 : const auto stc = SourceToCenterOfRotation{20 * volSize};
644 32 : const auto ctr = CenterOfRotationToDetector{volSize};
645 96 : auto volData = VolumeData2D{Size2D{volumeDims}};
646 96 : auto sinoData = SinogramData2D{Size2D{sinoDims}};
647 :
648 64 : std::vector<Geometry> geom;
649 :
650 32 : const real_t sqrt3r = std::sqrt(static_cast<real_t>(3));
651 32 : const data_t sqrt3d = std::sqrt(static_cast<data_t>(3));
652 32 : const data_t halfd = static_cast<data_t>(0.5);
653 32 : const data_t thirdd = static_cast<data_t>(1.0 / 3);
654 :
655 48 : GIVEN("An angle of -30 degrees")
656 : {
657 16 : const auto angle = Degree{-30};
658 :
659 20 : WHEN("Projections goes through center of volume")
660 : {
661 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData));
662 :
663 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
664 8 : DataContainer<data_t> sino(sinoDescriptor);
665 4 : sino = 0;
666 :
667 8 : TestType fast(volumeDescriptor, sinoDescriptor);
668 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
669 :
670 4 : real_t weight = 2 / sqrt3r;
671 :
672 8 : THEN("Ray intersects the correct pixels")
673 : {
674 4 : volume = 1;
675 4 : volume(3, 0) = 0;
676 4 : volume(1, 1) = 0;
677 4 : volume(1, 2) = 0;
678 4 : volume(1, 3) = 0;
679 :
680 4 : volume(2, 0) = 0;
681 4 : volume(2, 1) = 0;
682 4 : volume(2, 2) = 0;
683 4 : volume(0, 3) = 0;
684 :
685 4 : fast.apply(volume, sino);
686 8 : DataContainer<data_t> zero(sinoDescriptor);
687 4 : zero = 0;
688 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
689 :
690 4 : slow.apply(volume, sino);
691 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
692 :
693 8 : AND_THEN("The correct weighting is applied")
694 : {
695 4 : volume(3, 0) = 1;
696 4 : volume(1, 1) = 1;
697 4 : volume(1, 2) = 1;
698 4 : volume(1, 3) = 1;
699 :
700 4 : fast.apply(volume, sino);
701 4 : REQUIRE_EQ(sino[0], Approx(2 * weight));
702 :
703 4 : slow.apply(volume, sino);
704 4 : REQUIRE_EQ(sino[0], Approx(2 * weight));
705 :
706 4 : sino[0] = 1;
707 :
708 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
709 8 : slowExpected << 0, 0, (3 - sqrt3d) / 2, (sqrt3d - 1) / 2, 0,
710 4 : (sqrt3d - 1) / (2 * sqrt3d), (sqrt3d + 1) / (2 * sqrt3d), 0, 0,
711 4 : (sqrt3d + 1) / (2 * sqrt3d), (sqrt3d - 1) / (2 * sqrt3d), 0,
712 8 : (sqrt3d - 1) / 2, (3 - sqrt3d) / 2, 0, 0;
713 :
714 4 : slowExpected *= weight;
715 4 : slow.applyAdjoint(sino, volume);
716 4 : REQUIRE_UNARY(isCwiseApprox(
717 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
718 :
719 4 : fast.applyAdjoint(sino, volume);
720 20 : for (real_t i = 0.5; i < volSize; i += 1) {
721 80 : for (real_t j = 0.5; j < volSize; j += 1) {
722 : const real_t angle =
723 64 : std::abs(std::atan((sqrt3r * volSize * 10
724 64 : - static_cast<real_t>(volSize / 2.0) + j)
725 : / (volSize * 10
726 64 : + static_cast<real_t>(volSize / 2.0) - i))
727 : - pi_t / 3);
728 64 : const real_t len = volSize * 21 * std::tan(angle);
729 64 : if (len < 1) {
730 32 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
731 : Approx(1 - len).epsilon(0.005));
732 : } else {
733 32 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
734 : }
735 : }
736 : }
737 : }
738 : }
739 : }
740 :
741 20 : WHEN("Ray enters through the right border")
742 : {
743 : // In this case the ray exits through a border along the main ray direction, but
744 : // enters through a border not along the main direction First pixel should be
745 : // weighted differently
746 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
747 : PrincipalPointOffset{0}, RotationOffset2D{sqrt3r, 0});
748 :
749 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
750 8 : DataContainer<data_t> sino(sinoDescriptor);
751 4 : sino = 0;
752 :
753 8 : TestType fast(volumeDescriptor, sinoDescriptor);
754 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
755 :
756 : // TestType fast(volumeDescriptor, sinoDescriptor, geom);
757 : // TestType slow(volumeDescriptor, sinoDescriptor, geom, false);
758 :
759 8 : THEN("Ray intersects the correct pixels")
760 : {
761 4 : volume = 1;
762 4 : volume(3, 1) = 0;
763 4 : volume(3, 2) = 0;
764 4 : volume(3, 3) = 0;
765 4 : volume(2, 3) = 0;
766 4 : volume(2, 2) = 0;
767 :
768 4 : fast.apply(volume, sino);
769 8 : DataContainer<data_t> zero(sinoDescriptor);
770 4 : zero = 0;
771 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
772 :
773 4 : slow.apply(volume, sino);
774 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
775 :
776 8 : AND_THEN("The correct weighting is applied")
777 : {
778 4 : volume(3, 1) = 1;
779 4 : volume(2, 2) = 1;
780 4 : volume(2, 3) = 1;
781 :
782 4 : fast.apply(volume, sino);
783 4 : REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) * (sqrt3d - 1)
784 : + (2 / sqrt3d) * (3 - 8 * sqrt3d / 6))
785 : .epsilon(0.005));
786 :
787 4 : slow.apply(volume, sino);
788 4 : REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) * (sqrt3d - 1)
789 : + (2 / sqrt3d) * (3 - 8 * sqrt3d / 6))
790 : .epsilon(0.005));
791 :
792 4 : sino[0] = 1;
793 :
794 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
795 8 : slowExpected << 0, 0, 0, 0, 0, 0, 0, (4 - 2 * sqrt3d) * (sqrt3d - 1), 0, 0,
796 4 : (2 / sqrt3d) * (1 + halfd - 5 * sqrt3d / 6),
797 : (4 - 2 * sqrt3d) * (2 - sqrt3d)
798 4 : + (2 / sqrt3d) * (5 * sqrt3d / 6 - halfd),
799 4 : 0, 0, (2 / sqrt3d) * (1 + halfd - sqrt3d / 2),
800 4 : (2 / sqrt3d) * (sqrt3d / 2 - halfd);
801 :
802 4 : slow.applyAdjoint(sino, volume);
803 4 : REQUIRE_UNARY(isCwiseApprox(
804 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
805 :
806 4 : fast.applyAdjoint(sino, volume);
807 20 : for (real_t i = 0.5; i < volSize; i += 1) {
808 80 : for (real_t j = 0.5; j < volSize; j += 1) {
809 : const real_t angle =
810 64 : std::abs(std::atan((40 * sqrt3r - 2 + j) / (42 + sqrt3r - i))
811 : - pi_t / 3);
812 64 : const real_t len = 84 * std::tan(angle);
813 64 : if (len < 1) {
814 20 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
815 : Approx(1 - len).epsilon(0.01));
816 : } else {
817 44 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
818 : }
819 : }
820 : }
821 : }
822 : }
823 : }
824 :
825 20 : WHEN("Ray exits through the left border")
826 : {
827 : // In this case the ray enters through a border along the main ray direction, but
828 : // exits through a border not along the main direction Last pixel should be weighted
829 : // differently
830 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
831 : PrincipalPointOffset{0}, RotationOffset2D{-sqrt3r, 0});
832 :
833 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
834 8 : DataContainer<data_t> sino(sinoDescriptor);
835 4 : sino = 0;
836 :
837 8 : TestType fast(volumeDescriptor, sinoDescriptor);
838 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
839 :
840 8 : THEN("Ray intersects the correct pixels")
841 : {
842 4 : volume = 1;
843 4 : volume(0, 0) = 0;
844 4 : volume(1, 0) = 0;
845 4 : volume(0, 1) = 0;
846 4 : volume(1, 1) = 0;
847 4 : volume(0, 2) = 0;
848 :
849 4 : fast.apply(volume, sino);
850 8 : DataContainer<data_t> zero(sinoDescriptor);
851 4 : zero = 0;
852 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
853 :
854 4 : slow.apply(volume, sino);
855 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
856 :
857 8 : AND_THEN("The correct weighting is applied")
858 : {
859 4 : volume(1, 0) = 1;
860 4 : volume(0, 1) = 1;
861 :
862 4 : fast.apply(volume, sino);
863 4 : REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
864 : + (4 - 2 * sqrt3d) * (2 - sqrt3d))
865 : .epsilon(0.005));
866 :
867 4 : slow.apply(volume, sino);
868 4 : REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
869 : + (4 - 2 * sqrt3d) * (2 - sqrt3d))
870 : .epsilon(0.005));
871 :
872 4 : sino[0] = 1;
873 :
874 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
875 8 : slowExpected << 1 - 1 / sqrt3d, sqrt3d - 1, 0, 0,
876 4 : (5 * thirdd - 1 / sqrt3d) + (4 - 2 * sqrt3d) * (2 - sqrt3d),
877 4 : sqrt3d - 5 * thirdd, 0, 0, (sqrt3d - 1) * (4 - 2 * sqrt3d), 0, 0, 0, 0,
878 8 : 0, 0, 0;
879 :
880 4 : slow.applyAdjoint(sino, volume);
881 4 : REQUIRE_UNARY(isCwiseApprox(
882 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
883 :
884 4 : fast.applyAdjoint(sino, volume);
885 :
886 20 : for (real_t i = 0.5; i < volSize; i += 1) {
887 80 : for (real_t j = 0.5; j < volSize; j += 1) {
888 : const real_t angle =
889 64 : std::abs(std::atan((40 * sqrt3r - 2 + j) / (42 - sqrt3r - i))
890 : - pi_t / 3);
891 64 : const real_t len = 84 * std::tan(angle);
892 :
893 64 : if (len < 1) {
894 24 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
895 : Approx(1 - len).epsilon(0.002));
896 : } else {
897 40 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
898 : }
899 : }
900 : }
901 : }
902 : }
903 : }
904 :
905 20 : WHEN("Ray only intersects a single pixel")
906 : {
907 : // This is a special case that is handled separately in both forward and
908 : // backprojection
909 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
910 : PrincipalPointOffset{0}, RotationOffset2D{-2 - sqrt3r / 2, 0});
911 :
912 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
913 8 : DataContainer<data_t> sino(sinoDescriptor);
914 4 : sino = 0;
915 :
916 8 : TestType fast(volumeDescriptor, sinoDescriptor);
917 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
918 :
919 8 : THEN("Ray intersects the correct pixels")
920 : {
921 4 : volume = 1;
922 4 : volume(0, 0) = 0;
923 :
924 4 : fast.apply(volume, sino);
925 8 : DataContainer<data_t> zero(sinoDescriptor);
926 4 : zero = 0;
927 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
928 :
929 4 : slow.apply(volume, sino);
930 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
931 :
932 8 : AND_THEN("The correct weighting is applied")
933 : {
934 4 : volume(0, 0) = 1;
935 :
936 4 : fast.apply(volume, sino);
937 4 : REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
938 :
939 4 : slow.apply(volume, sino);
940 4 : REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
941 :
942 4 : sino[0] = 1;
943 :
944 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
945 4 : slowExpected << 1 / sqrt3d, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
946 :
947 4 : slow.applyAdjoint(sino, volume);
948 4 : REQUIRE_UNARY(isCwiseApprox(
949 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
950 :
951 4 : fast.applyAdjoint(sino, volume);
952 :
953 20 : for (real_t i = 0.5; i < volSize; i += 1) {
954 80 : for (real_t j = 0.5; j < volSize; j += 1) {
955 64 : const real_t angle = std::abs(
956 64 : std::atan((40 * sqrt3r - 2 + j) / (40 - sqrt3r / 2 - i))
957 : - pi_t / 3);
958 64 : const real_t len = 84 * std::tan(angle);
959 :
960 64 : if (len < 1) {
961 8 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
962 : Approx(1 - len).epsilon(0.002));
963 : } else {
964 56 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
965 : }
966 : }
967 : }
968 : }
969 : }
970 : }
971 : }
972 :
973 48 : GIVEN("An angle of -120 degrees")
974 : {
975 16 : const auto angle = Degree{-120};
976 :
977 20 : WHEN("Ray goes through center of volume")
978 : {
979 : // In this case the ray enters and exits the volume through the borders along the
980 : // main direction Weighting for all interpolated values should be the same
981 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData));
982 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
983 8 : DataContainer<data_t> sino(sinoDescriptor);
984 4 : sino = 0;
985 :
986 8 : TestType fast(volumeDescriptor, sinoDescriptor);
987 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
988 :
989 4 : real_t weight = 2 / sqrt3r;
990 8 : THEN("Ray intersects the correct pixels")
991 : {
992 4 : sino[0] = 1;
993 4 : slow.applyAdjoint(sino, volume);
994 :
995 4 : volume = 1;
996 4 : volume(0, 0) = 0;
997 4 : volume(0, 1) = 0;
998 4 : volume(1, 1) = 0;
999 4 : volume(1, 2) = 0;
1000 :
1001 4 : volume(2, 1) = 0;
1002 4 : volume(2, 2) = 0;
1003 4 : volume(3, 2) = 0;
1004 4 : volume(3, 3) = 0;
1005 :
1006 4 : fast.apply(volume, sino);
1007 8 : DataContainer<data_t> zero(sinoDescriptor);
1008 4 : zero = 0;
1009 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1010 :
1011 4 : slow.apply(volume, sino);
1012 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1013 :
1014 8 : AND_THEN("The correct weighting is applied")
1015 : {
1016 4 : volume(3, 3) = 1;
1017 4 : volume(0, 1) = 1;
1018 4 : volume(1, 1) = 1;
1019 4 : volume(2, 1) = 1;
1020 :
1021 4 : fast.apply(volume, sino);
1022 4 : REQUIRE_EQ(sino[0], Approx(2 * weight));
1023 :
1024 4 : slow.apply(volume, sino);
1025 4 : REQUIRE_EQ(sino[0], Approx(2 * weight));
1026 :
1027 4 : sino[0] = 1;
1028 :
1029 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
1030 :
1031 8 : slowExpected << (sqrt3d - 1) / 2, 0, 0, 0, (3 - sqrt3d) / 2,
1032 4 : (sqrt3d + 1) / (2 * sqrt3d), (sqrt3d - 1) / (2 * sqrt3d), 0, 0,
1033 4 : (sqrt3d - 1) / (2 * sqrt3d), (sqrt3d + 1) / (2 * sqrt3d),
1034 8 : (3 - sqrt3d) / 2, 0, 0, 0, (sqrt3d - 1) / 2;
1035 :
1036 4 : slowExpected *= weight;
1037 4 : slow.applyAdjoint(sino, volume);
1038 4 : REQUIRE_UNARY(isCwiseApprox(
1039 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
1040 :
1041 4 : fast.applyAdjoint(sino, volume);
1042 20 : for (real_t i = 0.5; i < volSize; i += 1) {
1043 80 : for (real_t j = 0.5; j < volSize; j += 1) {
1044 64 : const real_t angle = std::abs(
1045 64 : std::atan((sqrt3r * 40 + 2 - i) / (42 - j)) - pi_t / 3);
1046 64 : const real_t len = volSize * 21 * std::tan(angle);
1047 64 : if (len < 1) {
1048 32 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
1049 : Approx(1 - len).epsilon(0.005));
1050 : } else {
1051 32 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
1052 : }
1053 : }
1054 : }
1055 : }
1056 : }
1057 : }
1058 :
1059 20 : WHEN("Ray enters through the top border")
1060 : {
1061 : // In this case the ray exits through a border along the main ray direction, but
1062 : // enters through a border not along the main direction First pixel should be
1063 : // weighted differently
1064 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
1065 4 : PrincipalPointOffset{0}, RotationOffset2D{0, std::sqrt(3.f)});
1066 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1067 8 : DataContainer<data_t> sino(sinoDescriptor);
1068 4 : sino = 0;
1069 :
1070 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1071 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1072 :
1073 8 : THEN("Ray intersects the correct pixels")
1074 : {
1075 4 : volume = 1;
1076 4 : volume(0, 2) = 0;
1077 4 : volume(0, 3) = 0;
1078 4 : volume(1, 2) = 0;
1079 4 : volume(1, 3) = 0;
1080 4 : volume(2, 3) = 0;
1081 :
1082 4 : fast.apply(volume, sino);
1083 8 : DataContainer<data_t> zero(sinoDescriptor);
1084 4 : zero = 0;
1085 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1086 :
1087 4 : slow.apply(volume, sino);
1088 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1089 :
1090 8 : AND_THEN("The correct weighting is applied")
1091 : {
1092 4 : volume(2, 3) = 1;
1093 4 : volume(1, 2) = 1;
1094 4 : volume(1, 3) = 1;
1095 :
1096 4 : fast.apply(volume, sino);
1097 4 : REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) + (2 / sqrt3d)));
1098 :
1099 4 : slow.apply(volume, sino);
1100 4 : REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) + (2 / sqrt3d)));
1101 :
1102 4 : sino[0] = 1;
1103 :
1104 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
1105 :
1106 8 : slowExpected << 0, 0, 0, 0, 0, 0, 0, 0,
1107 4 : (2 / sqrt3d) * (1 + halfd - sqrt3d / 2),
1108 4 : (2 / sqrt3d) * (1 + halfd - 5 * sqrt3d / 6), 0, 0,
1109 4 : (2 / sqrt3d) * (sqrt3d / 2 - halfd),
1110 : (4 - 2 * sqrt3d) * (2 - sqrt3d)
1111 4 : + (2 / sqrt3d) * (5 * sqrt3d / 6 - halfd),
1112 8 : (4 - 2 * sqrt3d) * (sqrt3d - 1), 0;
1113 :
1114 4 : slow.applyAdjoint(sino, volume);
1115 4 : REQUIRE_UNARY(isCwiseApprox(
1116 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
1117 :
1118 4 : fast.applyAdjoint(sino, volume);
1119 20 : for (real_t i = 0.5; i < volSize; i += 1) {
1120 80 : for (real_t j = 0.5; j < volSize; j += 1) {
1121 : const real_t angle =
1122 64 : std::abs(std::atan((sqrt3r * 40 + 2 - i) / (42 + sqrt3r - j))
1123 : - pi_t / 3);
1124 64 : const real_t len = volSize * 21 * std::tan(angle);
1125 64 : if (len < 1) {
1126 20 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
1127 : Approx(1 - len).epsilon(0.01));
1128 : } else {
1129 44 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
1130 : }
1131 : }
1132 : }
1133 : }
1134 : }
1135 : }
1136 :
1137 20 : WHEN("Ray exits through the bottom border")
1138 : {
1139 : // In this case the ray enters through a border along the main ray direction, but
1140 : // exits through a border not along the main direction Last pixel should be weighted
1141 : // differently
1142 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
1143 4 : PrincipalPointOffset{0}, RotationOffset2D{0, -std::sqrt(3.f)});
1144 :
1145 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1146 8 : DataContainer<data_t> sino(sinoDescriptor);
1147 4 : sino = 0;
1148 :
1149 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1150 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1151 :
1152 8 : THEN("Ray intersects the correct pixels")
1153 : {
1154 4 : volume = 1;
1155 4 : volume(1, 0) = 0;
1156 4 : volume(2, 0) = 0;
1157 4 : volume(3, 0) = 0;
1158 4 : volume(2, 1) = 0;
1159 4 : volume(3, 1) = 0;
1160 :
1161 4 : fast.apply(volume, sino);
1162 8 : DataContainer<data_t> zero(sinoDescriptor);
1163 4 : zero = 0;
1164 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1165 :
1166 4 : slow.apply(volume, sino);
1167 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1168 :
1169 8 : AND_THEN("The correct weighting is applied")
1170 : {
1171 4 : volume(2, 0) = 1;
1172 4 : volume(3, 1) = 1;
1173 :
1174 4 : fast.apply(volume, sino);
1175 4 : REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
1176 : + (4 - 2 * sqrt3d) * (2 - sqrt3d))
1177 : .epsilon(0.005));
1178 :
1179 4 : slow.apply(volume, sino);
1180 4 : REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
1181 : + (4 - 2 * sqrt3d) * (2 - sqrt3d))
1182 : .epsilon(0.005));
1183 :
1184 4 : sino[0] = 1;
1185 :
1186 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
1187 :
1188 8 : slowExpected << 0, (sqrt3d - 1) * (4 - 2 * sqrt3d),
1189 4 : (5 * thirdd - 1 / sqrt3d) + (4 - 2 * sqrt3d) * (2 - sqrt3d),
1190 4 : 1 - 1 / sqrt3d, 0, 0, sqrt3d - 5 * thirdd, sqrt3d - 1, 0, 0, 0, 0, 0, 0,
1191 8 : 0, 0;
1192 :
1193 4 : slow.applyAdjoint(sino, volume);
1194 4 : REQUIRE_UNARY(isCwiseApprox(
1195 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
1196 :
1197 4 : fast.applyAdjoint(sino, volume);
1198 :
1199 20 : for (real_t i = 0.5; i < volSize; i += 1) {
1200 80 : for (real_t j = 0.5; j < volSize; j += 1) {
1201 : const real_t angle =
1202 64 : std::abs(std::atan((sqrt3r * 40 + 2 - i) / (42 - sqrt3r - j))
1203 : - pi_t / 3);
1204 64 : const real_t len = 84 * std::tan(angle);
1205 :
1206 64 : if (len < 1) {
1207 24 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
1208 : Approx(1 - len).epsilon(0.002));
1209 : } else {
1210 40 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
1211 : }
1212 : }
1213 : }
1214 : }
1215 : }
1216 : }
1217 :
1218 20 : WHEN("Ray only intersects a single pixel")
1219 : {
1220 : // This is a special case that is handled separately in both forward and
1221 : // backprojection
1222 4 : geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
1223 : PrincipalPointOffset{0},
1224 4 : RotationOffset2D{0, -2 - std::sqrt(3.f) / 2});
1225 :
1226 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1227 8 : DataContainer<data_t> sino(sinoDescriptor);
1228 4 : sino = 0;
1229 :
1230 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1231 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1232 :
1233 8 : THEN("Ray intersects the correct pixels")
1234 : {
1235 4 : volume = 1;
1236 4 : volume(3, 0) = 0;
1237 :
1238 4 : fast.apply(volume, sino);
1239 8 : DataContainer<data_t> zero(sinoDescriptor);
1240 4 : zero = 0;
1241 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1242 :
1243 4 : slow.apply(volume, sino);
1244 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1245 :
1246 8 : AND_THEN("The correct weighting is applied")
1247 : {
1248 4 : volume(3, 0) = 1;
1249 :
1250 4 : fast.apply(volume, sino);
1251 4 : REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
1252 :
1253 4 : slow.apply(volume, sino);
1254 4 : REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
1255 :
1256 4 : sino[0] = 1;
1257 :
1258 8 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
1259 4 : slowExpected << 0, 0, 0, 1 / sqrt3d, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
1260 :
1261 4 : slow.applyAdjoint(sino, volume);
1262 4 : REQUIRE_UNARY(isCwiseApprox(
1263 : volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
1264 :
1265 4 : fast.applyAdjoint(sino, volume);
1266 :
1267 20 : for (real_t i = 0.5; i < volSize; i += 1) {
1268 80 : for (real_t j = 0.5; j < volSize; j += 1) {
1269 64 : const real_t angle = std::abs(
1270 64 : std::atan((sqrt3r * 40 + 2 - i) / (40 - sqrt3r / 2 - j))
1271 : - pi_t / 3);
1272 64 : const real_t len = 84 * std::tan(angle);
1273 :
1274 64 : if (len < 1) {
1275 8 : REQUIRE_EQ(volume((index_t) i, (index_t) j),
1276 : Approx(1 - len).epsilon(0.002));
1277 : } else {
1278 56 : REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
1279 : }
1280 : }
1281 : }
1282 : }
1283 : }
1284 : }
1285 : }
1286 : }
1287 80 : }
1288 :
1289 34 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 2D setup with a multiple rays", TestType,
1290 : JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
1291 : {
1292 : // Turn logger of
1293 4 : Logger::setLevel(Logger::LogLevel::OFF);
1294 :
1295 : using data_t = decltype(return_data_t(std::declval<TestType>()));
1296 :
1297 8 : GIVEN("a 5x5 Volume and detector of size 1, with 4 angles")
1298 : {
1299 : // Domain setup
1300 4 : const index_t volSize = 5;
1301 :
1302 8 : IndexVector_t volumeDims(2);
1303 4 : volumeDims << volSize, volSize;
1304 :
1305 8 : VolumeDescriptor volumeDescriptor(volumeDims);
1306 8 : DataContainer<data_t> volume(volumeDescriptor);
1307 :
1308 : // range setup
1309 4 : const index_t detectorSize = 1;
1310 4 : const index_t numImgs = 4;
1311 :
1312 8 : IndexVector_t sinoDims(2);
1313 4 : sinoDims << detectorSize, numImgs;
1314 :
1315 : // Setup geometry
1316 4 : auto stc = SourceToCenterOfRotation{20 * volSize};
1317 4 : auto ctr = CenterOfRotationToDetector{volSize};
1318 12 : auto volData = VolumeData2D{Size2D{volumeDims}};
1319 12 : auto sinoData = SinogramData2D{Size2D{sinoDims}};
1320 :
1321 8 : std::vector<Geometry> geom;
1322 :
1323 8 : WHEN("Both x- and y-axis-aligned rays are present")
1324 : {
1325 4 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{volumeDims}},
1326 : SinogramData2D{Size2D{sinoDims}});
1327 4 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{volumeDims}},
1328 : SinogramData2D{Size2D{sinoDims}});
1329 4 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{volumeDims}},
1330 : SinogramData2D{Size2D{sinoDims}});
1331 4 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{volumeDims}},
1332 : SinogramData2D{Size2D{sinoDims}});
1333 :
1334 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1335 8 : DataContainer<data_t> sino(sinoDescriptor);
1336 4 : sino = 0;
1337 :
1338 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1339 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1340 :
1341 8 : THEN("Values are accumulated correctly along each ray's path")
1342 : {
1343 4 : volume = 0;
1344 :
1345 : // set only values along the rays' path to one to make sure interpolation is done
1346 : // correctly
1347 24 : for (index_t i = 0; i < volSize; i++) {
1348 20 : volume(i, volSize / 2) = 1;
1349 20 : volume(volSize / 2, i) = 1;
1350 : }
1351 :
1352 4 : slow.apply(volume, sino);
1353 20 : for (index_t i = 0; i < numImgs; i++)
1354 16 : REQUIRE_EQ(sino[i], Approx(5.0));
1355 :
1356 4 : fast.apply(volume, sino);
1357 20 : for (index_t i = 0; i < numImgs; i++)
1358 16 : REQUIRE_EQ(sino[i], Approx(5.0));
1359 :
1360 8 : AND_THEN("Both fast and slow backprojections yield the exact adjoint")
1361 : {
1362 4 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> cmp(volSize * volSize);
1363 :
1364 8 : cmp << 0, 0, 10, 0, 0, 0, 0, 10, 0, 0, 10, 10, 20, 10, 10, 0, 0, 10, 0, 0, 0, 0,
1365 8 : 10, 0, 0;
1366 :
1367 4 : slow.applyAdjoint(sino, volume);
1368 4 : REQUIRE_UNARY(
1369 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
1370 :
1371 4 : fast.applyAdjoint(sino, volume);
1372 4 : REQUIRE_UNARY(
1373 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
1374 : }
1375 : }
1376 : }
1377 : }
1378 4 : }
1379 74 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 3D setup with a single ray", TestType,
1380 : JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
1381 : {
1382 : // Turn logger of
1383 44 : Logger::setLevel(Logger::LogLevel::OFF);
1384 :
1385 : using data_t = decltype(return_data_t(std::declval<TestType>()));
1386 :
1387 80 : GIVEN("A 3x3x3 Volume")
1388 : {
1389 : // Domain setup
1390 36 : const index_t volSize = 3;
1391 :
1392 72 : IndexVector_t volumeDims(3);
1393 36 : volumeDims << volSize, volSize, volSize;
1394 :
1395 72 : VolumeDescriptor volumeDescriptor(volumeDims);
1396 72 : DataContainer<data_t> volume(volumeDescriptor);
1397 :
1398 : // range setup
1399 36 : const index_t detectorSize = 1;
1400 36 : const index_t numImgs = 1;
1401 :
1402 72 : IndexVector_t sinoDims(3);
1403 36 : sinoDims << detectorSize, detectorSize, numImgs;
1404 :
1405 : // Setup geometry
1406 36 : const auto stc = SourceToCenterOfRotation{20 * volSize};
1407 36 : const auto ctr = CenterOfRotationToDetector{volSize};
1408 108 : auto volData = VolumeData3D{Size3D{volumeDims}};
1409 108 : auto sinoData = SinogramData3D{Size3D{sinoDims}};
1410 :
1411 72 : std::vector<Geometry> geom;
1412 :
1413 44 : GIVEN("Single ray along z-axis")
1414 : {
1415 :
1416 8 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1417 : RotationAngles3D{Gamma{0}});
1418 :
1419 16 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1420 16 : DataContainer<data_t> sino(sinoDescriptor);
1421 8 : sino = 0;
1422 :
1423 16 : TestType fast(volumeDescriptor, sinoDescriptor);
1424 16 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1425 :
1426 12 : WHEN("Sinogram conatainer is not zero initialized and we project through an empty "
1427 : "volume")
1428 : {
1429 4 : volume = 0;
1430 4 : sino = 1;
1431 :
1432 8 : THEN("Result is zero")
1433 : {
1434 4 : fast.apply(volume, sino);
1435 4 : DataContainer<data_t> zero(sinoDescriptor);
1436 4 : zero = 0;
1437 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1438 :
1439 4 : sino = 1;
1440 4 : slow.apply(volume, sino);
1441 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
1442 : }
1443 : }
1444 :
1445 12 : WHEN("Volume container is not zero initialized and we backproject from an empty "
1446 : "sinogram")
1447 : {
1448 4 : sino = 0;
1449 4 : volume = 1;
1450 :
1451 8 : THEN("Result is zero")
1452 : {
1453 4 : fast.applyAdjoint(sino, volume);
1454 4 : DataContainer<data_t> zero(volumeDescriptor);
1455 4 : zero = 0;
1456 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
1457 :
1458 4 : volume = 1;
1459 4 : slow.applyAdjoint(sino, volume);
1460 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
1461 : }
1462 : }
1463 : }
1464 :
1465 36 : const index_t numCases = 6;
1466 36 : std::array<real_t, numCases> beta = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
1467 36 : std::array<real_t, numCases> gamma = {0.0, pi_t, pi_t / 2,
1468 : 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
1469 72 : std::array<std::string, numCases> al = {"z", "-z", "x", "-x", "y", "-y"};
1470 :
1471 504 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProjections[numCases];
1472 252 : for (auto& backPr : backProjections)
1473 216 : backPr.resize(volSize * volSize * volSize);
1474 :
1475 : // clang-format off
1476 72 : backProjections[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0,
1477 36 : 0, 0, 0, 0, 1, 0, 0, 0, 0,
1478 72 : 0, 0, 0, 0, 1, 0, 0, 0, 0;
1479 :
1480 72 : backProjections[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0,
1481 36 : 0, 0, 0, 1, 1, 1, 0, 0, 0,
1482 72 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1483 :
1484 72 : backProjections[2] << 0, 0, 0, 0, 0, 0, 0, 0, 0,
1485 36 : 0, 1, 0, 0, 1, 0, 0, 1, 0,
1486 72 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1487 : // clang-format on
1488 :
1489 252 : for (std::size_t i = 0; i < numCases; i++) {
1490 220 : WHEN("An axis-aligned ray passes through the center of a pixel")
1491 : {
1492 8 : INFO("A ", al[i], "-axis-aligned ray passes through the center of a pixel");
1493 4 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1494 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
1495 :
1496 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1497 8 : DataContainer<data_t> sino(sinoDescriptor);
1498 4 : sino = 0;
1499 :
1500 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1501 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1502 :
1503 8 : THEN("The result of projecting through a voxel is exactly the voxel value")
1504 : {
1505 16 : for (index_t j = 0; j < volSize; j++) {
1506 12 : volume = 0;
1507 12 : if (i / 2 == 0)
1508 12 : volume(volSize / 2, volSize / 2, j) = 1;
1509 0 : else if (i / 2 == 1)
1510 0 : volume(j, volSize / 2, volSize / 2) = 1;
1511 0 : else if (i / 2 == 2)
1512 0 : volume(volSize / 2, j, volSize / 2) = 1;
1513 :
1514 12 : fast.apply(volume, sino);
1515 : // Using doubles significantly increases interpolation accuracy
1516 : // For example: when using floats, points very near the border (less than
1517 : // ~1/500th of a pixel away from the border) are rounded to actually lie on
1518 : // the border. This then yields more accurate results when using floats in
1519 : // some of the axis-aligned test cases, despite the lower interpolation
1520 : // accuracy.
1521 : // => different requirements for floats and doubles, looser requirements
1522 : // for doubles
1523 12 : REQUIRE_EQ(sino[0], Approx(1.0));
1524 :
1525 12 : slow.apply(volume, sino);
1526 12 : REQUIRE_EQ(sino[0], Approx(1.0));
1527 : }
1528 :
1529 8 : AND_THEN("The backprojection sets the values of all hit voxels to the detector "
1530 : "value")
1531 : {
1532 4 : fast.applyAdjoint(sino, volume);
1533 4 : REQUIRE_UNARY(
1534 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor,
1535 : backProjections[i / 2])));
1536 :
1537 4 : slow.applyAdjoint(sino, volume);
1538 4 : REQUIRE_UNARY(
1539 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor,
1540 : backProjections[i / 2])));
1541 : }
1542 : }
1543 : }
1544 : }
1545 :
1546 36 : std::array<real_t, numCases> offsetx = {-0.25, -0.25, 0.0, 0.0, 0.0, 0.0};
1547 36 : std::array<real_t, numCases> offsety = {0.0, 0.0, -0.25, -0.25, 0.0, 0.0};
1548 36 : std::array<real_t, numCases> offsetz = {0.0, 0.0, 0.0, 0.0, -0.25, -0.25};
1549 :
1550 : // clang-format off
1551 72 : backProjections[0] << 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0,
1552 36 : 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0,
1553 72 : 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0;
1554 :
1555 72 : backProjections[1] << 0 , 0 , 0 , 0 , 0 , 0 , 0, 0, 0,
1556 36 : 0.25, 0.25, 0.25, 0.75, 0.75, 0.75, 0, 0, 0,
1557 72 : 0 , 0 , 0 , 0 , 0 , 0 , 0, 0, 0;
1558 :
1559 72 : backProjections[2] << 0, 0.25, 0, 0, 0.25, 0, 0, 0.25, 0,
1560 36 : 0, 0.75, 0, 0, 0.75, 0, 0, 0.75, 0,
1561 72 : 0, 0 , 0, 0, 0 , 0, 0, 0 , 0;
1562 : // clang-format on
1563 :
1564 252 : for (std::size_t i = 0; i < numCases; i++) {
1565 220 : WHEN("An axis-aligned ray not passing through the center of a pixel")
1566 : {
1567 8 : INFO("A ", al[i], "-axis-aligned ray not passing through the center of a pixel");
1568 : // x-ray source must be very far from the volume center to make testing of the fast
1569 : // backprojection simpler
1570 8 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, std::move(volData),
1571 4 : std::move(sinoData),
1572 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}},
1573 : PrincipalPointOffset2D{0, 0},
1574 : RotationOffset3D{offsetx[i], offsety[i], offsetz[i]});
1575 :
1576 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1577 8 : DataContainer<data_t> sino(sinoDescriptor);
1578 4 : sino = 0;
1579 :
1580 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1581 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1582 :
1583 8 : THEN("The result of projecting through a voxel is the interpolated value between "
1584 : "the four voxels nearest to the ray")
1585 : {
1586 16 : for (index_t j = 0; j < volSize; j++) {
1587 12 : volume = 0;
1588 12 : if (i / 2 == 0)
1589 12 : volume(volSize / 2, volSize / 2, j) = 1;
1590 0 : else if (i / 2 == 1)
1591 0 : volume(j, volSize / 2, volSize / 2) = 1;
1592 0 : else if (i / 2 == 2)
1593 0 : volume(volSize / 2, j, volSize / 2) = 1;
1594 :
1595 12 : fast.apply(volume, sino);
1596 :
1597 : // Using doubles significantly increases interpolation accuracy
1598 : // For example: when using floats, points very near the border (less than
1599 : // ~1/500th of a pixel away from the border) are rounded to actually lie on
1600 : // the border. This then yields more accurate results when using floats in
1601 : // some of the axis-aligned test cases, despite the lower interpolation
1602 : // accuracy.
1603 : //
1604 : // => different requirements for floats and doubles, looser requirements
1605 : // for doubles
1606 :
1607 12 : REQUIRE_EQ(sino[0], Approx(0.75));
1608 :
1609 12 : slow.apply(volume, sino);
1610 12 : REQUIRE_EQ(sino[0], Approx(0.75));
1611 : }
1612 :
1613 8 : AND_THEN("The slow backprojection yields the exact adjoint, the fast "
1614 : "backprojection is also exact for a very distant x-ray source")
1615 : {
1616 4 : sino[0] = 1;
1617 4 : fast.applyAdjoint(sino, volume);
1618 : // Using doubles significantly increases interpolation accuracy
1619 : // For example: when using floats, points very near the border (less than
1620 : // ~1/500th of a pixel away from the border) are rounded to actually lie on
1621 : // the border. This then yields more accurate results when using floats in
1622 : // some of the axis-aligned test cases, despite the lower interpolation
1623 : // accuracy.
1624 : //
1625 : // => different requirements for floats and doubles, looser requirements
1626 : // for doubles
1627 : if constexpr (std::is_same_v<data_t, float>)
1628 2 : REQUIRE_UNARY(isCwiseApprox(
1629 : volume,
1630 : DataContainer<data_t>(volumeDescriptor, backProjections[i / 2])));
1631 : else
1632 2 : REQUIRE_UNARY(isApprox(
1633 : volume,
1634 : DataContainer<data_t>(volumeDescriptor, backProjections[i / 2]),
1635 : static_cast<real_t>(0.005)));
1636 :
1637 4 : slow.applyAdjoint(sino, volume);
1638 4 : REQUIRE_UNARY(
1639 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor,
1640 : backProjections[i / 2])));
1641 : }
1642 : }
1643 : }
1644 : }
1645 :
1646 36 : offsetx[0] = -volSize / 2.0;
1647 36 : offsetx[1] = volSize / 2.0;
1648 36 : offsetx[2] = 0.0;
1649 36 : offsetx[3] = 0.0;
1650 36 : offsetx[4] = volSize / 2.0;
1651 36 : offsetx[5] = -volSize / 2.0;
1652 :
1653 36 : offsety[0] = 0.0;
1654 36 : offsety[1] = 0.0;
1655 36 : offsety[2] = -volSize / 2.0;
1656 36 : offsety[3] = volSize / 2.0;
1657 36 : offsety[4] = volSize / 2.0;
1658 36 : offsety[5] = -volSize / 2.0;
1659 :
1660 : // clang-format off
1661 72 : backProjections[0] << 0, 0, 0, 1, 0, 0, 0, 0, 0,
1662 36 : 0, 0, 0, 1, 0, 0, 0, 0, 0,
1663 72 : 0, 0, 0, 1, 0, 0, 0, 0, 0;
1664 :
1665 72 : backProjections[1] << 0, 0, 0, 0, 0, 1, 0, 0, 0,
1666 36 : 0, 0, 0, 0, 0, 1, 0, 0, 0,
1667 72 : 0, 0, 0, 0, 0, 1, 0, 0, 0;
1668 :
1669 72 : backProjections[2] << 0, 1, 0, 0, 0, 0, 0, 0, 0,
1670 36 : 0, 1, 0, 0, 0, 0, 0, 0, 0,
1671 72 : 0, 1, 0, 0, 0, 0, 0, 0, 0;
1672 :
1673 72 : backProjections[3] << 0, 0, 0, 0, 0, 0, 0, 1, 0,
1674 36 : 0, 0, 0, 0, 0, 0, 0, 1, 0,
1675 72 : 0, 0, 0, 0, 0, 0, 0, 1, 0;
1676 :
1677 72 : backProjections[4] << 0, 0, 0, 0, 0, 0, 0, 0, 1,
1678 36 : 0, 0, 0, 0, 0, 0, 0, 0, 1,
1679 72 : 0, 0, 0, 0, 0, 0, 0, 0, 1;
1680 :
1681 72 : backProjections[5] << 1, 0, 0, 0, 0, 0, 0, 0, 0,
1682 36 : 1, 0, 0, 0, 0, 0, 0, 0, 0,
1683 72 : 1, 0, 0, 0, 0, 0, 0, 0, 0;
1684 : // clang-format on
1685 :
1686 36 : al[0] = "left border";
1687 36 : al[1] = "right border";
1688 36 : al[2] = "top border";
1689 36 : al[3] = "bottom border";
1690 36 : al[4] = "top right edge";
1691 36 : al[5] = "bottom left edge";
1692 :
1693 252 : for (std::size_t i = 0; i < numCases; i++) {
1694 220 : WHEN("A z-axis-aligned ray runs along the corners and edges of the volume")
1695 : {
1696 8 : INFO("A z-axis-aligned ray runs along the ", al[i], " of the volume");
1697 : // x-ray source must be very far from the volume center to make testing of the fast
1698 : // backprojection simpler
1699 8 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, std::move(volData),
1700 4 : std::move(sinoData), RotationAngles3D{Gamma{0}},
1701 : PrincipalPointOffset2D{0, 0},
1702 : RotationOffset3D{offsetx[i], offsety[i], 0});
1703 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1704 8 : DataContainer<data_t> sino(sinoDescriptor);
1705 4 : sino = 0;
1706 :
1707 8 : TestType fast(volumeDescriptor, sinoDescriptor);
1708 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
1709 :
1710 8 : THEN("The result of projecting through a voxel is exactly the voxel's value (we "
1711 : "mirror values at the border for the purpose of interpolation)")
1712 : {
1713 16 : for (index_t j = 0; j < volSize; j++) {
1714 12 : volume = 0;
1715 12 : switch (i) {
1716 12 : case 0:
1717 12 : volume(0, volSize / 2, j) = 1;
1718 12 : break;
1719 0 : case 1:
1720 0 : volume(volSize - 1, volSize / 2, j) = 1;
1721 0 : break;
1722 0 : case 2:
1723 0 : volume(volSize / 2, 0, j) = 1;
1724 0 : break;
1725 0 : case 3:
1726 0 : volume(volSize / 2, volSize - 1, j) = 1;
1727 0 : break;
1728 0 : case 4:
1729 0 : volume(volSize - 1, volSize - 1, j) = 1;
1730 0 : break;
1731 0 : case 5:
1732 0 : volume(0, 0, j) = 1;
1733 0 : break;
1734 0 : default:
1735 0 : break;
1736 : }
1737 :
1738 12 : fast.apply(volume, sino);
1739 12 : REQUIRE_EQ(sino[0], Approx(1));
1740 :
1741 12 : slow.apply(volume, sino);
1742 12 : REQUIRE_EQ(sino[0], Approx(1));
1743 : }
1744 :
1745 8 : AND_THEN("The slow backprojection yields the exact adjoint, the fast "
1746 : "backprojection is also exact for a very distant x-ray source")
1747 : {
1748 4 : sino[0] = 1;
1749 4 : fast.applyAdjoint(sino, volume);
1750 :
1751 : // Using doubles significantly increases interpolation accuracy
1752 : // For example: when using floats, points very near the border (less than
1753 : // ~1/500th of a pixel away from the border) are rounded to actually lie on
1754 : // the border. This then yields more accurate results when using floats in
1755 : // some of the axis-aligned test cases, despite the lower interpolation
1756 : // accuracy.
1757 : //
1758 : // => different requirements for floats and doubles, looser requirements
1759 : // for doubles
1760 : if constexpr (std::is_same_v<data_t, float>)
1761 2 : REQUIRE_UNARY(isCwiseApprox(
1762 : volume, DataContainer<data_t>(
1763 : volumeDescriptor,
1764 : (backProjections[i] / (i > 3 ? 4 : 2)).eval())));
1765 : else
1766 2 : REQUIRE_UNARY(
1767 : isApprox(volume,
1768 : DataContainer<data_t>(
1769 : volumeDescriptor,
1770 : (backProjections[i] / (i > 3 ? 4 : 2)).eval()),
1771 : static_cast<real_t>(0.005)));
1772 :
1773 4 : slow.applyAdjoint(sino, volume);
1774 4 : REQUIRE_UNARY(isCwiseApprox(
1775 : volume, DataContainer<data_t>(volumeDescriptor, backProjections[i])));
1776 : }
1777 : }
1778 : }
1779 : }
1780 :
1781 36 : const data_t sqrt3d = std::sqrt(static_cast<data_t>(3));
1782 36 : const data_t thirdd = static_cast<data_t>(1.0 / 3);
1783 :
1784 72 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProj(volSize * volSize * volSize);
1785 :
1786 52 : GIVEN("A 30 degree rotate around z axis")
1787 : {
1788 16 : auto gamma = Gamma{Degree{30}};
1789 :
1790 20 : WHEN("A ray goes through the center of the volume")
1791 : {
1792 : // In this case the ray enters and exits the volume along the main direction
1793 4 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1794 : RotationAngles3D{gamma});
1795 :
1796 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1797 8 : DataContainer<data_t> sino(sinoDescriptor);
1798 :
1799 8 : TestType op(volumeDescriptor, sinoDescriptor, false);
1800 :
1801 8 : THEN("The ray intersects the correct voxels")
1802 : {
1803 4 : volume = 1;
1804 4 : volume(1, 1, 1) = 0;
1805 4 : volume(2, 1, 0) = 0;
1806 4 : volume(1, 1, 0) = 0;
1807 4 : volume(0, 1, 2) = 0;
1808 4 : volume(1, 1, 2) = 0;
1809 :
1810 4 : op.apply(volume, sino);
1811 4 : REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
1812 :
1813 8 : AND_THEN("The correct weighting is applied")
1814 : {
1815 4 : volume(1, 1, 1) = 1;
1816 4 : volume(2, 1, 0) = 3;
1817 4 : volume(1, 1, 2) = 2;
1818 :
1819 4 : op.apply(volume, sino);
1820 4 : REQUIRE_EQ(sino[0], Approx(6 / sqrt3d + 2.0 / 3).epsilon(0.001));
1821 :
1822 4 : sino[0] = 1;
1823 : // clang-format off
1824 8 : backProj << 0, 0, 0, 0 , 2 / sqrt3d - 2 * thirdd, 2 * thirdd, 0, 0, 0,
1825 4 : 0, 0, 0, 0 , 2 / sqrt3d , 0 , 0, 0, 0,
1826 8 : 0, 0, 0, 2 * thirdd, 2 / sqrt3d - 2 * thirdd, 0 , 0, 0, 0;
1827 : // clang-format on
1828 :
1829 4 : op.applyAdjoint(sino, volume);
1830 4 : REQUIRE_UNARY(isCwiseApprox(
1831 : volume, DataContainer<data_t>(volumeDescriptor, backProj)));
1832 : }
1833 : }
1834 : }
1835 :
1836 20 : WHEN("A ray enters through the right border")
1837 : {
1838 : // In this case the ray enters through a border orthogonal to a non-main direction
1839 4 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1840 : RotationAngles3D{gamma}, PrincipalPointOffset2D{0, 0},
1841 : RotationOffset3D{1, 0, 0});
1842 :
1843 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1844 8 : DataContainer<data_t> sino(sinoDescriptor);
1845 :
1846 8 : TestType op(volumeDescriptor, sinoDescriptor, false);
1847 :
1848 8 : THEN("The ray intersects the correct voxels")
1849 : {
1850 4 : volume = 1;
1851 4 : volume(2, 1, 1) = 0;
1852 4 : volume(2, 1, 0) = 0;
1853 4 : volume(2, 1, 2) = 0;
1854 4 : volume(1, 1, 2) = 0;
1855 :
1856 4 : op.apply(volume, sino);
1857 4 : REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
1858 :
1859 8 : AND_THEN("The correct weighting is applied")
1860 : {
1861 4 : volume(2, 1, 0) = 4;
1862 4 : volume(1, 1, 2) = 3;
1863 4 : volume(2, 1, 1) = 1;
1864 :
1865 4 : op.apply(volume, sino);
1866 4 : REQUIRE_EQ(sino[0], Approx((sqrt3d + 1) * (1 - 1 / sqrt3d) + 3 - sqrt3d / 2
1867 : + 2 / sqrt3d)
1868 : .epsilon(0.001));
1869 :
1870 4 : sino[0] = 1;
1871 : // clang-format off
1872 8 : backProj << 0, 0, 0, 0, 0, ((sqrt3d + 1) / 4) * (1 - 1 / sqrt3d), 0, 0, 0,
1873 4 : 0, 0, 0, 0, 0, 2 / sqrt3d + 1 - sqrt3d / 2, 0, 0, 0, 0, 0, 0, 0,
1874 8 : 2 * thirdd, 2 / sqrt3d - 2 * thirdd, 0, 0, 0;
1875 : // clang-format on
1876 :
1877 4 : op.applyAdjoint(sino, volume);
1878 4 : REQUIRE_UNARY(isCwiseApprox(
1879 : volume, DataContainer<data_t>(volumeDescriptor, backProj)));
1880 : }
1881 : }
1882 : }
1883 :
1884 20 : WHEN("A ray exits through the left border")
1885 : {
1886 : // In this case the ray exit through a border orthogonal to a non-main direction
1887 4 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1888 : RotationAngles3D{gamma}, PrincipalPointOffset2D{0, 0},
1889 : RotationOffset3D{-1, 0, 0});
1890 :
1891 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1892 8 : DataContainer<data_t> sino(sinoDescriptor);
1893 :
1894 8 : TestType op(volumeDescriptor, sinoDescriptor, false);
1895 :
1896 8 : THEN("The ray intersects the correct voxels")
1897 : {
1898 4 : volume = 1;
1899 4 : volume(0, 1, 0) = 0;
1900 4 : volume(1, 1, 0) = 0;
1901 4 : volume(0, 1, 1) = 0;
1902 4 : volume(0, 1, 2) = 0;
1903 :
1904 4 : op.apply(volume, sino);
1905 4 : REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
1906 :
1907 8 : AND_THEN("The correct weighting is applied")
1908 : {
1909 4 : volume(0, 1, 2) = 4;
1910 4 : volume(1, 1, 0) = 3;
1911 4 : volume(0, 1, 1) = 1;
1912 :
1913 4 : op.apply(volume, sino);
1914 4 : REQUIRE_EQ(sino[0], Approx((sqrt3d + 1) * (1 - 1 / sqrt3d) + 3 - sqrt3d / 2
1915 : + 2 / sqrt3d)
1916 : .epsilon(0.001));
1917 :
1918 4 : sino[0] = 1;
1919 : // clang-format off
1920 8 : backProj << 0, 0, 0, 2 / sqrt3d - 2 * thirdd, 2 * thirdd, 0, 0, 0, 0, 0, 0,
1921 4 : 0, 2 / sqrt3d + 1 - sqrt3d / 2, 0, 0, 0, 0, 0, 0, 0, 0,
1922 8 : ((sqrt3d + 1) / 4) * (1 - 1 / sqrt3d), 0, 0, 0, 0, 0;
1923 : // clang-format on
1924 :
1925 4 : op.applyAdjoint(sino, volume);
1926 4 : REQUIRE_UNARY(isCwiseApprox(
1927 : volume, DataContainer<data_t>(volumeDescriptor, backProj)));
1928 : }
1929 : }
1930 : }
1931 :
1932 20 : WHEN("A ray with an angle of 30 degrees only intersects a single voxel")
1933 : {
1934 : // special case - no interior voxels, entry and exit voxels are the same
1935 4 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
1936 : RotationAngles3D{gamma}, PrincipalPointOffset2D{0, 0},
1937 : RotationOffset3D{-2, 0, 0});
1938 :
1939 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
1940 8 : DataContainer<data_t> sino(sinoDescriptor);
1941 :
1942 8 : TestType op(volumeDescriptor, sinoDescriptor, false);
1943 :
1944 8 : THEN("The ray intersects the correct voxels")
1945 : {
1946 4 : volume = 1;
1947 4 : volume(0, 1, 0) = 0;
1948 :
1949 4 : op.apply(volume, sino);
1950 4 : REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
1951 :
1952 8 : AND_THEN("The correct weighting is applied")
1953 : {
1954 4 : volume(0, 1, 0) = 1;
1955 :
1956 4 : op.apply(volume, sino);
1957 4 : REQUIRE_EQ(sino[0], Approx(sqrt3d - 1));
1958 :
1959 4 : sino[0] = 1;
1960 : // clang-format off
1961 8 : backProj << 0, 0, 0, sqrt3d - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1962 8 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1963 : // clang-format off
1964 :
1965 4 : op.applyAdjoint(sino, volume);
1966 4 : REQUIRE_UNARY(
1967 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, backProj)));
1968 : }
1969 : }
1970 : }
1971 : }
1972 : }
1973 :
1974 52 : GIVEN("A 5x5x5 Volume")
1975 : {
1976 : // Domain setup
1977 8 : const index_t volSize = 5;
1978 :
1979 16 : IndexVector_t volumeDims(3);
1980 8 : volumeDims << volSize, volSize, volSize;
1981 :
1982 16 : VolumeDescriptor volumeDescriptor(volumeDims);
1983 16 : DataContainer<data_t> volume(volumeDescriptor);
1984 :
1985 : // range setup
1986 8 : const index_t detectorSize = 1;
1987 8 : const index_t numImgs = 1;
1988 :
1989 16 : IndexVector_t sinoDims(3);
1990 8 : sinoDims << detectorSize, detectorSize, numImgs;
1991 :
1992 : // Setup geometry
1993 8 : const auto stc = SourceToCenterOfRotation{20 * volSize};
1994 8 : const auto ctr = CenterOfRotationToDetector{volSize};
1995 24 : auto volData = VolumeData3D{Size3D{volumeDims}};
1996 24 : auto sinoData = SinogramData3D{Size3D{sinoDims}};
1997 :
1998 16 : std::vector<Geometry> geom;
1999 :
2000 16 : GIVEN("Rays not intersecting the volume")
2001 : {
2002 8 : constexpr index_t numCases = 9;
2003 8 : std::array<real_t, numCases> alpha = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2004 8 : std::array<real_t, numCases> beta = {0.0, 0.0, 0.0, 0.0, 0.0,
2005 : 0.0, pi_t / 2, pi_t / 2, pi_t / 2};
2006 8 : std::array<real_t, numCases> gamma = {0.0, 0.0, 0.0, pi_t / 2, pi_t / 2,
2007 : pi_t / 2, pi_t / 2, pi_t / 2, pi_t / 2};
2008 8 : std::array<real_t, numCases> offsetx = {-volSize, 0.0, -volSize, 0.0, 0.0,
2009 : 0.0, -volSize, 0.0, -volSize};
2010 8 : std::array<real_t, numCases> offsety = {0.0, -volSize, -volSize, -volSize, 0.0,
2011 : -volSize, 0.0, 0.0, 0.0};
2012 8 : std::array<real_t, numCases> offsetz = {0.0, 0.0, 0.0, 0.0, -volSize,
2013 : -volSize, 0.0, -volSize, -volSize};
2014 16 : std::array<std::string, numCases> neg = {"x", "y", "x and y", "y", "z",
2015 : "y and z", "x", "z", "x and z"};
2016 16 : std::array<std::string, numCases> ali = {"z", "z", "z", "x", "x", "x", "y", "y", "y"};
2017 :
2018 80 : for (std::size_t i = 0; i < numCases; i++) {
2019 76 : WHEN("Tracing along a fixed axis-aligned ray with negative coordinate of origin")
2020 : {
2021 8 : INFO("Tracing along a ", ali[i], "-axis-aligned ray with negative ", neg[i],
2022 : "-coodinate of origin");
2023 8 : geom.emplace_back(
2024 4 : stc, ctr, std::move(volData), std::move(sinoData),
2025 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}, Alpha{alpha[i]}},
2026 : PrincipalPointOffset2D{0, 0},
2027 4 : RotationOffset3D{-offsetx[i], -offsety[i], -offsetz[i]});
2028 :
2029 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
2030 8 : DataContainer<data_t> sino(sinoDescriptor);
2031 4 : sino = 0;
2032 :
2033 8 : TestType fast(volumeDescriptor, sinoDescriptor);
2034 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
2035 :
2036 8 : THEN("Result of forward projection is zero")
2037 : {
2038 4 : fast.apply(volume, sino);
2039 8 : DataContainer<data_t> zero(sinoDescriptor);
2040 4 : zero = 0;
2041 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
2042 :
2043 4 : slow.apply(volume, sino);
2044 4 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
2045 :
2046 8 : AND_THEN("Result of backprojection is zero")
2047 : {
2048 4 : fast.applyAdjoint(sino, volume);
2049 4 : DataContainer<data_t> zero(volumeDescriptor);
2050 4 : zero = 0;
2051 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
2052 :
2053 4 : slow.applyAdjoint(sino, volume);
2054 4 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
2055 : }
2056 : }
2057 : }
2058 : }
2059 : }
2060 : }
2061 44 : }
2062 34 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 3D setup with a multiple rays",TestType, JosephsMethodCUDA<float>,
2063 : JosephsMethodCUDA<double>)
2064 : {
2065 : // Turn logger of
2066 4 : Logger::setLevel(Logger::LogLevel::OFF);
2067 :
2068 : using data_t = decltype(return_data_t(std::declval<TestType>()));
2069 :
2070 8 : GIVEN("A 3x3x3 volume with 6 projection angles")
2071 : {
2072 : // Domain setup
2073 4 : const index_t volSize = 3;
2074 :
2075 8 : IndexVector_t volumeDims(3);
2076 4 : volumeDims << volSize, volSize, volSize;
2077 :
2078 8 : VolumeDescriptor volumeDescriptor(volumeDims);
2079 8 : DataContainer<data_t> volume(volumeDescriptor);
2080 :
2081 : // range setup
2082 4 : const index_t detectorSize = 1;
2083 4 : const index_t numImgs = 6;
2084 :
2085 8 : IndexVector_t sinoDims(3);
2086 4 : sinoDims << detectorSize, detectorSize, numImgs;
2087 :
2088 : // Setup geometry
2089 4 : const auto stc = SourceToCenterOfRotation{20 * volSize};
2090 4 : const auto ctr = CenterOfRotationToDetector{volSize};
2091 12 : auto volData = VolumeData3D{Size3D{volumeDims}};
2092 12 : auto sinoData = SinogramData3D{Size3D{sinoDims}};
2093 :
2094 8 : std::vector<Geometry> geom;
2095 :
2096 8 : WHEN("x-, y and z-axis-aligned rays are present")
2097 : {
2098 4 : real_t beta[numImgs] = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
2099 4 : real_t gamma[numImgs] = {0.0, pi_t, pi_t / 2, 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
2100 :
2101 28 : for (index_t i = 0; i < numImgs; i++)
2102 24 : geom.emplace_back(stc, ctr, VolumeData3D{Size3D{volumeDims}},
2103 : SinogramData3D{Size3D{sinoDims}},
2104 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
2105 :
2106 8 : PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
2107 8 : DataContainer<data_t> sino(sinoDescriptor);
2108 :
2109 8 : TestType slow(volumeDescriptor, sinoDescriptor, false);
2110 8 : TestType fast(volumeDescriptor, sinoDescriptor);
2111 :
2112 8 : THEN("Values are accumulated correctly along each ray's path")
2113 : {
2114 4 : volume = 0;
2115 :
2116 : // set only values along the rays' path to one to make sure interpolation is done
2117 : // correctly
2118 16 : for (index_t i = 0; i < volSize; i++) {
2119 12 : volume(i, volSize / 2, volSize / 2) = 1;
2120 12 : volume(volSize / 2, i, volSize / 2) = 1;
2121 12 : volume(volSize / 2, volSize / 2, i) = 1;
2122 : }
2123 :
2124 4 : slow.apply(volume, sino);
2125 28 : for (index_t i = 0; i < numImgs; i++)
2126 24 : REQUIRE_EQ(sino[i], Approx(3.0));
2127 :
2128 4 : fast.apply(volume, sino);
2129 28 : for (index_t i = 0; i < numImgs; i++)
2130 24 : REQUIRE_EQ(sino[i], Approx(3.0));
2131 :
2132 8 : AND_THEN("Both fast and slow backprojections yield the exact adjoint")
2133 : {
2134 4 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> cmp(volSize * volSize * volSize);
2135 :
2136 : // clang-format off
2137 8 : cmp << 0, 0, 0, 0, 6, 0, 0, 0, 0,
2138 4 : 0, 6, 0, 6, 18, 6, 0, 6, 0,
2139 8 : 0, 0, 0, 0, 6, 0, 0, 0, 0;
2140 : // clang-format on
2141 :
2142 4 : slow.applyAdjoint(sino, volume);
2143 4 : REQUIRE_UNARY(
2144 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
2145 :
2146 4 : fast.applyAdjoint(sino, volume);
2147 4 : REQUIRE_UNARY(
2148 : isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
2149 : }
2150 : }
2151 : }
2152 : }
2153 4 : }
|