QVision: Qt's Image, Video and Computer Vision Library

00001 /*
00002  *      Copyright (C) 2007, 2008, 2009, 2010, 2011, 2012. PARP Research Group.
00003  *      <http://perception.inf.um.es>
00004  *      University of Murcia, Spain.
00005  *
00006  *      This file is part of the QVision library.
00007  *
00008  *      QVision is free software: you can redistribute it and/or modify
00009  *      it under the terms of the GNU Lesser General Public License as
00010  *      published by the Free Software Foundation, version 3 of the License.
00011  *
00012  *      QVision is distributed in the hope that it will be useful,
00013  *      but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  *      GNU Lesser General Public License for more details.
00016  *
00017  *      You should have received a copy of the GNU Lesser General Public
00018  *      License along with QVision. If not, see <http://www.gnu.org/licenses/>.
00019  */
00020 
00021 #ifndef QVPROJECTIVE_H
00022 #define QVPROJECTIVE_H
00023 
00024 #include <qvmath.h>
00025 #include <QVMatrix>
00026 #include <QHash>
00027 #include <QVEuclideanMapping3>
00028 #include <QV3DPointF>
00029 #include <QVCameraPose>
00030 
00031 #ifdef QVIPP
00032 #include <QVImage>
00033 #include <qvipp.h>
00034 #endif // QVIPP
00035 
00036 #include <qvmath/qvepipolar.h>
00037 #include <qvmath/qvreprojectionerror.h>
00038 
00042 
00045 
00046 void homogenizePoints(const QList< QPointFMatching > &matchings,
00047         QVMatrix &premult, QVMatrix &postmult,
00048         QList< QPair<QPointF, QPointF> > &homogeneizedPairs);
00049 
00050 /*
00051 @brief Obtains an homography from two lists of corresponding points.
00052 @ingroup qvprojectivegeometry
00053 
00054 This function returns the homography that maps the points from a source position to a
00055 destination position, according to a projective transformation.
00056 
00057 @deprecated Use @ref ComputeProjectiveHomography instead
00058 @param sourcePoints list of source points.
00059 @param destinationPoints list of destination points.
00060 */
00061 //QVMatrix ComputeHomography(const QList<QPointF> &sourcePoints, const QList<QPointF> &destPoints);
00062 
00106 bool computeProjectiveHomography(const QList< QPointFMatching > &matchings, QVMatrix &H);
00107 
00108 #ifndef DOXYGEN_IGNORE_THIS
00109 bool computeProjectiveHomography2(const QList< QPointFMatching > &matchings, QVMatrix &H);
00110 #endif // DOXYGEN_IGNORE_THIS
00111 
00120 QVMatrix computeProjectiveHomography(const QList< QPointFMatching > &matchings);
00121 
00154 QVMatrix computeAffineHomography(const QList< QPointFMatching > &matchings);
00155 
00177 QVMatrix computeSimilarHomography(const QList< QPointFMatching > &matchings);
00178 
00179 /*
00180 @brief Obtains the fundamental matrix between two images using the <a href="http://en.wikipedia.org/wiki/Eight-point_algorithm">8-point algorithm</a>.
00181 
00182 This function obtains the fundamental matrix that models the epipolar geometry between two images from a set of point correspondences.
00183 
00184 This function performs a prior normalization of the point projections, to obtain valid and better fundamental matrices.
00185 This normalization includes a correction to move the mean of the point projections to the origin of coordinates, and a posterior whitening.
00186 
00187 The code for this function is based on the 8-point algorithm implementation contained in the function <i>cvFindfundamentalMat</i> from the OpenCV.
00188 
00189 @param matchings list of 8 or more point matchings
00190 @param normalize perform normalization of the coefficient matrix rows independently
00191 @ingroup qvprojectivegeometry
00192 */
00193 #ifndef DOXYGEN_IGNORE_THIS
00194 QVMatrix computeFundamentalMatrix(const QList<QPointFMatching> &matchings, const bool normalize = false);
00195 #endif // DOXYGEN_IGNORE_THIS
00196 
00209 #ifdef QVOPENCV
00210 QVMatrix cvFindFundamentalMat(const QList< QPointFMatching > &matchings);
00211 #endif // QVOPENCV
00212 
00230 QPointF applyHomography(const QVMatrix &homography, const QPointF &point);
00231 
00232 
00247 QList<QPointF> applyHomography(const QVMatrix &homography, const QList<QPointF> &sourcePoints);
00248 
00274 #ifdef QVIPP
00275 QVImage<uChar, 1> applyHomography(const QVMatrix &homography, const QVImage<uChar, 1> &image, const int interpolation = IPPI_INTER_CUBIC);
00276 #endif // QVIPP
00277 
00302 #ifdef QVIPP
00303 QVImage<uChar, 3> applyHomography(const QVMatrix &homography, const QVImage<uChar, 3> &image, const int interpolation = IPPI_INTER_CUBIC);
00304 #endif // QVIPP
00305 
00306 /*
00307 @brief Function to test if a 3x3 matrix corresponds to an homography.
00308 
00309 Every homography matrix corresponding to a perspective deformation from one plane to another should satisfy some constrains. The most important is that its two first columns
00310 should be perpendicular, and with a similar size, because for the matrix to correspond to an homography, they should be contained in a base of a rotated coordinate system.
00311 Given the following homography matrix:
00312 
00313 \f$ H = \left[ R_1 R_2 -CR^t \right] \f$
00314 
00315 The error value returned by this function for it will be:
00316 
00317 \f$ error = \left| \frac{ \|R_1\| - \|R_2\| }{ \|R_1\| + \|R_2\| } \right| + \left| \frac{ R_1 \cdot R_2}{ \|R_1\| \|R_2\| } \right| \f$
00318 
00319 That error is a measure of the difference of sizes between the norm of the two column vectors of the homography, corresponding to the two first columns of the rotation matrix, and their dot product. When both values are close to zero, the matrix corresponds to an homography, and it won't otherwise.
00320 
00321 A good method to prove that a matrix corresponds to an homography using this function, can be done testing the return value with a threshold of approximately 0.3. If the return value of this function for a matrix is lower than this threshold, that matrix is likely to correspond to an homography, and is not likely to correspond otherwise.
00322 
00323 @param homography a possible homography matrix.
00324 @return a value close to 1 when the matrix does not corresponds to an homography, and close to 0 when it is close to be an homography.
00325 @ingroup qvprojectivegeometry
00326 */
00327 #ifndef DOXYGEN_IGNORE_THIS
00328 double HomographyTestError(const QVMatrix &homography);
00329 #endif // DOXYGEN_IGNORE_THIS
00330 
00331 /*
00332 @todo write documentation for this function
00333 
00334 @ingroup qvprojectivegeometry
00335 */
00336 #ifndef DOXYGEN_IGNORE_THIS
00337 void GetExtrinsicCameraMatrixFromHomography(const QVMatrix &K, const QVMatrix &H, QVMatrix &M4x4);
00338 #endif // DOXYGEN_IGNORE_THIS
00339 
00340 /*
00341 @brief Diagonal intrinsic camera matrix calibration
00342 
00343 
00344 This function obtains a diagonal intrinsic camera matrix, consisting on the focal distance only. This matrix is such that
00345 
00346 \f$ H = K \left[ R | t\right]   \f$
00347 
00348 Where \f$ R \f$ is a rotation matrix. This function returns a direct approximation for the \f$ K \f$ matrix.
00349 
00350 @todo rename to GetDirectCameraIntrinsicsFromPlanarHomography
00351 @see GetIntrinsicCameraMatrixFromHomography
00352 @ingroup qvprojectivegeometry
00353 */
00354 //void GetDirectIntrinsicCameraMatrixFromHomography(const QVMatrix &H, QVMatrix &K);
00355 
00356 /*
00357 @brief Diagonal intrinsic camera matrix calibration
00358 
00359 This function obtains a diagonal intrinsic camera matrix, consisting on the focal distance only. This matrix is such that
00360 
00361 \f$ H = K \left[ R | t\right]   \f$
00362 
00363 Where \f$ R \f$ is a rotation matrix. This matrix is obtained with a minimization process, so its result is more precise than that obtained with @ref GetDirectIntrinsicCameraMatrixFromHomography function.
00364 
00365 @see GetDirectIntrinsicCameraMatrixFromHomography
00366 @ingroup qvprojectivegeometry
00367 */
00368 //void GetIntrinsicCameraMatrixFromHomography(const QVMatrix &H, QVMatrix &K,
00369 //      double focal = 3, const double maxFocal = 50, const int maxIterations = 100, const double maxError = 0.00001);
00370 
00371 /*
00372 @ingroup qvprojectivegeometry
00373 @todo write documentation for this function
00374 */
00375 //void CalibrateCameraFromPlanarHomography(const QVMatrix &H, QVMatrix &K, QVMatrix &Rt);
00376 
00377 /*
00378 @brief Obtains the intrinsic camera matrix from a planar homography
00379 
00380 This functions obtains the intrinsic calibration matrix \f$ K \f$ corresponding to a simple pinhole camera model. The intrinsic camera matrix has only one free parameter, related to the focal distance of the camera:
00381 
00382 \f$
00383 K = \left(\begin{array}{ccc}
00384 f       & 0     & 0     \\
00385 0       & f     & 0     \\
00386 0       & 0     & 1     \\
00387 \end{array}\right)
00388 \f$
00389 
00390 This function should receive a planar homography corresponding to the mapping of a set of know template points, to an image frame captured with the camera containing a view of that template.
00391 
00392 The following is an example of a full intrinsic camera calibration, knowing a set of point matchings between the template image and the input image read from the camera:
00393 
00394 @code
00395 [...]
00396 QList< QPointFMatching > matchings;
00397 matchings.append(QPointFMatching(QPointF(-171,109),     QPointF(-100,+100)));
00398 matchings.append(QPointFMatching(QPointF(-120,31),      QPointF(-100,-100)));
00399 matchings.append(QPointFMatching(QPointF(117,53),       QPointF(+100,-100)));
00400 matchings.append(QPointFMatching(QPointF(11,115),       QPointF(+100,+100)));
00401 
00402 QVMatrix H, K;
00403 H = ComputeProjectiveHomography(matchings);
00404 GetPinholeCameraIntrinsicsFromPlanarHomography(H, K);
00405 [...]
00406 @endcode
00407 
00408 For further understanding of the planar homography calibration process, see @ref ComputeProjectiveHomography function documentation.
00409 
00410 @param H                        Input planar homography.
00411 @param K                        Matrix to store the intrinsic camera matrix
00412 @param iterations               Cumber of iterations to perform camera calibration
00413 @param maxGradientNorm          Minimal value of the gradient size (norm 2) to stop the minimization when reached.
00414 @param step                     Corresponds to parameter <i>step</i> for the <i>gsl_multimin_fdfminimizer_set</i> function.
00415 @param tol                      Corresponds to parameter <i>tol</i> for the <i>gsl_multimin_fdfminimizer_set</i> function.
00416 
00417 @return a value close to 1 when the matrix does not correspond to an homography, and close to 0 when it is close to be an homography.
00418 @see ComputeProjectiveHomography
00419 @ingroup qvprojectivegeometry
00420 */
00421 //void GetPinholeCameraIntrinsicsFromPlanarHomography(  const QVMatrix &H, QVMatrix &K, const int iterations = 100,
00422 //                                                      const double maxGradientNorm = 1e-3, const double step = 0.01, const double tol = 1e-4);
00423 
00424 /*
00425 @brief Obtains the canonical matrices corresponding to an essential matrix
00426 
00427 See section 9.6.2 from <i>Multiple View Geometry in Computer Vision</i>.
00428 
00429 @param E The input essential matrix.
00430 @return A list containing the possible canonical camera matrices.
00431 @ingroup qvprojectivegeometry
00432 
00433 @deprecated Use @ref getCameraPosesFromEssentialMatrix instead.
00434 */
00435 #ifndef DOXYGEN_IGNORE_THIS
00436 QList<QVMatrix> getCanonicalCameraMatricesFromEssentialMatrix(const QVMatrix &E);
00437 #endif // DOXYGEN_IGNORE_THIS
00438 
00439 // TODO: document this function.
00440 #ifndef DOXYGEN_IGNORE_THIS
00441 bool getCorrectCameraPoseTestingCheirality(const QPointFMatching matching, const QVMatrix &R1, const QVMatrix &R2, const QV3DPointF t, bool &R1IsCorrect, bool &tIsPossitive);
00442 #endif // DOXYGEN_IGNORE_THIS
00443 
00471 void getCameraPosesFromEssentialMatrix(const QVMatrix &E, QVMatrix &R1, QVMatrix &R2, QV3DPointF &t);
00472 
00490 bool testCheiralityForCameraPoses(const QVMatrix &sourceRt, const QPointF &sourceProjection, const QVMatrix &destRt, const QPointF &destProjection);
00491 
00492 /*
00493 @brief Obtains the essential matrix corresponding to a canonical camera matrix
00494 
00495 This method estimates the essential matrix of a two view scenario, provided the second camera matrix \f$ P_2 \f$, considering the first camera matrix as the identity: \f$ P_1 = \left[ I | 0 \right] \f$
00496 See section 9.6.1 from <i>Multiple View Geometry in Computer Vision</i>.
00497 
00498 @param P The input canonical matrix.
00499 @return The essential matrix corresponding to P.
00500 @ingroup qvprojectivegeometry
00501 */
00502 //QVMatrix getEssentialMatrixFromCanonicalCameraMatrix(const QVMatrix &P);
00503 
00504 #ifndef DOXYGEN_IGNORE_THIS
00505 QVMatrix getCameraMatrixFrom2D3DPointCorrespondences(const QList<QPointF> &points2d, const QList<QV3DPointF> &points3d);
00506 #endif // DOXYGEN_IGNORE_THIS
00507 
00508 #ifndef DOXYGEN_IGNORE_THIS
00509 QV3DPointF triangulate3DPointFromNViews(const QList<QPointF> &points, const QList<QVMatrix> &Plist);
00510 #endif // DOXYGEN_IGNORE_THIS
00511 
00512 #ifndef DOXYGEN_IGNORE_THIS
00513 QV3DPointF triangulate3DPointFrom2Views(const QPointF &point1, const QVMatrix &P1, const QPointF &point2, const QVMatrix &P2);
00514 #endif // DOXYGEN_IGNORE_THIS
00515 
00516 /*
00517 @brief Eliminates errors in the rotation component of a canonical camera matrix using a QR decomposition
00518 
00519 Function @ref refineExtrinsicCameraMatrixWithPolarDecomposition(const QVMatrix &) obtains a valid canonical camera
00520 matrix <i>P = [R|t]</i> from a given initial approximation <i>P* = [R*|t*]</i>, using the Polar decomposition.
00521 
00522 This function obtains a computationally more efficient approximation of the <i>P</i> matrix using a QR decomposition.
00523 The rotation matrix <i>R</i> of the resulting <i>P</i> is not the closest rotation matrix to the initial <i>R*</i>
00524 regarding the Frobenius norm \f$ ||R-R*||_F \f$.
00525 
00526 Despite of that, this function can be used when a faster and slightly less accurate version of the best <i>P</i> is wanted.
00527 
00528 @ingroup qvprojectivegeometry
00529 */
00530 //QVMatrix refineExtrinsicCameraMatrixWithQRDecomposition(const QVMatrix &P);
00531 
00532 /*
00533 @brief Eliminates errors in the rotation component of a canonical camera matrix using a Polar decomposition
00534 
00535 Given a matrix <i>P* = [R*|t*]</i> which approximates a canonical camera matrix, this function obtains the closest
00536 valid canonical matrix \f$ P = [R|t]\f$. The rotation matrix of this valid canonical matrix must be an orthogonal matrix.
00537 Thus it must fulfill:
00538 
00539         \f$ R^T R = I \f$
00540 
00541 This function uses a correcting square 3x3 matrix <i>E<sup>-1</sup></i> that satisfies the following equation:
00542 
00543         \f$ E^-1 P* = P (1) \f$
00544 
00545 This matrix <i>E<sup>-1</sup></i> also satisfies that the matrix <i>R</i> is the rotation matrix closest to <i>R*</i>
00546 in a certain sense. Using the polar decomposition of <i>R*</i>, the function obtains the <i>E</i> which corrects it
00547 to the closest valid rotation matrix <i>R</i> regarding the Frobenius norm \f$ ||R-R*||_F \f$.
00548 
00549 The polar decomposition of the transpose of matrix <i>R*</i> is the following:
00550 
00551         \f$ R* = D^T U^T \f$
00552 
00553 Where <i>U</i> is a rotation matrix and <i>D</i> is a positive-semidefinite Hermitian matrix. The matrix
00554 <i>D<sup>T</sup></i> is used as the matrix <i>E</i>.
00555 
00556 @ingroup qvprojectivegeometry
00557 */
00558 #ifndef DOXYGEN_IGNORE_THIS
00559 QVMatrix refineExtrinsicCameraMatrixWithPolarDecomposition(const QVMatrix &P);
00560 #endif // DOXYGEN_IGNORE_THIS
00561 
00562 // ---------------------------------------------------------------------------
00563 
00580 QVMatrix linearCameraResection(const QList<QPointF> &points2d, const QList<QV3DPointF> &points3d);
00581 
00582 /*
00583 @brief Obtains the center of the camera from its rotation and a set of correspondences between 3D points and their respective image projections
00584 
00585 This method provides a full camera resection by estimating the camera center, provided its as a rotation matrix \f$ R \f$. 
00586 
00587 @param R Rotation matrix for the camera orientation.
00588 @param points2d List containing the points from the image.
00589 @param points3d List containing the 3D points.
00590 @return The camera center.
00591 @ingroup qvprojectivegeometry
00592 */
00593 #ifndef DOXYGEN_IGNORE_THIS
00594 QVVector linearCameraCenterResection(const QVMatrix &R, const QList<QPointF> &points2D, const QList<QV3DPointF> &points3D);
00595 #endif // DOXYGEN_IGNORE_THIS
00596 
00608 QV3DPointF linear3DPointTriangulation(const QVector<QVMatrix> &cameraMatrices, const QHash<int, QPointF> &projectionsOfAPoint, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00609 
00627 QV3DPointF linear3DPointTriangulation(const QList<QVMatrix> &cameraMatrices, const QList<QPointF> &projectionsOfAPoint, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00628 
00649 QV3DPointF linear3DPointTriangulation(const QPointF &point1, const QVMatrix &P1, const QPointF &point2, const QVMatrix &P2, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00650 
00669 QV3DPointF linear3DPointTriangulation(const QPointFMatching &matching, const QVMatrix &P1, const QVMatrix &P2, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00670 
00689 QV3DPointF linear3DPointTriangulation(const QPointFMatching &matching, const QVCameraPose &pose1, const QVCameraPose &pose2, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00690 
00703 QList<QV3DPointF> linear3DPointsTriangulation(const QList<QVEuclideanMapping3> &cameras, const QList<QHash<int, QPointF> > &pointProjections, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00704 
00717 QList<QV3DPointF> linear3DPointsTriangulation(const QList<QVEuclideanMapping3> &cameras, const QVector<QHash<int, QPointF> > &pointProjections, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00718 
00731 QList<QV3DPointF> linear3DPointsTriangulation(const QList<QVCameraPose> &cameras, const QList<QHash<int, QPointF> > &pointProjections, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00732 
00745 QList<QV3DPointF> linear3DPointsTriangulation(const QList<QVCameraPose> &cameras, const QVector<QHash<int, QPointF> > &pointProjections, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00746 
00790 #include <qvnumericalanalysis.h>
00791 #ifdef GSL_AVAILABLE
00792 bool getCameraFocals(   const QList<QPointFMatching> &matchings,
00793                         double &focal1, double &focal2,
00794                         const QPointF principalPoint1 = QPointF(0.0, 0.0),
00795                         const QPointF principalPoint2 = QPointF(0.0, 0.0),
00796                         const GSLMultiminFDFMinimizerType gslMinimizerAlgorithm = VectorBFGS,
00797                         const int optimizationIterations = 50);
00798 #endif // GSL_AVAILABLE
00799 
00800 
00834 double computeCameraFocalFromPlanarHomography(const QVMatrix &H, int w, int h, bool byzero = false);
00835 
00844 QVCameraPose getCameraPoseFromCalibratedHomography(const QVMatrix & K, const QVMatrix & H);
00845 
00846 #ifndef DOXYGEN_IGNORE_THIS
00847 
00848 typedef enum {
00849         GEA_DO_NOT_DECOMPOSE,           
00850         GEA_CHOLESKY_DECOMPOSITION,     
00851         GEA_EIGEN_DECOMPOSITION         
00852         } TGEA_decomposition_method;
00853 
00854 QVMatrix get8PointsCoefficientMatrix(const QList<QPointFMatching> &matchings, const bool normalize = true);
00855 QVMatrix getTransposeProductOf8PointsCoefficientMatrix(const QList<QPointFMatching> &matchings, const bool normalize = true);
00856 QVMatrix getReduced8PointsCoefficientsMatrix(   // Input matching list.
00857                                                 const QList<QPointFMatching> &matchingsList,
00858                                                 // Method for decomposition.
00859                                                 const TGEA_decomposition_method decomposition_method = GEA_CHOLESKY_DECOMPOSITION,
00860                                                 // Perform pre-normalization of the point matchings
00861                                                 const bool normalize = true,
00862                                                 // GSL should be faster, but more prone to errors.
00863                                                 const bool gsl= true,
00864                                                 // Value to increase diagonal elements prior to Cholesky decomposition.
00865                                                 // Important to avoid non-positive definite matrix errors.
00866                                                 const double choleskyLambda = 0.0);
00867 
00868 QList<QPointFMatching> applyHomographies(const QVMatrix &H1, const QVMatrix &H2, const QList<QPointFMatching> &matchings);
00869 QList<QPointFMatching> correctIntrinsics(const QVMatrix &K1, const QVMatrix &K2, const QList<QPointFMatching> &matchings);
00870 QList< QHash< int, QPointF> > correctIntrinsics(const QVMatrix &K, const QList< QHash< int, QPointF> > &pointsProjections);
00871 
00872 class QVCalibrationErrorFunction: public QVFunction<QVVector, double>
00873         {
00874         private:
00875                 const QVMatrix F;
00876                 double evaluate(const QVVector &x);
00877 
00878         public:
00879                 QVCalibrationErrorFunction(const QVMatrix &F);
00880         };
00881 #endif // DOXYGEN_IGNORE_THIS
00882 
00883 // Deprecated functions
00884 #ifndef DOXYGEN_IGNORE_THIS
00885 QPointF ApplyHomography(const QVMatrix &homography, const QPointF &point);
00886 QList<QPointF> ApplyHomography(const QVMatrix &homography, const QList<QPointF> &sourcePoints);
00887 
00888 #ifdef QVIPP
00889 QVImage<uChar, 1> ApplyHomography(const QVMatrix &homography, const QVImage<uChar, 1> &image, const int interpolation = IPPI_INTER_CUBIC);
00890 QVImage<uChar, 3> ApplyHomography(const QVMatrix &homography, const QVImage<uChar, 3> &image, const int interpolation = IPPI_INTER_CUBIC);
00891 #endif
00892 
00893 QVMatrix ComputeSimilarHomography(const QList< QPointFMatching > &matchings);
00894 QVMatrix ComputeAffineHomography(const QList< QPointFMatching > &matchings);
00895 QVMatrix ComputeProjectiveHomography(const QList< QPointFMatching > &matchings);
00896 QVMatrix ComputeEuclideanHomography(const QList< QPointFMatching > &matchings);
00897 #endif // QVIPP
00898 
00899 #endif //  QVPROJECTIVE_H
src/qvmath/qvprojective.h
