16 #include <mrpt/config.h> 17 #if MRPT_HAS_OPENCV && MRPT_OPENCV_VERSION_NUM<0x240 18 # undef MRPT_HAS_OPENCV 19 # define MRPT_HAS_OPENCV 0 35 fx = cam_intrinsic(0,0);
36 fy = cam_intrinsic(1,1);
37 cx = cam_intrinsic(0,2);
38 cy = cam_intrinsic(1,2);
39 init_inverse_parameters();
48 init_inverse_parameters();
54 if (cameraMatrix.depth() == CV_32F)
55 init_camera_parameters<float>(cameraMatrix);
57 init_camera_parameters<double>(cameraMatrix);
58 init_inverse_parameters();
63 double rotation_matrix[3][3], translation[3];
64 std::vector<double>
points;
65 if (opoints.depth() == ipoints.depth())
67 if (opoints.depth() == CV_32F)
68 extract_points<cv::Point3f,cv::Point2f>(opoints, ipoints,
points);
70 extract_points<cv::Point3d,cv::Point2d>(opoints, ipoints,
points);
72 else if (opoints.depth() == CV_32F)
73 extract_points<cv::Point3f,cv::Point2d>(opoints, ipoints,
points);
75 extract_points<cv::Point3d,cv::Point2f>(opoints, ipoints,
points);
80 cv::Mat(3, 1, CV_64F, translation).copyTo(tvec);
81 cv::Mat(3, 3, CV_64F, rotation_matrix).copyTo(
R);
89 double rotation_matrix[3][3], translation[3];
90 std::vector<double>
points;
91 extract_points(obj_pts, img_pts,
points);
95 R = Eigen::Map<Eigen::Matrix3d>((
double*)rotation_matrix);
96 t = Eigen::Map<Eigen::Vector3d>(translation);
101 double mu0,
double mv0,
double X0,
double Y0,
double Z0,
102 double mu1,
double mv1,
double X1,
double Y1,
double Z1,
103 double mu2,
double mv2,
double X2,
double Y2,
double Z2,
104 double mu3,
double mv3,
double X3,
double Y3,
double Z3)
106 double Rs[4][3][3], ts[4][3];
108 int n = solve(Rs, ts, mu0, mv0, X0, Y0, Z0, mu1, mv1, X1, Y1, Z1, mu2, mv2, X2, Y2, Z2);
114 double min_reproj = 0;
115 for(
int i = 0; i <
n; i++) {
116 double X3p = Rs[i][0][0] * X3 + Rs[i][0][1] * Y3 + Rs[i][0][2] * Z3 + ts[i][0];
117 double Y3p = Rs[i][1][0] * X3 + Rs[i][1][1] * Y3 + Rs[i][1][2] * Z3 + ts[i][1];
118 double Z3p = Rs[i][2][0] * X3 + Rs[i][2][1] * Y3 + Rs[i][2][2] * Z3 + ts[i][2];
119 double mu3p = cx + fx * X3p / Z3p;
120 double mv3p = cy + fy * Y3p / Z3p;
121 double reproj = (mu3p - mu3) * (mu3p - mu3) + (mv3p - mv3) * (mv3p - mv3);
122 if (i == 0 || min_reproj > reproj) {
128 for(
int i = 0; i < 3; i++) {
129 for(
int j = 0; j < 3; j++)
130 R[i][j] = Rs[ns][i][j];
138 double mu0,
double mv0,
double X0,
double Y0,
double Z0,
139 double mu1,
double mv1,
double X1,
double Y1,
double Z1,
140 double mu2,
double mv2,
double X2,
double Y2,
double Z2)
142 double mk0, mk1, mk2;
145 mu0 = inv_fx * mu0 - cx_fx;
146 mv0 = inv_fy * mv0 - cy_fy;
147 norm = sqrt(mu0 * mu0 + mv0 * mv0 + 1);
148 mk0 = 1. /
norm; mu0 *= mk0; mv0 *= mk0;
150 mu1 = inv_fx * mu1 - cx_fx;
151 mv1 = inv_fy * mv1 - cy_fy;
152 norm = sqrt(mu1 * mu1 + mv1 * mv1 + 1);
153 mk1 = 1. /
norm; mu1 *= mk1; mv1 *= mk1;
155 mu2 = inv_fx * mu2 - cx_fx;
156 mv2 = inv_fy * mv2 - cy_fy;
157 norm = sqrt(mu2 * mu2 + mv2 * mv2 + 1);
158 mk2 = 1. /
norm; mu2 *= mk2; mv2 *= mk2;
161 distances[0] = sqrt( (X1 - X2) * (X1 - X2) + (Y1 - Y2) * (Y1 - Y2) + (Z1 - Z2) * (Z1 - Z2) );
162 distances[1] = sqrt( (X0 - X2) * (X0 - X2) + (Y0 - Y2) * (Y0 - Y2) + (Z0 - Z2) * (Z0 - Z2) );
163 distances[2] = sqrt( (X0 - X1) * (X0 - X1) + (Y0 - Y1) * (Y0 - Y1) + (Z0 - Z1) * (Z0 - Z1) );
167 cosines[0] = mu1 * mu2 + mv1 * mv2 + mk1 * mk2;
168 cosines[1] = mu0 * mu2 + mv0 * mv2 + mk0 * mk2;
169 cosines[2] = mu0 * mu1 + mv0 * mv1 + mk0 * mk1;
171 double lengths[4][3];
172 int n = solve_for_lengths(lengths, distances, cosines);
174 int nb_solutions = 0;
175 for(
int i = 0; i <
n; i++) {
178 M_orig[0][0] = lengths[i][0] * mu0;
179 M_orig[0][1] = lengths[i][0] * mv0;
180 M_orig[0][2] = lengths[i][0] * mk0;
182 M_orig[1][0] = lengths[i][1] * mu1;
183 M_orig[1][1] = lengths[i][1] * mv1;
184 M_orig[1][2] = lengths[i][1] * mk1;
186 M_orig[2][0] = lengths[i][2] * mu2;
187 M_orig[2][1] = lengths[i][2] * mv2;
188 M_orig[2][2] = lengths[i][2] * mk2;
190 if (!align(M_orig, X0, Y0, Z0, X1, Y1, Z1, X2, Y2, Z2,
R[nb_solutions],
t[nb_solutions]))
213 double p = cosines[0] * 2;
214 double q = cosines[1] * 2;
215 double r = cosines[2] * 2;
217 double inv_d22 = 1. / (distances[2] * distances[2]);
218 double a = inv_d22 * (distances[0] * distances[0]);
219 double b = inv_d22 * (distances[1] * distances[1]);
221 double a2 =
a *
a,
b2 =
b *
b, p2 =
p *
p, q2 =
q *
q, r2 =
r *
r;
222 double pr =
p *
r, pqr =
q * pr;
225 if (p2 + q2 + r2 - pqr - 1 == 0)
228 double ab =
a *
b, a_2 = 2*
a;
230 double A = -2 *
b +
b2 +
a2 + 1 + ab*(2 - r2) - a_2;
233 if (A == 0)
return 0;
237 double B =
q*(-2*(ab +
a2 + 1 -
b) + r2*ab + a_4) + pr*(
b -
b2 + ab);
238 double C = q2 +
b2*(r2 + p2 - 2) -
b*(p2 + pqr) - ab*(r2 + pqr) + (
a2 - a_2)*(2 + q2) + 2;
239 double D = pr*(ab-
b2+
b) +
q*((p2-2)*
b + 2 * (ab -
a2) + a_4 - 2);
240 double E = 1 + 2*(
b -
a - ab) +
b2 -
b*p2 +
a2;
242 double temp = (p2*(
a-1+
b) + r2*(
a-1-
b) + pqr -
a*pqr);
243 double b0 =
b * temp * temp;
248 double real_roots[4];
249 int n =
solve_deg4(A, B, C, D, E, real_roots[0], real_roots[1], real_roots[2], real_roots[3]);
254 int nb_solutions = 0;
255 double r3 = r2*
r, pr2 =
p*r2, r3q = r3 *
q;
256 double inv_b0 = 1. /
b0;
259 for(
int i = 0; i <
n; i++) {
260 double x = real_roots[i];
269 ((1-
a-
b)*x2 + (
q*
a-
q)*
x + 1 -
a +
b) *
270 (((r3*(
a2 + ab*(2 - r2) - a_2 +
b2 - 2*
b + 1)) *
x +
272 (r3q*(2*(
b-
a2) + a_4 + ab*(r2 - 2) - 2) + pr2*(1 +
a2 + 2*(ab-
a-
b) + r2*(
b -
b2) +
b2))) * x2 +
274 (r3*(q2*(1-2*
a+
a2) + r2*(
b2-ab) - a_4 + 2*(
a2 -
b2) + 2) +
r*p2*(
b2 + 2*(ab -
b -
a) + 1 +
a2) + pr2*
q*(a_4 + 2*(
b - ab -
a2) - 2 - r2*
b)) *
x +
276 2*r3q*(a_2 -
b -
a2 + ab - 1) + pr2*(q2 - a_4 + 2*(
a2 -
b2) + r2*
b + q2*(
a2 - a_2) + 2) +
277 p2*(
p*(2*(ab -
a -
b) +
a2 +
b2 + 1) + 2*
q*
r*(
b + a_2 -
a2 - ab - 1)));
283 double y = inv_b0 *
b1;
284 double v = x2 +
y*
y -
x*
y*
r;
289 double Z = distances[2] / sqrt(
v);
293 lengths[nb_solutions][0] = X;
294 lengths[nb_solutions][1] = Y;
295 lengths[nb_solutions][2] = Z;
304 double X0,
double Y0,
double Z0,
305 double X1,
double Y1,
double Z1,
306 double X2,
double Y2,
double Z2,
307 double R[3][3],
double T[3])
310 double C_start[3], C_end[3];
311 for(
int i = 0; i < 3; i++) C_end[i] = (M_end[0][i] + M_end[1][i] + M_end[2][i]) / 3;
312 C_start[0] = (X0 + X1 + X2) / 3;
313 C_start[1] = (Y0 + Y1 + Y2) / 3;
314 C_start[2] = (Z0 + Z1 + Z2) / 3;
318 for(
int j = 0; j < 3; j++) {
319 s[0 * 3 + j] = (X0 * M_end[0][j] + X1 * M_end[1][j] + X2 * M_end[2][j]) / 3 - C_end[j] * C_start[0];
320 s[1 * 3 + j] = (Y0 * M_end[0][j] + Y1 * M_end[1][j] + Y2 * M_end[2][j]) / 3 - C_end[j] * C_start[1];
321 s[2 * 3 + j] = (Z0 * M_end[0][j] + Z1 * M_end[1][j] + Z2 * M_end[2][j]) / 3 - C_end[j] * C_start[2];
324 double Qs[16], evs[4], U[16];
326 Qs[0 * 4 + 0] =
s[0 * 3 + 0] +
s[1 * 3 + 1] +
s[2 * 3 + 2];
327 Qs[1 * 4 + 1] =
s[0 * 3 + 0] -
s[1 * 3 + 1] -
s[2 * 3 + 2];
328 Qs[2 * 4 + 2] =
s[1 * 3 + 1] -
s[2 * 3 + 2] -
s[0 * 3 + 0];
329 Qs[3 * 4 + 3] =
s[2 * 3 + 2] -
s[0 * 3 + 0] -
s[1 * 3 + 1];
331 Qs[1 * 4 + 0] = Qs[0 * 4 + 1] =
s[1 * 3 + 2] -
s[2 * 3 + 1];
332 Qs[2 * 4 + 0] = Qs[0 * 4 + 2] =
s[2 * 3 + 0] -
s[0 * 3 + 2];
333 Qs[3 * 4 + 0] = Qs[0 * 4 + 3] =
s[0 * 3 + 1] -
s[1 * 3 + 0];
334 Qs[2 * 4 + 1] = Qs[1 * 4 + 2] =
s[1 * 3 + 0] +
s[0 * 3 + 1];
335 Qs[3 * 4 + 1] = Qs[1 * 4 + 3] =
s[2 * 3 + 0] +
s[0 * 3 + 2];
336 Qs[3 * 4 + 2] = Qs[2 * 4 + 3] =
s[2 * 3 + 1] +
s[1 * 3 + 2];
338 jacobi_4x4(Qs, evs, U);
342 double ev_max = evs[i_ev];
343 for(
int i = 1; i < 4; i++)
345 ev_max = evs[i_ev = i];
349 for(
int i = 0; i < 4; i++)
350 q[i] = U[i * 4 + i_ev];
352 double q02 =
q[0] *
q[0], q12 =
q[1] *
q[1], q22 =
q[2] *
q[2], q32 =
q[3] *
q[3];
353 double q0_1 =
q[0] *
q[1], q0_2 =
q[0] *
q[2], q0_3 =
q[0] *
q[3];
354 double q1_2 =
q[1] *
q[2], q1_3 =
q[1] *
q[3];
355 double q2_3 =
q[2] *
q[3];
357 R[0][0] = q02 + q12 - q22 - q32;
358 R[0][1] = 2. * (q1_2 - q0_3);
359 R[0][2] = 2. * (q1_3 + q0_2);
361 R[1][0] = 2. * (q1_2 + q0_3);
362 R[1][1] = q02 + q22 - q12 - q32;
363 R[1][2] = 2. * (q2_3 - q0_1);
365 R[2][0] = 2. * (q1_3 - q0_2);
366 R[2][1] = 2. * (q2_3 + q0_1);
367 R[2][2] = q02 + q32 - q12 - q22;
369 for(
int i = 0; i < 3; i++)
370 T[i] = C_end[i] - (
R[i][0] * C_start[0] +
R[i][1] * C_start[1] +
R[i][2] * C_start[2]);
378 double Id[16] = {1., 0., 0., 0.,
383 memcpy(U, Id, 16 *
sizeof(
double));
385 B[0] = A[0]; B[1] = A[5]; B[2] = A[10]; B[3] = A[15];
386 memcpy(D, B, 4 *
sizeof(
double));
387 memset(Z, 0, 4 *
sizeof(
double));
389 for(
int iter = 0; iter < 50; iter++) {
390 double sum = fabs(A[1]) + fabs(A[2]) + fabs(A[3]) + fabs(A[6]) + fabs(A[7]) + fabs(A[11]);
395 double tresh = (iter < 3) ? 0.2 *
sum / 16. : 0.0;
396 for(
int i = 0; i < 3; i++) {
397 double * pAij = A + 5 * i + 1;
398 for(
int j = i + 1 ; j < 4; j++) {
400 double eps_machine = 100.0 * fabs(Aij);
402 if ( iter > 3 && fabs(D[i]) + eps_machine == fabs(D[i]) && fabs(D[j]) + eps_machine == fabs(D[j]) )
404 else if (fabs(Aij) > tresh) {
405 double hh = D[j] - D[i],
t;
406 if (fabs(hh) + eps_machine == fabs(hh))
409 double theta = 0.5 * hh / Aij;
410 t = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta));
411 if (theta < 0.0)
t = -
t;
421 double c = 1.0 / sqrt(1 +
t *
t);
423 double tau =
s / (1.0 +
c);
424 for(
int k = 0; k <= i - 1; k++) {
425 double g = A[k * 4 + i], h = A[k * 4 + j];
426 A[k * 4 + i] =
g -
s * (h +
g * tau);
427 A[k * 4 + j] = h +
s * (
g - h * tau);
429 for(
int k = i + 1; k <= j - 1; k++) {
430 double g = A[i * 4 + k], h = A[k * 4 + j];
431 A[i * 4 + k] =
g -
s * (h +
g * tau);
432 A[k * 4 + j] = h +
s * (
g - h * tau);
434 for(
int k = j + 1; k < 4; k++) {
435 double g = A[i * 4 + k], h = A[j * 4 + k];
436 A[i * 4 + k] =
g -
s * (h +
g * tau);
437 A[j * 4 + k] = h +
s * (
g - h * tau);
439 for(
int k = 0; k < 4; k++) {
440 double g = U[k * 4 + i], h = U[k * 4 + j];
441 U[k * 4 + i] =
g -
s * (h +
g * tau);
442 U[k * 4 + j] = h +
s * (
g - h * tau);
449 for(
int i = 0; i < 4; i++) B[i] += Z[i];
450 memcpy(D, B, 4 *
sizeof(
double));
451 memset(Z, 0, 4 *
sizeof(
double));
void BASE_IMPEXP memcpy(void *dest, size_t destSize, const void *src, size_t copyCount) MRPT_NO_THROWS
An OS and compiler independent version of "memcpy".
GLdouble GLdouble GLdouble GLdouble q
bool jacobi_4x4(double *A, double *D, double *U)
Function used to compute the SVD.
GLsizei const GLfloat * points
double inv_fy
Inverse of focal length x.
bool align(double M_start[3][3], double X0, double Y0, double Z0, double X1, double Y1, double Z1, double X2, double Y2, double Z2, double R[3][3], double T[3])
int solve_for_lengths(double lengths[4][3], double distances[3], double cosines[3])
Helper function to solve()
CONTAINER::Scalar sum(const CONTAINER &v)
Computes the sum of all the elements.
double cy_fy
Inverse of image center point x.
int solve_deg4(double a, double b, double c, double d, double e, double &x0, double &x1, double &x2, double &x3)
Reference : Eric W.
void init_inverse_parameters()
Function to compute inverse parameters of camera intrinsic matrix.
bool solve(Eigen::Ref< Eigen::Matrix3d > R, Eigen::Ref< Eigen::Vector3d > t, Eigen::MatrixXd obj_pts, Eigen::MatrixXd img_pts)
p3p(double fx, double fy, double cx, double cy)
Constructor for p3p class using C.
GLdouble GLdouble GLdouble r
double cx_fx
Inverse of focal length y.
GLubyte GLubyte GLubyte a
double inv_fx
Image center y.
CONTAINER::Scalar norm(const CONTAINER &v)
P3P Pose estimation Algorithm - Eigen Implementation.