Skip to content

Commit 60e1eda

Browse files
author
Alexey Spizhevoy
committed
modified focal estimation function in opencv_stitching
1 parent 34e2c78 commit 60e1eda

3 files changed

Lines changed: 55 additions & 118 deletions

File tree

‎modules/stitching/autocalib.cpp‎

Lines changed: 47 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -8,99 +8,63 @@ void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok, boo
88
{
99
CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3));
1010

11-
const double h[9] =
12-
{
13-
H.at<double>(0, 0), H.at<double>(0, 1), H.at<double>(0, 2),
14-
H.at<double>(1, 0), H.at<double>(1, 1), H.at<double>(1, 2),
15-
H.at<double>(2, 0), H.at<double>(2, 1), H.at<double>(2, 2)
16-
};
11+
const double* h = reinterpret_cast<const double*>(H.data);
12+
13+
double d1, d2; // Denominators
14+
double v1, v2; // Focal squares value candidates
1715

1816
f1_ok = true;
19-
double denom1 = h[6] * h[7];
20-
double denom2 = (h[7] - h[6]) * (h[7] + h[6]);
21-
if (max(abs(denom1), abs(denom2)) < 1e-5)
22-
f1_ok = false;
23-
else
24-
{
25-
double val1 = -(h[0] * h[1] + h[3] * h[4]) / denom1;
26-
double val2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / denom2;
27-
if (val1 < val2)
28-
swap(val1, val2);
29-
if (val1 > 0 && val2 > 0)
30-
f1 = sqrt(abs(denom1) > abs(denom2) ? val1 : val2);
31-
else if (val1 > 0)
32-
f1 = sqrt(val1);
33-
else
34-
f1_ok = false;
35-
}
17+
d1 = h[6] * h[7];
18+
d2 = (h[7] - h[6]) * (h[7] + h[6]);
19+
v1 = -(h[0] * h[1] + h[3] * h[4]) / d1;
20+
v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2;
21+
if (v1 < v2) swap(v1, v2);
22+
if (v1 > 0 && v2 > 0) f1 = sqrt(abs(d1) > abs(d2) ? v1 : v2);
23+
else if (v1 > 0) f1 = sqrt(v1);
24+
else f1_ok = false;
3625

3726
f0_ok = true;
38-
denom1 = h[0] * h[3] + h[1] * h[4];
39-
denom2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];
40-
if (max(abs(denom1), abs(denom2)) < 1e-5)
41-
f0_ok = false;
42-
else
43-
{
44-
double val1 = -h[2] * h[5] / denom1;
45-
double val2 = (h[5] * h[5] - h[2] * h[2]) / denom2;
46-
if (val1 < val2)
47-
swap(val1, val2);
48-
if (val1 > 0 && val2 > 0)
49-
f0 = sqrt(abs(denom1) > abs(denom2) ? val1 : val2);
50-
else if (val1 > 0)
51-
f0 = sqrt(val1);
52-
else
53-
f0_ok = false;
54-
}
27+
d1 = h[0] * h[3] + h[1] * h[4];
28+
d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];
29+
v1 = -h[2] * h[5] / d1;
30+
v2 = (h[5] * h[5] - h[2] * h[2]) / d2;
31+
if (v1 < v2) swap(v1, v2);
32+
if (v1 > 0 && v2 > 0) f0 = sqrt(abs(d1) > abs(d2) ? v1 : v2);
33+
else if (v1 > 0) f0 = sqrt(v1);
34+
else f0_ok = false;
5535
}
5636

5737

58-
bool focalsFromFundamental(const Mat &F, double &f0, double &f1)
38+
double estimateFocal(const vector<Mat> &images, const vector<ImageFeatures> &/*features*/,
39+
const vector<MatchesInfo> &pairwise_matches)
5940
{
60-
CV_Assert(F.type() == CV_64F);
61-
CV_Assert(F.size() == Size(3, 3));
62-
63-
Mat Ft = F.t();
64-
Mat k = Mat::zeros(3, 1, CV_64F);
65-
k.at<double>(2, 0) = 1.f;
41+
const int num_images = static_cast<int>(images.size());
6642

67-
// 1. Compute quantities
68-
double a = normL2sq(F*Ft*k) / normL2sq(Ft*k);
69-
double b = normL2sq(Ft*F*k) / normL2sq(F*k);
70-
double c = sqr(k.dot(F*k)) / (normL2sq(Ft*k) * normL2sq(F*k));
71-
double d = k.dot(F*Ft*F*k) / k.dot(F*k);
72-
double A = 1/c + a - 2*d;
73-
double B = 1/c + b - 2*d;
74-
double P = 2*(1/c - 2*d + 0.5*normL2sq(F));
75-
double Q = -(A + B)/c + 0.5*(normL2sq(F*Ft) - 0.5*sqr(normL2sq(F)));
76-
77-
// 2. Solve quadratic equation Z*Z*a_ + Z*b_ + c_ = 0
78-
double a_ = 1 + c*P;
79-
double b_ = -(c*P*P + 2*P + 4*c*Q);
80-
double c_ = P*P + 4*c*P*Q + 12*A*B;
81-
double D = b_*b_ - 4*a_*c_;
82-
if (abs(D) < 1e-5)
83-
D = 0;
84-
else if (D < 0)
85-
return false;
86-
double D_sqrt = sqrt(D);
87-
double Z0 = (-b_ - D_sqrt) / (2*a_);
88-
double Z1 = (-b_ + D_sqrt) / (2*a_);
89-
90-
// 3. Choose solution
91-
double w0 = abs(Z0*Z0*Z0 - 3*P*Z0*Z0 + 2*(P*P + 2*Q)*Z0 - 4*(P*Q + 4*A*B/c));
92-
double w1 = abs(Z1*Z1*Z1 - 3*P*Z1*Z1 + 2*(P*P + 2*Q)*Z1 - 4*(P*Q + 4*A*B/c));
93-
double Z = Z0;
94-
if (w1 < w0)
95-
Z = Z1;
43+
vector<double> focals;
44+
for (int src_idx = 0; src_idx < num_images; ++src_idx)
45+
{
46+
for (int dst_idx = 0; dst_idx < num_images; ++dst_idx)
47+
{
48+
const MatchesInfo &m = pairwise_matches[src_idx*num_images + dst_idx];
49+
if (m.H.empty())
50+
continue;
9651

97-
// 4.
98-
double X = -1/c*(1 + 2*B/(Z - P));
99-
double Y = -1/c*(1 + 2*A/(Z - P));
52+
double f0, f1;
53+
bool f0ok, f1ok;
54+
focalsFromHomography(m.H, f0, f1, f0ok, f1ok);
55+
if (f0ok && f1ok) focals.push_back(sqrt(f0*f1));
56+
}
57+
}
10058

101-
// 5. Compute focal lengths
102-
f0 = 1/sqrt(1 + X/normL2sq(Ft*k));
103-
f1 = 1/sqrt(1 + Y/normL2sq(F*k));
59+
if (focals.size() + 1 >= images.size())
60+
{
61+
nth_element(focals.begin(), focals.end(), focals.begin() + focals.size()/2);
62+
return focals[focals.size()/2];
63+
}
10464

105-
return true;
65+
LOGLN("Can't estimate focal length, will use naive approach");
66+
double focals_sum = 0;
67+
for (int i = 0; i < num_images; ++i)
68+
focals_sum += images[i].rows + images[i].cols;
69+
return focals_sum / num_images;
10670
}

‎modules/stitching/autocalib.hpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,15 @@
11
#ifndef __OPENCV_AUTOCALIB_HPP__
22
#define __OPENCV_AUTOCALIB_HPP__
33

4+
#include <vector>
45
#include <opencv2/core/core.hpp>
6+
#include "matchers.hpp"
57

68
// See "Construction of Panoramic Image Mosaics with Global and Local Alignment"
79
// by Heung-Yeung Shum and Richard Szeliski.
810
void focalsFromHomography(const cv::Mat &H, double &f0, double &f1, bool &f0_ok, bool &f1_ok);
911

10-
bool focalsFromFundamental(const cv::Mat &F, double &f0, double &f1);
12+
double estimateFocal(const std::vector<cv::Mat> &images, const std::vector<ImageFeatures> &features,
13+
const std::vector<MatchesInfo> &pairwise_matches);
1114

1215
#endif // __OPENCV_AUTOCALIB_HPP__

‎modules/stitching/motion_estimators.cpp‎

Lines changed: 4 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -64,46 +64,16 @@ struct CalcRotation
6464
};
6565

6666

67-
void HomographyBasedEstimator::estimate(const vector<Mat> &images, const vector<ImageFeatures> &/*features*/,
67+
void HomographyBasedEstimator::estimate(const vector<Mat> &images, const vector<ImageFeatures> &features,
6868
const vector<MatchesInfo> &pairwise_matches, vector<CameraParams> &cameras)
6969
{
7070
const int num_images = static_cast<int>(images.size());
7171

72-
// Find focals from pair-wise homographies
73-
vector<bool> is_focal_estimated(num_images, false);
74-
vector<double> focals;
75-
for (int i = 0; i < num_images; ++i)
76-
{
77-
for (int j = 0; j < num_images; ++j)
78-
{
79-
int pair_idx = i * num_images + j;
80-
if (pairwise_matches[pair_idx].H.empty())
81-
continue;
82-
83-
double f_to, f_from;
84-
bool f_to_ok, f_from_ok;
85-
focalsFromHomography(pairwise_matches[pair_idx].H.inv(), f_to, f_from, f_to_ok, f_from_ok);
86-
87-
if (f_from_ok) focals.push_back(f_from);
88-
if (f_to_ok) focals.push_back(f_to);
89-
90-
if (f_from_ok && f_to_ok)
91-
{
92-
is_focal_estimated[i] = true;
93-
is_focal_estimated[j] = true;
94-
}
95-
}
96-
}
97-
98-
is_focals_estimated_ = true;
99-
for (int i = 0; i < num_images; ++i)
100-
is_focals_estimated_ = is_focals_estimated_ && is_focal_estimated[i];
101-
102-
// Find focal median and use it as true focal length
103-
nth_element(focals.begin(), focals.end(), focals.begin() + focals.size() / 2);
72+
// Estimate focal length and set it for all cameras
73+
double focal = estimateFocal(images, features, pairwise_matches);
10474
cameras.resize(num_images);
10575
for (int i = 0; i < num_images; ++i)
106-
cameras[i].focal = focals[focals.size() / 2];
76+
cameras[i].focal = focal;
10777

10878
// Restore global motion
10979
Graph span_tree;

0 commit comments

Comments
 (0)