fixed bug in opencv_stitching (corrected resize images step), added matches…

fixed bug in opencv_stitching (corrected resize images step), added matches checking (both 1->2 and 2->1 must be presented)

modules/stitching/blenders.cpp

@@ -38,312 +38,312 @@
 // 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 "blenders.hpp"
-#include "util.hpp"
-
-using namespace std;
-using namespace cv;
-
-static const float WEIGHT_EPS = 1e-5f;
-
-Ptr<Blender> Blender::createDefault(int type)
-{
-    if (type == NO)
-        return new Blender();
-    if (type == FEATHER)
-        return new FeatherBlender();
-    if (type == MULTI_BAND)
-        return new MultiBandBlender();
-    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)
-{
-    normalize(dst_weight_map_, dst_);
-    dst_mask_ = dst_weight_map_ > WEIGHT_EPS;
-    Blender::blend(dst, dst_mask);
-}
-
-
-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
-    vector<Mat> src_pyr_gauss(num_bands_ + 1);
-    copyMakeBorder(img, src_pyr_gauss[0], top, bottom, left, right, 
-                   BORDER_REFLECT);
-    for (int i = 0; i < num_bands_; ++i)
-        pyrDown(src_pyr_gauss[i], src_pyr_gauss[i + 1]);
-    vector<Mat> src_pyr_laplace;
-    createLaplacePyr(src_pyr_gauss, src_pyr_laplace);
-    src_pyr_gauss.clear();
-
-    // 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)
-        normalize(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 normalize(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 vector<Mat> &pyr_gauss, vector<Mat> &pyr_laplace)
-{
-    if (pyr_gauss.size() == 0)
-        return;
-    pyr_laplace.resize(pyr_gauss.size());
-    Mat tmp;
-    for (size_t i = 0; i < pyr_laplace.size() - 1; ++i)
-    {
-        pyrUp(pyr_gauss[i + 1], tmp, pyr_gauss[i].size());
-        subtract(pyr_gauss[i], tmp, pyr_laplace[i]);
-    }
-    pyr_laplace[pyr_laplace.size() - 1] = pyr_gauss[pyr_laplace.size() - 1].clone();
-}
-
-
-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]);
-    }
-}
-
-
+//M*/
+#include "blenders.hpp"
+#include "util.hpp"
+
+using namespace std;
+using namespace cv;
+
+static const float WEIGHT_EPS = 1e-5f;
+
+Ptr<Blender> Blender::createDefault(int type)
+{
+    if (type == NO)
+        return new Blender();
+    if (type == FEATHER)
+        return new FeatherBlender();
+    if (type == MULTI_BAND)
+        return new MultiBandBlender();
+    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)
+{
+    normalize(dst_weight_map_, dst_);
+    dst_mask_ = dst_weight_map_ > WEIGHT_EPS;
+    Blender::blend(dst, dst_mask);
+}
+
+
+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
+    vector<Mat> src_pyr_gauss(num_bands_ + 1);
+    copyMakeBorder(img, src_pyr_gauss[0], top, bottom, left, right, 
+                   BORDER_REFLECT);
+    for (int i = 0; i < num_bands_; ++i)
+        pyrDown(src_pyr_gauss[i], src_pyr_gauss[i + 1]);
+    vector<Mat> src_pyr_laplace;
+    createLaplacePyr(src_pyr_gauss, src_pyr_laplace);
+    src_pyr_gauss.clear();
+
+    // 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)
+        normalize(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 normalize(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 vector<Mat> &pyr_gauss, vector<Mat> &pyr_laplace)
+{
+    if (pyr_gauss.size() == 0)
+        return;
+    pyr_laplace.resize(pyr_gauss.size());
+    Mat tmp;
+    for (size_t i = 0; i < pyr_laplace.size() - 1; ++i)
+    {
+        pyrUp(pyr_gauss[i + 1], tmp, pyr_gauss[i].size());
+        subtract(pyr_gauss[i], tmp, pyr_laplace[i]);
+    }
+    pyr_laplace[pyr_laplace.size() - 1] = pyr_gauss[pyr_laplace.size() - 1].clone();
+}
+
+
+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]);
+    }
+}
+
+

modules/stitching/main.cpp

@@ -75,7 +75,7 @@ void printUsage()
         "  --work_megapix <float>\n"
         "      Resolution for image registration step. The default is 0.6 Mpx.\n"
         "  --match_conf <float>\n"
-        "      Confidence for feature matching step. The default is 0.7.\n"
+        "      Confidence for feature matching step. The default is 0.65.\n"
         "  --conf_thresh <float>\n"
         "      Threshold for two images are from the same panorama confidence.\n"
         "      The default is 1.0.\n"
@@ -320,11 +320,14 @@ int main(int argc, char* argv[])
     Mat full_img, img;
 
     vector<Mat> images(num_images);
+    vector<Size> full_img_sizes(num_images);
     double seam_work_aspect = 1;
 
     for (int i = 0; i < num_images; ++i)
     {
         full_img = imread(img_names[i]);
+        full_img_sizes[i] = full_img.size();
+
         if (full_img.empty())
         {
             LOGLN("Can't open image " << img_names[i]);
@@ -376,14 +379,17 @@ int main(int argc, char* argv[])
     vector<int> indices = leaveBiggestComponent(features, pairwise_matches, conf_thresh);
     vector<Mat> img_subset;
     vector<string> img_names_subset;
+    vector<Size> full_img_sizes_subset;
     for (size_t i = 0; i < indices.size(); ++i)
     {
         img_names_subset.push_back(img_names[indices[i]]);
         img_subset.push_back(images[indices[i]]);
+        full_img_sizes_subset.push_back(full_img_sizes[indices[i]]);
     }
 
     images = img_subset;
     img_names = img_names_subset;
+    full_img_sizes = full_img_sizes_subset;
 
     // Check if we still have enough images
     num_images = static_cast<int>(img_names.size());
@@ -519,16 +525,21 @@ int main(int argc, char* argv[])
             warper = Warper::createByCameraFocal(warped_image_scale, warp_type);
 
             // Update corners and sizes
-            Rect dst_roi = resultRoi(corners, sizes);
             for (int i = 0; i < num_images; ++i)
             {
                 // Update camera focal
                 cameras[i].focal *= compose_work_aspect;
 
                 // Update corner and size
-                corners[i] = dst_roi.tl() + (corners[i] - dst_roi.tl()) * compose_seam_aspect;
-                sizes[i] = Size(static_cast<int>((sizes[i].width + 1) * compose_seam_aspect), 
-                                static_cast<int>((sizes[i].height + 1) * compose_seam_aspect));
+                Size sz = full_img_sizes[i];
+                if (abs(compose_scale - 1) > 1e-1)
+                {
+                    sz.width = cvRound(full_img_sizes[i].width * compose_scale);
+                    sz.height = cvRound(full_img_sizes[i].height * compose_scale);
+                }
+                Rect roi = warper->warpRoi(sz, static_cast<float>(cameras[i].focal), cameras[i].R);
+                corners[i] = roi.tl();
+                sizes[i] = roi.size();
             }
         }
         if (abs(compose_scale - 1) > 1e-1)
@@ -539,7 +550,7 @@ int main(int argc, char* argv[])
         Size img_size = img.size();
 
         // Warp the current image
-        warper->warp(img, static_cast<float>(cameras[img_idx].focal), cameras[img_idx].R, 
+        warper->warp(img, static_cast<float>(cameras[img_idx].focal), cameras[img_idx].R,
                      img_warped);
 
         // Warp the current image mask
@@ -587,7 +598,7 @@ int main(int argc, char* argv[])
     }
    
     Mat result, result_mask;
-    blender->blend(result, result_mask);
+    blender->blend(result, result_mask);    
 
     LOGLN("Compositing, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
 

modules/stitching/matchers.cpp

@@ -311,7 +311,7 @@ namespace
             const DMatch& m1 = pair_matches[i][1];
             if (m0.distance < (1.f - match_conf_) * m1.distance)
             {
-                matches_info.matches.push_back(m0);
+                //matches_info.matches.push_back(m0);
                 matches.insert(make_pair(m0.queryIdx, m0.trainIdx));
             }
         }
@@ -326,7 +326,7 @@ namespace
             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())
+                if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) != matches.end())
                     matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance));
         }
     }
@@ -352,7 +352,7 @@ namespace
             const DMatch& m1 = pair_matches[i][1];
             if (m0.distance < (1.f - match_conf_) * m1.distance)
             {
-                matches_info.matches.push_back(m0);
+                //matches_info.matches.push_back(m0);
                 matches.insert(make_pair(m0.queryIdx, m0.trainIdx));
             }
         }
@@ -368,7 +368,7 @@ namespace
             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())
+                if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) != matches.end())
                     matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance));
         }
     }

modules/stitching/motion_estimators.cpp

@@ -38,478 +38,478 @@
 // 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 <algorithm>
-#include "autocalib.hpp"
-#include "motion_estimators.hpp"
-#include "util.hpp"
-
-using namespace std;
-using 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;
-}
-
-
-//////////////////////////////////////////////////////////////////////////////
-
-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)
-{
-    const int num_images = static_cast<int>(features.size());
-
-    // 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));
-}
-
-
-//////////////////////////////////////////////////////////////////////////////
-
-void BundleAdjuster::estimate(const vector<ImageFeatures> &features, const vector<MatchesInfo> &pairwise_matches, 
-                              vector<CameraParams> &cameras)
-{
-    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;
-}
-
-
-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;
-            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;
-    }
-}
-
-
-//////////////////////////////////////////////////////////////////////////////
-
-void waveCorrect(vector<Mat> &rmats)
-{
-    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];
-}
-
-
-//////////////////////////////////////////////////////////////////////////////
-
-vector<int> leaveBiggestComponent(vector<ImageFeatures> &features,  vector<MatchesInfo> &pairwise_matches, 
-                                  float conf_threshold)
-{
-    const int num_images = static_cast<int>(features.size());
-
-    DjSets 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.find(i);
-            int comp2 = comps.find(j);
-            if (comp1 != comp2) 
-                comps.merge(comp1, comp2);
-        }
-    }
-
-    int max_comp = 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.find(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 = i;
-            pairwise_matches_subset.back().dst_img_idx = 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(")");
-
-    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));
-        }
-    }
-
-    DjSets 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.find(edges[i].from);
-        int comp2 = comps.find(edges[i].to);
-        if (comp1 != comp2)
-        {
-            comps.merge(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);
-}
+//M*/
+#include <algorithm>
+#include "autocalib.hpp"
+#include "motion_estimators.hpp"
+#include "util.hpp"
+
+using namespace std;
+using 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;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+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)
+{
+    const int num_images = static_cast<int>(features.size());
+
+    // 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));
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+void BundleAdjuster::estimate(const vector<ImageFeatures> &features, const vector<MatchesInfo> &pairwise_matches, 
+                              vector<CameraParams> &cameras)
+{
+    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;
+}
+
+
+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;
+            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;
+    }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+void waveCorrect(vector<Mat> &rmats)
+{
+    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];
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+vector<int> leaveBiggestComponent(vector<ImageFeatures> &features,  vector<MatchesInfo> &pairwise_matches, 
+                                  float conf_threshold)
+{
+    const int num_images = static_cast<int>(features.size());
+
+    DjSets 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.find(i);
+            int comp2 = comps.find(j);
+            if (comp1 != comp2) 
+                comps.merge(comp1, comp2);
+        }
+    }
+
+    int max_comp = 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.find(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 = i;
+            pairwise_matches_subset.back().dst_img_idx = 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 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));
+        }
+    }
+
+    DjSets 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.find(edges[i].from);
+        int comp2 = comps.find(edges[i].to);
+        if (comp1 != comp2)
+        {
+            comps.merge(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);
+}

modules/stitching/seam_finders.cpp

@@ -38,365 +38,365 @@
 // 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 "seam_finders.hpp"
-#include "util.hpp"
-
-using namespace std;
-using 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)
-{
-    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);
-        }
-    }
-}
-
-
-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:
+//M*/
+#include "seam_finders.hpp"
+#include "util.hpp"
+
+using namespace std;
+using 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)
+{
+    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);
+        }
+    }
+}
+
+
+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);
+                              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);
-}
-
-
+                                  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 = 1e-3f;
-    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);
-            }
-        }
-    }
-}
-
-
+                                                    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 = 1e-3f;
-    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)
-            {
+        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)
-            {
+                             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);
-}
+                             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);
+}

modules/stitching/warpers.hpp

@@ -38,119 +38,122 @@
 // 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_WARPERS_HPP__
-#define __OPENCV_WARPERS_HPP__
-
-#include "precomp.hpp"
-
-class Warper
-{
-public:
-    enum { PLANE, CYLINDRICAL, SPHERICAL };
-    static cv::Ptr<Warper> createByCameraFocal(float focal, int type);
-
-    virtual ~Warper() {}
-    virtual cv::Point warp(const cv::Mat &src, float focal, const cv::Mat& R, cv::Mat &dst,
-                           int interp_mode = cv::INTER_LINEAR, int border_mode = cv::BORDER_REFLECT) = 0;
-};
-
-
-struct ProjectorBase
-{
-    void setTransformation(const cv::Mat& R);
-
-    cv::Size size;
-    float focal;
-    float r[9];
-    float rinv[9];
-    float scale;
-};
-
-
-template <class P>
-class WarperBase : public Warper
-{   
-public:
-    cv::Point warp(const cv::Mat &src, float focal, const cv::Mat &R, cv::Mat &dst,
-                   int interp_mode, int border_mode);
-
-protected:
-    // Detects ROI of the destination image. It's correct for any projection.
-    virtual void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br);
-
-    // Detects ROI of the destination image by walking over image border.
-    // Correctness for any projection isn't guaranteed.
-    void detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br);
-
-    cv::Size src_size_;
-    P projector_;
-};
-
-
-struct PlaneProjector : ProjectorBase
-{
-    void mapForward(float x, float y, float &u, float &v);
-    void mapBackward(float u, float v, float &x, float &y);
-
-    float plane_dist;
-};
-
-
-// Projects image onto z = plane_dist plane
-class PlaneWarper : public WarperBase<PlaneProjector>
-{
-public:
-    PlaneWarper(float plane_dist = 1.f, float scale = 1.f)
-    {
-        projector_.plane_dist = plane_dist;
-        projector_.scale = scale;
-    }
-
-private:
-    void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br);
-};
-
-
-struct SphericalProjector : ProjectorBase
-{
-    void mapForward(float x, float y, float &u, float &v);
-    void mapBackward(float u, float v, float &x, float &y);
-};
-
-
-// Projects image onto unit sphere with origin at (0, 0, 0).
-// Poles are located at (0, -1, 0) and (0, 1, 0) points.
-class SphericalWarper : public WarperBase<SphericalProjector>
-{
-public:
-    SphericalWarper(float scale = 300.f) { projector_.scale = scale; }
-
-private:  
-    void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br);
-};
-
-
-struct CylindricalProjector : ProjectorBase
-{
-    void mapForward(float x, float y, float &u, float &v);
-    void mapBackward(float u, float v, float &x, float &y);
-};
-
-
-// Projects image onto x * x + z * z = 1 cylinder
-class CylindricalWarper : public WarperBase<CylindricalProjector>
-{
-public:
-    CylindricalWarper(float scale = 300.f) { projector_.scale = scale; }
-
-private:
-    void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br)
-    {
-        WarperBase<CylindricalProjector>::detectResultRoiByBorder(dst_tl, dst_br);
-    }
-};
-
-#include "warpers_inl.hpp"
-
-#endif // __OPENCV_WARPERS_HPP__
+//M*/
+#ifndef __OPENCV_WARPERS_HPP__
+#define __OPENCV_WARPERS_HPP__
+
+#include "precomp.hpp"
+
+class Warper
+{
+public:
+    enum { PLANE, CYLINDRICAL, SPHERICAL };
+    static cv::Ptr<Warper> createByCameraFocal(float focal, int type);
+
+    virtual ~Warper() {}
+    virtual cv::Point warp(const cv::Mat &src, float focal, const cv::Mat& R, cv::Mat &dst,
+                           int interp_mode = cv::INTER_LINEAR, int border_mode = cv::BORDER_REFLECT) = 0;
+    virtual cv::Rect warpRoi(const cv::Size &sz, float focal, const cv::Mat &R) = 0;
+};
+
+
+struct ProjectorBase
+{
+    void setTransformation(const cv::Mat& R);
+
+    cv::Size size;
+    float focal;
+    float r[9];
+    float rinv[9];
+    float scale;
+};
+
+
+template <class P>
+class WarperBase : public Warper
+{   
+public:
+    cv::Point warp(const cv::Mat &src, float focal, const cv::Mat &R, cv::Mat &dst,
+                   int interp_mode, int border_mode);
+
+    cv::Rect warpRoi(const cv::Size &sz, float focal, const cv::Mat &R);
+
+protected:
+    // Detects ROI of the destination image. It's correct for any projection.
+    virtual void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br);
+
+    // Detects ROI of the destination image by walking over image border.
+    // Correctness for any projection isn't guaranteed.
+    void detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br);
+
+    cv::Size src_size_;
+    P projector_;
+};
+
+
+struct PlaneProjector : ProjectorBase
+{
+    void mapForward(float x, float y, float &u, float &v);
+    void mapBackward(float u, float v, float &x, float &y);
+
+    float plane_dist;
+};
+
+
+// Projects image onto z = plane_dist plane
+class PlaneWarper : public WarperBase<PlaneProjector>
+{
+public:
+    PlaneWarper(float plane_dist = 1.f, float scale = 1.f)
+    {
+        projector_.plane_dist = plane_dist;
+        projector_.scale = scale;
+    }
+
+private:
+    void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br);
+};
+
+
+struct SphericalProjector : ProjectorBase
+{
+    void mapForward(float x, float y, float &u, float &v);
+    void mapBackward(float u, float v, float &x, float &y);
+};
+
+
+// Projects image onto unit sphere with origin at (0, 0, 0).
+// Poles are located at (0, -1, 0) and (0, 1, 0) points.
+class SphericalWarper : public WarperBase<SphericalProjector>
+{
+public:
+    SphericalWarper(float scale = 300.f) { projector_.scale = scale; }
+
+private:  
+    void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br);
+};
+
+
+struct CylindricalProjector : ProjectorBase
+{
+    void mapForward(float x, float y, float &u, float &v);
+    void mapBackward(float u, float v, float &x, float &y);
+};
+
+
+// Projects image onto x * x + z * z = 1 cylinder
+class CylindricalWarper : public WarperBase<CylindricalProjector>
+{
+public:
+    CylindricalWarper(float scale = 300.f) { projector_.scale = scale; }
+
+private:
+    void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br)
+    {
+        WarperBase<CylindricalProjector>::detectResultRoiByBorder(dst_tl, dst_br);
+    }
+};
+
+#include "warpers_inl.hpp"
+
+#endif // __OPENCV_WARPERS_HPP__

modules/stitching/warpers_inl.hpp

@@ -38,202 +38,218 @@
 // 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_WARPERS_INL_HPP__
-#define __OPENCV_WARPERS_INL_HPP__
-
-#include "warpers.hpp" // Make your IDE see declarations
-
-template <class P>
-cv::Point WarperBase<P>::warp(const cv::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);
-
-    cv::Mat xmap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
-    cv::Mat ymap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
-
-    float x, y;
-    for (int v = dst_tl.y; v <= dst_br.y; ++v)
-    {
-        for (int u = dst_tl.x; u <= dst_br.x; ++u)
-        {
-            projector_.mapBackward(static_cast<float>(u), static_cast<float>(v), x, y);
-            xmap.at<float>(v - dst_tl.y, u - dst_tl.x) = x;
-            ymap.at<float>(v - dst_tl.y, u - dst_tl.x) = y;
-        }
-    }
-
-    dst.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type());
-    remap(src, dst, xmap, ymap, interp_mode, border_mode);
-
-    return dst_tl;
-}
-
-
-template <class P>
-void WarperBase<P>::detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br)
-{
-    float tl_uf = std::numeric_limits<float>::max();
-    float tl_vf = std::numeric_limits<float>::max();
-    float br_uf = -std::numeric_limits<float>::max();
-    float br_vf = -std::numeric_limits<float>::max();
-
-    float u, v;
-    for (int y = 0; y < src_size_.height; ++y)
-    {
-        for (int x = 0; x < src_size_.width; ++x)
-        {
-            projector_.mapForward(static_cast<float>(x), static_cast<float>(y), u, v);
-            tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
-            br_uf = std::max(br_uf, u); br_vf = std::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);
-}
-
-
-template <class P>
-void WarperBase<P>::detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br)
-{
-    float tl_uf = std::numeric_limits<float>::max();
-    float tl_vf = std::numeric_limits<float>::max();
-    float br_uf = -std::numeric_limits<float>::max();
-    float br_vf = -std::numeric_limits<float>::max();
-
-    float u, v;
-    for (float x = 0; x < src_size_.width; ++x)
-    {
-        projector_.mapForward(static_cast<float>(x), 0, u, v);
-        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
-        br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v);
-
-        projector_.mapForward(static_cast<float>(x), static_cast<float>(src_size_.height - 1), u, v);
-        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
-        br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v);
-    }
-    for (int y = 0; y < src_size_.height; ++y)
-    {
-        projector_.mapForward(0, static_cast<float>(y), u, v);
-        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
-        br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v);
-
-        projector_.mapForward(static_cast<float>(src_size_.width - 1), static_cast<float>(y), u, v);
-        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
-        br_uf = std::max(br_uf, u); br_vf = std::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);
-}
-
-
-inline
-void PlaneProjector::mapForward(float x, float y, float &u, float &v)
-{
-    x -= size.width * 0.5f;
-    y -= size.height * 0.5f;
-
-    float x_ = r[0] * x + r[1] * y + r[2] * focal;
-    float y_ = r[3] * x + r[4] * y + r[5] * focal;
-    float z_ = r[6] * x + r[7] * y + r[8] * focal;
-
-    u = scale * x_ / z_ * plane_dist;
-    v = scale * y_ / z_ * plane_dist;
-}
-
-
-inline
-void PlaneProjector::mapBackward(float u, float v, float &x, float &y)
-{
-    float x_ = u / scale;
-    float y_ = v / scale;
-
-    float z;
-    x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * plane_dist;
-    y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * plane_dist;
-    z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * plane_dist;
-
-    x = focal * x / z + size.width * 0.5f;
-    y = focal * y / z + size.height * 0.5f;
-}
-
-
-inline
-void SphericalProjector::mapForward(float x, float y, float &u, float &v)
-{
-    x -= size.width * 0.5f;
-    y -= size.height * 0.5f;
-
-    float x_ = r[0] * x + r[1] * y + r[2] * focal;
-    float y_ = r[3] * x + r[4] * y + r[5] * focal;
-    float z_ = r[6] * x + r[7] * y + r[8] * focal;
-
-    u = scale * atan2f(x_, z_);
-    v = scale * (static_cast<float>(CV_PI) - acosf(y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_)));
-}
-
-
-inline
-void SphericalProjector::mapBackward(float u, float v, float &x, float &y)
-{
-    float sinv = sinf(static_cast<float>(CV_PI) - v / scale);
-    float x_ = sinv * sinf(u / scale);
-    float y_ = cosf(static_cast<float>(CV_PI) - v / scale);
-    float z_ = sinv * cosf(u / scale);
-
-    float z;
-    x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_;
-    y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_;
-    z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_;
-
-    x = focal * x / z + size.width * 0.5f;
-    y = focal * y / z + size.height * 0.5f;
-}
-
-
-inline
-void CylindricalProjector::mapForward(float x, float y, float &u, float &v)
-{
-    x -= size.width * 0.5f;
-    y -= size.height * 0.5f;
-
-    float x_ = r[0] * x + r[1] * y + r[2] * focal;
-    float y_ = r[3] * x + r[4] * y + r[5] * focal;
-    float z_ = r[6] * x + r[7] * y + r[8] * focal;
-
-    u = scale * atan2f(x_, z_);
-    v = scale * y_ / sqrtf(x_ * x_ + z_ * z_);
-}
-
-
-inline
-void CylindricalProjector::mapBackward(float u, float v, float &x, float &y)
-{
-    float x_ = sinf(u / scale);
-    float y_ = v / scale;
-    float z_ = cosf(u / scale);
-
-    float z;
-    x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_;
-    y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_;
-    z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_;
-
-    x = focal * x / z + size.width * 0.5f;
-    y = focal * y / z + size.height * 0.5f;
-}
-
-#endif // __OPENCV_WARPERS_INL_HPP__
+//M*/
+#ifndef __OPENCV_WARPERS_INL_HPP__
+#define __OPENCV_WARPERS_INL_HPP__
+
+#include "warpers.hpp" // Make your IDE see declarations
+
+template <class P>
+cv::Point WarperBase<P>::warp(const cv::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);
+
+    cv::Mat xmap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
+    cv::Mat ymap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
+
+    float x, y;
+    for (int v = dst_tl.y; v <= dst_br.y; ++v)
+    {
+        for (int u = dst_tl.x; u <= dst_br.x; ++u)
+        {
+            projector_.mapBackward(static_cast<float>(u), static_cast<float>(v), x, y);
+            xmap.at<float>(v - dst_tl.y, u - dst_tl.x) = x;
+            ymap.at<float>(v - dst_tl.y, u - dst_tl.x) = y;
+        }
+    }
+
+    dst.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type());
+    remap(src, dst, xmap, ymap, interp_mode, border_mode);
+
+    return dst_tl;
+}
+
+
+template <class P>
+cv::Rect WarperBase<P>::warpRoi(const cv::Size &sz, float focal, const cv::Mat &R)
+{
+    src_size_ = sz;
+
+    projector_.size = sz;
+    projector_.focal = focal;
+    projector_.setTransformation(R);
+
+    cv::Point dst_tl, dst_br;
+    detectResultRoi(dst_tl, dst_br);
+
+    return cv::Rect(dst_tl, cv::Point(dst_br.x + 1, dst_br.y + 1));
+}
+
+
+template <class P>
+void WarperBase<P>::detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br)
+{
+    float tl_uf = std::numeric_limits<float>::max();
+    float tl_vf = std::numeric_limits<float>::max();
+    float br_uf = -std::numeric_limits<float>::max();
+    float br_vf = -std::numeric_limits<float>::max();
+
+    float u, v;
+    for (int y = 0; y < src_size_.height; ++y)
+    {
+        for (int x = 0; x < src_size_.width; ++x)
+        {
+            projector_.mapForward(static_cast<float>(x), static_cast<float>(y), u, v);
+            tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
+            br_uf = std::max(br_uf, u); br_vf = std::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);
+}
+
+
+template <class P>
+void WarperBase<P>::detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br)
+{
+    float tl_uf = std::numeric_limits<float>::max();
+    float tl_vf = std::numeric_limits<float>::max();
+    float br_uf = -std::numeric_limits<float>::max();
+    float br_vf = -std::numeric_limits<float>::max();
+
+    float u, v;
+    for (float x = 0; x < src_size_.width; ++x)
+    {
+        projector_.mapForward(static_cast<float>(x), 0, u, v);
+        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
+        br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v);
+
+        projector_.mapForward(static_cast<float>(x), static_cast<float>(src_size_.height - 1), u, v);
+        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
+        br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v);
+    }
+    for (int y = 0; y < src_size_.height; ++y)
+    {
+        projector_.mapForward(0, static_cast<float>(y), u, v);
+        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
+        br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v);
+
+        projector_.mapForward(static_cast<float>(src_size_.width - 1), static_cast<float>(y), u, v);
+        tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v);
+        br_uf = std::max(br_uf, u); br_vf = std::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);
+}
+
+
+inline
+void PlaneProjector::mapForward(float x, float y, float &u, float &v)
+{
+    x -= size.width * 0.5f;
+    y -= size.height * 0.5f;
+
+    float x_ = r[0] * x + r[1] * y + r[2] * focal;
+    float y_ = r[3] * x + r[4] * y + r[5] * focal;
+    float z_ = r[6] * x + r[7] * y + r[8] * focal;
+
+    u = scale * x_ / z_ * plane_dist;
+    v = scale * y_ / z_ * plane_dist;
+}
+
+
+inline
+void PlaneProjector::mapBackward(float u, float v, float &x, float &y)
+{
+    float x_ = u / scale;
+    float y_ = v / scale;
+
+    float z;
+    x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * plane_dist;
+    y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * plane_dist;
+    z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * plane_dist;
+
+    x = focal * x / z + size.width * 0.5f;
+    y = focal * y / z + size.height * 0.5f;
+}
+
+
+inline
+void SphericalProjector::mapForward(float x, float y, float &u, float &v)
+{
+    x -= size.width * 0.5f;
+    y -= size.height * 0.5f;
+
+    float x_ = r[0] * x + r[1] * y + r[2] * focal;
+    float y_ = r[3] * x + r[4] * y + r[5] * focal;
+    float z_ = r[6] * x + r[7] * y + r[8] * focal;
+
+    u = scale * atan2f(x_, z_);
+    v = scale * (static_cast<float>(CV_PI) - acosf(y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_)));
+}
+
+
+inline
+void SphericalProjector::mapBackward(float u, float v, float &x, float &y)
+{
+    float sinv = sinf(static_cast<float>(CV_PI) - v / scale);
+    float x_ = sinv * sinf(u / scale);
+    float y_ = cosf(static_cast<float>(CV_PI) - v / scale);
+    float z_ = sinv * cosf(u / scale);
+
+    float z;
+    x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_;
+    y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_;
+    z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_;
+
+    x = focal * x / z + size.width * 0.5f;
+    y = focal * y / z + size.height * 0.5f;
+}
+
+
+inline
+void CylindricalProjector::mapForward(float x, float y, float &u, float &v)
+{
+    x -= size.width * 0.5f;
+    y -= size.height * 0.5f;
+
+    float x_ = r[0] * x + r[1] * y + r[2] * focal;
+    float y_ = r[3] * x + r[4] * y + r[5] * focal;
+    float z_ = r[6] * x + r[7] * y + r[8] * focal;
+
+    u = scale * atan2f(x_, z_);
+    v = scale * y_ / sqrtf(x_ * x_ + z_ * z_);
+}
+
+
+inline
+void CylindricalProjector::mapBackward(float u, float v, float &x, float &y)
+{
+    float x_ = sinf(u / scale);
+    float y_ = v / scale;
+    float z_ = cosf(u / scale);
+
+    float z;
+    x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_;
+    y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_;
+    z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_;
+
+    x = focal * x / z + size.width * 0.5f;
+    y = focal * y / z + size.height * 0.5f;
+}
+
+#endif // __OPENCV_WARPERS_INL_HPP__

samples/gpu/driver_api_stereo_multi.cpp

@@ -32,7 +32,11 @@ int main()
 
 #include <cuda.h>
 #include <cuda_runtime.h>
+#include <GL/gl.h>
+#include <cudaGL.h>
 #include "opencv2/core/internal.hpp" // For TBB wrappers
+#include "tbb/tbb.h"
+#include "tbb/mutex.h"
 
 using namespace std;
 using namespace cv;
@@ -54,7 +58,7 @@ inline void safeCall_(int code, const char* expr, const char* file, int line)
 }
 
 // Each GPU is associated with its own context
-CUcontext contexts[2];
+CUcontext contexts[/*2*/1];
 
 void inline contextOn(int id) 
 {
@@ -76,6 +80,10 @@ GpuMat d_result[2];
 // CPU result
 Mat result;
 
+int some[2];
+
+tbb::mutex mutex;
+
 int main(int argc, char** argv)
 {
     if (argc < 3)
@@ -85,11 +93,11 @@ int main(int argc, char** argv)
     }
 
     int num_devices = getCudaEnabledDeviceCount();
-    if (num_devices < 2)
-    {
-        std::cout << "Two or more GPUs are required\n";
-        return -1;
-    }
+//    if (num_devices < 2)
+//    {
+//        std::cout << "Two or more GPUs are required\n";
+//        return -1;
+//    }
 
     for (int i = 0; i < num_devices; ++i)
     {
@@ -123,13 +131,14 @@ int main(int argc, char** argv)
     // Create context for GPU #0
     CUdevice device;
     safeCall(cuDeviceGet(&device, 0));
-    safeCall(cuCtxCreate(&contexts[0], 0, device));
+    safeCall(cuGLCtxCreate(&contexts[0], 0, device));
+    //safeCall(cuCtxCreate(&contexts[0], 0, device));
     contextOff();
 
-    // Create context for GPU #1
-    safeCall(cuDeviceGet(&device, 1));
-    safeCall(cuCtxCreate(&contexts[1], 0, device));
-    contextOff();
+//    // Create context for GPU #1
+//    safeCall(cuDeviceGet(&device, 0));
+//    safeCall(cuCtxCreate(&contexts[1], 0, device));
+//    contextOff();
 
     // Split source images for processing on GPU #0
     contextOn(0);
@@ -139,15 +148,20 @@ int main(int argc, char** argv)
     contextOff();
 
     // Split source images for processing on the GPU #1
-    contextOn(1);
+    contextOn(0);
     d_left[1].upload(left.rowRange(left.rows / 2, left.rows));
     d_right[1].upload(right.rowRange(right.rows / 2, right.rows));
     bm[1] = new StereoBM_GPU();
     contextOff();
 
+    some[0] = some[1] = 0;
     // Execute calculation in two threads using two GPUs
-    int devices[] = {0, 1};
-    parallel_do(devices, devices + 2, Worker());
+    vector<int> devices;
+    for (int i = 0; i < 4; ++i)
+        devices.push_back(rand()%2);
+    tbb::parallel_do(&devices[0], &devices[devices.size() - 1], Worker());
+
+    cout << some[0] << " " << some[1] << endl;
 
     // Release the first GPU resources
     contextOn(0);
@@ -159,7 +173,7 @@ int main(int argc, char** argv)
     contextOff();
 
     // Release the second GPU resources
-    contextOn(1);
+    contextOn(0);
     imshow("GPU #1 result", Mat(d_result[1]));
     d_left[1].release();
     d_right[1].release();
@@ -175,7 +189,9 @@ int main(int argc, char** argv)
 
 void Worker::operator()(int device_id) const
 {
-    contextOn(device_id);
+    mutex.lock();
+
+    contextOn(0);
 
     bm[device_id]->operator()(d_left[device_id], d_right[device_id],
                               d_result[device_id]);
@@ -184,13 +200,16 @@ void Worker::operator()(int device_id) const
         << "): finished\n";
 
     contextOff();
+
+    mutex.unlock();
 }
 
 
 void destroyContexts()
 {
     safeCall(cuCtxDestroy(contexts[0]));
-    safeCall(cuCtxDestroy(contexts[1]));
+    //safeCall(cuCtxDestroy(contexts[1]));
 }
 
 #endif
+