Line data Source code
1 : /**
2 : * @file test_CurvedDetectorDescriptor.cpp
3 : *
4 : * @brief Test for CurvedDetectorDescriptor
5 : *
6 : * @author David Frank - initial code for PlanarDetectorDescriptor
7 : * @author Julia Spindler, Robert Imschweiler - adapt for CurvedDetectorDescriptor
8 : */
9 :
10 : #include "StrongTypes.h"
11 : #include "doctest/doctest.h"
12 :
13 : #include "CurvedDetectorDescriptor.h"
14 : #include "VolumeDescriptor.h"
15 : #include <iostream>
16 :
17 : using namespace elsa;
18 : using namespace elsa::geometry;
19 : using namespace doctest;
20 :
21 : TEST_SUITE_BEGIN("core");
22 :
23 : const geometry::Radian angleFakePlanar{static_cast<real_t>(1e-10)};
24 :
25 : TEST_CASE("CurvedDetectorDescriptor: Testing 2D CurvedDetectorDescriptor as fake "
26 : "PlanarDetectorDescriptor")
27 7 : {
28 7 : GIVEN("Given a 5x5 Volume and a single 5 wide detector pose")
29 7 : {
30 7 : IndexVector_t volSize(2);
31 7 : volSize << 5, 5;
32 7 : VolumeDescriptor ddVol(volSize);
33 :
34 7 : IndexVector_t sinoSize(2);
35 7 : sinoSize << 5, 1;
36 :
37 7 : real_t s2c = 10;
38 7 : real_t c2d = 4;
39 :
40 7 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Radian{0},
41 7 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
42 :
43 7 : CurvedDetectorDescriptor desc(sinoSize, {g}, angleFakePlanar, s2c + c2d);
44 :
45 7 : WHEN("Retreiving the single geometry pose")
46 7 : {
47 1 : auto geom = desc.getGeometryAt(0);
48 1 : auto geomList = desc.getGeometry();
49 :
50 1 : CHECK_EQ(desc.getNumberOfGeometryPoses(), 1);
51 1 : CHECK_EQ(geomList.size(), 1);
52 :
53 1 : THEN("Geometry is equal")
54 1 : {
55 1 : CHECK_EQ((geom), g);
56 1 : CHECK_EQ(geomList[0], g);
57 1 : }
58 1 : }
59 :
60 7 : WHEN("Generating rays for detecor pixels 0, 2 and 4")
61 7 : {
62 3 : for (real_t detPixel : std::initializer_list<real_t>{0, 2, 4}) {
63 3 : RealVector_t pixel(1);
64 3 : pixel << detPixel + 0.5f;
65 :
66 : // Check that ray for IndexVector_t is equal to previous one
67 3 : auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
68 :
69 : // Create variables, which make typing quicker
70 3 : auto ro = ray.origin();
71 3 : auto rd = ray.direction();
72 :
73 : // Check that ray origin is camera center
74 3 : auto c = g.getCameraCenter();
75 3 : CHECK_EQ((ro - c).sum(), Approx(0));
76 :
77 : // compute intersection manually
78 3 : real_t t = Approx(rd[0]) == 0 ? (s2c + c2d) : ((pixel[0] - ro[0]) / rd[0]);
79 :
80 3 : auto detCoordY = ro[1] + t * rd[1];
81 :
82 3 : CHECK_EQ(detCoordY, Approx(ddVol.getLocationOfOrigin()[1] + c2d));
83 3 : }
84 1 : }
85 :
86 7 : GIVEN("Given a 5x5 Volume and a multiple 5 wide detector pose")
87 7 : {
88 5 : IndexVector_t volSize(2);
89 5 : volSize << 5, 5;
90 5 : VolumeDescriptor ddVol(volSize);
91 :
92 5 : IndexVector_t sinoSize(2);
93 5 : sinoSize << 5, 4;
94 :
95 5 : Geometry g1(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{0},
96 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
97 5 : Geometry g2(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{90},
98 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
99 5 : Geometry g3(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{180},
100 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
101 5 : Geometry g4(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{270},
102 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
103 :
104 5 : CurvedDetectorDescriptor desc(sinoSize, {g1, g2, g3, g4}, angleFakePlanar, s2c + c2d);
105 :
106 5 : WHEN("Retreiving geometry poses")
107 5 : {
108 4 : auto geom = desc.getGeometryAt(0);
109 4 : auto geomList = desc.getGeometry();
110 :
111 4 : CHECK_EQ(desc.getNumberOfGeometryPoses(), 4);
112 4 : CHECK_EQ(geomList.size(), 4);
113 :
114 4 : THEN("Geometry is equal")
115 4 : {
116 1 : CHECK_EQ((desc.getGeometryAt(0)), g1);
117 1 : CHECK_EQ(geomList[0], g1);
118 1 : }
119 4 : THEN("Geometry is equal")
120 4 : {
121 1 : CHECK_EQ((desc.getGeometryAt(1)), g2);
122 1 : CHECK_EQ(geomList[1], g2);
123 1 : }
124 4 : THEN("Geometry is equal")
125 4 : {
126 1 : CHECK_EQ((desc.getGeometryAt(2)), g3);
127 1 : CHECK_EQ(geomList[2], g3);
128 1 : }
129 4 : THEN("Geometry is equal")
130 4 : {
131 1 : CHECK_EQ((desc.getGeometryAt(3)), g4);
132 1 : CHECK_EQ(geomList[3], g4);
133 1 : }
134 4 : }
135 :
136 5 : WHEN("Check for multiple poses, that all the overloads compute the same rays")
137 5 : {
138 4 : for (index_t pose : {0, 1, 2, 3}) {
139 :
140 12 : for (index_t detPixel : {0, 2, 4}) {
141 12 : IndexVector_t pixel(2);
142 12 : pixel << detPixel, pose;
143 :
144 12 : RealVector_t pixelReal(1);
145 12 : pixelReal << static_cast<real_t>(detPixel) + 0.5f;
146 :
147 12 : auto ray1 =
148 12 : desc.computeRayFromDetectorCoord(desc.getIndexFromCoordinate(pixel));
149 12 : auto ray2 = desc.computeRayFromDetectorCoord(pixel);
150 12 : auto ray3 = desc.computeRayFromDetectorCoord(pixelReal, pose);
151 :
152 12 : auto ro1 = ray1.origin();
153 12 : auto rd1 = ray1.direction();
154 :
155 12 : auto ro2 = ray2.origin();
156 12 : auto rd2 = ray2.direction();
157 :
158 12 : auto ro3 = ray3.origin();
159 12 : auto rd3 = ray3.direction();
160 :
161 12 : CHECK_EQ(ro1, ro2);
162 12 : CHECK_EQ(ro1, ro3);
163 :
164 12 : CHECK_EQ(rd1, rd2);
165 12 : CHECK_EQ(rd1, rd3);
166 :
167 : // Shouldn't be necessary, but whatever
168 12 : CHECK_EQ(ro2, ro3);
169 12 : CHECK_EQ(rd2, rd3);
170 12 : }
171 4 : }
172 1 : }
173 5 : }
174 7 : }
175 7 : }
176 :
177 : TEST_CASE("CurvedDetectorDescriptor: Testing 3D CurvedDetectorDescriptor as fake "
178 : "PlanarDetectorDescriptor")
179 2 : {
180 2 : GIVEN("Given a 5x5x5 Volume and a single 5x5 wide detector pose")
181 2 : {
182 2 : IndexVector_t volSize(3);
183 2 : volSize << 5, 5, 5;
184 2 : VolumeDescriptor ddVol(volSize);
185 :
186 2 : IndexVector_t sinoSize(3);
187 2 : sinoSize << 5, 5, 1;
188 :
189 2 : real_t s2c = 10;
190 2 : real_t c2d = 4;
191 :
192 2 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
193 2 : VolumeData3D{Size3D{volSize}}, SinogramData3D{Size3D{sinoSize}},
194 2 : RotationAngles3D{Gamma{0}});
195 :
196 2 : CurvedDetectorDescriptor desc(sinoSize, {g}, angleFakePlanar, s2c + c2d);
197 :
198 2 : WHEN("Retreiving the single geometry pose")
199 2 : {
200 1 : auto geom = desc.getGeometryAt(0);
201 1 : auto geomList = desc.getGeometry();
202 :
203 1 : CHECK_EQ(desc.getNumberOfGeometryPoses(), 1);
204 1 : CHECK_EQ(geomList.size(), 1);
205 :
206 1 : THEN("Geometry is equal")
207 1 : {
208 1 : CHECK_EQ((geom), g);
209 1 : CHECK_EQ(geomList[0], g);
210 1 : }
211 1 : }
212 :
213 2 : WHEN("Generating rays for detecor pixels 0, 2 and 4 for each dim")
214 2 : {
215 3 : for (index_t detPixel1 : {0, 2, 4}) {
216 9 : for (index_t detPixel2 : {0, 2, 4}) {
217 9 : RealVector_t pixel(2);
218 9 : pixel << static_cast<real_t>(detPixel1) + 0.5f,
219 9 : static_cast<real_t>(detPixel2) + 0.5f;
220 :
221 : // Check that ray for IndexVector_t is equal to previous one
222 9 : auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
223 :
224 : // Create variables, which make typing quicker
225 9 : auto ro = ray.origin();
226 9 : auto rd = ray.direction();
227 :
228 : // Check that ray origin is camera center
229 9 : auto c = g.getCameraCenter();
230 9 : CHECK_EQ((ro - c).sum(), Approx(0));
231 :
232 9 : auto o = ddVol.getLocationOfOrigin();
233 9 : RealVector_t detCoordWorld(3);
234 9 : detCoordWorld << pixel[0] - o[0], pixel[1] - o[1], c2d;
235 9 : RealVector_t rotD = g.getRotationMatrix().transpose() * detCoordWorld + o;
236 :
237 9 : real_t factor = 0;
238 9 : if (std::abs(rd[0]) > 0)
239 6 : factor = (rotD[0] - ro[0]) / rd[0];
240 3 : else if (std::abs(rd[1]) > 0)
241 2 : factor = (rotD[1] - ro[1]) / rd[1];
242 1 : else if (std::abs(rd[2]) > 0)
243 1 : factor = (rotD[2] - ro[2] / rd[2]);
244 :
245 9 : CHECK_EQ((ro[0] + factor * rd[0]), Approx(rotD[0]));
246 9 : CHECK_EQ((ro[1] + factor * rd[1]), Approx(rotD[1]));
247 9 : CHECK_EQ((ro[2] + factor * rd[2]), Approx(rotD[2]));
248 9 : }
249 3 : }
250 1 : }
251 2 : }
252 2 : }
253 :
254 : TEST_CASE("CurvedDetectorDescriptor: Testing 2D CurvedDetectorDescriptor translation to planar "
255 : "coordinates")
256 2 : {
257 2 : real_t angle{2.0};
258 2 : real_t s2d{14.0};
259 :
260 2 : GIVEN("Given a single 5 wide detector pose")
261 2 : {
262 1 : IndexVector_t sinoSize(2);
263 1 : sinoSize << 5, 1;
264 :
265 1 : CurvedDetectorDescriptor desc(sinoSize, {}, geometry::Radian(angle), s2d);
266 :
267 1 : const std::vector<RealVector_t>& planarCoords{desc.getPlanarCoords()};
268 :
269 1 : real_t radius{static_cast<real_t>(sinoSize[0]) / angle};
270 1 : CHECK_EQ(radius, desc.getRadius());
271 :
272 : /*
273 : * The detector has 5 pixels
274 : * -> the pixel in the center has index 2
275 : * -> its center has the x-coordinate 2.5
276 : */
277 1 : CHECK_EQ(planarCoords[2][0], 2.5);
278 :
279 : // Towards the edges the pixel size gets smaller
280 1 : CHECK_EQ(planarCoords[0][0], Approx(0.60392));
281 1 : CHECK_EQ(planarCoords[1][0], Approx(1.51253));
282 1 : CHECK_EQ(planarCoords[3][0], Approx(3.48747));
283 1 : CHECK_EQ(planarCoords[4][0], Approx(4.39608));
284 :
285 : // Check that the spacing is symmetrical
286 1 : CHECK_EQ(2.5 - planarCoords[0][0], Approx(abs(2.5 - planarCoords[4][0])));
287 1 : CHECK_EQ(2.5 - planarCoords[1][0], Approx(abs(2.5 - planarCoords[3][0])));
288 1 : }
289 :
290 2 : GIVEN("Given a single 6 wide detector pose")
291 2 : {
292 1 : IndexVector_t sinoSize(2);
293 1 : sinoSize << 6, 1;
294 :
295 1 : CurvedDetectorDescriptor desc(sinoSize, {}, geometry::Radian(angle), s2d);
296 :
297 1 : const std::vector<RealVector_t>& planarCoords{desc.getPlanarCoords()};
298 :
299 1 : real_t radius{static_cast<real_t>(sinoSize[0]) / angle};
300 1 : CHECK_EQ(radius, desc.getRadius());
301 :
302 : /*
303 : * The detector has 6 pixels
304 : * -> there is no pixel in the center
305 : * -> the center is at 3
306 : */
307 :
308 : // Towards the outer pixels, the pixel size gets smaller.
309 1 : CHECK_EQ(planarCoords[0][0], Approx(0.61183));
310 1 : CHECK_EQ(planarCoords[1][0], Approx(1.52298));
311 1 : CHECK_EQ(planarCoords[2][0], Approx(2.50083));
312 1 : CHECK_EQ(planarCoords[3][0], Approx(3.49917));
313 1 : CHECK_EQ(planarCoords[4][0], Approx(4.47702));
314 1 : CHECK_EQ(planarCoords[5][0], Approx(5.38817));
315 :
316 : // check that the spacing is symmetrical
317 1 : CHECK_EQ(3.0 - planarCoords[0][0], Approx(abs(3.0 - planarCoords[5][0])));
318 1 : CHECK_EQ(3.0 - planarCoords[1][0], Approx(abs(3.0 - planarCoords[4][0])));
319 1 : CHECK_EQ(3.0 - planarCoords[2][0], Approx(abs(3.0 - planarCoords[3][0])));
320 1 : }
321 2 : }
322 :
323 : TEST_CASE("CurvedDetectorDescriptor: Testing 3D CurvedDetectorDescriptor translation to planar "
324 : "coordinates")
325 2 : {
326 2 : real_t angle{2.0};
327 2 : real_t s2d{14.0};
328 :
329 2 : GIVEN("Given a single 5x5 wide detector pose")
330 2 : {
331 1 : IndexVector_t sinoSize(3);
332 1 : sinoSize << 5, 5, 1;
333 :
334 1 : CurvedDetectorDescriptor desc(sinoSize, {}, geometry::Radian(angle), s2d);
335 :
336 1 : const std::vector<RealVector_t>& planarCoords{desc.getPlanarCoords()};
337 :
338 1 : real_t radius{static_cast<real_t>(sinoSize[0]) / angle};
339 1 : CHECK_EQ(radius, desc.getRadius());
340 :
341 : /*
342 : * The detector has 5x5 pixels
343 : * -> the column in the center starts at index 2
344 : * -> x_value stays the same, y_value counts upwards
345 : * (y_value denotes the width of the detector)
346 : */
347 1 : float x_value = 2.5;
348 1 : float y_value_start = 0.5f;
349 1 : CHECK_EQ(planarCoords[2][0], Approx(x_value));
350 1 : CHECK_EQ(planarCoords[2][1], Approx(y_value_start));
351 1 : CHECK_EQ(planarCoords[7][0], Approx(x_value));
352 1 : CHECK_EQ(planarCoords[7][1], Approx(y_value_start + 1));
353 1 : CHECK_EQ(planarCoords[12][0], Approx(x_value));
354 1 : CHECK_EQ(planarCoords[12][1], Approx(y_value_start + 2));
355 1 : CHECK_EQ(planarCoords[17][0], Approx(x_value));
356 1 : CHECK_EQ(planarCoords[17][1], Approx(y_value_start + 3));
357 1 : CHECK_EQ(planarCoords[22][0], Approx(x_value));
358 1 : CHECK_EQ(planarCoords[22][1], Approx(y_value_start + 4));
359 :
360 : // check for symmetry in the last row and last column
361 1 : CHECK_EQ(2.5 - planarCoords[20][0], Approx(abs(2.5 - planarCoords[24][0])));
362 1 : CHECK_EQ(2.5 - planarCoords[21][0], Approx(abs(2.5 - planarCoords[23][0])));
363 1 : CHECK_EQ(2.5 - planarCoords[9][1], Approx(abs(2.5 - planarCoords[19][1])));
364 1 : CHECK_EQ(2.5 - planarCoords[4][1], Approx(abs(2.5 - planarCoords[24][1])));
365 :
366 : // check that the x-value deltas get smaller towards the edges as the source is far away
367 1 : CHECK_LT(planarCoords[1][0] - planarCoords[0][0], planarCoords[2][0] - planarCoords[1][0]);
368 1 : CHECK_LT(planarCoords[4][0] - planarCoords[3][0], planarCoords[2][0] - planarCoords[1][0]);
369 :
370 : // check that the y-value deltas get bigger towards the edges
371 1 : CHECK_LT(planarCoords[0][1], planarCoords[1][1]);
372 1 : CHECK_GT(planarCoords[20][1], planarCoords[21][1]);
373 1 : }
374 :
375 2 : GIVEN("Given a single 6x8 wide detector pose")
376 2 : {
377 1 : IndexVector_t sinoSize(3);
378 1 : sinoSize << 6, 8, 1;
379 :
380 1 : CurvedDetectorDescriptor desc(sinoSize, {}, geometry::Radian(angle), s2d);
381 :
382 1 : const std::vector<RealVector_t>& planarCoords{desc.getPlanarCoords()};
383 :
384 1 : real_t radius{static_cast<real_t>(sinoSize[0]) / angle};
385 1 : CHECK_EQ(radius, desc.getRadius());
386 :
387 : /*
388 : * The detector has 6 pixels
389 : * -> there is no pixel in the center
390 : */
391 :
392 : // check for symmetry in the last row and last column
393 1 : CHECK_EQ(3 - planarCoords[42][0], Approx(abs(3 - planarCoords[47][0])));
394 1 : CHECK_EQ(3 - planarCoords[43][0], Approx(abs(3 - planarCoords[46][0])));
395 1 : CHECK_EQ(3 - planarCoords[44][0], Approx(abs(3 - planarCoords[45][0])));
396 1 : CHECK_EQ(4 - planarCoords[5][1], Approx(abs(4 - planarCoords[47][1])));
397 1 : CHECK_EQ(4 - planarCoords[11][1], Approx(abs(4 - planarCoords[41][1])));
398 1 : CHECK_EQ(4 - planarCoords[17][1], Approx(abs(4 - planarCoords[35][1])));
399 1 : CHECK_EQ(4 - planarCoords[23][1], Approx(abs(4 - planarCoords[29][1])));
400 :
401 : // check that the x-value deltas get smaller towards the edges
402 1 : CHECK_LT(planarCoords[1][0] - planarCoords[0][0], planarCoords[3][0] - planarCoords[2][0]);
403 1 : CHECK_LT(planarCoords[5][0] - planarCoords[4][0], planarCoords[3][0] - planarCoords[2][0]);
404 :
405 : // check that the y-value deltas get bigger towards the edges
406 1 : CHECK_LT(planarCoords[0][1], planarCoords[1][1]);
407 1 : CHECK_GT(planarCoords[42][1], planarCoords[43][1]);
408 1 : }
409 2 : }
410 :
411 : TEST_SUITE_END();
|
