Blender  V3.3
panography.cc
Go to the documentation of this file.
1 // Copyright (c) 2009 libmv authors.
2 //
3 // Permission is hereby granted, free of charge, to any person obtaining a copy
4 // of this software and associated documentation files (the "Software"), to
5 // deal in the Software without restriction, including without limitation the
6 // rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
7 // sell copies of the Software, and to permit persons to whom the Software is
8 // furnished to do so, subject to the following conditions:
9 //
10 // The above copyright notice and this permission notice shall be included in
11 // all copies or substantial portions of the Software.
12 //
13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
18 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
19 // IN THE SOFTWARE.
20 //
21 
23 
24 namespace libmv {
25 
27  const Mat& x1,
28  const Mat& x2,
29  double* P) { // P must be a double[4]
30  assert(2 == x1.rows());
31  assert(2 == x1.cols());
32  assert(x1.rows() == x2.rows());
33  assert(x1.cols() == x2.cols());
34 
35  // Setup the variable of the input problem:
36  Vec xx1 = (x1.col(0)).transpose();
37  Vec yx1 = (x1.col(1)).transpose();
38 
39  double a12 = xx1.dot(yx1);
40  Vec xx2 = (x2.col(0)).transpose();
41  Vec yx2 = (x2.col(1)).transpose();
42  double b12 = xx2.dot(yx2);
43 
44  double a1 = xx1.squaredNorm();
45  double a2 = yx1.squaredNorm();
46 
47  double b1 = xx2.squaredNorm();
48  double b2 = yx2.squaredNorm();
49 
50  // Build the 3rd degre polynomial in F^2.
51  //
52  // f^6 * p + f^4 * q + f^2* r + s = 0;
53  //
54  // Coefficients in ascending powers of alpha, i.e. P[N]*x^N.
55  // Run panography_coeffs.py to get the below coefficients.
56  P[0] = b1 * b2 * a12 * a12 - a1 * a2 * b12 * b12;
57  P[1] = -2 * a1 * a2 * b12 + 2 * a12 * b1 * b2 + b1 * a12 * a12 +
58  b2 * a12 * a12 - a1 * b12 * b12 - a2 * b12 * b12;
59  P[2] = b1 * b2 - a1 * a2 - 2 * a1 * b12 - 2 * a2 * b12 + 2 * a12 * b1 +
60  2 * a12 * b2 + a12 * a12 - b12 * b12;
61  P[3] = b1 + b2 - 2 * b12 - a1 - a2 + 2 * a12;
62 
63  // If P[3] equal to 0 we get ill conditionned data
64  return (P[3] != 0.0);
65 }
66 
67 // This implements a minimal solution (2 points) for panoramic stitching:
68 //
69 // http://www.cs.ubc.ca/~mbrown/minimal/minimal.html
70 //
71 // [1] M. Brown and R. Hartley and D. Nister. Minimal Solutions for Panoramic
72 // Stitching. CVPR07.
74  const Mat& x2,
75  vector<double>* fs) {
76  // Build Polynomial factor to get squared focal value.
77  double P[4];
79 
80  // Solve it by using F = f^2 and a Cubic polynomial solver
81  //
82  // F^3 * p + F^2 * q + F^1 * r + s = 0
83  //
84  double roots[3];
85  int num_roots = SolveCubicPolynomial(P, roots);
86  for (int i = 0; i < num_roots; ++i) {
87  if (roots[i] > 0.0) {
88  fs->push_back(sqrt(roots[i]));
89  }
90  }
91 }
92 
93 // Compute the 3x3 rotation matrix that fits two 3D point clouds in the least
94 // square sense. The method is from:
95 //
96 // K. Arun,T. Huand and D. Blostein. Least-squares fitting of 2 3-D point
97 // sets. IEEE Transactions on Pattern Analysis and Machine Intelligence,
98 // 9:698-700, 1987.
99 void GetR_FixedCameraCenter(const Mat& x1,
100  const Mat& x2,
101  const double focal,
102  Mat3* R) {
103  assert(3 == x1.rows());
104  assert(2 <= x1.cols());
105  assert(x1.rows() == x2.rows());
106  assert(x1.cols() == x2.cols());
107 
108  // Build simplified K matrix
109  Mat3 K(Mat3::Identity() * 1.0 / focal);
110  K(2, 2) = 1.0;
111 
112  // Build the correlation matrix; equation (22) in [1].
113  Mat3 C = Mat3::Zero();
114  for (int i = 0; i < x1.cols(); ++i) {
115  Mat r1i = (K * x1.col(i)).normalized();
116  Mat r2i = (K * x2.col(i)).normalized();
117  C += r2i * r1i.transpose();
118  }
119 
120  // Solve for rotation. Equations (24) and (25) in [1].
121  Eigen::JacobiSVD<Mat> svd(C, Eigen::ComputeThinU | Eigen::ComputeThinV);
122  Mat3 scale = Mat3::Identity();
123  scale(2, 2) =
124  ((svd.matrixU() * svd.matrixV().transpose()).determinant() > 0.0) ? 1.0
125  : -1.0;
126 
127  (*R) = svd.matrixU() * scale * svd.matrixV().transpose();
128 }
129 
130 } // namespace libmv
sqrt(x)+1/max(0
#define K(key)
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble x2
#define C
Definition: RandGen.cpp:25
btMatrix3x3 transpose() const
Return the transpose of the matrix.
btScalar determinant() const
Return the determinant of the matrix.
Definition: btMatrix3x3.h:1022
SIMD_FORCE_INLINE btVector3 normalized() const
Return a normalized version of this vector.
static float P(float k)
Definition: math_interp.c:25
#define R
Eigen::VectorXd Vec
Definition: numeric.h:61
int SolveCubicPolynomial(Real a, Real b, Real c, Real *x0, Real *x1, Real *x2)
Definition: poly.h:39
Eigen::Matrix< double, 3, 3 > Mat3
Definition: numeric.h:72
static bool Build_Minimal2Point_PolynomialFactor(const Mat &x1, const Mat &x2, double *P)
Definition: panography.cc:26
void F_FromCorrespondance_2points(const Mat &x1, const Mat &x2, vector< double > *fs)
Definition: panography.cc:73
Eigen::MatrixXd Mat
Definition: numeric.h:60
void GetR_FixedCameraCenter(const Mat &x1, const Mat &x2, const double focal, Mat3 *R)
Definition: panography.cc:99