Mercurial > hg > orthanc-stone
comparison UnitTestsSources/UnitTestsMain.cpp @ 162:77715c340767 wasm
unit testing FiniteProjectiveCamera
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 14 Feb 2018 12:24:56 +0100 |
parents | 197a5ddaf68c |
children | 8c5b24892ed2 |
comparison
equal
deleted
inserted
replaced
161:197a5ddaf68c | 162:77715c340767 |
---|---|
145 91, 210, 162, 95, | 145 91, 210, 162, 95, |
146 51, 190, 80, 92)); | 146 51, 190, 80, 92)); |
147 } | 147 } |
148 | 148 |
149 | 149 |
150 TEST(FiniteProjectiveCamera, Decomposition) | 150 static bool CompareMatrix(const OrthancStone::Matrix& a, |
151 const OrthancStone::Matrix& b, | |
152 double threshold = 0.00000001) | |
153 { | |
154 if (a.size1() != b.size1() || | |
155 a.size2() != b.size2()) | |
156 { | |
157 return false; | |
158 } | |
159 | |
160 for (size_t i = 0; i < a.size1(); i++) | |
161 { | |
162 for (size_t j = 0; j < a.size2(); j++) | |
163 { | |
164 if (fabs(a(i, j) - b(i, j)) > threshold) | |
165 { | |
166 LOG(ERROR) << "Too large difference in component (" | |
167 << i << "," << j << "): " << a(i,j) << " != " << b(i,j); | |
168 return false; | |
169 } | |
170 } | |
171 } | |
172 | |
173 return true; | |
174 } | |
175 | |
176 | |
177 static bool CompareVector(const OrthancStone::Vector& a, | |
178 const OrthancStone::Vector& b, | |
179 double threshold = 0.00000001) | |
180 { | |
181 if (a.size() != b.size()) | |
182 { | |
183 return false; | |
184 } | |
185 | |
186 for (size_t i = 0; i < a.size(); i++) | |
187 { | |
188 if (fabs(a(i) - b(i)) > threshold) | |
189 { | |
190 LOG(ERROR) << "Too large difference in component " | |
191 << i << ": " << a(i) << " != " << b(i); | |
192 return false; | |
193 } | |
194 } | |
195 | |
196 return true; | |
197 } | |
198 | |
199 | |
200 | |
201 TEST(FiniteProjectiveCamera, Decomposition1) | |
151 { | 202 { |
152 // Example 6.2 of "Multiple View Geometry in Computer Vision - 2nd | 203 // Example 6.2 of "Multiple View Geometry in Computer Vision - 2nd |
153 // edition" (page 163) | 204 // edition" (page 163) |
154 const double p[12] = { | 205 const double p[12] = { |
155 3.53553e+2, 3.39645e+2, 2.77744e+2, -1.44946e+6, | 206 3.53553e+2, 3.39645e+2, 2.77744e+2, -1.44946e+6, |
194 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | 245 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); |
195 | 246 |
196 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), | 247 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), |
197 camera.GetRotation(), | 248 camera.GetRotation(), |
198 camera.GetCenter()); | 249 camera.GetCenter()); |
199 ASSERT_EQ(3u, camera2.GetMatrix().size1()); | 250 |
200 ASSERT_EQ(4u, camera2.GetMatrix().size2()); | 251 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); |
201 ASSERT_EQ(3u, camera2.GetIntrinsicParameters().size1()); | 252 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); |
202 ASSERT_EQ(3u, camera2.GetIntrinsicParameters().size2()); | 253 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); |
203 ASSERT_EQ(3u, camera2.GetRotation().size1()); | 254 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); |
204 ASSERT_EQ(3u, camera2.GetRotation().size2()); | 255 } |
205 ASSERT_EQ(3u, camera2.GetCenter().size()); | 256 |
206 | 257 |
207 for (size_t i = 0; i < 3; i++) | 258 TEST(FiniteProjectiveCamera, Decomposition2) |
208 { | 259 { |
209 for (size_t j = 0; j < 4; j++) | 260 const double p[] = { 1188.111986, 580.205341, -808.445330, 128000.000000, -366.466264, 1446.510501, 418.499736, 128000.000000, -0.487118, 0.291726, -0.823172, 500.000000 }; |
210 { | 261 const double k[] = { -1528.494743, 0.000000, 256.000000, 0.000000, 1528.494743, 256.000000, 0.000000, 0.000000, 1.000000 }; |
211 ASSERT_NEAR(camera.GetMatrix() (i, j), | 262 const double r[] = { -0.858893, -0.330733, 0.391047, -0.158171, 0.897503, 0.411668, -0.487118, 0.291726, -0.823172 }; |
212 camera2.GetMatrix() (i, j), 0.00000001); | 263 const double c[] = { 243.558936, -145.863085, 411.585964 }; |
213 } | 264 |
265 OrthancStone::FiniteProjectiveCamera camera(p); | |
266 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
267 | |
268 OrthancStone::FiniteProjectiveCamera camera2(k, r, c); | |
269 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix(), 1)); | |
270 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters(), 0.001)); | |
271 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation(), 0.000001)); | |
272 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter(), 0.0001)); | |
273 } | |
274 | |
275 | |
276 TEST(FiniteProjectiveCamera, Decomposition3) | |
277 { | |
278 const double p[] = { 10, 0, 0, 0, | |
279 0, 20, 0, 0, | |
280 0, 0, 30, 0 }; | |
281 | |
282 OrthancStone::FiniteProjectiveCamera camera(p); | |
283 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
284 ASSERT_DOUBLE_EQ(10, camera.GetIntrinsicParameters() (0, 0)); | |
285 ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); | |
286 ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); | |
287 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); | |
288 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); | |
289 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); | |
290 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); | |
291 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); | |
292 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); | |
293 } | |
294 | |
295 | |
296 TEST(FiniteProjectiveCamera, Decomposition4) | |
297 { | |
298 const double p[] = { 1, 0, 0, 10, | |
299 0, 1, 0, 20, | |
300 0, 0, 1, 30 }; | |
301 | |
302 OrthancStone::FiniteProjectiveCamera camera(p); | |
303 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
304 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (0, 0)); | |
305 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (1, 1)); | |
306 ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (2, 2)); | |
307 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); | |
308 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); | |
309 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); | |
310 ASSERT_DOUBLE_EQ(-10, camera.GetCenter() (0)); | |
311 ASSERT_DOUBLE_EQ(-20, camera.GetCenter() (1)); | |
312 ASSERT_DOUBLE_EQ(-30, camera.GetCenter() (2)); | |
313 } | |
314 | |
315 | |
316 TEST(FiniteProjectiveCamera, Decomposition5) | |
317 { | |
318 const double p[] = { 0, 0, 10, 0, | |
319 0, 20, 0, 0, | |
320 30, 0, 0, 0 }; | |
321 | |
322 OrthancStone::FiniteProjectiveCamera camera(p); | |
323 ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); | |
324 ASSERT_DOUBLE_EQ(-10, camera.GetIntrinsicParameters() (0, 0)); | |
325 ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); | |
326 ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); | |
327 ASSERT_DOUBLE_EQ(-1, camera.GetRotation() (0, 2)); | |
328 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); | |
329 ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 0)); | |
330 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); | |
331 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); | |
332 ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); | |
333 | |
334 OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), | |
335 camera.GetRotation(), | |
336 camera.GetCenter()); | |
337 ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); | |
338 ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); | |
339 ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); | |
340 ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); | |
341 } | |
342 | |
343 | |
344 static double GetCosAngle(const OrthancStone::Vector& a, | |
345 const OrthancStone::Vector& b) | |
346 { | |
347 // Returns the cosine of the angle between two vectors | |
348 // https://en.wikipedia.org/wiki/Dot_product#Geometric_definition | |
349 return boost::numeric::ublas::inner_prod(a, b) / | |
350 (boost::numeric::ublas::norm_2(a) * boost::numeric::ublas::norm_2(b)); | |
351 } | |
352 | |
353 | |
354 TEST(FiniteProjectiveCamera, Ray) | |
355 { | |
356 const double pp[] = { -1499.650894, 2954.618773, -259.737419, 637891.819097, | |
357 -2951.517707, -1501.019129, -285.785281, 637891.819097, | |
358 0.008528, 0.003067, -0.999959, 2491.764918 }; | |
359 | |
360 const OrthancStone::FiniteProjectiveCamera camera(pp); | |
361 | |
362 ASSERT_NEAR(-21.2492, camera.GetCenter() (0), 0.0001); | |
363 ASSERT_NEAR(-7.64234, camera.GetCenter() (1), 0.00001); | |
364 ASSERT_NEAR(2491.66, camera.GetCenter() (2), 0.01); | |
365 | |
366 // Image plane that led to these parameters, with principal point at | |
367 // (256,256). The image has dimensions 512x512. | |
368 OrthancStone::Vector o, ax, ay; | |
369 OrthancStone::LinearAlgebra::AssignVector(o, 7.009620, 2.521030, -821.942000); | |
370 OrthancStone::LinearAlgebra::AssignVector(ax, -0.453219, 0.891399, -0.001131 ); | |
371 OrthancStone::LinearAlgebra::AssignVector(ay, -0.891359, -0.453210, -0.008992); | |
372 OrthancStone::CoordinateSystem3D imagePlane(o, ax, ay); | |
373 | |
374 // Back-projection of the principal point | |
375 { | |
376 OrthancStone::Vector ray = camera.GetRayDirection(256, 256); | |
377 | |
378 // The principal axis vector is orthogonal to the image plane | |
379 // (i.e. parallel to the plane normal), in the opposite direction | |
380 // ("-1" corresponds to "cos(pi)"). | |
381 ASSERT_NEAR(-1, GetCosAngle(ray, imagePlane.GetNormal()), 0.0000001); | |
382 | |
383 // Forward projection of principal axis, resulting in the principal point | |
384 double x, y; | |
385 camera.ApplyFinite(x, y, camera.GetCenter() - ray); | |
386 | |
387 ASSERT_NEAR(256, x, 0.00001); | |
388 ASSERT_NEAR(256, y, 0.00001); | |
389 } | |
390 | |
391 // Back-projection of the 4 corners of the image | |
392 std::vector<double> cx, cy; | |
393 cx.push_back(0); | |
394 cy.push_back(0); | |
395 cx.push_back(512); | |
396 cy.push_back(0); | |
397 cx.push_back(512); | |
398 cy.push_back(512); | |
399 cx.push_back(0); | |
400 cy.push_back(512); | |
401 | |
402 bool first = true; | |
403 double angle; | |
404 | |
405 for (size_t i = 0; i < cx.size(); i++) | |
406 { | |
407 OrthancStone::Vector ray = camera.GetRayDirection(cx[i], cy[i]); | |
408 | |
409 // Check that the angle wrt. principal axis is the same for all | |
410 // the 4 corners | |
411 double a = GetCosAngle(ray, imagePlane.GetNormal()); | |
412 if (first) | |
413 { | |
414 first = false; | |
415 angle = a; | |
416 } | |
417 else | |
418 { | |
419 ASSERT_NEAR(angle, a, 0.000001); | |
420 } | |
421 | |
422 // Forward projection of the ray, going back to the original point | |
423 double x, y; | |
424 camera.ApplyFinite(x, y, camera.GetCenter() - ray); | |
425 | |
426 ASSERT_NEAR(cx[i], x, 0.00001); | |
427 ASSERT_NEAR(cy[i], y, 0.00001); | |
428 | |
429 // Alternative construction, by computing the intersection of the | |
430 // ray with the image plane | |
431 OrthancStone::Vector p; | |
432 ASSERT_TRUE(imagePlane.IntersectLine(p, camera.GetCenter(), -ray)); | |
433 imagePlane.ProjectPoint(x, y, p); | |
434 ASSERT_NEAR(cx[i], x + 256, 0.01); | |
435 ASSERT_NEAR(cy[i], y + 256, 0.01); | |
214 } | 436 } |
215 } | 437 } |
216 | 438 |
217 | 439 |
218 int main(int argc, char **argv) | 440 int main(int argc, char **argv) |