Line data Source code
1 : #include "doctest/doctest.h"
2 :
3 : #include "VoxelProjector.h"
4 : #include "PlanarDetectorDescriptor.h"
5 :
6 : #include "PrettyPrint/Eigen.h"
7 : #include "PrettyPrint/Stl.h"
8 : #include "spdlog/fmt/bundled/core.h"
9 :
10 : using namespace elsa;
11 : using namespace elsa::geometry;
12 : using namespace doctest;
13 :
14 : TEST_SUITE_BEGIN("projectors");
15 :
16 : Eigen::IOFormat vecfmt(10, 0, ", ", ", ", "", "", "[", "]");
17 : Eigen::IOFormat matfmt(10, 0, ", ", "\n", "\t\t[", "]");
18 :
19 : TYPE_TO_STRING(BlobVoxelProjector<float>);
20 :
21 : // Redefine GIVEN such that it's nicely usable inside an loop
22 : #undef GIVEN
23 137813 : #define GIVEN(...) DOCTEST_SUBCASE((std::string(" Given: ") + std::string(__VA_ARGS__)).c_str())
24 :
25 : TEST_CASE_TEMPLATE("BlobVoxelProjector: Testing simple volume 3D with one detector pixel", data_t,
26 : float, double)
27 540 : {
28 540 : const IndexVector_t sizeDomain({{5, 5, 5}});
29 540 : const IndexVector_t sizeRange({{1, 1, 1}});
30 :
31 540 : auto stc = SourceToCenterOfRotation{100};
32 540 : auto ctr = CenterOfRotationToDetector{5};
33 :
34 2160 : for (int i = 1; i < 4; i++) {
35 1620 : real_t scaling = i / 2.0f;
36 1620 : const RealVector_t spacing{{scaling, scaling, scaling}};
37 1620 : GIVEN("Spacing of " + std::to_string(scaling))
38 49140 : for (int i = 0; i < 360; i += 4) {
39 48600 : GIVEN("Ray of angle " + std::to_string(i))
40 48600 : {
41 540 : auto domain = VolumeDescriptor(sizeDomain, spacing);
42 540 : auto x = DataContainer<data_t>(domain);
43 540 : x = 0;
44 : // set center voxel to 1
45 540 : x(2, 2, 2) = 1;
46 540 : std::vector<Geometry> geom;
47 540 : geom.emplace_back(stc, ctr, VolumeData3D{Size3D{sizeDomain}, Spacing3D{spacing}},
48 540 : SinogramData3D{Size3D{sizeRange}},
49 540 : RotationAngles3D{Gamma{static_cast<real_t>(i)}});
50 540 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
51 540 : auto op = BlobVoxelProjector<data_t>(domain, range);
52 :
53 540 : auto Ax = DataContainer<data_t>(range);
54 540 : Ax = 0;
55 :
56 540 : WHEN("projecting forward and only the center voxel is set to 1")
57 540 : {
58 540 : op.apply(x, Ax);
59 :
60 540 : const auto weight = op.weight(0);
61 540 : CAPTURE(weight);
62 :
63 540 : THEN("The detector value is the weight for distance 0")
64 540 : {
65 540 : CAPTURE(Ax[0]);
66 540 : CHECK_EQ(weight, Approx(Ax[0]).epsilon(0.005));
67 540 : }
68 540 : }
69 540 : }
70 48600 : }
71 1620 : }
72 540 : }
73 :
74 : TEST_CASE_TEMPLATE("BlobVoxelProjector: Testing simple volume 2D with one detector Pixel", data_t,
75 : float, double)
76 540 : {
77 540 : const IndexVector_t sizeDomain({{5, 5}});
78 540 : const IndexVector_t sizeRange({{1, 1}});
79 :
80 540 : auto stc = SourceToCenterOfRotation{100};
81 540 : auto ctr = CenterOfRotationToDetector{5};
82 :
83 2160 : for (int i = 1; i < 4; i++) {
84 1620 : real_t scaling = i / 2.0f;
85 1620 : const RealVector_t spacing{{scaling, scaling}};
86 1620 : GIVEN("Spacing of " + std::to_string(scaling))
87 49140 : for (int i = 0; i < 360; i += 4) {
88 48600 : GIVEN("Ray of angle " + std::to_string(i))
89 48600 : {
90 540 : auto domain = VolumeDescriptor(sizeDomain, spacing);
91 540 : auto x = DataContainer<data_t>(domain);
92 540 : x = 0;
93 540 : std::vector<Geometry> geom;
94 540 : geom.emplace_back(stc, ctr, Degree{static_cast<real_t>(i)},
95 540 : VolumeData2D{Size2D{sizeDomain}, Spacing2D{spacing}},
96 540 : SinogramData2D{Size2D{sizeRange}});
97 540 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
98 540 : auto op = BlobVoxelProjector<data_t>(domain, range);
99 :
100 540 : auto Ax = DataContainer<data_t>(range);
101 540 : Ax = 0;
102 :
103 540 : WHEN("projecting forward and only the center voxel is set to 1")
104 540 : {
105 : // set center voxel to 1
106 540 : x(2, 2) = 1;
107 :
108 540 : op.apply(x, Ax);
109 :
110 540 : const auto weight = op.weight(0);
111 540 : CAPTURE(weight);
112 :
113 540 : THEN("The detector value is the weight for distance 0")
114 540 : {
115 540 : CAPTURE(Ax[0]);
116 540 : CHECK_EQ(weight, Approx(Ax[0]));
117 540 : }
118 540 : }
119 540 : }
120 48600 : }
121 1620 : }
122 540 : }
123 :
124 : TEST_CASE_TEMPLATE("BlobVoxelProjector: Testing simple volume 2D with two detector pixels", data_t,
125 : float, double)
126 180 : {
127 180 : const IndexVector_t sizeDomain({{5, 5}});
128 180 : const IndexVector_t sizeRange({{2, 1}});
129 :
130 180 : auto domain = VolumeDescriptor(sizeDomain);
131 180 : auto x = DataContainer<data_t>(domain);
132 180 : x = 0;
133 :
134 180 : auto stc = SourceToCenterOfRotation{100};
135 180 : auto ctr = CenterOfRotationToDetector{5};
136 180 : auto volData = VolumeData2D{Size2D{sizeDomain}};
137 180 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
138 :
139 16380 : for (int i = 0; i < 360; i += 4) {
140 16200 : GIVEN("Ray of angle " + std::to_string(i))
141 16200 : {
142 180 : std::vector<Geometry> geom;
143 180 : geom.emplace_back(stc, ctr, Degree{static_cast<real_t>(i)},
144 180 : VolumeData2D{Size2D{sizeDomain}}, SinogramData2D{Size2D{sizeRange}});
145 180 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
146 180 : auto op = BlobVoxelProjector<data_t>(domain, range);
147 :
148 180 : auto Ax = DataContainer<data_t>(range);
149 180 : Ax = 0;
150 :
151 180 : WHEN("only the center voxel is set to 1")
152 180 : {
153 : // set center voxel to 1
154 180 : x(2, 2) = 1;
155 :
156 180 : op.apply(x, Ax);
157 :
158 180 : const auto weight_dist_half = op.weight(0.4761878);
159 180 : CAPTURE(weight_dist_half);
160 :
161 180 : THEN("Detector values are symmetric")
162 180 : {
163 180 : CAPTURE(Ax[0]);
164 180 : CAPTURE(Ax[1]);
165 180 : CHECK_EQ(weight_dist_half, Approx(Ax[0]).epsilon(0.01));
166 180 : CHECK_EQ(weight_dist_half, Approx(Ax[1]).epsilon(0.01));
167 180 : }
168 180 : }
169 180 : }
170 16200 : }
171 180 : }
172 :
173 : TEST_CASE_TEMPLATE("BlobVoxelProjector: Testing simple volume 2D with three detector pixels",
174 : data_t, float, double)
175 288 : {
176 288 : const IndexVector_t sizeDomain({{5, 5}});
177 288 : const IndexVector_t sizeRange({{3, 1}});
178 :
179 288 : auto domain = VolumeDescriptor(sizeDomain);
180 288 : auto x = DataContainer<data_t>(domain);
181 288 : x = 0;
182 :
183 288 : auto stc = SourceToCenterOfRotation{100};
184 288 : auto ctr = CenterOfRotationToDetector{5};
185 288 : auto volData = VolumeData2D{Size2D{sizeDomain}};
186 288 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
187 :
188 21024 : for (int i = 0; i < 360; i += 5) {
189 20736 : GIVEN("Ray of angle " + std::to_string(i))
190 20736 : {
191 288 : std::vector<Geometry> geom;
192 288 : geom.emplace_back(stc, ctr, Degree{static_cast<real_t>(i)},
193 288 : VolumeData2D{Size2D{sizeDomain}}, SinogramData2D{Size2D{sizeRange}});
194 288 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
195 288 : auto op = BlobVoxelProjector<data_t>(domain, range);
196 :
197 288 : auto Ax = DataContainer<data_t>(range);
198 288 : Ax = 0;
199 :
200 288 : WHEN("only the center voxel is set to 1")
201 288 : {
202 : // set center voxel to 1
203 288 : x(2, 2) = 1;
204 :
205 288 : op.apply(x, Ax);
206 :
207 288 : const auto weight_dist0 = op.weight(0);
208 288 : const auto weight_dist1 = op.weight(0.95233786);
209 288 : CAPTURE(weight_dist0);
210 288 : CAPTURE(weight_dist1);
211 :
212 288 : THEN("the center detector pixel is correct")
213 288 : {
214 144 : CAPTURE(Ax[1]);
215 144 : CHECK_EQ(weight_dist0, Approx(Ax[1]));
216 144 : }
217 :
218 288 : THEN("the outer detector pixels are the same")
219 288 : {
220 144 : CAPTURE(Ax[0]);
221 144 : CAPTURE(Ax[2]);
222 144 : CHECK_EQ(weight_dist1, Approx(Ax[0]).epsilon(0.01));
223 144 : CHECK_EQ(weight_dist1, Approx(Ax[2]).epsilon(0.01));
224 144 : }
225 288 : }
226 288 : }
227 20736 : }
228 288 : }
229 :
230 : TEST_CASE("BlobVoxelProjector: Test single detector pixel")
231 22 : {
232 22 : const IndexVector_t sizeDomain({{5, 5}});
233 22 : const IndexVector_t sizeRange({{1, 1}});
234 :
235 22 : auto domain = VolumeDescriptor(sizeDomain);
236 22 : auto x = DataContainer(domain);
237 22 : x = 0;
238 :
239 22 : auto stc = SourceToCenterOfRotation{100};
240 22 : auto ctr = CenterOfRotationToDetector{5};
241 22 : auto volData = VolumeData2D{Size2D{sizeDomain}};
242 22 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
243 :
244 22 : GIVEN("a single detector of size 1, at 0 degree")
245 22 : {
246 3 : std::vector<Geometry> geom;
247 3 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{sizeDomain}},
248 3 : SinogramData2D{Size2D{sizeRange}});
249 :
250 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
251 3 : auto op = BlobVoxelProjector(domain, range);
252 :
253 3 : auto Ax = DataContainer(range);
254 3 : Ax = 0;
255 :
256 3 : WHEN("Setting only the voxels directly on the ray to 1")
257 3 : {
258 : // set all voxels on the ray to 1
259 1 : x(2, 0) = 1;
260 1 : x(2, 1) = 1;
261 1 : x(2, 2) = 1;
262 1 : x(2, 3) = 1;
263 1 : x(2, 4) = 1;
264 :
265 1 : op.apply(x, Ax);
266 :
267 1 : const auto weight = op.weight(0);
268 1 : CAPTURE(weight);
269 :
270 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
271 1 : {
272 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
273 1 : }
274 1 : }
275 :
276 3 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
277 3 : {
278 : // set all voxels above the ray to 1, i.e the visited neighbours
279 1 : x(3, 0) = 1;
280 1 : x(3, 1) = 1;
281 1 : x(3, 2) = 1;
282 1 : x(3, 3) = 1;
283 1 : x(3, 4) = 1;
284 :
285 1 : op.apply(x, Ax);
286 :
287 1 : THEN("The detector value is equal to the weight at the scaled distances")
288 1 : {
289 1 : real_t total_weight = 0;
290 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
291 5 : RealVector_t voxCoord{{3.f, slice}};
292 5 : voxCoord = voxCoord.array() + 0.5f;
293 5 : auto [detectorCoordShifted, scaling] =
294 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
295 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
296 5 : total_weight += op.weight(detectorCoord / scaling);
297 5 : }
298 :
299 1 : CHECK_EQ(Ax[0], Approx(total_weight));
300 1 : }
301 1 : }
302 :
303 3 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
304 3 : {
305 : // set all voxels above the ray to 1, i.e the visited neighbours
306 1 : x(1, 0) = 1;
307 1 : x(1, 1) = 1;
308 1 : x(1, 2) = 1;
309 1 : x(1, 3) = 1;
310 1 : x(1, 4) = 1;
311 :
312 1 : op.apply(x, Ax);
313 :
314 1 : THEN("The detector value is equal to the weight at the scaled distances")
315 1 : {
316 1 : real_t total_weight = 0;
317 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
318 5 : RealVector_t voxCoord{{1.f, slice}};
319 5 : voxCoord = voxCoord.array() + 0.5f;
320 5 : auto [detectorCoordShifted, scaling] =
321 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
322 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
323 5 : total_weight += op.weight(detectorCoord / scaling);
324 5 : }
325 :
326 1 : CHECK_EQ(Ax[0], Approx(total_weight));
327 1 : }
328 1 : }
329 3 : }
330 :
331 22 : GIVEN("a single detector of size 1, at 45 degree")
332 22 : {
333 5 : std::vector<Geometry> geom;
334 5 : geom.emplace_back(stc, ctr, Degree{45}, VolumeData2D{Size2D{sizeDomain}},
335 5 : SinogramData2D{Size2D{sizeRange}});
336 :
337 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
338 5 : auto op = BlobVoxelProjector(domain, range);
339 :
340 5 : auto Ax = DataContainer(range);
341 5 : Ax = 0;
342 :
343 5 : WHEN("Setting only the voxels directly on the ray to 1")
344 5 : {
345 : // set all voxels on the ray to 1
346 1 : x(0, 0) = 1;
347 1 : x(1, 1) = 1;
348 1 : x(2, 2) = 1;
349 1 : x(3, 3) = 1;
350 1 : x(4, 4) = 1;
351 :
352 1 : op.apply(x, Ax);
353 :
354 1 : const auto weight = op.weight(0);
355 1 : CAPTURE(weight);
356 :
357 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
358 1 : {
359 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
360 1 : }
361 1 : }
362 :
363 5 : WHEN("Setting only the voxels directly above the ray to 1")
364 5 : {
365 : // set all voxels on the ray to 1
366 1 : x(0, 1) = 1;
367 1 : x(1, 2) = 1;
368 1 : x(2, 3) = 1;
369 1 : x(3, 4) = 1;
370 :
371 1 : op.apply(x, Ax);
372 :
373 1 : THEN("The detector value is equal to the weight at the scaled distances")
374 1 : {
375 1 : real_t total_weight = 0;
376 4 : for (real_t slice : {0, 1, 2, 3}) {
377 4 : RealVector_t voxCoord{{slice, slice + 1}};
378 4 : voxCoord = voxCoord.array() + 0.5f;
379 4 : auto [detectorCoordShifted, scaling] =
380 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
381 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
382 4 : total_weight += op.weight(detectorCoord / scaling);
383 4 : }
384 :
385 1 : CHECK_EQ(Ax[0], Approx(total_weight));
386 1 : }
387 1 : }
388 :
389 5 : WHEN("Setting only the voxels directly below the ray to 1")
390 5 : {
391 : // set all voxels on the ray to 1
392 1 : x(1, 0) = 1;
393 1 : x(2, 1) = 1;
394 1 : x(3, 2) = 1;
395 1 : x(4, 3) = 1;
396 :
397 1 : op.apply(x, Ax);
398 :
399 1 : THEN("The detector value is equal to the weight at the scaled distances")
400 1 : {
401 1 : real_t total_weight = 0;
402 4 : for (real_t slice : {0, 1, 2, 3}) {
403 4 : RealVector_t voxCoord{{slice + 1, slice}};
404 4 : voxCoord = voxCoord.array() + 0.5f;
405 4 : auto [detectorCoordShifted, scaling] =
406 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
407 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
408 4 : total_weight += op.weight(detectorCoord / scaling);
409 4 : }
410 :
411 1 : CHECK_EQ(Ax[0], Approx(total_weight));
412 1 : }
413 1 : }
414 :
415 5 : WHEN("Setting only the voxels two steps above the ray")
416 5 : {
417 : // set all voxels on the ray to 1
418 1 : x(0, 2) = 1;
419 1 : x(1, 3) = 1;
420 1 : x(2, 4) = 1;
421 :
422 1 : op.apply(x, Ax);
423 :
424 1 : THEN("The detector value is equal to the weight at the scaled distances")
425 1 : {
426 1 : real_t total_weight = 0;
427 3 : for (real_t slice : {0, 1, 2}) {
428 3 : RealVector_t voxCoord{{slice, slice + 2}};
429 3 : voxCoord = voxCoord.array() + 0.5f;
430 3 : auto [detectorCoordShifted, scaling] =
431 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
432 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
433 3 : total_weight += op.weight(detectorCoord / scaling);
434 3 : }
435 :
436 1 : CHECK_EQ(Ax[0], Approx(total_weight));
437 1 : }
438 1 : }
439 :
440 5 : WHEN("Setting only the voxels two steps below the ray")
441 5 : {
442 : // set all voxels on the ray to 1
443 1 : x(2, 0) = 1;
444 1 : x(3, 1) = 1;
445 1 : x(4, 2) = 1;
446 :
447 1 : op.apply(x, Ax);
448 :
449 1 : THEN("The detector value is equal to the weight at the scaled distances")
450 1 : {
451 1 : real_t total_weight = 0;
452 3 : for (real_t slice : {0, 1, 2}) {
453 3 : RealVector_t voxCoord{{slice + 2, slice}};
454 3 : voxCoord = voxCoord.array() + 0.5f;
455 3 : auto [detectorCoordShifted, scaling] =
456 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
457 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
458 3 : total_weight += op.weight(detectorCoord / scaling);
459 3 : }
460 :
461 1 : CHECK_EQ(Ax[0], Approx(total_weight));
462 1 : }
463 1 : }
464 5 : }
465 :
466 22 : GIVEN("a single detector of size 1, at 90 degree")
467 22 : {
468 3 : std::vector<Geometry> geom;
469 3 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{sizeDomain}},
470 3 : SinogramData2D{Size2D{sizeRange}});
471 :
472 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
473 3 : auto op = BlobVoxelProjector(domain, range);
474 :
475 3 : auto Ax = DataContainer(range);
476 3 : Ax = 0;
477 :
478 3 : WHEN("Applying the BlobVoxelProjector to the volume")
479 3 : {
480 1 : x = 0;
481 : // set all voxels on the ray to 1
482 1 : x(0, 2) = 1;
483 1 : x(1, 2) = 1;
484 1 : x(2, 2) = 1;
485 1 : x(3, 2) = 1;
486 1 : x(4, 2) = 1;
487 :
488 1 : op.apply(x, Ax);
489 :
490 1 : const auto weight = op.weight(0);
491 1 : CAPTURE(weight);
492 :
493 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
494 1 : {
495 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
496 1 : }
497 1 : }
498 :
499 3 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
500 3 : {
501 : // set all voxels above the ray to 1, i.e the visited neighbours
502 1 : x(0, 3) = 1;
503 1 : x(1, 3) = 1;
504 1 : x(2, 3) = 1;
505 1 : x(3, 3) = 1;
506 1 : x(4, 3) = 1;
507 :
508 1 : op.apply(x, Ax);
509 :
510 1 : THEN("The detector value is equal to the weight at the scaled distances")
511 1 : {
512 1 : real_t total_weight = 0;
513 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
514 5 : RealVector_t voxCoord{{slice, 3.f}};
515 5 : voxCoord = voxCoord.array() + 0.5f;
516 5 : auto [detectorCoordShifted, scaling] =
517 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
518 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
519 5 : total_weight += op.weight(detectorCoord / scaling);
520 5 : }
521 :
522 1 : CHECK_EQ(Ax[0], Approx(total_weight));
523 1 : }
524 1 : }
525 :
526 3 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
527 3 : {
528 : // set all voxels above the ray to 1, i.e the visited neighbours
529 1 : x(0, 1) = 1;
530 1 : x(1, 1) = 1;
531 1 : x(2, 1) = 1;
532 1 : x(3, 1) = 1;
533 1 : x(4, 1) = 1;
534 :
535 1 : op.apply(x, Ax);
536 :
537 1 : THEN("The detector value is equal to the weight at the scaled distances")
538 1 : {
539 1 : real_t total_weight = 0;
540 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
541 5 : RealVector_t voxCoord{{slice, 1.f}};
542 5 : voxCoord = voxCoord.array() + 0.5f;
543 5 : auto [detectorCoordShifted, scaling] =
544 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
545 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
546 5 : total_weight += op.weight(detectorCoord / scaling);
547 5 : }
548 :
549 1 : CHECK_EQ(Ax[0], Approx(total_weight));
550 1 : }
551 1 : }
552 3 : }
553 :
554 22 : GIVEN("a single detector of size 1, at 135 degree")
555 22 : {
556 5 : std::vector<Geometry> geom;
557 5 : geom.emplace_back(stc, ctr, Degree{135}, VolumeData2D{Size2D{sizeDomain}},
558 5 : SinogramData2D{Size2D{sizeRange}});
559 :
560 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
561 5 : auto op = BlobVoxelProjector(domain, range);
562 :
563 5 : auto Ax = DataContainer(range);
564 5 : Ax = 0;
565 :
566 5 : WHEN("Applying the BlobVoxelProjector to the volume")
567 5 : {
568 : // set all voxels on the ray to 1
569 1 : x(0, 4) = 1;
570 1 : x(1, 3) = 1;
571 1 : x(2, 2) = 1;
572 1 : x(3, 1) = 1;
573 1 : x(4, 0) = 1;
574 :
575 1 : op.apply(x, Ax);
576 :
577 1 : const auto weight = op.weight(0);
578 1 : CAPTURE(weight);
579 :
580 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
581 1 : {
582 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
583 1 : }
584 1 : }
585 :
586 5 : WHEN("Setting only the voxels directly above on the ray to 1")
587 5 : {
588 : // set all voxels on the ray to 1
589 1 : x(1, 4) = 1;
590 1 : x(2, 3) = 1;
591 1 : x(3, 2) = 1;
592 1 : x(4, 1) = 1;
593 :
594 1 : op.apply(x, Ax);
595 :
596 1 : THEN("The detector value is equal to the weight at the scaled distances")
597 1 : {
598 1 : real_t total_weight = 0;
599 4 : for (real_t slice : {1, 2, 3, 4}) {
600 4 : RealVector_t voxCoord{{slice, 5 - slice}};
601 4 : voxCoord = voxCoord.array() + 0.5f;
602 4 : auto [detectorCoordShifted, scaling] =
603 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
604 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
605 4 : total_weight += op.weight(detectorCoord / scaling);
606 4 : }
607 :
608 1 : CHECK_EQ(Ax[0], Approx(total_weight));
609 1 : }
610 1 : }
611 :
612 5 : WHEN("Setting only the voxels directly above on the ray to 1")
613 5 : {
614 : // set all voxels on the ray to 1
615 1 : x(0, 3) = 1;
616 1 : x(1, 2) = 1;
617 1 : x(2, 1) = 1;
618 1 : x(3, 0) = 1;
619 :
620 1 : op.apply(x, Ax);
621 :
622 1 : THEN("The detector value is equal to the weight at the scaled distances")
623 1 : {
624 1 : real_t total_weight = 0;
625 4 : for (real_t slice : {0, 1, 2, 3}) {
626 4 : RealVector_t voxCoord{{slice, 3 - slice}};
627 4 : voxCoord = voxCoord.array() + 0.5f;
628 4 : auto [detectorCoordShifted, scaling] =
629 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
630 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
631 4 : total_weight += op.weight(detectorCoord / scaling);
632 4 : }
633 :
634 1 : CHECK_EQ(Ax[0], Approx(total_weight));
635 1 : }
636 1 : }
637 :
638 5 : WHEN("Setting only the voxels two above on the ray")
639 5 : {
640 : // set all voxels on the ray to 1
641 1 : x(2, 4) = 1;
642 1 : x(3, 3) = 1;
643 1 : x(4, 2) = 1;
644 :
645 1 : op.apply(x, Ax);
646 :
647 1 : THEN("The detector value is equal to the weight at the scaled distances")
648 1 : {
649 1 : real_t total_weight = 0;
650 3 : for (real_t slice : {2, 3, 4}) {
651 3 : RealVector_t voxCoord{{slice, 6 - slice}};
652 3 : voxCoord = voxCoord.array() + 0.5f;
653 3 : auto [detectorCoordShifted, scaling] =
654 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
655 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
656 3 : total_weight += op.weight(detectorCoord / scaling);
657 3 : }
658 :
659 1 : CHECK_EQ(Ax[0], Approx(total_weight));
660 1 : }
661 1 : }
662 :
663 5 : WHEN("Setting only the voxels directly above on the ray to 1")
664 5 : {
665 : // set all voxels on the ray to 1
666 1 : x(0, 2) = 1;
667 1 : x(1, 1) = 1;
668 1 : x(2, 0) = 1;
669 :
670 1 : op.apply(x, Ax);
671 :
672 1 : THEN("The detector value is equal to the weight at the scaled distances")
673 1 : {
674 1 : real_t total_weight = 0;
675 3 : for (real_t slice : {0, 1, 2}) {
676 3 : RealVector_t voxCoord{{slice, 2 - slice}};
677 3 : voxCoord = voxCoord.array() + 0.5f;
678 3 : auto [detectorCoordShifted, scaling] =
679 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
680 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
681 3 : total_weight += op.weight(detectorCoord / scaling);
682 3 : }
683 :
684 1 : CHECK_EQ(Ax[0], Approx(total_weight));
685 1 : }
686 1 : }
687 5 : }
688 :
689 22 : GIVEN("a single detector of size 1, at 180 degree")
690 22 : {
691 3 : std::vector<Geometry> geom;
692 3 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{sizeDomain}},
693 3 : SinogramData2D{Size2D{sizeRange}});
694 :
695 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
696 3 : auto op = BlobVoxelProjector(domain, range);
697 :
698 3 : auto Ax = DataContainer(range);
699 3 : Ax = 0;
700 :
701 3 : WHEN("Applying the BlobVoxelProjector to the volume")
702 3 : {
703 1 : x = 0;
704 : // set all voxels on the ray to 1
705 1 : x(2, 0) = 1;
706 1 : x(2, 1) = 1;
707 1 : x(2, 2) = 1;
708 1 : x(2, 3) = 1;
709 1 : x(2, 4) = 1;
710 :
711 1 : op.apply(x, Ax);
712 :
713 1 : const auto weight = op.weight(0);
714 1 : CAPTURE(weight);
715 :
716 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
717 1 : {
718 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
719 1 : }
720 1 : }
721 :
722 3 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
723 3 : {
724 : // set all voxels above the ray to 1, i.e the visited neighbours
725 1 : x(3, 0) = 1;
726 1 : x(3, 1) = 1;
727 1 : x(3, 2) = 1;
728 1 : x(3, 3) = 1;
729 1 : x(3, 4) = 1;
730 :
731 1 : op.apply(x, Ax);
732 :
733 1 : THEN("The detector value is equal to the weight at the scaled distances")
734 1 : {
735 1 : real_t total_weight = 0;
736 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
737 5 : RealVector_t voxCoord{{3.f, slice}};
738 5 : voxCoord = voxCoord.array() + 0.5f;
739 5 : auto [detectorCoordShifted, scaling] =
740 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
741 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
742 5 : total_weight += op.weight(detectorCoord / scaling);
743 5 : }
744 :
745 1 : CHECK_EQ(Ax[0], Approx(total_weight));
746 1 : }
747 1 : }
748 :
749 3 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
750 3 : {
751 : // set all voxels above the ray to 1, i.e the visited neighbours
752 1 : x(1, 0) = 1;
753 1 : x(1, 1) = 1;
754 1 : x(1, 2) = 1;
755 1 : x(1, 3) = 1;
756 1 : x(1, 4) = 1;
757 :
758 1 : op.apply(x, Ax);
759 :
760 1 : THEN("The detector value is equal to the weight at the scaled distances")
761 1 : {
762 1 : real_t total_weight = 0;
763 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
764 5 : RealVector_t voxCoord{{1.f, slice}};
765 5 : voxCoord = voxCoord.array() + 0.5f;
766 5 : auto [detectorCoordShifted, scaling] =
767 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
768 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
769 5 : total_weight += op.weight(detectorCoord / scaling);
770 5 : }
771 :
772 1 : CHECK_EQ(Ax[0], Approx(total_weight));
773 1 : }
774 1 : }
775 3 : }
776 :
777 22 : GIVEN("a single detector of size 1, at 225 degree")
778 22 : {
779 1 : std::vector<Geometry> geom;
780 1 : geom.emplace_back(stc, ctr, Degree{225}, VolumeData2D{Size2D{sizeDomain}},
781 1 : SinogramData2D{Size2D{sizeRange}});
782 :
783 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
784 1 : auto op = BlobVoxelProjector(domain, range);
785 :
786 1 : auto Ax = DataContainer(range);
787 1 : Ax = 0;
788 :
789 1 : WHEN("Applying the BlobVoxelProjector to the volume")
790 1 : {
791 : // set all voxels on the ray to 1
792 1 : x(0, 0) = 1;
793 1 : x(1, 1) = 1;
794 1 : x(2, 2) = 1;
795 1 : x(3, 3) = 1;
796 1 : x(4, 4) = 1;
797 :
798 1 : op.apply(x, Ax);
799 :
800 1 : const auto weight = op.weight(0);
801 1 : CAPTURE(weight);
802 :
803 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
804 1 : {
805 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
806 1 : }
807 1 : }
808 1 : }
809 :
810 22 : GIVEN("a single detector of size 1, at 270 degree")
811 22 : {
812 1 : std::vector<Geometry> geom;
813 1 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{sizeDomain}},
814 1 : SinogramData2D{Size2D{sizeRange}});
815 :
816 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
817 1 : auto op = BlobVoxelProjector(domain, range);
818 :
819 1 : auto Ax = DataContainer(range);
820 1 : Ax = 0;
821 :
822 1 : WHEN("Applying the BlobVoxelProjector to the volume")
823 1 : {
824 1 : x = 0;
825 : // set all voxels on the ray to 1
826 1 : x(0, 2) = 1;
827 1 : x(1, 2) = 1;
828 1 : x(2, 2) = 1;
829 1 : x(3, 2) = 1;
830 1 : x(4, 2) = 1;
831 :
832 1 : op.apply(x, Ax);
833 :
834 1 : const auto weight = op.weight(0);
835 1 : CAPTURE(weight);
836 :
837 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
838 1 : {
839 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
840 1 : }
841 1 : }
842 1 : }
843 :
844 22 : GIVEN("a single detector of size 1, at 315 degree")
845 22 : {
846 1 : std::vector<Geometry> geom;
847 1 : geom.emplace_back(stc, ctr, Degree{315}, VolumeData2D{Size2D{sizeDomain}},
848 1 : SinogramData2D{Size2D{sizeRange}});
849 :
850 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
851 1 : auto op = BlobVoxelProjector(domain, range);
852 :
853 1 : auto Ax = DataContainer(range);
854 1 : Ax = 0;
855 :
856 1 : WHEN("Applying the BlobVoxelProjector to the volume")
857 1 : {
858 : // set all voxels on the ray to 1
859 1 : x(0, 4) = 1;
860 1 : x(1, 3) = 1;
861 1 : x(2, 2) = 1;
862 1 : x(3, 1) = 1;
863 1 : x(4, 0) = 1;
864 :
865 1 : op.apply(x, Ax);
866 :
867 1 : const auto weight = op.weight(0);
868 1 : CAPTURE(weight);
869 :
870 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
871 1 : {
872 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
873 1 : }
874 1 : }
875 1 : }
876 22 : }
877 :
878 : TEST_CASE("BlobVoxelProjector: Test single detector pixel")
879 22 : {
880 22 : const IndexVector_t sizeDomain({{5, 5}});
881 22 : const IndexVector_t sizeRange({{1, 1}});
882 :
883 22 : auto domain = VolumeDescriptor(sizeDomain);
884 22 : auto x = DataContainer(domain);
885 22 : x = 0;
886 :
887 22 : auto stc = SourceToCenterOfRotation{100};
888 22 : auto ctr = CenterOfRotationToDetector{5};
889 22 : auto volData = VolumeData2D{Size2D{sizeDomain}};
890 22 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
891 :
892 22 : GIVEN("a single detector of size 1, at 0 degree")
893 22 : {
894 3 : std::vector<Geometry> geom;
895 3 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{sizeDomain}},
896 3 : SinogramData2D{Size2D{sizeRange}});
897 :
898 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
899 3 : auto op = BlobVoxelProjector(domain, range);
900 :
901 3 : auto Ax = DataContainer(range);
902 3 : Ax = 0;
903 :
904 3 : WHEN("Setting only the voxels directly on the ray to 1")
905 3 : {
906 : // set all voxels on the ray to 1
907 1 : x(2, 0) = 1;
908 1 : x(2, 1) = 1;
909 1 : x(2, 2) = 1;
910 1 : x(2, 3) = 1;
911 1 : x(2, 4) = 1;
912 :
913 1 : op.apply(x, Ax);
914 :
915 1 : const auto weight = op.weight(0);
916 1 : CAPTURE(weight);
917 :
918 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
919 1 : {
920 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
921 1 : }
922 1 : }
923 :
924 3 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
925 3 : {
926 : // set all voxels above the ray to 1, i.e the visited neighbours
927 1 : x(3, 0) = 1;
928 1 : x(3, 1) = 1;
929 1 : x(3, 2) = 1;
930 1 : x(3, 3) = 1;
931 1 : x(3, 4) = 1;
932 :
933 1 : op.apply(x, Ax);
934 :
935 1 : THEN("The detector value is equal to the weight at the scaled distances")
936 1 : {
937 1 : real_t total_weight = 0;
938 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
939 5 : RealVector_t voxCoord{{3.f, slice}};
940 5 : voxCoord = voxCoord.array() + 0.5f;
941 5 : auto [detectorCoordShifted, scaling] =
942 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
943 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
944 5 : total_weight += op.weight(detectorCoord / scaling);
945 5 : }
946 :
947 1 : CHECK_EQ(Ax[0], Approx(total_weight));
948 1 : }
949 1 : }
950 :
951 3 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
952 3 : {
953 : // set all voxels above the ray to 1, i.e the visited neighbours
954 1 : x(1, 0) = 1;
955 1 : x(1, 1) = 1;
956 1 : x(1, 2) = 1;
957 1 : x(1, 3) = 1;
958 1 : x(1, 4) = 1;
959 :
960 1 : op.apply(x, Ax);
961 :
962 1 : THEN("The detector value is equal to the weight at the scaled distances")
963 1 : {
964 1 : real_t total_weight = 0;
965 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
966 5 : RealVector_t voxCoord{{1.f, slice}};
967 5 : voxCoord = voxCoord.array() + 0.5f;
968 5 : auto [detectorCoordShifted, scaling] =
969 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
970 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
971 5 : total_weight += op.weight(detectorCoord / scaling);
972 5 : }
973 :
974 1 : CHECK_EQ(Ax[0], Approx(total_weight));
975 1 : }
976 1 : }
977 3 : }
978 :
979 22 : GIVEN("a single detector of size 1, at 45 degree")
980 22 : {
981 5 : std::vector<Geometry> geom;
982 5 : geom.emplace_back(stc, ctr, Degree{45}, VolumeData2D{Size2D{sizeDomain}},
983 5 : SinogramData2D{Size2D{sizeRange}});
984 :
985 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
986 5 : auto op = BlobVoxelProjector(domain, range);
987 :
988 5 : auto Ax = DataContainer(range);
989 5 : Ax = 0;
990 :
991 5 : WHEN("Setting only the voxels directly on the ray to 1")
992 5 : {
993 : // set all voxels on the ray to 1
994 1 : x(0, 0) = 1;
995 1 : x(1, 1) = 1;
996 1 : x(2, 2) = 1;
997 1 : x(3, 3) = 1;
998 1 : x(4, 4) = 1;
999 :
1000 1 : op.apply(x, Ax);
1001 :
1002 1 : const auto weight = op.weight(0);
1003 1 : CAPTURE(weight);
1004 :
1005 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1006 1 : {
1007 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1008 1 : }
1009 1 : }
1010 :
1011 5 : WHEN("Setting only the voxels directly above the ray to 1")
1012 5 : {
1013 : // set all voxels on the ray to 1
1014 1 : x(0, 1) = 1;
1015 1 : x(1, 2) = 1;
1016 1 : x(2, 3) = 1;
1017 1 : x(3, 4) = 1;
1018 :
1019 1 : op.apply(x, Ax);
1020 :
1021 1 : THEN("The detector value is equal to the weight at the scaled distances")
1022 1 : {
1023 1 : real_t total_weight = 0;
1024 4 : for (real_t slice : {0, 1, 2, 3}) {
1025 4 : RealVector_t voxCoord{{slice, slice + 1}};
1026 4 : voxCoord = voxCoord.array() + 0.5f;
1027 4 : auto [detectorCoordShifted, scaling] =
1028 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1029 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1030 4 : total_weight += op.weight(detectorCoord / scaling);
1031 4 : }
1032 :
1033 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1034 1 : }
1035 1 : }
1036 :
1037 5 : WHEN("Setting only the voxels directly below the ray to 1")
1038 5 : {
1039 : // set all voxels on the ray to 1
1040 1 : x(1, 0) = 1;
1041 1 : x(2, 1) = 1;
1042 1 : x(3, 2) = 1;
1043 1 : x(4, 3) = 1;
1044 :
1045 1 : op.apply(x, Ax);
1046 :
1047 1 : THEN("The detector value is equal to the weight at the scaled distances")
1048 1 : {
1049 1 : real_t total_weight = 0;
1050 4 : for (real_t slice : {0, 1, 2, 3}) {
1051 4 : RealVector_t voxCoord{{slice + 1, slice}};
1052 4 : voxCoord = voxCoord.array() + 0.5f;
1053 4 : auto [detectorCoordShifted, scaling] =
1054 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1055 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1056 4 : total_weight += op.weight(detectorCoord / scaling);
1057 4 : }
1058 :
1059 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1060 1 : }
1061 1 : }
1062 :
1063 5 : WHEN("Setting only the voxels two steps above the ray")
1064 5 : {
1065 : // set all voxels on the ray to 1
1066 1 : x(0, 2) = 1;
1067 1 : x(1, 3) = 1;
1068 1 : x(2, 4) = 1;
1069 :
1070 1 : op.apply(x, Ax);
1071 :
1072 1 : THEN("The detector value is equal to the weight at the scaled distances")
1073 1 : {
1074 1 : real_t total_weight = 0;
1075 3 : for (real_t slice : {0, 1, 2}) {
1076 3 : RealVector_t voxCoord{{slice, slice + 2}};
1077 3 : voxCoord = voxCoord.array() + 0.5f;
1078 3 : auto [detectorCoordShifted, scaling] =
1079 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1080 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1081 3 : total_weight += op.weight(detectorCoord / scaling);
1082 3 : }
1083 :
1084 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1085 1 : }
1086 1 : }
1087 :
1088 5 : WHEN("Setting only the voxels two steps below the ray")
1089 5 : {
1090 : // set all voxels on the ray to 1
1091 1 : x(2, 0) = 1;
1092 1 : x(3, 1) = 1;
1093 1 : x(4, 2) = 1;
1094 :
1095 1 : op.apply(x, Ax);
1096 :
1097 1 : THEN("The detector value is equal to the weight at the scaled distances")
1098 1 : {
1099 1 : real_t total_weight = 0;
1100 3 : for (real_t slice : {0, 1, 2}) {
1101 3 : RealVector_t voxCoord{{slice + 2, slice}};
1102 3 : voxCoord = voxCoord.array() + 0.5f;
1103 3 : auto [detectorCoordShifted, scaling] =
1104 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1105 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1106 3 : total_weight += op.weight(detectorCoord / scaling);
1107 3 : }
1108 :
1109 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1110 1 : }
1111 1 : }
1112 5 : }
1113 :
1114 22 : GIVEN("a single detector of size 1, at 90 degree")
1115 22 : {
1116 3 : std::vector<Geometry> geom;
1117 3 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{sizeDomain}},
1118 3 : SinogramData2D{Size2D{sizeRange}});
1119 :
1120 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1121 3 : auto op = BlobVoxelProjector(domain, range);
1122 :
1123 3 : auto Ax = DataContainer(range);
1124 3 : Ax = 0;
1125 :
1126 3 : WHEN("Applying the BlobVoxelProjector to the volume")
1127 3 : {
1128 1 : x = 0;
1129 : // set all voxels on the ray to 1
1130 1 : x(0, 2) = 1;
1131 1 : x(1, 2) = 1;
1132 1 : x(2, 2) = 1;
1133 1 : x(3, 2) = 1;
1134 1 : x(4, 2) = 1;
1135 :
1136 1 : op.apply(x, Ax);
1137 :
1138 1 : const auto weight = op.weight(0);
1139 1 : CAPTURE(weight);
1140 :
1141 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1142 1 : {
1143 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1144 1 : }
1145 1 : }
1146 :
1147 3 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
1148 3 : {
1149 : // set all voxels above the ray to 1, i.e the visited neighbours
1150 1 : x(0, 3) = 1;
1151 1 : x(1, 3) = 1;
1152 1 : x(2, 3) = 1;
1153 1 : x(3, 3) = 1;
1154 1 : x(4, 3) = 1;
1155 :
1156 1 : op.apply(x, Ax);
1157 :
1158 1 : THEN("The detector value is equal to the weight at the scaled distances")
1159 1 : {
1160 1 : real_t total_weight = 0;
1161 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
1162 5 : RealVector_t voxCoord{{slice, 3.f}};
1163 5 : voxCoord = voxCoord.array() + 0.5f;
1164 5 : auto [detectorCoordShifted, scaling] =
1165 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1166 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1167 5 : total_weight += op.weight(detectorCoord / scaling);
1168 5 : }
1169 :
1170 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1171 1 : }
1172 1 : }
1173 :
1174 3 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
1175 3 : {
1176 : // set all voxels above the ray to 1, i.e the visited neighbours
1177 1 : x(0, 1) = 1;
1178 1 : x(1, 1) = 1;
1179 1 : x(2, 1) = 1;
1180 1 : x(3, 1) = 1;
1181 1 : x(4, 1) = 1;
1182 :
1183 1 : op.apply(x, Ax);
1184 :
1185 1 : THEN("The detector value is equal to the weight at the scaled distances")
1186 1 : {
1187 1 : real_t total_weight = 0;
1188 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
1189 5 : RealVector_t voxCoord{{slice, 1.f}};
1190 5 : voxCoord = voxCoord.array() + 0.5f;
1191 5 : auto [detectorCoordShifted, scaling] =
1192 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1193 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1194 5 : total_weight += op.weight(detectorCoord / scaling);
1195 5 : }
1196 :
1197 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1198 1 : }
1199 1 : }
1200 3 : }
1201 :
1202 22 : GIVEN("a single detector of size 1, at 135 degree")
1203 22 : {
1204 5 : std::vector<Geometry> geom;
1205 5 : geom.emplace_back(stc, ctr, Degree{135}, VolumeData2D{Size2D{sizeDomain}},
1206 5 : SinogramData2D{Size2D{sizeRange}});
1207 :
1208 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1209 5 : auto op = BlobVoxelProjector(domain, range);
1210 :
1211 5 : auto Ax = DataContainer(range);
1212 5 : Ax = 0;
1213 :
1214 5 : WHEN("Applying the BlobVoxelProjector to the volume")
1215 5 : {
1216 : // set all voxels on the ray to 1
1217 1 : x(0, 4) = 1;
1218 1 : x(1, 3) = 1;
1219 1 : x(2, 2) = 1;
1220 1 : x(3, 1) = 1;
1221 1 : x(4, 0) = 1;
1222 :
1223 1 : op.apply(x, Ax);
1224 :
1225 1 : const auto weight = op.weight(0);
1226 1 : CAPTURE(weight);
1227 :
1228 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1229 1 : {
1230 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1231 1 : }
1232 1 : }
1233 :
1234 5 : WHEN("Setting only the voxels directly above on the ray to 1")
1235 5 : {
1236 : // set all voxels on the ray to 1
1237 1 : x(1, 4) = 1;
1238 1 : x(2, 3) = 1;
1239 1 : x(3, 2) = 1;
1240 1 : x(4, 1) = 1;
1241 :
1242 1 : op.apply(x, Ax);
1243 :
1244 1 : THEN("The detector value is equal to the weight at the scaled distances")
1245 1 : {
1246 1 : real_t total_weight = 0;
1247 4 : for (real_t slice : {1, 2, 3, 4}) {
1248 4 : RealVector_t voxCoord{{slice, 5 - slice}};
1249 4 : voxCoord = voxCoord.array() + 0.5f;
1250 4 : auto [detectorCoordShifted, scaling] =
1251 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1252 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1253 4 : total_weight += op.weight(detectorCoord / scaling);
1254 4 : }
1255 :
1256 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1257 1 : }
1258 1 : }
1259 :
1260 5 : WHEN("Setting only the voxels directly above on the ray to 1")
1261 5 : {
1262 : // set all voxels on the ray to 1
1263 1 : x(0, 3) = 1;
1264 1 : x(1, 2) = 1;
1265 1 : x(2, 1) = 1;
1266 1 : x(3, 0) = 1;
1267 :
1268 1 : op.apply(x, Ax);
1269 :
1270 1 : THEN("The detector value is equal to the weight at the scaled distances")
1271 1 : {
1272 1 : real_t total_weight = 0;
1273 4 : for (real_t slice : {0, 1, 2, 3}) {
1274 4 : RealVector_t voxCoord{{slice, 3 - slice}};
1275 4 : voxCoord = voxCoord.array() + 0.5f;
1276 4 : auto [detectorCoordShifted, scaling] =
1277 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1278 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1279 4 : total_weight += op.weight(detectorCoord / scaling);
1280 4 : }
1281 :
1282 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1283 1 : }
1284 1 : }
1285 :
1286 5 : WHEN("Setting only the voxels two above on the ray")
1287 5 : {
1288 : // set all voxels on the ray to 1
1289 1 : x(2, 4) = 1;
1290 1 : x(3, 3) = 1;
1291 1 : x(4, 2) = 1;
1292 :
1293 1 : op.apply(x, Ax);
1294 :
1295 1 : THEN("The detector value is equal to the weight at the scaled distances")
1296 1 : {
1297 1 : real_t total_weight = 0;
1298 3 : for (real_t slice : {2, 3, 4}) {
1299 3 : RealVector_t voxCoord{{slice, 6 - slice}};
1300 3 : voxCoord = voxCoord.array() + 0.5f;
1301 3 : auto [detectorCoordShifted, scaling] =
1302 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1303 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1304 3 : total_weight += op.weight(detectorCoord / scaling);
1305 3 : }
1306 :
1307 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1308 1 : }
1309 1 : }
1310 :
1311 5 : WHEN("Setting only the voxels directly above on the ray to 1")
1312 5 : {
1313 : // set all voxels on the ray to 1
1314 1 : x(0, 2) = 1;
1315 1 : x(1, 1) = 1;
1316 1 : x(2, 0) = 1;
1317 :
1318 1 : op.apply(x, Ax);
1319 :
1320 1 : THEN("The detector value is equal to the weight at the scaled distances")
1321 1 : {
1322 1 : real_t total_weight = 0;
1323 3 : for (real_t slice : {0, 1, 2}) {
1324 3 : RealVector_t voxCoord{{slice, 2 - slice}};
1325 3 : voxCoord = voxCoord.array() + 0.5f;
1326 3 : auto [detectorCoordShifted, scaling] =
1327 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1328 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1329 3 : total_weight += op.weight(detectorCoord / scaling);
1330 3 : }
1331 :
1332 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1333 1 : }
1334 1 : }
1335 5 : }
1336 :
1337 22 : GIVEN("a single detector of size 1, at 180 degree")
1338 22 : {
1339 3 : std::vector<Geometry> geom;
1340 3 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{sizeDomain}},
1341 3 : SinogramData2D{Size2D{sizeRange}});
1342 :
1343 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1344 3 : auto op = BlobVoxelProjector(domain, range);
1345 :
1346 3 : auto Ax = DataContainer(range);
1347 3 : Ax = 0;
1348 :
1349 3 : WHEN("Applying the BlobVoxelProjector to the volume")
1350 3 : {
1351 1 : x = 0;
1352 : // set all voxels on the ray to 1
1353 1 : x(2, 0) = 1;
1354 1 : x(2, 1) = 1;
1355 1 : x(2, 2) = 1;
1356 1 : x(2, 3) = 1;
1357 1 : x(2, 4) = 1;
1358 :
1359 1 : op.apply(x, Ax);
1360 :
1361 1 : const auto weight = op.weight(0);
1362 1 : CAPTURE(weight);
1363 :
1364 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1365 1 : {
1366 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1367 1 : }
1368 1 : }
1369 :
1370 3 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
1371 3 : {
1372 : // set all voxels above the ray to 1, i.e the visited neighbours
1373 1 : x(3, 0) = 1;
1374 1 : x(3, 1) = 1;
1375 1 : x(3, 2) = 1;
1376 1 : x(3, 3) = 1;
1377 1 : x(3, 4) = 1;
1378 :
1379 1 : op.apply(x, Ax);
1380 :
1381 1 : THEN("The detector value is equal to the weight at the scaled distances")
1382 1 : {
1383 1 : real_t total_weight = 0;
1384 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
1385 5 : RealVector_t voxCoord{{3.f, slice}};
1386 5 : voxCoord = voxCoord.array() + 0.5f;
1387 5 : auto [detectorCoordShifted, scaling] =
1388 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1389 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1390 5 : total_weight += op.weight(detectorCoord / scaling);
1391 5 : }
1392 :
1393 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1394 1 : }
1395 1 : }
1396 :
1397 3 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
1398 3 : {
1399 : // set all voxels above the ray to 1, i.e the visited neighbours
1400 1 : x(1, 0) = 1;
1401 1 : x(1, 1) = 1;
1402 1 : x(1, 2) = 1;
1403 1 : x(1, 3) = 1;
1404 1 : x(1, 4) = 1;
1405 :
1406 1 : op.apply(x, Ax);
1407 :
1408 1 : THEN("The detector value is equal to the weight at the scaled distances")
1409 1 : {
1410 1 : real_t total_weight = 0;
1411 5 : for (real_t slice : {0, 1, 2, 3, 4}) {
1412 5 : RealVector_t voxCoord{{1.f, slice}};
1413 5 : voxCoord = voxCoord.array() + 0.5f;
1414 5 : auto [detectorCoordShifted, scaling] =
1415 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1416 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1417 5 : total_weight += op.weight(detectorCoord / scaling);
1418 5 : }
1419 :
1420 1 : CHECK_EQ(Ax[0], Approx(total_weight));
1421 1 : }
1422 1 : }
1423 3 : }
1424 :
1425 22 : GIVEN("a single detector of size 1, at 225 degree")
1426 22 : {
1427 1 : std::vector<Geometry> geom;
1428 1 : geom.emplace_back(stc, ctr, Degree{225}, VolumeData2D{Size2D{sizeDomain}},
1429 1 : SinogramData2D{Size2D{sizeRange}});
1430 :
1431 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1432 1 : auto op = BlobVoxelProjector(domain, range);
1433 :
1434 1 : auto Ax = DataContainer(range);
1435 1 : Ax = 0;
1436 :
1437 1 : WHEN("Applying the BlobVoxelProjector to the volume")
1438 1 : {
1439 : // set all voxels on the ray to 1
1440 1 : x(0, 0) = 1;
1441 1 : x(1, 1) = 1;
1442 1 : x(2, 2) = 1;
1443 1 : x(3, 3) = 1;
1444 1 : x(4, 4) = 1;
1445 :
1446 1 : op.apply(x, Ax);
1447 :
1448 1 : const auto weight = op.weight(0);
1449 1 : CAPTURE(weight);
1450 :
1451 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1452 1 : {
1453 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1454 1 : }
1455 1 : }
1456 1 : }
1457 :
1458 22 : GIVEN("a single detector of size 1, at 270 degree")
1459 22 : {
1460 1 : std::vector<Geometry> geom;
1461 1 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{sizeDomain}},
1462 1 : SinogramData2D{Size2D{sizeRange}});
1463 :
1464 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1465 1 : auto op = BlobVoxelProjector(domain, range);
1466 :
1467 1 : auto Ax = DataContainer(range);
1468 1 : Ax = 0;
1469 :
1470 1 : WHEN("Applying the BlobVoxelProjector to the volume")
1471 1 : {
1472 1 : x = 0;
1473 : // set all voxels on the ray to 1
1474 1 : x(0, 2) = 1;
1475 1 : x(1, 2) = 1;
1476 1 : x(2, 2) = 1;
1477 1 : x(3, 2) = 1;
1478 1 : x(4, 2) = 1;
1479 :
1480 1 : op.apply(x, Ax);
1481 :
1482 1 : const auto weight = op.weight(0);
1483 1 : CAPTURE(weight);
1484 :
1485 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1486 1 : {
1487 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1488 1 : }
1489 1 : }
1490 1 : }
1491 :
1492 22 : GIVEN("a single detector of size 1, at 315 degree")
1493 22 : {
1494 1 : std::vector<Geometry> geom;
1495 1 : geom.emplace_back(stc, ctr, Degree{315}, VolumeData2D{Size2D{sizeDomain}},
1496 1 : SinogramData2D{Size2D{sizeRange}});
1497 :
1498 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1499 1 : auto op = BlobVoxelProjector(domain, range);
1500 :
1501 1 : auto Ax = DataContainer(range);
1502 1 : Ax = 0;
1503 :
1504 1 : WHEN("Applying the BlobVoxelProjector to the volume")
1505 1 : {
1506 : // set all voxels on the ray to 1
1507 1 : x(0, 4) = 1;
1508 1 : x(1, 3) = 1;
1509 1 : x(2, 2) = 1;
1510 1 : x(3, 1) = 1;
1511 1 : x(4, 0) = 1;
1512 :
1513 1 : op.apply(x, Ax);
1514 :
1515 1 : const auto weight = op.weight(0);
1516 1 : CAPTURE(weight);
1517 :
1518 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
1519 1 : {
1520 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
1521 1 : }
1522 1 : }
1523 1 : }
1524 22 : }
1525 :
1526 : TEST_CASE("BlobVoxelProjector: Test backward projection")
1527 17 : {
1528 17 : const IndexVector_t sizeDomain({{5, 5}});
1529 17 : const IndexVector_t sizeRange({{1, 1}});
1530 :
1531 17 : auto domain = VolumeDescriptor(sizeDomain);
1532 17 : auto x = DataContainer(domain);
1533 :
1534 17 : auto stc = SourceToCenterOfRotation{100};
1535 17 : auto ctr = CenterOfRotationToDetector{5};
1536 17 : auto volData = VolumeData2D{Size2D{sizeDomain}};
1537 17 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
1538 :
1539 17 : GIVEN("a single detector of size 1, at 0 degree")
1540 17 : {
1541 3 : std::vector<Geometry> geom;
1542 3 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{sizeDomain}},
1543 3 : SinogramData2D{Size2D{sizeRange}});
1544 :
1545 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1546 3 : auto op = BlobVoxelProjector(domain, range);
1547 :
1548 3 : auto Ax = DataContainer(range);
1549 3 : Ax = 0;
1550 3 : x = 0;
1551 :
1552 3 : WHEN("Backward projecting the operator")
1553 3 : {
1554 3 : Ax[0] = 1;
1555 :
1556 3 : op.applyAdjoint(Ax, x);
1557 :
1558 3 : THEN("The main direction of the ray is equal to the weight")
1559 3 : {
1560 1 : const auto weight = op.weight(0);
1561 1 : CAPTURE(weight);
1562 :
1563 1 : CHECK_EQ(x(2, 0), Approx(weight));
1564 1 : CHECK_EQ(x(2, 1), Approx(weight));
1565 1 : CHECK_EQ(x(2, 2), Approx(weight));
1566 1 : CHECK_EQ(x(2, 3), Approx(weight));
1567 1 : CHECK_EQ(x(2, 4), Approx(weight));
1568 1 : }
1569 :
1570 3 : THEN("The pixel 1 above the main direction of the ray have the projected weight")
1571 3 : {
1572 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1573 5 : RealVector_t voxCoord{{3.f, static_cast<float>(slice)}};
1574 5 : voxCoord = voxCoord.array() + 0.5f;
1575 5 : auto [detectorCoordShifted, scaling] =
1576 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1577 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1578 5 : auto weight = op.weight(detectorCoord / scaling);
1579 5 : CHECK_EQ(x(3, slice), Approx(weight));
1580 5 : }
1581 1 : }
1582 3 : THEN("The pixel 1 below the main direction of the ray have the projected weight")
1583 3 : {
1584 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1585 5 : RealVector_t voxCoord{{1.f, static_cast<float>(slice)}};
1586 5 : voxCoord = voxCoord.array() + 0.5f;
1587 5 : auto [detectorCoordShifted, scaling] =
1588 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1589 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1590 5 : auto weight = op.weight(detectorCoord / scaling);
1591 5 : CHECK_EQ(x(1, slice), Approx(weight));
1592 5 : }
1593 1 : }
1594 3 : }
1595 3 : }
1596 :
1597 17 : GIVEN("a single detector of size 1, at 180 degree")
1598 17 : {
1599 3 : std::vector<Geometry> geom;
1600 3 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{sizeDomain}},
1601 3 : SinogramData2D{Size2D{sizeRange}});
1602 :
1603 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1604 3 : auto op = BlobVoxelProjector(domain, range);
1605 :
1606 3 : auto Ax = DataContainer(range);
1607 3 : Ax = 0;
1608 3 : x = 0;
1609 :
1610 3 : WHEN("Backward projecting the operator")
1611 3 : {
1612 3 : Ax[0] = 1;
1613 :
1614 3 : op.applyAdjoint(Ax, x);
1615 :
1616 3 : THEN("The main direction of the ray is equal to the weight")
1617 3 : {
1618 1 : const auto weight = op.weight(0);
1619 1 : CAPTURE(weight);
1620 :
1621 1 : CHECK_EQ(x(2, 0), Approx(weight));
1622 1 : CHECK_EQ(x(2, 1), Approx(weight));
1623 1 : CHECK_EQ(x(2, 2), Approx(weight));
1624 1 : CHECK_EQ(x(2, 3), Approx(weight));
1625 1 : CHECK_EQ(x(2, 4), Approx(weight));
1626 1 : }
1627 :
1628 3 : THEN("The pixel 1 above the main direction of the ray have the projected weight")
1629 3 : {
1630 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1631 5 : RealVector_t voxCoord{{3.f, static_cast<float>(slice)}};
1632 5 : voxCoord = voxCoord.array() + 0.5f;
1633 5 : auto [detectorCoordShifted, scaling] =
1634 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1635 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1636 5 : auto weight = op.weight(detectorCoord / scaling);
1637 5 : CHECK_EQ(x(3, slice), Approx(weight));
1638 5 : }
1639 1 : }
1640 3 : THEN("The pixel 1 below the main direction of the ray have the projected weight")
1641 3 : {
1642 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1643 5 : RealVector_t voxCoord{{1.f, static_cast<float>(slice)}};
1644 5 : voxCoord = voxCoord.array() + 0.5f;
1645 5 : auto [detectorCoordShifted, scaling] =
1646 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1647 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1648 5 : auto weight = op.weight(detectorCoord / scaling);
1649 5 : CHECK_EQ(x(1, slice), Approx(weight));
1650 5 : }
1651 1 : }
1652 3 : }
1653 3 : }
1654 :
1655 17 : GIVEN("a single detector of size 1, at 90 degree")
1656 17 : {
1657 3 : std::vector<Geometry> geom;
1658 3 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{sizeDomain}},
1659 3 : SinogramData2D{Size2D{sizeRange}});
1660 :
1661 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1662 3 : auto op = BlobVoxelProjector(domain, range);
1663 :
1664 3 : auto Ax = DataContainer(range);
1665 3 : Ax = 0;
1666 3 : x = 0;
1667 :
1668 3 : WHEN("Backward projecting the operator")
1669 3 : {
1670 3 : Ax[0] = 1;
1671 :
1672 3 : op.applyAdjoint(Ax, x);
1673 :
1674 3 : THEN("The main direction of the ray is equal to the weight")
1675 3 : {
1676 1 : const auto weight = op.weight(0);
1677 1 : CAPTURE(weight);
1678 :
1679 1 : CHECK_EQ(x(0, 2), Approx(weight));
1680 1 : CHECK_EQ(x(1, 2), Approx(weight));
1681 1 : CHECK_EQ(x(2, 2), Approx(weight));
1682 1 : CHECK_EQ(x(3, 2), Approx(weight));
1683 1 : CHECK_EQ(x(4, 2), Approx(weight));
1684 1 : }
1685 :
1686 3 : THEN("The pixel 1 above the main direction of the ray have the projected weight")
1687 3 : {
1688 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1689 5 : RealVector_t voxCoord{{static_cast<float>(slice), 3.f}};
1690 5 : voxCoord = voxCoord.array() + 0.5f;
1691 5 : auto [detectorCoordShifted, scaling] =
1692 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1693 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1694 5 : auto weight = op.weight(detectorCoord / scaling);
1695 5 : CHECK_EQ(x(slice, 3), Approx(weight));
1696 5 : }
1697 1 : }
1698 3 : THEN("The pixel 1 below the main direction of the ray have the projected weight")
1699 3 : {
1700 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1701 5 : RealVector_t voxCoord{{static_cast<float>(slice), 1.f}};
1702 5 : voxCoord = voxCoord.array() + 0.5f;
1703 5 : auto [detectorCoordShifted, scaling] =
1704 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1705 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1706 5 : auto weight = op.weight(detectorCoord / scaling);
1707 5 : CHECK_EQ(x(slice, 1), Approx(weight));
1708 5 : }
1709 1 : }
1710 3 : }
1711 3 : }
1712 :
1713 17 : GIVEN("a single detector of size 1, at 270 degree")
1714 17 : {
1715 3 : std::vector<Geometry> geom;
1716 3 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{sizeDomain}},
1717 3 : SinogramData2D{Size2D{sizeRange}});
1718 :
1719 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1720 3 : auto op = BlobVoxelProjector(domain, range);
1721 :
1722 3 : auto Ax = DataContainer(range);
1723 3 : Ax = 0;
1724 3 : x = 0;
1725 :
1726 3 : WHEN("Backward projecting the operator")
1727 3 : {
1728 3 : Ax[0] = 1;
1729 :
1730 3 : op.applyAdjoint(Ax, x);
1731 :
1732 3 : THEN("The main direction of the ray is equal to the weight")
1733 3 : {
1734 1 : const auto weight = op.weight(0);
1735 1 : CAPTURE(weight);
1736 :
1737 1 : CHECK_EQ(x(0, 2), Approx(weight));
1738 1 : CHECK_EQ(x(1, 2), Approx(weight));
1739 1 : CHECK_EQ(x(2, 2), Approx(weight));
1740 1 : CHECK_EQ(x(3, 2), Approx(weight));
1741 1 : CHECK_EQ(x(4, 2), Approx(weight));
1742 1 : }
1743 3 : THEN("The pixel 1 above the main direction of the ray have the projected weight")
1744 3 : {
1745 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1746 5 : RealVector_t voxCoord{{static_cast<float>(slice), 3.f}};
1747 5 : voxCoord = voxCoord.array() + 0.5f;
1748 5 : auto [detectorCoordShifted, scaling] =
1749 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1750 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1751 5 : auto weight = op.weight(detectorCoord / scaling);
1752 5 : CHECK_EQ(x(slice, 3), Approx(weight));
1753 5 : }
1754 1 : }
1755 3 : THEN("The pixel 1 below the main direction of the ray have the projected weight")
1756 3 : {
1757 5 : for (index_t slice : {0, 1, 2, 3, 4}) {
1758 5 : RealVector_t voxCoord{{static_cast<float>(slice), 1.f}};
1759 5 : voxCoord = voxCoord.array() + 0.5f;
1760 5 : auto [detectorCoordShifted, scaling] =
1761 5 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1762 5 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1763 5 : auto weight = op.weight(detectorCoord / scaling);
1764 5 : CHECK_EQ(x(slice, 1), Approx(weight));
1765 5 : }
1766 1 : }
1767 3 : }
1768 3 : }
1769 :
1770 17 : GIVEN("a single detector of size 1, at 45 degree")
1771 17 : {
1772 5 : std::vector<Geometry> geom;
1773 5 : geom.emplace_back(stc, ctr, Degree{45}, VolumeData2D{Size2D{sizeDomain}},
1774 5 : SinogramData2D{Size2D{sizeRange}});
1775 :
1776 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1777 5 : auto op = BlobVoxelProjector(domain, range);
1778 :
1779 5 : auto Ax = DataContainer(range);
1780 5 : Ax = 0;
1781 5 : x = 0;
1782 :
1783 5 : WHEN("Backward projecting the operator")
1784 5 : {
1785 5 : Ax[0] = 1;
1786 :
1787 5 : op.applyAdjoint(Ax, x);
1788 :
1789 5 : THEN("The main direction of the ray is equal to the weight")
1790 5 : {
1791 1 : const auto weight = op.weight(0);
1792 1 : CAPTURE(weight);
1793 :
1794 1 : CHECK_EQ(x(0, 0), Approx(weight));
1795 1 : CHECK_EQ(x(1, 1), Approx(weight));
1796 1 : CHECK_EQ(x(2, 2), Approx(weight));
1797 1 : CHECK_EQ(x(3, 3), Approx(weight));
1798 1 : CHECK_EQ(x(4, 4), Approx(weight));
1799 1 : }
1800 :
1801 5 : THEN("The pixel 1 above the main direction of the ray have the projected weight")
1802 5 : {
1803 4 : for (index_t slice : {0, 1, 2, 3}) {
1804 4 : auto sliceF = static_cast<float>(slice);
1805 4 : RealVector_t voxCoord{{sliceF, sliceF + 1}};
1806 4 : voxCoord = voxCoord.array() + 0.5f;
1807 4 : auto [detectorCoordShifted, scaling] =
1808 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1809 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1810 4 : auto weight = op.weight(detectorCoord / scaling);
1811 4 : CHECK_EQ(x(slice, slice + 1), Approx(weight));
1812 4 : }
1813 1 : }
1814 :
1815 5 : THEN("The pixel 1 below the main direction of the ray have the projected weight")
1816 5 : {
1817 4 : for (index_t slice : {0, 1, 2, 3}) {
1818 4 : auto sliceF = static_cast<float>(slice);
1819 4 : RealVector_t voxCoord{{sliceF + 1, sliceF}};
1820 4 : voxCoord = voxCoord.array() + 0.5f;
1821 4 : auto [detectorCoordShifted, scaling] =
1822 4 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1823 4 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1824 4 : auto weight = op.weight(detectorCoord / scaling);
1825 4 : CHECK_EQ(x(slice + 1, slice), Approx(weight));
1826 4 : }
1827 1 : }
1828 :
1829 5 : THEN("The pixel 2 above the main direction of the ray have the projected weight")
1830 5 : {
1831 3 : for (index_t slice : {0, 1, 2}) {
1832 3 : auto sliceF = static_cast<float>(slice);
1833 3 : RealVector_t voxCoord{{sliceF, 2 - sliceF}};
1834 3 : voxCoord = voxCoord.array() + 0.5f;
1835 3 : auto [detectorCoordShifted, scaling] =
1836 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1837 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1838 3 : auto weight = op.weight(detectorCoord / scaling);
1839 3 : CHECK_EQ(x(slice, 2 - slice), Approx(weight));
1840 3 : }
1841 1 : }
1842 :
1843 5 : THEN("The pixel 2 above the main direction of the ray have the projected weight")
1844 5 : {
1845 3 : for (index_t slice : {0, 1, 2}) {
1846 3 : auto sliceF = static_cast<float>(slice);
1847 3 : RealVector_t voxCoord{{2 - sliceF, sliceF}};
1848 3 : voxCoord = voxCoord.array() + 0.5f;
1849 3 : auto [detectorCoordShifted, scaling] =
1850 3 : range.projectAndScaleVoxelOnDetector(voxCoord, 0);
1851 3 : real_t detectorCoord = detectorCoordShifted[0] - 0.5f;
1852 3 : auto weight = op.weight(detectorCoord / scaling);
1853 3 : CHECK_EQ(x(2 - slice, slice), Approx(weight));
1854 3 : }
1855 1 : }
1856 5 : }
1857 5 : }
1858 17 : }
1859 :
1860 : TEST_CASE_TEMPLATE("BlobVoxelProjector: Test weights", data_t, float, double)
1861 2 : {
1862 2 : const auto a = 2;
1863 2 : const auto alpha = 10.83;
1864 2 : const auto m = 2;
1865 :
1866 2 : std::array<double, 51> expected{
1867 2 : 1.3671064952680276, 1.3635202864368146, 1.3528128836429958,
1868 2 : 1.3351368521026497, 1.3107428112196733, 1.2799740558068384,
1869 2 : 1.2432592326454344, 1.2011032712842662, 1.1540768137691881,
1870 2 : 1.1028044260488254, 1.0479519029788988, 0.9902129983165086,
1871 2 : 0.930295920364307, 0.868909932853028, 0.80675238945173,
1872 2 : 0.7444965095220905, 0.6827801732549599, 0.6221959772930218,
1873 2 : 0.5632827487344612, 0.5065186675909519, 0.4523160970195944,
1874 2 : 0.40101816870068707, 0.35289711932154216, 0.30815432491147815,
1875 2 : 0.2669219342875695, 0.22926596246745282, 0.1951906707135796,
1876 2 : 0.1646440327638897, 0.13752406736650655, 0.11368580576665363,
1877 2 : 0.09294865929450916, 0.0751039563916854, 0.05992242974801804,
1878 2 : 0.04716145192475515, 0.036571840947646775, 0.027904084748497336,
1879 2 : 0.020913863801848218, 0.015366783581446854, 0.01104226128798172,
1880 2 : 0.007736543464620423, 0.005264861504239613, 0.003462759678911652,
1881 2 : 0.0021866543689359847, 0.0013137030013652602, 0.0007410763873099119,
1882 2 : 0.0003847384204545552, 0.0001778423521946198, 6.885294293780879e-05,
1883 2 : 1.9497773431607567e-05, 2.632583006403572e-06, 0};
1884 :
1885 2 : const IndexVector_t sizeDomain({{5, 5}});
1886 2 : const IndexVector_t sizeRange({{1, 1}});
1887 :
1888 2 : auto domain = VolumeDescriptor(sizeDomain);
1889 2 : auto x = DataContainer(domain);
1890 2 : x = 1;
1891 :
1892 2 : auto stc = SourceToCenterOfRotation{100};
1893 2 : auto ctr = CenterOfRotationToDetector{5};
1894 2 : auto volData = VolumeData2D{Size2D{sizeDomain}};
1895 2 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
1896 :
1897 2 : std::vector<Geometry> geom;
1898 2 : geom.emplace_back(stc, ctr, Degree{0}, std::move(volData), std::move(sinoData));
1899 :
1900 2 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1901 2 : auto op = BlobVoxelProjector(domain, range);
1902 :
1903 102 : for (int i = 0; i < 50; ++i) {
1904 100 : const auto distance = i / 25.;
1905 :
1906 100 : CAPTURE(i);
1907 100 : CAPTURE(distance);
1908 100 : CHECK_EQ(Approx(op.weight(distance)), expected[i]);
1909 100 : }
1910 2 : }
1911 :
1912 : TEST_SUITE_END();
|
