1/*M/////////////////////////////////////////////////////////////////////////////////////// 2// 3// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4// 5// By downloading, copying, installing or using the software you agree to this license. 6// If you do not agree to this license, do not download, install, 7// copy or use the software. 8// 9// 10// License Agreement 11// For Open Source Computer Vision Library 12// 13// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. 14// Copyright (C) 2009, Willow Garage Inc., all rights reserved. 15// Third party copyrights are property of their respective owners. 16// 17// Redistribution and use in source and binary forms, with or without modification, 18// are permitted provided that the following conditions are met: 19// 20// * Redistribution's of source code must retain the above copyright notice, 21// this list of conditions and the following disclaimer. 22// 23// * Redistribution's in binary form must reproduce the above copyright notice, 24// this list of conditions and the following disclaimer in the documentation 25// and/or other materials provided with the distribution. 26// 27// * The name of the copyright holders may not be used to endorse or promote products 28// derived from this software without specific prior written permission. 29// 30// This software is provided by the copyright holders and contributors "as is" and 31// any express or implied warranties, including, but not limited to, the implied 32// warranties of merchantability and fitness for a particular purpose are disclaimed. 33// In no event shall the Intel Corporation or contributors be liable for any direct, 34// indirect, incidental, special, exemplary, or consequential damages 35// (including, but not limited to, procurement of substitute goods or services; 36// loss of use, data, or profits; or business interruption) however caused 37// and on any theory of liability, whether in contract, strict liability, 38// or tort (including negligence or otherwise) arising in any way out of 39// the use of this software, even if advised of the possibility of such damage. 40// 41//M*/ 42 43#include "precomp.hpp" 44 45using namespace cv; 46 47namespace { 48 49template<typename _Tp> static inline bool 50decomposeCholesky(_Tp* A, size_t astep, int m) 51{ 52 if (!hal::Cholesky(A, astep, m, 0, 0, 0)) 53 return false; 54 astep /= sizeof(A[0]); 55 for (int i = 0; i < m; ++i) 56 A[i*astep + i] = (_Tp)(1./A[i*astep + i]); 57 return true; 58} 59 60} // namespace 61 62 63namespace cv { 64namespace detail { 65 66void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok, bool &f1_ok) 67{ 68 CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3)); 69 70 const double* h = H.ptr<double>(); 71 72 double d1, d2; // Denominators 73 double v1, v2; // Focal squares value candidates 74 75 f1_ok = true; 76 d1 = h[6] * h[7]; 77 d2 = (h[7] - h[6]) * (h[7] + h[6]); 78 v1 = -(h[0] * h[1] + h[3] * h[4]) / d1; 79 v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2; 80 if (v1 < v2) std::swap(v1, v2); 81 if (v1 > 0 && v2 > 0) f1 = std::sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2); 82 else if (v1 > 0) f1 = std::sqrt(v1); 83 else f1_ok = false; 84 85 f0_ok = true; 86 d1 = h[0] * h[3] + h[1] * h[4]; 87 d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4]; 88 v1 = -h[2] * h[5] / d1; 89 v2 = (h[5] * h[5] - h[2] * h[2]) / d2; 90 if (v1 < v2) std::swap(v1, v2); 91 if (v1 > 0 && v2 > 0) f0 = std::sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2); 92 else if (v1 > 0) f0 = std::sqrt(v1); 93 else f0_ok = false; 94} 95 96 97void estimateFocal(const std::vector<ImageFeatures> &features, const std::vector<MatchesInfo> &pairwise_matches, 98 std::vector<double> &focals) 99{ 100 const int num_images = static_cast<int>(features.size()); 101 focals.resize(num_images); 102 103 std::vector<double> all_focals; 104 105 for (int i = 0; i < num_images; ++i) 106 { 107 for (int j = 0; j < num_images; ++j) 108 { 109 const MatchesInfo &m = pairwise_matches[i*num_images + j]; 110 if (m.H.empty()) 111 continue; 112 double f0, f1; 113 bool f0ok, f1ok; 114 focalsFromHomography(m.H, f0, f1, f0ok, f1ok); 115 if (f0ok && f1ok) 116 all_focals.push_back(std::sqrt(f0 * f1)); 117 } 118 } 119 120 if (static_cast<int>(all_focals.size()) >= num_images - 1) 121 { 122 double median; 123 124 std::sort(all_focals.begin(), all_focals.end()); 125 if (all_focals.size() % 2 == 1) 126 median = all_focals[all_focals.size() / 2]; 127 else 128 median = (all_focals[all_focals.size() / 2 - 1] + all_focals[all_focals.size() / 2]) * 0.5; 129 130 for (int i = 0; i < num_images; ++i) 131 focals[i] = median; 132 } 133 else 134 { 135 LOGLN("Can't estimate focal length, will use naive approach"); 136 double focals_sum = 0; 137 for (int i = 0; i < num_images; ++i) 138 focals_sum += features[i].img_size.width + features[i].img_size.height; 139 for (int i = 0; i < num_images; ++i) 140 focals[i] = focals_sum / num_images; 141 } 142} 143 144 145bool calibrateRotatingCamera(const std::vector<Mat> &Hs, Mat &K) 146{ 147 int m = static_cast<int>(Hs.size()); 148 CV_Assert(m >= 1); 149 150 std::vector<Mat> Hs_(m); 151 for (int i = 0; i < m; ++i) 152 { 153 CV_Assert(Hs[i].size() == Size(3, 3) && Hs[i].type() == CV_64F); 154 Hs_[i] = Hs[i] / std::pow(determinant(Hs[i]), 1./3.); 155 } 156 157 const int idx_map[3][3] = {{0, 1, 2}, {1, 3, 4}, {2, 4, 5}}; 158 Mat_<double> A(6*m, 6); 159 A.setTo(0); 160 161 int eq_idx = 0; 162 for (int k = 0; k < m; ++k) 163 { 164 Mat_<double> H(Hs_[k]); 165 for (int i = 0; i < 3; ++i) 166 { 167 for (int j = i; j < 3; ++j, ++eq_idx) 168 { 169 for (int l = 0; l < 3; ++l) 170 { 171 for (int s = 0; s < 3; ++s) 172 { 173 int idx = idx_map[l][s]; 174 A(eq_idx, idx) += H(i,l) * H(j,s); 175 } 176 } 177 A(eq_idx, idx_map[i][j]) -= 1; 178 } 179 } 180 } 181 182 Mat_<double> wcoef; 183 SVD::solveZ(A, wcoef); 184 185 Mat_<double> W(3,3); 186 for (int i = 0; i < 3; ++i) 187 for (int j = i; j < 3; ++j) 188 W(i,j) = W(j,i) = wcoef(idx_map[i][j], 0) / wcoef(5,0); 189 if (!decomposeCholesky(W.ptr<double>(), W.step, 3)) 190 return false; 191 W(0,1) = W(0,2) = W(1,2) = 0; 192 K = W.t(); 193 return true; 194} 195 196} // namespace detail 197} // namespace cv 198