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