implemented rotating-only cameras calibration

modules/stitching/autocalib.cpp

@@ -115,3 +115,70 @@ void estimateFocal(const vector<ImageFeatures> &features, const vector<MatchesIn
             focals[i] = focals_sum / num_images;
     }
 }
+
+
+namespace
+{
+    template<typename _Tp> static inline bool
+    decomposeCholesky(_Tp* A, size_t astep, int m)
+    {
+        if (!Cholesky(A, astep, m, 0, 0, 0))
+            return false;
+        astep /= sizeof(A[0]);
+        for (int i = 0; i < m; ++i)
+            A[i*astep + i] = (_Tp)(1./A[i*astep + i]);
+        return true;
+    }
+} // namespace
+
+
+bool calibrateRotatingCamera(const vector<Mat> &Hs, Mat &K)
+{
+    int m = static_cast<int>(Hs.size());
+    CV_Assert(m >= 1);
+
+    vector<Mat> Hs_(m);
+    for (int i = 0; i < m; ++i)
+    {
+        CV_Assert(Hs[i].size() == Size(3, 3) && Hs[i].type() == CV_64F);
+        Hs_[i] = Hs[i] / pow(determinant(Hs[i]), 1./3.);
+    }
+
+    const int idx_map[3][3] = {{0, 1, 2}, {1, 3, 4}, {2, 4, 5}};
+    Mat_<double> A(6*m, 6);
+    A.setTo(0);
+
+    int eq_idx = 0;
+    for (int k = 0; k < m; ++k)
+    {
+        Mat_<double> H(Hs_[k]);
+        for (int i = 0; i < 3; ++i)
+        {
+            for (int j = i; j < 3; ++j, ++eq_idx)
+            {
+                for (int l = 0; l < 3; ++l)
+                {
+                    for (int s = 0; s < 3; ++s)
+                    {
+                        int idx = idx_map[l][s];
+                        A(eq_idx, idx) += H(i,l) * H(j,s);
+                    }
+                }
+                A(eq_idx, idx_map[i][j]) -= 1;
+            }
+        }
+    }
+
+    Mat_<double> wcoef;
+    SVD::solveZ(A, wcoef);
+
+    Mat_<double> W(3,3);
+    for (int i = 0; i < 3; ++i)
+        for (int j = i; j < 3; ++j)
+            W(i,j) = W(j,i) = wcoef(idx_map[i][j], 0) / wcoef(5,0);
+    if (!decomposeCholesky(W.ptr<double>(), W.step, 3))
+        return false;
+    W(0,1) = W(0,2) = W(1,2) = 0;
+    K = W.t();
+    return true;
+}

modules/stitching/autocalib.hpp

@@ -38,18 +38,20 @@
 // or tort (including negligence or otherwise) arising in any way out of
 // the use of this software, even if advised of the possibility of such damage.
 //
-//M*/
-#ifndef __OPENCV_AUTOCALIB_HPP__
-#define __OPENCV_AUTOCALIB_HPP__
-
-#include "precomp.hpp"
-#include "matchers.hpp"
-
-// See "Construction of Panoramic Image Mosaics with Global and Local Alignment"
-// by Heung-Yeung Shum and Richard Szeliski.
-void focalsFromHomography(const cv::Mat &H, double &f0, double &f1, bool &f0_ok, bool &f1_ok);
-
-void estimateFocal(const std::vector<ImageFeatures> &features, const std::vector<MatchesInfo> &pairwise_matches, 
-                   std::vector<double> &focals);
-
-#endif // __OPENCV_AUTOCALIB_HPP__
+//M*/
+#ifndef __OPENCV_AUTOCALIB_HPP__
+#define __OPENCV_AUTOCALIB_HPP__
+
+#include "precomp.hpp"
+#include "matchers.hpp"
+
+// See "Construction of Panoramic Image Mosaics with Global and Local Alignment"
+// by Heung-Yeung Shum and Richard Szeliski.
+void focalsFromHomography(const cv::Mat &H, double &f0, double &f1, bool &f0_ok, bool &f1_ok);
+
+void estimateFocal(const std::vector<ImageFeatures> &features, const std::vector<MatchesInfo> &pairwise_matches, 
+                   std::vector<double> &focals);
+
+bool calibrateRotatingCamera(const std::vector<cv::Mat> &Hs, cv::Mat &K);
+
+#endif // __OPENCV_AUTOCALIB_HPP__

modules/stitching/main.cpp

@@ -100,7 +100,7 @@ void printUsage()
         "  --blend_strength <float>\n"
         "      Blending strength from [0,100] range. The default is 5.\n"
         "  --output <result_img>\n"
-        "      The default is 'result.png'.\n";
+        "      The default is 'result.jpg'.\n";
 }
 
 
@@ -120,7 +120,7 @@ float match_conf = 0.65f;
 int seam_find_type = SeamFinder::GC_COLOR;
 int blend_type = Blender::MULTI_BAND;
 float blend_strength = 5;
-string result_name = "result.png";
+string result_name = "result.jpg";
 
 int parseCmdArgs(int argc, char** argv)
 {

modules/stitching/motion_estimators.cpp

@@ -108,6 +108,27 @@ void HomographyBasedEstimator::estimate(const vector<ImageFeatures> &features, c
 {
     const int num_images = static_cast<int>(features.size());
 
+#if 0
+    // Robustly estimate focal length from rotating cameras
+    vector<Mat> Hs;
+    for (int iter = 0; iter < 100; ++iter)
+    {
+        int len = 2 + rand()%(pairwise_matches.size() - 1);
+        vector<int> subset;
+        selectRandomSubset(len, pairwise_matches.size(), subset);
+        Hs.clear();
+        for (size_t i = 0; i < subset.size(); ++i)
+            if (!pairwise_matches[subset[i]].H.empty())
+                Hs.push_back(pairwise_matches[subset[i]].H);
+        Mat_<double> K;
+        if (Hs.size() >= 2)
+        {
+            if (calibrateRotatingCamera(Hs, K))
+                cin.get();
+        }
+    }
+#endif
+
     // Estimate focal length and set it for all cameras
     vector<double> focals;
     estimateFocal(features, pairwise_matches, focals);

modules/stitching/util.cpp

@@ -38,113 +38,126 @@
 // or tort (including negligence or otherwise) arising in any way out of
 // the use of this software, even if advised of the possibility of such damage.
 //
-//M*/
-#include "util.hpp"
-
-using namespace std;
-using namespace cv;
-
-void DjSets::create(int n) 
-{
-    rank_.assign(n, 0);
-    size.assign(n, 1);
-    parent.resize(n);
-    for (int i = 0; i < n; ++i)
-        parent[i] = i;
-}
-
-
-int DjSets::find(int elem) 
-{
-    int set = elem;
-    while (set != parent[set])
-        set = parent[set];
-    int next;
-    while (elem != parent[elem]) 
-    {
-        next = parent[elem];
-        parent[elem] = set;
-        elem = next;
-    }
-    return set;
-}
-
-
-int DjSets::merge(int set1, int set2) 
-{
-    if (rank_[set1] < rank_[set2]) 
-    {
-        parent[set1] = set2;
-        size[set2] += size[set1];
-        return set2;
-    }
-    if (rank_[set2] < rank_[set1]) 
-    {
-        parent[set2] = set1;
-        size[set1] += size[set2];
-        return set1;
-    }
-    parent[set1] = set2;
-    rank_[set2]++;
-    size[set2] += size[set1];
-    return set2;
-}
-
-
-void Graph::addEdge(int from, int to, float weight)
-{
-    edges_[from].push_back(GraphEdge(from, to, weight));
-}
-
-
-bool overlapRoi(Point tl1, Point tl2, Size sz1, Size sz2, Rect &roi)
-{
-    int x_tl = max(tl1.x, tl2.x);
-    int y_tl = max(tl1.y, tl2.y);
-    int x_br = min(tl1.x + sz1.width, tl2.x + sz2.width);
-    int y_br = min(tl1.y + sz1.height, tl2.y + sz2.height);
-    if (x_tl < x_br && y_tl < y_br)
-    {
-        roi = Rect(x_tl, y_tl, x_br - x_tl, y_br - y_tl);
-        return true;
-    }
-    return false;
-}
-
-
-Rect resultRoi(const vector<Point> &corners, const vector<Mat> &images)
-{
-    vector<Size> sizes(images.size());
-    for (size_t i = 0; i < images.size(); ++i)
-        sizes[i] = images[i].size();
-    return resultRoi(corners, sizes);
-}
-
-
-Rect resultRoi(const vector<Point> &corners, const vector<Size> &sizes)
-{
-    CV_Assert(sizes.size() == corners.size());
-    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
-    Point br(numeric_limits<int>::min(), numeric_limits<int>::min());
-    for (size_t i = 0; i < corners.size(); ++i)
-    {
-        tl.x = min(tl.x, corners[i].x);
-        tl.y = min(tl.y, corners[i].y);
-        br.x = max(br.x, corners[i].x + sizes[i].width);
-        br.y = max(br.y, corners[i].y + sizes[i].height);
-    }
-    return Rect(tl, br);
-}
-
-
-Point resultTl(const vector<Point> &corners)
-{
-    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
-    for (size_t i = 0; i < corners.size(); ++i)
-    {
-        tl.x = min(tl.x, corners[i].x);
-        tl.y = min(tl.y, corners[i].y);
-    }
-    return tl;
-}
-
+//M*/
+#include "util.hpp"
+
+using namespace std;
+using namespace cv;
+
+void DjSets::create(int n) 
+{
+    rank_.assign(n, 0);
+    size.assign(n, 1);
+    parent.resize(n);
+    for (int i = 0; i < n; ++i)
+        parent[i] = i;
+}
+
+
+int DjSets::find(int elem) 
+{
+    int set = elem;
+    while (set != parent[set])
+        set = parent[set];
+    int next;
+    while (elem != parent[elem]) 
+    {
+        next = parent[elem];
+        parent[elem] = set;
+        elem = next;
+    }
+    return set;
+}
+
+
+int DjSets::merge(int set1, int set2) 
+{
+    if (rank_[set1] < rank_[set2]) 
+    {
+        parent[set1] = set2;
+        size[set2] += size[set1];
+        return set2;
+    }
+    if (rank_[set2] < rank_[set1]) 
+    {
+        parent[set2] = set1;
+        size[set1] += size[set2];
+        return set1;
+    }
+    parent[set1] = set2;
+    rank_[set2]++;
+    size[set2] += size[set1];
+    return set2;
+}
+
+
+void Graph::addEdge(int from, int to, float weight)
+{
+    edges_[from].push_back(GraphEdge(from, to, weight));
+}
+
+
+bool overlapRoi(Point tl1, Point tl2, Size sz1, Size sz2, Rect &roi)
+{
+    int x_tl = max(tl1.x, tl2.x);
+    int y_tl = max(tl1.y, tl2.y);
+    int x_br = min(tl1.x + sz1.width, tl2.x + sz2.width);
+    int y_br = min(tl1.y + sz1.height, tl2.y + sz2.height);
+    if (x_tl < x_br && y_tl < y_br)
+    {
+        roi = Rect(x_tl, y_tl, x_br - x_tl, y_br - y_tl);
+        return true;
+    }
+    return false;
+}
+
+
+Rect resultRoi(const vector<Point> &corners, const vector<Mat> &images)
+{
+    vector<Size> sizes(images.size());
+    for (size_t i = 0; i < images.size(); ++i)
+        sizes[i] = images[i].size();
+    return resultRoi(corners, sizes);
+}
+
+
+Rect resultRoi(const vector<Point> &corners, const vector<Size> &sizes)
+{
+    CV_Assert(sizes.size() == corners.size());
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
+    Point br(numeric_limits<int>::min(), numeric_limits<int>::min());
+    for (size_t i = 0; i < corners.size(); ++i)
+    {
+        tl.x = min(tl.x, corners[i].x);
+        tl.y = min(tl.y, corners[i].y);
+        br.x = max(br.x, corners[i].x + sizes[i].width);
+        br.y = max(br.y, corners[i].y + sizes[i].height);
+    }
+    return Rect(tl, br);
+}
+
+
+Point resultTl(const vector<Point> &corners)
+{
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
+    for (size_t i = 0; i < corners.size(); ++i)
+    {
+        tl.x = min(tl.x, corners[i].x);
+        tl.y = min(tl.y, corners[i].y);
+    }
+    return tl;
+}
+
+
+void selectRandomSubset(int count, int size, vector<int> &subset)
+{
+    subset.clear();
+    for (int i = 0; i < size; ++i)
+    {
+        if (randu<int>() % (size - i) < count)
+        {
+            subset.push_back(i);
+            count--;
+        }
+    }
+}

modules/stitching/util.hpp

@@ -108,6 +108,7 @@ bool overlapRoi(cv::Point tl1, cv::Point tl2, cv::Size sz1, cv::Size sz2, cv::Re
 cv::Rect resultRoi(const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &images);
 cv::Rect resultRoi(const std::vector<cv::Point> &corners, const std::vector<cv::Size> &sizes);
 cv::Point resultTl(const std::vector<cv::Point> &corners);
+void selectRandomSubset(int count, int size, std::vector<int> &subset);
 
 
 #include "util_inl.hpp"