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 : TEST_CASE("BinaryMethod: Testing with only one ray")
28 6 : {
29 : // Turn logger of
30 6 : Logger::setLevel(Logger::LogLevel::OFF);
31 :
32 6 : IndexVector_t sizeDomain(2);
33 6 : sizeDomain << 5, 5;
34 :
35 6 : IndexVector_t sizeRange(2);
36 6 : sizeRange << 1, 1;
37 :
38 6 : auto domain = VolumeDescriptor(sizeDomain);
39 : // auto range = VolumeDescriptor(sizeRange);
40 :
41 6 : auto stc = SourceToCenterOfRotation{100};
42 6 : auto ctr = CenterOfRotationToDetector{5};
43 6 : auto volData = VolumeData2D{Size2D{sizeDomain}};
44 6 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
45 :
46 6 : Eigen::IOFormat CommaInitFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ", ", "", "",
47 6 : " << ", ";");
48 :
49 6 : GIVEN("A BinaryMethod for 2D and a domain data with all ones")
50 6 : {
51 6 : std::vector<Geometry> geom;
52 :
53 6 : auto dataDomain = DataContainer(domain);
54 6 : dataDomain = 1;
55 :
56 6 : WHEN("We have a single ray with 0 degrees")
57 6 : {
58 1 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData));
59 :
60 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
61 1 : auto op = BinaryMethod(domain, range);
62 :
63 1 : auto dataRange = DataContainer(range);
64 1 : dataRange = 0;
65 :
66 1 : THEN("A^t A x should be close to the original data")
67 1 : {
68 1 : 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 1 : 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 1 : }
82 1 : }
83 :
84 6 : WHEN("We have a single ray with 180 degrees")
85 6 : {
86 1 : geom.emplace_back(stc, ctr, Degree{180}, std::move(volData), std::move(sinoData));
87 :
88 : // auto op = BinaryMethod(domain, range, geom);
89 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
90 1 : auto op = BinaryMethod(domain, range);
91 1 : auto dataRange = DataContainer(range);
92 1 : dataRange = 0;
93 :
94 1 : THEN("A^t A x should be close to the original data")
95 1 : {
96 1 : 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 1 : 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 1 : }
110 1 : }
111 :
112 6 : WHEN("We have a single ray with 90 degrees")
113 6 : {
114 1 : geom.emplace_back(stc, ctr, Degree{90}, std::move(volData), std::move(sinoData));
115 :
116 : // auto op = BinaryMethod(domain, range, geom);
117 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
118 1 : auto op = BinaryMethod(domain, range);
119 1 : auto dataRange = DataContainer(range);
120 1 : dataRange = 0;
121 :
122 1 : THEN("A^t A x should be close to the original data")
123 1 : {
124 1 : 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 1 : 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 1 : }
138 1 : }
139 :
140 6 : WHEN("We have a single ray with 270 degrees")
141 6 : {
142 1 : geom.emplace_back(stc, ctr, Degree{270}, std::move(volData), std::move(sinoData));
143 :
144 : // auto op = BinaryMethod(domain, range, geom);
145 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
146 1 : auto op = BinaryMethod(domain, range);
147 1 : auto dataRange = DataContainer(range);
148 1 : dataRange = 0;
149 :
150 1 : THEN("A^t A x should be close to the original data")
151 1 : {
152 1 : 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 1 : 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 1 : }
166 1 : }
167 :
168 6 : WHEN("We have a single ray with 45 degrees")
169 6 : {
170 1 : geom.emplace_back(stc, ctr, Degree{45}, std::move(volData), std::move(sinoData));
171 :
172 : // auto op = BinaryMethod(domain, range, geom);
173 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
174 1 : auto op = BinaryMethod(domain, range);
175 1 : auto dataRange = DataContainer(range);
176 1 : dataRange = 0;
177 :
178 1 : THEN("A^t A x should be close to the original data")
179 1 : {
180 1 : auto AtAx = DataContainer(domain);
181 :
182 1 : op.apply(dataDomain, dataRange);
183 1 : op.applyAdjoint(dataRange, AtAx);
184 :
185 1 : 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 1 : }
191 1 : }
192 :
193 6 : WHEN("We have a single ray with 225 degrees")
194 6 : {
195 1 : geom.emplace_back(stc, ctr, Degree{225}, std::move(volData), std::move(sinoData));
196 :
197 : // auto op = BinaryMethod(domain, range, geom);
198 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
199 1 : auto op = BinaryMethod(domain, range);
200 :
201 1 : 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 1 : THEN("A^t A x should be close to the original data")
209 1 : {
210 1 : 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 1 : 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 1 : }
223 1 : }
224 6 : }
225 6 : }
226 :
227 : TEST_CASE("BinaryMethod: Testing with only 1 rays for 4 angles")
228 2 : {
229 : // Turn logger of
230 2 : Logger::setLevel(Logger::LogLevel::OFF);
231 :
232 2 : IndexVector_t sizeDomain(2);
233 2 : sizeDomain << 5, 5;
234 :
235 2 : IndexVector_t sizeRange(2);
236 2 : sizeRange << 1, 4;
237 :
238 2 : auto domain = VolumeDescriptor(sizeDomain);
239 : // auto range = VolumeDescriptor(sizeRange);
240 :
241 2 : auto stc = SourceToCenterOfRotation{100};
242 2 : auto ctr = CenterOfRotationToDetector{5};
243 2 : auto volData = VolumeData2D{Size2D{sizeDomain}};
244 2 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
245 :
246 2 : Eigen::IOFormat CommaInitFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ", ", "", "",
247 2 : " << ", ";");
248 :
249 2 : GIVEN("A BinaryMethod for 2D and a domain data with all ones")
250 2 : {
251 2 : std::vector<Geometry> geom;
252 :
253 2 : auto dataDomain = DataContainer(domain);
254 2 : dataDomain = 1;
255 :
256 2 : WHEN("We have a single ray with 0, 90, 180, 270 degrees")
257 2 : {
258 2 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{sizeDomain}},
259 2 : SinogramData2D{Size2D{sizeRange}});
260 2 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{sizeDomain}},
261 2 : SinogramData2D{Size2D{sizeRange}});
262 2 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{sizeDomain}},
263 2 : SinogramData2D{Size2D{sizeRange}});
264 2 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{sizeDomain}},
265 2 : SinogramData2D{Size2D{sizeRange}});
266 :
267 : // auto op = BinaryMethod(domain, range, geom);
268 2 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
269 2 : auto op = BinaryMethod(domain, range);
270 :
271 2 : auto dataRange = DataContainer(range);
272 2 : dataRange = 0;
273 :
274 2 : THEN("A^t A x should be close to the original data")
275 2 : {
276 1 : 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 1 : auto cmp = RealVector_t(sizeDomain.prod());
284 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,
285 1 : 0, 0;
286 1 : DataContainer tmpCmp(domain, cmp);
287 :
288 1 : REQUIRE_UNARY(isCwiseApprox(tmpCmp, AtAx));
289 1 : }
290 :
291 2 : WHEN("A clone is created from the projector")
292 2 : {
293 1 : auto cloneOp = op.clone();
294 :
295 1 : THEN("The results will stay the same")
296 1 : {
297 1 : auto AtAx = DataContainer(domain);
298 :
299 1 : cloneOp->apply(dataDomain, dataRange);
300 1 : cloneOp->applyAdjoint(dataRange, AtAx);
301 :
302 1 : auto cmp = RealVector_t(sizeDomain.prod());
303 1 : cmp << 0, 0, 10, 0, 0, 0, 0, 10, 0, 0, 10, 10, 20, 10, 10, 0, 0, 10, 0, 0, 0, 0,
304 1 : 10, 0, 0;
305 1 : DataContainer tmpCmp(domain, cmp);
306 :
307 1 : REQUIRE_UNARY(isCwiseApprox(tmpCmp, AtAx));
308 1 : }
309 1 : }
310 2 : }
311 2 : }
312 2 : }
313 :
314 : TEST_CASE("BinaryMethod: Testing different setup")
315 6 : {
316 : // Turn logger of
317 6 : Logger::setLevel(Logger::LogLevel::OFF);
318 :
319 6 : IndexVector_t sizeDomain(2);
320 6 : sizeDomain << 5, 5;
321 :
322 6 : IndexVector_t sizeRange(2);
323 6 : sizeRange << 5, 1;
324 :
325 6 : auto domain = VolumeDescriptor(sizeDomain);
326 : // auto range = VolumeDescriptor(sizeRange);
327 :
328 6 : auto stc = SourceToCenterOfRotation{100};
329 6 : auto ctr = CenterOfRotationToDetector{5};
330 6 : auto volData = VolumeData2D{Size2D{sizeDomain}};
331 6 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
332 :
333 6 : Eigen::IOFormat CommaInitFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ", ", "", "",
334 6 : " << ", ";");
335 :
336 6 : GIVEN("A BinaryMethod with 1 angle at 0 degree")
337 6 : {
338 3 : 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 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
343 3 : auto op = BinaryMethod(domain, range);
344 :
345 : // THEN("It is not spd")
346 : // {
347 : // REQUIRE_FALSE(op.isSpd());
348 : // }
349 :
350 3 : THEN("Domain descriptor is still the same")
351 3 : {
352 1 : auto& retDescriptor = op.getDomainDescriptor();
353 :
354 1 : CHECK_EQ(retDescriptor.getNumberOfCoefficientsPerDimension()(0), 5);
355 1 : CHECK_EQ(retDescriptor.getNumberOfCoefficientsPerDimension()(1), 5);
356 1 : }
357 :
358 3 : THEN("Domain descriptor is still the same")
359 3 : {
360 1 : auto& retDescriptor = op.getRangeDescriptor();
361 :
362 1 : CHECK_EQ(retDescriptor.getNumberOfCoefficientsPerDimension()(0), 5);
363 1 : CHECK_EQ(retDescriptor.getNumberOfCoefficientsPerDimension()(1), 1);
364 1 : }
365 :
366 3 : WHEN("We have domain data with only ones")
367 3 : {
368 1 : auto dataDomain = DataContainer(domain);
369 1 : dataDomain = 1;
370 :
371 1 : auto dataRange = DataContainer(range);
372 1 : dataRange = 0;
373 :
374 1 : THEN("A^t A x should be close to the original data")
375 1 : {
376 1 : auto AtAx = DataContainer(domain);
377 :
378 1 : op.apply(dataDomain, dataRange);
379 :
380 1 : RealVector_t res = RealVector_t::Constant(sizeRange.prod(), 1, 5);
381 1 : DataContainer tmpRes(range, res);
382 1 : REQUIRE_UNARY(isCwiseApprox(tmpRes, dataRange));
383 :
384 1 : op.applyAdjoint(dataRange, AtAx);
385 :
386 1 : 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 1 : }
392 1 : }
393 3 : }
394 :
395 6 : GIVEN("A traversal with 5 rays at 180 degrees")
396 6 : {
397 1 : 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 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
402 1 : auto op = BinaryMethod(domain, range);
403 :
404 1 : THEN("A^t A x should be close to the original data")
405 1 : {
406 1 : auto dataDomain = DataContainer(domain);
407 1 : dataDomain = 1;
408 :
409 1 : auto dataRange = DataContainer(range);
410 1 : dataRange = 0;
411 :
412 1 : auto AtAx = DataContainer(domain);
413 :
414 1 : op.apply(dataDomain, dataRange);
415 :
416 1 : RealVector_t res = RealVector_t::Constant(sizeRange.prod(), 1, 5);
417 1 : DataContainer tmpRes(range, res);
418 1 : REQUIRE_UNARY(isCwiseApprox(tmpRes, dataRange));
419 :
420 1 : op.applyAdjoint(dataRange, AtAx);
421 :
422 1 : 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 1 : }
428 1 : }
429 :
430 6 : GIVEN("A traversal with 5 rays at 90 degrees")
431 6 : {
432 1 : 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 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
437 1 : auto op = BinaryMethod(domain, range);
438 :
439 1 : THEN("A^t A x should be close to the original data")
440 1 : {
441 1 : auto dataDomain = DataContainer(domain);
442 1 : dataDomain = 1;
443 :
444 1 : auto dataRange = DataContainer(range);
445 1 : dataRange = 0;
446 :
447 1 : auto AtAx = DataContainer(domain);
448 :
449 1 : op.apply(dataDomain, dataRange);
450 :
451 1 : RealVector_t res = RealVector_t::Constant(sizeRange.prod(), 1, 5);
452 1 : DataContainer tmpRes(range, res);
453 1 : REQUIRE_UNARY(isApprox(tmpRes, dataRange));
454 :
455 1 : op.applyAdjoint(dataRange, AtAx);
456 :
457 1 : 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 1 : }
463 1 : }
464 :
465 6 : GIVEN("A traversal with 5 rays at 270 degrees")
466 6 : {
467 1 : 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 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
472 1 : auto op = BinaryMethod(domain, range);
473 :
474 1 : THEN("A^t A x should be close to the original data")
475 1 : {
476 1 : auto dataDomain = DataContainer(domain);
477 1 : dataDomain = 1;
478 :
479 1 : auto dataRange = DataContainer(range);
480 1 : dataRange = 0;
481 :
482 1 : auto AtAx = DataContainer(domain);
483 :
484 1 : op.apply(dataDomain, dataRange);
485 :
486 1 : RealVector_t res = RealVector_t::Constant(sizeRange.prod(), 1, 5);
487 1 : DataContainer tmpRes(range, res);
488 1 : REQUIRE_UNARY(isApprox(tmpRes, dataRange));
489 :
490 1 : op.applyAdjoint(dataRange, AtAx);
491 :
492 1 : 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 1 : }
498 1 : }
499 6 : }
500 :
501 : TEST_CASE("BinaryMethod: Calls to functions of super class")
502 1 : {
503 : // Turn logger of
504 1 : Logger::setLevel(Logger::LogLevel::OFF);
505 :
506 1 : GIVEN("A projector")
507 1 : {
508 1 : 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 1 : VolumeDescriptor volumeDescriptor(volumeDims);
515 1 : VolumeDescriptor sinoDescriptor(sizeRange);
516 :
517 1 : RealVector_t randomStuff(volumeDescriptor.getNumberOfCoefficients());
518 1 : randomStuff.setRandom();
519 1 : DataContainer volume(volumeDescriptor, randomStuff);
520 1 : DataContainer sino(sinoDescriptor);
521 :
522 1 : auto stc = SourceToCenterOfRotation{20 * volSize};
523 1 : auto ctr = CenterOfRotationToDetector{volSize};
524 :
525 1 : 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 50 : geom.emplace_back(stc, ctr, Radian{angle}, VolumeData2D{Size2D{volumeDims}},
529 50 : SinogramData2D{Size2D{sizeRange}});
530 50 : }
531 :
532 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
533 1 : auto op = BinaryMethod(volumeDescriptor, range);
534 :
535 1 : WHEN("Projector is cloned")
536 1 : {
537 1 : auto opClone = op.clone();
538 1 : auto sinoClone = sino;
539 1 : auto volumeClone = volume;
540 :
541 1 : THEN("Results do not change (may still be slightly different due to summation being "
542 1 : "performed in a different order)")
543 1 : {
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 1 : }
552 1 : }
553 1 : }
554 1 : }
555 :
556 : TEST_CASE("BinaryMethod: Output DataContainer is not zero initialized")
557 4 : {
558 : // Turn logger of
559 4 : Logger::setLevel(Logger::LogLevel::OFF);
560 :
561 4 : GIVEN("A 2D setting")
562 4 : {
563 2 : 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 2 : VolumeDescriptor volumeDescriptor(volumeDims);
570 2 : VolumeDescriptor sinoDescriptor(sizeRange);
571 :
572 2 : DataContainer volume(volumeDescriptor);
573 2 : DataContainer sino(sinoDescriptor);
574 :
575 2 : auto stc = SourceToCenterOfRotation{20 * volSize};
576 2 : auto ctr = CenterOfRotationToDetector{volSize};
577 2 : auto volData = VolumeData2D{Size2D{volumeDims}};
578 2 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
579 :
580 2 : std::vector<Geometry> geom;
581 2 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData));
582 :
583 2 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
584 2 : auto op = BinaryMethod(volumeDescriptor, range);
585 :
586 2 : WHEN("Sinogram conatainer is not zero initialized and we project through an empty volume")
587 2 : {
588 1 : volume = 0;
589 1 : sino = 1;
590 :
591 1 : THEN("Result is zero")
592 1 : {
593 1 : op.apply(volume, sino);
594 :
595 1 : DataContainer zero(sinoDescriptor);
596 1 : zero = 0;
597 1 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
598 1 : }
599 1 : }
600 :
601 2 : WHEN("Volume container is not zero initialized and we backproject from an empty sinogram")
602 2 : {
603 1 : sino = 0;
604 1 : volume = 1;
605 :
606 1 : THEN("Result is zero")
607 1 : {
608 1 : op.applyAdjoint(sino, volume);
609 :
610 1 : DataContainer zero(volumeDescriptor);
611 1 : zero = 0;
612 1 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
613 1 : }
614 1 : }
615 2 : }
616 :
617 4 : GIVEN("A 3D setting")
618 4 : {
619 2 : 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 2 : VolumeDescriptor volumeDescriptor(volumeDims);
626 2 : VolumeDescriptor sinoDescriptor(sizeRange);
627 :
628 2 : auto stc = SourceToCenterOfRotation{20 * volSize};
629 2 : auto ctr = CenterOfRotationToDetector{volSize};
630 2 : auto volData = VolumeData3D{Size3D{volumeDims}};
631 2 : auto sinoData = SinogramData3D{Size3D{sizeRange}};
632 :
633 2 : DataContainer volume(volumeDescriptor);
634 2 : DataContainer sino(sinoDescriptor);
635 :
636 2 : std::vector<Geometry> geom;
637 2 : geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
638 2 : RotationAngles3D{Gamma{0}});
639 :
640 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
641 2 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
642 2 : auto op = BinaryMethod(volumeDescriptor, range);
643 :
644 2 : WHEN("Sinogram conatainer is not zero initialized and we project through an empty volume")
645 2 : {
646 1 : volume = 0;
647 1 : sino = 1;
648 :
649 1 : THEN("Result is zero")
650 1 : {
651 1 : op.apply(volume, sino);
652 :
653 1 : DataContainer zero(sinoDescriptor);
654 1 : zero = 0;
655 1 : REQUIRE_UNARY(isCwiseApprox(sino, zero));
656 1 : }
657 1 : }
658 :
659 2 : WHEN("Volume container is not zero initialized and we backproject from an empty sinogram")
660 2 : {
661 1 : sino = 0;
662 1 : volume = 1;
663 :
664 1 : THEN("Result is zero")
665 1 : {
666 1 : op.applyAdjoint(sino, volume);
667 :
668 1 : DataContainer zero(volumeDescriptor);
669 1 : zero = 0;
670 1 : REQUIRE_UNARY(isCwiseApprox(volume, zero));
671 1 : }
672 1 : }
673 2 : }
674 4 : }
675 :
676 : TEST_CASE("BinaryMethod: Rays not intersecting the bounding box are present")
677 6 : {
678 : // Turn logger of
679 6 : Logger::setLevel(Logger::LogLevel::OFF);
680 :
681 6 : GIVEN("A 2D setting")
682 6 : {
683 4 : 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 4 : VolumeDescriptor volumeDescriptor(volumeDims);
690 4 : VolumeDescriptor sinoDescriptor(sizeRange);
691 :
692 4 : auto stc = SourceToCenterOfRotation{20 * volSize};
693 4 : auto ctr = CenterOfRotationToDetector{volSize};
694 4 : auto volData = VolumeData2D{Size2D{volumeDims}};
695 4 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
696 :
697 4 : DataContainer volume(volumeDescriptor);
698 4 : DataContainer sino(sinoDescriptor);
699 4 : volume = 1;
700 4 : sino = 0;
701 :
702 4 : std::vector<Geometry> geom;
703 :
704 4 : WHEN("Tracing along a y-axis-aligned ray with a negative x-coordinate of origin")
705 4 : {
706 1 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
707 1 : PrincipalPointOffset{}, RotationOffset2D{-volSize, 0});
708 :
709 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
710 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
711 1 : auto op = BinaryMethod(volumeDescriptor, range);
712 :
713 1 : THEN("Result of forward projection is zero")
714 1 : {
715 1 : op.apply(volume, sino);
716 :
717 1 : DataContainer zeroSino(sinoDescriptor);
718 1 : zeroSino = 0;
719 1 : REQUIRE_UNARY(isCwiseApprox(sino, zeroSino));
720 :
721 1 : AND_THEN("Result of backprojection is zero")
722 1 : {
723 1 : op.applyAdjoint(sino, volume);
724 :
725 1 : DataContainer zeroVolume(volumeDescriptor);
726 1 : zeroVolume = 0;
727 1 : REQUIRE_UNARY(isCwiseApprox(volume, zeroVolume));
728 1 : }
729 1 : }
730 1 : }
731 :
732 4 : WHEN("Tracing along a y-axis-aligned ray with a x-coordinate of origin beyond the bounding "
733 4 : "box")
734 4 : {
735 1 : geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
736 1 : PrincipalPointOffset{}, RotationOffset2D{volSize, 0});
737 :
738 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
739 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
740 1 : auto op = BinaryMethod(volumeDescriptor, range);
741 :
742 1 : THEN("Result of forward projection is zero")
743 1 : {
744 1 : op.apply(volume, sino);
745 :
746 1 : DataContainer zeroSino(sinoDescriptor);
747 1 : zeroSino = 0;
748 1 : REQUIRE_UNARY(isCwiseApprox(sino, zeroSino));
749 :
750 1 : AND_THEN("Result of backprojection is zero")
751 1 : {
752 1 : op.applyAdjoint(sino, volume);
753 :
754 1 : DataContainer zeroVolume(volumeDescriptor);
755 1 : zeroVolume = 0;
756 1 : REQUIRE_UNARY(isCwiseApprox(volume, zeroVolume));
757 1 : }
758 1 : }
759 1 : }
760 :
761 4 : WHEN("Tracing along a x-axis-aligned ray with a negative y-coordinate of origin")
762 4 : {
763 1 : geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData), std::move(sinoData),
764 1 : PrincipalPointOffset{}, RotationOffset2D{0, -volSize});
765 :
766 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
767 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
768 1 : auto op = BinaryMethod(volumeDescriptor, range);
769 :
770 1 : THEN("Result of forward projection is zero")
771 1 : {
772 1 : op.apply(volume, sino);
773 :
774 1 : DataContainer zeroSino(sinoDescriptor);
775 1 : zeroSino = 0;
776 1 : REQUIRE_UNARY(isCwiseApprox(sino, zeroSino));
777 :
778 1 : AND_THEN("Result of backprojection is zero")
779 1 : {
780 1 : op.applyAdjoint(sino, volume);
781 :
782 1 : DataContainer zeroVolume(volumeDescriptor);
783 1 : zeroVolume = 0;
784 1 : REQUIRE_UNARY(isCwiseApprox(volume, zeroVolume));
785 1 : }
786 1 : }
787 1 : }
788 :
789 4 : WHEN("Tracing along a x-axis-aligned ray with a y-coordinate of origin beyond the bounding "
790 4 : "box")
791 4 : {
792 1 : geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData), std::move(sinoData),
793 1 : PrincipalPointOffset{}, RotationOffset2D{0, volSize});
794 :
795 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
796 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
797 1 : auto op = BinaryMethod(volumeDescriptor, range);
798 :
799 1 : THEN("Result of forward projection is zero")
800 1 : {
801 1 : op.apply(volume, sino);
802 :
803 1 : DataContainer zeroSino(sinoDescriptor);
804 1 : zeroSino = 0;
805 1 : REQUIRE_UNARY(isCwiseApprox(sino, zeroSino));
806 :
807 1 : AND_THEN("Result of backprojection is zero")
808 1 : {
809 1 : op.applyAdjoint(sino, volume);
810 :
811 1 : DataContainer zeroVolume(volumeDescriptor);
812 1 : zeroVolume = 0;
813 1 : REQUIRE_UNARY(isCwiseApprox(volume, zeroVolume));
814 1 : }
815 1 : }
816 1 : }
817 4 : }
818 :
819 6 : GIVEN("A 3D setting")
820 6 : {
821 2 : 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 2 : VolumeDescriptor volumeDescriptor(volumeDims);
828 2 : VolumeDescriptor sinoDescriptor(sizeRange);
829 :
830 2 : DataContainer volume(volumeDescriptor);
831 2 : 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 2 : auto volData = VolumeData3D{Size3D{volumeDims}};
838 2 : auto sinoData = SinogramData3D{Size3D{sizeRange}};
839 :
840 2 : 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 2 : 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 2 : std::string neg[numCases] = {"x", "y", "x and y", "y", "z", "y and z", "x", "z", "x and z"};
851 2 : std::string ali[numCases] = {"z", "z", "z", "x", "x", "x", "y", "y", "y"};
852 :
853 20 : for (int i = 0; i < numCases; i++) {
854 18 : WHEN("Tracing rays along different axis")
855 18 : {
856 1 : INFO("Tracing along a ", ali[i], "-axis-aligned ray with negative ", neg[i],
857 1 : "-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 1 : PrincipalPointOffset2D{0, 0},
861 1 : RotationOffset3D{-offsetx[i], -offsety[i], -offsetz[i]});
862 :
863 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
864 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
865 1 : auto op = BinaryMethod(volumeDescriptor, range);
866 :
867 1 : THEN("Result of forward projection is zero")
868 1 : {
869 1 : op.apply(volume, sino);
870 :
871 1 : DataContainer zeroSino(sinoDescriptor);
872 1 : zeroSino = 0;
873 1 : REQUIRE_UNARY(isCwiseApprox(sino, zeroSino));
874 :
875 1 : AND_THEN("Result of backprojection is zero")
876 1 : {
877 1 : op.applyAdjoint(sino, volume);
878 :
879 1 : DataContainer zeroVolume(volumeDescriptor);
880 1 : zeroVolume = 0;
881 1 : REQUIRE_UNARY(isCwiseApprox(volume, zeroVolume));
882 1 : }
883 1 : }
884 1 : }
885 18 : }
886 2 : }
887 6 : }
888 :
889 : TEST_CASE("BinaryMethod: Axis-aligned rays are present")
890 10 : {
891 : // Turn logger of
892 10 : Logger::setLevel(Logger::LogLevel::OFF);
893 :
894 10 : GIVEN("A 2D setting with a single ray")
895 10 : {
896 4 : 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 4 : VolumeDescriptor volumeDescriptor(volumeDims);
903 4 : VolumeDescriptor sinoDescriptor(sizeRange);
904 :
905 4 : DataContainer volume(volumeDescriptor);
906 4 : DataContainer sino(sinoDescriptor);
907 :
908 4 : auto stc = SourceToCenterOfRotation{20 * volSize};
909 4 : auto ctr = CenterOfRotationToDetector{volSize};
910 4 : auto volData = VolumeData2D{Size2D{volumeDims}};
911 4 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
912 :
913 4 : 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 4 : 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 16 : WHEN("Axis-aligned ray through the center of the pixel")
926 16 : {
927 1 : INFO("An axis-aligned ray with an angle of ", angles[i],
928 1 : " radians passes through the center of a pixel");
929 :
930 1 : geom.emplace_back(stc, ctr, Radian{angles[i]}, std::move(volData),
931 1 : std::move(sinoData));
932 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
933 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
934 1 : auto op = BinaryMethod(volumeDescriptor, range);
935 1 : THEN("The result of projecting through a pixel is exactly the pixel value")
936 1 : {
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 5 : } else {
944 0 : IndexVector_t coord(2);
945 0 : coord << j, volSize / 2;
946 0 : volume[volumeDescriptor.getIndexFromCoordinate(coord)] = 1;
947 0 : }
948 :
949 5 : op.apply(volume, sino);
950 5 : REQUIRE_EQ(sino[0], Approx(1));
951 5 : }
952 :
953 1 : AND_THEN("The backprojection sets the values of all hit pixels to the detector "
954 1 : "value")
955 1 : {
956 1 : op.applyAdjoint(sino, volume);
957 :
958 1 : DataContainer res(volumeDescriptor, backProj[i % 2]);
959 1 : REQUIRE_UNARY(isCwiseApprox(volume, res));
960 1 : }
961 1 : }
962 1 : }
963 16 : }
964 :
965 4 : WHEN("A y-axis-aligned ray runs along a voxel boundary")
966 4 : {
967 1 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, Radian{0},
968 1 : std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
969 1 : RotationOffset2D{-0.5, 0});
970 :
971 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
972 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
973 1 : auto op = BinaryMethod(volumeDescriptor, range);
974 1 : THEN("The result of projecting through a pixel is the value of the pixel with the "
975 1 : "higher index")
976 1 : {
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 5 : }
986 :
987 1 : AND_THEN("The backprojection yields the exact adjoint")
988 1 : {
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 1 : }
995 1 : }
996 1 : }
997 :
998 4 : WHEN("A y-axis-aligned ray runs along the right volume boundary")
999 4 : {
1000 : // For Siddon's values in the range [0,boxMax) are considered, a ray running along
1001 : // boxMax should be ignored
1002 :
1003 1 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, Radian{0},
1004 1 : std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
1005 1 : RotationOffset2D{volSize * 0.5, 0});
1006 :
1007 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1008 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1009 1 : auto op = BinaryMethod(volumeDescriptor, range);
1010 :
1011 1 : THEN("The result of projecting is zero")
1012 1 : {
1013 1 : volume = 0;
1014 1 : op.apply(volume, sino);
1015 1 : REQUIRE_EQ(sino[0], Approx(0.0));
1016 :
1017 1 : AND_THEN("The result of backprojection is also zero")
1018 1 : {
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 1 : }
1027 1 : }
1028 1 : }
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 4 : WHEN("A y-axis-aligned ray runs along the left volume boundary")
1033 4 : {
1034 1 : geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, Radian{0},
1035 1 : std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
1036 1 : RotationOffset2D{-volSize / 2.0, 0});
1037 :
1038 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1039 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1040 1 : auto op = BinaryMethod(volumeDescriptor, range);
1041 1 : THEN("The result of projecting through a pixel is exactly the pixel's value")
1042 1 : {
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 5 : }
1050 :
1051 1 : AND_THEN("The backprojection yields the exact adjoint")
1052 1 : {
1053 1 : sino[0] = 1;
1054 1 : op.applyAdjoint(sino, volume);
1055 1 : REQUIRE_UNARY(
1056 1 : isCwiseApprox(volume, DataContainer(volumeDescriptor, backProj[0])));
1057 1 : }
1058 1 : }
1059 1 : }
1060 4 : }
1061 :
1062 10 : GIVEN("A 3D setting with a single ray")
1063 10 : {
1064 4 : 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 4 : VolumeDescriptor volumeDescriptor(volumeDims);
1071 4 : VolumeDescriptor sinoDescriptor(sizeRange);
1072 :
1073 4 : DataContainer volume(volumeDescriptor);
1074 4 : DataContainer sino(sinoDescriptor);
1075 :
1076 4 : auto stc = SourceToCenterOfRotation{20 * volSize};
1077 4 : auto ctr = CenterOfRotationToDetector{volSize};
1078 4 : auto volData = VolumeData3D{Size3D{volumeDims}};
1079 4 : auto sinoData = SinogramData3D{Size3D{sizeRange}};
1080 :
1081 4 : 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 4 : std::string al[numCases] = {"z", "-z", "x", "-x", "y", "-y"};
1087 :
1088 4 : RealVector_t backProj[numCases];
1089 4 : for (auto& backPr : backProj)
1090 24 : backPr.resize(volSize * volSize * volSize);
1091 :
1092 4 : 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 4 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1097 :
1098 4 : 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 4 : 0, 0, 0, 0, 0, 0, 0, 0, 0;
1103 :
1104 4 : 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 4 : 0, 0, 0, 0, 1, 0, 0, 0, 0;
1109 :
1110 28 : for (index_t i = 0; i < numCases; i++) {
1111 24 : WHEN("Tracing an axis-aligned ray trough the pixel center")
1112 24 : {
1113 1 : 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 1 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
1116 :
1117 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1118 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1119 1 : auto op = BinaryMethod(volumeDescriptor, range);
1120 1 : THEN("The result of projecting through a voxel is exactly the voxel value")
1121 1 : {
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 3 : } 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 0 : }
1137 :
1138 3 : op.apply(volume, sino);
1139 3 : REQUIRE_EQ(sino[0], Approx(1));
1140 3 : }
1141 :
1142 1 : AND_THEN("The backprojection sets the values of all hit voxels to the detector "
1143 1 : "value")
1144 1 : {
1145 1 : op.applyAdjoint(sino, volume);
1146 :
1147 1 : DataContainer res(volumeDescriptor, backProj[i / 2]);
1148 1 : REQUIRE_UNARY(isCwiseApprox(volume, res));
1149 1 : }
1150 1 : }
1151 1 : }
1152 24 : }
1153 :
1154 4 : real_t offsetx[numCases];
1155 4 : 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 4 : 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 4 : 0, 0, 0, 1, 0, 0, 0, 0, 0;
1176 :
1177 4 : 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 4 : 0, 1, 0, 0, 0, 0, 0, 0, 0;
1182 :
1183 4 : 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 4 : 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 12 : WHEN("A z-axis-aligned ray runs along the corners and edges of the volume")
1198 12 : {
1199 1 : 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 1 : std::move(sinoData), RotationAngles3D{Gamma{0}},
1204 1 : PrincipalPointOffset2D{0, 0},
1205 1 : RotationOffset3D{-offsetx[i], -offsety[i], 0});
1206 :
1207 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1208 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1209 1 : auto op = BinaryMethod(volumeDescriptor, range);
1210 1 : THEN("The result of projecting through a voxel is exactly the voxel's value")
1211 1 : {
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 0 : break;
1222 0 : case 2:
1223 0 : volume(0, 0, j) = 1;
1224 0 : break;
1225 0 : default:
1226 0 : break;
1227 3 : }
1228 :
1229 3 : op.apply(volume, sino);
1230 3 : REQUIRE_EQ(sino[0], Approx(1));
1231 3 : }
1232 :
1233 1 : AND_THEN("The backprojection yields the exact adjoints")
1234 1 : {
1235 1 : sino[0] = 1;
1236 1 : op.applyAdjoint(sino, volume);
1237 :
1238 1 : REQUIRE_UNARY(
1239 1 : isCwiseApprox(volume, DataContainer(volumeDescriptor, backProj[i])));
1240 1 : }
1241 1 : }
1242 1 : }
1243 12 : }
1244 :
1245 16 : for (index_t i = numCases / 2; i < numCases; i++) {
1246 12 : WHEN("A z-axis-aligned ray runs along the edges and corners of the volume")
1247 12 : {
1248 1 : 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 1 : std::move(sinoData), RotationAngles3D{Gamma{0}},
1253 1 : PrincipalPointOffset2D{0, 0},
1254 1 : RotationOffset3D{-offsetx[i], -offsety[i], 0});
1255 :
1256 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1257 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1258 1 : auto op = BinaryMethod(volumeDescriptor, range);
1259 1 : THEN("The result of projecting is zero")
1260 1 : {
1261 1 : volume = 1;
1262 :
1263 1 : op.apply(volume, sino);
1264 1 : REQUIRE_EQ(sino[0], Approx(0));
1265 :
1266 1 : AND_THEN("The result of backprojection is also zero")
1267 1 : {
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 1 : }
1274 1 : }
1275 1 : }
1276 12 : }
1277 4 : }
1278 :
1279 10 : GIVEN("A 2D setting with multiple projection angles")
1280 10 : {
1281 1 : 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 1 : VolumeDescriptor volumeDescriptor(volumeDims);
1288 1 : VolumeDescriptor sinoDescriptor(sizeRange);
1289 1 : DataContainer volume(volumeDescriptor);
1290 1 : DataContainer sino(sinoDescriptor);
1291 :
1292 1 : auto stc = SourceToCenterOfRotation{20 * volSize};
1293 1 : auto ctr = CenterOfRotationToDetector{volSize};
1294 :
1295 1 : std::vector<Geometry> geom;
1296 :
1297 1 : WHEN("Both x- and y-axis-aligned rays are present")
1298 1 : {
1299 1 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{volumeDims}},
1300 1 : SinogramData2D{Size2D{sizeRange}});
1301 1 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{volumeDims}},
1302 1 : SinogramData2D{Size2D{sizeRange}});
1303 1 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{volumeDims}},
1304 1 : SinogramData2D{Size2D{sizeRange}});
1305 1 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{volumeDims}},
1306 1 : SinogramData2D{Size2D{sizeRange}});
1307 :
1308 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1309 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1310 1 : auto op = BinaryMethod(volumeDescriptor, range);
1311 :
1312 1 : THEN("Values are accumulated correctly along each ray's path")
1313 1 : {
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 5 : }
1325 :
1326 1 : op.apply(volume, sino);
1327 5 : for (index_t i = 0; i < numImgs; i++)
1328 1 : REQUIRE_EQ(sino[i], Approx(5.0));
1329 :
1330 1 : AND_THEN("Backprojection yields the exact adjoint")
1331 1 : {
1332 1 : RealVector_t cmp(volSize * volSize);
1333 :
1334 1 : cmp << 0, 0, 10, 0, 0, 0, 0, 10, 0, 0, 10, 10, 20, 10, 10, 0, 0, 10, 0, 0, 0, 0,
1335 1 : 10, 0, 0;
1336 1 : DataContainer tmpCmp(volumeDescriptor, cmp);
1337 :
1338 1 : op.applyAdjoint(sino, volume);
1339 1 : REQUIRE_UNARY(isCwiseApprox(volume, tmpCmp));
1340 1 : }
1341 1 : }
1342 1 : }
1343 1 : }
1344 :
1345 10 : GIVEN("A 3D setting with multiple projection angles")
1346 10 : {
1347 1 : 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 1 : VolumeDescriptor volumeDescriptor(volumeDims);
1354 1 : VolumeDescriptor sinoDescriptor(sizeRange);
1355 1 : DataContainer volume(volumeDescriptor);
1356 1 : DataContainer sino(sinoDescriptor);
1357 :
1358 1 : auto stc = SourceToCenterOfRotation{20 * volSize};
1359 1 : auto ctr = CenterOfRotationToDetector{volSize};
1360 :
1361 1 : std::vector<Geometry> geom;
1362 :
1363 1 : WHEN("x-, y and z-axis-aligned rays are present")
1364 1 : {
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 6 : geom.emplace_back(stc, ctr, VolumeData3D{Size3D{volumeDims}},
1370 6 : SinogramData3D{Size3D{sizeRange}},
1371 6 : RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
1372 :
1373 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1374 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1375 1 : auto op = BinaryMethod(volumeDescriptor, range);
1376 :
1377 1 : THEN("Values are accumulated correctly along each ray's path")
1378 1 : {
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 3 : }
1392 :
1393 1 : op.apply(volume, sino);
1394 7 : for (index_t i = 0; i < numImgs; i++)
1395 1 : REQUIRE_EQ(sino[i], Approx(3.0));
1396 :
1397 1 : AND_THEN("Backprojection yields the exact adjoint")
1398 1 : {
1399 1 : RealVector_t cmp(volSize * volSize * volSize);
1400 :
1401 1 : cmp << 0, 0, 0, 0, 6, 0, 0, 0, 0,
1402 :
1403 1 : 0, 6, 0, 6, 18, 6, 0, 6, 0,
1404 :
1405 1 : 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 1 : }
1411 1 : }
1412 1 : }
1413 1 : }
1414 10 : }
1415 :
1416 : TEST_CASE("BinaryMethod: Projection under an angle")
1417 12 : {
1418 : // Turn logger of
1419 12 : Logger::setLevel(Logger::LogLevel::OFF);
1420 :
1421 12 : GIVEN("A 2D setting with a single ray")
1422 12 : {
1423 8 : 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 8 : VolumeDescriptor volumeDescriptor(volumeDims);
1430 8 : VolumeDescriptor sinoDescriptor(sizeRange);
1431 :
1432 8 : DataContainer volume(volumeDescriptor);
1433 8 : DataContainer sino(sinoDescriptor);
1434 :
1435 8 : auto stc = SourceToCenterOfRotation{20 * volSize};
1436 8 : auto ctr = CenterOfRotationToDetector{volSize};
1437 8 : auto volData = VolumeData2D{Size2D{volumeDims}};
1438 8 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
1439 :
1440 8 : std::vector<Geometry> geom;
1441 :
1442 8 : WHEN("Projecting under an angle of 30 degrees and ray goes through center of volume")
1443 8 : {
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 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1449 1 : auto op = BinaryMethod(volumeDescriptor, range);
1450 :
1451 1 : THEN("Ray intersects the correct pixels")
1452 1 : {
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 1 : DataContainer sZero(sinoDescriptor);
1471 1 : sZero = 0;
1472 1 : CHECK_LE(std::abs(sino[0]), Approx(0.0001f).epsilon(epsilon));
1473 :
1474 1 : AND_THEN("The correct weighting is applied")
1475 1 : {
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 1 : 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 1 : }
1502 1 : }
1503 1 : }
1504 :
1505 8 : WHEN("Projecting under an angle of 30 degrees and ray enters through the right border")
1506 8 : {
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 1 : PrincipalPointOffset{0}, RotationOffset2D{std::sqrt(3.f), 0});
1511 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1512 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1513 1 : auto op = BinaryMethod(volumeDescriptor, range);
1514 :
1515 1 : THEN("Ray intersects the correct pixels")
1516 1 : {
1517 1 : volume = 0;
1518 1 : 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 1 : DataContainer sZero(sinoDescriptor);
1530 1 : sZero = 0;
1531 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1532 :
1533 1 : AND_THEN("The correct weighting is applied")
1534 1 : {
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 1 : 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 1 : }
1556 1 : }
1557 1 : }
1558 :
1559 8 : WHEN("Projecting under an angle of 30 degrees and ray exits through the left border")
1560 8 : {
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 1 : PrincipalPointOffset{0}, RotationOffset2D{-std::sqrt(3.f), 0});
1565 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1566 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1567 1 : auto op = BinaryMethod(volumeDescriptor, range);
1568 :
1569 1 : THEN("Ray intersects the correct pixels")
1570 1 : {
1571 1 : volume = 1;
1572 1 : 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 1 : DataContainer sZero(sinoDescriptor);
1584 1 : sZero = 0;
1585 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1586 :
1587 1 : AND_THEN("The correct weighting is applied")
1588 1 : {
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 1 : 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 1 : }
1610 1 : }
1611 1 : }
1612 :
1613 8 : WHEN("Projecting under an angle of 30 degrees and ray only intersects a single pixel")
1614 8 : {
1615 1 : geom.emplace_back(stc, ctr, Radian{-pi_t / 6}, std::move(volData), std::move(sinoData),
1616 1 : PrincipalPointOffset{0},
1617 1 : RotationOffset2D{-2 - std::sqrt(3.f) / 2, 0});
1618 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1619 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1620 1 : auto op = BinaryMethod(volumeDescriptor, range);
1621 :
1622 1 : THEN("Ray intersects the correct pixels")
1623 1 : {
1624 1 : volume = 1;
1625 1 : IndexVector_t coord(2);
1626 1 : coord << 0, 0;
1627 1 : volume[volumeDescriptor.getIndexFromCoordinate(coord)] = 0;
1628 :
1629 1 : op.apply(volume, sino);
1630 1 : DataContainer sZero(sinoDescriptor);
1631 1 : sZero = 0;
1632 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1633 :
1634 1 : AND_THEN("The correct weighting is applied")
1635 1 : {
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 1 : 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 1 : }
1651 1 : }
1652 1 : }
1653 :
1654 8 : WHEN("Projecting under an angle of 120 degrees and ray goes through center of volume")
1655 8 : {
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 1 : std::move(sinoData));
1660 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1661 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1662 1 : auto op = BinaryMethod(volumeDescriptor, range);
1663 :
1664 1 : THEN("Ray intersects the correct pixels")
1665 1 : {
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 1 : DataContainer sZero(sinoDescriptor);
1677 1 : sZero = 0;
1678 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1679 :
1680 1 : AND_THEN("The correct weighting is applied")
1681 1 : {
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 1 : 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 1 : }
1709 1 : }
1710 1 : }
1711 :
1712 8 : WHEN("Projecting under an angle of 120 degrees and ray enters through the top border")
1713 8 : {
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 1 : RotationOffset2D{0, std::sqrt(3.f)});
1719 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1720 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1721 1 : auto op = BinaryMethod(volumeDescriptor, range);
1722 :
1723 1 : THEN("Ray intersects the correct pixels")
1724 1 : {
1725 1 : volume = 1;
1726 1 : 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 1 : DataContainer sZero(sinoDescriptor);
1738 1 : sZero = 0;
1739 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1740 :
1741 1 : AND_THEN("The correct weighting is applied")
1742 1 : {
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 1 : 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 1 : }
1765 1 : }
1766 1 : }
1767 :
1768 8 : WHEN("Projecting under an angle of 120 degrees and ray exits through the bottom border")
1769 8 : {
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 1 : RotationOffset2D{0, -std::sqrt(3.f)});
1776 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1777 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1778 1 : auto op = BinaryMethod(volumeDescriptor, range);
1779 :
1780 1 : THEN("Ray intersects the correct pixels")
1781 1 : {
1782 1 : volume = 1;
1783 1 : 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 1 : DataContainer sZero(sinoDescriptor);
1795 1 : sZero = 0;
1796 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1797 :
1798 1 : AND_THEN("The correct weighting is applied")
1799 1 : {
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 1 : 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 1 : }
1822 1 : }
1823 1 : }
1824 :
1825 8 : WHEN("Projecting under an angle of 120 degrees and ray only intersects a single pixel")
1826 8 : {
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 1 : RotationOffset2D{0, -2 - std::sqrt(3.f) / 2});
1832 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1833 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1834 1 : auto op = BinaryMethod(volumeDescriptor, range);
1835 :
1836 1 : THEN("Ray intersects the correct pixels")
1837 1 : {
1838 1 : volume = 1;
1839 1 : IndexVector_t coord(2);
1840 1 : coord << 3, 0;
1841 1 : volume[volumeDescriptor.getIndexFromCoordinate(coord)] = 0;
1842 :
1843 1 : op.apply(volume, sino);
1844 1 : DataContainer sZero(sinoDescriptor);
1845 1 : sZero = 0;
1846 1 : CHECK_UNARY(isApprox(sino, sZero, epsilon));
1847 :
1848 1 : AND_THEN("The correct weighting is applied")
1849 1 : {
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 1 : 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 1 : }
1865 1 : }
1866 1 : }
1867 8 : }
1868 :
1869 12 : GIVEN("A 3D setting with a single ray")
1870 12 : {
1871 4 : 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 4 : VolumeDescriptor volumeDescriptor(volumeDims);
1878 4 : VolumeDescriptor sinoDescriptor(sizeRange);
1879 4 : DataContainer volume(volumeDescriptor);
1880 4 : DataContainer sino(sinoDescriptor);
1881 :
1882 4 : auto stc = SourceToCenterOfRotation{20 * volSize};
1883 4 : auto ctr = CenterOfRotationToDetector{volSize};
1884 4 : auto volData = VolumeData3D{Size3D{volumeDims}};
1885 4 : auto sinoData = SinogramData3D{Size3D{sizeRange}};
1886 :
1887 4 : std::vector<Geometry> geom;
1888 :
1889 4 : RealVector_t backProj(volSize * volSize * volSize);
1890 :
1891 4 : WHEN("A ray with an angle of 30 degrees goes through the center of the volume")
1892 4 : {
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 1 : RotationAngles3D{Gamma{pi_t / 6}});
1896 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1897 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1898 1 : auto op = BinaryMethod(volumeDescriptor, range);
1899 :
1900 1 : THEN("The ray intersects the correct voxels")
1901 1 : {
1902 1 : volume = 1;
1903 1 : 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 1 : AND_THEN("The correct weighting is applied")
1919 1 : {
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 1 : backProj << 0, 0, 0, 0, 1, 1, 0, 0, 0,
1932 :
1933 1 : 0, 0, 0, 0, 1, 0, 0, 0, 0,
1934 :
1935 1 : 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 1 : }
1941 1 : }
1942 1 : }
1943 :
1944 4 : WHEN("A ray with an angle of 30 degrees enters through the right border")
1945 4 : {
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 1 : RotationOffset3D{1, 0, 0});
1951 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
1952 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1953 1 : auto op = BinaryMethod(volumeDescriptor, range);
1954 :
1955 1 : THEN("The ray intersects the correct voxels")
1956 1 : {
1957 1 : volume = 1;
1958 1 : 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 1 : AND_THEN("The correct weighting is applied")
1972 1 : {
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 1 : backProj << 0, 0, 0, 0, 0, 1, 0, 0, 0,
1985 :
1986 1 : 0, 0, 0, 0, 0, 1, 0, 0, 0,
1987 :
1988 1 : 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 1 : }
1994 1 : }
1995 1 : }
1996 :
1997 4 : WHEN("A ray with an angle of 30 degrees exits through the left border")
1998 4 : {
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 1 : RotationOffset3D{-1, 0, 0});
2003 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
2004 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
2005 1 : auto op = BinaryMethod(volumeDescriptor, range);
2006 :
2007 1 : THEN("The ray intersects the correct voxels")
2008 1 : {
2009 1 : volume = 1;
2010 1 : 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 1 : AND_THEN("The correct weighting is applied")
2024 1 : {
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 1 : backProj << 0, 0, 0, 1, 1, 0, 0, 0, 0,
2037 :
2038 1 : 0, 0, 0, 1, 0, 0, 0, 0, 0,
2039 :
2040 1 : 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 1 : }
2046 1 : }
2047 1 : }
2048 :
2049 4 : WHEN("A ray with an angle of 30 degrees only intersects a single voxel")
2050 4 : {
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 1 : RotationOffset3D{-2, 0, 0});
2055 : // BinaryMethod op(volumeDescriptor, sinoDescriptor, geom);
2056 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
2057 1 : auto op = BinaryMethod(volumeDescriptor, range);
2058 :
2059 1 : THEN("The ray intersects the correct voxels")
2060 1 : {
2061 1 : volume = 1;
2062 1 : 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 1 : AND_THEN("The correct weighting is applied")
2070 1 : {
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 1 : backProj << 0, 0, 0, 1, 0, 0, 0, 0, 0,
2079 :
2080 1 : 0, 0, 0, 0, 0, 0, 0, 0, 0,
2081 :
2082 1 : 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 1 : }
2088 1 : }
2089 1 : }
2090 4 : }
2091 12 : }
|
