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>(i).x;
28
p.at<double>(1,i) = opoints.at<OpointType>(i).y;
29
p.at<double>(2,i) = opoints.at<OpointType>(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>(i).x, 2) +
38
std::pow(ipoints.at<IpointType>(i).y, 2) + (double)1;
41
z.at<double>(0,i) = ipoints.at<IpointType>(i).x / sr;
42
z.at<double>(1,i) = ipoints.at<IpointType>(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--) {
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
double vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
533
double 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; }