implemented rotating-only cameras calibration
authorAlexey Spizhevoy <no@email>
Mon, 15 Aug 2011 06:15:06 +0000 (06:15 +0000)
committerAlexey Spizhevoy <no@email>
Mon, 15 Aug 2011 06:15:06 +0000 (06:15 +0000)
modules/stitching/autocalib.cpp
modules/stitching/autocalib.hpp
modules/stitching/main.cpp
modules/stitching/motion_estimators.cpp
modules/stitching/util.cpp
modules/stitching/util.hpp

index fbd2967..03abb63 100644 (file)
@@ -115,3 +115,70 @@ void estimateFocal(const vector<ImageFeatures> &features, const vector<MatchesIn
             focals[i] = focals_sum / num_images;\r
     }\r
 }\r
+\r
+\r
+namespace\r
+{\r
+    template<typename _Tp> static inline bool\r
+    decomposeCholesky(_Tp* A, size_t astep, int m)\r
+    {\r
+        if (!Cholesky(A, astep, m, 0, 0, 0))\r
+            return false;\r
+        astep /= sizeof(A[0]);\r
+        for (int i = 0; i < m; ++i)\r
+            A[i*astep + i] = (_Tp)(1./A[i*astep + i]);\r
+        return true;\r
+    }\r
+} // namespace\r
+\r
+\r
+bool calibrateRotatingCamera(const vector<Mat> &Hs, Mat &K)\r
+{\r
+    int m = static_cast<int>(Hs.size());\r
+    CV_Assert(m >= 1);\r
+\r
+    vector<Mat> Hs_(m);\r
+    for (int i = 0; i < m; ++i)\r
+    {\r
+        CV_Assert(Hs[i].size() == Size(3, 3) && Hs[i].type() == CV_64F);\r
+        Hs_[i] = Hs[i] / pow(determinant(Hs[i]), 1./3.);\r
+    }\r
+\r
+    const int idx_map[3][3] = {{0, 1, 2}, {1, 3, 4}, {2, 4, 5}};\r
+    Mat_<double> A(6*m, 6);\r
+    A.setTo(0);\r
+\r
+    int eq_idx = 0;\r
+    for (int k = 0; k < m; ++k)\r
+    {\r
+        Mat_<double> H(Hs_[k]);\r
+        for (int i = 0; i < 3; ++i)\r
+        {\r
+            for (int j = i; j < 3; ++j, ++eq_idx)\r
+            {\r
+                for (int l = 0; l < 3; ++l)\r
+                {\r
+                    for (int s = 0; s < 3; ++s)\r
+                    {\r
+                        int idx = idx_map[l][s];\r
+                        A(eq_idx, idx) += H(i,l) * H(j,s);\r
+                    }\r
+                }\r
+                A(eq_idx, idx_map[i][j]) -= 1;\r
+            }\r
+        }\r
+    }\r
+\r
+    Mat_<double> wcoef;\r
+    SVD::solveZ(A, wcoef);\r
+\r
+    Mat_<double> W(3,3);\r
+    for (int i = 0; i < 3; ++i)\r
+        for (int j = i; j < 3; ++j)\r
+            W(i,j) = W(j,i) = wcoef(idx_map[i][j], 0) / wcoef(5,0);\r
+    if (!decomposeCholesky(W.ptr<double>(), W.step, 3))\r
+        return false;\r
+    W(0,1) = W(0,2) = W(1,2) = 0;\r
+    K = W.t();\r
+    return true;\r
+}\r
index 53ecc4d..0807ee7 100644 (file)
 // or tort (including negligence or otherwise) arising in any way out of\r
 // the use of this software, even if advised of the possibility of such damage.\r
 //\r
-//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*/\r
+#ifndef __OPENCV_AUTOCALIB_HPP__\r
+#define __OPENCV_AUTOCALIB_HPP__\r
+\r
+#include "precomp.hpp"\r
+#include "matchers.hpp"\r
+\r
+// See "Construction of Panoramic Image Mosaics with Global and Local Alignment"\r
+// by Heung-Yeung Shum and Richard Szeliski.\r
+void focalsFromHomography(const cv::Mat &H, double &f0, double &f1, bool &f0_ok, bool &f1_ok);\r
+\r
+void estimateFocal(const std::vector<ImageFeatures> &features, const std::vector<MatchesInfo> &pairwise_matches, \r
+                   std::vector<double> &focals);\r
+\r
+bool calibrateRotatingCamera(const std::vector<cv::Mat> &Hs, cv::Mat &K);\r
+\r
+#endif // __OPENCV_AUTOCALIB_HPP__\r
index fbd4657..76a4c44 100644 (file)
@@ -100,7 +100,7 @@ void printUsage()
         "  --blend_strength <float>\n"\r
         "      Blending strength from [0,100] range. The default is 5.\n"\r
         "  --output <result_img>\n"\r
-        "      The default is 'result.png'.\n";\r
+        "      The default is 'result.jpg'.\n";\r
 }\r
 \r
 \r
@@ -120,7 +120,7 @@ float match_conf = 0.65f;
 int seam_find_type = SeamFinder::GC_COLOR;\r
 int blend_type = Blender::MULTI_BAND;\r
 float blend_strength = 5;\r
-string result_name = "result.png";\r
+string result_name = "result.jpg";\r
 \r
 int parseCmdArgs(int argc, char** argv)\r
 {\r
index 08d1aae..d354456 100644 (file)
@@ -108,6 +108,27 @@ void HomographyBasedEstimator::estimate(const vector<ImageFeatures> &features, c
 {\r
     const int num_images = static_cast<int>(features.size());\r
 \r
+#if 0\r
+    // Robustly estimate focal length from rotating cameras\r
+    vector<Mat> Hs;\r
+    for (int iter = 0; iter < 100; ++iter)\r
+    {\r
+        int len = 2 + rand()%(pairwise_matches.size() - 1);\r
+        vector<int> subset;\r
+        selectRandomSubset(len, pairwise_matches.size(), subset);\r
+        Hs.clear();\r
+        for (size_t i = 0; i < subset.size(); ++i)\r
+            if (!pairwise_matches[subset[i]].H.empty())\r
+                Hs.push_back(pairwise_matches[subset[i]].H);\r
+        Mat_<double> K;\r
+        if (Hs.size() >= 2)\r
+        {\r
+            if (calibrateRotatingCamera(Hs, K))\r
+                cin.get();\r
+        }\r
+    }\r
+#endif\r
+\r
     // Estimate focal length and set it for all cameras\r
     vector<double> focals;\r
     estimateFocal(features, pairwise_matches, focals);\r
index 1bcb7f8..6cba9e3 100644 (file)
 // or tort (including negligence or otherwise) arising in any way out of\r
 // the use of this software, even if advised of the possibility of such damage.\r
 //\r
-//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*/\r
+#include "util.hpp"\r
+\r
+using namespace std;\r
+using namespace cv;\r
+\r
+void DjSets::create(int n) \r
+{\r
+    rank_.assign(n, 0);\r
+    size.assign(n, 1);\r
+    parent.resize(n);\r
+    for (int i = 0; i < n; ++i)\r
+        parent[i] = i;\r
+}\r
+\r
+\r
+int DjSets::find(int elem) \r
+{\r
+    int set = elem;\r
+    while (set != parent[set])\r
+        set = parent[set];\r
+    int next;\r
+    while (elem != parent[elem]) \r
+    {\r
+        next = parent[elem];\r
+        parent[elem] = set;\r
+        elem = next;\r
+    }\r
+    return set;\r
+}\r
+\r
+\r
+int DjSets::merge(int set1, int set2) \r
+{\r
+    if (rank_[set1] < rank_[set2]) \r
+    {\r
+        parent[set1] = set2;\r
+        size[set2] += size[set1];\r
+        return set2;\r
+    }\r
+    if (rank_[set2] < rank_[set1]) \r
+    {\r
+        parent[set2] = set1;\r
+        size[set1] += size[set2];\r
+        return set1;\r
+    }\r
+    parent[set1] = set2;\r
+    rank_[set2]++;\r
+    size[set2] += size[set1];\r
+    return set2;\r
+}\r
+\r
+\r
+void Graph::addEdge(int from, int to, float weight)\r
+{\r
+    edges_[from].push_back(GraphEdge(from, to, weight));\r
+}\r
+\r
+\r
+bool overlapRoi(Point tl1, Point tl2, Size sz1, Size sz2, Rect &roi)\r
+{\r
+    int x_tl = max(tl1.x, tl2.x);\r
+    int y_tl = max(tl1.y, tl2.y);\r
+    int x_br = min(tl1.x + sz1.width, tl2.x + sz2.width);\r
+    int y_br = min(tl1.y + sz1.height, tl2.y + sz2.height);\r
+    if (x_tl < x_br && y_tl < y_br)\r
+    {\r
+        roi = Rect(x_tl, y_tl, x_br - x_tl, y_br - y_tl);\r
+        return true;\r
+    }\r
+    return false;\r
+}\r
+\r
+\r
+Rect resultRoi(const vector<Point> &corners, const vector<Mat> &images)\r
+{\r
+    vector<Size> sizes(images.size());\r
+    for (size_t i = 0; i < images.size(); ++i)\r
+        sizes[i] = images[i].size();\r
+    return resultRoi(corners, sizes);\r
+}\r
+\r
+\r
+Rect resultRoi(const vector<Point> &corners, const vector<Size> &sizes)\r
+{\r
+    CV_Assert(sizes.size() == corners.size());\r
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());\r
+    Point br(numeric_limits<int>::min(), numeric_limits<int>::min());\r
+    for (size_t i = 0; i < corners.size(); ++i)\r
+    {\r
+        tl.x = min(tl.x, corners[i].x);\r
+        tl.y = min(tl.y, corners[i].y);\r
+        br.x = max(br.x, corners[i].x + sizes[i].width);\r
+        br.y = max(br.y, corners[i].y + sizes[i].height);\r
+    }\r
+    return Rect(tl, br);\r
+}\r
+\r
+\r
+Point resultTl(const vector<Point> &corners)\r
+{\r
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());\r
+    for (size_t i = 0; i < corners.size(); ++i)\r
+    {\r
+        tl.x = min(tl.x, corners[i].x);\r
+        tl.y = min(tl.y, corners[i].y);\r
+    }\r
+    return tl;\r
+}\r
+\r
+\r
+void selectRandomSubset(int count, int size, vector<int> &subset)\r
+{\r
+    subset.clear();\r
+    for (int i = 0; i < size; ++i)\r
+    {\r
+        if (randu<int>() % (size - i) < count)\r
+        {\r
+            subset.push_back(i);\r
+            count--;\r
+        }\r
+    }\r
+}\r
index aa3cfb2..6bec221 100644 (file)
@@ -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);\r
 cv::Rect resultRoi(const std::vector<cv::Point> &corners, const std::vector<cv::Size> &sizes);\r
 cv::Point resultTl(const std::vector<cv::Point> &corners);\r
+void selectRandomSubset(int count, int size, std::vector<int> &subset);\r
 \r
 \r
 #include "util_inl.hpp"\r