39 dls(
const cv::Mat& opoints,
const cv::Mat& ipoints);
43 bool compute_pose(cv::Mat& R, cv::Mat&
t);
52 template <
typename Opo
intType,
typename Ipo
intType>
53 void init_points(
const cv::Mat& opoints,
const cv::Mat& ipoints)
55 for(
int i = 0; i < N; i++)
57 p.at<
double>(0,i) = opoints.at<OpointType>(i).x;
58 p.at<
double>(1,i) = opoints.at<OpointType>(i).y;
59 p.at<
double>(2,i) = opoints.at<OpointType>(i).z;
62 mn.at<
double>(0) += p.at<
double>(0,i);
63 mn.at<
double>(1) += p.at<
double>(1,i);
64 mn.at<
double>(2) += p.at<
double>(2,i);
67 double sr =
std::pow(ipoints.at<IpointType>(i).x, 2) +
68 std::pow(ipoints.at<IpointType>(i).y, 2) + (double)1;
71 z.at<
double>(0,i) = ipoints.at<IpointType>(i).x / sr;
72 z.at<
double>(1,i) = ipoints.at<IpointType>(i).y / sr;
73 z.at<
double>(2,i) = (
double)1 / sr;
76 mn.at<
double>(0) /= (
double)N;
77 mn.at<
double>(1) /= (
double)N;
78 mn.at<
double>(2) /= (
double)N;
86 cv::Mat LeftMultVec(
const cv::Mat& v);
92 void run_kernel(
const cv::Mat& pp);
100 void build_coeff_matrix(
const cv::Mat& pp, cv::Mat& Mtilde, cv::Mat& D);
110 void compute_eigenvec(
const cv::Mat& Mtilde, cv::Mat& eigenval_real, cv::Mat& eigenval_imag,
111 cv::Mat& eigenvec_real, cv::Mat& eigenvec_imag);
117 void fill_coeff(
const cv::Mat * D);
127 cv::Mat cayley_LS_M(
const std::vector<double>& a,
const std::vector<double>& b,
128 const std::vector<double>& c,
const std::vector<double>& u);
135 cv::Mat Hessian(
const double s[]);
142 cv::Mat cayley2rotbar(
const cv::Mat& s);
149 cv::Mat skewsymm(
const cv::Mat * X1);
156 cv::Mat rotx(
const double t);
163 cv::Mat roty(
const double t);
170 cv::Mat rotz(
const double t);
177 cv::Mat
mean(
const cv::Mat& M);
184 bool is_empty(
const cv::Mat * v);
191 bool positive_eigenvalues(
const cv::Mat * eigenvalues);
195 std::vector<double> f1coeff, f2coeff, f3coeff, cost_;
196 std::vector<cv::Mat> C_est_, t_est_;
197 cv::Mat C_est__, t_est__;
205 class EigenvalueDecomposition {
217 cv::Mat _eigenvalues;
219 cv::Mat _eigenvectors;
226 template<
typename _Tp>
227 _Tp *alloc_1d(
int m) {
237 template<
typename _Tp>
238 _Tp *alloc_1d(
int m, _Tp val) {
239 _Tp *arr = alloc_1d<_Tp> (m);
240 for (
int i = 0; i < m; i++)
251 template<
typename _Tp>
252 _Tp **alloc_2d(
int m,
int _n) {
253 _Tp **arr =
new _Tp*[m];
254 for (
int i = 0; i < m; i++)
255 arr[i] =
new _Tp[_n];
266 template<
typename _Tp>
267 _Tp **alloc_2d(
int m,
int _n, _Tp val) {
268 _Tp **arr = alloc_2d<_Tp> (m, _n);
269 for (
int i = 0; i < m; i++) {
270 for (
int j = 0; j < _n; j++) {
284 void cdiv(
double xr,
double xi,
double yr,
double yi) {
286 if (std::abs(yr) > std::abs(yi)) {
289 cdivr = (xr + r * xi) / dv;
290 cdivi = (xi - r * xr) / dv;
294 cdivr = (r * xr + xi) / dv;
295 cdivi = (r * xi - xr) / dv;
314 double exshift = 0.0;
315 double p = 0, q = 0, r = 0, s = 0,
z = 0,
t, w,
x, y;
320 for (
int i = 0; i < nn; i++) {
321 if (i < low || i > high) {
325 for (
int j =
std::max(i - 1, 0); j < nn; j++) {
326 norm = norm + std::abs(H[i][j]);
337 s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]);
341 if (std::abs(H[l][l - 1]) < eps * s) {
351 H[n1][n1] = H[n1][n1] + exshift;
359 }
else if (l == n1 - 1) {
360 w = H[n1][n1 - 1] * H[n1 - 1][n1];
361 p = (H[n1 - 1][n1 - 1] - H[n1][n1]) / 2.0;
363 z = std::sqrt(std::abs(q));
364 H[n1][n1] = H[n1][n1] + exshift;
365 H[n1 - 1][n1 - 1] = H[n1 - 1][n1 - 1] + exshift;
384 s = std::abs(x) + std::abs(z);
387 r = std::sqrt(p * p + q * q);
393 for (
int j = n1 - 1; j < nn; j++) {
395 H[n1 - 1][j] = q * z + p * H[n1][j];
396 H[n1][j] = q * H[n1][j] - p *
z;
401 for (
int i = 0; i <= n1; i++) {
403 H[i][n1 - 1] = q * z + p * H[i][n1];
404 H[i][n1] = q * H[i][n1] - p *
z;
409 for (
int i = low; i <= high; i++) {
411 V[i][n1 - 1] = q * z + p * V[i][n1];
412 V[i][n1] = q * V[i][n1] - p *
z;
436 y = H[n1 - 1][n1 - 1];
437 w = H[n1][n1 - 1] * H[n1 - 1][n1];
444 for (
int i = low; i <= n1; i++) {
447 s = std::abs(H[n1][n1 - 1]) + std::abs(H[n1 - 1][n1 - 2]);
462 s = x - w / ((y -
x) / 2.0 + s);
463 for (
int i = low; i <= n1; i++) {
479 p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
480 q = H[m + 1][m + 1] - z - r - s;
482 s = std::abs(p) + std::abs(q) + std::abs(r);
489 if (std::abs(H[m][m - 1]) * (std::abs(q) + std::abs(r)) < eps * (std::abs(p)
490 * (std::abs(H[m - 1][m - 1]) + std::abs(z) + std::abs(
491 H[m + 1][m + 1])))) {
497 for (
int i = m + 2; i <= n1; i++) {
506 for (
int k = m; k <= n1 - 1; k++) {
507 bool notlast = (k != n1 - 1);
511 r = (notlast ? H[k + 2][k - 1] : 0.0);
512 x = std::abs(p) + std::abs(q) + std::abs(r);
522 s = std::sqrt(p * p + q * q + r * r);
528 H[k][k - 1] = -s *
x;
530 H[k][k - 1] = -H[k][k - 1];
541 for (
int j = k; j < nn; j++) {
542 p = H[k][j] + q * H[k + 1][j];
544 p = p + r * H[k + 2][j];
545 H[k + 2][j] = H[k + 2][j] - p *
z;
547 H[k][j] = H[k][j] - p *
x;
548 H[k + 1][j] = H[k + 1][j] - p * y;
553 for (
int i = 0; i <=
std::min(n1, k + 3); i++) {
554 p = x * H[i][k] + y * H[i][k + 1];
556 p = p + z * H[i][k + 2];
557 H[i][k + 2] = H[i][k + 2] - p * r;
559 H[i][k] = H[i][k] - p;
560 H[i][k + 1] = H[i][k + 1] - p * q;
565 for (
int i = low; i <= high; i++) {
566 p = x * V[i][k] + y * V[i][k + 1];
568 p = p + z * V[i][k + 2];
569 V[i][k + 2] = V[i][k + 2] - p * r;
571 V[i][k] = V[i][k] - p;
572 V[i][k + 1] = V[i][k + 1] - p * q;
585 for (n1 = nn - 1; n1 >= 0; n1--) {
594 for (
int i = n1 - 1; i >= 0; i--) {
597 for (
int j = l; j <= n1; j++) {
598 r = r + H[i][j] * H[j][n1];
609 H[i][n1] = -r / (eps *
norm);
617 q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
618 t = (x * s - z * r) / q;
620 if (std::abs(x) > std::abs(z)) {
621 H[i + 1][n1] = (-r - w *
t) / x;
623 H[i + 1][n1] = (-s - y *
t) / z;
629 t = std::abs(H[i][n1]);
630 if ((eps *
t) * t > 1) {
631 for (
int j = i; j <= n1; j++) {
632 H[j][n1] = H[j][n1] /
t;
643 if (std::abs(H[n1][n1 - 1]) > std::abs(H[n1 - 1][n1])) {
644 H[n1 - 1][n1 - 1] = q / H[n1][n1 - 1];
645 H[n1 - 1][n1] = -(H[n1][n1] - p) / H[n1][n1 - 1];
647 cdiv(0.0, -H[n1 - 1][n1], H[n1 - 1][n1 - 1] - p, q);
648 H[n1 - 1][n1 - 1] = cdivr;
649 H[n1 - 1][n1] = cdivi;
653 for (
int i = n1 - 2; i >= 0; i--) {
657 for (
int j = l; j <= n1; j++) {
658 ra = ra + H[i][j] * H[j][n1 - 1];
659 sa = sa + H[i][j] * H[j][n1];
670 cdiv(-ra, -sa, w, q);
671 H[i][n1 - 1] = cdivr;
679 double vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
680 double vi = (d[i] - p) * 2.0 * q;
681 if (vr == 0.0 && vi == 0.0) {
682 vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x)
683 + std::abs(y) + std::abs(z));
685 cdiv(x * r - z * ra + q * sa,
686 x * s - z * sa - q * ra, vr, vi);
687 H[i][n1 - 1] = cdivr;
689 if (std::abs(x) > (std::abs(z) + std::abs(q))) {
690 H[i + 1][n1 - 1] = (-ra - w * H[i][n1 - 1] + q
692 H[i + 1][n1] = (-sa - w * H[i][n1] - q * H[i][n1
695 cdiv(-r - y * H[i][n1 - 1], -s - y * H[i][n1], z,
697 H[i + 1][n1 - 1] = cdivr;
698 H[i + 1][n1] = cdivi;
704 t =
std::max(std::abs(H[i][n1 - 1]), std::abs(H[i][n1]));
705 if ((eps *
t) * t > 1) {
706 for (
int j = i; j <= n1; j++) {
707 H[j][n1 - 1] = H[j][n1 - 1] /
t;
708 H[j][n1] = H[j][n1] /
t;
718 for (
int i = 0; i < nn; i++) {
719 if (i < low || i > high) {
720 for (
int j = i; j < nn; j++) {
728 for (
int j = nn - 1; j >= low; j--) {
729 for (
int i = low; i <= high; i++) {
731 for (
int k = low; k <=
std::min(j, high); k++) {
732 z = z + V[i][k] * H[k][j];
752 for (
int m = low + 1; m <= high - 1; m++) {
757 for (
int i = m; i <= high; i++) {
758 scale = scale + std::abs(H[i][m - 1]);
765 for (
int i = high; i >= m; i--) {
766 ort[i] = H[i][m - 1] / scale;
767 h += ort[i] * ort[i];
769 double g = std::sqrt(h);
779 for (
int j = m; j < n; j++) {
781 for (
int i = high; i >= m; i--) {
782 f += ort[i] * H[i][j];
785 for (
int i = m; i <= high; i++) {
786 H[i][j] -= f * ort[i];
790 for (
int i = 0; i <= high; i++) {
792 for (
int j = high; j >= m; j--) {
793 f += ort[j] * H[i][j];
796 for (
int j = m; j <= high; j++) {
797 H[i][j] -= f * ort[j];
800 ort[m] = scale * ort[m];
801 H[m][m - 1] = scale * g;
807 for (
int i = 0; i < n; i++) {
808 for (
int j = 0; j < n; j++) {
809 V[i][j] = (i == j ? 1.0 : 0.0);
813 for (
int m = high - 1; m >= low + 1; m--) {
814 if (H[m][m - 1] != 0.0) {
815 for (
int i = m + 1; i <= high; i++) {
816 ort[i] = H[i][m - 1];
818 for (
int j = m; j <= high; j++) {
820 for (
int i = m; i <= high; i++) {
821 g += ort[i] * V[i][j];
824 g = (g / ort[m]) / H[m][m - 1];
825 for (
int i = m; i <= high; i++) {
826 V[i][j] += g * ort[i];
841 for (
int i = 0; i < n; i++) {
854 V = alloc_2d<double> (n, n, 0.0);
855 d = alloc_1d<double> (n);
856 e = alloc_1d<double> (n);
857 ort = alloc_1d<double> (n);
863 _eigenvalues.create(1, n, CV_64FC1);
864 for (
int i = 0; i < n; i++) {
865 _eigenvalues.at<
double> (0, i) = d[i];
868 _eigenvectors.create(n, n, CV_64FC1);
869 for (
int i = 0; i < n; i++)
870 for (
int j = 0; j < n; j++)
871 _eigenvectors.at<
double> (i, j) = V[i][j];
879 EigenvalueDecomposition()
888 EigenvalueDecomposition(cv::InputArray src) {
897 void compute(cv::InputArray src)
907 src.getMat().convertTo(tmp, CV_64FC1);
911 this->H = alloc_2d<double> (n, n);
913 for (
int i = 0; i < tmp.rows; i++) {
914 for (
int j = 0; j < tmp.cols; j++) {
915 this->H[i][j] = tmp.at<
double>(i, j);
926 ~EigenvalueDecomposition() {}
932 cv::Mat eigenvalues() {
return _eigenvalues; }
938 cv::Mat eigenvectors() {
return _eigenvectors; }
950 #endif // OPENCV_Check bool compute(const std::string &expression_string, T &result)
EIGEN_STRONG_INLINE const AdjointReturnType t() const
Transpose.
T max(const T v0, const T v1)
T min(const T v0, const T v1)
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
x y t t *t x y t t t x y t t t x *y t *t t x *y t *t t x y t t t x y t t t x(y+z)
T pow(const T v0, const T v1)
EIGEN_STRONG_INLINE double mean() const
Computes the mean of the entire matrix.
CONTAINER::Scalar norm(const CONTAINER &v)