Line data Source code
1 : /**
2 : * @file test_PlanarDetectorDescriptor.cpp
3 : *
4 : * @brief Test for PlanarDetectorDescriptor
5 : *
6 : * @author David Frank - initial code
7 : */
8 :
9 : #include "doctest/doctest.h"
10 :
11 : #include "PlanarDetectorDescriptor.h"
12 : #include "VolumeDescriptor.h"
13 :
14 : using namespace elsa;
15 : using namespace elsa::geometry;
16 : using namespace doctest;
17 :
18 : TEST_SUITE_BEGIN("core");
19 :
20 : TEST_CASE("PlanarDetectorDescriptor: Testing 2D PlanarDetectorDescriptor")
21 10 : {
22 10 : GIVEN("Given a 5x5 Volume and a single 5 wide detector pose")
23 10 : {
24 10 : IndexVector_t volSize(2);
25 10 : volSize << 5, 5;
26 10 : VolumeDescriptor ddVol(volSize);
27 :
28 10 : IndexVector_t sinoSize(2);
29 10 : sinoSize << 5, 1;
30 :
31 10 : real_t s2c = 10;
32 10 : real_t c2d = 4;
33 :
34 10 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Radian{0},
35 10 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
36 :
37 10 : PlanarDetectorDescriptor desc(sinoSize, {g});
38 :
39 10 : WHEN("Retreiving the single geometry pose")
40 10 : {
41 1 : auto geom = desc.getGeometryAt(0);
42 1 : auto geomList = desc.getGeometry();
43 :
44 1 : CHECK_EQ(desc.getNumberOfGeometryPoses(), 1);
45 1 : CHECK_EQ(geomList.size(), 1);
46 :
47 1 : THEN("Geometry is equal")
48 1 : {
49 1 : CHECK_EQ((geom), g);
50 1 : CHECK_EQ(geomList[0], g);
51 1 : }
52 1 : }
53 :
54 10 : WHEN("Generating rays for detecor pixels 0, 2 and 4")
55 10 : {
56 3 : for (real_t detPixel : std::initializer_list<real_t>{0, 2, 4}) {
57 3 : RealVector_t pixel(1);
58 3 : pixel << detPixel + 0.5f;
59 :
60 : // Check that ray for IndexVector_t is equal to previous one
61 3 : auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
62 :
63 : // Create variables, which make typing quicker
64 3 : auto ro = ray.origin();
65 3 : auto rd = ray.direction();
66 :
67 : // Check that ray origin is camera center
68 3 : auto c = g.getCameraCenter();
69 3 : CHECK_EQ((ro - c).sum(), Approx(0));
70 :
71 : // compute intersection manually
72 3 : real_t t = Approx(rd[0]) == 0 ? (s2c + c2d) : ((pixel[0] - ro[0]) / rd[0]);
73 :
74 3 : auto detCoordY = ro[1] + t * rd[1];
75 :
76 3 : CHECK_EQ(detCoordY, Approx(ddVol.getLocationOfOrigin()[1] + c2d));
77 3 : }
78 1 : }
79 10 : WHEN("Computing detector coord from ray")
80 10 : {
81 3 : auto ro = g.getCameraCenter();
82 3 : auto rd = RealVector_t(2);
83 :
84 3 : THEN("The detector coord is correct")
85 3 : {
86 1 : rd << -0.141421f, 0.989949f;
87 1 : rd.normalize();
88 :
89 1 : auto detCoord = desc.computeDetectorCoordFromRay(RealRay_t(ro, rd), 0);
90 1 : CHECK_EQ(detCoord[0], Approx(0.5).epsilon(0.05));
91 1 : }
92 3 : THEN("The detector coord is correct asd")
93 3 : {
94 1 : rd << 0.f, 1.f;
95 1 : rd.normalize();
96 :
97 1 : auto detCoord = desc.computeDetectorCoordFromRay(RealRay_t(ro, rd), 0);
98 1 : CHECK_EQ(detCoord[0], Approx(2.5));
99 1 : }
100 :
101 3 : THEN("The detector coord is correct asd")
102 3 : {
103 1 : rd << 0.141421f, 0.989949f;
104 1 : rd.normalize();
105 :
106 1 : auto detCoord = desc.computeDetectorCoordFromRay(RealRay_t(ro, rd), 0);
107 1 : CHECK_EQ(detCoord[0], Approx(4.5).epsilon(0.05));
108 1 : }
109 3 : }
110 :
111 10 : GIVEN("Given a 5x5 Volume and a multiple 5 wide detector pose")
112 10 : {
113 5 : IndexVector_t volSize(2);
114 5 : volSize << 5, 5;
115 5 : VolumeDescriptor ddVol(volSize);
116 :
117 5 : IndexVector_t sinoSize(2);
118 5 : sinoSize << 5, 4;
119 :
120 5 : real_t s2c = 10;
121 5 : real_t c2d = 4;
122 :
123 5 : Geometry g1(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{0},
124 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
125 5 : Geometry g2(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{90},
126 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
127 5 : Geometry g3(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{180},
128 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
129 5 : Geometry g4(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{270},
130 5 : VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
131 :
132 5 : PlanarDetectorDescriptor desc(sinoSize, {g1, g2, g3, g4});
133 :
134 5 : WHEN("Retreiving geometry poses")
135 5 : {
136 4 : auto geom = desc.getGeometryAt(0);
137 4 : auto geomList = desc.getGeometry();
138 :
139 4 : CHECK_EQ(desc.getNumberOfGeometryPoses(), 4);
140 4 : CHECK_EQ(geomList.size(), 4);
141 :
142 4 : THEN("Geometry is equal")
143 4 : {
144 1 : CHECK_EQ((desc.getGeometryAt(0)), g1);
145 1 : CHECK_EQ(geomList[0], g1);
146 1 : }
147 4 : THEN("Geometry is equal")
148 4 : {
149 1 : CHECK_EQ((desc.getGeometryAt(1)), g2);
150 1 : CHECK_EQ(geomList[1], g2);
151 1 : }
152 4 : THEN("Geometry is equal")
153 4 : {
154 1 : CHECK_EQ((desc.getGeometryAt(2)), g3);
155 1 : CHECK_EQ(geomList[2], g3);
156 1 : }
157 4 : THEN("Geometry is equal")
158 4 : {
159 1 : CHECK_EQ((desc.getGeometryAt(3)), g4);
160 1 : CHECK_EQ(geomList[3], g4);
161 1 : }
162 4 : }
163 :
164 5 : WHEN("Check for multiple poses, that all the overloads compute the same rays")
165 5 : {
166 4 : for (index_t pose : {0, 1, 2, 3}) {
167 :
168 12 : for (index_t detPixel : {0, 2, 4}) {
169 12 : IndexVector_t pixel(2);
170 12 : pixel << detPixel, pose;
171 :
172 12 : RealVector_t pixelReal(1);
173 12 : pixelReal << static_cast<real_t>(detPixel) + 0.5f;
174 :
175 12 : auto ray1 =
176 12 : desc.computeRayFromDetectorCoord(desc.getIndexFromCoordinate(pixel));
177 12 : auto ray2 = desc.computeRayFromDetectorCoord(pixel);
178 12 : auto ray3 = desc.computeRayFromDetectorCoord(pixelReal, pose);
179 :
180 12 : auto ro1 = ray1.origin();
181 12 : auto rd1 = ray1.direction();
182 :
183 12 : auto ro2 = ray2.origin();
184 12 : auto rd2 = ray2.direction();
185 :
186 12 : auto ro3 = ray3.origin();
187 12 : auto rd3 = ray3.direction();
188 :
189 12 : CHECK_EQ(ro1, ro2);
190 12 : CHECK_EQ(ro1, ro3);
191 :
192 12 : CHECK_EQ(rd1, rd2);
193 12 : CHECK_EQ(rd1, rd3);
194 :
195 : // Shouldn't be necessary, but whatever
196 12 : CHECK_EQ(ro2, ro3);
197 12 : CHECK_EQ(rd2, rd3);
198 12 : }
199 4 : }
200 1 : }
201 5 : }
202 10 : }
203 10 : }
204 :
205 : TEST_CASE("PlanarDetectorDescriptor: Testing 3D PlanarDetectorDescriptor")
206 2 : {
207 2 : GIVEN("Given a 5x5x5 Volume and a single 5x5 wide detector pose")
208 2 : {
209 2 : IndexVector_t volSize(3);
210 2 : volSize << 5, 5, 5;
211 2 : VolumeDescriptor ddVol(volSize);
212 :
213 2 : IndexVector_t sinoSize(3);
214 2 : sinoSize << 5, 5, 1;
215 :
216 2 : real_t s2c = 10;
217 2 : real_t c2d = 4;
218 :
219 2 : Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
220 2 : VolumeData3D{Size3D{volSize}}, SinogramData3D{Size3D{sinoSize}},
221 2 : RotationAngles3D{Gamma{0}});
222 :
223 2 : PlanarDetectorDescriptor desc(sinoSize, {g});
224 :
225 2 : WHEN("Retreiving the single geometry pose")
226 2 : {
227 1 : auto geom = desc.getGeometryAt(0);
228 1 : auto geomList = desc.getGeometry();
229 :
230 1 : CHECK_EQ(desc.getNumberOfGeometryPoses(), 1);
231 1 : CHECK_EQ(geomList.size(), 1);
232 :
233 1 : THEN("Geometry is equal")
234 1 : {
235 1 : CHECK_EQ((geom), g);
236 1 : CHECK_EQ(geomList[0], g);
237 1 : }
238 1 : }
239 :
240 2 : WHEN("Generating rays for detecor pixels 0, 2 and 4 for each dim")
241 2 : {
242 3 : for (index_t detPixel1 : {0, 2, 4}) {
243 9 : for (index_t detPixel2 : {0, 2, 4}) {
244 9 : RealVector_t pixel(2);
245 9 : pixel << static_cast<real_t>(detPixel1) + 0.5f,
246 9 : static_cast<real_t>(detPixel2) + 0.5f;
247 :
248 : // Check that ray for IndexVector_t is equal to previous one
249 9 : auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
250 :
251 : // Create variables, which make typing quicker
252 9 : auto ro = ray.origin();
253 9 : auto rd = ray.direction();
254 :
255 : // Check that ray origin is camera center
256 9 : auto c = g.getCameraCenter();
257 9 : CHECK_EQ((ro - c).sum(), Approx(0));
258 :
259 9 : auto o = ddVol.getLocationOfOrigin();
260 9 : RealVector_t detCoordWorld(3);
261 9 : detCoordWorld << pixel[0] - o[0], pixel[1] - o[1], c2d;
262 9 : RealVector_t rotD = g.getRotationMatrix().transpose() * detCoordWorld + o;
263 :
264 9 : real_t factor = 0;
265 9 : if (std::abs(rd[0]) > 0)
266 6 : factor = (rotD[0] - ro[0]) / rd[0];
267 3 : else if (std::abs(rd[1]) > 0)
268 2 : factor = (rotD[1] - ro[1]) / rd[1];
269 1 : else if (std::abs(rd[2]) > 0)
270 1 : factor = (rotD[2] - ro[2] / rd[2]);
271 :
272 9 : CHECK_EQ((ro[0] + factor * rd[0]), Approx(rotD[0]));
273 9 : CHECK_EQ((ro[1] + factor * rd[1]), Approx(rotD[1]));
274 9 : CHECK_EQ((ro[2] + factor * rd[2]), Approx(rotD[2]));
275 9 : }
276 3 : }
277 1 : }
278 2 : }
279 2 : }
280 :
281 : TEST_SUITE_END();
|
