diff --git a/doc/tutorials/stitching/stitcher/images/boat.jpg b/doc/tutorials/stitching/stitcher/images/boat.jpg new file mode 100644 index 0000000000..da8dd48e62 Binary files /dev/null and b/doc/tutorials/stitching/stitcher/images/boat.jpg differ diff --git a/doc/tutorials/stitching/stitcher/images/budapest.jpg b/doc/tutorials/stitching/stitcher/images/budapest.jpg new file mode 100644 index 0000000000..a4cd0d47de Binary files /dev/null and b/doc/tutorials/stitching/stitcher/images/budapest.jpg differ diff --git a/doc/tutorials/stitching/stitcher/images/newspaper.jpg b/doc/tutorials/stitching/stitcher/images/newspaper.jpg new file mode 100644 index 0000000000..eff7c51f14 Binary files /dev/null and b/doc/tutorials/stitching/stitcher/images/newspaper.jpg differ diff --git a/doc/tutorials/stitching/stitcher/stitcher.markdown b/doc/tutorials/stitching/stitcher/stitcher.markdown new file mode 100644 index 0000000000..d28bd21d70 --- /dev/null +++ b/doc/tutorials/stitching/stitcher/stitcher.markdown @@ -0,0 +1,115 @@ +High level stitching API (Stitcher class) {#tutorial_stitcher} +========================================= + +Goal +---- + +In this tutorial you will learn how to: + +- use the high-level stitching API for stitching provided by + - @ref cv::Stitcher +- learn how to use preconfigured Stitcher configurations to stitch images + using different camera models. + +Code +---- + +This tutorial code's is shown lines below. You can also download it from +[here](https://github.com/opencv/opencv/tree/master/samples/cpp/samples/cpp/stitching.cpp). + +@include samples/cpp/stitching.cpp + +Explanation +----------- + +The most important code part is: + +@code{.cpp} +Mat pano; +Ptr stitcher = Stitcher::create(mode, try_use_gpu); +Stitcher::Status status = stitcher->stitch(imgs, pano); + +if (status != Stitcher::OK) +{ + cout << "Can't stitch images, error code = " << int(status) << endl; + return -1; +} +@endcode + +A new instance of stitcher is created and the @ref cv::Stitcher::stitch will +do all the hard work. + +@ref cv::Stitcher::create can create stitcher in one of the predefined +configurations (argument `mode`). See @ref cv::Stitcher::Mode for details. These +configurations will setup multiple stitcher properties to operate in one of +predefined scenarios. After you create stitcher in one of predefined +configurations you can adjust stitching by setting any of the stitcher +properties. + +If you have cuda device @ref cv::Stitcher can be configured to offload certain +operations to GPU. If you prefer this configuration set `try_use_gpu` to true. +OpenCL acceleration will be used transparently based on global OpenCV settings +regardless of this flag. + +Stitching might fail for several reasons, you should always check if +everything went good and resulting pano is stored in `pano`. See +@ref cv::Stitcher::Status documentation for possible error codes. + +Camera models +------------- + +There are currently 2 camera models implemented in stitching pipeline. + +- _Homography model_ expecting perspective transformations between images + implemented in @ref cv::detail::BestOf2NearestMatcher cv::detail::HomographyBasedEstimator + cv::detail::BundleAdjusterReproj cv::detail::BundleAdjusterRay +- _Affine model_ expecting affine transformation with 6 DOF or 4 DOF implemented in + @ref cv::detail::AffineBestOf2NearestMatcher cv::detail::AffineBasedEstimator + cv::detail::BundleAdjusterAffine cv::detail::BundleAdjusterAffinePartial cv::AffineWarper + +Homography model is useful for creating photo panoramas captured by camera, +while affine-based model can be used to stitch scans and object captured by +specialized devices. + +@note +Certain detailed settings of @ref cv::Stitcher might not make sense. Especially +you should not mix classes implementing affine model and classes implementing +Homography model, as they work with different transformations. + +Try it out +---------- + +If you enabled building samples you can found binary under +`build/bin/cpp-example-stitching`. This example is a console application, run it without +arguments to see help. `opencv_extra` provides some sample data for testing all available +configurations. + +to try panorama mode run: +``` +./cpp-example-stitching --mode panorama /testdata/stitching/boat* +``` +![](images/boat.jpg) + +to try scans mode run (dataset from home-grade scanner): +``` +./cpp-example-stitching --mode scans /testdata/stitching/newspaper* +``` +![](images/newspaper.jpg) + +or (dataset from professional book scanner): +``` +./cpp-example-stitching --mode scans /testdata/stitching/budapest* +``` +![](images/budapest.jpg) + +@note +Examples above expects POSIX platform, on windows you have to provide all files names explicitly +(e.g. `boat1.jpg` `boat2.jpg`...) as windows command line does not support `*` expansion. + +See also +-------- + +If you want to study internals of the stitching pipeline or you want to experiment with detailed +configuration see +[stitching_detailed.cpp](https://github.com/opencv/opencv/tree/master/samples/cpp/stitching_detailed.cpp) +in `opencv/samples/cpp` folder. diff --git a/doc/tutorials/stitching/table_of_content_stitching.markdown b/doc/tutorials/stitching/table_of_content_stitching.markdown new file mode 100644 index 0000000000..d85571cd7e --- /dev/null +++ b/doc/tutorials/stitching/table_of_content_stitching.markdown @@ -0,0 +1,15 @@ +Images stitching (stitching module) {#tutorial_table_of_content_stitching} +=================================== + +Sometimes a single image can't capture it all. Here you will learn how to join +more images together to create a large pano. Doesn't matter if you want to +create a photo panorama or you want to stitch scans. + +- @subpage tutorial_stitcher + + *Compatibility:* \>= OpenCV 3.2 + + *Author:* Jiri Horner + + You will use high level stitching api to create a photo panorama. You will + learn about Stitcher class and its configurations. diff --git a/doc/tutorials/tutorials.markdown b/doc/tutorials/tutorials.markdown index 9319217772..1f1a7e9f98 100644 --- a/doc/tutorials/tutorials.markdown +++ b/doc/tutorials/tutorials.markdown @@ -68,6 +68,10 @@ As always, we would be happy to hear your comments and receive your contribution Use OpenCV for advanced photo processing. +- @subpage tutorial_table_of_content_stitching + + Learn how to create beautiful photo panoramas and more with OpenCV stitching pipeline. + - @subpage tutorial_table_of_content_gpu Squeeze out every diff --git a/modules/calib3d/include/opencv2/calib3d.hpp b/modules/calib3d/include/opencv2/calib3d.hpp index e5d3b1cc74..6c61e6a854 100644 --- a/modules/calib3d/include/opencv2/calib3d.hpp +++ b/modules/calib3d/include/opencv2/calib3d.hpp @@ -360,8 +360,9 @@ determined up to a scale. Thus, it is normalized so that \f$h_{33}=1\f$. Note th cannot be estimated, an empty one will be returned. @sa - getAffineTransform, getPerspectiveTransform, estimateRigidTransform, warpPerspective, - perspectiveTransform +getAffineTransform, estimateAffine2D, estimateAffinePartial2D, getPerspectiveTransform, warpPerspective, +perspectiveTransform + @note - A example on calculating a homography for image matching can be found at @@ -1583,6 +1584,93 @@ CV_EXPORTS_W int estimateAffine3D(InputArray src, InputArray dst, OutputArray out, OutputArray inliers, double ransacThreshold = 3, double confidence = 0.99); +/** @brief Computes an optimal affine transformation between two 2D point sets. + +@param from First input 2D point set. +@param to Second input 2D point set. +@param inliers Output vector indicating which points are inliers. +@param method Robust method used to compute tranformation. The following methods are possible: +- cv::RANSAC - RANSAC-based robust method +- cv::LMEDS - Least-Median robust method +RANSAC is the default method. +@param ransacReprojThreshold Maximum reprojection error in the RANSAC algorithm to consider +a point as an inlier. Applies only to RANSAC. +@param maxIters The maximum number of robust method iterations, 2000 is the maximum it can be. +@param confidence Confidence level, between 0 and 1, for the estimated transformation. Anything +between 0.95 and 0.99 is usually good enough. Values too close to 1 can slow down the estimation +significantly. Values lower than 0.8-0.9 can result in an incorrectly estimated transformation. +@param refineIters Maximum number of iterations of refining algorithm (Levenberg-Marquardt). +Passing 0 will disable refining, so the output matrix will be output of robust method. + +@return Output 2D affine transformation matrix \f$2 \times 3\f$ or empty matrix if transformation +could not be estimated. + +The function estimates an optimal 2D affine transformation between two 2D point sets using the +selected robust algorithm. + +The computed transformation is then refined further (using only inliers) with the +Levenberg-Marquardt method to reduce the re-projection error even more. + +@note +The RANSAC method can handle practically any ratio of outliers but need a threshold to +distinguish inliers from outliers. The method LMeDS does not need any threshold but it works +correctly only when there are more than 50% of inliers. + +@sa estimateAffinePartial2D, getAffineTransform +*/ +CV_EXPORTS_W cv::Mat estimateAffine2D(InputArray from, InputArray to, OutputArray inliers = noArray(), + int method = RANSAC, double ransacReprojThreshold = 3, + size_t maxIters = 2000, double confidence = 0.99, + size_t refineIters = 10); + +/** @brief Computes an optimal limited affine transformation with 4 degrees of freedom between +two 2D point sets. + +@param from First input 2D point set. +@param to Second input 2D point set. +@param inliers Output vector indicating which points are inliers. +@param method Robust method used to compute tranformation. The following methods are possible: +- cv::RANSAC - RANSAC-based robust method +- cv::LMEDS - Least-Median robust method +RANSAC is the default method. +@param ransacReprojThreshold Maximum reprojection error in the RANSAC algorithm to consider +a point as an inlier. Applies only to RANSAC. +@param maxIters The maximum number of robust method iterations, 2000 is the maximum it can be. +@param confidence Confidence level, between 0 and 1, for the estimated transformation. Anything +between 0.95 and 0.99 is usually good enough. Values too close to 1 can slow down the estimation +significantly. Values lower than 0.8-0.9 can result in an incorrectly estimated transformation. +@param refineIters Maximum number of iterations of refining algorithm (Levenberg-Marquardt). +Passing 0 will disable refining, so the output matrix will be output of robust method. + +@return Output 2D affine transformation (4 degrees of freedom) matrix \f$2 \times 3\f$ or +empty matrix if transformation could not be estimated. + +The function estimates an optimal 2D affine transformation with 4 degrees of freedom limited to +combinations of translation, rotation, and uniform scaling. Uses the selected algorithm for robust +estimation. + +The computed transformation is then refined further (using only inliers) with the +Levenberg-Marquardt method to reduce the re-projection error even more. + +Estimated transformation matrix is: +\f[ \begin{bmatrix} \cos(\theta)s & -\sin(\theta)s & tx \\ + \sin(\theta)s & \cos(\theta)s & ty +\end{bmatrix} \f] +Where \f$ \theta \f$ is the rotation angle, \f$ s \f$ the scaling factor and \f$ tx, ty \f$ are +translations in \f$ x, y \f$ axes respectively. + +@note +The RANSAC method can handle practically any ratio of outliers but need a threshold to +distinguish inliers from outliers. The method LMeDS does not need any threshold but it works +correctly only when there are more than 50% of inliers. + +@sa estimateAffine2D, getAffineTransform +*/ +CV_EXPORTS_W cv::Mat estimateAffinePartial2D(InputArray from, InputArray to, OutputArray inliers = noArray(), + int method = RANSAC, double ransacReprojThreshold = 3, + size_t maxIters = 2000, double confidence = 0.99, + size_t refineIters = 10); + /** @brief Decompose a homography matrix to rotation(s), translation(s) and plane normal(s). @param H The input homography matrix between two images. diff --git a/modules/calib3d/perf/perf_affine2d.cpp b/modules/calib3d/perf/perf_affine2d.cpp new file mode 100644 index 0000000000..b893c39840 --- /dev/null +++ b/modules/calib3d/perf/perf_affine2d.cpp @@ -0,0 +1,170 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// (3-clause BSD License) +// +// Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the names of the copyright holders nor the names of the contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall copyright holders or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#include "perf_precomp.hpp" +#include +#include + +namespace cvtest +{ + +using std::tr1::tuple; +using std::tr1::get; +using namespace perf; +using namespace testing; +using namespace cv; + +CV_ENUM(Method, RANSAC, LMEDS) +typedef tuple AffineParams; +typedef TestBaseWithParam EstimateAffine; +#define ESTIMATE_PARAMS Combine(Values(100000, 5000, 100), Values(0.99, 0.95, 0.9), Method::all(), Values(10, 0)) + +static float rngIn(float from, float to) { return from + (to-from) * (float)theRNG(); } + +static Mat rngPartialAffMat() { + double theta = rngIn(0, (float)CV_PI*2.f); + double scale = rngIn(0, 3); + double tx = rngIn(-2, 2); + double ty = rngIn(-2, 2); + double aff[2*3] = { std::cos(theta) * scale, -std::sin(theta) * scale, tx, + std::sin(theta) * scale, std::cos(theta) * scale, ty }; + return Mat(2, 3, CV_64F, aff).clone(); +} + +PERF_TEST_P( EstimateAffine, EstimateAffine2D, ESTIMATE_PARAMS ) +{ + AffineParams params = GetParam(); + const int n = get<0>(params); + const double confidence = get<1>(params); + const int method = get<2>(params); + const size_t refining = get<3>(params); + + Mat aff(2, 3, CV_64F); + cv::randu(aff, -2., 2.); + + // LMEDS can't handle more than 50% outliers (by design) + int m; + if (method == LMEDS) + m = 3*n/5; + else + m = 2*n/5; + const float shift_outl = 15.f; + const float noise_level = 20.f; + + Mat fpts(1, n, CV_32FC2); + Mat tpts(1, n, CV_32FC2); + + randu(fpts, 0., 100.); + transform(fpts, tpts, aff); + + /* adding noise to some points */ + Mat outliers = tpts.colRange(m, n); + outliers.reshape(1) += shift_outl; + + Mat noise (outliers.size(), outliers.type()); + randu(noise, 0., noise_level); + outliers += noise; + + Mat aff_est; + vector inliers (n); + + warmup(inliers, WARMUP_WRITE); + warmup(fpts, WARMUP_READ); + warmup(tpts, WARMUP_READ); + + TEST_CYCLE() + { + aff_est = estimateAffine2D(fpts, tpts, inliers, method, 3, 2000, confidence, refining); + } + + // we already have accuracy tests + SANITY_CHECK_NOTHING(); +} + +PERF_TEST_P( EstimateAffine, EstimateAffinePartial2D, ESTIMATE_PARAMS ) +{ + AffineParams params = GetParam(); + const int n = get<0>(params); + const double confidence = get<1>(params); + const int method = get<2>(params); + const size_t refining = get<3>(params); + + Mat aff = rngPartialAffMat(); + + int m; + // LMEDS can't handle more than 50% outliers (by design) + if (method == LMEDS) + m = 3*n/5; + else + m = 2*n/5; + const float shift_outl = 15.f; const float noise_level = 20.f; + + Mat fpts(1, n, CV_32FC2); + Mat tpts(1, n, CV_32FC2); + + randu(fpts, 0., 100.); + transform(fpts, tpts, aff); + + /* adding noise*/ + Mat outliers = tpts.colRange(m, n); + outliers.reshape(1) += shift_outl; + + Mat noise (outliers.size(), outliers.type()); + randu(noise, 0., noise_level); + outliers += noise; + + Mat aff_est; + vector inliers (n); + + warmup(inliers, WARMUP_WRITE); + warmup(fpts, WARMUP_READ); + warmup(tpts, WARMUP_READ); + + TEST_CYCLE() + { + aff_est = estimateAffinePartial2D(fpts, tpts, inliers, method, 3, 2000, confidence, refining); + } + + // we already have accuracy tests + SANITY_CHECK_NOTHING(); +} + +} // namespace cvtest diff --git a/modules/calib3d/src/fundam.cpp b/modules/calib3d/src/fundam.cpp index cdcbeadd77..9af70dd210 100644 --- a/modules/calib3d/src/fundam.cpp +++ b/modules/calib3d/src/fundam.cpp @@ -47,29 +47,6 @@ namespace cv { -static bool haveCollinearPoints( const Mat& m, int count ) -{ - int j, k, i = count-1; - const Point2f* ptr = m.ptr(); - - // check that the i-th selected point does not belong - // to a line connecting some previously selected points - for( j = 0; j < i; j++ ) - { - double dx1 = ptr[j].x - ptr[i].x; - double dy1 = ptr[j].y - ptr[i].y; - for( k = 0; k < j; k++ ) - { - double dx2 = ptr[k].x - ptr[i].x; - double dy2 = ptr[k].y - ptr[i].y; - if( fabs(dx2*dy1 - dy2*dx1) <= FLT_EPSILON*(fabs(dx1) + fabs(dy1) + fabs(dx2) + fabs(dy2))) - return true; - } - } - return false; -} - - class HomographyEstimatorCallback : public PointSetRegistrator::Callback { public: diff --git a/modules/calib3d/src/precomp.hpp b/modules/calib3d/src/precomp.hpp index 0005f234e2..3208680c5f 100644 --- a/modules/calib3d/src/precomp.hpp +++ b/modules/calib3d/src/precomp.hpp @@ -115,8 +115,31 @@ template inline int compressElems( T* ptr, const uchar* mask, int ms return j; } +static inline bool haveCollinearPoints( const Mat& m, int count ) +{ + int j, k, i = count-1; + const Point2f* ptr = m.ptr(); + + // check that the i-th selected point does not belong + // to a line connecting some previously selected points + // also checks that points are not too close to each other + for( j = 0; j < i; j++ ) + { + double dx1 = ptr[j].x - ptr[i].x; + double dy1 = ptr[j].y - ptr[i].y; + for( k = 0; k < j; k++ ) + { + double dx2 = ptr[k].x - ptr[i].x; + double dy2 = ptr[k].y - ptr[i].y; + if( fabs(dx2*dy1 - dy2*dx1) <= FLT_EPSILON*(fabs(dx1) + fabs(dy1) + fabs(dx2) + fabs(dy2))) + return true; + } + } + return false; } +} // namespace cv + int checkChessboard(const cv::Mat & img, const cv::Size & size); int checkChessboardBinary(const cv::Mat & img, const cv::Size & size); diff --git a/modules/calib3d/src/ptsetreg.cpp b/modules/calib3d/src/ptsetreg.cpp index 5e652d4fd5..cbf8175d48 100644 --- a/modules/calib3d/src/ptsetreg.cpp +++ b/modules/calib3d/src/ptsetreg.cpp @@ -499,11 +499,274 @@ public: } }; -} +class Affine2DEstimatorCallback : public PointSetRegistrator::Callback +{ +public: + int runKernel( InputArray _m1, InputArray _m2, OutputArray _model ) const + { + Mat m1 = _m1.getMat(), m2 = _m2.getMat(); + const Point2f* from = m1.ptr(); + const Point2f* to = m2.ptr(); + _model.create(2, 3, CV_64F); + Mat M_mat = _model.getMat(); + double *M = M_mat.ptr(); -int cv::estimateAffine3D(InputArray _from, InputArray _to, - OutputArray _out, OutputArray _inliers, - double param1, double param2) + // we need 3 points to estimate affine transform + double x1 = from[0].x; + double y1 = from[0].y; + double x2 = from[1].x; + double y2 = from[1].y; + double x3 = from[2].x; + double y3 = from[2].y; + + double X1 = to[0].x; + double Y1 = to[0].y; + double X2 = to[1].x; + double Y2 = to[1].y; + double X3 = to[2].x; + double Y3 = to[2].y; + + /* + We want to solve AX = B + + | x1 y1 1 0 0 0 | + | 0 0 0 x1 y1 1 | + | x2 y2 1 0 0 0 | + A = | 0 0 0 x2 y2 1 | + | x3 y3 1 0 0 0 | + | 0 0 0 x3 y3 1 | + B = (X1, Y1, X2, Y2, X3, Y3).t() + X = (a, b, c, d, e, f).t() + + As the estimate of (a, b, c) only depends on the Xi, and (d, e, f) only + depends on the Yi, we do the *trick* to solve each one analytically. + + | X1 | | x1 y1 1 | | a | + | X2 | = | x2 y2 1 | * | b | + | X3 | | x3 y3 1 | | c | + + | Y1 | | x1 y1 1 | | d | + | Y2 | = | x2 y2 1 | * | e | + | Y3 | | x3 y3 1 | | f | + */ + + double d = 1. / ( x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2) ); + + M[0] = d * ( X1*(y2-y3) + X2*(y3-y1) + X3*(y1-y2) ); + M[1] = d * ( X1*(x3-x2) + X2*(x1-x3) + X3*(x2-x1) ); + M[2] = d * ( X1*(x2*y3 - x3*y2) + X2*(x3*y1 - x1*y3) + X3*(x1*y2 - x2*y1) ); + + M[3] = d * ( Y1*(y2-y3) + Y2*(y3-y1) + Y3*(y1-y2) ); + M[4] = d * ( Y1*(x3-x2) + Y2*(x1-x3) + Y3*(x2-x1) ); + M[5] = d * ( Y1*(x2*y3 - x3*y2) + Y2*(x3*y1 - x1*y3) + Y3*(x1*y2 - x2*y1) ); + return 1; + } + + void computeError( InputArray _m1, InputArray _m2, InputArray _model, OutputArray _err ) const + { + Mat m1 = _m1.getMat(), m2 = _m2.getMat(), model = _model.getMat(); + const Point2f* from = m1.ptr(); + const Point2f* to = m2.ptr(); + const double* F = model.ptr(); + + int count = m1.checkVector(2); + CV_Assert( count > 0 ); + + _err.create(count, 1, CV_32F); + Mat err = _err.getMat(); + float* errptr = err.ptr(); + // transform matrix to floats + float F0 = (float)F[0], F1 = (float)F[1], F2 = (float)F[2]; + float F3 = (float)F[3], F4 = (float)F[4], F5 = (float)F[5]; + + for(int i = 0; i < count; i++ ) + { + const Point2f& f = from[i]; + const Point2f& t = to[i]; + + float a = F0*f.x + F1*f.y + F2 - t.x; + float b = F3*f.x + F4*f.y + F5 - t.y; + + errptr[i] = a*a + b*b; + } + } + + bool checkSubset( InputArray _ms1, InputArray, int count ) const + { + Mat ms1 = _ms1.getMat(); + // check colinearity and also check that points are too close + // only ms1 affects actual estimation stability + return !haveCollinearPoints(ms1, count); + } +}; + +class AffinePartial2DEstimatorCallback : public Affine2DEstimatorCallback +{ +public: + int runKernel( InputArray _m1, InputArray _m2, OutputArray _model ) const + { + Mat m1 = _m1.getMat(), m2 = _m2.getMat(); + const Point2f* from = m1.ptr(); + const Point2f* to = m2.ptr(); + _model.create(2, 3, CV_64F); + Mat M_mat = _model.getMat(); + double *M = M_mat.ptr(); + + // we need only 2 points to estimate transform + double x1 = from[0].x; + double y1 = from[0].y; + double x2 = from[1].x; + double y2 = from[1].y; + + double X1 = to[0].x; + double Y1 = to[0].y; + double X2 = to[1].x; + double Y2 = to[1].y; + + /* + we are solving AS = B + | x1 -y1 1 0 | + | y1 x1 0 1 | + A = | x2 -y2 1 0 | + | y2 x2 0 1 | + B = (X1, Y1, X2, Y2).t() + we solve that analytically + */ + double d = 1./((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + + // solution vector + double S0 = d * ( (X1-X2)*(x1-x2) + (Y1-Y2)*(y1-y2) ); + double S1 = d * ( (Y1-Y2)*(x1-x2) - (X1-X2)*(y1-y2) ); + double S2 = d * ( (Y1-Y2)*(x1*y2 - x2*y1) - (X1*y2 - X2*y1)*(y1-y2) - (X1*x2 - X2*x1)*(x1-x2) ); + double S3 = d * (-(X1-X2)*(x1*y2 - x2*y1) - (Y1*x2 - Y2*x1)*(x1-x2) - (Y1*y2 - Y2*y1)*(y1-y2) ); + + // set model, rotation part is antisymmetric + M[0] = M[4] = S0; + M[1] = -S1; + M[2] = S2; + M[3] = S1; + M[5] = S3; + return 1; + } +}; + +class Affine2DRefineCallback : public LMSolver::Callback +{ +public: + Affine2DRefineCallback(InputArray _src, InputArray _dst) + { + src = _src.getMat(); + dst = _dst.getMat(); + } + + bool compute(InputArray _param, OutputArray _err, OutputArray _Jac) const + { + int i, count = src.checkVector(2); + Mat param = _param.getMat(); + _err.create(count*2, 1, CV_64F); + Mat err = _err.getMat(), J; + if( _Jac.needed()) + { + _Jac.create(count*2, param.rows, CV_64F); + J = _Jac.getMat(); + CV_Assert( J.isContinuous() && J.cols == 6 ); + } + + const Point2f* M = src.ptr(); + const Point2f* m = dst.ptr(); + const double* h = param.ptr(); + double* errptr = err.ptr(); + double* Jptr = J.data ? J.ptr() : 0; + + for( i = 0; i < count; i++ ) + { + double Mx = M[i].x, My = M[i].y; + double xi = h[0]*Mx + h[1]*My + h[2]; + double yi = h[3]*Mx + h[4]*My + h[5]; + errptr[i*2] = xi - m[i].x; + errptr[i*2+1] = yi - m[i].y; + + /* + Jacobian should be: + {x, y, 1, 0, 0, 0} + {0, 0, 0, x, y, 1} + */ + if( Jptr ) + { + Jptr[0] = Mx; Jptr[1] = My; Jptr[2] = 1.; + Jptr[3] = Jptr[4] = Jptr[5] = 0.; + Jptr[6] = Jptr[7] = Jptr[8] = 0.; + Jptr[9] = Mx; Jptr[10] = My; Jptr[11] = 1.; + + Jptr += 6*2; + } + } + + return true; + } + + Mat src, dst; +}; + +class AffinePartial2DRefineCallback : public LMSolver::Callback +{ +public: + AffinePartial2DRefineCallback(InputArray _src, InputArray _dst) + { + src = _src.getMat(); + dst = _dst.getMat(); + } + + bool compute(InputArray _param, OutputArray _err, OutputArray _Jac) const + { + int i, count = src.checkVector(2); + Mat param = _param.getMat(); + _err.create(count*2, 1, CV_64F); + Mat err = _err.getMat(), J; + if( _Jac.needed()) + { + _Jac.create(count*2, param.rows, CV_64F); + J = _Jac.getMat(); + CV_Assert( J.isContinuous() && J.cols == 4 ); + } + + const Point2f* M = src.ptr(); + const Point2f* m = dst.ptr(); + const double* h = param.ptr(); + double* errptr = err.ptr(); + double* Jptr = J.data ? J.ptr() : 0; + + for( i = 0; i < count; i++ ) + { + double Mx = M[i].x, My = M[i].y; + double xi = h[0]*Mx - h[1]*My + h[2]; + double yi = h[1]*Mx + h[0]*My + h[3]; + errptr[i*2] = xi - m[i].x; + errptr[i*2+1] = yi - m[i].y; + + /* + Jacobian should be: + {x, -y, 1, 0} + {y, x, 0, 1} + */ + if( Jptr ) + { + Jptr[0] = Mx; Jptr[1] = -My; Jptr[2] = 1.; Jptr[3] = 0.; + Jptr[4] = My; Jptr[5] = Mx; Jptr[6] = 0.; Jptr[7] = 1.; + + Jptr += 4*2; + } + } + + return true; + } + + Mat src, dst; +}; + +int estimateAffine3D(InputArray _from, InputArray _to, + OutputArray _out, OutputArray _inliers, + double param1, double param2) { CV_INSTRUMENT_REGION() @@ -524,3 +787,152 @@ int cv::estimateAffine3D(InputArray _from, InputArray _to, return createRANSACPointSetRegistrator(makePtr(), 4, param1, param2)->run(dFrom, dTo, _out, _inliers); } + +Mat estimateAffine2D(InputArray _from, InputArray _to, OutputArray _inliers, + const int method, const double ransacReprojThreshold, + const size_t maxIters, const double confidence, + const size_t refineIters) +{ + Mat from = _from.getMat(), to = _to.getMat(); + int count = from.checkVector(2); + bool result = false; + Mat H; + + CV_Assert( count >= 0 && to.checkVector(2) == count ); + + if (from.type() != CV_32FC2 || to.type() != CV_32FC2) + { + Mat tmp; + from.convertTo(tmp, CV_32FC2); + from = tmp; + to.convertTo(tmp, CV_32FC2); + to = tmp; + } + // convert to N x 1 vectors + from = from.reshape(2, count); + to = to.reshape(2, count); + + Mat inliers; + if(_inliers.needed()) + { + _inliers.create(count, 1, CV_8U, -1, true); + inliers = _inliers.getMat(); + } + + // run robust method + Ptr cb = makePtr(); + if( method == RANSAC ) + result = createRANSACPointSetRegistrator(cb, 3, ransacReprojThreshold, confidence, static_cast(maxIters))->run(from, to, H, inliers); + else if( method == LMEDS ) + result = createLMeDSPointSetRegistrator(cb, 3, confidence, static_cast(maxIters))->run(from, to, H, inliers); + else + CV_Error(Error::StsBadArg, "Unknown or unsupported robust estimation method"); + + if(result && count > 3 && refineIters) + { + // reorder to start with inliers + compressElems(from.ptr(), inliers.ptr(), 1, count); + int inliers_count = compressElems(to.ptr(), inliers.ptr(), 1, count); + if(inliers_count > 0) + { + Mat src = from.rowRange(0, inliers_count); + Mat dst = to.rowRange(0, inliers_count); + Mat Hvec = H.reshape(1, 6); + createLMSolver(makePtr(src, dst), static_cast(refineIters))->run(Hvec); + } + } + + if (!result) + { + H.release(); + if(_inliers.needed()) + { + inliers = Mat::zeros(count, 1, CV_8U); + inliers.copyTo(_inliers); + } + } + + return H; +} + +Mat estimateAffinePartial2D(InputArray _from, InputArray _to, OutputArray _inliers, + const int method, const double ransacReprojThreshold, + const size_t maxIters, const double confidence, + const size_t refineIters) +{ + Mat from = _from.getMat(), to = _to.getMat(); + const int count = from.checkVector(2); + bool result = false; + Mat H; + + CV_Assert( count >= 0 && to.checkVector(2) == count ); + + if (from.type() != CV_32FC2 || to.type() != CV_32FC2) + { + Mat tmp; + from.convertTo(tmp, CV_32FC2); + from = tmp; + to.convertTo(tmp, CV_32FC2); + to = tmp; + } + // convert to N x 1 vectors + from = from.reshape(2, count); + to = to.reshape(2, count); + + Mat inliers; + if(_inliers.needed()) + { + _inliers.create(count, 1, CV_8U, -1, true); + inliers = _inliers.getMat(); + } + + // run robust estimation + Ptr cb = makePtr(); + if( method == RANSAC ) + result = createRANSACPointSetRegistrator(cb, 2, ransacReprojThreshold, confidence, static_cast(maxIters))->run(from, to, H, inliers); + else if( method == LMEDS ) + result = createLMeDSPointSetRegistrator(cb, 2, confidence, static_cast(maxIters))->run(from, to, H, inliers); + else + CV_Error(Error::StsBadArg, "Unknown or unsupported robust estimation method"); + + if(result && count > 2 && refineIters) + { + // reorder to start with inliers + compressElems(from.ptr(), inliers.ptr(), 1, count); + int inliers_count = compressElems(to.ptr(), inliers.ptr(), 1, count); + if(inliers_count > 0) + { + Mat src = from.rowRange(0, inliers_count); + Mat dst = to.rowRange(0, inliers_count); + // H is + // a -b tx + // b a ty + // Hvec model for LevMarq is + // (a, b, tx, ty) + double *Hptr = H.ptr(); + double Hvec_buf[4] = {Hptr[0], Hptr[3], Hptr[2], Hptr[5]}; + Mat Hvec (4, 1, CV_64F, Hvec_buf); + createLMSolver(makePtr(src, dst), static_cast(refineIters))->run(Hvec); + // update H with refined parameters + Hptr[0] = Hptr[4] = Hvec_buf[0]; + Hptr[1] = -Hvec_buf[1]; + Hptr[2] = Hvec_buf[2]; + Hptr[3] = Hvec_buf[1]; + Hptr[5] = Hvec_buf[3]; + } + } + + if (!result) + { + H.release(); + if(_inliers.needed()) + { + inliers = Mat::zeros(count, 1, CV_8U); + inliers.copyTo(_inliers); + } + } + + return H; +} + +} // namespace cv diff --git a/modules/calib3d/test/test_affine2d_estimator.cpp b/modules/calib3d/test/test_affine2d_estimator.cpp new file mode 100644 index 0000000000..de9f7003e1 --- /dev/null +++ b/modules/calib3d/test/test_affine2d_estimator.cpp @@ -0,0 +1,130 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// (3-clause BSD License) +// +// Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the names of the copyright holders nor the names of the contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall copyright holders or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#include "test_precomp.hpp" + +using namespace cv; +using namespace std; +using namespace testing; + +#include +#include + +CV_ENUM(Method, RANSAC, LMEDS) +typedef TestWithParam EstimateAffine2D; + +static float rngIn(float from, float to) { return from + (to-from) * (float)theRNG(); } + +TEST_P(EstimateAffine2D, test3Points) +{ + // try more transformations + for (size_t i = 0; i < 500; ++i) + { + Mat aff(2, 3, CV_64F); + cv::randu(aff, 1., 3.); + + Mat fpts(1, 3, CV_32FC2); + Mat tpts(1, 3, CV_32FC2); + + // setting points that are not in the same line + fpts.at(0) = Point2f( rngIn(1,2), rngIn(5,6) ); + fpts.at(1) = Point2f( rngIn(3,4), rngIn(3,4) ); + fpts.at(2) = Point2f( rngIn(1,2), rngIn(3,4) ); + + transform(fpts, tpts, aff); + + vector inliers; + Mat aff_est = estimateAffine2D(fpts, tpts, inliers, GetParam() /* method */); + + EXPECT_NEAR(0., cvtest::norm(aff_est, aff, NORM_INF), 1e-3); + + // all must be inliers + EXPECT_EQ(countNonZero(inliers), 3); + } +} + +TEST_P(EstimateAffine2D, testNPoints) +{ + // try more transformations + for (size_t i = 0; i < 500; ++i) + { + Mat aff(2, 3, CV_64F); + cv::randu(aff, -2., 2.); + const int method = GetParam(); + const int n = 100; + int m; + // LMEDS can't handle more than 50% outliers (by design) + if (method == LMEDS) + m = 3*n/5; + else + m = 2*n/5; + const float shift_outl = 15.f; + const float noise_level = 20.f; + + Mat fpts(1, n, CV_32FC2); + Mat tpts(1, n, CV_32FC2); + + randu(fpts, 0., 100.); + transform(fpts, tpts, aff); + + /* adding noise to some points */ + Mat outliers = tpts.colRange(m, n); + outliers.reshape(1) += shift_outl; + + Mat noise (outliers.size(), outliers.type()); + randu(noise, 0., noise_level); + outliers += noise; + + vector inliers; + Mat aff_est = estimateAffine2D(fpts, tpts, inliers, method); + + EXPECT_FALSE(aff_est.empty()) << "estimation failed, unable to estimate transform"; + + EXPECT_NEAR(0., cvtest::norm(aff_est, aff, NORM_INF), 1e-4); + + bool inliers_good = count(inliers.begin(), inliers.end(), 1) == m && + m == accumulate(inliers.begin(), inliers.begin() + m, 0); + + EXPECT_TRUE(inliers_good); + } +} + +INSTANTIATE_TEST_CASE_P(Calib3d, EstimateAffine2D, Method::all()); diff --git a/modules/calib3d/test/test_affine3d_estimator.cpp b/modules/calib3d/test/test_affine3d_estimator.cpp index 26088157d9..aa41bf39e5 100644 --- a/modules/calib3d/test/test_affine3d_estimator.cpp +++ b/modules/calib3d/test/test_affine3d_estimator.cpp @@ -194,4 +194,4 @@ void CV_Affine3D_EstTest::run( int /* start_from */) ts->set_failed_test_info(cvtest::TS::OK); } -TEST(Calib3d_EstimateAffineTransform, accuracy) { CV_Affine3D_EstTest test; test.safe_run(); } +TEST(Calib3d_EstimateAffine3D, accuracy) { CV_Affine3D_EstTest test; test.safe_run(); } diff --git a/modules/calib3d/test/test_affine_partial2d_estimator.cpp b/modules/calib3d/test/test_affine_partial2d_estimator.cpp new file mode 100644 index 0000000000..dde7d7da1c --- /dev/null +++ b/modules/calib3d/test/test_affine_partial2d_estimator.cpp @@ -0,0 +1,140 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// (3-clause BSD License) +// +// Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the names of the copyright holders nor the names of the contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall copyright holders or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#include "test_precomp.hpp" + +using namespace cv; +using namespace std; +using namespace testing; + +#include +#include + +CV_ENUM(Method, RANSAC, LMEDS) +typedef TestWithParam EstimateAffinePartial2D; + +static float rngIn(float from, float to) { return from + (to-from) * (float)theRNG(); } + +// get random matrix of affine transformation limited to combinations of translation, +// rotation, and uniform scaling +static Mat rngPartialAffMat() { + double theta = rngIn(0, (float)CV_PI*2.f); + double scale = rngIn(0, 3); + double tx = rngIn(-2, 2); + double ty = rngIn(-2, 2); + double aff[2*3] = { std::cos(theta) * scale, -std::sin(theta) * scale, tx, + std::sin(theta) * scale, std::cos(theta) * scale, ty }; + return Mat(2, 3, CV_64F, aff).clone(); +} + +TEST_P(EstimateAffinePartial2D, test2Points) +{ + // try more transformations + for (size_t i = 0; i < 500; ++i) + { + Mat aff = rngPartialAffMat(); + + // setting points that are no in the same line + Mat fpts(1, 2, CV_32FC2); + Mat tpts(1, 2, CV_32FC2); + + fpts.at(0) = Point2f( rngIn(1,2), rngIn(5,6) ); + fpts.at(1) = Point2f( rngIn(3,4), rngIn(3,4) ); + + transform(fpts, tpts, aff); + + vector inliers; + Mat aff_est = estimateAffinePartial2D(fpts, tpts, inliers, GetParam() /* method */); + + EXPECT_NEAR(0., cvtest::norm(aff_est, aff, NORM_INF), 1e-3); + + // all must be inliers + EXPECT_EQ(countNonZero(inliers), 2); + } +} + +TEST_P(EstimateAffinePartial2D, testNPoints) +{ + // try more transformations + for (size_t i = 0; i < 500; ++i) + { + Mat aff = rngPartialAffMat(); + + const int method = GetParam(); + const int n = 100; + int m; + // LMEDS can't handle more than 50% outliers (by design) + if (method == LMEDS) + m = 3*n/5; + else + m = 2*n/5; + const float shift_outl = 15.f; + const float noise_level = 20.f; + + Mat fpts(1, n, CV_32FC2); + Mat tpts(1, n, CV_32FC2); + + randu(fpts, 0., 100.); + transform(fpts, tpts, aff); + + /* adding noise to some points */ + Mat outliers = tpts.colRange(m, n); + outliers.reshape(1) += shift_outl; + + Mat noise (outliers.size(), outliers.type()); + randu(noise, 0., noise_level); + outliers += noise; + + vector inliers; + Mat aff_est = estimateAffinePartial2D(fpts, tpts, inliers, method); + + EXPECT_FALSE(aff_est.empty()); + + EXPECT_NEAR(0., cvtest::norm(aff_est, aff, NORM_INF), 1e-4); + + bool inliers_good = count(inliers.begin(), inliers.end(), 1) == m && + m == accumulate(inliers.begin(), inliers.begin() + m, 0); + + EXPECT_TRUE(inliers_good); + } +} + +INSTANTIATE_TEST_CASE_P(Calib3d, EstimateAffinePartial2D, Method::all()); diff --git a/modules/core/include/opencv2/core.hpp b/modules/core/include/opencv2/core.hpp index 4a21aab66e..97a13f04ec 100644 --- a/modules/core/include/opencv2/core.hpp +++ b/modules/core/include/opencv2/core.hpp @@ -1650,7 +1650,7 @@ m.cols or m.cols-1. @param dst output array of the same size and depth as src; it has as many channels as m.rows. @param m transformation 2x2 or 2x3 floating-point matrix. -@sa perspectiveTransform, getAffineTransform, estimateRigidTransform, warpAffine, warpPerspective +@sa perspectiveTransform, getAffineTransform, estimateAffine2D, warpAffine, warpPerspective */ CV_EXPORTS_W void transform(InputArray src, OutputArray dst, InputArray m ); diff --git a/modules/stitching/include/opencv2/stitching.hpp b/modules/stitching/include/opencv2/stitching.hpp index c78be1d665..387e1dec78 100644 --- a/modules/stitching/include/opencv2/stitching.hpp +++ b/modules/stitching/include/opencv2/stitching.hpp @@ -69,7 +69,29 @@ one can combine and use them separately. The implemented stitching pipeline is very similar to the one proposed in @cite BL07 . -![image](StitchingPipeline.jpg) +![stitching pipeline](StitchingPipeline.jpg) + +Camera models +------------- + +There are currently 2 camera models implemented in stitching pipeline. + +- _Homography model_ expecting perspective transformations between images + implemented in @ref cv::detail::BestOf2NearestMatcher cv::detail::HomographyBasedEstimator + cv::detail::BundleAdjusterReproj cv::detail::BundleAdjusterRay +- _Affine model_ expecting affine transformation with 6 DOF or 4 DOF implemented in + @ref cv::detail::AffineBestOf2NearestMatcher cv::detail::AffineBasedEstimator + cv::detail::BundleAdjusterAffine cv::detail::BundleAdjusterAffinePartial cv::AffineWarper + +Homography model is useful for creating photo panoramas captured by camera, +while affine-based model can be used to stitch scans and object captured by +specialized devices. Use @ref cv::Stitcher::create to get preconfigured pipeline for one +of those models. + +@note +Certain detailed settings of @ref cv::Stitcher might not make sense. Especially +you should not mix classes implementing affine model and classes implementing +Homography model, as they work with different transformations. @{ @defgroup stitching_match Features Finding and Images Matching @@ -110,6 +132,22 @@ public: ERR_HOMOGRAPHY_EST_FAIL = 2, ERR_CAMERA_PARAMS_ADJUST_FAIL = 3 }; + enum Mode + { + /** Mode for creating photo panoramas. Expects images under perspective + transformation and projects resulting pano to sphere. + + @sa detail::BestOf2NearestMatcher SphericalWarper + */ + PANORAMA = 0, + /** Mode for composing scans. Expects images under affine transformation does + not compensate exposure by default. + + @sa detail::AffineBestOf2NearestMatcher AffineWarper + */ + SCANS = 1, + + }; // Stitcher() {} /** @brief Creates a stitcher with the default parameters. @@ -118,6 +156,15 @@ public: @return Stitcher class instance. */ static Stitcher createDefault(bool try_use_gpu = false); + /** @brief Creates a Stitcher configured in one of the stitching modes. + + @param mode Scenario for stitcher operation. This is usually determined by source of images + to stitch and their transformation. Default parameters will be chosen for operation in given + scenario. + @param try_use_gpu Flag indicating whether GPU should be used whenever it's possible. + @return Stitcher class instance. + */ + static Ptr create(Mode mode = PANORAMA, bool try_use_gpu = false); CV_WRAP double registrationResol() const { return registr_resol_; } CV_WRAP void setRegistrationResol(double resol_mpx) { registr_resol_ = resol_mpx; } @@ -159,6 +206,13 @@ public: void setBundleAdjuster(Ptr bundle_adjuster) { bundle_adjuster_ = bundle_adjuster; } + /* TODO OpenCV ABI 4.x + Ptr estimator() { return estimator_; } + const Ptr estimator() const { return estimator_; } + void setEstimator(Ptr estimator) + { estimator_ = estimator; } + */ + Ptr warper() { return warper_; } const Ptr warper() const { return warper_; } void setWarper(Ptr creator) { warper_ = creator; } @@ -233,6 +287,9 @@ private: Ptr features_matcher_; cv::UMat matching_mask_; Ptr bundle_adjuster_; + /* TODO OpenCV ABI 4.x + Ptr estimator_; + */ bool do_wave_correct_; detail::WaveCorrectKind wave_correct_kind_; Ptr warper_; diff --git a/modules/stitching/include/opencv2/stitching/detail/matchers.hpp b/modules/stitching/include/opencv2/stitching/detail/matchers.hpp index a46efc5760..bc81a84629 100644 --- a/modules/stitching/include/opencv2/stitching/detail/matchers.hpp +++ b/modules/stitching/include/opencv2/stitching/detail/matchers.hpp @@ -83,9 +83,27 @@ public: @sa detail::ImageFeatures, Rect_ */ void operator ()(InputArray image, ImageFeatures &features, const std::vector &rois); + /** @brief Finds features in the given images in parallel. + + @param images Source images + @param features Found features for each image + @param rois Regions of interest for each image + + @sa detail::ImageFeatures, Rect_ + */ + void operator ()(InputArrayOfArrays images, std::vector &features, + const std::vector > &rois); + /** @overload */ + void operator ()(InputArrayOfArrays images, std::vector &features); /** @brief Frees unused memory allocated before if there is any. */ virtual void collectGarbage() {} + /* TODO OpenCV ABI 4.x + reimplement this as public method similar to FeaturesMatcher and remove private function hack + @return True, if it's possible to use the same finder instance in parallel, false otherwise + bool isThreadSafe() const { return is_thread_safe_; } + */ + protected: /** @brief This method must implement features finding logic in order to make the wrappers detail::FeaturesFinder::operator()_ work. @@ -95,6 +113,10 @@ protected: @sa detail::ImageFeatures */ virtual void find(InputArray image, ImageFeatures &features) = 0; + /** @brief uses dynamic_cast to determine thread-safety + @return True, if it's possible to use the same finder instance in parallel, false otherwise + */ + bool isThreadSafe() const; }; /** @brief SURF features finder. @@ -152,7 +174,6 @@ private: Ptr akaze; }; - #ifdef HAVE_OPENCV_XFEATURES2D class CV_EXPORTS SurfFeaturesFinderGpu : public FeaturesFinder { @@ -177,7 +198,10 @@ private: /** @brief Structure containing information about matches between two images. -It's assumed that there is a homography between those images. +It's assumed that there is a transformation between those images. Transformation may be +homography or affine transformation based on selected matcher. + +@sa detail::FeaturesMatcher */ struct CV_EXPORTS MatchesInfo { @@ -189,7 +213,7 @@ struct CV_EXPORTS MatchesInfo std::vector matches; std::vector inliers_mask; //!< Geometrically consistent matches mask int num_inliers; //!< Number of geometrically consistent matches - Mat H; //!< Estimated homography + Mat H; //!< Estimated transformation double confidence; //!< Confidence two images are from the same panorama }; @@ -288,6 +312,41 @@ protected: int range_width_; }; +/** @brief Features matcher similar to cv::detail::BestOf2NearestMatcher which +finds two best matches for each feature and leaves the best one only if the +ratio between descriptor distances is greater than the threshold match_conf. + +Unlike cv::detail::BestOf2NearestMatcher this matcher uses affine +transformation (affine trasformation estimate will be placed in matches_info). + +@sa cv::detail::FeaturesMatcher cv::detail::BestOf2NearestMatcher + */ +class CV_EXPORTS AffineBestOf2NearestMatcher : public BestOf2NearestMatcher +{ +public: + /** @brief Constructs a "best of 2 nearest" matcher that expects affine trasformation + between images + + @param full_affine whether to use full affine transformation with 6 degress of freedom or reduced + transformation with 4 degrees of freedom using only rotation, translation and uniform scaling + @param try_use_gpu Should try to use GPU or not + @param match_conf Match distances ration threshold + @param num_matches_thresh1 Minimum number of matches required for the 2D affine transform + estimation used in the inliers classification step + + @sa cv::estimateAffine2D cv::estimateAffinePartial2D + */ + AffineBestOf2NearestMatcher(bool full_affine = false, bool try_use_gpu = false, + float match_conf = 0.3f, int num_matches_thresh1 = 6) : + BestOf2NearestMatcher(try_use_gpu, match_conf, num_matches_thresh1, num_matches_thresh1), + full_affine_(full_affine) {} + +protected: + void match(const ImageFeatures &features1, const ImageFeatures &features2, MatchesInfo &matches_info); + + bool full_affine_; +}; + //! @} stitching_match } // namespace detail diff --git a/modules/stitching/include/opencv2/stitching/detail/motion_estimators.hpp b/modules/stitching/include/opencv2/stitching/detail/motion_estimators.hpp index 9f7b8bde1a..5276fd1d98 100644 --- a/modules/stitching/include/opencv2/stitching/detail/motion_estimators.hpp +++ b/modules/stitching/include/opencv2/stitching/detail/motion_estimators.hpp @@ -109,6 +109,21 @@ private: bool is_focals_estimated_; }; +/** @brief Affine transformation based estimator. + +This estimator uses pairwise tranformations estimated by matcher to estimate +final transformation for each camera. + +@sa cv::detail::HomographyBasedEstimator + */ +class CV_EXPORTS AffineBasedEstimator : public Estimator +{ +private: + virtual bool estimate(const std::vector &features, + const std::vector &pairwise_matches, + std::vector &cameras); +}; + /** @brief Base class for all camera parameters refinement methods. */ class CV_EXPORTS BundleAdjusterBase : public Estimator @@ -195,6 +210,26 @@ protected: }; +/** @brief Stub bundle adjuster that does nothing. + */ +class CV_EXPORTS NoBundleAdjuster : public BundleAdjusterBase +{ +public: + NoBundleAdjuster() : BundleAdjusterBase(0, 0) {} + +private: + bool estimate(const std::vector &, const std::vector &, + std::vector &) + { + return true; + } + void setUpInitialCameraParams(const std::vector &) {} + void obtainRefinedCameraParams(std::vector &) const {} + void calcError(Mat &) {} + void calcJacobian(Mat &) {} +}; + + /** @brief Implementation of the camera parameters refinement algorithm which minimizes sum of the reprojection error squares @@ -236,6 +271,54 @@ private: }; +/** @brief Bundle adjuster that expects affine transformation +represented in homogeneous coordinates in R for each camera param. Implements +camera parameters refinement algorithm which minimizes sum of the reprojection +error squares + +It estimates all transformation parameters. Refinement mask is ignored. + +@sa AffineBasedEstimator AffineBestOf2NearestMatcher BundleAdjusterAffinePartial + */ +class CV_EXPORTS BundleAdjusterAffine : public BundleAdjusterBase +{ +public: + BundleAdjusterAffine() : BundleAdjusterBase(6, 2) {} + +private: + void setUpInitialCameraParams(const std::vector &cameras); + void obtainRefinedCameraParams(std::vector &cameras) const; + void calcError(Mat &err); + void calcJacobian(Mat &jac); + + Mat err1_, err2_; +}; + + +/** @brief Bundle adjuster that expects affine transformation with 4 DOF +represented in homogeneous coordinates in R for each camera param. Implements +camera parameters refinement algorithm which minimizes sum of the reprojection +error squares + +It estimates all transformation parameters. Refinement mask is ignored. + +@sa AffineBasedEstimator AffineBestOf2NearestMatcher BundleAdjusterAffine + */ +class CV_EXPORTS BundleAdjusterAffinePartial : public BundleAdjusterBase +{ +public: + BundleAdjusterAffinePartial() : BundleAdjusterBase(4, 2) {} + +private: + void setUpInitialCameraParams(const std::vector &cameras); + void obtainRefinedCameraParams(std::vector &cameras) const; + void calcError(Mat &err); + void calcJacobian(Mat &jac); + + Mat err1_, err2_; +}; + + enum WaveCorrectKind { WAVE_CORRECT_HORIZ, diff --git a/modules/stitching/include/opencv2/stitching/detail/warpers.hpp b/modules/stitching/include/opencv2/stitching/detail/warpers.hpp index f848c49360..1515d76260 100644 --- a/modules/stitching/include/opencv2/stitching/detail/warpers.hpp +++ b/modules/stitching/include/opencv2/stitching/detail/warpers.hpp @@ -205,6 +205,34 @@ protected: }; +/** @brief Affine warper that uses rotations and translations + + Uses affine transformation in homogeneous coordinates to represent both rotation and + translation in camera rotation matrix. + */ +class CV_EXPORTS AffineWarper : public PlaneWarper +{ +public: + /** @brief Construct an instance of the affine warper class. + + @param scale Projected image scale multiplier + */ + AffineWarper(float scale = 1.f) : PlaneWarper(scale) {} + + Point2f warpPoint(const Point2f &pt, InputArray K, InputArray R); + Rect buildMaps(Size src_size, InputArray K, InputArray R, OutputArray xmap, OutputArray ymap); + Point warp(InputArray src, InputArray K, InputArray R, + int interp_mode, int border_mode, OutputArray dst); + Rect warpRoi(Size src_size, InputArray K, InputArray R); + +protected: + /** @brief Extracts rotation and translation matrices from matrix H representing + affine transformation in homogeneous coordinates + */ + void getRTfromHomogeneous(InputArray H, Mat &R, Mat &T); +}; + + struct CV_EXPORTS SphericalProjector : ProjectorBase { void mapForward(float x, float y, float &u, float &v); diff --git a/modules/stitching/include/opencv2/stitching/warpers.hpp b/modules/stitching/include/opencv2/stitching/warpers.hpp index 23f69516ac..139e0522b9 100644 --- a/modules/stitching/include/opencv2/stitching/warpers.hpp +++ b/modules/stitching/include/opencv2/stitching/warpers.hpp @@ -68,6 +68,15 @@ public: Ptr create(float scale) const { return makePtr(scale); } }; +/** @brief Affine warper factory class. + @sa detail::AffineWarper + */ +class AffineWarper : public WarperCreator +{ +public: + Ptr create(float scale) const { return makePtr(scale); } +}; + /** @brief Cylindrical warper factory class. @sa detail::CylindricalWarper */ diff --git a/modules/stitching/perf/opencl/perf_stitch.cpp b/modules/stitching/perf/opencl/perf_stitch.cpp index ce7c3a9f11..8b25c50940 100644 --- a/modules/stitching/perf/opencl/perf_stitch.cpp +++ b/modules/stitching/perf/opencl/perf_stitch.cpp @@ -7,9 +7,13 @@ #include "../perf_precomp.hpp" #include "opencv2/ts/ocl_perf.hpp" +#ifdef HAVE_OPENCL + +namespace cvtest { +namespace ocl { + using namespace cv; using namespace perf; -using namespace cvtest::ocl; using namespace std; using namespace std::tr1; @@ -19,7 +23,7 @@ using namespace std::tr1; typedef TestBaseWithParam stitch; -#ifdef HAVE_OPENCV_NONFREE_TODO_FIND_WHY_SURF_IS_NOT_ABLE_TO_STITCH_PANOS +#ifdef HAVE_OPENCV_XFEATURES2D #define TEST_DETECTORS testing::Values("surf", "orb") #else #define TEST_DETECTORS testing::Values("orb") @@ -142,3 +146,7 @@ OCL_PERF_TEST_P(stitch, boat, TEST_DETECTORS) SANITY_CHECK_NOTHING(); } + +} } // namespace cvtest::ocl + +#endif // HAVE_OPENCL diff --git a/modules/stitching/perf/opencl/perf_warpers.cpp b/modules/stitching/perf/opencl/perf_warpers.cpp index 57ca9a602d..1aa738cd9f 100644 --- a/modules/stitching/perf/opencl/perf_warpers.cpp +++ b/modules/stitching/perf/opencl/perf_warpers.cpp @@ -54,7 +54,8 @@ enum { SphericalWarperType = 0, CylindricalWarperType = 1, - PlaneWarperType = 2 + PlaneWarperType = 2, + AffineWarperType = 3, }; class WarperBase @@ -69,6 +70,8 @@ public: creator = makePtr(); else if (type == PlaneWarperType) creator = makePtr(); + else if (type == AffineWarperType) + creator = makePtr(); CV_Assert(!creator.empty()); K = Mat::eye(3, 3, CV_32FC1); @@ -98,7 +101,7 @@ private: Mat K, R; }; -CV_ENUM(WarperType, SphericalWarperType, CylindricalWarperType, PlaneWarperType) +CV_ENUM(WarperType, SphericalWarperType, CylindricalWarperType, PlaneWarperType, AffineWarperType) typedef tuple StitchingWarpersParams; typedef TestBaseWithParam StitchingWarpersFixture; diff --git a/modules/stitching/perf/perf_estimators.cpp b/modules/stitching/perf/perf_estimators.cpp new file mode 100644 index 0000000000..7de470c0e5 --- /dev/null +++ b/modules/stitching/perf/perf_estimators.cpp @@ -0,0 +1,100 @@ +#include "perf_precomp.hpp" +#include "opencv2/imgcodecs.hpp" +#include "opencv2/opencv_modules.hpp" + +using namespace std; +using namespace cv; +using namespace perf; +using std::tr1::tuple; +using std::tr1::get; + +typedef TestBaseWithParam > bundleAdjuster; + +#ifdef HAVE_OPENCV_XFEATURES2D +#define TEST_DETECTORS testing::Values("surf", "orb") +#else +#define TEST_DETECTORS testing::Values("orb") +#endif +#define WORK_MEGAPIX 0.6 +#define AFFINE_FUNCTIONS testing::Values("affinePartial", "affine") + +PERF_TEST_P(bundleAdjuster, affine, testing::Combine(TEST_DETECTORS, AFFINE_FUNCTIONS)) +{ + Mat img1, img1_full = imread(getDataPath("stitching/s1.jpg")); + Mat img2, img2_full = imread(getDataPath("stitching/s2.jpg")); + float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); + float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); + resize(img1_full, img1, Size(), scale1, scale1); + resize(img2_full, img2, Size(), scale2, scale2); + + string detector = get<0>(GetParam()); + string affine_fun = get<1>(GetParam()); + + Ptr finder; + Ptr matcher; + Ptr bundle_adjuster; + if (detector == "surf") + finder = makePtr(); + else if (detector == "orb") + finder = makePtr(); + if (affine_fun == "affinePartial") + { + matcher = makePtr(false); + bundle_adjuster = makePtr(); + } + else if (affine_fun == "affine") + { + matcher = makePtr(true); + bundle_adjuster = makePtr(); + } + Ptr estimator = makePtr(); + + std::vector images; + images.push_back(img1), images.push_back(img2); + std::vector features; + std::vector pairwise_matches; + std::vector cameras; + std::vector cameras2; + + (*finder)(images, features); + (*matcher)(features, pairwise_matches); + if (!(*estimator)(features, pairwise_matches, cameras)) + FAIL() << "estimation failed. this should never happen."; + // this is currently required + for (size_t i = 0; i < cameras.size(); ++i) + { + Mat R; + cameras[i].R.convertTo(R, CV_32F); + cameras[i].R = R; + } + + cameras2 = cameras; + bool success = true; + while(next()) + { + cameras = cameras2; // revert cameras back to original initial guess + startTimer(); + success = (*bundle_adjuster)(features, pairwise_matches, cameras); + stopTimer(); + } + + EXPECT_TRUE(success); + EXPECT_TRUE(cameras.size() == 2); + + // fist camera should be just identity + Mat &first = cameras[0].R; + SANITY_CHECK(first, 1e-3, ERROR_ABSOLUTE); + // second camera should be the estimated transform between images + // separate rotation and translation in transform matrix + Mat T_second (cameras[1].R, Range(0, 2), Range(2, 3)); + Mat R_second (cameras[1].R, Range(0, 2), Range(0, 2)); + Mat h (cameras[1].R, Range(2, 3), Range::all()); + SANITY_CHECK(T_second, 5, ERROR_ABSOLUTE); // allow 5 pixels diff in translations + SANITY_CHECK(R_second, .01, ERROR_ABSOLUTE); // rotations must be more precise + // last row should be precisely (0, 0, 1) as it is just added for representation in homogeneous + // coordinates + EXPECT_TRUE(h.type() == CV_32F); + EXPECT_FLOAT_EQ(h.at(0), 0.f); + EXPECT_FLOAT_EQ(h.at(1), 0.f); + EXPECT_FLOAT_EQ(h.at(2), 1.f); +} diff --git a/modules/stitching/perf/perf_matchers.cpp b/modules/stitching/perf/perf_matchers.cpp new file mode 100644 index 0000000000..4e5b03d430 --- /dev/null +++ b/modules/stitching/perf/perf_matchers.cpp @@ -0,0 +1,301 @@ +#include "perf_precomp.hpp" +#include "opencv2/imgcodecs.hpp" +#include "opencv2/opencv_modules.hpp" +#include "opencv2/flann.hpp" + +using namespace std; +using namespace cv; +using namespace perf; +using std::tr1::make_tuple; +using std::tr1::get; + +typedef TestBaseWithParam FeaturesFinderVec; +typedef TestBaseWithParam match; +typedef std::tr1::tuple matchVector_t; +typedef TestBaseWithParam matchVector; + +#define NUMBER_IMAGES testing::Values(1, 5, 20) +#define SURF_MATCH_CONFIDENCE 0.65f +#define ORB_MATCH_CONFIDENCE 0.3f +#define WORK_MEGAPIX 0.6 + +#ifdef HAVE_OPENCV_XFEATURES2D +#define TEST_DETECTORS testing::Values("surf", "orb") +#else +#define TEST_DETECTORS testing::Values("orb") +#endif + +PERF_TEST_P(FeaturesFinderVec, ParallelFeaturesFinder, NUMBER_IMAGES) +{ + Mat img = imread( getDataPath("stitching/a1.png") ); + vector imgs(GetParam(), img); + vector features(imgs.size()); + + Ptr featuresFinder = makePtr(); + + TEST_CYCLE() + { + (*featuresFinder)(imgs, features); + } + + SANITY_CHECK_NOTHING(); +} + +PERF_TEST_P(FeaturesFinderVec, SerialFeaturesFinder, NUMBER_IMAGES) +{ + Mat img = imread( getDataPath("stitching/a1.png") ); + vector imgs(GetParam(), img); + vector features(imgs.size()); + + Ptr featuresFinder = makePtr(); + + TEST_CYCLE() + { + for (size_t i = 0; i < imgs.size(); ++i) + (*featuresFinder)(imgs[i], features[i]); + } + + SANITY_CHECK_NOTHING(); +} + +PERF_TEST_P( match, bestOf2Nearest, TEST_DETECTORS) +{ + Mat img1, img1_full = imread( getDataPath("stitching/boat1.jpg") ); + Mat img2, img2_full = imread( getDataPath("stitching/boat2.jpg") ); + float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); + float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); + resize(img1_full, img1, Size(), scale1, scale1); + resize(img2_full, img2, Size(), scale2, scale2); + + Ptr finder; + Ptr matcher; + if (GetParam() == "surf") + { + finder = makePtr(); + matcher = makePtr(false, SURF_MATCH_CONFIDENCE); + } + else if (GetParam() == "orb") + { + finder = makePtr(); + matcher = makePtr(false, ORB_MATCH_CONFIDENCE); + } + else + { + FAIL() << "Unknown 2D features type: " << GetParam(); + } + + detail::ImageFeatures features1, features2; + (*finder)(img1, features1); + (*finder)(img2, features2); + + detail::MatchesInfo pairwise_matches; + + declare.in(features1.descriptors, features2.descriptors); + + while(next()) + { + cvflann::seed_random(42);//for predictive FlannBasedMatcher + startTimer(); + (*matcher)(features1, features2, pairwise_matches); + stopTimer(); + matcher->collectGarbage(); + } + + Mat dist (pairwise_matches.H, Range::all(), Range(2, 3)); + Mat R (pairwise_matches.H, Range::all(), Range(0, 2)); + // separate transform matrix, use lower error on rotations + SANITY_CHECK(dist, 1., ERROR_ABSOLUTE); + SANITY_CHECK(R, .015, ERROR_ABSOLUTE); +} + +PERF_TEST_P( matchVector, bestOf2NearestVectorFeatures, testing::Combine( + TEST_DETECTORS, + testing::Values(2, 4, 8)) + ) +{ + Mat img1, img1_full = imread( getDataPath("stitching/boat1.jpg") ); + Mat img2, img2_full = imread( getDataPath("stitching/boat2.jpg") ); + float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); + float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); + resize(img1_full, img1, Size(), scale1, scale1); + resize(img2_full, img2, Size(), scale2, scale2); + + Ptr finder; + Ptr matcher; + string detectorName = get<0>(GetParam()); + int featuresVectorSize = get<1>(GetParam()); + if (detectorName == "surf") + { + finder = makePtr(); + matcher = makePtr(false, SURF_MATCH_CONFIDENCE); + } + else if (detectorName == "orb") + { + finder = makePtr(); + matcher = makePtr(false, ORB_MATCH_CONFIDENCE); + } + else + { + FAIL() << "Unknown 2D features type: " << get<0>(GetParam()); + } + + detail::ImageFeatures features1, features2; + (*finder)(img1, features1); + (*finder)(img2, features2); + vector features; + vector pairwise_matches; + for(int i = 0; i < featuresVectorSize/2; i++) + { + features.push_back(features1); + features.push_back(features2); + } + + declare.time(200); + while(next()) + { + cvflann::seed_random(42);//for predictive FlannBasedMatcher + startTimer(); + (*matcher)(features, pairwise_matches); + stopTimer(); + matcher->collectGarbage(); + } + + size_t matches_count = 0; + for (size_t i = 0; i < pairwise_matches.size(); ++i) + { + if (pairwise_matches[i].src_img_idx < 0) + continue; + + EXPECT_TRUE(pairwise_matches[i].matches.size() > 100); + EXPECT_FALSE(pairwise_matches[i].H.empty()); + ++matches_count; + } + + EXPECT_TRUE(matches_count > 0); + + SANITY_CHECK_NOTHING(); +} + +PERF_TEST_P( match, affineBestOf2Nearest, TEST_DETECTORS) +{ + Mat img1, img1_full = imread( getDataPath("stitching/s1.jpg") ); + Mat img2, img2_full = imread( getDataPath("stitching/s2.jpg") ); + float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); + float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); + resize(img1_full, img1, Size(), scale1, scale1); + resize(img2_full, img2, Size(), scale2, scale2); + + Ptr finder; + Ptr matcher; + if (GetParam() == "surf") + { + finder = makePtr(); + matcher = makePtr(false, false, SURF_MATCH_CONFIDENCE); + } + else if (GetParam() == "orb") + { + finder = makePtr(); + matcher = makePtr(false, false, ORB_MATCH_CONFIDENCE); + } + else + { + FAIL() << "Unknown 2D features type: " << GetParam(); + } + + detail::ImageFeatures features1, features2; + (*finder)(img1, features1); + (*finder)(img2, features2); + + detail::MatchesInfo pairwise_matches; + + declare.in(features1.descriptors, features2.descriptors); + + while(next()) + { + cvflann::seed_random(42);//for predictive FlannBasedMatcher + startTimer(); + (*matcher)(features1, features2, pairwise_matches); + stopTimer(); + matcher->collectGarbage(); + } + + // separate rotation and translation in transform matrix + Mat T (pairwise_matches.H, Range(0, 2), Range(2, 3)); + Mat R (pairwise_matches.H, Range(0, 2), Range(0, 2)); + Mat h (pairwise_matches.H, Range(2, 3), Range::all()); + SANITY_CHECK(T, 5, ERROR_ABSOLUTE); // allow 5 pixels diff in translations + SANITY_CHECK(R, .01, ERROR_ABSOLUTE); // rotations must be more precise + // last row should be precisely (0, 0, 1) as it is just added for representation in homogeneous + // coordinates + EXPECT_DOUBLE_EQ(h.at(0), 0.); + EXPECT_DOUBLE_EQ(h.at(1), 0.); + EXPECT_DOUBLE_EQ(h.at(2), 1.); +} + +PERF_TEST_P( matchVector, affineBestOf2NearestVectorFeatures, testing::Combine( + TEST_DETECTORS, + testing::Values(2, 4, 8)) + ) +{ + Mat img1, img1_full = imread( getDataPath("stitching/s1.jpg") ); + Mat img2, img2_full = imread( getDataPath("stitching/s2.jpg") ); + float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); + float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); + resize(img1_full, img1, Size(), scale1, scale1); + resize(img2_full, img2, Size(), scale2, scale2); + + Ptr finder; + Ptr matcher; + string detectorName = get<0>(GetParam()); + int featuresVectorSize = get<1>(GetParam()); + if (detectorName == "surf") + { + finder = makePtr(); + matcher = makePtr(false, false, SURF_MATCH_CONFIDENCE); + } + else if (detectorName == "orb") + { + finder = makePtr(); + matcher = makePtr(false, false, ORB_MATCH_CONFIDENCE); + } + else + { + FAIL() << "Unknown 2D features type: " << get<0>(GetParam()); + } + + detail::ImageFeatures features1, features2; + (*finder)(img1, features1); + (*finder)(img2, features2); + vector features; + vector pairwise_matches; + for(int i = 0; i < featuresVectorSize/2; i++) + { + features.push_back(features1); + features.push_back(features2); + } + + declare.time(200); + while(next()) + { + cvflann::seed_random(42);//for predictive FlannBasedMatcher + startTimer(); + (*matcher)(features, pairwise_matches); + stopTimer(); + matcher->collectGarbage(); + } + + size_t matches_count = 0; + for (size_t i = 0; i < pairwise_matches.size(); ++i) + { + if (pairwise_matches[i].src_img_idx < 0) + continue; + + EXPECT_TRUE(pairwise_matches[i].matches.size() > 400); + EXPECT_FALSE(pairwise_matches[i].H.empty()); + ++matches_count; + } + + EXPECT_TRUE(matches_count > 0); + + SANITY_CHECK_NOTHING(); +} diff --git a/modules/stitching/perf/perf_stich.cpp b/modules/stitching/perf/perf_stich.cpp index 74fd1ccaf9..5a6d0237e9 100644 --- a/modules/stitching/perf/perf_stich.cpp +++ b/modules/stitching/perf/perf_stich.cpp @@ -1,12 +1,11 @@ #include "perf_precomp.hpp" #include "opencv2/imgcodecs.hpp" -#include "opencv2/flann.hpp" #include "opencv2/opencv_modules.hpp" using namespace std; using namespace cv; using namespace perf; -using std::tr1::make_tuple; +using std::tr1::tuple; using std::tr1::get; #define SURF_MATCH_CONFIDENCE 0.65f @@ -14,15 +13,14 @@ using std::tr1::get; #define WORK_MEGAPIX 0.6 typedef TestBaseWithParam stitch; -typedef TestBaseWithParam match; -typedef std::tr1::tuple matchVector_t; -typedef TestBaseWithParam matchVector; +typedef TestBaseWithParam > stitchDatasets; -#ifdef HAVE_OPENCV_XFEATURES2D_TODO_FIND_WHY_SURF_IS_NOT_ABLE_TO_STITCH_PANOS +#ifdef HAVE_OPENCV_XFEATURES2D #define TEST_DETECTORS testing::Values("surf", "orb") #else #define TEST_DETECTORS testing::Values("orb") #endif +#define AFFINE_DATASETS testing::Values("s", "budapest", "newspaper", "prague") PERF_TEST_P(stitch, a123, TEST_DETECTORS) { @@ -93,118 +91,101 @@ PERF_TEST_P(stitch, b12, TEST_DETECTORS) stopTimer(); } + EXPECT_NEAR(pano.size().width, 1117, 50); + EXPECT_NEAR(pano.size().height, 642, 30); + Mat pano_small; - if (!pano.empty()) - resize(pano, pano_small, Size(320, 240), 0, 0, INTER_AREA); + resize(pano, pano_small, Size(320, 240), 0, 0, INTER_AREA); - SANITY_CHECK(pano_small, 5); + // results from orb and surf are slightly different (but both of them are good). Regression data + // are for orb. + /* transformations are: + orb: + [0.99213386, 0.062509097, -351.83731; + -0.073042989, 1.0615162, -89.869858; + 0.0005330033, -4.0937066e-05, 1] + surf: + [1.0034728, 0.022535477, -352.76849; + -0.080653802, 1.0742083, -89.602058; + 0.0004876224, 0.00012311155, 1] + */ + if (GetParam() == "orb") + SANITY_CHECK(pano_small, 5); + else + SANITY_CHECK_NOTHING(); } -PERF_TEST_P( match, bestOf2Nearest, TEST_DETECTORS) +PERF_TEST_P(stitchDatasets, affine, testing::Combine(AFFINE_DATASETS, TEST_DETECTORS)) { - Mat img1, img1_full = imread( getDataPath("stitching/b1.png") ); - Mat img2, img2_full = imread( getDataPath("stitching/b2.png") ); - float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); - float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); - resize(img1_full, img1, Size(), scale1, scale1); - resize(img2_full, img2, Size(), scale2, scale2); + string dataset = get<0>(GetParam()); + string detector = get<1>(GetParam()); - Ptr finder; - Ptr matcher; - if (GetParam() == "surf") - { - finder = makePtr(); - matcher = makePtr(false, SURF_MATCH_CONFIDENCE); - } - else if (GetParam() == "orb") - { - finder = makePtr(); - matcher = makePtr(false, ORB_MATCH_CONFIDENCE); - } + Mat pano; + vector imgs; + int width, height, allowed_diff = 10; + Ptr featuresFinder; + + if(detector == "orb") + featuresFinder = makePtr(); else + featuresFinder = makePtr(); + + if(dataset == "budapest") { - FAIL() << "Unknown 2D features type: " << GetParam(); + imgs.push_back(imread(getDataPath("stitching/budapest1.jpg"))); + imgs.push_back(imread(getDataPath("stitching/budapest2.jpg"))); + imgs.push_back(imread(getDataPath("stitching/budapest3.jpg"))); + imgs.push_back(imread(getDataPath("stitching/budapest4.jpg"))); + imgs.push_back(imread(getDataPath("stitching/budapest5.jpg"))); + imgs.push_back(imread(getDataPath("stitching/budapest6.jpg"))); + width = 2313; + height = 1158; + // this dataset is big, the results between surf and orb differ slightly, + // but both are still good + allowed_diff = 27; + } + else if (dataset == "newspaper") + { + imgs.push_back(imread(getDataPath("stitching/newspaper1.jpg"))); + imgs.push_back(imread(getDataPath("stitching/newspaper2.jpg"))); + imgs.push_back(imread(getDataPath("stitching/newspaper3.jpg"))); + imgs.push_back(imread(getDataPath("stitching/newspaper4.jpg"))); + width = 1791; + height = 1136; + // we need to boost ORB number of features to be able to stitch this dataset + // SURF works just fine with default settings + if(detector == "orb") + featuresFinder = makePtr(Size(3,1), 3000); + } + else if (dataset == "prague") + { + imgs.push_back(imread(getDataPath("stitching/prague1.jpg"))); + imgs.push_back(imread(getDataPath("stitching/prague2.jpg"))); + width = 983; + height = 1759; + } + else // dataset == "s" + { + imgs.push_back(imread(getDataPath("stitching/s1.jpg"))); + imgs.push_back(imread(getDataPath("stitching/s2.jpg"))); + width = 1815; + height = 700; } - detail::ImageFeatures features1, features2; - (*finder)(img1, features1); - (*finder)(img2, features2); - - detail::MatchesInfo pairwise_matches; - - declare.in(features1.descriptors, features2.descriptors); + declare.time(30 * 20).iterations(20); while(next()) { - cvflann::seed_random(42);//for predictive FlannBasedMatcher + Ptr stitcher = Stitcher::create(Stitcher::SCANS, false); + stitcher->setFeaturesFinder(featuresFinder); + startTimer(); - (*matcher)(features1, features2, pairwise_matches); + stitcher->stitch(imgs, pano); stopTimer(); - matcher->collectGarbage(); } - std::vector& matches = pairwise_matches.matches; - if (GetParam() == "orb") matches.resize(0); - for(size_t q = 0; q < matches.size(); ++q) - if (matches[q].imgIdx < 0) { matches.resize(q); break;} - SANITY_CHECK_MATCHES(matches); -} - -PERF_TEST_P( matchVector, bestOf2NearestVectorFeatures, testing::Combine( - TEST_DETECTORS, - testing::Values(2, 4, 8)) - ) -{ - Mat img1, img1_full = imread( getDataPath("stitching/b1.png") ); - Mat img2, img2_full = imread( getDataPath("stitching/b2.png") ); - float scale1 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img1_full.total())); - float scale2 = (float)std::min(1.0, sqrt(WORK_MEGAPIX * 1e6 / img2_full.total())); - resize(img1_full, img1, Size(), scale1, scale1); - resize(img2_full, img2, Size(), scale2, scale2); - - Ptr finder; - Ptr matcher; - string detectorName = get<0>(GetParam()); - int featuresVectorSize = get<1>(GetParam()); - if (detectorName == "surf") - { - finder = makePtr(); - matcher = makePtr(false, SURF_MATCH_CONFIDENCE); - } - else if (detectorName == "orb") - { - finder = makePtr(); - matcher = makePtr(false, ORB_MATCH_CONFIDENCE); - } - else - { - FAIL() << "Unknown 2D features type: " << get<0>(GetParam()); - } - - detail::ImageFeatures features1, features2; - (*finder)(img1, features1); - (*finder)(img2, features2); - vector features; - vector pairwise_matches; - for(int i = 0; i < featuresVectorSize/2; i++) - { - features.push_back(features1); - features.push_back(features2); - } - - declare.time(200); - while(next()) - { - cvflann::seed_random(42);//for predictive FlannBasedMatcher - startTimer(); - (*matcher)(features, pairwise_matches); - stopTimer(); - matcher->collectGarbage(); - } - - - std::vector& matches = pairwise_matches[detectorName == "surf" ? 1 : 0].matches; - for(size_t q = 0; q < matches.size(); ++q) - if (matches[q].imgIdx < 0) { matches.resize(q); break;} - SANITY_CHECK_MATCHES(matches); + EXPECT_NEAR(pano.size().width, width, allowed_diff); + EXPECT_NEAR(pano.size().height, height, allowed_diff); + + SANITY_CHECK_NOTHING(); } diff --git a/modules/stitching/src/matchers.cpp b/modules/stitching/src/matchers.cpp index df424a797a..edd6b61931 100644 --- a/modules/stitching/src/matchers.cpp +++ b/modules/stitching/src/matchers.cpp @@ -107,6 +107,35 @@ private: }; +struct FindFeaturesBody : ParallelLoopBody +{ + FindFeaturesBody(FeaturesFinder &finder, InputArrayOfArrays images, + std::vector &features, const std::vector > *rois) + : finder_(finder), images_(images), features_(features), rois_(rois) {} + + void operator ()(const Range &r) const + { + for (int i = r.start; i < r.end; ++i) + { + Mat image = images_.getMat(i); + if (rois_) + finder_(image, features_[i], (*rois_)[i]); + else + finder_(image, features_[i]); + } + } + +private: + FeaturesFinder &finder_; + InputArrayOfArrays images_; + std::vector &features_; + const std::vector > *rois_; + + // to cease visual studio warning + void operator =(const FindFeaturesBody&); +}; + + ////////////////////////////////////////////////////////////////////////////// typedef std::set > MatchesSet; @@ -324,6 +353,55 @@ void FeaturesFinder::operator ()(InputArray image, ImageFeatures &features, cons } +void FeaturesFinder::operator ()(InputArrayOfArrays images, std::vector &features) +{ + size_t count = images.total(); + features.resize(count); + + FindFeaturesBody body(*this, images, features, NULL); + if (isThreadSafe()) + parallel_for_(Range(0, static_cast(count)), body); + else + body(Range(0, static_cast(count))); +} + + +void FeaturesFinder::operator ()(InputArrayOfArrays images, std::vector &features, + const std::vector > &rois) +{ + CV_Assert(rois.size() == images.total()); + size_t count = images.total(); + features.resize(count); + + FindFeaturesBody body(*this, images, features, &rois); + if (isThreadSafe()) + parallel_for_(Range(0, static_cast(count)), body); + else + body(Range(0, static_cast(count))); +} + + +bool FeaturesFinder::isThreadSafe() const +{ + if (ocl::useOpenCL()) + { + return false; + } + if (dynamic_cast(this)) + { + return true; + } + else if (dynamic_cast(this)) + { + return true; + } + else + { + return false; + } +} + + SurfFeaturesFinder::SurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers, int num_octaves_descr, int num_layers_descr) { @@ -722,5 +800,57 @@ void BestOf2NearestRangeMatcher::operator ()(const std::vector &f } +void AffineBestOf2NearestMatcher::match(const ImageFeatures &features1, const ImageFeatures &features2, + MatchesInfo &matches_info) +{ + (*impl_)(features1, features2, matches_info); + + // Check if it makes sense to find transform + if (matches_info.matches.size() < static_cast(num_matches_thresh1_)) + return; + + // Construct point-point correspondences for transform estimation + Mat src_points(1, static_cast(matches_info.matches.size()), CV_32FC2); + Mat dst_points(1, static_cast(matches_info.matches.size()), CV_32FC2); + for (size_t i = 0; i < matches_info.matches.size(); ++i) + { + const cv::DMatch &m = matches_info.matches[i]; + src_points.at(0, static_cast(i)) = features1.keypoints[m.queryIdx].pt; + dst_points.at(0, static_cast(i)) = features2.keypoints[m.trainIdx].pt; + } + + // Find pair-wise motion + if (full_affine_) + matches_info.H = estimateAffine2D(src_points, dst_points, matches_info.inliers_mask); + else + matches_info.H = estimateAffinePartial2D(src_points, dst_points, matches_info.inliers_mask); + + if (matches_info.H.empty()) { + // could not find transformation + matches_info.confidence = 0; + matches_info.num_inliers = 0; + return; + } + + // Find number of inliers + matches_info.num_inliers = 0; + for (size_t i = 0; i < matches_info.inliers_mask.size(); ++i) + if (matches_info.inliers_mask[i]) + matches_info.num_inliers++; + + // These coeffs are from paper M. Brown and D. Lowe. "Automatic Panoramic + // Image Stitching using Invariant Features" + matches_info.confidence = + matches_info.num_inliers / (8 + 0.3 * matches_info.matches.size()); + + /* should we remove matches between too close images? */ + // matches_info.confidence = matches_info.confidence > 3. ? 0. : matches_info.confidence; + + // extend H to represent linear tranformation in homogeneous coordinates + matches_info.H.push_back(Mat::zeros(1, 3, CV_64F)); + matches_info.H.at(2, 2) = 1; +} + + } // namespace detail } // namespace cv diff --git a/modules/stitching/src/motion_estimators.cpp b/modules/stitching/src/motion_estimators.cpp index dd87725bb6..f76309f777 100644 --- a/modules/stitching/src/motion_estimators.cpp +++ b/modules/stitching/src/motion_estimators.cpp @@ -88,6 +88,28 @@ struct CalcRotation }; +/** + * @brief Functor calculating final tranformation by chaining linear transformations + */ +struct CalcAffineTransform +{ + CalcAffineTransform(int _num_images, + const std::vector &_pairwise_matches, + std::vector &_cameras) + : num_images(_num_images), pairwise_matches(&_pairwise_matches[0]), cameras(&_cameras[0]) {} + + void operator()(const GraphEdge &edge) + { + int pair_idx = edge.from * num_images + edge.to; + cameras[edge.to].R = cameras[edge.from].R * pairwise_matches[pair_idx].H; + } + + int num_images; + const MatchesInfo *pairwise_matches; + CameraParams *cameras; +}; + + ////////////////////////////////////////////////////////////////////////////// void calcDeriv(const Mat &err1, const Mat &err2, double h, Mat res) @@ -171,6 +193,31 @@ bool HomographyBasedEstimator::estimate( } +////////////////////////////////////////////////////////////////////////////// + +bool AffineBasedEstimator::estimate(const std::vector &features, + const std::vector &pairwise_matches, + std::vector &cameras) +{ + cameras.resize(features.size()); + const int num_images = static_cast(features.size()); + + // find maximum spaning tree on pairwise matches + cv::detail::Graph span_tree; + std::vector span_tree_centers; + // uses number of inliers as weights + findMaxSpanningTree(num_images, pairwise_matches, span_tree, + span_tree_centers); + + // compute final transform by chaining H together + span_tree.walkBreadthFirst( + span_tree_centers[0], + CalcAffineTransform(num_images, pairwise_matches, cameras)); + // this estimator never fails + return true; +} + + ////////////////////////////////////////////////////////////////////////////// bool BundleAdjusterBase::estimate(const std::vector &features, @@ -598,6 +645,243 @@ void BundleAdjusterRay::calcJacobian(Mat &jac) } } +////////////////////////////////////////////////////////////////////////////// + +void BundleAdjusterAffine::setUpInitialCameraParams(const std::vector &cameras) +{ + cam_params_.create(num_images_ * 6, 1, CV_64F); + for (size_t i = 0; i < static_cast(num_images_); ++i) + { + CV_Assert(cameras[i].R.type() == CV_32F); + // cameras[i].R is + // a b tx + // c d ty + // 0 0 1. (optional) + // cam_params_ model for LevMarq is + // (a, b, tx, c, d, ty) + Mat params (2, 3, CV_64F, cam_params_.ptr() + i * 6); + cameras[i].R.rowRange(0, 2).convertTo(params, CV_64F); + } +} + + +void BundleAdjusterAffine::obtainRefinedCameraParams(std::vector &cameras) const +{ + for (int i = 0; i < num_images_; ++i) + { + // cameras[i].R will be + // a b tx + // c d ty + // 0 0 1 + cameras[i].R = Mat::eye(3, 3, CV_32F); + Mat params = cam_params_.rowRange(i * 6, i * 6 + 6).reshape(1, 2); + params.convertTo(cameras[i].R.rowRange(0, 2), CV_32F); + } +} + + +void BundleAdjusterAffine::calcError(Mat &err) +{ + err.create(total_num_matches_ * 2, 1, CV_64F); + + int match_idx = 0; + for (size_t edge_idx = 0; edge_idx < edges_.size(); ++edge_idx) + { + size_t i = edges_[edge_idx].first; + size_t j = edges_[edge_idx].second; + + const ImageFeatures& features1 = features_[i]; + const ImageFeatures& features2 = features_[j]; + const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j]; + + Mat H1 (2, 3, CV_64F, cam_params_.ptr() + i * 6); + Mat H2 (2, 3, CV_64F, cam_params_.ptr() + j * 6); + + // invert H1 + Mat H1_inv; + invertAffineTransform(H1, H1_inv); + + // convert to representation in homogeneous coordinates + Mat last_row = Mat::zeros(1, 3, CV_64F); + last_row.at(2) = 1.; + H1_inv.push_back(last_row); + H2.push_back(last_row); + + Mat_ H = H1_inv * H2; + + for (size_t k = 0; k < matches_info.matches.size(); ++k) + { + if (!matches_info.inliers_mask[k]) + continue; + + const DMatch& m = matches_info.matches[k]; + const Point2f& p1 = features1.keypoints[m.queryIdx].pt; + const Point2f& p2 = features2.keypoints[m.trainIdx].pt; + + double x = H(0,0)*p1.x + H(0,1)*p1.y + H(0,2); + double y = H(1,0)*p1.x + H(1,1)*p1.y + H(1,2); + + err.at(2 * match_idx + 0, 0) = p2.x - x; + err.at(2 * match_idx + 1, 0) = p2.y - y; + + ++match_idx; + } + } +} + + +void BundleAdjusterAffine::calcJacobian(Mat &jac) +{ + jac.create(total_num_matches_ * 2, num_images_ * 6, CV_64F); + + double val; + const double step = 1e-4; + + for (int i = 0; i < num_images_; ++i) + { + for (int j = 0; j < 6; ++j) + { + val = cam_params_.at(i * 6 + j, 0); + cam_params_.at(i * 6 + j, 0) = val - step; + calcError(err1_); + cam_params_.at(i * 6 + j, 0) = val + step; + calcError(err2_); + calcDeriv(err1_, err2_, 2 * step, jac.col(i * 6 + j)); + cam_params_.at(i * 6 + j, 0) = val; + } + } +} + + +////////////////////////////////////////////////////////////////////////////// + +void BundleAdjusterAffinePartial::setUpInitialCameraParams(const std::vector &cameras) +{ + cam_params_.create(num_images_ * 4, 1, CV_64F); + for (size_t i = 0; i < static_cast(num_images_); ++i) + { + CV_Assert(cameras[i].R.type() == CV_32F); + // cameras[i].R is + // a -b tx + // b a ty + // 0 0 1. (optional) + // cam_params_ model for LevMarq is + // (a, b, tx, ty) + double *params = cam_params_.ptr() + i * 4; + params[0] = cameras[i].R.at(0, 0); + params[1] = cameras[i].R.at(1, 0); + params[2] = cameras[i].R.at(0, 2); + params[3] = cameras[i].R.at(1, 2); + } +} + + +void BundleAdjusterAffinePartial::obtainRefinedCameraParams(std::vector &cameras) const +{ + for (size_t i = 0; i < static_cast(num_images_); ++i) + { + // cameras[i].R will be + // a -b tx + // b a ty + // 0 0 1 + // cam_params_ model for LevMarq is + // (a, b, tx, ty) + const double *params = cam_params_.ptr() + i * 4; + double transform_buf[9] = + { + params[0], -params[1], params[2], + params[1], params[0], params[3], + 0., 0., 1. + }; + Mat transform(3, 3, CV_64F, transform_buf); + transform.convertTo(cameras[i].R, CV_32F); + } +} + + +void BundleAdjusterAffinePartial::calcError(Mat &err) +{ + err.create(total_num_matches_ * 2, 1, CV_64F); + + int match_idx = 0; + for (size_t edge_idx = 0; edge_idx < edges_.size(); ++edge_idx) + { + size_t i = edges_[edge_idx].first; + size_t j = edges_[edge_idx].second; + + const ImageFeatures& features1 = features_[i]; + const ImageFeatures& features2 = features_[j]; + const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j]; + + const double *H1_ptr = cam_params_.ptr() + i * 4; + double H1_buf[9] = + { + H1_ptr[0], -H1_ptr[1], H1_ptr[2], + H1_ptr[1], H1_ptr[0], H1_ptr[3], + 0., 0., 1. + }; + Mat H1 (3, 3, CV_64F, H1_buf); + const double *H2_ptr = cam_params_.ptr() + j * 4; + double H2_buf[9] = + { + H2_ptr[0], -H2_ptr[1], H2_ptr[2], + H2_ptr[1], H2_ptr[0], H2_ptr[3], + 0., 0., 1. + }; + Mat H2 (3, 3, CV_64F, H2_buf); + + // invert H1 + Mat H1_aff (H1, Range(0, 2)); + double H1_inv_buf[6]; + Mat H1_inv (2, 3, CV_64F, H1_inv_buf); + invertAffineTransform(H1_aff, H1_inv); + H1_inv.copyTo(H1_aff); + + Mat_ H = H1 * H2; + + for (size_t k = 0; k < matches_info.matches.size(); ++k) + { + if (!matches_info.inliers_mask[k]) + continue; + + const DMatch& m = matches_info.matches[k]; + const Point2f& p1 = features1.keypoints[m.queryIdx].pt; + const Point2f& p2 = features2.keypoints[m.trainIdx].pt; + + double x = H(0,0)*p1.x + H(0,1)*p1.y + H(0,2); + double y = H(1,0)*p1.x + H(1,1)*p1.y + H(1,2); + + err.at(2 * match_idx + 0, 0) = p2.x - x; + err.at(2 * match_idx + 1, 0) = p2.y - y; + + ++match_idx; + } + } +} + + +void BundleAdjusterAffinePartial::calcJacobian(Mat &jac) +{ + jac.create(total_num_matches_ * 2, num_images_ * 4, CV_64F); + + double val; + const double step = 1e-4; + + for (int i = 0; i < num_images_; ++i) + { + for (int j = 0; j < 4; ++j) + { + val = cam_params_.at(i * 4 + j, 0); + cam_params_.at(i * 4 + j, 0) = val - step; + calcError(err1_); + cam_params_.at(i * 4 + j, 0) = val + step; + calcError(err2_); + calcDeriv(err1_, err2_, 2 * step, jac.col(i * 4 + j)); + cam_params_.at(i * 4 + j, 0) = val; + } + } +} + ////////////////////////////////////////////////////////////////////////////// diff --git a/modules/stitching/src/stitcher.cpp b/modules/stitching/src/stitcher.cpp index 19d3924337..905e5a938a 100644 --- a/modules/stitching/src/stitcher.cpp +++ b/modules/stitching/src/stitcher.cpp @@ -91,6 +91,34 @@ Stitcher Stitcher::createDefault(bool try_use_gpu) } +Ptr Stitcher::create(Mode mode, bool try_use_gpu) +{ + Stitcher stit = createDefault(try_use_gpu); + Ptr stitcher = makePtr(stit); + + switch (mode) + { + case PANORAMA: // PANORAMA is the default + // already setup + break; + + case SCANS: + stitcher->setWaveCorrection(false); + stitcher->setFeaturesMatcher(makePtr(false, try_use_gpu)); + stitcher->setBundleAdjuster(makePtr()); + stitcher->setWarper(makePtr()); + stitcher->setExposureCompensator(makePtr()); + break; + + default: + CV_Error(Error::StsBadArg, "Invalid stitching mode. Must be one of Stitcher::Mode"); + break; + } + + return stitcher; +} + + Stitcher::Status Stitcher::estimateTransform(InputArrayOfArrays images) { CV_INSTRUMENT_REGION() @@ -416,6 +444,9 @@ Stitcher::Status Stitcher::matchImages() int64 t = getTickCount(); #endif + std::vector feature_find_imgs(imgs_.size()); + std::vector > feature_find_rois(rois_.size()); + for (size_t i = 0; i < imgs_.size(); ++i) { full_img = imgs_[i]; @@ -444,17 +475,17 @@ Stitcher::Status Stitcher::matchImages() } if (rois_.empty()) - (*features_finder_)(img, features_[i]); + feature_find_imgs[i] = img; else { - std::vector rois(rois_[i].size()); + feature_find_rois[i].resize(rois_[i].size()); for (size_t j = 0; j < rois_[i].size(); ++j) { Point tl(cvRound(rois_[i][j].x * work_scale_), cvRound(rois_[i][j].y * work_scale_)); Point br(cvRound(rois_[i][j].br().x * work_scale_), cvRound(rois_[i][j].br().y * work_scale_)); - rois[j] = Rect(tl, br); + feature_find_rois[i][j] = Rect(tl, br); } - (*features_finder_)(img, features_[i], rois); + feature_find_imgs[i] = img; } features_[i].img_idx = (int)i; LOGLN("Features in image #" << i+1 << ": " << features_[i].keypoints.size()); @@ -463,10 +494,18 @@ Stitcher::Status Stitcher::matchImages() seam_est_imgs_[i] = img.clone(); } + // find features possibly in parallel + if (rois_.empty()) + (*features_finder_)(feature_find_imgs, features_); + else + (*features_finder_)(feature_find_imgs, features_, feature_find_rois); + // Do it to save memory features_finder_->collectGarbage(); full_img.release(); img.release(); + feature_find_imgs.clear(); + feature_find_rois.clear(); LOGLN("Finding features, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); @@ -505,8 +544,16 @@ Stitcher::Status Stitcher::matchImages() Stitcher::Status Stitcher::estimateCameraParams() { - detail::HomographyBasedEstimator estimator; - if (!estimator(features_, pairwise_matches_, cameras_)) + /* TODO OpenCV ABI 4.x + get rid of this dynamic_cast hack and use estimator_ + */ + Ptr estimator; + if (dynamic_cast(features_matcher_.get())) + estimator = makePtr(); + else + estimator = makePtr(); + + if (!(*estimator)(features_, pairwise_matches_, cameras_)) return ERR_HOMOGRAPHY_EST_FAIL; for (size_t i = 0; i < cameras_.size(); ++i) diff --git a/modules/stitching/src/warpers.cpp b/modules/stitching/src/warpers.cpp index 3163fa4414..96fe7f7cb5 100644 --- a/modules/stitching/src/warpers.cpp +++ b/modules/stitching/src/warpers.cpp @@ -222,6 +222,58 @@ void PlaneWarper::detectResultRoi(Size src_size, Point &dst_tl, Point &dst_br) } +Point2f AffineWarper::warpPoint(const Point2f &pt, InputArray K, InputArray H) +{ + Mat R, T; + getRTfromHomogeneous(H, R, T); + return PlaneWarper::warpPoint(pt, K, R, T); +} + + +Rect AffineWarper::buildMaps(Size src_size, InputArray K, InputArray H, OutputArray xmap, OutputArray ymap) +{ + Mat R, T; + getRTfromHomogeneous(H, R, T); + return PlaneWarper::buildMaps(src_size, K, R, T, xmap, ymap); +} + + +Point AffineWarper::warp(InputArray src, InputArray K, InputArray H, + int interp_mode, int border_mode, OutputArray dst) +{ + Mat R, T; + getRTfromHomogeneous(H, R, T); + return PlaneWarper::warp(src, K, R, T, interp_mode, border_mode, dst); +} + + +Rect AffineWarper::warpRoi(Size src_size, InputArray K, InputArray H) +{ + Mat R, T; + getRTfromHomogeneous(H, R, T); + return PlaneWarper::warpRoi(src_size, K, R, T); +} + + +void AffineWarper::getRTfromHomogeneous(InputArray H_, Mat &R, Mat &T) +{ + Mat H = H_.getMat(); + CV_Assert(H.size() == Size(3, 3) && H.type() == CV_32F); + + T = Mat::zeros(3, 1, CV_32F); + R = H.clone(); + + T.at(0,0) = R.at(0,2); + T.at(1,0) = R.at(1,2); + R.at(0,2) = 0.f; + R.at(1,2) = 0.f; + + // we want to compensate transform to fit into plane warper + R = R.t(); + T = (R * T) * -1; +} + + void SphericalWarper::detectResultRoi(Size src_size, Point &dst_tl, Point &dst_br) { detectResultRoiByBorder(src_size, dst_tl, dst_br); diff --git a/modules/stitching/test/ocl/test_warpers.cpp b/modules/stitching/test/ocl/test_warpers.cpp index a42ebb36c2..2b74372414 100644 --- a/modules/stitching/test/ocl/test_warpers.cpp +++ b/modules/stitching/test/ocl/test_warpers.cpp @@ -143,6 +143,27 @@ OCL_TEST_F(PlaneWarperTest, Mat) } } +typedef WarperTestBase AffineWarperTest; + +OCL_TEST_F(AffineWarperTest, Mat) +{ + for (int j = 0; j < test_loop_times; j++) + { + generateTestData(); + + Ptr creator = makePtr(); + Ptr warper = creator->create(1.0); + + OCL_OFF(warper->buildMaps(src.size(), K, R, xmap, ymap)); + OCL_ON(warper->buildMaps(usrc.size(), K, R, uxmap, uymap)); + + OCL_OFF(warper->warp(src, K, R, INTER_LINEAR, BORDER_REPLICATE, dst)); + OCL_ON(warper->warp(usrc, K, R, INTER_LINEAR, BORDER_REPLICATE, udst)); + + Near(1.5e-4); + } +} + } } // namespace cvtest::ocl #endif // HAVE_OPENCL diff --git a/modules/stitching/test/test_matchers.cpp b/modules/stitching/test/test_matchers.cpp index 6deed7b75f..61d86a6a0f 100644 --- a/modules/stitching/test/test_matchers.cpp +++ b/modules/stitching/test/test_matchers.cpp @@ -42,11 +42,11 @@ #include "test_precomp.hpp" #include "opencv2/opencv_modules.hpp" -#ifdef HAVE_OPENCV_XFEATURES2D - using namespace cv; using namespace std; +#ifdef HAVE_OPENCV_XFEATURES2D + TEST(SurfFeaturesFinder, CanFindInROIs) { Ptr finder = makePtr(); @@ -75,4 +75,27 @@ TEST(SurfFeaturesFinder, CanFindInROIs) ASSERT_EQ(bad_count, 0); } -#endif +#endif // HAVE_OPENCV_XFEATURES2D + +TEST(ParallelFeaturesFinder, IsSameWithSerial) +{ + Ptr para_finder = makePtr(); + Ptr serial_finder = makePtr(); + Mat img = imread(string(cvtest::TS::ptr()->get_data_path()) + "stitching/a3.png", IMREAD_GRAYSCALE); + + vector imgs(50, img); + detail::ImageFeatures serial_features; + vector para_features(imgs.size()); + + (*serial_finder)(img, serial_features); + (*para_finder)(imgs, para_features); + + // results must be the same + for(size_t i = 0; i < para_features.size(); ++i) + { + Mat diff_descriptors = serial_features.descriptors.getMat(ACCESS_READ) != para_features[i].descriptors.getMat(ACCESS_READ); + ASSERT_EQ(countNonZero(diff_descriptors), 0); + ASSERT_EQ(serial_features.img_size, para_features[i].img_size); + ASSERT_EQ(serial_features.keypoints.size(), para_features[i].keypoints.size()); + } +} diff --git a/modules/video/include/opencv2/video/tracking.hpp b/modules/video/include/opencv2/video/tracking.hpp index d4ccd532aa..c6ead3a42c 100644 --- a/modules/video/include/opencv2/video/tracking.hpp +++ b/modules/video/include/opencv2/video/tracking.hpp @@ -245,7 +245,7 @@ where src[i] and dst[i] are the i-th points in src and dst, respectively when fullAffine=false. @sa -getAffineTransform, getPerspectiveTransform, findHomography +estimateAffine2D, estimateAffinePartial2D, getAffineTransform, getPerspectiveTransform, findHomography */ CV_EXPORTS_W Mat estimateRigidTransform( InputArray src, InputArray dst, bool fullAffine ); @@ -306,7 +306,7 @@ sample image_alignment.cpp that demonstrates the use of the function. Note that an exception if algorithm does not converges. @sa -estimateRigidTransform, findHomography +estimateAffine2D, estimateAffinePartial2D, findHomography */ CV_EXPORTS_W double findTransformECC( InputArray templateImage, InputArray inputImage, InputOutputArray warpMatrix, int motionType = MOTION_AFFINE, diff --git a/samples/cpp/stitching.cpp b/samples/cpp/stitching.cpp index 5b4437ac9e..d0b38b4e5c 100644 --- a/samples/cpp/stitching.cpp +++ b/samples/cpp/stitching.cpp @@ -50,6 +50,7 @@ using namespace std; using namespace cv; bool try_use_gpu = false; +Stitcher::Mode mode = Stitcher::PANORAMA; vector imgs; string result_name = "result.jpg"; @@ -62,8 +63,8 @@ int main(int argc, char* argv[]) if (retval) return -1; Mat pano; - Stitcher stitcher = Stitcher::createDefault(try_use_gpu); - Stitcher::Status status = stitcher.stitch(imgs, pano); + Ptr stitcher = Stitcher::create(mode, try_use_gpu); + Stitcher::Status status = stitcher->stitch(imgs, pano); if (status != Stitcher::OK) { @@ -79,12 +80,16 @@ int main(int argc, char* argv[]) void printUsage() { cout << - "Rotation model images stitcher.\n\n" + "Images stitcher.\n\n" "stitching img1 img2 [...imgN]\n\n" "Flags:\n" " --try_use_gpu (yes|no)\n" " Try to use GPU. The default value is 'no'. All default values\n" " are for CPU mode.\n" + " --mode (panorama|scans)\n" + " Determines configuration of stitcher. The default is 'panorama',\n" + " mode suitable for creating photo panoramas. Option 'scans' is suitable\n" + " for stitching materials under affine transformation, such as scans.\n" " --output \n" " The default is 'result.jpg'.\n"; } @@ -122,6 +127,19 @@ int parseCmdArgs(int argc, char** argv) result_name = argv[i + 1]; i++; } + else if (string(argv[i]) == "--mode") + { + if (string(argv[i + 1]) == "panorama") + mode = Stitcher::PANORAMA; + else if (string(argv[i + 1]) == "scans") + mode = Stitcher::SCANS; + else + { + cout << "Bad --mode flag value\n"; + return -1; + } + i++; + } else { Mat img = imread(argv[i]); diff --git a/samples/cpp/stitching_detailed.cpp b/samples/cpp/stitching_detailed.cpp index f333bc6e36..df412e14fa 100644 --- a/samples/cpp/stitching_detailed.cpp +++ b/samples/cpp/stitching_detailed.cpp @@ -84,12 +84,16 @@ static void printUsage() " Resolution for image registration step. The default is 0.6 Mpx.\n" " --features (surf|orb)\n" " Type of features used for images matching. The default is surf.\n" + " --matcher (homography|affine)\n" + " Matcher used for pairwise image matching.\n" + " --estimator (homography|affine)\n" + " Type of estimator used for transformation estimation.\n" " --match_conf \n" " Confidence for feature matching step. The default is 0.65 for surf and 0.3 for orb.\n" " --conf_thresh \n" " Threshold for two images are from the same panorama confidence.\n" " The default is 1.0.\n" - " --ba (reproj|ray)\n" + " --ba (no|reproj|ray|affine)\n" " Bundle adjustment cost function. The default is ray.\n" " --ba_refine_mask (mask)\n" " Set refinement mask for bundle adjustment. It looks like 'x_xxx',\n" @@ -105,7 +109,7 @@ static void printUsage() " Labels description: Nm is number of matches, Ni is number of inliers,\n" " C is confidence.\n" "\nCompositing Flags:\n" - " --warp (plane|cylindrical|spherical|fisheye|stereographic|compressedPlaneA2B1|compressedPlaneA1.5B1|compressedPlanePortraitA2B1|compressedPlanePortraitA1.5B1|paniniA2B1|paniniA1.5B1|paniniPortraitA2B1|paniniPortraitA1.5B1|mercator|transverseMercator)\n" + " --warp (affine|plane|cylindrical|spherical|fisheye|stereographic|compressedPlaneA2B1|compressedPlaneA1.5B1|compressedPlanePortraitA2B1|compressedPlanePortraitA1.5B1|paniniA2B1|paniniA1.5B1|paniniPortraitA2B1|paniniPortraitA1.5B1|mercator|transverseMercator)\n" " Warp surface type. The default is 'spherical'.\n" " --seam_megapix \n" " Resolution for seam estimation step. The default is 0.1 Mpx.\n" @@ -138,6 +142,8 @@ double seam_megapix = 0.1; double compose_megapix = -1; float conf_thresh = 1.f; string features_type = "surf"; +string matcher_type = "homography"; +string estimator_type = "homography"; string ba_cost_func = "ray"; string ba_refine_mask = "xxxxx"; bool do_wave_correct = true; @@ -214,6 +220,28 @@ static int parseCmdArgs(int argc, char** argv) match_conf = 0.3f; i++; } + else if (string(argv[i]) == "--matcher") + { + if (string(argv[i + 1]) == "homography" || string(argv[i + 1]) == "affine") + matcher_type = argv[i + 1]; + else + { + cout << "Bad --matcher flag value\n"; + return -1; + } + i++; + } + else if (string(argv[i]) == "--estimator") + { + if (string(argv[i + 1]) == "homography" || string(argv[i + 1]) == "affine") + estimator_type = argv[i + 1]; + else + { + cout << "Bad --estimator flag value\n"; + return -1; + } + i++; + } else if (string(argv[i]) == "--match_conf") { match_conf = static_cast(atof(argv[i + 1])); @@ -465,18 +493,16 @@ int main(int argc, char* argv[]) t = getTickCount(); #endif vector pairwise_matches; - if (range_width==-1) - { - BestOf2NearestMatcher matcher(try_cuda, match_conf); - matcher(features, pairwise_matches); - matcher.collectGarbage(); - } + Ptr matcher; + if (matcher_type == "affine") + matcher = makePtr(false, try_cuda, match_conf); + else if (range_width==-1) + matcher = makePtr(try_cuda, match_conf); else - { - BestOf2NearestRangeMatcher matcher(range_width, try_cuda, match_conf); - matcher(features, pairwise_matches); - matcher.collectGarbage(); - } + matcher = makePtr(range_width, try_cuda, match_conf); + + (*matcher)(features, pairwise_matches); + matcher->collectGarbage(); LOGLN("Pairwise matching, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); @@ -512,9 +538,14 @@ int main(int argc, char* argv[]) return -1; } - HomographyBasedEstimator estimator; + Ptr estimator; + if (estimator_type == "affine") + estimator = makePtr(); + else + estimator = makePtr(); + vector cameras; - if (!estimator(features, pairwise_matches, cameras)) + if (!(*estimator)(features, pairwise_matches, cameras)) { cout << "Homography estimation failed.\n"; return -1; @@ -525,12 +556,14 @@ int main(int argc, char* argv[]) Mat R; cameras[i].R.convertTo(R, CV_32F); cameras[i].R = R; - LOGLN("Initial intrinsics #" << indices[i]+1 << ":\n" << cameras[i].K()); + LOGLN("Initial camera intrinsics #" << indices[i]+1 << ":\nK:\n" << cameras[i].K() << "\nR:\n" << cameras[i].R); } Ptr adjuster; if (ba_cost_func == "reproj") adjuster = makePtr(); else if (ba_cost_func == "ray") adjuster = makePtr(); + else if (ba_cost_func == "affine") adjuster = makePtr(); + else if (ba_cost_func == "no") adjuster = makePtr(); else { cout << "Unknown bundle adjustment cost function: '" << ba_cost_func << "'.\n"; @@ -555,7 +588,7 @@ int main(int argc, char* argv[]) vector focals; for (size_t i = 0; i < cameras.size(); ++i) { - LOGLN("Camera #" << indices[i]+1 << ":\n" << cameras[i].K()); + LOGLN("Camera #" << indices[i]+1 << ":\nK:\n" << cameras[i].K() << "\nR:\n" << cameras[i].R); focals.push_back(cameras[i].focal); } @@ -612,6 +645,8 @@ int main(int argc, char* argv[]) { if (warp_type == "plane") warper_creator = makePtr(); + else if (warp_type == "affine") + warper_creator = makePtr(); else if (warp_type == "cylindrical") warper_creator = makePtr(); else if (warp_type == "spherical")