14 dls(const cv::Mat& opoints, const cv::Mat& ipoints);
17 bool compute_pose(cv::Mat& R, cv::Mat& t);
22 template <typename OpointType, typename IpointType>
23 void init_points(const cv::Mat& opoints, const cv::Mat& ipoints)
25 for(int i = 0; i < N; i++)
27 p.at<double>(0,i) = opoints.at<OpointType>(0,i).x;
28 p.at<double>(1,i) = opoints.at<OpointType>(0,i).y;
29 p.at<double>(2,i) = opoints.at<OpointType>(0,i).z;
31 // compute mean of object points
32 mn.at<double>(0) += p.at<double>(0,i);
33 mn.at<double>(1) += p.at<double>(1,i);
34 mn.at<double>(2) += p.at<double>(2,i);
36 // make z into unit vectors from normalized pixel coords
37 double sr = std::pow(ipoints.at<IpointType>(0,i).x, 2) +
38 std::pow(ipoints.at<IpointType>(0,i).y, 2) + (double)1;
41 z.at<double>(0,i) = ipoints.at<IpointType>(0,i).x / sr;
42 z.at<double>(1,i) = ipoints.at<IpointType>(0,i).y / sr;
43 z.at<double>(2,i) = (double)1 / sr;
46 mn.at<double>(0) /= (double)N;
47 mn.at<double>(1) /= (double)N;
48 mn.at<double>(2) /= (double)N;
52 cv::Mat LeftMultVec(const cv::Mat& v);
53 void run_kernel(const cv::Mat& pp);
54 void build_coeff_matrix(const cv::Mat& pp, cv::Mat& Mtilde, cv::Mat& D);
55 void compute_eigenvec(const cv::Mat& Mtilde, cv::Mat& eigenval_real, cv::Mat& eigenval_imag,
56 cv::Mat& eigenvec_real, cv::Mat& eigenvec_imag);
57 void fill_coeff(const cv::Mat * D);
60 cv::Mat cayley_LS_M(const std::vector<double>& a, const std::vector<double>& b,
61 const std::vector<double>& c, const std::vector<double>& u);
62 cv::Mat Hessian(const double s[]);
63 cv::Mat cayley2rotbar(const cv::Mat& s);
64 cv::Mat skewsymm(const cv::Mat * X1);
67 cv::Mat rotx(const double t);
68 cv::Mat roty(const double t);
69 cv::Mat rotz(const double t);
70 cv::Mat mean(const cv::Mat& M);
71 bool is_empty(const cv::Mat * v);
72 bool positive_eigenvalues(const cv::Mat * eigenvalues);
74 cv::Mat p, z, mn; // object-image points
75 int N; // number of input points
76 std::vector<double> f1coeff, f2coeff, f3coeff, cost_; // coefficient for coefficients matrix
77 std::vector<cv::Mat> C_est_, t_est_; // optimal candidates
78 cv::Mat C_est__, t_est__; // optimal found solution
79 double cost__; // optimal found solution
82 class EigenvalueDecomposition {
85 // Holds the data dimension.
88 // Stores real/imag part of a complex division.
91 // Pointer to internal memory.
95 // Holds the computed eigenvalues.
98 // Holds the computed eigenvectors.
102 template<typename _Tp>
103 _Tp *alloc_1d(int m) {
108 template<typename _Tp>
109 _Tp *alloc_1d(int m, _Tp val) {
110 _Tp *arr = alloc_1d<_Tp> (m);
111 for (int i = 0; i < m; i++)
117 template<typename _Tp>
118 _Tp **alloc_2d(int m, int _n) {
119 _Tp **arr = new _Tp*[m];
120 for (int i = 0; i < m; i++)
121 arr[i] = new _Tp[_n];
126 template<typename _Tp>
127 _Tp **alloc_2d(int m, int _n, _Tp val) {
128 _Tp **arr = alloc_2d<_Tp> (m, _n);
129 for (int i = 0; i < m; i++) {
130 for (int j = 0; j < _n; j++) {
137 void cdiv(double xr, double xi, double yr, double yi) {
139 if (std::abs(yr) > std::abs(yi)) {
142 cdivr = (xr + r * xi) / dv;
143 cdivi = (xi - r * xr) / dv;
147 cdivr = (r * xr + xi) / dv;
148 cdivi = (r * xi - xr) / dv;
152 // Nonsymmetric reduction from Hessenberg to real Schur form.
156 // This is derived from the Algol procedure hqr2,
157 // by Martin and Wilkinson, Handbook for Auto. Comp.,
158 // Vol.ii-Linear Algebra, and the corresponding
159 // Fortran subroutine in EISPACK.
166 double eps = std::pow(2.0, -52.0);
167 double exshift = 0.0;
168 double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
170 // Store roots isolated by balanc and compute matrix norm
173 for (int i = 0; i < nn; i++) {
174 if (i < low || i > high) {
178 for (int j = std::max(i - 1, 0); j < nn; j++) {
179 norm = norm + std::abs(H[i][j]);
183 // Outer loop over eigenvalue index
187 // Look for single small sub-diagonal element
190 s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]);
194 if (std::abs(H[l][l - 1]) < eps * s) {
200 // Check for convergence
204 H[n1][n1] = H[n1][n1] + exshift;
212 } else if (l == n1 - 1) {
213 w = H[n1][n1 - 1] * H[n1 - 1][n1];
214 p = (H[n1 - 1][n1 - 1] - H[n1][n1]) / 2.0;
216 z = std::sqrt(std::abs(q));
217 H[n1][n1] = H[n1][n1] + exshift;
218 H[n1 - 1][n1 - 1] = H[n1 - 1][n1 - 1] + exshift;
237 s = std::abs(x) + std::abs(z);
240 r = std::sqrt(p * p + q * q);
246 for (int j = n1 - 1; j < nn; j++) {
248 H[n1 - 1][j] = q * z + p * H[n1][j];
249 H[n1][j] = q * H[n1][j] - p * z;
252 // Column modification
254 for (int i = 0; i <= n1; i++) {
256 H[i][n1 - 1] = q * z + p * H[i][n1];
257 H[i][n1] = q * H[i][n1] - p * z;
260 // Accumulate transformations
262 for (int i = low; i <= high; i++) {
264 V[i][n1 - 1] = q * z + p * V[i][n1];
265 V[i][n1] = q * V[i][n1] - p * z;
279 // No convergence yet
289 y = H[n1 - 1][n1 - 1];
290 w = H[n1][n1 - 1] * H[n1 - 1][n1];
293 // Wilkinson's original ad hoc shift
297 for (int i = low; i <= n1; i++) {
300 s = std::abs(H[n1][n1 - 1]) + std::abs(H[n1 - 1][n1 - 2]);
305 // MATLAB's new ad hoc shift
315 s = x - w / ((y - x) / 2.0 + s);
316 for (int i = low; i <= n1; i++) {
324 iter = iter + 1; // (Could check iteration count here.)
326 // Look for two consecutive small sub-diagonal elements
332 p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
333 q = H[m + 1][m + 1] - z - r - s;
335 s = std::abs(p) + std::abs(q) + std::abs(r);
342 if (std::abs(H[m][m - 1]) * (std::abs(q) + std::abs(r)) < eps * (std::abs(p)
343 * (std::abs(H[m - 1][m - 1]) + std::abs(z) + std::abs(
344 H[m + 1][m + 1])))) {
350 for (int i = m + 2; i <= n1; i++) {
357 // Double QR step involving rows l:n and columns m:n
359 for (int k = m; k <= n1 - 1; k++) {
360 bool notlast = (k != n1 - 1);
364 r = (notlast ? H[k + 2][k - 1] : 0.0);
365 x = std::abs(p) + std::abs(q) + std::abs(r);
375 s = std::sqrt(p * p + q * q + r * r);
381 H[k][k - 1] = -s * x;
383 H[k][k - 1] = -H[k][k - 1];
394 for (int j = k; j < nn; j++) {
395 p = H[k][j] + q * H[k + 1][j];
397 p = p + r * H[k + 2][j];
398 H[k + 2][j] = H[k + 2][j] - p * z;
400 H[k][j] = H[k][j] - p * x;
401 H[k + 1][j] = H[k + 1][j] - p * y;
404 // Column modification
406 for (int i = 0; i <= std::min(n1, k + 3); i++) {
407 p = x * H[i][k] + y * H[i][k + 1];
409 p = p + z * H[i][k + 2];
410 H[i][k + 2] = H[i][k + 2] - p * r;
412 H[i][k] = H[i][k] - p;
413 H[i][k + 1] = H[i][k + 1] - p * q;
416 // Accumulate transformations
418 for (int i = low; i <= high; i++) {
419 p = x * V[i][k] + y * V[i][k + 1];
421 p = p + z * V[i][k + 2];
422 V[i][k + 2] = V[i][k + 2] - p * r;
424 V[i][k] = V[i][k] - p;
425 V[i][k + 1] = V[i][k + 1] - p * q;
429 } // check convergence
430 } // while (n1 >= low)
432 // Backsubstitute to find vectors of upper triangular form
438 for (n1 = nn - 1; n1 >= 0; n1--) {
447 for (int i = n1 - 1; i >= 0; i--) {
450 for (int j = l; j <= n1; j++) {
451 r = r + H[i][j] * H[j][n1];
462 H[i][n1] = -r / (eps * norm);
465 // Solve real equations
470 q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
471 t = (x * s - z * r) / q;
473 if (std::abs(x) > std::abs(z)) {
474 H[i + 1][n1] = (-r - w * t) / x;
476 H[i + 1][n1] = (-s - y * t) / z;
482 t = std::abs(H[i][n1]);
483 if ((eps * t) * t > 1) {
484 for (int j = i; j <= n1; j++) {
485 H[j][n1] = H[j][n1] / t;
494 // Last vector component imaginary so matrix is triangular
496 if (std::abs(H[n1][n1 - 1]) > std::abs(H[n1 - 1][n1])) {
497 H[n1 - 1][n1 - 1] = q / H[n1][n1 - 1];
498 H[n1 - 1][n1] = -(H[n1][n1] - p) / H[n1][n1 - 1];
500 cdiv(0.0, -H[n1 - 1][n1], H[n1 - 1][n1 - 1] - p, q);
501 H[n1 - 1][n1 - 1] = cdivr;
502 H[n1 - 1][n1] = cdivi;
506 for (int i = n1 - 2; i >= 0; i--) {
507 double ra, sa, vr, vi;
510 for (int j = l; j <= n1; j++) {
511 ra = ra + H[i][j] * H[j][n1 - 1];
512 sa = sa + H[i][j] * H[j][n1];
523 cdiv(-ra, -sa, w, q);
524 H[i][n1 - 1] = cdivr;
528 // Solve complex equations
532 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
533 vi = (d[i] - p) * 2.0 * q;
534 if (vr == 0.0 && vi == 0.0) {
535 vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x)
536 + std::abs(y) + std::abs(z));
538 cdiv(x * r - z * ra + q * sa,
539 x * s - z * sa - q * ra, vr, vi);
540 H[i][n1 - 1] = cdivr;
542 if (std::abs(x) > (std::abs(z) + std::abs(q))) {
543 H[i + 1][n1 - 1] = (-ra - w * H[i][n1 - 1] + q
545 H[i + 1][n1] = (-sa - w * H[i][n1] - q * H[i][n1
548 cdiv(-r - y * H[i][n1 - 1], -s - y * H[i][n1], z,
550 H[i + 1][n1 - 1] = cdivr;
551 H[i + 1][n1] = cdivi;
557 t = std::max(std::abs(H[i][n1 - 1]), std::abs(H[i][n1]));
558 if ((eps * t) * t > 1) {
559 for (int j = i; j <= n1; j++) {
560 H[j][n1 - 1] = H[j][n1 - 1] / t;
561 H[j][n1] = H[j][n1] / t;
569 // Vectors of isolated roots
571 for (int i = 0; i < nn; i++) {
572 if (i < low || i > high) {
573 for (int j = i; j < nn; j++) {
579 // Back transformation to get eigenvectors of original matrix
581 for (int j = nn - 1; j >= low; j--) {
582 for (int i = low; i <= high; i++) {
584 for (int k = low; k <= std::min(j, high); k++) {
585 z = z + V[i][k] * H[k][j];
592 // Nonsymmetric reduction to Hessenberg form.
594 // This is derived from the Algol procedures orthes and ortran,
595 // by Martin and Wilkinson, Handbook for Auto. Comp.,
596 // Vol.ii-Linear Algebra, and the corresponding
597 // Fortran subroutines in EISPACK.
601 for (int m = low + 1; m <= high - 1; m++) {
606 for (int i = m; i <= high; i++) {
607 scale = scale + std::abs(H[i][m - 1]);
611 // Compute Householder transformation.
614 for (int i = high; i >= m; i--) {
615 ort[i] = H[i][m - 1] / scale;
616 h += ort[i] * ort[i];
618 double g = std::sqrt(h);
625 // Apply Householder similarity transformation
626 // H = (I-u*u'/h)*H*(I-u*u')/h)
628 for (int j = m; j < n; j++) {
630 for (int i = high; i >= m; i--) {
631 f += ort[i] * H[i][j];
634 for (int i = m; i <= high; i++) {
635 H[i][j] -= f * ort[i];
639 for (int i = 0; i <= high; i++) {
641 for (int j = high; j >= m; j--) {
642 f += ort[j] * H[i][j];
645 for (int j = m; j <= high; j++) {
646 H[i][j] -= f * ort[j];
649 ort[m] = scale * ort[m];
650 H[m][m - 1] = scale * g;
654 // Accumulate transformations (Algol's ortran).
656 for (int i = 0; i < n; i++) {
657 for (int j = 0; j < n; j++) {
658 V[i][j] = (i == j ? 1.0 : 0.0);
662 for (int m = high - 1; m >= low + 1; m--) {
663 if (H[m][m - 1] != 0.0) {
664 for (int i = m + 1; i <= high; i++) {
665 ort[i] = H[i][m - 1];
667 for (int j = m; j <= high; j++) {
669 for (int i = m; i <= high; i++) {
670 g += ort[i] * V[i][j];
672 // Double division avoids possible underflow
673 g = (g / ort[m]) / H[m][m - 1];
674 for (int i = m; i <= high; i++) {
675 V[i][j] += g * ort[i];
682 // Releases all internal working memory.
684 // releases the working data
688 for (int i = 0; i < n; i++) {
696 // Computes the Eigenvalue Decomposition for a matrix given in H.
698 // Allocate memory for the working data.
699 V = alloc_2d<double> (n, n, 0.0);
700 d = alloc_1d<double> (n);
701 e = alloc_1d<double> (n);
702 ort = alloc_1d<double> (n);
703 // Reduce to Hessenberg form.
705 // Reduce Hessenberg to real Schur form.
707 // Copy eigenvalues to OpenCV Matrix.
708 _eigenvalues.create(1, n, CV_64FC1);
709 for (int i = 0; i < n; i++) {
710 _eigenvalues.at<double> (0, i) = d[i];
712 // Copy eigenvectors to OpenCV Matrix.
713 _eigenvectors.create(n, n, CV_64FC1);
714 for (int i = 0; i < n; i++)
715 for (int j = 0; j < n; j++)
716 _eigenvectors.at<double> (i, j) = V[i][j];
717 // Deallocate the memory by releasing all internal working data.
722 EigenvalueDecomposition()
725 // Initializes & computes the Eigenvalue Decomposition for a general matrix
726 // given in src. This function is a port of the EigenvalueSolver in JAMA,
727 // which has been released to public domain by The MathWorks and the
728 // National Institute of Standards and Technology (NIST).
729 EigenvalueDecomposition(InputArray src) {
733 // This function computes the Eigenvalue Decomposition for a general matrix
734 // given in src. This function is a port of the EigenvalueSolver in JAMA,
735 // which has been released to public domain by The MathWorks and the
736 // National Institute of Standards and Technology (NIST).
737 void compute(InputArray src)
739 /*if(isSymmetric(src)) {
740 // Fall back to OpenCV for a symmetric matrix!
741 cv::eigen(src, _eigenvalues, _eigenvectors);
744 // Convert the given input matrix to double. Is there any way to
745 // prevent allocating the temporary memory? Only used for copying
746 // into working memory and deallocated after.
747 src.getMat().convertTo(tmp, CV_64FC1);
748 // Get dimension of the matrix.
750 // Allocate the matrix data to work on.
751 this->H = alloc_2d<double> (n, n);
752 // Now safely copy the data.
753 for (int i = 0; i < tmp.rows; i++) {
754 for (int j = 0; j < tmp.cols; j++) {
755 this->H[i][j] = tmp.at<double>(i, j);
758 // Deallocates the temporary matrix before computing.
760 // Performs the eigenvalue decomposition of H.
765 ~EigenvalueDecomposition() {}
767 // Returns the eigenvalues of the Eigenvalue Decomposition.
768 Mat eigenvalues() { return _eigenvalues; }
769 // Returns the eigenvectors of the Eigenvalue Decomposition.
770 Mat eigenvectors() { return _eigenvectors; }