turned opencv_stitching application to module and sample

modules/stitching/src/autocalib.cpp

@@ -0,0 +1,187 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok, bool &f1_ok)
+{
+    CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3));
+
+    const double* h = reinterpret_cast<const double*>(H.data);
+
+    double d1, d2; // Denominators
+    double v1, v2; // Focal squares value candidates
+
+    f1_ok = true;
+    d1 = h[6] * h[7];
+    d2 = (h[7] - h[6]) * (h[7] + h[6]);
+    v1 = -(h[0] * h[1] + h[3] * h[4]) / d1;
+    v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2;
+    if (v1 < v2) std::swap(v1, v2);
+    if (v1 > 0 && v2 > 0) f1 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
+    else if (v1 > 0) f1 = sqrt(v1);
+    else f1_ok = false;
+
+    f0_ok = true;
+    d1 = h[0] * h[3] + h[1] * h[4];
+    d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];
+    v1 = -h[2] * h[5] / d1;
+    v2 = (h[5] * h[5] - h[2] * h[2]) / d2;
+    if (v1 < v2) std::swap(v1, v2);
+    if (v1 > 0 && v2 > 0) f0 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
+    else if (v1 > 0) f0 = sqrt(v1);
+    else f0_ok = false;
+}
+
+
+void estimateFocal(const vector<ImageFeatures> &features, const vector<MatchesInfo> &pairwise_matches, 
+                   vector<double> &focals)
+{
+    const int num_images = static_cast<int>(features.size());
+    focals.resize(num_images);
+
+    vector<double> all_focals;
+
+    for (int i = 0; i < num_images; ++i)
+    {        
+        for (int j = 0; j < num_images; ++j)
+        {
+            const MatchesInfo &m = pairwise_matches[i*num_images + j];
+            if (m.H.empty())
+                continue;
+            double f0, f1;
+            bool f0ok, f1ok;
+            focalsFromHomography(m.H, f0, f1, f0ok, f1ok);
+            if (f0ok && f1ok) 
+                all_focals.push_back(sqrt(f0 * f1));
+        }
+    }
+
+    if (static_cast<int>(all_focals.size()) >= num_images - 1)
+    {
+        nth_element(all_focals.begin(), all_focals.begin() + all_focals.size()/2, all_focals.end());
+        for (int i = 0; i < num_images; ++i)
+            focals[i] = all_focals[all_focals.size()/2];
+    }
+    else
+    {
+        LOGLN("Can't estimate focal length, will use naive approach");
+        double focals_sum = 0;
+        for (int i = 0; i < num_images; ++i)
+            focals_sum += features[i].img_size.width + features[i].img_size.height;
+        for (int i = 0; i < num_images; ++i)
+            focals[i] = focals_sum / num_images;
+    }
+}
+
+
+namespace
+{
+    template<typename _Tp> static inline bool
+    decomposeCholesky(_Tp* A, size_t astep, int m)
+    {
+        if (!Cholesky(A, astep, m, 0, 0, 0))
+            return false;
+        astep /= sizeof(A[0]);
+        for (int i = 0; i < m; ++i)
+            A[i*astep + i] = (_Tp)(1./A[i*astep + i]);
+        return true;
+    }
+} // namespace
+
+
+bool calibrateRotatingCamera(const vector<Mat> &Hs, Mat &K)
+{
+    int m = static_cast<int>(Hs.size());
+    CV_Assert(m >= 1);
+
+    vector<Mat> Hs_(m);
+    for (int i = 0; i < m; ++i)
+    {
+        CV_Assert(Hs[i].size() == Size(3, 3) && Hs[i].type() == CV_64F);
+        Hs_[i] = Hs[i] / pow(determinant(Hs[i]), 1./3.);
+    }
+
+    const int idx_map[3][3] = {{0, 1, 2}, {1, 3, 4}, {2, 4, 5}};
+    Mat_<double> A(6*m, 6);
+    A.setTo(0);
+
+    int eq_idx = 0;
+    for (int k = 0; k < m; ++k)
+    {
+        Mat_<double> H(Hs_[k]);
+        for (int i = 0; i < 3; ++i)
+        {
+            for (int j = i; j < 3; ++j, ++eq_idx)
+            {
+                for (int l = 0; l < 3; ++l)
+                {
+                    for (int s = 0; s < 3; ++s)
+                    {
+                        int idx = idx_map[l][s];
+                        A(eq_idx, idx) += H(i,l) * H(j,s);
+                    }
+                }
+                A(eq_idx, idx_map[i][j]) -= 1;
+            }
+        }
+    }
+
+    Mat_<double> wcoef;
+    SVD::solveZ(A, wcoef);
+
+    Mat_<double> W(3,3);
+    for (int i = 0; i < 3; ++i)
+        for (int j = i; j < 3; ++j)
+            W(i,j) = W(j,i) = wcoef(idx_map[i][j], 0) / wcoef(5,0);
+    if (!decomposeCholesky(W.ptr<double>(), W.step, 3))
+        return false;
+    W(0,1) = W(0,2) = W(1,2) = 0;
+    K = W.t();
+    return true;
+}
+
+} // namespace cv

modules/stitching/src/blenders.cpp

@@ -0,0 +1,378 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+static const float WEIGHT_EPS = 1e-5f;
+
+Ptr<Blender> Blender::createDefault(int type, bool try_gpu)
+{
+    if (type == NO)
+        return new Blender();
+    if (type == FEATHER)
+        return new FeatherBlender();
+    if (type == MULTI_BAND)
+        return new MultiBandBlender(try_gpu);
+    CV_Error(CV_StsBadArg, "unsupported blending method");
+    return NULL;
+}
+
+
+void Blender::prepare(const vector<Point> &corners, const vector<Size> &sizes)
+{
+    prepare(resultRoi(corners, sizes));
+}
+
+
+void Blender::prepare(Rect dst_roi)
+{
+    dst_.create(dst_roi.size(), CV_16SC3);
+    dst_.setTo(Scalar::all(0));
+    dst_mask_.create(dst_roi.size(), CV_8U);
+    dst_mask_.setTo(Scalar::all(0));
+    dst_roi_ = dst_roi;
+}
+
+
+void Blender::feed(const Mat &img, const Mat &mask, Point tl) 
+{
+    CV_Assert(img.type() == CV_16SC3);
+    CV_Assert(mask.type() == CV_8U);
+    int dx = tl.x - dst_roi_.x;
+    int dy = tl.y - dst_roi_.y;
+
+    for (int y = 0; y < img.rows; ++y)
+    {
+        const Point3_<short> *src_row = img.ptr<Point3_<short> >(y);
+        Point3_<short> *dst_row = dst_.ptr<Point3_<short> >(dy + y);
+        const uchar *mask_row = mask.ptr<uchar>(y);
+        uchar *dst_mask_row = dst_mask_.ptr<uchar>(dy + y);
+
+        for (int x = 0; x < img.cols; ++x)
+        {
+            if (mask_row[x]) 
+                dst_row[dx + x] = src_row[x];
+            dst_mask_row[dx + x] |= mask_row[x];
+        }
+    }
+}
+
+
+void Blender::blend(Mat &dst, Mat &dst_mask)
+{
+    dst_.setTo(Scalar::all(0), dst_mask_ == 0);
+    dst = dst_;
+    dst_mask = dst_mask_;
+    dst_.release();
+    dst_mask_.release();
+}
+
+
+void FeatherBlender::prepare(Rect dst_roi)
+{
+    Blender::prepare(dst_roi);
+    dst_weight_map_.create(dst_roi.size(), CV_32F);
+    dst_weight_map_.setTo(0);
+}
+
+
+void FeatherBlender::feed(const Mat &img, const Mat &mask, Point tl)
+{
+    CV_Assert(img.type() == CV_16SC3);
+    CV_Assert(mask.type() == CV_8U);
+
+    createWeightMap(mask, sharpness_, weight_map_);
+    int dx = tl.x - dst_roi_.x;
+    int dy = tl.y - dst_roi_.y;
+
+    for (int y = 0; y < img.rows; ++y)
+    {
+        const Point3_<short>* src_row = img.ptr<Point3_<short> >(y);
+        Point3_<short>* dst_row = dst_.ptr<Point3_<short> >(dy + y);
+        const float* weight_row = weight_map_.ptr<float>(y);
+        float* dst_weight_row = dst_weight_map_.ptr<float>(dy + y);
+
+        for (int x = 0; x < img.cols; ++x)               
+        {
+            dst_row[dx + x].x += static_cast<short>(src_row[x].x * weight_row[x]);
+            dst_row[dx + x].y += static_cast<short>(src_row[x].y * weight_row[x]);
+            dst_row[dx + x].z += static_cast<short>(src_row[x].z * weight_row[x]);
+            dst_weight_row[dx + x] += weight_row[x];
+        }
+    }
+}
+
+
+void FeatherBlender::blend(Mat &dst, Mat &dst_mask)
+{
+    normalizeUsingWeightMap(dst_weight_map_, dst_);
+    dst_mask_ = dst_weight_map_ > WEIGHT_EPS;
+    Blender::blend(dst, dst_mask);
+}
+
+
+MultiBandBlender::MultiBandBlender(int try_gpu, int num_bands)
+{
+    setNumBands(num_bands);
+    can_use_gpu_ = try_gpu && gpu::getCudaEnabledDeviceCount();
+}
+
+
+void MultiBandBlender::prepare(Rect dst_roi)
+{
+    dst_roi_final_ = dst_roi;
+
+    // Crop unnecessary bands
+    double max_len = static_cast<double>(max(dst_roi.width, dst_roi.height));
+    num_bands_ = min(actual_num_bands_, static_cast<int>(ceil(log(max_len) / log(2.0))));
+
+    // Add border to the final image, to ensure sizes are divided by (1 << num_bands_)
+    dst_roi.width += ((1 << num_bands_) - dst_roi.width % (1 << num_bands_)) % (1 << num_bands_);
+    dst_roi.height += ((1 << num_bands_) - dst_roi.height % (1 << num_bands_)) % (1 << num_bands_);
+
+    Blender::prepare(dst_roi);
+
+    dst_pyr_laplace_.resize(num_bands_ + 1);
+    dst_pyr_laplace_[0] = dst_;
+
+    dst_band_weights_.resize(num_bands_ + 1);
+    dst_band_weights_[0].create(dst_roi.size(), CV_32F);
+    dst_band_weights_[0].setTo(0);
+
+    for (int i = 1; i <= num_bands_; ++i)
+    {
+        dst_pyr_laplace_[i].create((dst_pyr_laplace_[i - 1].rows + 1) / 2, 
+                                   (dst_pyr_laplace_[i - 1].cols + 1) / 2, CV_16SC3);
+        dst_band_weights_[i].create((dst_band_weights_[i - 1].rows + 1) / 2,
+                                    (dst_band_weights_[i - 1].cols + 1) / 2, CV_32F);
+        dst_pyr_laplace_[i].setTo(Scalar::all(0));
+        dst_band_weights_[i].setTo(0);
+    }
+}
+
+
+void MultiBandBlender::feed(const Mat &img, const Mat &mask, Point tl)
+{
+    CV_Assert(img.type() == CV_16SC3);
+    CV_Assert(mask.type() == CV_8U);
+
+    // Keep source image in memory with small border
+    int gap = 3 * (1 << num_bands_);
+    Point tl_new(max(dst_roi_.x, tl.x - gap), 
+                 max(dst_roi_.y, tl.y - gap));
+    Point br_new(min(dst_roi_.br().x, tl.x + img.cols + gap), 
+                 min(dst_roi_.br().y, tl.y + img.rows + gap));
+
+    // Ensure coordinates of top-left, bottom-right corners are divided by (1 << num_bands_). 
+    // After that scale between layers is exactly 2.
+    //
+    // We do it to avoid interpolation problems when keeping sub-images only. There is no such problem when 
+    // image is bordered to have size equal to the final image size, but this is too memory hungry approach.
+    tl_new.x = dst_roi_.x + (((tl_new.x - dst_roi_.x) >> num_bands_) << num_bands_);
+    tl_new.y = dst_roi_.y + (((tl_new.y - dst_roi_.y) >> num_bands_) << num_bands_);
+    int width = br_new.x - tl_new.x;
+    int height = br_new.y - tl_new.y;
+    width += ((1 << num_bands_) - width % (1 << num_bands_)) % (1 << num_bands_);
+    height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_bands_);
+    br_new.x = tl_new.x + width;
+    br_new.y = tl_new.y + height;
+    int dy = max(br_new.y - dst_roi_.br().y, 0);
+    int dx = max(br_new.x - dst_roi_.br().x, 0);
+    tl_new.x -= dx; br_new.x -= dx;
+    tl_new.y -= dy; br_new.y -= dy;
+
+    int top = tl.y - tl_new.y;
+    int left = tl.x - tl_new.x;
+    int bottom = br_new.y - tl.y - img.rows;
+    int right = br_new.x - tl.x - img.cols;
+
+    // Create the source image Laplacian pyramid
+    Mat img_with_border;
+    copyMakeBorder(img, img_with_border, top, bottom, left, right,
+                   BORDER_REFLECT);
+    vector<Mat> src_pyr_laplace;
+    if (can_use_gpu_)
+        createLaplacePyrGpu(img_with_border, num_bands_, src_pyr_laplace);
+    else
+        createLaplacePyr(img_with_border, num_bands_, src_pyr_laplace);
+
+    // Create the weight map Gaussian pyramid
+    Mat weight_map;
+    mask.convertTo(weight_map, CV_32F, 1./255.);
+    vector<Mat> weight_pyr_gauss(num_bands_ + 1);
+    copyMakeBorder(weight_map, weight_pyr_gauss[0], top, bottom, left, right, 
+                   BORDER_CONSTANT);
+    for (int i = 0; i < num_bands_; ++i)
+        pyrDown(weight_pyr_gauss[i], weight_pyr_gauss[i + 1]);
+
+    int y_tl = tl_new.y - dst_roi_.y;
+    int y_br = br_new.y - dst_roi_.y;
+    int x_tl = tl_new.x - dst_roi_.x;
+    int x_br = br_new.x - dst_roi_.x;
+
+    // Add weighted layer of the source image to the final Laplacian pyramid layer
+    for (int i = 0; i <= num_bands_; ++i)
+    {
+        for (int y = y_tl; y < y_br; ++y)
+        {
+            int y_ = y - y_tl;
+            const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Point3_<short> >(y_);
+            Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(y);
+            const float* weight_row = weight_pyr_gauss[i].ptr<float>(y_);
+            float* dst_weight_row = dst_band_weights_[i].ptr<float>(y);
+
+            for (int x = x_tl; x < x_br; ++x)               
+            {
+                int x_ = x - x_tl;
+                dst_row[x].x += static_cast<short>(src_row[x_].x * weight_row[x_]);
+                dst_row[x].y += static_cast<short>(src_row[x_].y * weight_row[x_]);
+                dst_row[x].z += static_cast<short>(src_row[x_].z * weight_row[x_]);
+                dst_weight_row[x] += weight_row[x_];
+            }
+        }
+        x_tl /= 2; y_tl /= 2; 
+        x_br /= 2; y_br /= 2;
+    }
+}
+
+
+void MultiBandBlender::blend(Mat &dst, Mat &dst_mask)
+{
+    for (int i = 0; i <= num_bands_; ++i)
+        normalizeUsingWeightMap(dst_band_weights_[i], dst_pyr_laplace_[i]);
+
+    restoreImageFromLaplacePyr(dst_pyr_laplace_);
+
+    dst_ = dst_pyr_laplace_[0];
+    dst_ = dst_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width));
+    dst_mask_ = dst_band_weights_[0] > WEIGHT_EPS;
+    dst_mask_ = dst_mask_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width));
+    dst_pyr_laplace_.clear();
+    dst_band_weights_.clear();
+
+    Blender::blend(dst, dst_mask);
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+// Auxiliary functions
+
+void normalizeUsingWeightMap(const Mat& weight, Mat& src)
+{
+    CV_Assert(weight.type() == CV_32F);
+    CV_Assert(src.type() == CV_16SC3);
+    for (int y = 0; y < src.rows; ++y)
+    {
+        Point3_<short> *row = src.ptr<Point3_<short> >(y);
+        const float *weight_row = weight.ptr<float>(y);
+
+        for (int x = 0; x < src.cols; ++x)
+        {
+            row[x].x = static_cast<short>(row[x].x / (weight_row[x] + WEIGHT_EPS));
+            row[x].y = static_cast<short>(row[x].y / (weight_row[x] + WEIGHT_EPS));
+            row[x].z = static_cast<short>(row[x].z / (weight_row[x] + WEIGHT_EPS));
+        }
+    }
+}
+
+
+void createWeightMap(const Mat &mask, float sharpness, Mat &weight)
+{
+    CV_Assert(mask.type() == CV_8U);
+    distanceTransform(mask, weight, CV_DIST_L1, 3);
+    threshold(weight * sharpness, weight, 1.f, 1.f, THRESH_TRUNC);
+}
+
+
+void createLaplacePyr(const Mat &img, int num_levels, vector<Mat> &pyr)
+{
+    pyr.resize(num_levels + 1);
+    pyr[0] = img;
+    for (int i = 0; i < num_levels; ++i)
+        pyrDown(pyr[i], pyr[i + 1]);
+    Mat tmp;
+    for (int i = 0; i < num_levels; ++i)
+    {
+        pyrUp(pyr[i + 1], tmp, pyr[i].size());
+        subtract(pyr[i], tmp, pyr[i]);
+    }
+}
+
+
+void createLaplacePyrGpu(const Mat &img, int num_levels, vector<Mat> &pyr)
+{
+    pyr.resize(num_levels + 1);
+
+    vector<gpu::GpuMat> gpu_pyr(num_levels + 1);
+    gpu_pyr[0] = img;
+    for (int i = 0; i < num_levels; ++i)
+        gpu::pyrDown(gpu_pyr[i], gpu_pyr[i + 1]);
+
+    gpu::GpuMat tmp;
+    for (int i = 0; i < num_levels; ++i)
+    {
+        gpu::pyrUp(gpu_pyr[i + 1], tmp);
+        gpu::subtract(gpu_pyr[i], tmp, gpu_pyr[i]);
+        pyr[i] = gpu_pyr[i];
+    }
+
+    pyr[num_levels] = gpu_pyr[num_levels];
+}
+
+
+void restoreImageFromLaplacePyr(vector<Mat> &pyr)
+{
+    if (pyr.size() == 0)
+        return;
+    Mat tmp;
+    for (size_t i = pyr.size() - 1; i > 0; --i)
+    {
+        pyrUp(pyr[i], tmp, pyr[i - 1].size());
+        add(tmp, pyr[i - 1], pyr[i - 1]);
+    }
+}
+
+} // namespace cv

modules/stitching/src/camera.cpp

@@ -0,0 +1,20 @@
+#include "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+CameraParams::CameraParams() : focal(1), R(Mat::eye(3, 3, CV_64F)), t(Mat::zeros(3, 1, CV_64F)) {}
+
+CameraParams::CameraParams(const CameraParams &other) { *this = other; }
+
+const CameraParams& CameraParams::operator =(const CameraParams &other)
+{
+    focal = other.focal;
+    R = other.R.clone();
+    t = other.t.clone();
+    return *this;
+}
+
+} // namespace cv

modules/stitching/src/exposure_compensate.cpp

@@ -0,0 +1,246 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+using namespace cv::gpu;
+
+namespace cv
+{
+
+Ptr<ExposureCompensator> ExposureCompensator::createDefault(int type)
+{
+    if (type == NO)
+        return new NoExposureCompensator();
+    if (type == GAIN)
+        return new GainCompensator();
+    if (type == GAIN_BLOCKS)
+        return new BlocksGainCompensator();
+    CV_Error(CV_StsBadArg, "unsupported exposure compensation method");
+    return NULL;
+}
+
+
+void ExposureCompensator::feed(const vector<Point> &corners, const vector<Mat> &images, 
+                               const vector<Mat> &masks)
+{
+    vector<pair<Mat,uchar> > level_masks;
+    for (size_t i = 0; i < masks.size(); ++i)
+        level_masks.push_back(make_pair(masks[i], 255));
+    feed(corners, images, level_masks);
+}
+
+
+void GainCompensator::feed(const vector<Point> &corners, const vector<Mat> &images, 
+                           const vector<pair<Mat,uchar> > &masks)
+{
+    LOGLN("Exposure compensation...");
+    int64 t = getTickCount();
+
+    CV_Assert(corners.size() == images.size() && images.size() == masks.size());
+
+    const int num_images = static_cast<int>(images.size());
+    Mat_<int> N(num_images, num_images); N.setTo(0);
+    Mat_<double> I(num_images, num_images); I.setTo(0);
+
+    Rect dst_roi = resultRoi(corners, images);
+    Mat subimg1, subimg2;
+    Mat_<uchar> submask1, submask2, intersect;
+
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = i; j < num_images; ++j)
+        {
+            Rect roi;
+            if (overlapRoi(corners[i], corners[j], images[i].size(), images[j].size(), roi))
+            {
+                subimg1 = images[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));
+                subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));
+
+                submask1 = masks[i].first(Rect(roi.tl() - corners[i], roi.br() - corners[i]));
+                submask2 = masks[j].first(Rect(roi.tl() - corners[j], roi.br() - corners[j]));
+                intersect = (submask1 == masks[i].second) & (submask2 == masks[j].second);
+
+                N(i, j) = N(j, i) = max(1, countNonZero(intersect));
+
+                double Isum1 = 0, Isum2 = 0;
+                for (int y = 0; y < roi.height; ++y)
+                {
+                    const Point3_<uchar>* r1 = subimg1.ptr<Point3_<uchar> >(y);
+                    const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(y);
+                    for (int x = 0; x < roi.width; ++x)
+                    {
+                        if (intersect(y, x))
+                        {
+                            Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z)));
+                            Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) + sqr(r2[x].y) + sqr(r2[x].z)));
+                        }
+                    }
+                }
+                I(i, j) = Isum1 / N(i, j);
+                I(j, i) = Isum2 / N(i, j);
+            }
+        }
+    }
+
+    double alpha = 0.01;
+    double beta = 100;
+
+    Mat_<double> A(num_images, num_images); A.setTo(0);
+    Mat_<double> b(num_images, 1); b.setTo(0);
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = 0; j < num_images; ++j)
+        {
+            b(i, 0) += beta * N(i, j);
+            A(i, i) += beta * N(i, j);
+            if (j == i) continue;
+            A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);
+            A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);
+        }
+    }
+
+    solve(A, b, gains_);
+
+    LOGLN("Exposure compensation, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
+}
+
+
+void GainCompensator::apply(int index, Point /*corner*/, Mat &image, const Mat &/*mask*/)
+{
+    image *= gains_(index, 0);
+}
+
+
+vector<double> GainCompensator::gains() const
+{
+    vector<double> gains_vec(gains_.rows);
+    for (int i = 0; i < gains_.rows; ++i)
+        gains_vec[i] = gains_(i, 0);
+    return gains_vec;
+}
+
+
+void BlocksGainCompensator::feed(const vector<Point> &corners, const vector<Mat> &images, 
+                                const vector<pair<Mat,uchar> > &masks)
+{
+    CV_Assert(corners.size() == images.size() && images.size() == masks.size());
+
+    const int num_images = static_cast<int>(images.size());
+
+    vector<Size> bl_per_imgs(num_images);
+    vector<Point> block_corners;
+    vector<Mat> block_images;
+    vector<pair<Mat,uchar> > block_masks;
+
+    // Construct blocks for gain compensator
+    for (int img_idx = 0; img_idx < num_images; ++img_idx)
+    {
+        Size bl_per_img((images[img_idx].cols + bl_width_ - 1) / bl_width_,
+                        (images[img_idx].rows + bl_height_ - 1) / bl_height_);
+        int bl_width = (images[img_idx].cols + bl_per_img.width - 1) / bl_per_img.width;
+        int bl_height = (images[img_idx].rows + bl_per_img.height - 1) / bl_per_img.height;
+        bl_per_imgs[img_idx] = bl_per_img;
+        for (int by = 0; by < bl_per_img.height; ++by)
+        {
+            for (int bx = 0; bx < bl_per_img.width; ++bx)
+            {
+                Point bl_tl(bx * bl_width, by * bl_height);
+                Point bl_br(min(bl_tl.x + bl_width, images[img_idx].cols),
+                            min(bl_tl.y + bl_height, images[img_idx].rows));
+
+                block_corners.push_back(corners[img_idx] + bl_tl);
+                block_images.push_back(images[img_idx](Rect(bl_tl, bl_br)));
+                block_masks.push_back(make_pair(masks[img_idx].first(Rect(bl_tl, bl_br)), 
+                                                masks[img_idx].second));
+            }
+        }
+    }
+
+    GainCompensator compensator;
+    compensator.feed(block_corners, block_images, block_masks);
+    vector<double> gains = compensator.gains();
+    gain_maps_.resize(num_images);
+
+    Mat_<float> ker(1, 3); 
+    ker(0,0) = 0.25; ker(0,1) = 0.5; ker(0,2) = 0.25;
+
+    int bl_idx = 0;
+    for (int img_idx = 0; img_idx < num_images; ++img_idx)
+    {
+        Size bl_per_img = bl_per_imgs[img_idx];
+        gain_maps_[img_idx].create(bl_per_img);
+
+        for (int by = 0; by < bl_per_img.height; ++by)
+            for (int bx = 0; bx < bl_per_img.width; ++bx, ++bl_idx)
+                gain_maps_[img_idx](by, bx) = static_cast<float>(gains[bl_idx]);
+        
+        sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker);
+        sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker);
+    }
+}
+
+
+void BlocksGainCompensator::apply(int index, Point /*corner*/, Mat &image, const Mat &/*mask*/)
+{
+    CV_Assert(image.type() == CV_8UC3);
+
+    Mat_<float> gain_map;
+    if (gain_maps_[index].size() == image.size())
+        gain_map = gain_maps_[index];
+    else
+        resize(gain_maps_[index], gain_map, image.size(), 0, 0, INTER_LINEAR);
+
+    for (int y = 0; y < image.rows; ++y)
+    {
+        const float* gain_row = gain_map.ptr<float>(y);
+        Point3_<uchar>* row = image.ptr<Point3_<uchar> >(y);
+        for (int x = 0; x < image.cols; ++x)
+        {
+            row[x].x = saturate_cast<uchar>(row[x].x * gain_row[x]);
+            row[x].y = saturate_cast<uchar>(row[x].y * gain_row[x]);
+            row[x].z = saturate_cast<uchar>(row[x].z * gain_row[x]);
+        }
+    }
+}
+
+} // namespace cv

modules/stitching/src/matchers.cpp

@@ -0,0 +1,510 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+using namespace cv::gpu;
+
+namespace cv
+{
+
+void FeaturesFinder::operator ()(const Mat &image, ImageFeatures &features) 
+{ 
+    find(image, features);
+    features.img_size = image.size();
+    //features.img = image.clone();
+}
+
+} // namespace cv
+
+
+namespace
+{
+    using namespace cv;
+
+    class CpuSurfFeaturesFinder : public FeaturesFinder
+    {
+    public:
+        CpuSurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers, 
+                              int num_octaves_descr, int num_layers_descr) 
+        {
+            detector_ = new SurfFeatureDetector(hess_thresh, num_octaves, num_layers);
+            extractor_ = new SurfDescriptorExtractor(num_octaves_descr, num_layers_descr);
+        }
+
+    protected:
+        void find(const Mat &image, ImageFeatures &features);
+
+    private:
+        Ptr<FeatureDetector> detector_;
+        Ptr<DescriptorExtractor> extractor_;
+    };
+
+
+    class GpuSurfFeaturesFinder : public FeaturesFinder
+    {
+    public:
+        GpuSurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers, 
+                              int num_octaves_descr, int num_layers_descr) 
+        {
+            surf_.keypointsRatio = 0.1f;
+            surf_.hessianThreshold = hess_thresh;
+            surf_.extended = false;
+            num_octaves_ = num_octaves;
+            num_layers_ = num_layers;
+            num_octaves_descr_ = num_octaves_descr;
+            num_layers_descr_ = num_layers_descr;
+        }
+
+        void releaseMemory();
+
+    protected:
+        void find(const Mat &image, ImageFeatures &features);
+
+    private:
+        GpuMat image_;
+        GpuMat gray_image_;
+        SURF_GPU surf_;
+        GpuMat keypoints_;
+        GpuMat descriptors_;
+        int num_octaves_, num_layers_;
+        int num_octaves_descr_, num_layers_descr_;
+    };
+
+
+    void CpuSurfFeaturesFinder::find(const Mat &image, ImageFeatures &features)
+    {
+        Mat gray_image;
+        CV_Assert(image.depth() == CV_8U);
+        cvtColor(image, gray_image, CV_BGR2GRAY);
+        detector_->detect(gray_image, features.keypoints);
+        extractor_->compute(gray_image, features.keypoints, features.descriptors);
+    }
+  
+
+    void GpuSurfFeaturesFinder::find(const Mat &image, ImageFeatures &features)
+    {
+        CV_Assert(image.depth() == CV_8U);
+
+        ensureSizeIsEnough(image.size(), image.type(), image_);
+        image_.upload(image);
+
+        ensureSizeIsEnough(image.size(), CV_8UC1, gray_image_);
+        cvtColor(image_, gray_image_, CV_BGR2GRAY);
+
+        surf_.nOctaves = num_octaves_;
+        surf_.nOctaveLayers = num_layers_;
+        surf_(gray_image_, GpuMat(), keypoints_);
+
+        surf_.nOctaves = num_octaves_descr_;
+        surf_.nOctaveLayers = num_layers_descr_;
+        surf_.upright = true;
+        surf_(gray_image_, GpuMat(), keypoints_, descriptors_, true);
+        surf_.downloadKeypoints(keypoints_, features.keypoints);
+
+        descriptors_.download(features.descriptors);
+    }
+
+    void GpuSurfFeaturesFinder::releaseMemory()
+    {
+        surf_.releaseMemory();
+        image_.release();
+        gray_image_.release();
+        keypoints_.release();
+        descriptors_.release();
+    }
+} // namespace
+
+
+namespace cv
+{
+
+SurfFeaturesFinder::SurfFeaturesFinder(bool try_use_gpu, double hess_thresh, int num_octaves, int num_layers, 
+                                       int num_octaves_descr, int num_layers_descr)
+{
+    if (try_use_gpu && getCudaEnabledDeviceCount() > 0)
+        impl_ = new GpuSurfFeaturesFinder(hess_thresh, num_octaves, num_layers, num_octaves_descr, num_layers_descr);
+    else
+        impl_ = new CpuSurfFeaturesFinder(hess_thresh, num_octaves, num_layers, num_octaves_descr, num_layers_descr);
+}
+
+
+void SurfFeaturesFinder::find(const Mat &image, ImageFeatures &features)
+{
+    (*impl_)(image, features);
+}
+
+void SurfFeaturesFinder::releaseMemory()
+{
+    impl_->releaseMemory();
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+MatchesInfo::MatchesInfo() : src_img_idx(-1), dst_img_idx(-1), num_inliers(0), confidence(0) {}
+
+MatchesInfo::MatchesInfo(const MatchesInfo &other) { *this = other; }
+
+const MatchesInfo& MatchesInfo::operator =(const MatchesInfo &other)
+{
+    src_img_idx = other.src_img_idx;
+    dst_img_idx = other.dst_img_idx;
+    matches = other.matches;
+    inliers_mask = other.inliers_mask;
+    num_inliers = other.num_inliers;
+    H = other.H.clone();
+    confidence = other.confidence;
+    return *this;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+struct DistIdxPair
+{
+    bool operator<(const DistIdxPair &other) const { return dist < other.dist; }
+    double dist;
+    int idx;
+};
+
+
+struct MatchPairsBody
+{
+    MatchPairsBody(const MatchPairsBody& other)
+            : matcher(other.matcher), features(other.features), 
+              pairwise_matches(other.pairwise_matches), near_pairs(other.near_pairs) {}
+
+    MatchPairsBody(FeaturesMatcher &matcher, const vector<ImageFeatures> &features, 
+                   vector<MatchesInfo> &pairwise_matches, vector<pair<int,int> > &near_pairs)
+            : matcher(matcher), features(features), 
+              pairwise_matches(pairwise_matches), near_pairs(near_pairs) {}
+
+    void operator ()(const BlockedRange &r) const 
+    {
+        const int num_images = static_cast<int>(features.size());
+        for (int i = r.begin(); i < r.end(); ++i)
+        {
+            int from = near_pairs[i].first;
+            int to = near_pairs[i].second;
+            int pair_idx = from*num_images + to;
+
+            matcher(features[from], features[to], pairwise_matches[pair_idx]);
+            pairwise_matches[pair_idx].src_img_idx = from;
+            pairwise_matches[pair_idx].dst_img_idx = to;
+
+            size_t dual_pair_idx = to*num_images + from;
+
+            pairwise_matches[dual_pair_idx] = pairwise_matches[pair_idx];
+            pairwise_matches[dual_pair_idx].src_img_idx = to;
+            pairwise_matches[dual_pair_idx].dst_img_idx = from;
+
+            if (!pairwise_matches[pair_idx].H.empty())
+                pairwise_matches[dual_pair_idx].H = pairwise_matches[pair_idx].H.inv();
+
+            for (size_t j = 0; j < pairwise_matches[dual_pair_idx].matches.size(); ++j)
+                std::swap(pairwise_matches[dual_pair_idx].matches[j].queryIdx,
+                          pairwise_matches[dual_pair_idx].matches[j].trainIdx);
+            LOG(".");
+        }
+    }
+
+    FeaturesMatcher &matcher;
+    const vector<ImageFeatures> &features;
+    vector<MatchesInfo> &pairwise_matches;
+    vector<pair<int,int> > &near_pairs;
+
+private:
+    void operator =(const MatchPairsBody&);
+};
+
+
+void FeaturesMatcher::operator ()(const vector<ImageFeatures> &features, vector<MatchesInfo> &pairwise_matches)
+{
+    const int num_images = static_cast<int>(features.size());
+
+    vector<pair<int,int> > near_pairs;
+    for (int i = 0; i < num_images - 1; ++i)
+        for (int j = i + 1; j < num_images; ++j)
+            near_pairs.push_back(make_pair(i, j));
+
+    pairwise_matches.resize(num_images * num_images);
+    MatchPairsBody body(*this, features, pairwise_matches, near_pairs);
+
+    if (is_thread_safe_)
+        parallel_for(BlockedRange(0, static_cast<int>(near_pairs.size())), body);
+    else
+        body(BlockedRange(0, static_cast<int>(near_pairs.size())));
+    LOGLN("");
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+namespace 
+{
+    typedef set<pair<int,int> > MatchesSet;
+
+    // These two classes are aimed to find features matches only, not to 
+    // estimate homography
+
+    class CpuMatcher : public FeaturesMatcher
+    {
+    public:
+        CpuMatcher(float match_conf) : FeaturesMatcher(true), match_conf_(match_conf) {}
+        void match(const ImageFeatures &features1, const ImageFeatures &features2, MatchesInfo& matches_info);
+
+    private:
+        float match_conf_;
+    };
+
+
+    class GpuMatcher : public FeaturesMatcher
+    {
+    public:
+        GpuMatcher(float match_conf) : match_conf_(match_conf) {}
+        void match(const ImageFeatures &features1, const ImageFeatures &features2, MatchesInfo& matches_info);
+
+        void releaseMemory();
+
+    private:
+        float match_conf_;
+        GpuMat descriptors1_, descriptors2_;
+        GpuMat train_idx_, distance_, all_dist_;
+        vector< vector<DMatch> > pair_matches;
+    };
+
+
+    void CpuMatcher::match(const ImageFeatures &features1, const ImageFeatures &features2, MatchesInfo& matches_info)
+    {
+        matches_info.matches.clear();
+        FlannBasedMatcher matcher;
+        vector< vector<DMatch> > pair_matches;        
+        MatchesSet matches;
+
+        // Find 1->2 matches
+        matcher.knnMatch(features1.descriptors, features2.descriptors, pair_matches, 2);
+        for (size_t i = 0; i < pair_matches.size(); ++i)
+        {
+            if (pair_matches[i].size() < 2)
+                continue;
+            const DMatch& m0 = pair_matches[i][0];
+            const DMatch& m1 = pair_matches[i][1];
+            if (m0.distance < (1.f - match_conf_) * m1.distance)
+            {
+                matches_info.matches.push_back(m0);
+                matches.insert(make_pair(m0.queryIdx, m0.trainIdx));
+            }
+        }
+
+        // Find 2->1 matches
+        pair_matches.clear();
+        matcher.knnMatch(features2.descriptors, features1.descriptors, pair_matches, 2);
+        for (size_t i = 0; i < pair_matches.size(); ++i)
+        {
+            if (pair_matches[i].size() < 2)
+                continue;
+            const DMatch& m0 = pair_matches[i][0];
+            const DMatch& m1 = pair_matches[i][1];
+            if (m0.distance < (1.f - match_conf_) * m1.distance)
+                if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) == matches.end())
+                    matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance));
+        }
+    }
+       
+
+    void GpuMatcher::match(const ImageFeatures &features1, const ImageFeatures &features2, MatchesInfo& matches_info)
+    {
+        matches_info.matches.clear(); 
+
+        ensureSizeIsEnough(features1.descriptors.size(), features1.descriptors.type(), descriptors1_);
+        ensureSizeIsEnough(features2.descriptors.size(), features2.descriptors.type(), descriptors2_);
+
+        descriptors1_.upload(features1.descriptors);
+        descriptors2_.upload(features2.descriptors);
+
+        BruteForceMatcher_GPU< L2<float> > matcher;
+        MatchesSet matches;
+
+        // Find 1->2 matches
+        pair_matches.clear();
+        matcher.knnMatch(descriptors1_, descriptors2_, train_idx_, distance_, all_dist_, 2);
+        matcher.knnMatchDownload(train_idx_, distance_, pair_matches);
+        for (size_t i = 0; i < pair_matches.size(); ++i)
+        {
+            if (pair_matches[i].size() < 2)
+                continue;
+            const DMatch& m0 = pair_matches[i][0];
+            const DMatch& m1 = pair_matches[i][1];
+            if (m0.distance < (1.f - match_conf_) * m1.distance)
+            {
+                matches_info.matches.push_back(m0);
+                matches.insert(make_pair(m0.queryIdx, m0.trainIdx));
+            }
+        }
+
+        // Find 2->1 matches
+        pair_matches.clear();
+        matcher.knnMatch(descriptors2_, descriptors1_, train_idx_, distance_, all_dist_, 2);
+        matcher.knnMatchDownload(train_idx_, distance_, pair_matches);
+        for (size_t i = 0; i < pair_matches.size(); ++i)
+        {
+            if (pair_matches[i].size() < 2)
+                continue;
+            const DMatch& m0 = pair_matches[i][0];
+            const DMatch& m1 = pair_matches[i][1];
+            if (m0.distance < (1.f - match_conf_) * m1.distance)
+                if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) == matches.end())
+                    matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance));
+        }
+    }
+
+    void GpuMatcher::releaseMemory()
+    {
+        descriptors1_.release();
+        descriptors2_.release();
+        train_idx_.release();
+        distance_.release();
+        all_dist_.release();
+        vector< vector<DMatch> >().swap(pair_matches);
+    }
+
+} // anonymous namespace
+
+
+BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu, float match_conf, int num_matches_thresh1, int num_matches_thresh2)
+{
+    if (try_use_gpu && getCudaEnabledDeviceCount() > 0)
+        impl_ = new GpuMatcher(match_conf);
+    else
+        impl_ = new CpuMatcher(match_conf);
+
+    is_thread_safe_ = impl_->isThreadSafe();
+    num_matches_thresh1_ = num_matches_thresh1;
+    num_matches_thresh2_ = num_matches_thresh2;
+}
+
+
+void BestOf2NearestMatcher::match(const ImageFeatures &features1, const ImageFeatures &features2,
+                                  MatchesInfo &matches_info)
+{
+    (*impl_)(features1, features2, matches_info);
+
+    //Mat out;
+    //drawMatches(features1.img, features1.keypoints, features2.img, features2.keypoints, matches_info.matches, out);
+    //stringstream ss;
+    //ss << features1.img_idx << features2.img_idx << ".png";
+    //imwrite(ss.str(), out);
+
+    // Check if it makes sense to find homography
+    if (matches_info.matches.size() < static_cast<size_t>(num_matches_thresh1_))
+        return;
+
+    // Construct point-point correspondences for homography estimation
+    Mat src_points(1, static_cast<int>(matches_info.matches.size()), CV_32FC2);
+    Mat dst_points(1, static_cast<int>(matches_info.matches.size()), CV_32FC2);
+    for (size_t i = 0; i < matches_info.matches.size(); ++i)
+    {
+        const DMatch& m = matches_info.matches[i];
+
+        Point2f p = features1.keypoints[m.queryIdx].pt;
+        p.x -= features1.img_size.width * 0.5f;
+        p.y -= features1.img_size.height * 0.5f;
+        src_points.at<Point2f>(0, static_cast<int>(i)) = p;
+
+        p = features2.keypoints[m.trainIdx].pt;
+        p.x -= features2.img_size.width * 0.5f;
+        p.y -= features2.img_size.height * 0.5f;
+        dst_points.at<Point2f>(0, static_cast<int>(i)) = p;
+    }
+
+    // Find pair-wise motion
+    matches_info.H = findHomography(src_points, dst_points, matches_info.inliers_mask, CV_RANSAC);
+    if (std::abs(determinant(matches_info.H)) < numeric_limits<double>::epsilon())
+        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++;
+
+    matches_info.confidence = matches_info.num_inliers / (8 + 0.3*matches_info.matches.size());
+
+    // Check if we should try to refine motion
+    if (matches_info.num_inliers < num_matches_thresh2_)
+        return;
+
+    // Construct point-point correspondences for inliers only
+    src_points.create(1, matches_info.num_inliers, CV_32FC2);
+    dst_points.create(1, matches_info.num_inliers, CV_32FC2);
+    int inlier_idx = 0;
+    for (size_t i = 0; i < matches_info.matches.size(); ++i)
+    {
+        if (!matches_info.inliers_mask[i])
+            continue;
+
+        const DMatch& m = matches_info.matches[i];
+
+        Point2f p = features1.keypoints[m.queryIdx].pt;
+        p.x -= features1.img_size.width * 0.5f;
+        p.y -= features1.img_size.height * 0.5f;
+        src_points.at<Point2f>(0, inlier_idx) = p;
+
+        p = features2.keypoints[m.trainIdx].pt;
+        p.x -= features2.img_size.width * 0.5f;
+        p.y -= features2.img_size.height * 0.5f;
+        dst_points.at<Point2f>(0, inlier_idx) = p;
+
+        inlier_idx++;
+    }
+
+    // Rerun motion estimation on inliers only
+    matches_info.H = findHomography(src_points, dst_points, CV_RANSAC);
+}
+
+void BestOf2NearestMatcher::releaseMemory()
+{
+    impl_->releaseMemory();
+}
+
+} // namespace cv

modules/stitching/src/motion_estimators.cpp

@@ -0,0 +1,604 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+struct IncDistance
+{
+    IncDistance(vector<int> &dists) : dists(&dists[0]) {}
+    void operator ()(const GraphEdge &edge) { dists[edge.to] = dists[edge.from] + 1; }
+    int* dists;
+};
+
+
+struct CalcRotation
+{
+    CalcRotation(int num_images, const vector<MatchesInfo> &pairwise_matches, vector<CameraParams> &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;
+
+        double f_from = cameras[edge.from].focal;
+        double f_to = cameras[edge.to].focal;
+
+        Mat K_from = Mat::eye(3, 3, CV_64F);
+        K_from.at<double>(0, 0) = f_from;
+        K_from.at<double>(1, 1) = f_from;
+
+        Mat K_to = Mat::eye(3, 3, CV_64F);
+        K_to.at<double>(0, 0) = f_to;
+        K_to.at<double>(1, 1) = f_to;
+
+        Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to;
+        cameras[edge.to].R = cameras[edge.from].R * R;
+    }
+
+    int num_images;
+    const MatchesInfo* pairwise_matches;
+    CameraParams* cameras;
+};
+
+
+void HomographyBasedEstimator::estimate(const vector<ImageFeatures> &features, const vector<MatchesInfo> &pairwise_matches, 
+                                        vector<CameraParams> &cameras)
+{
+    LOGLN("Estimating rotations...");
+    int64 t = getTickCount();
+
+    const int num_images = static_cast<int>(features.size());
+
+#if 0
+    // Robustly estimate focal length from rotating cameras
+    vector<Mat> Hs;
+    for (int iter = 0; iter < 100; ++iter)
+    {
+        int len = 2 + rand()%(pairwise_matches.size() - 1);
+        vector<int> subset;
+        selectRandomSubset(len, pairwise_matches.size(), subset);
+        Hs.clear();
+        for (size_t i = 0; i < subset.size(); ++i)
+            if (!pairwise_matches[subset[i]].H.empty())
+                Hs.push_back(pairwise_matches[subset[i]].H);
+        Mat_<double> K;
+        if (Hs.size() >= 2)
+        {
+            if (calibrateRotatingCamera(Hs, K))
+                cin.get();
+        }
+    }
+#endif
+
+    // Estimate focal length and set it for all cameras
+    vector<double> focals;
+    estimateFocal(features, pairwise_matches, focals);
+    cameras.resize(num_images);
+    for (int i = 0; i < num_images; ++i)
+        cameras[i].focal = focals[i];
+
+    // Restore global motion
+    Graph span_tree;
+    vector<int> span_tree_centers;
+    findMaxSpanningTree(num_images, pairwise_matches, span_tree, span_tree_centers);
+    span_tree.walkBreadthFirst(span_tree_centers[0], CalcRotation(num_images, pairwise_matches, cameras));
+
+    LOGLN("Estimating rotations, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+void BundleAdjuster::estimate(const vector<ImageFeatures> &features, const vector<MatchesInfo> &pairwise_matches, 
+                              vector<CameraParams> &cameras)
+{
+    if (cost_space_ == NO)
+        return;
+
+    LOG("Bundle adjustment");
+    int64 t = getTickCount();
+
+    num_images_ = static_cast<int>(features.size());
+    features_ = &features[0];
+    pairwise_matches_ = &pairwise_matches[0];
+
+    // Prepare focals and rotations
+    cameras_.create(num_images_ * 4, 1, CV_64F);
+    SVD svd;
+    for (int i = 0; i < num_images_; ++i)
+    {
+        cameras_.at<double>(i * 4, 0) = cameras[i].focal;
+
+        svd(cameras[i].R, SVD::FULL_UV);
+        Mat R = svd.u * svd.vt;
+        if (determinant(R) < 0) 
+            R *= -1;
+
+        Mat rvec;
+        Rodrigues(R, rvec); CV_Assert(rvec.type() == CV_32F);
+        cameras_.at<double>(i * 4 + 1, 0) = rvec.at<float>(0, 0);
+        cameras_.at<double>(i * 4 + 2, 0) = rvec.at<float>(1, 0);
+        cameras_.at<double>(i * 4 + 3, 0) = rvec.at<float>(2, 0);
+    }
+
+    // Select only consistent image pairs for futher adjustment
+    edges_.clear();
+    for (int i = 0; i < num_images_ - 1; ++i)
+    {
+        for (int j = i + 1; j < num_images_; ++j)
+        {
+            const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j];
+            if (matches_info.confidence > conf_thresh_)
+                edges_.push_back(make_pair(i, j));
+        }
+    }
+
+    // Compute number of correspondences
+    total_num_matches_ = 0;
+    for (size_t i = 0; i < edges_.size(); ++i)
+        total_num_matches_ += static_cast<int>(pairwise_matches[edges_[i].first * num_images_ + edges_[i].second].num_inliers);
+
+    CvLevMarq solver(num_images_ * 4, total_num_matches_ * 3,
+                     cvTermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 1000, DBL_EPSILON));
+
+    CvMat matParams = cameras_;
+    cvCopy(&matParams, solver.param);
+
+    int count = 0;
+    for(;;)
+    {
+        const CvMat* _param = 0;
+        CvMat* _J = 0;
+        CvMat* _err = 0;
+
+        bool proceed = solver.update(_param, _J, _err);
+
+        cvCopy( _param, &matParams );
+
+        if( !proceed || !_err )
+            break;
+
+        if( _J )
+        {
+            calcJacobian();
+            CvMat matJ = J_;
+            cvCopy( &matJ, _J );
+        }
+
+        if (_err)
+        {
+            calcError(err_);
+            LOG(".");
+            count++;
+            CvMat matErr = err_;
+            cvCopy( &matErr, _err );
+        }
+    }
+    LOGLN("");
+    LOGLN("Bundle adjustment, final error: " << sqrt(err_.dot(err_)));
+    LOGLN("Bundle adjustment, iterations done: " << count);
+
+    // Obtain global motion
+    for (int i = 0; i < num_images_; ++i)
+    {
+        cameras[i].focal = cameras_.at<double>(i * 4, 0);
+        Mat rvec(3, 1, CV_64F);
+        rvec.at<double>(0, 0) = cameras_.at<double>(i * 4 + 1, 0);
+        rvec.at<double>(1, 0) = cameras_.at<double>(i * 4 + 2, 0);
+        rvec.at<double>(2, 0) = cameras_.at<double>(i * 4 + 3, 0);
+        Rodrigues(rvec, cameras[i].R);
+        Mat Mf;
+        cameras[i].R.convertTo(Mf, CV_32F);
+        cameras[i].R = Mf;
+    }
+
+    // Normalize motion to center image
+    Graph span_tree;
+    vector<int> span_tree_centers;
+    findMaxSpanningTree(num_images_, pairwise_matches, span_tree, span_tree_centers);
+    Mat R_inv = cameras[span_tree_centers[0]].R.inv();
+    for (int i = 0; i < num_images_; ++i)
+        cameras[i].R = R_inv * cameras[i].R;
+
+    LOGLN("Bundle adjustment, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
+}
+
+
+void BundleAdjuster::calcError(Mat &err)
+{
+    err.create(total_num_matches_ * 3, 1, CV_64F);
+
+    int match_idx = 0;
+    for (size_t edge_idx = 0; edge_idx < edges_.size(); ++edge_idx)
+    {
+        int i = edges_[edge_idx].first;
+        int j = edges_[edge_idx].second;
+        double f1 = cameras_.at<double>(i * 4, 0);
+        double f2 = cameras_.at<double>(j * 4, 0);
+        double R1[9], R2[9];
+        Mat R1_(3, 3, CV_64F, R1), R2_(3, 3, CV_64F, R2);
+        Mat rvec(3, 1, CV_64F);
+        rvec.at<double>(0, 0) = cameras_.at<double>(i * 4 + 1, 0);
+        rvec.at<double>(1, 0) = cameras_.at<double>(i * 4 + 2, 0);
+        rvec.at<double>(2, 0) = cameras_.at<double>(i * 4 + 3, 0);
+        Rodrigues(rvec, R1_); CV_Assert(R1_.type() == CV_64F);
+        rvec.at<double>(0, 0) = cameras_.at<double>(j * 4 + 1, 0);
+        rvec.at<double>(1, 0) = cameras_.at<double>(j * 4 + 2, 0);
+        rvec.at<double>(2, 0) = cameras_.at<double>(j * 4 + 3, 0);
+        Rodrigues(rvec, R2_); CV_Assert(R2_.type() == CV_64F);
+
+        const ImageFeatures& features1 = features_[i];
+        const ImageFeatures& features2 = features_[j];
+        const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j];
+
+        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];
+
+            Point2d kp1 = features1.keypoints[m.queryIdx].pt;
+            kp1.x -= 0.5 * features1.img_size.width;
+            kp1.y -= 0.5 * features1.img_size.height;
+            Point2d kp2 = features2.keypoints[m.trainIdx].pt;
+            kp2.x -= 0.5 * features2.img_size.width;
+            kp2.y -= 0.5 * features2.img_size.height;
+            double len1 = sqrt(kp1.x * kp1.x + kp1.y * kp1.y + f1 * f1);
+            double len2 = sqrt(kp2.x * kp2.x + kp2.y * kp2.y + f2 * f2);
+            Point3d p1(kp1.x / len1, kp1.y / len1, f1 / len1);
+            Point3d p2(kp2.x / len2, kp2.y / len2, f2 / len2);
+
+            Point3d d1(p1.x * R1[0] + p1.y * R1[1] + p1.z * R1[2],
+                       p1.x * R1[3] + p1.y * R1[4] + p1.z * R1[5],
+                       p1.x * R1[6] + p1.y * R1[7] + p1.z * R1[8]);
+            Point3d d2(p2.x * R2[0] + p2.y * R2[1] + p2.z * R2[2],
+                       p2.x * R2[3] + p2.y * R2[4] + p2.z * R2[5],
+                       p2.x * R2[6] + p2.y * R2[7] + p2.z * R2[8]);
+
+            double mult = 1; // For cost_space_ == RAY_SPACE
+            if (cost_space_ == FOCAL_RAY_SPACE)
+                mult = sqrt(f1 * f2);
+            err.at<double>(3 * match_idx, 0) = mult * (d1.x - d2.x);
+            err.at<double>(3 * match_idx + 1, 0) = mult * (d1.y - d2.y);
+            err.at<double>(3 * match_idx + 2, 0) = mult * (d1.z - d2.z);
+            match_idx++;
+        }
+    }
+}
+
+
+void calcDeriv(const Mat &err1, const Mat &err2, double h, Mat res)
+{
+    for (int i = 0; i < err1.rows; ++i)
+        res.at<double>(i, 0) = (err2.at<double>(i, 0) - err1.at<double>(i, 0)) / h;
+}
+
+
+void BundleAdjuster::calcJacobian()
+{
+    J_.create(total_num_matches_ * 3, num_images_ * 4, CV_64F);
+
+    double f, r;
+    const double df = 0.001; // Focal length step
+    const double dr = 0.001; // Angle step
+
+    for (int i = 0; i < num_images_; ++i)
+    {
+        f = cameras_.at<double>(i * 4, 0);
+        cameras_.at<double>(i * 4, 0) = f - df;
+        calcError(err1_);
+        cameras_.at<double>(i * 4, 0) = f + df;
+        calcError(err2_);
+        calcDeriv(err1_, err2_, 2 * df, J_.col(i * 4));
+        cameras_.at<double>(i * 4, 0) = f;
+
+        r = cameras_.at<double>(i * 4 + 1, 0);
+        cameras_.at<double>(i * 4 + 1, 0) = r - dr;
+        calcError(err1_);
+        cameras_.at<double>(i * 4 + 1, 0) = r + dr;
+        calcError(err2_);
+        calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 1));
+        cameras_.at<double>(i * 4 + 1, 0) = r;
+
+        r = cameras_.at<double>(i * 4 + 2, 0);
+        cameras_.at<double>(i * 4 + 2, 0) = r - dr;
+        calcError(err1_);
+        cameras_.at<double>(i * 4 + 2, 0) = r + dr;
+        calcError(err2_);
+        calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 2));
+        cameras_.at<double>(i * 4 + 2, 0) = r;
+
+        r = cameras_.at<double>(i * 4 + 3, 0);
+        cameras_.at<double>(i * 4 + 3, 0) = r - dr;
+        calcError(err1_);
+        cameras_.at<double>(i * 4 + 3, 0) = r + dr;
+        calcError(err2_);
+        calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 3));
+        cameras_.at<double>(i * 4 + 3, 0) = r;
+    }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+// TODO replace SVD with eigen
+void waveCorrect(vector<Mat> &rmats)
+{
+    LOGLN("Wave correcting...");
+    int64 t = getTickCount();
+
+    float data[9];
+    Mat r0(1, 3, CV_32F, data);
+    Mat r1(1, 3, CV_32F, data + 3);
+    Mat r2(1, 3, CV_32F, data + 6);
+    Mat R(3, 3, CV_32F, data);
+
+    Mat cov = Mat::zeros(3, 3, CV_32F);
+    for (size_t i = 0; i < rmats.size(); ++i)
+    {   
+        Mat r0 = rmats[i].col(0);
+        cov += r0 * r0.t();
+    }
+
+    SVD svd;
+    svd(cov, SVD::FULL_UV);
+    svd.vt.row(2).copyTo(r1);
+    if (determinant(svd.vt) < 0) r1 *= -1;
+
+    Mat avgz = Mat::zeros(3, 1, CV_32F);
+    for (size_t i = 0; i < rmats.size(); ++i)
+        avgz += rmats[i].col(2);
+    r1.cross(avgz.t()).copyTo(r0);
+    normalize(r0, r0);
+
+    r1.cross(r0).copyTo(r2);
+    if (determinant(R) < 0) R *= -1;
+
+    for (size_t i = 0; i < rmats.size(); ++i)
+        rmats[i] = R * rmats[i];
+
+    LOGLN("Wave correcting, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+string matchesGraphAsString(vector<string> &pathes, vector<MatchesInfo> &pairwise_matches,
+                            float conf_threshold)
+{
+    stringstream str;
+    str << "graph matches_graph{\n";
+
+    const int num_images = static_cast<int>(pathes.size());
+    set<pair<int,int> > span_tree_edges;
+    DisjointSets comps(num_images);
+
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = 0; j < num_images; ++j)
+        {
+            if (pairwise_matches[i*num_images + j].confidence < conf_threshold)
+                continue;
+            int comp1 = comps.findSetByElem(i);
+            int comp2 = comps.findSetByElem(j);
+            if (comp1 != comp2)
+            {
+                comps.mergeSets(comp1, comp2);
+                span_tree_edges.insert(make_pair(i, j));
+            }
+        }
+    }
+
+    for (set<pair<int,int> >::const_iterator itr = span_tree_edges.begin();
+         itr != span_tree_edges.end(); ++itr)
+    {
+        pair<int,int> edge = *itr;
+        if (span_tree_edges.find(edge) != span_tree_edges.end())
+        {
+            string name_src = pathes[edge.first];
+            size_t prefix_len = name_src.find_last_of("/\\");
+            if (prefix_len != string::npos) prefix_len++; else prefix_len = 0;
+            name_src = name_src.substr(prefix_len, name_src.size() - prefix_len);
+
+            string name_dst = pathes[edge.second];
+            prefix_len = name_dst.find_last_of("/\\");
+            if (prefix_len != string::npos) prefix_len++; else prefix_len = 0;
+            name_dst = name_dst.substr(prefix_len, name_dst.size() - prefix_len);
+
+            int pos = edge.first*num_images + edge.second;
+            str << "\"" << name_src << "\" -- \"" << name_dst << "\""
+                << "[label=\"Nm=" << pairwise_matches[pos].matches.size()
+                << ", Ni=" << pairwise_matches[pos].num_inliers
+                << ", C=" << pairwise_matches[pos].confidence << "\"];\n";
+        }
+    }
+
+    for (size_t i = 0; i < comps.size.size(); ++i)
+    {
+        if (comps.size[comps.findSetByElem(i)] == 1)
+        {
+            string name = pathes[i];
+            size_t prefix_len = name.find_last_of("/\\");
+            if (prefix_len != string::npos) prefix_len++; else prefix_len = 0;
+            name = name.substr(prefix_len, name.size() - prefix_len);
+            str << "\"" << name << "\";\n";
+        }
+    }
+
+    str << "}";
+    return str.str();
+}
+
+vector<int> leaveBiggestComponent(vector<ImageFeatures> &features,  vector<MatchesInfo> &pairwise_matches, 
+                                  float conf_threshold)
+{
+    const int num_images = static_cast<int>(features.size());
+
+    DisjointSets comps(num_images);
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = 0; j < num_images; ++j)
+        {
+            if (pairwise_matches[i*num_images + j].confidence < conf_threshold)
+                continue;
+            int comp1 = comps.findSetByElem(i);
+            int comp2 = comps.findSetByElem(j);
+            if (comp1 != comp2) 
+                comps.mergeSets(comp1, comp2);
+        }
+    }
+
+    int max_comp = static_cast<int>(max_element(comps.size.begin(), comps.size.end()) - comps.size.begin());
+
+    vector<int> indices;
+    vector<int> indices_removed;
+    for (int i = 0; i < num_images; ++i)
+        if (comps.findSetByElem(i) == max_comp)
+            indices.push_back(i);    
+        else
+            indices_removed.push_back(i);
+
+    vector<ImageFeatures> features_subset;
+    vector<MatchesInfo> pairwise_matches_subset;
+    for (size_t i = 0; i < indices.size(); ++i)
+    {
+        features_subset.push_back(features[indices[i]]);
+        for (size_t j = 0; j < indices.size(); ++j)
+        {
+            pairwise_matches_subset.push_back(pairwise_matches[indices[i]*num_images + indices[j]]);
+            pairwise_matches_subset.back().src_img_idx = static_cast<int>(i);
+            pairwise_matches_subset.back().dst_img_idx = static_cast<int>(j);
+        }
+    }
+
+    if (static_cast<int>(features_subset.size()) == num_images)
+        return indices;
+
+    LOG("Removed some images, because can't match them: (");
+    LOG(indices_removed[0]+1);
+    for (size_t i = 1; i < indices_removed.size(); ++i) 
+        LOG(", " << indices_removed[i]+1);
+    LOGLN("). Try to decrease --match_conf value.");
+
+    features = features_subset;
+    pairwise_matches = pairwise_matches_subset;
+
+    return indices;
+}
+
+
+void findMaxSpanningTree(int num_images, const vector<MatchesInfo> &pairwise_matches,
+                         Graph &span_tree, vector<int> &centers)
+{
+    Graph graph(num_images);
+    vector<GraphEdge> edges;
+
+    // Construct images graph and remember its edges
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = 0; j < num_images; ++j)
+        {
+            if (pairwise_matches[i * num_images + j].H.empty())
+                continue;
+            float conf = static_cast<float>(pairwise_matches[i * num_images + j].num_inliers);
+            graph.addEdge(i, j, conf);
+            edges.push_back(GraphEdge(i, j, conf));
+        }
+    }
+
+    DisjointSets comps(num_images);
+    span_tree.create(num_images);
+    vector<int> span_tree_powers(num_images, 0);
+
+    // Find maximum spanning tree
+    sort(edges.begin(), edges.end(), greater<GraphEdge>());
+    for (size_t i = 0; i < edges.size(); ++i)
+    {
+        int comp1 = comps.findSetByElem(edges[i].from);
+        int comp2 = comps.findSetByElem(edges[i].to);
+        if (comp1 != comp2)
+        {
+            comps.mergeSets(comp1, comp2);
+            span_tree.addEdge(edges[i].from, edges[i].to, edges[i].weight);
+            span_tree.addEdge(edges[i].to, edges[i].from, edges[i].weight);
+            span_tree_powers[edges[i].from]++;
+            span_tree_powers[edges[i].to]++;
+        }
+    }
+
+    // Find spanning tree leafs
+    vector<int> span_tree_leafs;
+    for (int i = 0; i < num_images; ++i)
+        if (span_tree_powers[i] == 1)
+            span_tree_leafs.push_back(i);
+
+    // Find maximum distance from each spanning tree vertex
+    vector<int> max_dists(num_images, 0);
+    vector<int> cur_dists;
+    for (size_t i = 0; i < span_tree_leafs.size(); ++i)
+    {
+        cur_dists.assign(num_images, 0);
+        span_tree.walkBreadthFirst(span_tree_leafs[i], IncDistance(cur_dists));
+        for (int j = 0; j < num_images; ++j)
+            max_dists[j] = max(max_dists[j], cur_dists[j]);
+    }
+
+    // Find min-max distance
+    int min_max_dist = max_dists[0];
+    for (int i = 1; i < num_images; ++i)
+        if (min_max_dist > max_dists[i])
+            min_max_dist = max_dists[i];
+
+    // Find spanning tree centers
+    centers.clear();
+    for (int i = 0; i < num_images; ++i)
+        if (max_dists[i] == min_max_dist)
+            centers.push_back(i);
+    CV_Assert(centers.size() > 0 && centers.size() <= 2);
+}
+
+} // namespace cv

modules/stitching/src/precomp.cpp

@@ -0,0 +1,43 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"

modules/stitching/src/precomp.hpp

@@ -0,0 +1,72 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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*/
+#ifndef __OPENCV_STITCHING_PRECOMP_H__
+#define __OPENCV_STITCHING_PRECOMP_H__
+
+#ifdef HAVE_CVCONFIG_H 
+#include "cvconfig.h"
+#endif
+
+#include <vector>
+#include <algorithm>
+#include <utility>
+#include <set>
+#include <functional>
+#include <sstream>
+#include "opencv2/stitching/autocalib.hpp"
+#include "opencv2/stitching/blenders.hpp"
+#include "opencv2/stitching/camera.hpp"
+#include "opencv2/stitching/exposure_compensate.hpp"
+#include "opencv2/stitching/matchers.hpp"
+#include "opencv2/stitching/motion_estimators.hpp"
+#include "opencv2/stitching/seam_finders.hpp"
+#include "opencv2/stitching/util.hpp"
+#include "opencv2/stitching/warpers.hpp"
+#include "opencv2/core/core.hpp"
+#include "opencv2/core/internal.hpp"
+#include "opencv2/imgproc/imgproc.hpp"
+#include "opencv2/features2d/features2d.hpp"
+#include "opencv2/calib3d/calib3d.hpp"
+#include "opencv2/gpu/gpu.hpp"
+#include "gcgraph.hpp"
+
+#endif

modules/stitching/src/seam_finders.cpp

@@ -0,0 +1,409 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+Ptr<SeamFinder> SeamFinder::createDefault(int type)
+{
+    if (type == NO)
+        return new NoSeamFinder();
+    if (type == VORONOI)
+        return new VoronoiSeamFinder();
+    if (type == GC_COLOR)
+        return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR);
+    if (type == GC_COLOR_GRAD)
+        return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR_GRAD);
+    CV_Error(CV_StsBadArg, "unsupported seam finding method");
+    return NULL;
+}
+
+
+void PairwiseSeamFinder::find(const vector<Mat> &src, const vector<Point> &corners,
+                              vector<Mat> &masks)
+{
+    LOGLN("Finding seams...");
+    int64 t = getTickCount();
+
+    if (src.size() == 0)
+        return;
+
+    images_ = src;
+    corners_ = corners;
+    masks_ = masks;
+
+    for (size_t i = 0; i < src.size() - 1; ++i)
+    {
+        for (size_t j = i + 1; j < src.size(); ++j)
+        {
+            Rect roi;
+            if (overlapRoi(corners[i], corners[j], src[i].size(), src[j].size(), roi))
+                findInPair(i, j, roi);
+        }
+    }
+
+    LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
+}
+
+
+void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi)
+{
+    const int gap = 10;
+    Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
+    Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
+
+    Mat img1 = images_[first], img2 = images_[second];
+    Mat mask1 = masks_[first], mask2 = masks_[second];
+    Point tl1 = corners_[first], tl2 = corners_[second];
+
+    // Cut submasks with some gap
+    for (int y = -gap; y < roi.height + gap; ++y)
+    {
+        for (int x = -gap; x < roi.width + gap; ++x)
+        {
+            int y1 = roi.y - tl1.y + y;
+            int x1 = roi.x - tl1.x + x;
+            if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols)
+                submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1);
+            else
+                submask1.at<uchar>(y + gap, x + gap) = 0;
+
+            int y2 = roi.y - tl2.y + y;
+            int x2 = roi.x - tl2.x + x;
+            if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols)
+                submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2);
+            else
+                submask2.at<uchar>(y + gap, x + gap) = 0;
+        }
+    }
+
+    Mat collision = (submask1 != 0) & (submask2 != 0);
+    Mat unique1 = submask1.clone(); unique1.setTo(0, collision);
+    Mat unique2 = submask2.clone(); unique2.setTo(0, collision);
+
+    Mat dist1, dist2;
+    distanceTransform(unique1 == 0, dist1, CV_DIST_L1, 3);
+    distanceTransform(unique2 == 0, dist2, CV_DIST_L1, 3);
+
+    Mat seam = dist1 < dist2;
+
+    for (int y = 0; y < roi.height; ++y)
+    {
+        for (int x = 0; x < roi.width; ++x)
+        {
+            if (seam.at<uchar>(y + gap, x + gap))
+                mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0;
+            else
+                mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0;
+        }
+    }
+}
+
+
+class GraphCutSeamFinder::Impl : public PairwiseSeamFinder
+{
+public:
+    Impl(int cost_type, float terminal_cost, float bad_region_penalty)
+        : cost_type_(cost_type), terminal_cost_(terminal_cost), bad_region_penalty_(bad_region_penalty) {}
+
+    void find(const vector<Mat> &src, const vector<Point> &corners, vector<Mat> &masks);
+    void findInPair(size_t first, size_t second, Rect roi);
+
+private:
+    void setGraphWeightsColor(const Mat &img1, const Mat &img2, 
+                              const Mat &mask1, const Mat &mask2, GCGraph<float> &graph);
+    void setGraphWeightsColorGrad(const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, 
+                                  const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, 
+                                  GCGraph<float> &graph);
+
+    vector<Mat> dx_, dy_;
+    int cost_type_;
+    float terminal_cost_;
+    float bad_region_penalty_;
+};
+
+
+void GraphCutSeamFinder::Impl::find(const vector<Mat> &src, const vector<Point> &corners, 
+                                    vector<Mat> &masks)
+{
+    // Compute gradients
+    dx_.resize(src.size());
+    dy_.resize(src.size());
+    Mat dx, dy;
+    for (size_t i = 0; i < src.size(); ++i)
+    {
+        CV_Assert(src[i].channels() == 3);
+        Sobel(src[i], dx, CV_32F, 1, 0);
+        Sobel(src[i], dy, CV_32F, 0, 1);
+        dx_[i].create(src[i].size(), CV_32F);
+        dy_[i].create(src[i].size(), CV_32F);
+        for (int y = 0; y < src[i].rows; ++y)
+        {
+            const Point3f* dx_row = dx.ptr<Point3f>(y);
+            const Point3f* dy_row = dy.ptr<Point3f>(y);
+            float* dx_row_ = dx_[i].ptr<float>(y);
+            float* dy_row_ = dy_[i].ptr<float>(y);
+            for (int x = 0; x < src[i].cols; ++x)
+            {
+                dx_row_[x] = normL2(dx_row[x]);
+                dy_row_[x] = normL2(dy_row[x]);
+            }
+        }
+    }
+    PairwiseSeamFinder::find(src, corners, masks);
+}
+
+
+void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &img2, 
+                                                    const Mat &mask1, const Mat &mask2, GCGraph<float> &graph)
+{
+    const Size img_size = img1.size();
+
+    // Set terminal weights
+    for (int y = 0; y < img_size.height; ++y)
+    {
+        for (int x = 0; x < img_size.width; ++x)
+        {
+            int v = graph.addVtx();
+            graph.addTermWeights(v, mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f,
+                                    mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f);
+        }
+    }
+
+    // Set regular edge weights
+    const float weight_eps = 1.f;
+    for (int y = 0; y < img_size.height; ++y)
+    {
+        for (int x = 0; x < img_size.width; ++x)
+        {
+            int v = y * img_size.width + x;
+            if (x < img_size.width - 1)
+            {
+                float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
+                               normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1)) +
+                               weight_eps;
+                if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) ||
+                    !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1))
+                    weight += bad_region_penalty_;
+                graph.addEdges(v, v + 1, weight, weight);
+            }
+            if (y < img_size.height - 1)
+            {
+                float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
+                               normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x)) +
+                               weight_eps;
+                if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) ||
+                    !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x))
+                    weight += bad_region_penalty_;
+                graph.addEdges(v, v + img_size.width, weight, weight);
+            }
+        }
+    }
+}
+
+
+void GraphCutSeamFinder::Impl::setGraphWeightsColorGrad(
+        const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, 
+        const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, 
+        GCGraph<float> &graph)
+{
+    const Size img_size = img1.size();
+
+    // Set terminal weights
+    for (int y = 0; y < img_size.height; ++y)
+    {
+        for (int x = 0; x < img_size.width; ++x)
+        {
+            int v = graph.addVtx();
+            graph.addTermWeights(v, mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f,
+                                    mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f);
+        }
+    }
+
+    // Set regular edge weights
+    const float weight_eps = 1.f;
+    for (int y = 0; y < img_size.height; ++y)
+    {
+        for (int x = 0; x < img_size.width; ++x)
+        {
+            int v = y * img_size.width + x;
+            if (x < img_size.width - 1)
+            {
+                float grad = dx1.at<float>(y, x) + dx1.at<float>(y, x + 1) +
+                             dx2.at<float>(y, x) + dx2.at<float>(y, x + 1) + weight_eps;
+                float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
+                                normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1))) / grad + 
+                               weight_eps;
+                if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) ||
+                    !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1))
+                    weight += bad_region_penalty_;
+                graph.addEdges(v, v + 1, weight, weight);
+            }
+            if (y < img_size.height - 1)
+            {
+                float grad = dy1.at<float>(y, x) + dy1.at<float>(y + 1, x) + 
+                             dy2.at<float>(y, x) + dy2.at<float>(y + 1, x) + weight_eps;
+                float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 
+                                normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x))) / grad + 
+                               weight_eps;
+                if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) ||
+                    !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x))
+                    weight += bad_region_penalty_;
+                graph.addEdges(v, v + img_size.width, weight, weight);
+            }
+        }
+    }
+}
+
+
+void GraphCutSeamFinder::Impl::findInPair(size_t first, size_t second, Rect roi)
+{
+    Mat img1 = images_[first], img2 = images_[second];
+    Mat dx1 = dx_[first], dx2 = dx_[second];
+    Mat dy1 = dy_[first], dy2 = dy_[second];
+    Mat mask1 = masks_[first], mask2 = masks_[second];
+    Point tl1 = corners_[first], tl2 = corners_[second];
+
+    const int gap = 10;
+    Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3);
+    Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3);
+    Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
+    Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
+    Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
+    Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
+    Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
+    Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
+
+    // Cut subimages and submasks with some gap
+    for (int y = -gap; y < roi.height + gap; ++y)
+    {
+        for (int x = -gap; x < roi.width + gap; ++x)
+        {
+            int y1 = roi.y - tl1.y + y;
+            int x1 = roi.x - tl1.x + x;
+            if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols)
+            {
+                subimg1.at<Point3f>(y + gap, x + gap) = img1.at<Point3f>(y1, x1);
+                submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1);
+                subdx1.at<float>(y + gap, x + gap) = dx1.at<float>(y1, x1);
+                subdy1.at<float>(y + gap, x + gap) = dy1.at<float>(y1, x1);
+            }
+            else
+            {
+                subimg1.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0);
+                submask1.at<uchar>(y + gap, x + gap) = 0;
+                subdx1.at<float>(y + gap, x + gap) = 0.f;
+                subdy1.at<float>(y + gap, x + gap) = 0.f;
+            }
+
+            int y2 = roi.y - tl2.y + y;
+            int x2 = roi.x - tl2.x + x;
+            if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols)
+            {
+                subimg2.at<Point3f>(y + gap, x + gap) = img2.at<Point3f>(y2, x2);
+                submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2);
+                subdx2.at<float>(y + gap, x + gap) = dx2.at<float>(y2, x2);
+                subdy2.at<float>(y + gap, x + gap) = dy2.at<float>(y2, x2);
+            }
+            else
+            {
+                subimg2.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0);
+                submask2.at<uchar>(y + gap, x + gap) = 0;
+                subdx2.at<float>(y + gap, x + gap) = 0.f;
+                subdy2.at<float>(y + gap, x + gap) = 0.f;
+            }
+        }
+    }
+
+    const int vertex_count = (roi.height + 2 * gap) * (roi.width + 2 * gap);
+    const int edge_count = (roi.height - 1 + 2 * gap) * (roi.width + 2 * gap) +
+                           (roi.width - 1 + 2 * gap) * (roi.height + 2 * gap);
+    GCGraph<float> graph(vertex_count, edge_count);
+
+    switch (cost_type_)
+    {
+    case GraphCutSeamFinder::COST_COLOR:
+        setGraphWeightsColor(subimg1, subimg2, submask1, submask2, graph);
+        break;
+    case GraphCutSeamFinder::COST_COLOR_GRAD:
+        setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2, 
+                                 submask1, submask2, graph);
+        break;
+    default:
+        CV_Error(CV_StsBadArg, "unsupported pixel similarity measure");
+    }
+
+    graph.maxFlow();
+
+    for (int y = 0; y < roi.height; ++y)
+    {
+        for (int x = 0; x < roi.width; ++x)
+        {
+            if (graph.inSourceSegment((y + gap) * (roi.width + 2 * gap) + x + gap))
+            {
+                if (mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x))
+                    mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0;
+            }
+            else
+            {
+                if (mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x))
+                    mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0;
+            }
+        }
+    }
+}
+
+
+GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float bad_region_penalty)
+    : impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {}
+
+
+void GraphCutSeamFinder::find(const vector<Mat> &src, const vector<Point> &corners,
+                              vector<Mat> &masks)
+{
+    impl_->find(src, corners, masks);
+}
+
+} // namespace cv

modules/stitching/src/util.cpp

@@ -0,0 +1,167 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+void DisjointSets::createOneElemSets(int n)
+{
+    rank_.assign(n, 0);
+    size.assign(n, 1);
+    parent.resize(n);
+    for (int i = 0; i < n; ++i)
+        parent[i] = i;
+}
+
+
+int DisjointSets::findSetByElem(int elem)
+{
+    int set = elem;
+    while (set != parent[set])
+        set = parent[set];
+    int next;
+    while (elem != parent[elem]) 
+    {
+        next = parent[elem];
+        parent[elem] = set;
+        elem = next;
+    }
+    return set;
+}
+
+
+int DisjointSets::mergeSets(int set1, int set2)
+{
+    if (rank_[set1] < rank_[set2]) 
+    {
+        parent[set1] = set2;
+        size[set2] += size[set1];
+        return set2;
+    }
+    if (rank_[set2] < rank_[set1]) 
+    {
+        parent[set2] = set1;
+        size[set1] += size[set2];
+        return set1;
+    }
+    parent[set1] = set2;
+    rank_[set2]++;
+    size[set2] += size[set1];
+    return set2;
+}
+
+
+void Graph::addEdge(int from, int to, float weight)
+{
+    edges_[from].push_back(GraphEdge(from, to, weight));
+}
+
+
+bool overlapRoi(Point tl1, Point tl2, Size sz1, Size sz2, Rect &roi)
+{
+    int x_tl = max(tl1.x, tl2.x);
+    int y_tl = max(tl1.y, tl2.y);
+    int x_br = min(tl1.x + sz1.width, tl2.x + sz2.width);
+    int y_br = min(tl1.y + sz1.height, tl2.y + sz2.height);
+    if (x_tl < x_br && y_tl < y_br)
+    {
+        roi = Rect(x_tl, y_tl, x_br - x_tl, y_br - y_tl);
+        return true;
+    }
+    return false;
+}
+
+
+Rect resultRoi(const vector<Point> &corners, const vector<Mat> &images)
+{
+    vector<Size> sizes(images.size());
+    for (size_t i = 0; i < images.size(); ++i)
+        sizes[i] = images[i].size();
+    return resultRoi(corners, sizes);
+}
+
+
+Rect resultRoi(const vector<Point> &corners, const vector<Size> &sizes)
+{
+    CV_Assert(sizes.size() == corners.size());
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
+    Point br(numeric_limits<int>::min(), numeric_limits<int>::min());
+    for (size_t i = 0; i < corners.size(); ++i)
+    {
+        tl.x = min(tl.x, corners[i].x);
+        tl.y = min(tl.y, corners[i].y);
+        br.x = max(br.x, corners[i].x + sizes[i].width);
+        br.y = max(br.y, corners[i].y + sizes[i].height);
+    }
+    return Rect(tl, br);
+}
+
+
+Point resultTl(const vector<Point> &corners)
+{
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
+    for (size_t i = 0; i < corners.size(); ++i)
+    {
+        tl.x = min(tl.x, corners[i].x);
+        tl.y = min(tl.y, corners[i].y);
+    }
+    return tl;
+}
+
+
+void selectRandomSubset(int count, int size, vector<int> &subset)
+{
+    subset.clear();
+    for (int i = 0; i < size; ++i)
+    {
+        if (randu<int>() % (size - i) < count)
+        {
+            subset.push_back(i);
+            count--;
+        }
+    }
+}
+
+} // namespace cv

modules/stitching/src/warpers.cpp

@@ -0,0 +1,233 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  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
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., 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:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's 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.
+//
+//   * The name of the copyright holders may not 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 the Intel Corporation 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 "precomp.hpp"
+
+using namespace std;
+
+namespace cv
+{
+
+Ptr<Warper> Warper::createByCameraFocal(float focal, int type, bool try_gpu)
+{
+    bool can_use_gpu = try_gpu && gpu::getCudaEnabledDeviceCount();
+    if (type == PLANE)
+        return !can_use_gpu ? new PlaneWarper(focal) : new PlaneWarperGpu(focal);
+    if (type == CYLINDRICAL)
+        return !can_use_gpu ? new CylindricalWarper(focal) : new CylindricalWarperGpu(focal);
+    if (type == SPHERICAL)
+        return !can_use_gpu ? new SphericalWarper(focal) : new SphericalWarperGpu(focal);
+    CV_Error(CV_StsBadArg, "unsupported warping type");
+    return NULL;
+}
+
+
+void ProjectorBase::setTransformation(const Mat &R)
+{
+    CV_Assert(R.size() == Size(3, 3));
+    CV_Assert(R.type() == CV_32F);
+    r[0] = R.at<float>(0, 0); r[1] = R.at<float>(0, 1); r[2] = R.at<float>(0, 2);
+    r[3] = R.at<float>(1, 0); r[4] = R.at<float>(1, 1); r[5] = R.at<float>(1, 2);
+    r[6] = R.at<float>(2, 0); r[7] = R.at<float>(2, 1); r[8] = R.at<float>(2, 2);
+
+    Mat Rinv = R.inv();
+    rinv[0] = Rinv.at<float>(0, 0); rinv[1] = Rinv.at<float>(0, 1); rinv[2] = Rinv.at<float>(0, 2);
+    rinv[3] = Rinv.at<float>(1, 0); rinv[4] = Rinv.at<float>(1, 1); rinv[5] = Rinv.at<float>(1, 2);
+    rinv[6] = Rinv.at<float>(2, 0); rinv[7] = Rinv.at<float>(2, 1); rinv[8] = Rinv.at<float>(2, 2);
+}
+
+
+void PlaneWarper::detectResultRoi(Point &dst_tl, Point &dst_br)
+{
+    float tl_uf = numeric_limits<float>::max();
+    float tl_vf = numeric_limits<float>::max();
+    float br_uf = -numeric_limits<float>::max();
+    float br_vf = -numeric_limits<float>::max();
+
+    float u, v;
+
+    projector_.mapForward(0, 0, u, v);
+    tl_uf = min(tl_uf, u); tl_vf = min(tl_vf, v);
+    br_uf = max(br_uf, u); br_vf = max(br_vf, v);
+
+    projector_.mapForward(0, static_cast<float>(src_size_.height - 1), u, v);
+    tl_uf = min(tl_uf, u); tl_vf = min(tl_vf, v);
+    br_uf = max(br_uf, u); br_vf = max(br_vf, v);
+
+    projector_.mapForward(static_cast<float>(src_size_.width - 1), 0, u, v);
+    tl_uf = min(tl_uf, u); tl_vf = min(tl_vf, v);
+    br_uf = max(br_uf, u); br_vf = max(br_vf, v);
+
+    projector_.mapForward(static_cast<float>(src_size_.width - 1), static_cast<float>(src_size_.height - 1), u, v);
+    tl_uf = min(tl_uf, u); tl_vf = min(tl_vf, v);
+    br_uf = max(br_uf, u); br_vf = max(br_vf, v);
+
+    dst_tl.x = static_cast<int>(tl_uf);
+    dst_tl.y = static_cast<int>(tl_vf);
+    dst_br.x = static_cast<int>(br_uf);
+    dst_br.y = static_cast<int>(br_vf);
+}
+
+
+Point PlaneWarperGpu::warp(const Mat &src, float focal, const cv::Mat &R, cv::Mat &dst, int interp_mode, int border_mode)
+{
+    src_size_ = src.size();
+    projector_.size = src.size();
+    projector_.focal = focal;
+    projector_.setTransformation(R);
+
+    cv::Point dst_tl, dst_br;
+    detectResultRoi(dst_tl, dst_br);
+
+    gpu::buildWarpPlaneMaps(src.size(), Rect(dst_tl, Point(dst_br.x+1, dst_br.y+1)),
+                            R, focal, projector_.scale, projector_.plane_dist, d_xmap_, d_ymap_);
+
+    gpu::ensureSizeIsEnough(src.size(), src.type(), d_src_);
+    d_src_.upload(src);
+
+    gpu::ensureSizeIsEnough(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type(), d_dst_);
+
+    gpu::remap(d_src_, d_dst_, d_xmap_, d_ymap_, interp_mode, border_mode);
+
+    d_dst_.download(dst);
+
+    return dst_tl;
+}
+
+
+void SphericalWarper::detectResultRoi(Point &dst_tl, Point &dst_br)
+{
+    detectResultRoiByBorder(dst_tl, dst_br);
+
+    float tl_uf = static_cast<float>(dst_tl.x);
+    float tl_vf = static_cast<float>(dst_tl.y);
+    float br_uf = static_cast<float>(dst_br.x);
+    float br_vf = static_cast<float>(dst_br.y);
+
+    float x = projector_.rinv[1];
+    float y = projector_.rinv[4];
+    float z = projector_.rinv[7];
+    if (y > 0.f)
+    {
+        x = projector_.focal * x / z + src_size_.width * 0.5f;
+        y = projector_.focal * y / z + src_size_.height * 0.5f;
+        if (x > 0.f && x < src_size_.width && y > 0.f && y < src_size_.height)
+        {
+            tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(CV_PI * projector_.scale));
+            br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(CV_PI * projector_.scale));
+        }
+    }
+
+    x = projector_.rinv[1];
+    y = -projector_.rinv[4];
+    z = projector_.rinv[7];
+    if (y > 0.f)
+    {
+        x = projector_.focal * x / z + src_size_.width * 0.5f;
+        y = projector_.focal * y / z + src_size_.height * 0.5f;
+        if (x > 0.f && x < src_size_.width && y > 0.f && y < src_size_.height)
+        {
+            tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(0));
+            br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(0));
+        }
+    }
+
+    dst_tl.x = static_cast<int>(tl_uf);
+    dst_tl.y = static_cast<int>(tl_vf);
+    dst_br.x = static_cast<int>(br_uf);
+    dst_br.y = static_cast<int>(br_vf);
+}
+
+
+Point SphericalWarperGpu::warp(const Mat &src, float focal, const Mat &R, Mat &dst,
+                               int interp_mode, int border_mode)
+{
+    src_size_ = src.size();
+    projector_.size = src.size();
+    projector_.focal = focal;
+    projector_.setTransformation(R);
+
+    cv::Point dst_tl, dst_br;
+    detectResultRoi(dst_tl, dst_br);
+
+    gpu::buildWarpSphericalMaps(src.size(), Rect(dst_tl, Point(dst_br.x+1, dst_br.y+1)),
+                                R, focal, projector_.scale, d_xmap_, d_ymap_);
+
+    gpu::ensureSizeIsEnough(src.size(), src.type(), d_src_);
+    d_src_.upload(src);
+
+    gpu::ensureSizeIsEnough(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type(), d_dst_);
+
+    gpu::remap(d_src_, d_dst_, d_xmap_, d_ymap_, interp_mode, border_mode);
+
+    d_dst_.download(dst);
+
+    return dst_tl;
+}
+
+
+Point CylindricalWarperGpu::warp(const Mat &src, float focal, const Mat &R, Mat &dst,
+                                 int interp_mode, int border_mode)
+{
+    src_size_ = src.size();
+    projector_.size = src.size();
+    projector_.focal = focal;
+    projector_.setTransformation(R);
+
+    cv::Point dst_tl, dst_br;
+    detectResultRoi(dst_tl, dst_br);
+
+    gpu::buildWarpCylindricalMaps(src.size(), Rect(dst_tl, Point(dst_br.x+1, dst_br.y+1)),
+                                  R, focal, projector_.scale, d_xmap_, d_ymap_);
+
+    gpu::ensureSizeIsEnough(src.size(), src.type(), d_src_);
+    d_src_.upload(src);
+
+    gpu::ensureSizeIsEnough(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type(), d_dst_);
+
+    gpu::remap(d_src_, d_dst_, d_xmap_, d_ymap_, interp_mode, border_mode);
+
+    d_dst_.download(dst);
+
+    return dst_tl;
+}
+
+} // namespace cv