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