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 : TEST_CASE("LimitedAngleTrajectoryGenerator: Create a mirrored Limited Angle Trajectory")
20 2 : {
21 2 : 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 2 : GIVEN("A 2D descriptor and 256 angles")
29 2 : {
30 2 : index_t numberOfAngles = 256;
31 2 : IndexVector_t volSize(2);
32 2 : volSize << s, s;
33 2 : VolumeDescriptor desc{volSize};
34 :
35 2 : WHEN("We create a half circular limited angle trajectory for this scenario")
36 2 : {
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 1 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
43 1 : numberOfAngles, missingWedgeAngles, desc, halfCircular, diffCenterSource,
44 1 : diffCenterDetector);
45 :
46 : // check that the detector size is correct
47 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
48 1 : expectedDetectorSize);
49 :
50 1 : THEN("Every geomList in our list has the same camera center and the same projection "
51 1 : "matrix")
52 1 : {
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 1 : 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 769 : 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 1 : }
68 :
69 768 : if (!((angle.to_degree() >= missingWedgeAngles.first
70 768 : && angle.to_degree() <= missingWedgeAngles.second)
71 768 : || (angle.to_degree() >= (missingWedgeAngles.first + 180)
72 512 : && angle.to_degree() <= (missingWedgeAngles.second + 180)))) {
73 512 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
74 512 : CenterOfRotationToDetector{centerToDetector},
75 512 : Radian{angle}, VolumeData2D{volSize},
76 512 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
77 512 : sinoDescriptor->getLocationOfOrigin()});
78 :
79 512 : auto geom = sinoDescriptor->getGeometryAt(j++);
80 512 : CHECK(geom);
81 :
82 512 : const auto centerNorm =
83 512 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
84 512 : const auto projMatNorm =
85 512 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
86 512 : 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 512 : }
93 768 : }
94 1 : }
95 1 : }
96 :
97 2 : WHEN("We create a full circular limited angle trajectory for this scenario")
98 2 : {
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 1 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
105 1 : numberOfAngles, missingWedgeAngles, desc, fullyCircular, diffCenterSource,
106 1 : diffCenterDetector);
107 :
108 : // check that the detector size is correct
109 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
110 1 : expectedDetectorSize);
111 :
112 1 : THEN("Every geomList in our list has the same camera center and the same projection "
113 1 : "matrix")
114 1 : {
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 1 : 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 331 : 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 1 : }
130 :
131 330 : if (!((angle.to_degree() >= missingWedgeAngles.first
132 330 : && angle.to_degree() <= missingWedgeAngles.second)
133 330 : || (angle.to_degree() >= (missingWedgeAngles.first + 180)
134 293 : && angle.to_degree() <= (missingWedgeAngles.second + 180)))) {
135 256 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
136 256 : CenterOfRotationToDetector{centerToDetector},
137 256 : Radian{angle}, VolumeData2D{volSize},
138 256 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
139 256 : sinoDescriptor->getLocationOfOrigin()});
140 :
141 256 : auto geom = sinoDescriptor->getGeometryAt(j++);
142 256 : CHECK(geom);
143 :
144 256 : const auto centerNorm =
145 256 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
146 256 : const auto projMatNorm =
147 256 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
148 256 : 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 256 : }
155 330 : }
156 1 : }
157 1 : }
158 2 : }
159 2 : }
160 :
161 : TEST_CASE("LimitedAngleTrajectoryGenerator: Create a non-mirrored Limited Angle Trajectory")
162 2 : {
163 2 : 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 2 : GIVEN("A 2D descriptor and 256 angles")
171 2 : {
172 2 : index_t numberOfAngles = 256;
173 2 : IndexVector_t volSize(2);
174 2 : volSize << s, s;
175 2 : VolumeDescriptor desc{volSize};
176 :
177 2 : WHEN("We create a half circular limited angle trajectory for this scenario")
178 2 : {
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 1 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
185 1 : numberOfAngles, missingWedgeAngles, desc, halfCircular, diffCenterSource,
186 1 : diffCenterDetector, false);
187 :
188 : // check that the detector size is correct
189 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
190 1 : expectedDetectorSize);
191 :
192 1 : THEN("Every geomList in our list has the same camera center and the same projection "
193 1 : "matrix")
194 1 : {
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 1 : 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 265 : 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 1 : }
210 :
211 264 : if (!(angle.to_degree() >= missingWedgeAngles.first
212 264 : && angle.to_degree() <= missingWedgeAngles.second)) {
213 257 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
214 257 : CenterOfRotationToDetector{centerToDetector},
215 257 : Radian{angle}, VolumeData2D{volSize},
216 257 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
217 257 : sinoDescriptor->getLocationOfOrigin()});
218 :
219 257 : auto geom = sinoDescriptor->getGeometryAt(j++);
220 257 : CHECK(geom);
221 :
222 257 : const auto centerNorm =
223 257 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
224 257 : const auto projMatNorm =
225 257 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
226 257 : 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 257 : }
233 264 : }
234 1 : }
235 1 : }
236 :
237 2 : WHEN("We create a full circular limited angle trajectory for this scenario")
238 2 : {
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 1 : auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
245 1 : numberOfAngles, missingWedgeAngles, desc, fullyCircular, diffCenterSource,
246 1 : diffCenterDetector, false);
247 :
248 : // check that the detector size is correct
249 1 : REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
250 1 : expectedDetectorSize);
251 :
252 1 : THEN("Every geomList in our list has the same camera center and the same projection "
253 1 : "matrix")
254 1 : {
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 1 : 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 314 : 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 1 : }
270 :
271 313 : if (!(angle.to_degree() >= missingWedgeAngles.first
272 313 : && angle.to_degree() <= missingWedgeAngles.second)) {
273 256 : Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
274 256 : CenterOfRotationToDetector{centerToDetector},
275 256 : Radian{angle}, VolumeData2D{volSize},
276 256 : SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
277 256 : sinoDescriptor->getLocationOfOrigin()});
278 :
279 256 : auto geom = sinoDescriptor->getGeometryAt(j++);
280 256 : CHECK(geom);
281 :
282 256 : const auto centerNorm =
283 256 : (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
284 256 : const auto projMatNorm =
285 256 : (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
286 256 : 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 256 : }
293 313 : }
294 1 : }
295 1 : }
296 2 : }
297 2 : }
|
