From a39bce204d05e382f99cfb8819bf78b83b5f15c7 Mon Sep 17 00:00:00 2001 From: "alexey.spizhevoy" Date: Thu, 2 Aug 2012 11:37:47 +0400 Subject: [PATCH] implemented DP-based seam estimation method --- .../opencv2/stitching/detail/seam_finders.hpp | 109 +++ modules/stitching/src/seam_finders.cpp | 850 +++++++++++++++++++++ samples/cpp/stitching_detailed.cpp | 8 +- 3 files changed, 966 insertions(+), 1 deletion(-) diff --git a/modules/stitching/include/opencv2/stitching/detail/seam_finders.hpp b/modules/stitching/include/opencv2/stitching/detail/seam_finders.hpp index bab29a0..009f853 100644 --- a/modules/stitching/include/opencv2/stitching/detail/seam_finders.hpp +++ b/modules/stitching/include/opencv2/stitching/detail/seam_finders.hpp @@ -43,6 +43,7 @@ #ifndef __OPENCV_STITCHING_SEAM_FINDERS_HPP__ #define __OPENCV_STITCHING_SEAM_FINDERS_HPP__ +#include #include "opencv2/core/core.hpp" #include "opencv2/opencv_modules.hpp" @@ -92,6 +93,114 @@ private: }; +class CV_EXPORTS DpSeamFinder : public SeamFinder +{ +public: + enum CostFunction { COLOR, COLOR_GRAD }; + + DpSeamFinder(CostFunction costFunc = COLOR_GRAD); + + CostFunction costFunction() const { return costFunc_; } + void setCostFunction(CostFunction val) { costFunc_ = val; } + +private: + enum ComponentState + { + FIRST = 1, SECOND = 2, INTERS = 4, + INTERS_FIRST = INTERS | FIRST, + INTERS_SECOND = INTERS | SECOND + }; + + class ImagePairLess + { + public: + ImagePairLess(const std::vector &images, const std::vector &corners) + : src_(&images[0]), corners_(&corners[0]) {} + + bool operator() (const std::pair &l, const std::pair &r) const + { + Point c1 = corners_[l.first] + Point(src_[l.first].cols / 2, src_[l.first].rows / 2); + Point c2 = corners_[l.second] + Point(src_[l.second].cols / 2, src_[l.second].rows / 2); + int d1 = (c1 - c2).dot(c1 - c2); + + c1 = corners_[r.first] + Point(src_[r.first].cols / 2, src_[r.first].rows / 2); + c2 = corners_[r.second] + Point(src_[r.second].cols / 2, src_[r.second].rows / 2); + int d2 = (c1 - c2).dot(c1 - c2); + + return d1 < d2; + } + + private: + const Mat *src_; + const Point *corners_; + }; + + class ClosePoints + { + public: + ClosePoints(int minDist) : minDist_(minDist) {} + + bool operator() (const Point &p1, const Point &p2) const + { + int dist2 = (p1.x-p2.x) * (p1.x-p2.x) + (p1.y-p2.y) * (p1.y-p2.y); + return dist2 < minDist_ * minDist_; + } + + private: + int minDist_; + }; + + virtual void find(const std::vector &src, const std::vector &corners, + std::vector &masks); + + void process(const Mat &image1, const Mat &image2, Point tl1, Point tl2, + Mat &mask1, Mat &mask2); + + void findComponents(); + + void findEdges(); + + void resolveConflicts(const Mat &image1, const Mat &image2, + Point tl1, Point tl2, Mat &mask1, Mat &mask2); + + void computeGradients(const Mat &image1, const Mat &image2); + + bool hasOnlyOneNeighbor(int c); + + bool closeToContour(int y, int x, const Mat_ &contourMask); + + bool getSeamTips(int c1, int c2, Point &p1, Point &p2); + + void computeCosts(const Mat &image1, const Mat &image2, Point tl1, Point tl2, + int c, Mat_ &costV, Mat_ &costH); + + bool estimateSeam( + const Mat &image1, const Mat &image2, Point tl1, Point tl2, int c, + Point p1, Point p2, std::vector &seam, bool &isHorizontal); + + void updateLabelsUsingSeam(int c1, int c2, const std::vector &seam, + bool isHorizontalSeam); + + CostFunction costFunc_; + + // processing images pair data + Point unionTl_, unionBr_; + Size unionSize_; + Mat_ mask1_, mask2_; + Mat_ contour1mask_, contour2mask_; + Mat_ gradx1_, grady1_; + Mat_ gradx2_, grady2_; + + // components data + int ncomps_; + Mat_ labels_; + std::vector states_; + std::vector tls_, brs_; + std::vector > contours_; + std::set > edges_; +}; + + class CV_EXPORTS GraphCutSeamFinderBase { public: diff --git a/modules/stitching/src/seam_finders.cpp b/modules/stitching/src/seam_finders.cpp index 2f6c78c..c155d07 100644 --- a/modules/stitching/src/seam_finders.cpp +++ b/modules/stitching/src/seam_finders.cpp @@ -41,6 +41,9 @@ //M*/ #include "precomp.hpp" +#include + +using namespace std; namespace cv { namespace detail { @@ -152,6 +155,853 @@ void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi) } +DpSeamFinder::DpSeamFinder(CostFunction costFunc) : costFunc_(costFunc) {} + + +void DpSeamFinder::find(const vector &src, const vector &corners, vector &masks) +{ + LOGLN("Finding seams..."); + int64 t = getTickCount(); + + if (src.size() == 0) + return; + + vector > pairs; + + for (size_t i = 0; i+1 < src.size(); ++i) + for (size_t j = i+1; j < src.size(); ++j) + pairs.push_back(make_pair(i, j)); + + sort(pairs.begin(), pairs.end(), ImagePairLess(src, corners)); + reverse(pairs.begin(), pairs.end()); + + for (size_t i = 0; i < pairs.size(); ++i) + { + int i0 = pairs[i].first, i1 = pairs[i].second; + process(src[i0], src[i1], corners[i0], corners[i1], masks[i0], masks[i1]); + } + + LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); +} + + +void DpSeamFinder::process( + const Mat &image1, const Mat &image2, Point tl1, Point tl2, + Mat &mask1, Mat &mask2) +{ + Point intersectTl(std::max(tl1.x, tl2.x), std::max(tl1.y, tl2.y)); + + Point intersectBr(std::min(tl1.x + image1.cols, tl2.x + image2.cols), + std::min(tl1.y + image1.rows, tl2.y + image2.rows)); + + if (intersectTl.x >= intersectBr.x || intersectTl.y >= intersectBr.y) + return; // there are no conflicts + + unionTl_ = Point(std::min(tl1.x, tl2.x), std::min(tl1.y, tl2.y)); + + unionBr_ = Point(std::max(tl1.x + image1.cols, tl2.x + image2.cols), + std::max(tl1.y + image1.rows, tl2.y + image2.rows)); + + unionSize_ = Size(unionBr_.x - unionTl_.x, unionBr_.y - unionTl_.y); + + mask1_ = Mat::zeros(unionSize_, CV_8U); + mask2_ = Mat::zeros(unionSize_, CV_8U); + + Mat tmp = mask1_(Rect(tl1.x - unionTl_.x, tl1.y - unionTl_.y, mask1.cols, mask1.rows)); + mask1.copyTo(tmp); + + tmp = mask2_(Rect(tl2.x - unionTl_.x, tl2.y - unionTl_.y, mask2.cols, mask2.rows)); + mask2.copyTo(tmp); + + // find both images contour masks + + contour1mask_ = Mat::zeros(unionSize_, CV_8U); + contour2mask_ = Mat::zeros(unionSize_, CV_8U); + + for (int y = 0; y < unionSize_.height; ++y) + { + for (int x = 0; x < unionSize_.width; ++x) + { + if (mask1_(y, x) && + ((x == 0 || !mask1_(y, x-1)) || (x == unionSize_.width-1 || !mask1_(y, x+1)) || + (y == 0 || !mask1_(y-1, x)) || (y == unionSize_.height-1 || !mask1_(y+1, x)))) + { + contour1mask_(y, x) = 255; + } + + if (mask2_(y, x) && + ((x == 0 || !mask2_(y, x-1)) || (x == unionSize_.width-1 || !mask2_(y, x+1)) || + (y == 0 || !mask2_(y-1, x)) || (y == unionSize_.height-1 || !mask2_(y+1, x)))) + { + contour2mask_(y, x) = 255; + } + } + } + + findComponents(); + + findEdges(); + + resolveConflicts(image1, image2, tl1, tl2, mask1, mask2); +} + + +void DpSeamFinder::findComponents() +{ + // label all connected components and get information about them + + ncomps_ = 0; + labels_.create(unionSize_); + states_.clear(); + tls_.clear(); + brs_.clear(); + contours_.clear(); + + for (int y = 0; y < unionSize_.height; ++y) + { + for (int x = 0; x < unionSize_.width; ++x) + { + if (mask1_(y, x) && mask2_(y, x)) + labels_(y, x) = numeric_limits::max(); + else if (mask1_(y, x)) + labels_(y, x) = numeric_limits::max()-1; + else if (mask2_(y, x)) + labels_(y, x) = numeric_limits::max()-2; + else + labels_(y, x) = 0; + } + } + + for (int y = 0; y < unionSize_.height; ++y) + { + for (int x = 0; x < unionSize_.width; ++x) + { + if (labels_(y, x) >= numeric_limits::max()-2) + { + if (labels_(y, x) == numeric_limits::max()) + states_.push_back(INTERS); + else if (labels_(y, x) == numeric_limits::max()-1) + states_.push_back(FIRST); + else if (labels_(y, x) == numeric_limits::max()-2) + states_.push_back(SECOND); + + floodFill(labels_, Point(x, y), ++ncomps_); + tls_.push_back(Point(x, y)); + brs_.push_back(Point(x+1, y+1)); + contours_.push_back(vector()); + } + + if (labels_(y, x)) + { + int l = labels_(y, x); + int ci = l-1; + + tls_[ci].x = std::min(tls_[ci].x, x); + tls_[ci].y = std::min(tls_[ci].y, y); + brs_[ci].x = std::max(brs_[ci].x, x+1); + brs_[ci].y = std::max(brs_[ci].y, y+1); + + if ((x == 0 || labels_(y, x-1) != l) || (x == unionSize_.width-1 || labels_(y, x+1) != l) || + (y == 0 || labels_(y-1, x) != l) || (y == unionSize_.height-1 || labels_(y+1, x) != l)) + { + contours_[ci].push_back(Point(x, y)); + } + } + } + } +} + + +void DpSeamFinder::findEdges() +{ + // find edges between components + + map, int> wedges; // weighted edges + + for (int ci = 0; ci < ncomps_-1; ++ci) + { + for (int cj = ci+1; cj < ncomps_; ++cj) + { + wedges[make_pair(ci, cj)] = 0; + wedges[make_pair(cj, ci)] = 0; + } + } + + for (int ci = 0; ci < ncomps_; ++ci) + { + for (size_t i = 0; i < contours_[ci].size(); ++i) + { + int x = contours_[ci][i].x; + int y = contours_[ci][i].y; + int l = ci + 1; + + if (x > 0 && labels_(y, x-1) && labels_(y, x-1) != l) + { + wedges[make_pair(ci, labels_(y, x-1)-1)]++; + wedges[make_pair(labels_(y, x-1)-1, ci)]++; + } + + if (y > 0 && labels_(y-1, x) && labels_(y-1, x) != l) + { + wedges[make_pair(ci, labels_(y-1, x)-1)]++; + wedges[make_pair(labels_(y-1, x)-1, ci)]++; + } + + if (x < unionSize_.width-1 && labels_(y, x+1) && labels_(y, x+1) != l) + { + wedges[make_pair(ci, labels_(y, x+1)-1)]++; + wedges[make_pair(labels_(y, x+1)-1, ci)]++; + } + + if (y < unionSize_.height-1 && labels_(y+1, x) && labels_(y+1, x) != l) + { + wedges[make_pair(ci, labels_(y+1, x)-1)]++; + wedges[make_pair(labels_(y+1, x)-1, ci)]++; + } + } + } + + edges_.clear(); + + for (int ci = 0; ci < ncomps_-1; ++ci) + { + for (int cj = ci+1; cj < ncomps_; ++cj) + { + map, int>::iterator itr = wedges.find(make_pair(ci, cj)); + if (itr != wedges.end() && itr->second > 0) + edges_.insert(itr->first); + + itr = wedges.find(make_pair(cj, ci)); + if (itr != wedges.end() && itr->second > 0) + edges_.insert(itr->first); + } + } +} + + +void DpSeamFinder::resolveConflicts(const Mat &image1, const Mat &image2, + Point tl1, Point tl2, Mat &mask1, Mat &mask2) +{ + if (costFunc_ == COLOR_GRAD) + computeGradients(image1, image2); + + // resolve conflicts between components + + bool hasConflict = true; + while (hasConflict) + { + int c1, c2; + hasConflict = false; + + for (set >::iterator itr = edges_.begin(); itr != edges_.end(); ++itr) + { + c1 = itr->first; + c2 = itr->second; + + if ((states_[c1] & INTERS) && (states_[c1] & (~INTERS)) != states_[c2]) + { + hasConflict = true; + break; + } + } + + if (hasConflict) + { + int l1 = c1+1, l2 = c2+1; + + if (hasOnlyOneNeighbor(c1)) + { + // if the first components has only one adjacent component + + for (int y = tls_[c1].y; y < brs_[c1].y; ++y) + for (int x = tls_[c1].x; x < brs_[c1].x; ++x) + if (labels_(y, x) == l1) + labels_(y, x) = l2; + + states_[c1] = states_[c2] == FIRST ? SECOND : FIRST; + } + else + { + // if the first component has more than one adjacent component + + Point p1, p2; + if (getSeamTips(c1, c2, p1, p2)) + { + vector seam; + bool isHorizontalSeam; + + if (estimateSeam(image1, image2, tl1, tl2, c1, p1, p2, seam, isHorizontalSeam)) + updateLabelsUsingSeam(c1, c2, seam, isHorizontalSeam); + } + + states_[c1] = states_[c2] == FIRST ? INTERS_SECOND : INTERS_FIRST; + } + + const int c[] = {c1, c2}; + const int l[] = {l1, l2}; + + for (int i = 0; i < 2; ++i) + { + // update information about the (i+1)-th component + + int x0 = tls_[c[i]].x, x1 = brs_[c[i]].x; + int y0 = tls_[c[i]].y, y1 = brs_[c[i]].y; + + tls_[c[i]] = Point(numeric_limits::max(), numeric_limits::max()); + brs_[c[i]] = Point(numeric_limits::min(), numeric_limits::min()); + contours_[c[i]].clear(); + + for (int y = y0; y < y1; ++y) + { + for (int x = x0; x < x1; ++x) + { + if (labels_(y, x) == l[i]) + { + tls_[c[i]].x = std::min(tls_[c[i]].x, x); + tls_[c[i]].y = std::min(tls_[c[i]].y, y); + brs_[c[i]].x = std::max(brs_[c[i]].x, x+1); + brs_[c[i]].y = std::max(brs_[c[i]].y, y+1); + + if ((x == 0 || labels_(y, x-1) != l[i]) || (x == unionSize_.width-1 || labels_(y, x+1) != l[i]) || + (y == 0 || labels_(y-1, x) != l[i]) || (y == unionSize_.height-1 || labels_(y+1, x) != l[i])) + { + contours_[c[i]].push_back(Point(x, y)); + } + } + } + } + } + + // remove edges + + edges_.erase(make_pair(c1, c2)); + edges_.erase(make_pair(c2, c1)); + } + } + + // update masks + + int dx1 = unionTl_.x - tl1.x, dy1 = unionTl_.y - tl1.y; + int dx2 = unionTl_.x - tl2.x, dy2 = unionTl_.y - tl2.y; + + for (int y = 0; y < mask2.rows; ++y) + { + for (int x = 0; x < mask2.cols; ++x) + { + int l = labels_(y - dy2, x - dx2); + if ((states_[l-1] & FIRST) && mask1.at(y - dy2 + dy1, x - dx2 + dx1)) + mask2.at(y, x) = 0; + } + } + + for (int y = 0; y < mask1.rows; ++y) + { + for (int x = 0; x < mask1.cols; ++x) + { + int l = labels_(y - dy1, x - dx1); + if ((states_[l-1] & SECOND) && mask2.at(y - dy1 + dy2, x - dx1 + dx2)) + mask1.at(y, x) = 0; + } + } +} + + +void DpSeamFinder::computeGradients(const Mat &image1, const Mat &image2) +{ + CV_Assert(costFunction() == COLOR_GRAD); + CV_Assert(image1.type() == CV_32FC3); + CV_Assert(image2.type() == CV_32FC3); + + Mat gray; + cvtColor(image1, gray, CV_BGR2GRAY); + Sobel(gray, gradx1_, CV_32F, 1, 0); + Sobel(gray, grady1_, CV_32F, 0, 1); + + cvtColor(image2, gray, CV_BGR2GRAY); + Sobel(gray, gradx2_, CV_32F, 1, 0); + Sobel(gray, grady2_, CV_32F, 0, 1); +} + + +bool DpSeamFinder::hasOnlyOneNeighbor(int c) +{ + set >::iterator begin, end; + begin = lower_bound(edges_.begin(), edges_.end(), make_pair(c, numeric_limits::min())); + end = upper_bound(edges_.begin(), edges_.end(), make_pair(c, numeric_limits::max())); + return ++begin == end; +} + + +bool DpSeamFinder::closeToContour(int y, int x, const Mat_ &contourMask) +{ + const int rad = 2; + + for (int dy = -rad; dy <= rad; ++dy) + { + if (y + dy >= 0 && y + dy < unionSize_.height) + { + for (int dx = -rad; dx <= rad; ++dx) + { + if (x + dx >= 0 && x + dx < unionSize_.width && + contourMask(y + dy, x + dx)) + { + return true; + } + } + } + } + + return false; +} + + +bool DpSeamFinder::getSeamTips(int c1, int c2, Point &p1, Point &p2) +{ + CV_Assert(states_[c1] & INTERS); + + // find special points + + vector specialPoints; + int l2 = c2+1; + + for (size_t i = 0; i < contours_[c1].size(); ++i) + { + int x = contours_[c1][i].x; + int y = contours_[c1][i].y; + + if (closeToContour(y, x, contour1mask_) && + closeToContour(y, x, contour2mask_) && + ((x > 0 && labels_(y, x-1) == l2) || + (y > 0 && labels_(y-1, x) == l2) || + (x < unionSize_.width-1 && labels_(y, x+1) == l2) || + (y < unionSize_.height-1 && labels_(y+1, x) == l2))) + { + specialPoints.push_back(Point(x, y)); + } + } + + if (specialPoints.size() < 2) + return false; + + // find clusters + + vector labels; + cv::partition(specialPoints, labels, ClosePoints(10)); + + int nlabels = *max_element(labels.begin(), labels.end()) + 1; + if (nlabels < 2) + return false; + + vector sum(nlabels); + vector > points(nlabels); + + for (size_t i = 0; i < specialPoints.size(); ++i) + { + sum[labels[i]] += specialPoints[i]; + points[labels[i]].push_back(specialPoints[i]); + } + + // select two most distant clusters + + int idx[2]; + + double maxDist = -numeric_limits::max(); + + for (int i = 0; i < nlabels-1; ++i) + { + for (int j = i+1; j < nlabels; ++j) + { + double size1 = points[i].size(), size2 = points[j].size(); + double cx1 = cvRound(sum[i].x / size1), cy1 = cvRound(sum[i].y / size1); + double cx2 = cvRound(sum[j].x / size2), cy2 = cvRound(sum[j].y / size1); + + double dist = (cx1 - cx2) * (cx1 - cx2) + (cy1 - cy2) * (cy1 - cy2); + if (dist > maxDist) + { + maxDist = dist; + idx[0] = i; + idx[1] = j; + } + } + } + + // select two points closest to the clusters' centers + + Point p[2]; + + for (int i = 0; i < 2; ++i) + { + double size = points[idx[i]].size(); + double cx = cvRound(sum[idx[i]].x / size); + double cy = cvRound(sum[idx[i]].y / size); + + int closest = -1; + double minDist = numeric_limits::max(); + + for (size_t j = 0; j < points[idx[i]].size(); ++j) + { + double dist = (points[idx[i]][j].x - cx) * (points[idx[i]][j].x - cx) + + (points[idx[i]][j].y - cy) * (points[idx[i]][j].y - cy); + if (dist < minDist) + { + minDist = dist; + closest = j; + } + } + + p[i] = points[idx[i]][closest]; + } + + p1 = p[0]; + p2 = p[1]; + return true; +} + + +void DpSeamFinder::computeCosts(const Mat &image1, const Mat &image2, Point tl1, Point tl2, + int c, Mat_ &costV, Mat_ &costH) +{ + CV_Assert(image1.type() == CV_32FC3); + CV_Assert(image2.type() == CV_32FC3); + CV_Assert(states_[c] & INTERS); + + // compute costs + + int l = c+1; + Rect roi(tls_[c], brs_[c]); + + int dx1 = unionTl_.x - tl1.x, dy1 = unionTl_.y - tl1.y; + int dx2 = unionTl_.x - tl2.x, dy2 = unionTl_.y - tl2.y; + + const float badRegionCost = normL2(Point3f(255.f, 255.f, 255.f), + Point3f(0.f, 0.f, 0.f)); + + costV.create(roi.height, roi.width+1); + + for (int y = roi.y; y < roi.br().y; ++y) + { + for (int x = roi.x; x < roi.br().x+1; ++x) + { + if (labels_(y, x) == l && x > 0 && labels_(y, x-1) == l) + { + const Point3f &pr1 = image1.at(y + dy1, x + dx1); + const Point3f &pl1 = image1.at(y + dy1, x + dx1 - 1); + const Point3f &pr2 = image2.at(y + dy2, x + dx2); + const Point3f &pl2 = image2.at(y + dy2, x + dx2 - 1); + + float costColor = (normL2(pl1, pr2) + normL2(pl2, pr1)) / 2; + if (costFunc_ == COLOR) + costV(y - roi.y, x - roi.x) = costColor; + else if (costFunc_ == COLOR_GRAD) + { + float costGrad = fabs(gradx1_(y + dy1, x + dx1)) + fabs(gradx1_(y + dy1, x + dx1 - 1)) + + fabs(gradx2_(y + dy2, x + dx2)) + fabs(gradx2_(y + dy2, x + dx2 - 1)) + 1.f; + costV(y - roi.y, x - roi.x) = costColor / costGrad; + } + } + else + costV(y - roi.y, x - roi.x) = badRegionCost; + } + } + + costH.create(roi.height+1, roi.width); + + for (int y = roi.y; y < roi.br().y+1; ++y) + { + for (int x = roi.x; x < roi.br().x; ++x) + { + if (labels_(y, x) == l && y > 0 && labels_(y-1, x) == l) + { + const Point3f &pd1 = image1.at(y + dy1, x + dx1); + const Point3f &pu1 = image1.at(y + dy1 - 1, x + dx1); + const Point3f &pd2 = image2.at(y + dy2, x + dx2); + const Point3f &pu2 = image2.at(y + dy2 - 1, x + dx2); + + float costColor = (normL2(pu1, pd2) + normL2(pu2, pd1)) / 2; + if (costFunc_ == COLOR) + costH(y - roi.y, x - roi.x) = costColor; + else if (costFunc_ == COLOR_GRAD) + { + float costGrad = fabs(grady1_(y + dy1, x + dx1)) + fabs(grady1_(y + dy1 - 1, x + dx1)) + + fabs(grady2_(y + dy2, x + dx2)) + fabs(grady2_(y + dy2 - 1, x + dx2)) + 1.f; + costH(y - roi.y, x - roi.x) = costColor / costGrad; + } + } + else + costH(y - roi.y, x - roi.x) = badRegionCost; + } + } +} + + +bool DpSeamFinder::estimateSeam( + const Mat &image1, const Mat &image2, Point tl1, Point tl2, int c, + Point p1, Point p2, vector &seam, bool &isHorizontal) +{ + CV_Assert(states_[c] & INTERS); + + Mat_ costV, costH; + computeCosts(image1, image2, tl1, tl2, c, costV, costH); + + Rect roi(tls_[c], brs_[c]); + Point src = p1 - roi.tl(); + Point dst = p2 - roi.tl(); + int l = c+1; + + // estimate seam direction + + bool swapped = false; + isHorizontal = std::abs(dst.x - src.x) > std::abs(dst.y - src.y); + + if (isHorizontal) + { + if (src.x > dst.x) + { + std::swap(src, dst); + swapped = true; + } + } + else if (src.y > dst.y) + { + swapped = true; + std::swap(src, dst); + } + + // find optimal control + + Mat_ control = Mat::zeros(roi.size(), CV_8U); + Mat_ reachable = Mat::zeros(roi.size(), CV_8U); + Mat_ cost = Mat::zeros(roi.size(), CV_32F); + + reachable(src) = 1; + cost(src) = 0.f; + + int nsteps; + pair steps[3]; + + if (isHorizontal) + { + for (int x = src.x+1; x <= dst.x; ++x) + { + for (int y = 0; y < roi.height; ++y) + { + // seam follows along upper side of pixels + + nsteps = 0; + + if (labels_(y + roi.y, x + roi.x) == l) + { + if (reachable(y, x-1)) + steps[nsteps++] = make_pair(cost(y, x-1) + costH(y, x-1), 1); + if (y > 0 && reachable(y-1, x-1)) + steps[nsteps++] = make_pair(cost(y-1, x-1) + costH(y-1, x-1) + costV(y-1, x), 2); + if (y < roi.height-1 && reachable(y+1, x-1)) + steps[nsteps++] = make_pair(cost(y+1, x-1) + costH(y+1, x-1) + costV(y, x), 3); + } + + if (nsteps) + { + pair opt = *min_element(steps, steps + nsteps); + cost(y, x) = opt.first; + control(y, x) = opt.second; + reachable(y, x) = 255; + } + } + } + } + else + { + for (int y = src.y+1; y <= dst.y; ++y) + { + for (int x = 0; x < roi.width; ++x) + { + // seam follows along left side of pixels + + nsteps = 0; + + if (labels_(y + roi.y, x + roi.x) == l) + { + if (reachable(y-1, x)) + steps[nsteps++] = make_pair(cost(y-1, x) + costV(y-1, x), 1); + if (x > 0 && reachable(y-1, x-1)) + steps[nsteps++] = make_pair(cost(y-1, x-1) + costV(y-1, x-1) + costH(y, x-1), 2); + if (x < roi.width-1 && reachable(y-1, x+1)) + steps[nsteps++] = make_pair(cost(y-1, x+1) + costV(y-1, x+1) + costH(y, x), 3); + } + + if (nsteps) + { + pair opt = *min_element(steps, steps + nsteps); + cost(y, x) = opt.first; + control(y, x) = opt.second; + reachable(y, x) = 255; + } + } + } + } + + if (!reachable(dst)) + return false; + + // restore seam + + Point p = dst; + seam.clear(); + seam.push_back(p + roi.tl()); + + if (isHorizontal) + { + for (; p.x != src.x; seam.push_back(p + roi.tl())) + { + if (control(p) == 2) p.y--; + else if (control(p) == 3) p.y++; + p.x--; + } + } + else + { + for (; p.y != src.y; seam.push_back(p + roi.tl())) + { + if (control(p) == 2) p.x--; + else if (control(p) == 3) p.x++; + p.y--; + } + } + + if (!swapped) + reverse(seam.begin(), seam.end()); + + CV_Assert(seam.front() == p1); + CV_Assert(seam.back() == p2); + return true; +} + + +void DpSeamFinder::updateLabelsUsingSeam(int c1, int c2, const vector &seam, bool isHorizontalSeam) +{ + Mat_ mask = Mat::zeros(brs_[c1].y - tls_[c1].y, brs_[c1].x - tls_[c1].x, CV_32S); + + for (size_t i = 0; i < contours_[c1].size(); ++i) + mask(contours_[c1][i] - tls_[c1]) = 255; + + for (size_t i = 0; i < seam.size(); ++i) + mask(seam[i] - tls_[c1]) = 255; + + // find connected components after seam carving + + int l1 = c1+1, l2 = c2+1; + + int ncomps = 0; + + for (int y = 0; y < mask.rows; ++y) + for (int x = 0; x < mask.cols; ++x) + if (!mask(y, x) && labels_(y + tls_[c1].y, x + tls_[c1].x) == l1) + floodFill(mask, Point(x, y), ++ncomps); + + for (size_t i = 0; i < contours_[c1].size(); ++i) + { + int x = contours_[c1][i].x - tls_[c1].x; + int y = contours_[c1][i].y - tls_[c1].y; + + bool ok = false; + static const int dx[] = {-1, +1, 0, 0, -1, +1, -1, +1}; + static const int dy[] = {0, 0, -1, +1, -1, -1, +1, +1}; + + for (int j = 0; j < 8; ++j) + { + int c = x + dx[j]; + int r = y + dy[j]; + + if (c > 0 && c < mask.cols && r > 0 && r < mask.rows && + mask(r, c) && mask(r, c) != 255) + { + ok = true; + mask(y, x) = mask(r, c); + } + } + + if (!ok) + mask(y, x) = 0; + } + + if (isHorizontalSeam) + { + for (size_t i = 0; i < seam.size(); ++i) + { + int x = seam[i].x - tls_[c1].x; + int y = seam[i].y - tls_[c1].y; + + if (y < mask.rows-1 && mask(y+1, x) && mask(y+1, x) != 255) + mask(y, x) = mask(y+1, x); + else + mask(y, x) = 0; + } + } + else + { + for (size_t i = 0; i < seam.size(); ++i) + { + int x = seam[i].x - tls_[c1].x; + int y = seam[i].y - tls_[c1].y; + + if (x < mask.cols-1 && mask(y, x+1) && mask(y, x+1) != 255) + mask(y, x) = mask(y, x+1); + else + mask(y, x) = 0; + } + } + + // find new components connected with the second component and + // with other components except the ones we are working with + + map connect2; + map connectOther; + + for (int i = 1; i <= ncomps; ++i) + { + connect2.insert(make_pair(i, 0)); + connectOther.insert(make_pair(i, 0)); + } + + for (size_t i = 0; i < contours_[c1].size(); ++i) + { + int x = contours_[c1][i].x; + int y = contours_[c1][i].y; + + if ((x > 0 && labels_(y, x-1) == l2) || + (y > 0 && labels_(y-1, x) == l2) || + (x < unionSize_.width-1 && labels_(y, x+1) == l2) || + (y < unionSize_.height-1 && labels_(y+1, x) == l2)) + { + connect2[mask(y - tls_[c1].y, x - tls_[c1].x)]++; + } + + if ((x > 0 && labels_(y, x-1) != l1 && labels_(y, x-1) != l2) || + (y > 0 && labels_(y-1, x) != l1 && labels_(y-1, x) != l2) || + (x < unionSize_.width-1 && labels_(y, x+1) != l1 && labels_(y, x+1) != l2) || + (y < unionSize_.height-1 && labels_(y+1, x) != l1 && labels_(y+1, x) != l2)) + { + connectOther[mask(y - tls_[c1].y, x - tls_[c1].x)]++; + } + } + + vector isAdjComp(ncomps + 1, 0); + + for (map::iterator itr = connect2.begin(); itr != connect2.end(); ++itr) + { + double len = contours_[c1].size(); + isAdjComp[itr->first] = itr->second / len > 0.05 && connectOther.find(itr->first)->second / len < 0.1; + } + + // update labels + + for (int y = 0; y < mask.rows; ++y) + for (int x = 0; x < mask.cols; ++x) + if (mask(y, x) && isAdjComp[mask(y, x)]) + labels_(y + tls_[c1].y, x + tls_[c1].x) = l2; +} + + class GraphCutSeamFinder::Impl : public PairwiseSeamFinder { public: diff --git a/samples/cpp/stitching_detailed.cpp b/samples/cpp/stitching_detailed.cpp index 220d466..1e46117 100644 --- a/samples/cpp/stitching_detailed.cpp +++ b/samples/cpp/stitching_detailed.cpp @@ -276,7 +276,9 @@ static int parseCmdArgs(int argc, char** argv) if (string(argv[i + 1]) == "no" || string(argv[i + 1]) == "voronoi" || string(argv[i + 1]) == "gc_color" || - string(argv[i + 1]) == "gc_colorgrad") + string(argv[i + 1]) == "gc_colorgrad" || + string(argv[i + 1]) == "dp_color" || + string(argv[i + 1]) == "dp_colorgrad") seam_find_type = argv[i + 1]; else { @@ -612,6 +614,10 @@ int main(int argc, char* argv[]) #endif seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR_GRAD); } + else if (seam_find_type == "dp_color") + seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR); + else if (seam_find_type == "dp_colorgrad") + seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR_GRAD); if (seam_finder.empty()) { cout << "Can't create the following seam finder '" << seam_find_type << "'\n"; -- 2.7.4