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