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