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
// 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
// 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