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