Mercurial > hg > orthanc-stone

comparison UnitTestsSources/UnitTestsMain.cpp @ 162:77715c340767 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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)