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