Line data Source code
1 : #include "doctest/doctest.h"
2 :
3 : #include "LutProjector.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(BlobProjector<float>);
20 :
21 : // Redefine GIVEN such that it's nicely usable inside an loop
22 : #undef GIVEN
23 69635 : #define GIVEN(...) DOCTEST_SUBCASE((std::string(" Given: ") + std::string(__VA_ARGS__)).c_str())
24 :
25 : TEST_CASE_TEMPLATE("BlobProjector: Testing rays going through the center of the volume", data_t,
26 : float, double)
27 2 : {
28 : // const IndexVector_t sizeDomain({{5, 5, 5}});
29 : // const IndexVector_t sizeRange({{1, 1, 1}});
30 2 : const IndexVector_t sizeDomain({{1, 1}});
31 2 : const IndexVector_t sizeRange({{1, 1}});
32 :
33 2 : auto domain = VolumeDescriptor(sizeDomain);
34 2 : auto x = DataContainer<data_t>(domain);
35 2 : x = 0;
36 :
37 2 : auto stc = SourceToCenterOfRotation{100};
38 2 : auto ctr = CenterOfRotationToDetector{0};
39 :
40 2 : const auto theta = 0;
41 2 : const auto thetad = Degree{theta}.to_radian();
42 2 : const auto beta = 0;
43 2 : const auto betad = Degree{beta}.to_radian();
44 2 : const auto alpha = 0;
45 2 : const auto alphad = Degree{alpha}.to_radian();
46 2 : std::vector<Geometry> geom;
47 : // geom.emplace_back(stc, ctr, VolumeData3D{Size3D{sizeDomain}},
48 : // SinogramData3D{Size3D{sizeRange}},
49 : // RotationAngles3D{Gamma{theta}, Beta{beta}, Alpha{alpha}});
50 2 : geom.emplace_back(stc, ctr, Degree{theta}, VolumeData2D{Size2D{sizeDomain}},
51 2 : SinogramData2D{Size2D{sizeRange}}, PrincipalPointOffset{0},
52 2 : RotationOffset2D{-0.5, -0.5});
53 2 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
54 :
55 2 : Geometry& g = geom[0];
56 :
57 2 : auto projInvMatrix = g.getInverseProjectionMatrix();
58 :
59 2 : MESSAGE(fmt::format("camera center: {}", g.getCameraCenter().format(vecfmt)));
60 2 : MESSAGE(fmt::format("Projection Matrix:\n{}", g.getProjectionMatrix().format(matfmt)));
61 :
62 2 : const RealVector_t principalpoint =
63 2 : (projInvMatrix * RealVector_t({{0.5, 1}})).head(2).normalized();
64 2 : MESSAGE(fmt::format("principal point: {}", principalpoint.format(vecfmt)));
65 :
66 2 : const auto projMatrix = g.getProjectionMatrix();
67 2 : const auto proj = RealMatrix_t({{std::cos(thetad), std::sin(thetad)}});
68 2 : MESSAGE(fmt::format("My Projection Matrix: {}", proj.format(vecfmt)));
69 :
70 : // const auto xaxis = RealMatrix_t({{1, 0, 0}});
71 : // const auto yaxis = RealMatrix_t({{0, 1, 0}});
72 : // const auto zaxis = RealMatrix_t({{0, 0, 1}});
73 : // const auto xaxish = RealMatrix_t({{1, 0, 0, 1}});
74 : // const auto yaxish = RealMatrix_t({{0, 1, 0, 1}});
75 : // const auto zaxish = RealMatrix_t({{0, 0, 1, 1}});
76 :
77 2 : const auto xaxis = RealVector_t({{1, 0}});
78 2 : const auto yaxis = RealVector_t({{0, 1}});
79 2 : const auto xaxish = RealVector_t({{1, 0, 1}});
80 2 : const auto yaxish = RealVector_t({{0, 1, 1}});
81 :
82 2 : const RealVector_t projx = (projMatrix.block(0, 0, 2, 2) * xaxis).head(1);
83 2 : const RealVector_t projy = (projMatrix.block(0, 0, 2, 2) * yaxis).head(1);
84 : // const RealVector_t projx = (projMatrix * xaxish).hnormalized();
85 : // const RealVector_t projy = (projMatrix * yaxish).hnormalized();
86 : // const RealVector_t projz = (projMatrix * zaxish);
87 :
88 2 : MESSAGE(fmt::format("project xaxis: {}", projx.format(vecfmt)));
89 2 : MESSAGE(fmt::format("project yaxis: {}", projy.format(vecfmt)));
90 :
91 2 : const RealVector_t myprojx = proj * xaxis;
92 2 : const RealVector_t myprojy = proj * yaxis;
93 2 : MESSAGE(fmt::format("project xaxis: {}", myprojx.format(vecfmt)));
94 2 : MESSAGE(fmt::format("project yaxis: {}", myprojy.format(vecfmt)));
95 2 : }
96 :
97 : TEST_CASE_TEMPLATE("BlobProjector: Testing rays going through the center of the volume", data_t,
98 : float, double)
99 180 : {
100 180 : const IndexVector_t sizeDomain({{5, 5, 5}});
101 180 : const IndexVector_t sizeRange({{1, 1, 1}});
102 :
103 180 : auto domain = VolumeDescriptor(sizeDomain);
104 180 : auto x = DataContainer<data_t>(domain);
105 180 : x = 0;
106 :
107 180 : auto stc = SourceToCenterOfRotation{100};
108 180 : auto ctr = CenterOfRotationToDetector{5};
109 :
110 : // set center voxel to 1
111 180 : x(2, 2, 2) = 1;
112 :
113 16380 : for (int i = 0; i < 360; i += 4) {
114 16200 : GIVEN("Ray of angle " + std::to_string(i))
115 16200 : {
116 180 : std::vector<Geometry> geom;
117 180 : geom.emplace_back(stc, ctr, VolumeData3D{Size3D{sizeDomain}},
118 180 : SinogramData3D{Size3D{sizeRange}},
119 180 : RotationAngles3D{Gamma{static_cast<real_t>(i)}});
120 180 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
121 180 : auto op = BlobProjector<data_t>(domain, range);
122 :
123 180 : auto Ax = DataContainer<data_t>(range);
124 180 : Ax = 0;
125 :
126 180 : WHEN("projecting forward and only the center voxel is set to 1")
127 180 : {
128 180 : op.apply(x, Ax);
129 :
130 180 : const auto weight = op.weight(0);
131 180 : CAPTURE(weight);
132 :
133 180 : THEN("The detector value is the weight for distance 0")
134 180 : {
135 180 : CAPTURE(Ax[0]);
136 180 : CHECK_EQ(weight, Approx(Ax[0]));
137 180 : }
138 180 : }
139 180 : }
140 16200 : }
141 180 : }
142 :
143 : TEST_CASE_TEMPLATE("BlobProjector: Testing rays going through the center of the volume", data_t,
144 : float, double)
145 180 : {
146 180 : const IndexVector_t sizeDomain({{5, 5}});
147 180 : const IndexVector_t sizeRange({{1, 1}});
148 :
149 180 : auto domain = VolumeDescriptor(sizeDomain);
150 180 : auto x = DataContainer<data_t>(domain);
151 180 : x = 0;
152 :
153 180 : auto stc = SourceToCenterOfRotation{100};
154 180 : auto ctr = CenterOfRotationToDetector{5};
155 180 : auto volData = VolumeData2D{Size2D{sizeDomain}};
156 180 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
157 :
158 16380 : for (int i = 0; i < 360; i += 4) {
159 16200 : GIVEN("Ray of angle " + std::to_string(i))
160 16200 : {
161 180 : std::vector<Geometry> geom;
162 180 : geom.emplace_back(stc, ctr, Degree{static_cast<real_t>(i)},
163 180 : VolumeData2D{Size2D{sizeDomain}}, SinogramData2D{Size2D{sizeRange}});
164 180 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
165 180 : auto op = BlobProjector<data_t>(domain, range);
166 :
167 180 : auto Ax = DataContainer<data_t>(range);
168 180 : Ax = 0;
169 :
170 180 : WHEN("projecting forward and only the center voxel is set to 1")
171 180 : {
172 : // set center voxel to 1
173 180 : x(2, 2) = 1;
174 :
175 180 : op.apply(x, Ax);
176 :
177 180 : const auto weight = op.weight(0);
178 180 : CAPTURE(weight);
179 :
180 180 : THEN("The detector value is the weight for distance 0")
181 180 : {
182 180 : CAPTURE(Ax[0]);
183 180 : CHECK_EQ(weight, Approx(Ax[0]));
184 180 : }
185 180 : }
186 180 : }
187 16200 : }
188 180 : }
189 :
190 : TEST_CASE_TEMPLATE("BlobProjector: Testing 2 rays going through the center of the volume", data_t,
191 : float, double)
192 180 : {
193 180 : const IndexVector_t sizeDomain({{5, 5}});
194 180 : const IndexVector_t sizeRange({{2, 1}});
195 :
196 180 : auto domain = VolumeDescriptor(sizeDomain);
197 180 : auto x = DataContainer<data_t>(domain);
198 180 : x = 0;
199 :
200 180 : auto stc = SourceToCenterOfRotation{100};
201 180 : auto ctr = CenterOfRotationToDetector{5};
202 180 : auto volData = VolumeData2D{Size2D{sizeDomain}};
203 180 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
204 :
205 16380 : for (int i = 0; i < 360; i += 4) {
206 16200 : GIVEN("Ray of angle " + std::to_string(i))
207 16200 : {
208 180 : std::vector<Geometry> geom;
209 180 : geom.emplace_back(stc, ctr, Degree{static_cast<real_t>(i)},
210 180 : VolumeData2D{Size2D{sizeDomain}}, SinogramData2D{Size2D{sizeRange}});
211 180 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
212 180 : auto op = BlobProjector<data_t>(domain, range);
213 :
214 180 : auto Ax = DataContainer<data_t>(range);
215 180 : Ax = 0;
216 :
217 180 : WHEN("only the center voxel is set to 1")
218 180 : {
219 : // set center voxel to 1
220 180 : x(2, 2) = 1;
221 :
222 180 : op.apply(x, Ax);
223 :
224 180 : const auto weight_dist_half = op.weight(0.4761878);
225 180 : CAPTURE(weight_dist_half);
226 :
227 180 : THEN("Detector values are symmetric")
228 180 : {
229 180 : CAPTURE(Ax[0]);
230 180 : CAPTURE(Ax[1]);
231 180 : CHECK_EQ(weight_dist_half, Approx(Ax[0]).epsilon(0.01));
232 180 : CHECK_EQ(weight_dist_half, Approx(Ax[1]).epsilon(0.01));
233 180 : }
234 180 : }
235 180 : }
236 16200 : }
237 180 : }
238 :
239 : TEST_CASE_TEMPLATE("BlobProjector: Testing 3 rays going through the center of the volume", data_t,
240 : float, double)
241 288 : {
242 288 : const IndexVector_t sizeDomain({{5, 5}});
243 288 : const IndexVector_t sizeRange({{3, 1}});
244 :
245 288 : auto domain = VolumeDescriptor(sizeDomain);
246 288 : auto x = DataContainer<data_t>(domain);
247 288 : x = 0;
248 :
249 288 : auto stc = SourceToCenterOfRotation{100};
250 288 : auto ctr = CenterOfRotationToDetector{5};
251 288 : auto volData = VolumeData2D{Size2D{sizeDomain}};
252 288 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
253 :
254 21024 : for (int i = 0; i < 360; i += 5) {
255 20736 : GIVEN("Ray of angle " + std::to_string(i))
256 20736 : {
257 288 : std::vector<Geometry> geom;
258 288 : geom.emplace_back(stc, ctr, Degree{static_cast<real_t>(i)},
259 288 : VolumeData2D{Size2D{sizeDomain}}, SinogramData2D{Size2D{sizeRange}});
260 288 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
261 288 : auto op = BlobProjector<data_t>(domain, range);
262 :
263 288 : auto Ax = DataContainer<data_t>(range);
264 288 : Ax = 0;
265 :
266 288 : WHEN("only the center voxel is set to 1")
267 288 : {
268 : // set center voxel to 1
269 288 : x(2, 2) = 1;
270 :
271 288 : op.apply(x, Ax);
272 :
273 288 : const auto weight_dist0 = op.weight(0);
274 288 : const auto weight_dist1 = op.weight(0.95233786);
275 288 : CAPTURE(weight_dist0);
276 288 : CAPTURE(weight_dist1);
277 :
278 288 : THEN("the center detector pixel is correct")
279 288 : {
280 144 : CAPTURE(Ax[1]);
281 144 : CHECK_EQ(weight_dist0, Approx(Ax[1]));
282 144 : }
283 :
284 288 : THEN("the outer detector pixels are the same")
285 288 : {
286 144 : CAPTURE(Ax[0]);
287 144 : CAPTURE(Ax[2]);
288 144 : CHECK_EQ(weight_dist1, Approx(Ax[0]).epsilon(0.01));
289 144 : CHECK_EQ(weight_dist1, Approx(Ax[2]).epsilon(0.01));
290 144 : }
291 288 : }
292 288 : }
293 20736 : }
294 288 : }
295 :
296 : TEST_CASE("BlobProjector: Test stepping of single rays")
297 28 : {
298 28 : const IndexVector_t sizeDomain({{5, 5}});
299 28 : const IndexVector_t sizeRange({{1, 1}});
300 :
301 28 : auto domain = VolumeDescriptor(sizeDomain);
302 28 : auto x = DataContainer(domain);
303 28 : x = 0;
304 :
305 28 : auto stc = SourceToCenterOfRotation{100};
306 28 : auto ctr = CenterOfRotationToDetector{5};
307 28 : auto volData = VolumeData2D{Size2D{sizeDomain}};
308 28 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
309 :
310 28 : GIVEN("a single detector of size 1, at 0 degree")
311 28 : {
312 5 : std::vector<Geometry> geom;
313 5 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{sizeDomain}},
314 5 : SinogramData2D{Size2D{sizeRange}});
315 :
316 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
317 5 : auto op = BlobProjector(domain, range);
318 :
319 5 : auto Ax = DataContainer(range);
320 5 : Ax = 0;
321 :
322 5 : WHEN("Setting only the voxels directly on the ray to 1")
323 5 : {
324 : // set all voxels on the ray to 1
325 1 : x(2, 0) = 1;
326 1 : x(2, 1) = 1;
327 1 : x(2, 2) = 1;
328 1 : x(2, 3) = 1;
329 1 : x(2, 4) = 1;
330 :
331 1 : op.apply(x, Ax);
332 :
333 1 : const auto weight = op.weight(0);
334 1 : CAPTURE(weight);
335 :
336 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
337 1 : {
338 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
339 1 : }
340 1 : }
341 :
342 5 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
343 5 : {
344 : // set all voxels above the ray to 1, i.e the visited neighbours
345 1 : x(3, 0) = 1;
346 1 : x(3, 1) = 1;
347 1 : x(3, 2) = 1;
348 1 : x(3, 3) = 1;
349 1 : x(3, 4) = 1;
350 :
351 1 : op.apply(x, Ax);
352 :
353 1 : const auto weight = op.weight(1);
354 1 : CAPTURE(weight);
355 :
356 1 : THEN("Each detector value is equal to 5 * the weight of distance 1")
357 1 : {
358 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
359 1 : }
360 1 : }
361 :
362 5 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
363 5 : {
364 : // set all voxels above the ray to 1, i.e the visited neighbours
365 1 : x(1, 0) = 1;
366 1 : x(1, 1) = 1;
367 1 : x(1, 2) = 1;
368 1 : x(1, 3) = 1;
369 1 : x(1, 4) = 1;
370 :
371 1 : op.apply(x, Ax);
372 :
373 1 : const auto weight = op.weight(1);
374 1 : CAPTURE(weight);
375 :
376 1 : THEN("Each detector value is equal to 5 * the weight of distance 1")
377 1 : {
378 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
379 1 : }
380 1 : }
381 :
382 5 : WHEN("Setting only the voxels 2 away neighbours above the ray to 1")
383 5 : {
384 : // set all voxels above the ray to 1, i.e the visited neighbours
385 1 : x(4, 0) = 1;
386 1 : x(4, 1) = 1;
387 1 : x(4, 2) = 1;
388 1 : x(4, 3) = 1;
389 1 : x(4, 4) = 1;
390 :
391 1 : op.apply(x, Ax);
392 :
393 1 : const auto weight = op.weight(2);
394 1 : CAPTURE(weight);
395 :
396 1 : THEN("Each detector value is equal to 5 * the weight of distance 2")
397 1 : {
398 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
399 1 : }
400 1 : }
401 :
402 5 : WHEN("Setting only the voxels 2 away neighbours below the ray to 1")
403 5 : {
404 : // set all voxels above the ray to 1, i.e the visited neighbours
405 1 : x(0, 0) = 1;
406 1 : x(0, 1) = 1;
407 1 : x(0, 2) = 1;
408 1 : x(0, 3) = 1;
409 1 : x(0, 4) = 1;
410 :
411 1 : op.apply(x, Ax);
412 :
413 1 : const auto weight = op.weight(2);
414 1 : CAPTURE(weight);
415 :
416 1 : THEN("Each detector value is equal to 5 * the weight of distance 2")
417 1 : {
418 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
419 1 : }
420 1 : }
421 5 : }
422 :
423 28 : GIVEN("a single detector of size 1, at 45 degree")
424 28 : {
425 5 : std::vector<Geometry> geom;
426 5 : geom.emplace_back(stc, ctr, Degree{45}, VolumeData2D{Size2D{sizeDomain}},
427 5 : SinogramData2D{Size2D{sizeRange}});
428 :
429 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
430 5 : auto op = BlobProjector(domain, range);
431 :
432 5 : auto Ax = DataContainer(range);
433 5 : Ax = 0;
434 :
435 5 : WHEN("Setting only the voxels directly on the ray to 1")
436 5 : {
437 : // set all voxels on the ray to 1
438 1 : x(0, 0) = 1;
439 1 : x(1, 1) = 1;
440 1 : x(2, 2) = 1;
441 1 : x(3, 3) = 1;
442 1 : x(4, 4) = 1;
443 :
444 1 : op.apply(x, Ax);
445 :
446 1 : const auto weight = op.weight(0);
447 1 : CAPTURE(weight);
448 :
449 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
450 1 : {
451 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
452 1 : }
453 1 : }
454 :
455 5 : WHEN("Setting only the voxels directly above on the ray to 1")
456 5 : {
457 : // set all voxels on the ray to 1
458 1 : x(0, 1) = 1;
459 1 : x(1, 2) = 1;
460 1 : x(2, 3) = 1;
461 1 : x(3, 4) = 1;
462 :
463 1 : op.apply(x, Ax);
464 :
465 1 : const auto weight = op.weight(0.707106769084930);
466 1 : CAPTURE(weight);
467 :
468 1 : THEN("Each detector value is equal to 4 * the weight of distance 0.7071")
469 1 : {
470 1 : CHECK_EQ(Ax[0], Approx(4 * weight).epsilon(0.01));
471 1 : }
472 1 : }
473 :
474 5 : WHEN("Setting only the voxels directly above on the ray to 1")
475 5 : {
476 : // set all voxels on the ray to 1
477 1 : x(1, 0) = 1;
478 1 : x(2, 1) = 1;
479 1 : x(3, 2) = 1;
480 1 : x(4, 3) = 1;
481 :
482 1 : op.apply(x, Ax);
483 :
484 1 : const auto weight = op.weight(0.707106769084930);
485 1 : CAPTURE(weight);
486 :
487 1 : THEN("Each detector value is equal to 5 * the weight of distance 0.7071")
488 1 : {
489 1 : CHECK_EQ(Ax[0], Approx(4 * weight).epsilon(0.01));
490 1 : }
491 1 : }
492 :
493 5 : WHEN("Setting only the voxels two above on the ray")
494 5 : {
495 : // set all voxels on the ray to 1
496 1 : x(0, 2) = 1;
497 1 : x(1, 3) = 1;
498 1 : x(2, 4) = 1;
499 :
500 1 : op.apply(x, Ax);
501 :
502 1 : const auto distance = 2 * 0.707106769084930;
503 1 : const auto weight = op.weight(distance);
504 1 : CAPTURE(distance);
505 1 : CAPTURE(weight);
506 :
507 1 : THEN("Each detector value is equal to 3 * the weight of distance 0.7071")
508 1 : {
509 1 : CHECK_EQ(Ax[0], Approx(3 * weight));
510 1 : }
511 1 : }
512 :
513 5 : WHEN("Setting only the voxels directly above on the ray to 1")
514 5 : {
515 : // set all voxels on the ray to 1
516 1 : x(2, 0) = 1;
517 1 : x(3, 1) = 1;
518 1 : x(4, 2) = 1;
519 :
520 1 : op.apply(x, Ax);
521 :
522 1 : const auto distance = 2 * 0.707106769084930;
523 1 : const auto weight = op.weight(distance);
524 1 : CAPTURE(distance);
525 1 : CAPTURE(weight);
526 :
527 1 : THEN("Each detector value is equal to 3 * the weight of distance 0.7071")
528 1 : {
529 1 : CHECK_EQ(Ax[0], Approx(3 * weight));
530 1 : }
531 1 : }
532 5 : }
533 :
534 28 : GIVEN("a single detector of size 1, at 90 degree")
535 28 : {
536 5 : std::vector<Geometry> geom;
537 5 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{sizeDomain}},
538 5 : SinogramData2D{Size2D{sizeRange}});
539 :
540 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
541 5 : auto op = BlobProjector(domain, range);
542 :
543 5 : auto Ax = DataContainer(range);
544 5 : Ax = 0;
545 :
546 5 : WHEN("Applying the BlobProjector to the volume")
547 5 : {
548 1 : x = 0;
549 : // set all voxels on the ray to 1
550 1 : x(0, 2) = 1;
551 1 : x(1, 2) = 1;
552 1 : x(2, 2) = 1;
553 1 : x(3, 2) = 1;
554 1 : x(4, 2) = 1;
555 :
556 1 : op.apply(x, Ax);
557 :
558 1 : const auto weight = op.weight(0);
559 1 : CAPTURE(weight);
560 :
561 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
562 1 : {
563 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
564 1 : }
565 1 : }
566 :
567 5 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
568 5 : {
569 : // set all voxels above the ray to 1, i.e the visited neighbours
570 1 : x(0, 3) = 1;
571 1 : x(1, 3) = 1;
572 1 : x(2, 3) = 1;
573 1 : x(3, 3) = 1;
574 1 : x(4, 3) = 1;
575 :
576 1 : op.apply(x, Ax);
577 :
578 1 : const auto weight = op.weight(1);
579 1 : CAPTURE(weight);
580 :
581 1 : THEN("Each detector value is equal to 5 * the weight of distance 1")
582 1 : {
583 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
584 1 : }
585 1 : }
586 :
587 5 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
588 5 : {
589 : // set all voxels above the ray to 1, i.e the visited neighbours
590 1 : x(0, 1) = 1;
591 1 : x(1, 1) = 1;
592 1 : x(2, 1) = 1;
593 1 : x(3, 1) = 1;
594 1 : x(4, 1) = 1;
595 :
596 1 : op.apply(x, Ax);
597 :
598 1 : const auto weight = op.weight(1);
599 1 : CAPTURE(weight);
600 :
601 1 : THEN("Each detector value is equal to 5 * the weight of distance 1")
602 1 : {
603 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
604 1 : }
605 1 : }
606 :
607 5 : WHEN("Setting only the voxels 2 away neighbours above the ray to 1")
608 5 : {
609 : // set all voxels above the ray to 1, i.e the visited neighbours
610 1 : x(0, 4) = 1;
611 1 : x(1, 4) = 1;
612 1 : x(2, 4) = 1;
613 1 : x(3, 4) = 1;
614 1 : x(4, 4) = 1;
615 :
616 1 : op.apply(x, Ax);
617 :
618 1 : const auto weight = op.weight(2);
619 1 : CAPTURE(weight);
620 :
621 1 : THEN("Each detector value is equal to 5 * the weight of distance 2")
622 1 : {
623 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
624 1 : }
625 1 : }
626 :
627 5 : WHEN("Setting only the voxels 2 away neighbours below the ray to 1")
628 5 : {
629 : // set all voxels above the ray to 1, i.e the visited neighbours
630 1 : x(0, 0) = 1;
631 1 : x(1, 0) = 1;
632 1 : x(2, 0) = 1;
633 1 : x(3, 0) = 1;
634 1 : x(4, 0) = 1;
635 :
636 1 : op.apply(x, Ax);
637 :
638 1 : const auto weight = op.weight(2);
639 1 : CAPTURE(weight);
640 :
641 1 : THEN("Each detector value is equal to 5 * the weight of distance 2")
642 1 : {
643 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
644 1 : }
645 1 : }
646 5 : }
647 :
648 28 : GIVEN("a single detector of size 1, at 135 degree")
649 28 : {
650 5 : std::vector<Geometry> geom;
651 5 : geom.emplace_back(stc, ctr, Degree{135}, VolumeData2D{Size2D{sizeDomain}},
652 5 : SinogramData2D{Size2D{sizeRange}});
653 :
654 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
655 5 : auto op = BlobProjector(domain, range);
656 :
657 5 : auto Ax = DataContainer(range);
658 5 : Ax = 0;
659 :
660 5 : WHEN("Applying the BlobProjector to the volume")
661 5 : {
662 : // set all voxels on the ray to 1
663 1 : x(0, 4) = 1;
664 1 : x(1, 3) = 1;
665 1 : x(2, 2) = 1;
666 1 : x(3, 1) = 1;
667 1 : x(4, 0) = 1;
668 :
669 1 : op.apply(x, Ax);
670 :
671 1 : const auto weight = op.weight(0);
672 1 : CAPTURE(weight);
673 :
674 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
675 1 : {
676 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
677 1 : }
678 1 : }
679 :
680 5 : WHEN("Setting only the voxels directly above on the ray to 1")
681 5 : {
682 : // set all voxels on the ray to 1
683 1 : x(1, 4) = 1;
684 1 : x(2, 3) = 1;
685 1 : x(3, 2) = 1;
686 1 : x(4, 1) = 1;
687 :
688 1 : op.apply(x, Ax);
689 :
690 1 : const auto distance = 0.707106769084930;
691 1 : const auto weight = op.weight(distance);
692 1 : CAPTURE(distance);
693 1 : CAPTURE(weight);
694 :
695 1 : THEN("Each detector value is equal to 4 * the weight of distance 0.7071")
696 1 : {
697 1 : CHECK_EQ(Ax[0], Approx(4 * weight).epsilon(0.01));
698 1 : }
699 1 : }
700 :
701 5 : WHEN("Setting only the voxels directly above on the ray to 1")
702 5 : {
703 : // set all voxels on the ray to 1
704 1 : x(0, 3) = 1;
705 1 : x(1, 2) = 1;
706 1 : x(2, 1) = 1;
707 1 : x(3, 0) = 1;
708 :
709 1 : op.apply(x, Ax);
710 :
711 1 : const auto distance = 0.707106769084930;
712 1 : const auto weight = op.weight(distance);
713 1 : CAPTURE(distance);
714 1 : CAPTURE(weight);
715 :
716 1 : THEN("Each detector value is equal to 5 * the weight of distance 0.7071")
717 1 : {
718 1 : CHECK_EQ(Ax[0], Approx(4 * weight).epsilon(0.01));
719 1 : }
720 1 : }
721 :
722 5 : WHEN("Setting only the voxels two above on the ray")
723 5 : {
724 : // set all voxels on the ray to 1
725 1 : x(2, 4) = 1;
726 1 : x(3, 3) = 1;
727 1 : x(4, 2) = 1;
728 :
729 1 : op.apply(x, Ax);
730 :
731 1 : const auto distance = 2 * 0.707106769084930;
732 1 : const auto weight = op.weight(distance);
733 1 : CAPTURE(distance);
734 1 : CAPTURE(weight);
735 :
736 1 : THEN("Each detector value is equal to 3 * the weight of distance 0.7071")
737 1 : {
738 1 : CHECK_EQ(Ax[0], Approx(3 * weight));
739 1 : }
740 1 : }
741 :
742 5 : WHEN("Setting only the voxels directly above on the ray to 1")
743 5 : {
744 : // set all voxels on the ray to 1
745 1 : x(0, 2) = 1;
746 1 : x(1, 1) = 1;
747 1 : x(2, 0) = 1;
748 :
749 1 : op.apply(x, Ax);
750 :
751 1 : const auto distance = 2 * 0.707106769084930;
752 1 : const auto weight = op.weight(distance);
753 1 : CAPTURE(distance);
754 1 : CAPTURE(weight);
755 :
756 1 : THEN("Each detector value is equal to 3 * the weight of distance 0.7071")
757 1 : {
758 1 : CHECK_EQ(Ax[0], Approx(3 * weight));
759 1 : }
760 1 : }
761 5 : }
762 :
763 28 : GIVEN("a single detector of size 1, at 180 degree")
764 28 : {
765 5 : std::vector<Geometry> geom;
766 5 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{sizeDomain}},
767 5 : SinogramData2D{Size2D{sizeRange}});
768 :
769 5 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
770 5 : auto op = BlobProjector(domain, range);
771 :
772 5 : auto Ax = DataContainer(range);
773 5 : Ax = 0;
774 :
775 5 : WHEN("Applying the BlobProjector to the volume")
776 5 : {
777 1 : x = 0;
778 : // set all voxels on the ray to 1
779 1 : x(2, 0) = 1;
780 1 : x(2, 1) = 1;
781 1 : x(2, 2) = 1;
782 1 : x(2, 3) = 1;
783 1 : x(2, 4) = 1;
784 :
785 1 : op.apply(x, Ax);
786 :
787 1 : const auto weight = op.weight(0);
788 1 : CAPTURE(weight);
789 :
790 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
791 1 : {
792 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
793 1 : }
794 1 : }
795 :
796 5 : WHEN("Setting only the voxels direct neighbours above the ray to 1")
797 5 : {
798 : // set all voxels above the ray to 1, i.e the visited neighbours
799 1 : x(3, 0) = 1;
800 1 : x(3, 1) = 1;
801 1 : x(3, 2) = 1;
802 1 : x(3, 3) = 1;
803 1 : x(3, 4) = 1;
804 :
805 1 : op.apply(x, Ax);
806 :
807 1 : const auto weight = op.weight(1);
808 1 : CAPTURE(weight);
809 :
810 1 : THEN("Each detector value is equal to 5 * the weight of distance 1")
811 1 : {
812 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
813 1 : }
814 1 : }
815 :
816 5 : WHEN("Setting only the voxels direct neighbours below the ray to 1")
817 5 : {
818 : // set all voxels above the ray to 1, i.e the visited neighbours
819 1 : x(1, 0) = 1;
820 1 : x(1, 1) = 1;
821 1 : x(1, 2) = 1;
822 1 : x(1, 3) = 1;
823 1 : x(1, 4) = 1;
824 :
825 1 : op.apply(x, Ax);
826 :
827 1 : const auto weight = op.weight(1);
828 1 : CAPTURE(weight);
829 :
830 1 : THEN("Each detector value is equal to 5 * the weight of distance 1")
831 1 : {
832 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
833 1 : }
834 1 : }
835 :
836 5 : WHEN("Setting only the voxels 2 away neighbours above the ray to 1")
837 5 : {
838 : // set all voxels above the ray to 1, i.e the visited neighbours
839 1 : x(4, 0) = 1;
840 1 : x(4, 1) = 1;
841 1 : x(4, 2) = 1;
842 1 : x(4, 3) = 1;
843 1 : x(4, 4) = 1;
844 :
845 1 : op.apply(x, Ax);
846 :
847 1 : const auto weight = op.weight(2);
848 1 : CAPTURE(weight);
849 :
850 1 : THEN("Each detector value is equal to 5 * the weight of distance 2")
851 1 : {
852 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
853 1 : }
854 1 : }
855 :
856 5 : WHEN("Setting only the voxels 2 away neighbours below the ray to 1")
857 5 : {
858 : // set all voxels above the ray to 1, i.e the visited neighbours
859 1 : x(0, 0) = 1;
860 1 : x(0, 1) = 1;
861 1 : x(0, 2) = 1;
862 1 : x(0, 3) = 1;
863 1 : x(0, 4) = 1;
864 :
865 1 : op.apply(x, Ax);
866 :
867 1 : const auto weight = op.weight(2);
868 1 : CAPTURE(weight);
869 :
870 1 : THEN("Each detector value is equal to 5 * the weight of distance 2")
871 1 : {
872 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
873 1 : }
874 1 : }
875 5 : }
876 :
877 28 : GIVEN("a single detector of size 1, at 225 degree")
878 28 : {
879 1 : std::vector<Geometry> geom;
880 1 : geom.emplace_back(stc, ctr, Degree{225}, VolumeData2D{Size2D{sizeDomain}},
881 1 : SinogramData2D{Size2D{sizeRange}});
882 :
883 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
884 1 : auto op = BlobProjector(domain, range);
885 :
886 1 : auto Ax = DataContainer(range);
887 1 : Ax = 0;
888 :
889 1 : WHEN("Applying the BlobProjector to the volume")
890 1 : {
891 : // set all voxels on the ray to 1
892 1 : x(0, 0) = 1;
893 1 : x(1, 1) = 1;
894 1 : x(2, 2) = 1;
895 1 : x(3, 3) = 1;
896 1 : x(4, 4) = 1;
897 :
898 1 : op.apply(x, Ax);
899 :
900 1 : const auto weight = op.weight(0);
901 1 : CAPTURE(weight);
902 :
903 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
904 1 : {
905 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
906 1 : }
907 1 : }
908 1 : }
909 :
910 28 : GIVEN("a single detector of size 1, at 270 degree")
911 28 : {
912 1 : std::vector<Geometry> geom;
913 1 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{sizeDomain}},
914 1 : SinogramData2D{Size2D{sizeRange}});
915 :
916 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
917 1 : auto op = BlobProjector(domain, range);
918 :
919 1 : auto Ax = DataContainer(range);
920 1 : Ax = 0;
921 :
922 1 : WHEN("Applying the BlobProjector to the volume")
923 1 : {
924 1 : x = 0;
925 : // set all voxels on the ray to 1
926 1 : x(0, 2) = 1;
927 1 : x(1, 2) = 1;
928 1 : x(2, 2) = 1;
929 1 : x(3, 2) = 1;
930 1 : x(4, 2) = 1;
931 :
932 1 : op.apply(x, Ax);
933 :
934 1 : const auto weight = op.weight(0);
935 1 : CAPTURE(weight);
936 :
937 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
938 1 : {
939 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
940 1 : }
941 1 : }
942 1 : }
943 :
944 28 : GIVEN("a single detector of size 1, at 315 degree")
945 28 : {
946 1 : std::vector<Geometry> geom;
947 1 : geom.emplace_back(stc, ctr, Degree{315}, VolumeData2D{Size2D{sizeDomain}},
948 1 : SinogramData2D{Size2D{sizeRange}});
949 :
950 1 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
951 1 : auto op = BlobProjector(domain, range);
952 :
953 1 : auto Ax = DataContainer(range);
954 1 : Ax = 0;
955 :
956 1 : WHEN("Applying the BlobProjector to the volume")
957 1 : {
958 : // set all voxels on the ray to 1
959 1 : x(0, 4) = 1;
960 1 : x(1, 3) = 1;
961 1 : x(2, 2) = 1;
962 1 : x(3, 1) = 1;
963 1 : x(4, 0) = 1;
964 :
965 1 : op.apply(x, Ax);
966 :
967 1 : const auto weight = op.weight(0);
968 1 : CAPTURE(weight);
969 :
970 1 : THEN("Each detector value is equal to 5 * the weight of distance 0")
971 1 : {
972 1 : CHECK_EQ(Ax[0], Approx(5 * weight));
973 1 : }
974 1 : }
975 1 : }
976 28 : }
977 :
978 : TEST_CASE("LutProjector: Test single rays backward projection")
979 15 : {
980 15 : const IndexVector_t sizeDomain({{5, 5}});
981 15 : const IndexVector_t sizeRange({{1, 1}});
982 :
983 15 : auto domain = VolumeDescriptor(sizeDomain);
984 15 : auto x = DataContainer(domain);
985 :
986 15 : auto stc = SourceToCenterOfRotation{100};
987 15 : auto ctr = CenterOfRotationToDetector{5};
988 15 : auto volData = VolumeData2D{Size2D{sizeDomain}};
989 15 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
990 :
991 15 : GIVEN("a single detector of size 1, at 0 degree")
992 15 : {
993 3 : std::vector<Geometry> geom;
994 3 : geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{sizeDomain}},
995 3 : SinogramData2D{Size2D{sizeRange}});
996 :
997 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
998 3 : auto op = BlobProjector(domain, range);
999 :
1000 3 : auto Ax = DataContainer(range);
1001 3 : Ax = 0;
1002 3 : x = 0;
1003 :
1004 3 : WHEN("Backward projecting the operator")
1005 3 : {
1006 3 : Ax[0] = 1;
1007 :
1008 3 : op.applyAdjoint(Ax, x);
1009 :
1010 3 : THEN("The main direction of the ray is equal to the weight")
1011 3 : {
1012 1 : const auto weight = op.weight(0);
1013 1 : CAPTURE(weight);
1014 :
1015 1 : CHECK_EQ(x(2, 0), Approx(weight));
1016 1 : CHECK_EQ(x(2, 1), Approx(weight));
1017 1 : CHECK_EQ(x(2, 2), Approx(weight));
1018 1 : CHECK_EQ(x(2, 3), Approx(weight));
1019 1 : CHECK_EQ(x(2, 4), Approx(weight));
1020 1 : }
1021 :
1022 3 : THEN("The main direction of the ray is equal to the weight")
1023 3 : {
1024 1 : const auto weight = op.weight(1);
1025 1 : CAPTURE(weight);
1026 :
1027 1 : CHECK_EQ(x(3, 0), Approx(weight));
1028 1 : CHECK_EQ(x(3, 1), Approx(weight));
1029 1 : CHECK_EQ(x(3, 2), Approx(weight));
1030 1 : CHECK_EQ(x(3, 3), Approx(weight));
1031 1 : CHECK_EQ(x(3, 4), Approx(weight));
1032 :
1033 1 : CHECK_EQ(x(1, 0), Approx(weight));
1034 1 : CHECK_EQ(x(1, 1), Approx(weight));
1035 1 : CHECK_EQ(x(1, 2), Approx(weight));
1036 1 : CHECK_EQ(x(1, 3), Approx(weight));
1037 1 : CHECK_EQ(x(1, 4), Approx(weight));
1038 1 : }
1039 :
1040 3 : THEN("The main direction of the ray is equal to the weight")
1041 3 : {
1042 1 : const auto weight = op.weight(2);
1043 1 : CAPTURE(weight);
1044 :
1045 1 : CHECK_EQ(x(4, 0), Approx(weight));
1046 1 : CHECK_EQ(x(4, 1), Approx(weight));
1047 1 : CHECK_EQ(x(4, 2), Approx(weight));
1048 1 : CHECK_EQ(x(4, 3), Approx(weight));
1049 1 : CHECK_EQ(x(4, 4), Approx(weight));
1050 :
1051 1 : CHECK_EQ(x(0, 0), Approx(weight));
1052 1 : CHECK_EQ(x(0, 1), Approx(weight));
1053 1 : CHECK_EQ(x(0, 2), Approx(weight));
1054 1 : CHECK_EQ(x(0, 3), Approx(weight));
1055 1 : CHECK_EQ(x(0, 4), Approx(weight));
1056 1 : }
1057 3 : }
1058 3 : }
1059 :
1060 15 : GIVEN("a single detector of size 1, at 180 degree")
1061 15 : {
1062 3 : std::vector<Geometry> geom;
1063 3 : geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{sizeDomain}},
1064 3 : SinogramData2D{Size2D{sizeRange}});
1065 :
1066 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1067 3 : auto op = BlobProjector(domain, range);
1068 :
1069 3 : auto Ax = DataContainer(range);
1070 3 : Ax = 0;
1071 3 : x = 0;
1072 :
1073 3 : WHEN("Backward projecting the operator")
1074 3 : {
1075 3 : Ax[0] = 1;
1076 :
1077 3 : op.applyAdjoint(Ax, x);
1078 :
1079 3 : THEN("The main direction of the ray is equal to the weight")
1080 3 : {
1081 1 : const auto weight = op.weight(0);
1082 1 : CAPTURE(weight);
1083 :
1084 1 : CHECK_EQ(x(2, 0), Approx(weight));
1085 1 : CHECK_EQ(x(2, 1), Approx(weight));
1086 1 : CHECK_EQ(x(2, 2), Approx(weight));
1087 1 : CHECK_EQ(x(2, 3), Approx(weight));
1088 1 : CHECK_EQ(x(2, 4), Approx(weight));
1089 1 : }
1090 :
1091 3 : THEN("The main direction of the ray is equal to the weight")
1092 3 : {
1093 1 : const auto weight = op.weight(1);
1094 1 : CAPTURE(weight);
1095 :
1096 1 : CHECK_EQ(x(3, 0), Approx(weight));
1097 1 : CHECK_EQ(x(3, 1), Approx(weight));
1098 1 : CHECK_EQ(x(3, 2), Approx(weight));
1099 1 : CHECK_EQ(x(3, 3), Approx(weight));
1100 1 : CHECK_EQ(x(3, 4), Approx(weight));
1101 :
1102 1 : CHECK_EQ(x(1, 0), Approx(weight));
1103 1 : CHECK_EQ(x(1, 1), Approx(weight));
1104 1 : CHECK_EQ(x(1, 2), Approx(weight));
1105 1 : CHECK_EQ(x(1, 3), Approx(weight));
1106 1 : CHECK_EQ(x(1, 4), Approx(weight));
1107 1 : }
1108 :
1109 3 : THEN("The main direction of the ray is equal to the weight")
1110 3 : {
1111 1 : const auto weight = op.weight(2);
1112 1 : CAPTURE(weight);
1113 :
1114 1 : CHECK_EQ(x(4, 0), Approx(weight));
1115 1 : CHECK_EQ(x(4, 1), Approx(weight));
1116 1 : CHECK_EQ(x(4, 2), Approx(weight));
1117 1 : CHECK_EQ(x(4, 3), Approx(weight));
1118 1 : CHECK_EQ(x(4, 4), Approx(weight));
1119 :
1120 1 : CHECK_EQ(x(0, 0), Approx(weight));
1121 1 : CHECK_EQ(x(0, 1), Approx(weight));
1122 1 : CHECK_EQ(x(0, 2), Approx(weight));
1123 1 : CHECK_EQ(x(0, 3), Approx(weight));
1124 1 : CHECK_EQ(x(0, 4), Approx(weight));
1125 1 : }
1126 3 : }
1127 3 : }
1128 :
1129 15 : GIVEN("a single detector of size 1, at 90 degree")
1130 15 : {
1131 3 : std::vector<Geometry> geom;
1132 3 : geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{sizeDomain}},
1133 3 : SinogramData2D{Size2D{sizeRange}});
1134 :
1135 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1136 3 : auto op = BlobProjector(domain, range);
1137 :
1138 3 : auto Ax = DataContainer(range);
1139 3 : Ax = 0;
1140 3 : x = 0;
1141 :
1142 3 : WHEN("Backward projecting the operator")
1143 3 : {
1144 3 : Ax[0] = 1;
1145 :
1146 3 : op.applyAdjoint(Ax, x);
1147 :
1148 3 : THEN("The main direction of the ray is equal to the weight")
1149 3 : {
1150 1 : const auto weight = op.weight(0);
1151 1 : CAPTURE(weight);
1152 :
1153 1 : CHECK_EQ(x(0, 2), Approx(weight));
1154 1 : CHECK_EQ(x(1, 2), Approx(weight));
1155 1 : CHECK_EQ(x(2, 2), Approx(weight));
1156 1 : CHECK_EQ(x(3, 2), Approx(weight));
1157 1 : CHECK_EQ(x(4, 2), Approx(weight));
1158 1 : }
1159 :
1160 3 : THEN("The main direction of the ray is equal to the weight")
1161 3 : {
1162 1 : const auto weight = op.weight(1);
1163 1 : CAPTURE(weight);
1164 :
1165 1 : CHECK_EQ(x(0, 3), Approx(weight));
1166 1 : CHECK_EQ(x(1, 3), Approx(weight));
1167 1 : CHECK_EQ(x(2, 3), Approx(weight));
1168 1 : CHECK_EQ(x(3, 3), Approx(weight));
1169 1 : CHECK_EQ(x(4, 3), Approx(weight));
1170 :
1171 1 : CHECK_EQ(x(0, 1), Approx(weight));
1172 1 : CHECK_EQ(x(1, 1), Approx(weight));
1173 1 : CHECK_EQ(x(2, 1), Approx(weight));
1174 1 : CHECK_EQ(x(3, 1), Approx(weight));
1175 1 : CHECK_EQ(x(4, 1), Approx(weight));
1176 1 : }
1177 :
1178 3 : THEN("The main direction of the ray is equal to the weight")
1179 3 : {
1180 1 : const auto weight = op.weight(2);
1181 1 : CAPTURE(weight);
1182 :
1183 1 : CHECK_EQ(x(0, 4), Approx(weight));
1184 1 : CHECK_EQ(x(1, 4), Approx(weight));
1185 1 : CHECK_EQ(x(2, 4), Approx(weight));
1186 1 : CHECK_EQ(x(3, 4), Approx(weight));
1187 1 : CHECK_EQ(x(4, 4), Approx(weight));
1188 :
1189 1 : CHECK_EQ(x(0, 0), Approx(weight));
1190 1 : CHECK_EQ(x(1, 0), Approx(weight));
1191 1 : CHECK_EQ(x(2, 0), Approx(weight));
1192 1 : CHECK_EQ(x(3, 0), Approx(weight));
1193 1 : CHECK_EQ(x(4, 0), Approx(weight));
1194 1 : }
1195 3 : }
1196 3 : }
1197 :
1198 15 : GIVEN("a single detector of size 1, at 270 degree")
1199 15 : {
1200 3 : std::vector<Geometry> geom;
1201 3 : geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{sizeDomain}},
1202 3 : SinogramData2D{Size2D{sizeRange}});
1203 :
1204 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1205 3 : auto op = BlobProjector(domain, range);
1206 :
1207 3 : auto Ax = DataContainer(range);
1208 3 : Ax = 0;
1209 3 : x = 0;
1210 :
1211 3 : WHEN("Backward projecting the operator")
1212 3 : {
1213 3 : Ax[0] = 1;
1214 :
1215 3 : op.applyAdjoint(Ax, x);
1216 :
1217 3 : THEN("The main direction of the ray is equal to the weight")
1218 3 : {
1219 1 : const auto weight = op.weight(0);
1220 1 : CAPTURE(weight);
1221 :
1222 1 : CHECK_EQ(x(0, 2), Approx(weight));
1223 1 : CHECK_EQ(x(1, 2), Approx(weight));
1224 1 : CHECK_EQ(x(2, 2), Approx(weight));
1225 1 : CHECK_EQ(x(3, 2), Approx(weight));
1226 1 : CHECK_EQ(x(4, 2), Approx(weight));
1227 1 : }
1228 :
1229 3 : THEN("The main direction of the ray is equal to the weight")
1230 3 : {
1231 1 : const auto weight = op.weight(1);
1232 1 : CAPTURE(weight);
1233 :
1234 1 : CHECK_EQ(x(0, 3), Approx(weight));
1235 1 : CHECK_EQ(x(1, 3), Approx(weight));
1236 1 : CHECK_EQ(x(2, 3), Approx(weight));
1237 1 : CHECK_EQ(x(3, 3), Approx(weight));
1238 1 : CHECK_EQ(x(4, 3), Approx(weight));
1239 :
1240 1 : CHECK_EQ(x(0, 1), Approx(weight));
1241 1 : CHECK_EQ(x(1, 1), Approx(weight));
1242 1 : CHECK_EQ(x(2, 1), Approx(weight));
1243 1 : CHECK_EQ(x(3, 1), Approx(weight));
1244 1 : CHECK_EQ(x(4, 1), Approx(weight));
1245 1 : }
1246 :
1247 3 : THEN("The main direction of the ray is equal to the weight")
1248 3 : {
1249 1 : const auto weight = op.weight(2);
1250 1 : CAPTURE(weight);
1251 :
1252 1 : CHECK_EQ(x(0, 4), Approx(weight));
1253 1 : CHECK_EQ(x(1, 4), Approx(weight));
1254 1 : CHECK_EQ(x(2, 4), Approx(weight));
1255 1 : CHECK_EQ(x(3, 4), Approx(weight));
1256 1 : CHECK_EQ(x(4, 4), Approx(weight));
1257 :
1258 1 : CHECK_EQ(x(0, 0), Approx(weight));
1259 1 : CHECK_EQ(x(1, 0), Approx(weight));
1260 1 : CHECK_EQ(x(2, 0), Approx(weight));
1261 1 : CHECK_EQ(x(3, 0), Approx(weight));
1262 1 : CHECK_EQ(x(4, 0), Approx(weight));
1263 1 : }
1264 3 : }
1265 3 : }
1266 :
1267 15 : GIVEN("a single detector of size 1, at 45 degree")
1268 15 : {
1269 3 : std::vector<Geometry> geom;
1270 3 : geom.emplace_back(stc, ctr, Degree{45}, VolumeData2D{Size2D{sizeDomain}},
1271 3 : SinogramData2D{Size2D{sizeRange}});
1272 :
1273 3 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1274 3 : auto op = BlobProjector(domain, range);
1275 :
1276 3 : auto Ax = DataContainer(range);
1277 3 : Ax = 0;
1278 3 : x = 0;
1279 :
1280 3 : WHEN("Backward projecting the operator")
1281 3 : {
1282 3 : Ax[0] = 1;
1283 :
1284 3 : op.applyAdjoint(Ax, x);
1285 :
1286 3 : THEN("The main direction of the ray is equal to the weight")
1287 3 : {
1288 1 : const auto weight = op.weight(0);
1289 1 : CAPTURE(weight);
1290 :
1291 1 : CHECK_EQ(x(0, 0), Approx(weight));
1292 1 : CHECK_EQ(x(1, 1), Approx(weight));
1293 1 : CHECK_EQ(x(2, 2), Approx(weight));
1294 1 : CHECK_EQ(x(3, 3), Approx(weight));
1295 1 : CHECK_EQ(x(4, 4), Approx(weight));
1296 1 : }
1297 :
1298 3 : THEN("The main direction of the ray is equal to the weight")
1299 3 : {
1300 1 : const auto distance = 0.707106769084930;
1301 1 : const auto weight = op.weight(distance);
1302 1 : CAPTURE(distance);
1303 1 : CAPTURE(weight);
1304 :
1305 1 : CHECK_EQ(x(0, 1), Approx(weight));
1306 1 : CHECK_EQ(x(1, 2), Approx(weight));
1307 1 : CHECK_EQ(x(2, 3), Approx(weight));
1308 1 : CHECK_EQ(x(3, 4), Approx(weight));
1309 :
1310 1 : CHECK_EQ(x(1, 0), Approx(weight));
1311 1 : CHECK_EQ(x(2, 1), Approx(weight));
1312 1 : CHECK_EQ(x(3, 2), Approx(weight));
1313 1 : CHECK_EQ(x(4, 3), Approx(weight));
1314 1 : }
1315 :
1316 3 : THEN("The main direction of the ray is equal to the weight")
1317 3 : {
1318 1 : const auto distance = 2 * 0.707106769084930;
1319 1 : const auto weight = op.weight(distance);
1320 1 : CAPTURE(weight);
1321 :
1322 1 : CHECK_EQ(x(0, 2), Approx(weight));
1323 1 : CHECK_EQ(x(1, 3), Approx(weight));
1324 1 : CHECK_EQ(x(2, 4), Approx(weight));
1325 :
1326 1 : CHECK_EQ(x(2, 0), Approx(weight));
1327 1 : CHECK_EQ(x(3, 1), Approx(weight));
1328 1 : CHECK_EQ(x(4, 2), Approx(weight));
1329 1 : }
1330 3 : }
1331 3 : }
1332 15 : }
1333 :
1334 : TEST_CASE_TEMPLATE("BlobProjector: Test weights", data_t, float, double)
1335 2 : {
1336 2 : const auto a = 2;
1337 2 : const auto alpha = 10.83;
1338 2 : const auto m = 2;
1339 :
1340 2 : std::array<double, 51> expected{
1341 2 : 1.3671064952680276, 1.3635202864368146, 1.3528128836429958,
1342 2 : 1.3351368521026497, 1.3107428112196733, 1.2799740558068384,
1343 2 : 1.2432592326454344, 1.2011032712842662, 1.1540768137691881,
1344 2 : 1.1028044260488254, 1.0479519029788988, 0.9902129983165086,
1345 2 : 0.930295920364307, 0.868909932853028, 0.80675238945173,
1346 2 : 0.7444965095220905, 0.6827801732549599, 0.6221959772930218,
1347 2 : 0.5632827487344612, 0.5065186675909519, 0.4523160970195944,
1348 2 : 0.40101816870068707, 0.35289711932154216, 0.30815432491147815,
1349 2 : 0.2669219342875695, 0.22926596246745282, 0.1951906707135796,
1350 2 : 0.1646440327638897, 0.13752406736650655, 0.11368580576665363,
1351 2 : 0.09294865929450916, 0.0751039563916854, 0.05992242974801804,
1352 2 : 0.04716145192475515, 0.036571840947646775, 0.027904084748497336,
1353 2 : 0.020913863801848218, 0.015366783581446854, 0.01104226128798172,
1354 2 : 0.007736543464620423, 0.005264861504239613, 0.003462759678911652,
1355 2 : 0.0021866543689359847, 0.0013137030013652602, 0.0007410763873099119,
1356 2 : 0.0003847384204545552, 0.0001778423521946198, 6.885294293780879e-05,
1357 2 : 1.9497773431607567e-05, 2.632583006403572e-06, 0};
1358 :
1359 2 : const IndexVector_t sizeDomain({{5, 5}});
1360 2 : const IndexVector_t sizeRange({{1, 1}});
1361 :
1362 2 : auto domain = VolumeDescriptor(sizeDomain);
1363 2 : auto x = DataContainer(domain);
1364 2 : x = 1;
1365 :
1366 2 : auto stc = SourceToCenterOfRotation{100};
1367 2 : auto ctr = CenterOfRotationToDetector{5};
1368 2 : auto volData = VolumeData2D{Size2D{sizeDomain}};
1369 2 : auto sinoData = SinogramData2D{Size2D{sizeRange}};
1370 :
1371 2 : std::vector<Geometry> geom;
1372 2 : geom.emplace_back(stc, ctr, Degree{0}, std::move(volData), std::move(sinoData));
1373 :
1374 2 : auto range = PlanarDetectorDescriptor(sizeRange, geom);
1375 2 : auto op = BlobProjector(domain, range);
1376 :
1377 102 : for (int i = 0; i < 50; ++i) {
1378 100 : const auto distance = i / 25.;
1379 :
1380 100 : CAPTURE(i);
1381 100 : CAPTURE(distance);
1382 100 : CHECK_EQ(Approx(op.weight(distance)), expected[i]);
1383 100 : }
1384 2 : }
1385 :
1386 : TEST_SUITE_END();
|
