Line data Source code
1 : /**
2 : * @file test_LimitedAngleTrajectoryGenerator.cpp
3 : *
4 : * @brief Tests for the LimitedAngleTrajectoryGenerator class
5 : *
6 : * @author Andi Braimllari
7 : */
8 :
9 : #include "LimitedAngleTrajectoryGenerator.h"
10 : #include "Logger.h"
11 : #include "VolumeDescriptor.h"
12 :
13 : #include "testHelpers.h"
14 : #include "doctest/doctest.h"
15 :
16 : using namespace elsa;
17 : using namespace doctest;
18 :
19 2 : TEST_CASE("LimitedAngleTrajectoryGenerator: Create a mirrored Limited Angle Trajectory")
20 : {
21 : using namespace geometry;
22 :
23 2 : const index_t s = 64;
24 :
25 : // Detector size is the volume size scaled by the square root of 2
26 2 : const auto expectedDetectorSize = static_cast<index_t>(s * std::sqrt(2));
27 :
28 4 : GIVEN("A 2D descriptor and 256 angles")
29 : {
30 2 : index_t numberOfAngles = 256;
31 4 : IndexVector_t volSize(2);
32 2 : volSize << s, s;
33 4 : VolumeDescriptor desc{volSize};
34 :
35 3 : WHEN("We create a half circular limited angle trajectory for this scenario")
36 : {
37 1 : index_t halfCircular = 180;
38 1 : real_t diffCenterSource{s * 100};
39 1 : real_t diffCenterDetector{s};
40 1 : std::pair missingWedgeAngles(geometry::Degree(110), geometry::Degree(170));
41 :
42 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
43 : numberOfAngles, missingWedgeAngles, desc, halfCircular, diffCenterSource,
44 2 : diffCenterDetector);
45 :
46 : // check that the detector size is correct
47 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
48 : expectedDetectorSize);
49 :
50 2 : THEN("Every geomList in our list has the same camera center and the same projection "
51 : "matrix")
52 : {
53 1 : const real_t sourceToCenter = diffCenterSource;
54 1 : const real_t centerToDetector = diffCenterDetector;
55 :
56 1 : real_t wedgeArc = 2 * (missingWedgeAngles.second - missingWedgeAngles.first);
57 :
58 2 : const real_t angleIncrement = (static_cast<real_t>(halfCircular) - wedgeArc)
59 1 : / (static_cast<real_t>(numberOfAngles));
60 :
61 1 : index_t j = 0;
62 1 : for (index_t i = 0;; ++i) {
63 769 : Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
64 :
65 769 : if (angle.to_degree() >= static_cast<real_t>(halfCircular)) {
66 1 : break;
67 : }
68 :
69 1536 : if (!((angle.to_degree() >= missingWedgeAngles.first
70 298 : && angle.to_degree() <= missingWedgeAngles.second)
71 512 : || (angle.to_degree() >= (missingWedgeAngles.first + 180)
72 0 : && angle.to_degree() <= (missingWedgeAngles.second + 180)))) {
73 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
74 : CenterOfRotationToDetector{centerToDetector},
75 1024 : Radian{angle}, VolumeData2D{volSize},
76 1536 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
77 2560 : sinoDescriptor->getLocationOfOrigin()});
78 :
79 512 : auto geom = sinoDescriptor->getGeometryAt(j++);
80 512 : CHECK(geom);
81 :
82 : const auto centerNorm =
83 512 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
84 : const auto projMatNorm =
85 512 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
86 1024 : const auto invProjMatNorm = (tmpGeom.getInverseProjectionMatrix()
87 512 : - geom->getInverseProjectionMatrix())
88 512 : .norm();
89 512 : REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
90 512 : REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
91 512 : REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
92 : }
93 768 : }
94 : }
95 : }
96 :
97 3 : WHEN("We create a full circular limited angle trajectory for this scenario")
98 : {
99 1 : index_t fullyCircular = 359;
100 1 : real_t diffCenterSource{s * 100};
101 1 : real_t diffCenterDetector{s};
102 1 : std::pair missingWedgeAngles(geometry::Degree(30), geometry::Degree(70));
103 :
104 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
105 : numberOfAngles, missingWedgeAngles, desc, fullyCircular, diffCenterSource,
106 2 : diffCenterDetector);
107 :
108 : // check that the detector size is correct
109 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
110 : expectedDetectorSize);
111 :
112 2 : THEN("Every geomList in our list has the same camera center and the same projection "
113 : "matrix")
114 : {
115 1 : const real_t sourceToCenter = diffCenterSource;
116 1 : const real_t centerToDetector = diffCenterDetector;
117 :
118 1 : real_t wedgeArc = 2 * (missingWedgeAngles.second - missingWedgeAngles.first);
119 :
120 2 : const real_t angleIncrement = (static_cast<real_t>(fullyCircular) - wedgeArc)
121 1 : / (static_cast<real_t>(numberOfAngles));
122 :
123 1 : index_t j = 0;
124 1 : for (index_t i = 0;; ++i) {
125 331 : Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
126 :
127 331 : if (angle.to_degree() > static_cast<real_t>(fullyCircular)) {
128 1 : break;
129 : }
130 :
131 660 : if (!((angle.to_degree() >= missingWedgeAngles.first
132 302 : && angle.to_degree() <= missingWedgeAngles.second)
133 293 : || (angle.to_degree() >= (missingWedgeAngles.first + 180)
134 137 : && angle.to_degree() <= (missingWedgeAngles.second + 180)))) {
135 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
136 : CenterOfRotationToDetector{centerToDetector},
137 512 : Radian{angle}, VolumeData2D{volSize},
138 768 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
139 1280 : sinoDescriptor->getLocationOfOrigin()});
140 :
141 256 : auto geom = sinoDescriptor->getGeometryAt(j++);
142 256 : CHECK(geom);
143 :
144 : const auto centerNorm =
145 256 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
146 : const auto projMatNorm =
147 256 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
148 512 : const auto invProjMatNorm = (tmpGeom.getInverseProjectionMatrix()
149 256 : - geom->getInverseProjectionMatrix())
150 256 : .norm();
151 256 : REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
152 256 : REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
153 256 : REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
154 : }
155 330 : }
156 : }
157 : }
158 : }
159 2 : }
160 :
161 2 : TEST_CASE("LimitedAngleTrajectoryGenerator: Create a non-mirrored Limited Angle Trajectory")
162 : {
163 : using namespace geometry;
164 :
165 2 : const index_t s = 64;
166 :
167 : // Detector size is the volume size scaled by the square root of 2
168 2 : const auto expectedDetectorSize = static_cast<index_t>(s * std::sqrt(2));
169 :
170 4 : GIVEN("A 2D descriptor and 256 angles")
171 : {
172 2 : index_t numberOfAngles = 256;
173 4 : IndexVector_t volSize(2);
174 2 : volSize << s, s;
175 4 : VolumeDescriptor desc{volSize};
176 :
177 3 : WHEN("We create a half circular limited angle trajectory for this scenario")
178 : {
179 1 : index_t halfCircular = 180;
180 1 : real_t diffCenterSource{s * 100};
181 1 : real_t diffCenterDetector{s};
182 1 : std::pair missingWedgeAngles(geometry::Degree(40), geometry::Degree(45));
183 :
184 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
185 : numberOfAngles, missingWedgeAngles, desc, halfCircular, diffCenterSource,
186 2 : diffCenterDetector, false);
187 :
188 : // check that the detector size is correct
189 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
190 : expectedDetectorSize);
191 :
192 2 : THEN("Every geomList in our list has the same camera center and the same projection "
193 : "matrix")
194 : {
195 1 : const real_t sourceToCenter = diffCenterSource;
196 1 : const real_t centerToDetector = diffCenterDetector;
197 :
198 1 : real_t wedgeArc = missingWedgeAngles.second - missingWedgeAngles.first;
199 :
200 2 : const real_t angleIncrement = (static_cast<real_t>(halfCircular) - wedgeArc)
201 1 : / (static_cast<real_t>(numberOfAngles));
202 :
203 1 : index_t j = 0;
204 1 : for (index_t i = 0;; ++i) {
205 265 : Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
206 :
207 265 : if (angle.to_degree() > static_cast<real_t>(halfCircular)) {
208 1 : break;
209 : }
210 :
211 469 : if (!(angle.to_degree() >= missingWedgeAngles.first
212 205 : && angle.to_degree() <= missingWedgeAngles.second)) {
213 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
214 : CenterOfRotationToDetector{centerToDetector},
215 514 : Radian{angle}, VolumeData2D{volSize},
216 771 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
217 1285 : sinoDescriptor->getLocationOfOrigin()});
218 :
219 257 : auto geom = sinoDescriptor->getGeometryAt(j++);
220 257 : CHECK(geom);
221 :
222 : const auto centerNorm =
223 257 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
224 : const auto projMatNorm =
225 257 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
226 514 : const auto invProjMatNorm = (tmpGeom.getInverseProjectionMatrix()
227 257 : - geom->getInverseProjectionMatrix())
228 257 : .norm();
229 257 : REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
230 257 : REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
231 257 : REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
232 : }
233 264 : }
234 : }
235 : }
236 :
237 3 : WHEN("We create a full circular limited angle trajectory for this scenario")
238 : {
239 1 : index_t fullyCircular = 359;
240 1 : real_t diffCenterSource{s * 100};
241 1 : real_t diffCenterDetector{s};
242 1 : std::pair missingWedgeAngles(geometry::Degree(20), geometry::Degree(85));
243 :
244 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
245 : numberOfAngles, missingWedgeAngles, desc, fullyCircular, diffCenterSource,
246 2 : diffCenterDetector, false);
247 :
248 : // check that the detector size is correct
249 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
250 : expectedDetectorSize);
251 :
252 2 : THEN("Every geomList in our list has the same camera center and the same projection "
253 : "matrix")
254 : {
255 1 : const real_t sourceToCenter = diffCenterSource;
256 1 : const real_t centerToDetector = diffCenterDetector;
257 :
258 1 : real_t wedgeArc = missingWedgeAngles.second - missingWedgeAngles.first;
259 :
260 2 : const real_t angleIncrement = (static_cast<real_t>(fullyCircular) - wedgeArc)
261 1 : / (static_cast<real_t>(numberOfAngles));
262 :
263 1 : index_t j = 0;
264 1 : for (index_t i = 0;; ++i) {
265 314 : Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
266 :
267 314 : if (angle.to_degree() > static_cast<real_t>(fullyCircular)) {
268 1 : break;
269 : }
270 :
271 608 : if (!(angle.to_degree() >= missingWedgeAngles.first
272 295 : && angle.to_degree() <= missingWedgeAngles.second)) {
273 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
274 : CenterOfRotationToDetector{centerToDetector},
275 512 : Radian{angle}, VolumeData2D{volSize},
276 768 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
277 1280 : sinoDescriptor->getLocationOfOrigin()});
278 :
279 256 : auto geom = sinoDescriptor->getGeometryAt(j++);
280 256 : CHECK(geom);
281 :
282 : const auto centerNorm =
283 256 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
284 : const auto projMatNorm =
285 256 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
286 512 : const auto invProjMatNorm = (tmpGeom.getInverseProjectionMatrix()
287 256 : - geom->getInverseProjectionMatrix())
288 256 : .norm();
289 256 : REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
290 256 : REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
291 256 : REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
292 : }
293 313 : }
294 : }
295 : }
296 : }
297 2 : }
|
