1/*M/////////////////////////////////////////////////////////////////////////////////////// 2// 3// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4// 5// By downloading, copying, installing or using the software you agree to this license. 6// If you do not agree to this license, do not download, install, 7// copy or use the software. 8// 9// 10// License Agreement 11// For Open Source Computer Vision Library 12// 13// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. 14// Copyright (C) 2009, Willow Garage Inc., all rights reserved. 15// Third party copyrights are property of their respective owners. 16// 17// Redistribution and use in source and binary forms, with or without modification, 18// are permitted provided that the following conditions are met: 19// 20// * Redistribution's of source code must retain the above copyright notice, 21// this list of conditions and the following disclaimer. 22// 23// * Redistribution's in binary form must reproduce the above copyright notice, 24// this list of conditions and the following disclaimer in the documentation 25// and/or other materials provided with the distribution. 26// 27// * The name of the copyright holders may not be used to endorse or promote products 28// derived from this software without specific prior written permission. 29// 30// This software is provided by the copyright holders and contributors "as is" and 31// any express or implied warranties, including, but not limited to, the implied 32// warranties of merchantability and fitness for a particular purpose are disclaimed. 33// In no event shall the Intel Corporation or contributors be liable for any direct, 34// indirect, incidental, special, exemplary, or consequential damages 35// (including, but not limited to, procurement of substitute goods or services; 36// loss of use, data, or profits; or business interruption) however caused 37// and on any theory of liability, whether in contract, strict liability, 38// or tort (including negligence or otherwise) arising in any way out of 39// the use of this software, even if advised of the possibility of such damage. 40// 41//M*/ 42 43#include "precomp.hpp" 44#include <map> 45 46namespace cv { 47namespace detail { 48 49void PairwiseSeamFinder::find(const std::vector<UMat> &src, const std::vector<Point> &corners, 50 std::vector<UMat> &masks) 51{ 52 LOGLN("Finding seams..."); 53 if (src.size() == 0) 54 return; 55 56#if ENABLE_LOG 57 int64 t = getTickCount(); 58#endif 59 60 images_ = src; 61 sizes_.resize(src.size()); 62 for (size_t i = 0; i < src.size(); ++i) 63 sizes_[i] = src[i].size(); 64 corners_ = corners; 65 masks_ = masks; 66 run(); 67 68 LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); 69} 70 71 72void PairwiseSeamFinder::run() 73{ 74 for (size_t i = 0; i < sizes_.size() - 1; ++i) 75 { 76 for (size_t j = i + 1; j < sizes_.size(); ++j) 77 { 78 Rect roi; 79 if (overlapRoi(corners_[i], corners_[j], sizes_[i], sizes_[j], roi)) 80 findInPair(i, j, roi); 81 } 82 } 83} 84 85void VoronoiSeamFinder::find(const std::vector<UMat> &src, const std::vector<Point> &corners, 86 std::vector<UMat> &masks) 87{ 88 PairwiseSeamFinder::find(src, corners, masks); 89} 90 91void VoronoiSeamFinder::find(const std::vector<Size> &sizes, const std::vector<Point> &corners, 92 std::vector<UMat> &masks) 93{ 94 LOGLN("Finding seams..."); 95 if (sizes.size() == 0) 96 return; 97 98#if ENABLE_LOG 99 int64 t = getTickCount(); 100#endif 101 102 sizes_ = sizes; 103 corners_ = corners; 104 masks_ = masks; 105 run(); 106 107 LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); 108} 109 110 111void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi) 112{ 113 const int gap = 10; 114 Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); 115 Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); 116 117 Size img1 = sizes_[first], img2 = sizes_[second]; 118 Mat mask1 = masks_[first].getMat(ACCESS_READ), mask2 = masks_[second].getMat(ACCESS_READ); 119 Point tl1 = corners_[first], tl2 = corners_[second]; 120 121 // Cut submasks with some gap 122 for (int y = -gap; y < roi.height + gap; ++y) 123 { 124 for (int x = -gap; x < roi.width + gap; ++x) 125 { 126 int y1 = roi.y - tl1.y + y; 127 int x1 = roi.x - tl1.x + x; 128 if (y1 >= 0 && x1 >= 0 && y1 < img1.height && x1 < img1.width) 129 submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1); 130 else 131 submask1.at<uchar>(y + gap, x + gap) = 0; 132 133 int y2 = roi.y - tl2.y + y; 134 int x2 = roi.x - tl2.x + x; 135 if (y2 >= 0 && x2 >= 0 && y2 < img2.height && x2 < img2.width) 136 submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2); 137 else 138 submask2.at<uchar>(y + gap, x + gap) = 0; 139 } 140 } 141 142 Mat collision = (submask1 != 0) & (submask2 != 0); 143 Mat unique1 = submask1.clone(); unique1.setTo(0, collision); 144 Mat unique2 = submask2.clone(); unique2.setTo(0, collision); 145 146 Mat dist1, dist2; 147 distanceTransform(unique1 == 0, dist1, DIST_L1, 3); 148 distanceTransform(unique2 == 0, dist2, DIST_L1, 3); 149 150 Mat seam = dist1 < dist2; 151 152 for (int y = 0; y < roi.height; ++y) 153 { 154 for (int x = 0; x < roi.width; ++x) 155 { 156 if (seam.at<uchar>(y + gap, x + gap)) 157 mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; 158 else 159 mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; 160 } 161 } 162} 163 164 165DpSeamFinder::DpSeamFinder(CostFunction costFunc) : costFunc_(costFunc) {} 166 167 168void DpSeamFinder::find(const std::vector<UMat> &src, const std::vector<Point> &corners, std::vector<UMat> &masks) 169{ 170 LOGLN("Finding seams..."); 171#if ENABLE_LOG 172 int64 t = getTickCount(); 173#endif 174 175 if (src.size() == 0) 176 return; 177 178 std::vector<std::pair<size_t, size_t> > pairs; 179 180 for (size_t i = 0; i+1 < src.size(); ++i) 181 for (size_t j = i+1; j < src.size(); ++j) 182 pairs.push_back(std::make_pair(i, j)); 183 184 { 185 std::vector<Mat> _src(src.size()); 186 for (size_t i = 0; i < src.size(); ++i) _src[i] = src[i].getMat(ACCESS_READ); 187 sort(pairs.begin(), pairs.end(), ImagePairLess(_src, corners)); 188 } 189 std::reverse(pairs.begin(), pairs.end()); 190 191 for (size_t i = 0; i < pairs.size(); ++i) 192 { 193 size_t i0 = pairs[i].first, i1 = pairs[i].second; 194 Mat mask0 = masks[i0].getMat(ACCESS_RW), mask1 = masks[i1].getMat(ACCESS_RW); 195 process(src[i0].getMat(ACCESS_READ), src[i1].getMat(ACCESS_READ), corners[i0], corners[i1], mask0, mask1); 196 } 197 198 LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); 199} 200 201 202void DpSeamFinder::process( 203 const Mat &image1, const Mat &image2, Point tl1, Point tl2, 204 Mat &mask1, Mat &mask2) 205{ 206 CV_Assert(image1.size() == mask1.size()); 207 CV_Assert(image2.size() == mask2.size()); 208 209 Point intersectTl(std::max(tl1.x, tl2.x), std::max(tl1.y, tl2.y)); 210 211 Point intersectBr(std::min(tl1.x + image1.cols, tl2.x + image2.cols), 212 std::min(tl1.y + image1.rows, tl2.y + image2.rows)); 213 214 if (intersectTl.x >= intersectBr.x || intersectTl.y >= intersectBr.y) 215 return; // there are no conflicts 216 217 unionTl_ = Point(std::min(tl1.x, tl2.x), std::min(tl1.y, tl2.y)); 218 219 unionBr_ = Point(std::max(tl1.x + image1.cols, tl2.x + image2.cols), 220 std::max(tl1.y + image1.rows, tl2.y + image2.rows)); 221 222 unionSize_ = Size(unionBr_.x - unionTl_.x, unionBr_.y - unionTl_.y); 223 224 mask1_ = Mat::zeros(unionSize_, CV_8U); 225 mask2_ = Mat::zeros(unionSize_, CV_8U); 226 227 Mat tmp = mask1_(Rect(tl1.x - unionTl_.x, tl1.y - unionTl_.y, mask1.cols, mask1.rows)); 228 mask1.copyTo(tmp); 229 230 tmp = mask2_(Rect(tl2.x - unionTl_.x, tl2.y - unionTl_.y, mask2.cols, mask2.rows)); 231 mask2.copyTo(tmp); 232 233 // find both images contour masks 234 235 contour1mask_ = Mat::zeros(unionSize_, CV_8U); 236 contour2mask_ = Mat::zeros(unionSize_, CV_8U); 237 238 for (int y = 0; y < unionSize_.height; ++y) 239 { 240 for (int x = 0; x < unionSize_.width; ++x) 241 { 242 if (mask1_(y, x) && 243 ((x == 0 || !mask1_(y, x-1)) || (x == unionSize_.width-1 || !mask1_(y, x+1)) || 244 (y == 0 || !mask1_(y-1, x)) || (y == unionSize_.height-1 || !mask1_(y+1, x)))) 245 { 246 contour1mask_(y, x) = 255; 247 } 248 249 if (mask2_(y, x) && 250 ((x == 0 || !mask2_(y, x-1)) || (x == unionSize_.width-1 || !mask2_(y, x+1)) || 251 (y == 0 || !mask2_(y-1, x)) || (y == unionSize_.height-1 || !mask2_(y+1, x)))) 252 { 253 contour2mask_(y, x) = 255; 254 } 255 } 256 } 257 258 findComponents(); 259 260 findEdges(); 261 262 resolveConflicts(image1, image2, tl1, tl2, mask1, mask2); 263} 264 265 266void DpSeamFinder::findComponents() 267{ 268 // label all connected components and get information about them 269 270 ncomps_ = 0; 271 labels_.create(unionSize_); 272 states_.clear(); 273 tls_.clear(); 274 brs_.clear(); 275 contours_.clear(); 276 277 for (int y = 0; y < unionSize_.height; ++y) 278 { 279 for (int x = 0; x < unionSize_.width; ++x) 280 { 281 if (mask1_(y, x) && mask2_(y, x)) 282 labels_(y, x) = std::numeric_limits<int>::max(); 283 else if (mask1_(y, x)) 284 labels_(y, x) = std::numeric_limits<int>::max()-1; 285 else if (mask2_(y, x)) 286 labels_(y, x) = std::numeric_limits<int>::max()-2; 287 else 288 labels_(y, x) = 0; 289 } 290 } 291 292 for (int y = 0; y < unionSize_.height; ++y) 293 { 294 for (int x = 0; x < unionSize_.width; ++x) 295 { 296 if (labels_(y, x) >= std::numeric_limits<int>::max()-2) 297 { 298 if (labels_(y, x) == std::numeric_limits<int>::max()) 299 states_.push_back(INTERS); 300 else if (labels_(y, x) == std::numeric_limits<int>::max()-1) 301 states_.push_back(FIRST); 302 else if (labels_(y, x) == std::numeric_limits<int>::max()-2) 303 states_.push_back(SECOND); 304 305 floodFill(labels_, Point(x, y), ++ncomps_); 306 tls_.push_back(Point(x, y)); 307 brs_.push_back(Point(x+1, y+1)); 308 contours_.push_back(std::vector<Point>()); 309 } 310 311 if (labels_(y, x)) 312 { 313 int l = labels_(y, x); 314 int ci = l-1; 315 316 tls_[ci].x = std::min(tls_[ci].x, x); 317 tls_[ci].y = std::min(tls_[ci].y, y); 318 brs_[ci].x = std::max(brs_[ci].x, x+1); 319 brs_[ci].y = std::max(brs_[ci].y, y+1); 320 321 if ((x == 0 || labels_(y, x-1) != l) || (x == unionSize_.width-1 || labels_(y, x+1) != l) || 322 (y == 0 || labels_(y-1, x) != l) || (y == unionSize_.height-1 || labels_(y+1, x) != l)) 323 { 324 contours_[ci].push_back(Point(x, y)); 325 } 326 } 327 } 328 } 329} 330 331 332void DpSeamFinder::findEdges() 333{ 334 // find edges between components 335 336 std::map<std::pair<int, int>, int> wedges; // weighted edges 337 338 for (int ci = 0; ci < ncomps_-1; ++ci) 339 { 340 for (int cj = ci+1; cj < ncomps_; ++cj) 341 { 342 wedges[std::make_pair(ci, cj)] = 0; 343 wedges[std::make_pair(cj, ci)] = 0; 344 } 345 } 346 347 for (int ci = 0; ci < ncomps_; ++ci) 348 { 349 for (size_t i = 0; i < contours_[ci].size(); ++i) 350 { 351 int x = contours_[ci][i].x; 352 int y = contours_[ci][i].y; 353 int l = ci + 1; 354 355 if (x > 0 && labels_(y, x-1) && labels_(y, x-1) != l) 356 { 357 wedges[std::make_pair(ci, labels_(y, x-1)-1)]++; 358 wedges[std::make_pair(labels_(y, x-1)-1, ci)]++; 359 } 360 361 if (y > 0 && labels_(y-1, x) && labels_(y-1, x) != l) 362 { 363 wedges[std::make_pair(ci, labels_(y-1, x)-1)]++; 364 wedges[std::make_pair(labels_(y-1, x)-1, ci)]++; 365 } 366 367 if (x < unionSize_.width-1 && labels_(y, x+1) && labels_(y, x+1) != l) 368 { 369 wedges[std::make_pair(ci, labels_(y, x+1)-1)]++; 370 wedges[std::make_pair(labels_(y, x+1)-1, ci)]++; 371 } 372 373 if (y < unionSize_.height-1 && labels_(y+1, x) && labels_(y+1, x) != l) 374 { 375 wedges[std::make_pair(ci, labels_(y+1, x)-1)]++; 376 wedges[std::make_pair(labels_(y+1, x)-1, ci)]++; 377 } 378 } 379 } 380 381 edges_.clear(); 382 383 for (int ci = 0; ci < ncomps_-1; ++ci) 384 { 385 for (int cj = ci+1; cj < ncomps_; ++cj) 386 { 387 std::map<std::pair<int, int>, int>::iterator itr = wedges.find(std::make_pair(ci, cj)); 388 if (itr != wedges.end() && itr->second > 0) 389 edges_.insert(itr->first); 390 391 itr = wedges.find(std::make_pair(cj, ci)); 392 if (itr != wedges.end() && itr->second > 0) 393 edges_.insert(itr->first); 394 } 395 } 396} 397 398 399void DpSeamFinder::resolveConflicts( 400 const Mat &image1, const Mat &image2, Point tl1, Point tl2, Mat &mask1, Mat &mask2) 401{ 402 if (costFunc_ == COLOR_GRAD) 403 computeGradients(image1, image2); 404 405 // resolve conflicts between components 406 407 bool hasConflict = true; 408 while (hasConflict) 409 { 410 int c1 = 0, c2 = 0; 411 hasConflict = false; 412 413 for (std::set<std::pair<int, int> >::iterator itr = edges_.begin(); itr != edges_.end(); ++itr) 414 { 415 c1 = itr->first; 416 c2 = itr->second; 417 418 if ((states_[c1] & INTERS) && (states_[c1] & (~INTERS)) != states_[c2]) 419 { 420 hasConflict = true; 421 break; 422 } 423 } 424 425 if (hasConflict) 426 { 427 int l1 = c1+1, l2 = c2+1; 428 429 if (hasOnlyOneNeighbor(c1)) 430 { 431 // if the first components has only one adjacent component 432 433 for (int y = tls_[c1].y; y < brs_[c1].y; ++y) 434 for (int x = tls_[c1].x; x < brs_[c1].x; ++x) 435 if (labels_(y, x) == l1) 436 labels_(y, x) = l2; 437 438 states_[c1] = states_[c2] == FIRST ? SECOND : FIRST; 439 } 440 else 441 { 442 // if the first component has more than one adjacent component 443 444 Point p1, p2; 445 if (getSeamTips(c1, c2, p1, p2)) 446 { 447 std::vector<Point> seam; 448 bool isHorizontalSeam; 449 450 if (estimateSeam(image1, image2, tl1, tl2, c1, p1, p2, seam, isHorizontalSeam)) 451 updateLabelsUsingSeam(c1, c2, seam, isHorizontalSeam); 452 } 453 454 states_[c1] = states_[c2] == FIRST ? INTERS_SECOND : INTERS_FIRST; 455 } 456 457 const int c[] = {c1, c2}; 458 const int l[] = {l1, l2}; 459 460 for (int i = 0; i < 2; ++i) 461 { 462 // update information about the (i+1)-th component 463 464 int x0 = tls_[c[i]].x, x1 = brs_[c[i]].x; 465 int y0 = tls_[c[i]].y, y1 = brs_[c[i]].y; 466 467 tls_[c[i]] = Point(std::numeric_limits<int>::max(), std::numeric_limits<int>::max()); 468 brs_[c[i]] = Point(std::numeric_limits<int>::min(), std::numeric_limits<int>::min()); 469 contours_[c[i]].clear(); 470 471 for (int y = y0; y < y1; ++y) 472 { 473 for (int x = x0; x < x1; ++x) 474 { 475 if (labels_(y, x) == l[i]) 476 { 477 tls_[c[i]].x = std::min(tls_[c[i]].x, x); 478 tls_[c[i]].y = std::min(tls_[c[i]].y, y); 479 brs_[c[i]].x = std::max(brs_[c[i]].x, x+1); 480 brs_[c[i]].y = std::max(brs_[c[i]].y, y+1); 481 482 if ((x == 0 || labels_(y, x-1) != l[i]) || (x == unionSize_.width-1 || labels_(y, x+1) != l[i]) || 483 (y == 0 || labels_(y-1, x) != l[i]) || (y == unionSize_.height-1 || labels_(y+1, x) != l[i])) 484 { 485 contours_[c[i]].push_back(Point(x, y)); 486 } 487 } 488 } 489 } 490 } 491 492 // remove edges 493 494 edges_.erase(std::make_pair(c1, c2)); 495 edges_.erase(std::make_pair(c2, c1)); 496 } 497 } 498 499 // update masks 500 501 int dx1 = unionTl_.x - tl1.x, dy1 = unionTl_.y - tl1.y; 502 int dx2 = unionTl_.x - tl2.x, dy2 = unionTl_.y - tl2.y; 503 504 for (int y = 0; y < mask2.rows; ++y) 505 { 506 for (int x = 0; x < mask2.cols; ++x) 507 { 508 int l = labels_(y - dy2, x - dx2); 509 if (l > 0 && (states_[l-1] & FIRST) && mask1.at<uchar>(y - dy2 + dy1, x - dx2 + dx1)) 510 mask2.at<uchar>(y, x) = 0; 511 } 512 } 513 514 for (int y = 0; y < mask1.rows; ++y) 515 { 516 for (int x = 0; x < mask1.cols; ++x) 517 { 518 int l = labels_(y - dy1, x - dx1); 519 if (l > 0 && (states_[l-1] & SECOND) && mask2.at<uchar>(y - dy1 + dy2, x - dx1 + dx2)) 520 mask1.at<uchar>(y, x) = 0; 521 } 522 } 523} 524 525 526void DpSeamFinder::computeGradients(const Mat &image1, const Mat &image2) 527{ 528 CV_Assert(image1.channels() == 3 || image1.channels() == 4); 529 CV_Assert(image2.channels() == 3 || image2.channels() == 4); 530 CV_Assert(costFunction() == COLOR_GRAD); 531 532 Mat gray; 533 534 if (image1.channels() == 3) 535 cvtColor(image1, gray, COLOR_BGR2GRAY); 536 else if (image1.channels() == 4) 537 cvtColor(image1, gray, COLOR_BGRA2GRAY); 538 539 Sobel(gray, gradx1_, CV_32F, 1, 0); 540 Sobel(gray, grady1_, CV_32F, 0, 1); 541 542 if (image2.channels() == 3) 543 cvtColor(image2, gray, COLOR_BGR2GRAY); 544 else if (image2.channels() == 4) 545 cvtColor(image2, gray, COLOR_BGRA2GRAY); 546 547 Sobel(gray, gradx2_, CV_32F, 1, 0); 548 Sobel(gray, grady2_, CV_32F, 0, 1); 549} 550 551 552bool DpSeamFinder::hasOnlyOneNeighbor(int comp) 553{ 554 std::set<std::pair<int, int> >::iterator begin, end; 555 begin = lower_bound(edges_.begin(), edges_.end(), std::make_pair(comp, std::numeric_limits<int>::min())); 556 end = upper_bound(edges_.begin(), edges_.end(), std::make_pair(comp, std::numeric_limits<int>::max())); 557 return ++begin == end; 558} 559 560 561bool DpSeamFinder::closeToContour(int y, int x, const Mat_<uchar> &contourMask) 562{ 563 const int rad = 2; 564 565 for (int dy = -rad; dy <= rad; ++dy) 566 { 567 if (y + dy >= 0 && y + dy < unionSize_.height) 568 { 569 for (int dx = -rad; dx <= rad; ++dx) 570 { 571 if (x + dx >= 0 && x + dx < unionSize_.width && 572 contourMask(y + dy, x + dx)) 573 { 574 return true; 575 } 576 } 577 } 578 } 579 580 return false; 581} 582 583 584bool DpSeamFinder::getSeamTips(int comp1, int comp2, Point &p1, Point &p2) 585{ 586 CV_Assert(states_[comp1] & INTERS); 587 588 // find special points 589 590 std::vector<Point> specialPoints; 591 int l2 = comp2+1; 592 593 for (size_t i = 0; i < contours_[comp1].size(); ++i) 594 { 595 int x = contours_[comp1][i].x; 596 int y = contours_[comp1][i].y; 597 598 if (closeToContour(y, x, contour1mask_) && 599 closeToContour(y, x, contour2mask_) && 600 ((x > 0 && labels_(y, x-1) == l2) || 601 (y > 0 && labels_(y-1, x) == l2) || 602 (x < unionSize_.width-1 && labels_(y, x+1) == l2) || 603 (y < unionSize_.height-1 && labels_(y+1, x) == l2))) 604 { 605 specialPoints.push_back(Point(x, y)); 606 } 607 } 608 609 if (specialPoints.size() < 2) 610 return false; 611 612 // find clusters 613 614 std::vector<int> labels; 615 cv::partition(specialPoints, labels, ClosePoints(10)); 616 617 int nlabels = *std::max_element(labels.begin(), labels.end()) + 1; 618 if (nlabels < 2) 619 return false; 620 621 std::vector<Point> sum(nlabels); 622 std::vector<std::vector<Point> > points(nlabels); 623 624 for (size_t i = 0; i < specialPoints.size(); ++i) 625 { 626 sum[labels[i]] += specialPoints[i]; 627 points[labels[i]].push_back(specialPoints[i]); 628 } 629 630 // select two most distant clusters 631 632 int idx[2] = {-1,-1}; 633 double maxDist = -std::numeric_limits<double>::max(); 634 635 for (int i = 0; i < nlabels-1; ++i) 636 { 637 for (int j = i+1; j < nlabels; ++j) 638 { 639 double size1 = static_cast<double>(points[i].size()), size2 = static_cast<double>(points[j].size()); 640 double cx1 = cvRound(sum[i].x / size1), cy1 = cvRound(sum[i].y / size1); 641 double cx2 = cvRound(sum[j].x / size2), cy2 = cvRound(sum[j].y / size1); 642 643 double dist = (cx1 - cx2) * (cx1 - cx2) + (cy1 - cy2) * (cy1 - cy2); 644 if (dist > maxDist) 645 { 646 maxDist = dist; 647 idx[0] = i; 648 idx[1] = j; 649 } 650 } 651 } 652 653 // select two points closest to the clusters' centers 654 655 Point p[2]; 656 657 for (int i = 0; i < 2; ++i) 658 { 659 double size = static_cast<double>(points[idx[i]].size()); 660 double cx = cvRound(sum[idx[i]].x / size); 661 double cy = cvRound(sum[idx[i]].y / size); 662 663 size_t closest = points[idx[i]].size(); 664 double minDist = std::numeric_limits<double>::max(); 665 666 for (size_t j = 0; j < points[idx[i]].size(); ++j) 667 { 668 double dist = (points[idx[i]][j].x - cx) * (points[idx[i]][j].x - cx) + 669 (points[idx[i]][j].y - cy) * (points[idx[i]][j].y - cy); 670 if (dist < minDist) 671 { 672 minDist = dist; 673 closest = j; 674 } 675 } 676 677 p[i] = points[idx[i]][closest]; 678 } 679 680 p1 = p[0]; 681 p2 = p[1]; 682 return true; 683} 684 685 686namespace 687{ 688 689template <typename T> 690float diffL2Square3(const Mat &image1, int y1, int x1, const Mat &image2, int y2, int x2) 691{ 692 const T *r1 = image1.ptr<T>(y1); 693 const T *r2 = image2.ptr<T>(y2); 694 return static_cast<float>(sqr(r1[3*x1] - r2[3*x2]) + sqr(r1[3*x1+1] - r2[3*x2+1]) + 695 sqr(r1[3*x1+2] - r2[3*x2+2])); 696} 697 698 699template <typename T> 700float diffL2Square4(const Mat &image1, int y1, int x1, const Mat &image2, int y2, int x2) 701{ 702 const T *r1 = image1.ptr<T>(y1); 703 const T *r2 = image2.ptr<T>(y2); 704 return static_cast<float>(sqr(r1[4*x1] - r2[4*x2]) + sqr(r1[4*x1+1] - r2[4*x2+1]) + 705 sqr(r1[4*x1+2] - r2[4*x2+2])); 706} 707 708} // namespace 709 710 711void DpSeamFinder::computeCosts( 712 const Mat &image1, const Mat &image2, Point tl1, Point tl2, 713 int comp, Mat_<float> &costV, Mat_<float> &costH) 714{ 715 CV_Assert(states_[comp] & INTERS); 716 717 // compute costs 718 719 float (*diff)(const Mat&, int, int, const Mat&, int, int) = 0; 720 if (image1.type() == CV_32FC3 && image2.type() == CV_32FC3) 721 diff = diffL2Square3<float>; 722 else if (image1.type() == CV_8UC3 && image2.type() == CV_8UC3) 723 diff = diffL2Square3<uchar>; 724 else if (image1.type() == CV_32FC4 && image2.type() == CV_32FC4) 725 diff = diffL2Square4<float>; 726 else if (image1.type() == CV_8UC4 && image2.type() == CV_8UC4) 727 diff = diffL2Square4<uchar>; 728 else 729 CV_Error(Error::StsBadArg, "both images must have CV_32FC3(4) or CV_8UC3(4) type"); 730 731 int l = comp+1; 732 Rect roi(tls_[comp], brs_[comp]); 733 734 int dx1 = unionTl_.x - tl1.x, dy1 = unionTl_.y - tl1.y; 735 int dx2 = unionTl_.x - tl2.x, dy2 = unionTl_.y - tl2.y; 736 737 const float badRegionCost = normL2(Point3f(255.f, 255.f, 255.f), 738 Point3f(0.f, 0.f, 0.f)); 739 740 costV.create(roi.height, roi.width+1); 741 742 for (int y = roi.y; y < roi.br().y; ++y) 743 { 744 for (int x = roi.x; x < roi.br().x+1; ++x) 745 { 746 if (labels_(y, x) == l && x > 0 && labels_(y, x-1) == l) 747 { 748 float costColor = (diff(image1, y + dy1, x + dx1 - 1, image2, y + dy2, x + dx2) + 749 diff(image1, y + dy1, x + dx1, image2, y + dy2, x + dx2 - 1)) / 2; 750 if (costFunc_ == COLOR) 751 costV(y - roi.y, x - roi.x) = costColor; 752 else if (costFunc_ == COLOR_GRAD) 753 { 754 float costGrad = std::abs(gradx1_(y + dy1, x + dx1)) + std::abs(gradx1_(y + dy1, x + dx1 - 1)) + 755 std::abs(gradx2_(y + dy2, x + dx2)) + std::abs(gradx2_(y + dy2, x + dx2 - 1)) + 1.f; 756 costV(y - roi.y, x - roi.x) = costColor / costGrad; 757 } 758 } 759 else 760 costV(y - roi.y, x - roi.x) = badRegionCost; 761 } 762 } 763 764 costH.create(roi.height+1, roi.width); 765 766 for (int y = roi.y; y < roi.br().y+1; ++y) 767 { 768 for (int x = roi.x; x < roi.br().x; ++x) 769 { 770 if (labels_(y, x) == l && y > 0 && labels_(y-1, x) == l) 771 { 772 float costColor = (diff(image1, y + dy1 - 1, x + dx1, image2, y + dy2, x + dx2) + 773 diff(image1, y + dy1, x + dx1, image2, y + dy2 - 1, x + dx2)) / 2; 774 if (costFunc_ == COLOR) 775 costH(y - roi.y, x - roi.x) = costColor; 776 else if (costFunc_ == COLOR_GRAD) 777 { 778 float costGrad = std::abs(grady1_(y + dy1, x + dx1)) + std::abs(grady1_(y + dy1 - 1, x + dx1)) + 779 std::abs(grady2_(y + dy2, x + dx2)) + std::abs(grady2_(y + dy2 - 1, x + dx2)) + 1.f; 780 costH(y - roi.y, x - roi.x) = costColor / costGrad; 781 } 782 } 783 else 784 costH(y - roi.y, x - roi.x) = badRegionCost; 785 } 786 } 787} 788 789 790bool DpSeamFinder::estimateSeam( 791 const Mat &image1, const Mat &image2, Point tl1, Point tl2, int comp, 792 Point p1, Point p2, std::vector<Point> &seam, bool &isHorizontal) 793{ 794 CV_Assert(states_[comp] & INTERS); 795 796 Mat_<float> costV, costH; 797 computeCosts(image1, image2, tl1, tl2, comp, costV, costH); 798 799 Rect roi(tls_[comp], brs_[comp]); 800 Point src = p1 - roi.tl(); 801 Point dst = p2 - roi.tl(); 802 int l = comp+1; 803 804 // estimate seam direction 805 806 bool swapped = false; 807 isHorizontal = std::abs(dst.x - src.x) > std::abs(dst.y - src.y); 808 809 if (isHorizontal) 810 { 811 if (src.x > dst.x) 812 { 813 std::swap(src, dst); 814 swapped = true; 815 } 816 } 817 else if (src.y > dst.y) 818 { 819 swapped = true; 820 std::swap(src, dst); 821 } 822 823 // find optimal control 824 825 Mat_<uchar> control = Mat::zeros(roi.size(), CV_8U); 826 Mat_<uchar> reachable = Mat::zeros(roi.size(), CV_8U); 827 Mat_<float> cost = Mat::zeros(roi.size(), CV_32F); 828 829 reachable(src) = 1; 830 cost(src) = 0.f; 831 832 int nsteps; 833 std::pair<float, int> steps[3]; 834 835 if (isHorizontal) 836 { 837 for (int x = src.x+1; x <= dst.x; ++x) 838 { 839 for (int y = 0; y < roi.height; ++y) 840 { 841 // seam follows along upper side of pixels 842 843 nsteps = 0; 844 845 if (labels_(y + roi.y, x + roi.x) == l) 846 { 847 if (reachable(y, x-1)) 848 steps[nsteps++] = std::make_pair(cost(y, x-1) + costH(y, x-1), 1); 849 if (y > 0 && reachable(y-1, x-1)) 850 steps[nsteps++] = std::make_pair(cost(y-1, x-1) + costH(y-1, x-1) + costV(y-1, x), 2); 851 if (y < roi.height-1 && reachable(y+1, x-1)) 852 steps[nsteps++] = std::make_pair(cost(y+1, x-1) + costH(y+1, x-1) + costV(y, x), 3); 853 } 854 855 if (nsteps) 856 { 857 std::pair<float, int> opt = *min_element(steps, steps + nsteps); 858 cost(y, x) = opt.first; 859 control(y, x) = (uchar)opt.second; 860 reachable(y, x) = 255; 861 } 862 } 863 } 864 } 865 else 866 { 867 for (int y = src.y+1; y <= dst.y; ++y) 868 { 869 for (int x = 0; x < roi.width; ++x) 870 { 871 // seam follows along left side of pixels 872 873 nsteps = 0; 874 875 if (labels_(y + roi.y, x + roi.x) == l) 876 { 877 if (reachable(y-1, x)) 878 steps[nsteps++] = std::make_pair(cost(y-1, x) + costV(y-1, x), 1); 879 if (x > 0 && reachable(y-1, x-1)) 880 steps[nsteps++] = std::make_pair(cost(y-1, x-1) + costV(y-1, x-1) + costH(y, x-1), 2); 881 if (x < roi.width-1 && reachable(y-1, x+1)) 882 steps[nsteps++] = std::make_pair(cost(y-1, x+1) + costV(y-1, x+1) + costH(y, x), 3); 883 } 884 885 if (nsteps) 886 { 887 std::pair<float, int> opt = *min_element(steps, steps + nsteps); 888 cost(y, x) = opt.first; 889 control(y, x) = (uchar)opt.second; 890 reachable(y, x) = 255; 891 } 892 } 893 } 894 } 895 896 if (!reachable(dst)) 897 return false; 898 899 // restore seam 900 901 Point p = dst; 902 seam.clear(); 903 seam.push_back(p + roi.tl()); 904 905 if (isHorizontal) 906 { 907 for (; p.x != src.x; seam.push_back(p + roi.tl())) 908 { 909 if (control(p) == 2) p.y--; 910 else if (control(p) == 3) p.y++; 911 p.x--; 912 } 913 } 914 else 915 { 916 for (; p.y != src.y; seam.push_back(p + roi.tl())) 917 { 918 if (control(p) == 2) p.x--; 919 else if (control(p) == 3) p.x++; 920 p.y--; 921 } 922 } 923 924 if (!swapped) 925 std::reverse(seam.begin(), seam.end()); 926 927 CV_Assert(seam.front() == p1); 928 CV_Assert(seam.back() == p2); 929 return true; 930} 931 932 933void DpSeamFinder::updateLabelsUsingSeam( 934 int comp1, int comp2, const std::vector<Point> &seam, bool isHorizontalSeam) 935{ 936 Mat_<int> mask = Mat::zeros(brs_[comp1].y - tls_[comp1].y, 937 brs_[comp1].x - tls_[comp1].x, CV_32S); 938 939 for (size_t i = 0; i < contours_[comp1].size(); ++i) 940 mask(contours_[comp1][i] - tls_[comp1]) = 255; 941 942 for (size_t i = 0; i < seam.size(); ++i) 943 mask(seam[i] - tls_[comp1]) = 255; 944 945 // find connected components after seam carving 946 947 int l1 = comp1+1, l2 = comp2+1; 948 949 int ncomps = 0; 950 951 for (int y = 0; y < mask.rows; ++y) 952 for (int x = 0; x < mask.cols; ++x) 953 if (!mask(y, x) && labels_(y + tls_[comp1].y, x + tls_[comp1].x) == l1) 954 floodFill(mask, Point(x, y), ++ncomps); 955 956 for (size_t i = 0; i < contours_[comp1].size(); ++i) 957 { 958 int x = contours_[comp1][i].x - tls_[comp1].x; 959 int y = contours_[comp1][i].y - tls_[comp1].y; 960 961 bool ok = false; 962 static const int dx[] = {-1, +1, 0, 0, -1, +1, -1, +1}; 963 static const int dy[] = {0, 0, -1, +1, -1, -1, +1, +1}; 964 965 for (int j = 0; j < 8; ++j) 966 { 967 int c = x + dx[j]; 968 int r = y + dy[j]; 969 970 if (c >= 0 && c < mask.cols && r >= 0 && r < mask.rows && 971 mask(r, c) && mask(r, c) != 255) 972 { 973 ok = true; 974 mask(y, x) = mask(r, c); 975 } 976 } 977 978 if (!ok) 979 mask(y, x) = 0; 980 } 981 982 if (isHorizontalSeam) 983 { 984 for (size_t i = 0; i < seam.size(); ++i) 985 { 986 int x = seam[i].x - tls_[comp1].x; 987 int y = seam[i].y - tls_[comp1].y; 988 989 if (y < mask.rows-1 && mask(y+1, x) && mask(y+1, x) != 255) 990 mask(y, x) = mask(y+1, x); 991 else 992 mask(y, x) = 0; 993 } 994 } 995 else 996 { 997 for (size_t i = 0; i < seam.size(); ++i) 998 { 999 int x = seam[i].x - tls_[comp1].x; 1000 int y = seam[i].y - tls_[comp1].y; 1001 1002 if (x < mask.cols-1 && mask(y, x+1) && mask(y, x+1) != 255) 1003 mask(y, x) = mask(y, x+1); 1004 else 1005 mask(y, x) = 0; 1006 } 1007 } 1008 1009 // find new components connected with the second component and 1010 // with other components except the ones we are working with 1011 1012 std::map<int, int> connect2; 1013 std::map<int, int> connectOther; 1014 1015 for (int i = 1; i <= ncomps; ++i) 1016 { 1017 connect2.insert(std::make_pair(i, 0)); 1018 connectOther.insert(std::make_pair(i, 0)); 1019 } 1020 1021 for (size_t i = 0; i < contours_[comp1].size(); ++i) 1022 { 1023 int x = contours_[comp1][i].x; 1024 int y = contours_[comp1][i].y; 1025 1026 if ((x > 0 && labels_(y, x-1) == l2) || 1027 (y > 0 && labels_(y-1, x) == l2) || 1028 (x < unionSize_.width-1 && labels_(y, x+1) == l2) || 1029 (y < unionSize_.height-1 && labels_(y+1, x) == l2)) 1030 { 1031 connect2[mask(y - tls_[comp1].y, x - tls_[comp1].x)]++; 1032 } 1033 1034 if ((x > 0 && labels_(y, x-1) != l1 && labels_(y, x-1) != l2) || 1035 (y > 0 && labels_(y-1, x) != l1 && labels_(y-1, x) != l2) || 1036 (x < unionSize_.width-1 && labels_(y, x+1) != l1 && labels_(y, x+1) != l2) || 1037 (y < unionSize_.height-1 && labels_(y+1, x) != l1 && labels_(y+1, x) != l2)) 1038 { 1039 connectOther[mask(y - tls_[comp1].y, x - tls_[comp1].x)]++; 1040 } 1041 } 1042 1043 std::vector<int> isAdjComp(ncomps + 1, 0); 1044 1045 for (std::map<int, int>::iterator itr = connect2.begin(); itr != connect2.end(); ++itr) 1046 { 1047 double len = static_cast<double>(contours_[comp1].size()); 1048 isAdjComp[itr->first] = itr->second / len > 0.05 && connectOther.find(itr->first)->second / len < 0.1; 1049 } 1050 1051 // update labels 1052 1053 for (int y = 0; y < mask.rows; ++y) 1054 for (int x = 0; x < mask.cols; ++x) 1055 if (mask(y, x) && isAdjComp[mask(y, x)]) 1056 labels_(y + tls_[comp1].y, x + tls_[comp1].x) = l2; 1057} 1058 1059 1060class GraphCutSeamFinder::Impl : public PairwiseSeamFinder 1061{ 1062public: 1063 Impl(int cost_type, float terminal_cost, float bad_region_penalty) 1064 : cost_type_(cost_type), terminal_cost_(terminal_cost), bad_region_penalty_(bad_region_penalty) {} 1065 1066 ~Impl() {} 1067 1068 void find(const std::vector<UMat> &src, const std::vector<Point> &corners, std::vector<UMat> &masks); 1069 void findInPair(size_t first, size_t second, Rect roi); 1070 1071private: 1072 void setGraphWeightsColor(const Mat &img1, const Mat &img2, 1073 const Mat &mask1, const Mat &mask2, GCGraph<float> &graph); 1074 void setGraphWeightsColorGrad(const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, 1075 const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, 1076 GCGraph<float> &graph); 1077 1078 std::vector<Mat> dx_, dy_; 1079 int cost_type_; 1080 float terminal_cost_; 1081 float bad_region_penalty_; 1082}; 1083 1084 1085void GraphCutSeamFinder::Impl::find(const std::vector<UMat> &src, const std::vector<Point> &corners, 1086 std::vector<UMat> &masks) 1087{ 1088 // Compute gradients 1089 dx_.resize(src.size()); 1090 dy_.resize(src.size()); 1091 Mat dx, dy; 1092 for (size_t i = 0; i < src.size(); ++i) 1093 { 1094 CV_Assert(src[i].channels() == 3); 1095 Sobel(src[i], dx, CV_32F, 1, 0); 1096 Sobel(src[i], dy, CV_32F, 0, 1); 1097 dx_[i].create(src[i].size(), CV_32F); 1098 dy_[i].create(src[i].size(), CV_32F); 1099 for (int y = 0; y < src[i].rows; ++y) 1100 { 1101 const Point3f* dx_row = dx.ptr<Point3f>(y); 1102 const Point3f* dy_row = dy.ptr<Point3f>(y); 1103 float* dx_row_ = dx_[i].ptr<float>(y); 1104 float* dy_row_ = dy_[i].ptr<float>(y); 1105 for (int x = 0; x < src[i].cols; ++x) 1106 { 1107 dx_row_[x] = normL2(dx_row[x]); 1108 dy_row_[x] = normL2(dy_row[x]); 1109 } 1110 } 1111 } 1112 PairwiseSeamFinder::find(src, corners, masks); 1113} 1114 1115 1116void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &img2, 1117 const Mat &mask1, const Mat &mask2, GCGraph<float> &graph) 1118{ 1119 const Size img_size = img1.size(); 1120 1121 // Set terminal weights 1122 for (int y = 0; y < img_size.height; ++y) 1123 { 1124 for (int x = 0; x < img_size.width; ++x) 1125 { 1126 int v = graph.addVtx(); 1127 graph.addTermWeights(v, mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f, 1128 mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f); 1129 } 1130 } 1131 1132 // Set regular edge weights 1133 const float weight_eps = 1.f; 1134 for (int y = 0; y < img_size.height; ++y) 1135 { 1136 for (int x = 0; x < img_size.width; ++x) 1137 { 1138 int v = y * img_size.width + x; 1139 if (x < img_size.width - 1) 1140 { 1141 float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1142 normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1)) + 1143 weight_eps; 1144 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) || 1145 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1)) 1146 weight += bad_region_penalty_; 1147 graph.addEdges(v, v + 1, weight, weight); 1148 } 1149 if (y < img_size.height - 1) 1150 { 1151 float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1152 normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x)) + 1153 weight_eps; 1154 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) || 1155 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x)) 1156 weight += bad_region_penalty_; 1157 graph.addEdges(v, v + img_size.width, weight, weight); 1158 } 1159 } 1160 } 1161} 1162 1163 1164void GraphCutSeamFinder::Impl::setGraphWeightsColorGrad( 1165 const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, 1166 const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, 1167 GCGraph<float> &graph) 1168{ 1169 const Size img_size = img1.size(); 1170 1171 // Set terminal weights 1172 for (int y = 0; y < img_size.height; ++y) 1173 { 1174 for (int x = 0; x < img_size.width; ++x) 1175 { 1176 int v = graph.addVtx(); 1177 graph.addTermWeights(v, mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f, 1178 mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f); 1179 } 1180 } 1181 1182 // Set regular edge weights 1183 const float weight_eps = 1.f; 1184 for (int y = 0; y < img_size.height; ++y) 1185 { 1186 for (int x = 0; x < img_size.width; ++x) 1187 { 1188 int v = y * img_size.width + x; 1189 if (x < img_size.width - 1) 1190 { 1191 float grad = dx1.at<float>(y, x) + dx1.at<float>(y, x + 1) + 1192 dx2.at<float>(y, x) + dx2.at<float>(y, x + 1) + weight_eps; 1193 float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1194 normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1))) / grad + 1195 weight_eps; 1196 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) || 1197 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1)) 1198 weight += bad_region_penalty_; 1199 graph.addEdges(v, v + 1, weight, weight); 1200 } 1201 if (y < img_size.height - 1) 1202 { 1203 float grad = dy1.at<float>(y, x) + dy1.at<float>(y + 1, x) + 1204 dy2.at<float>(y, x) + dy2.at<float>(y + 1, x) + weight_eps; 1205 float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1206 normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x))) / grad + 1207 weight_eps; 1208 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) || 1209 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x)) 1210 weight += bad_region_penalty_; 1211 graph.addEdges(v, v + img_size.width, weight, weight); 1212 } 1213 } 1214 } 1215} 1216 1217 1218void GraphCutSeamFinder::Impl::findInPair(size_t first, size_t second, Rect roi) 1219{ 1220 Mat img1 = images_[first].getMat(ACCESS_READ), img2 = images_[second].getMat(ACCESS_READ); 1221 Mat dx1 = dx_[first], dx2 = dx_[second]; 1222 Mat dy1 = dy_[first], dy2 = dy_[second]; 1223 Mat mask1 = masks_[first].getMat(ACCESS_RW), mask2 = masks_[second].getMat(ACCESS_RW); 1224 Point tl1 = corners_[first], tl2 = corners_[second]; 1225 1226 const int gap = 10; 1227 Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); 1228 Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); 1229 Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); 1230 Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); 1231 Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1232 Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1233 Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1234 Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1235 1236 // Cut subimages and submasks with some gap 1237 for (int y = -gap; y < roi.height + gap; ++y) 1238 { 1239 for (int x = -gap; x < roi.width + gap; ++x) 1240 { 1241 int y1 = roi.y - tl1.y + y; 1242 int x1 = roi.x - tl1.x + x; 1243 if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols) 1244 { 1245 subimg1.at<Point3f>(y + gap, x + gap) = img1.at<Point3f>(y1, x1); 1246 submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1); 1247 subdx1.at<float>(y + gap, x + gap) = dx1.at<float>(y1, x1); 1248 subdy1.at<float>(y + gap, x + gap) = dy1.at<float>(y1, x1); 1249 } 1250 else 1251 { 1252 subimg1.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0); 1253 submask1.at<uchar>(y + gap, x + gap) = 0; 1254 subdx1.at<float>(y + gap, x + gap) = 0.f; 1255 subdy1.at<float>(y + gap, x + gap) = 0.f; 1256 } 1257 1258 int y2 = roi.y - tl2.y + y; 1259 int x2 = roi.x - tl2.x + x; 1260 if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols) 1261 { 1262 subimg2.at<Point3f>(y + gap, x + gap) = img2.at<Point3f>(y2, x2); 1263 submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2); 1264 subdx2.at<float>(y + gap, x + gap) = dx2.at<float>(y2, x2); 1265 subdy2.at<float>(y + gap, x + gap) = dy2.at<float>(y2, x2); 1266 } 1267 else 1268 { 1269 subimg2.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0); 1270 submask2.at<uchar>(y + gap, x + gap) = 0; 1271 subdx2.at<float>(y + gap, x + gap) = 0.f; 1272 subdy2.at<float>(y + gap, x + gap) = 0.f; 1273 } 1274 } 1275 } 1276 1277 const int vertex_count = (roi.height + 2 * gap) * (roi.width + 2 * gap); 1278 const int edge_count = (roi.height - 1 + 2 * gap) * (roi.width + 2 * gap) + 1279 (roi.width - 1 + 2 * gap) * (roi.height + 2 * gap); 1280 GCGraph<float> graph(vertex_count, edge_count); 1281 1282 switch (cost_type_) 1283 { 1284 case GraphCutSeamFinder::COST_COLOR: 1285 setGraphWeightsColor(subimg1, subimg2, submask1, submask2, graph); 1286 break; 1287 case GraphCutSeamFinder::COST_COLOR_GRAD: 1288 setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2, 1289 submask1, submask2, graph); 1290 break; 1291 default: 1292 CV_Error(Error::StsBadArg, "unsupported pixel similarity measure"); 1293 } 1294 1295 graph.maxFlow(); 1296 1297 for (int y = 0; y < roi.height; ++y) 1298 { 1299 for (int x = 0; x < roi.width; ++x) 1300 { 1301 if (graph.inSourceSegment((y + gap) * (roi.width + 2 * gap) + x + gap)) 1302 { 1303 if (mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x)) 1304 mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; 1305 } 1306 else 1307 { 1308 if (mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x)) 1309 mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; 1310 } 1311 } 1312 } 1313} 1314 1315 1316GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float bad_region_penalty) 1317 : impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {} 1318 1319GraphCutSeamFinder::~GraphCutSeamFinder() {} 1320 1321 1322void GraphCutSeamFinder::find(const std::vector<UMat> &src, const std::vector<Point> &corners, 1323 std::vector<UMat> &masks) 1324{ 1325 impl_->find(src, corners, masks); 1326} 1327 1328 1329#ifdef HAVE_OPENCV_CUDALEGACY 1330void GraphCutSeamFinderGpu::find(const std::vector<UMat> &src, const std::vector<Point> &corners, 1331 std::vector<UMat> &masks) 1332{ 1333 // Compute gradients 1334 dx_.resize(src.size()); 1335 dy_.resize(src.size()); 1336 Mat dx, dy; 1337 for (size_t i = 0; i < src.size(); ++i) 1338 { 1339 CV_Assert(src[i].channels() == 3); 1340 Sobel(src[i], dx, CV_32F, 1, 0); 1341 Sobel(src[i], dy, CV_32F, 0, 1); 1342 dx_[i].create(src[i].size(), CV_32F); 1343 dy_[i].create(src[i].size(), CV_32F); 1344 for (int y = 0; y < src[i].rows; ++y) 1345 { 1346 const Point3f* dx_row = dx.ptr<Point3f>(y); 1347 const Point3f* dy_row = dy.ptr<Point3f>(y); 1348 float* dx_row_ = dx_[i].ptr<float>(y); 1349 float* dy_row_ = dy_[i].ptr<float>(y); 1350 for (int x = 0; x < src[i].cols; ++x) 1351 { 1352 dx_row_[x] = normL2(dx_row[x]); 1353 dy_row_[x] = normL2(dy_row[x]); 1354 } 1355 } 1356 } 1357 PairwiseSeamFinder::find(src, corners, masks); 1358} 1359 1360 1361void GraphCutSeamFinderGpu::findInPair(size_t first, size_t second, Rect roi) 1362{ 1363 Mat img1 = images_[first].getMat(ACCESS_READ), img2 = images_[second].getMat(ACCESS_READ); 1364 Mat dx1 = dx_[first], dx2 = dx_[second]; 1365 Mat dy1 = dy_[first], dy2 = dy_[second]; 1366 Mat mask1 = masks_[first].getMat(ACCESS_READ), mask2 = masks_[second].getMat(ACCESS_READ); 1367 Point tl1 = corners_[first], tl2 = corners_[second]; 1368 1369 const int gap = 10; 1370 Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); 1371 Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); 1372 Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); 1373 Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); 1374 Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1375 Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1376 Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1377 Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); 1378 1379 // Cut subimages and submasks with some gap 1380 for (int y = -gap; y < roi.height + gap; ++y) 1381 { 1382 for (int x = -gap; x < roi.width + gap; ++x) 1383 { 1384 int y1 = roi.y - tl1.y + y; 1385 int x1 = roi.x - tl1.x + x; 1386 if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols) 1387 { 1388 subimg1.at<Point3f>(y + gap, x + gap) = img1.at<Point3f>(y1, x1); 1389 submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1); 1390 subdx1.at<float>(y + gap, x + gap) = dx1.at<float>(y1, x1); 1391 subdy1.at<float>(y + gap, x + gap) = dy1.at<float>(y1, x1); 1392 } 1393 else 1394 { 1395 subimg1.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0); 1396 submask1.at<uchar>(y + gap, x + gap) = 0; 1397 subdx1.at<float>(y + gap, x + gap) = 0.f; 1398 subdy1.at<float>(y + gap, x + gap) = 0.f; 1399 } 1400 1401 int y2 = roi.y - tl2.y + y; 1402 int x2 = roi.x - tl2.x + x; 1403 if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols) 1404 { 1405 subimg2.at<Point3f>(y + gap, x + gap) = img2.at<Point3f>(y2, x2); 1406 submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2); 1407 subdx2.at<float>(y + gap, x + gap) = dx2.at<float>(y2, x2); 1408 subdy2.at<float>(y + gap, x + gap) = dy2.at<float>(y2, x2); 1409 } 1410 else 1411 { 1412 subimg2.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0); 1413 submask2.at<uchar>(y + gap, x + gap) = 0; 1414 subdx2.at<float>(y + gap, x + gap) = 0.f; 1415 subdy2.at<float>(y + gap, x + gap) = 0.f; 1416 } 1417 } 1418 } 1419 1420 Mat terminals, leftT, rightT, top, bottom; 1421 1422 switch (cost_type_) 1423 { 1424 case GraphCutSeamFinder::COST_COLOR: 1425 setGraphWeightsColor(subimg1, subimg2, submask1, submask2, 1426 terminals, leftT, rightT, top, bottom); 1427 break; 1428 case GraphCutSeamFinder::COST_COLOR_GRAD: 1429 setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2, 1430 submask1, submask2, terminals, leftT, rightT, top, bottom); 1431 break; 1432 default: 1433 CV_Error(Error::StsBadArg, "unsupported pixel similarity measure"); 1434 } 1435 1436 cuda::GpuMat terminals_d(terminals); 1437 cuda::GpuMat leftT_d(leftT); 1438 cuda::GpuMat rightT_d(rightT); 1439 cuda::GpuMat top_d(top); 1440 cuda::GpuMat bottom_d(bottom); 1441 cuda::GpuMat labels_d, buf_d; 1442 1443 cuda::graphcut(terminals_d, leftT_d, rightT_d, top_d, bottom_d, labels_d, buf_d); 1444 1445 Mat_<uchar> labels = (Mat)labels_d; 1446 for (int y = 0; y < roi.height; ++y) 1447 { 1448 for (int x = 0; x < roi.width; ++x) 1449 { 1450 if (labels(y + gap, x + gap)) 1451 { 1452 if (mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x)) 1453 mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; 1454 } 1455 else 1456 { 1457 if (mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x)) 1458 mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; 1459 } 1460 } 1461 } 1462} 1463 1464 1465void GraphCutSeamFinderGpu::setGraphWeightsColor(const Mat &img1, const Mat &img2, const Mat &mask1, const Mat &mask2, 1466 Mat &terminals, Mat &leftT, Mat &rightT, Mat &top, Mat &bottom) 1467{ 1468 const Size img_size = img1.size(); 1469 1470 terminals.create(img_size, CV_32S); 1471 leftT.create(Size(img_size.height, img_size.width), CV_32S); 1472 rightT.create(Size(img_size.height, img_size.width), CV_32S); 1473 top.create(img_size, CV_32S); 1474 bottom.create(img_size, CV_32S); 1475 1476 Mat_<int> terminals_(terminals); 1477 Mat_<int> leftT_(leftT); 1478 Mat_<int> rightT_(rightT); 1479 Mat_<int> top_(top); 1480 Mat_<int> bottom_(bottom); 1481 1482 // Set terminal weights 1483 for (int y = 0; y < img_size.height; ++y) 1484 { 1485 for (int x = 0; x < img_size.width; ++x) 1486 { 1487 float source = mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f; 1488 float sink = mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f; 1489 terminals_(y, x) = saturate_cast<int>((source - sink) * 255.f); 1490 } 1491 } 1492 1493 // Set regular edge weights 1494 const float weight_eps = 1.f; 1495 for (int y = 0; y < img_size.height; ++y) 1496 { 1497 for (int x = 0; x < img_size.width; ++x) 1498 { 1499 if (x > 0) 1500 { 1501 float weight = normL2(img1.at<Point3f>(y, x - 1), img2.at<Point3f>(y, x - 1)) + 1502 normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1503 weight_eps; 1504 if (!mask1.at<uchar>(y, x - 1) || !mask1.at<uchar>(y, x) || 1505 !mask2.at<uchar>(y, x - 1) || !mask2.at<uchar>(y, x)) 1506 weight += bad_region_penalty_; 1507 leftT_(x, y) = saturate_cast<int>(weight * 255.f); 1508 } 1509 else 1510 leftT_(x, y) = 0; 1511 1512 if (x < img_size.width - 1) 1513 { 1514 float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1515 normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1)) + 1516 weight_eps; 1517 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) || 1518 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1)) 1519 weight += bad_region_penalty_; 1520 rightT_(x, y) = saturate_cast<int>(weight * 255.f); 1521 } 1522 else 1523 rightT_(x, y) = 0; 1524 1525 if (y > 0) 1526 { 1527 float weight = normL2(img1.at<Point3f>(y - 1, x), img2.at<Point3f>(y - 1, x)) + 1528 normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1529 weight_eps; 1530 if (!mask1.at<uchar>(y - 1, x) || !mask1.at<uchar>(y, x) || 1531 !mask2.at<uchar>(y - 1, x) || !mask2.at<uchar>(y, x)) 1532 weight += bad_region_penalty_; 1533 top_(y, x) = saturate_cast<int>(weight * 255.f); 1534 } 1535 else 1536 top_(y, x) = 0; 1537 1538 if (y < img_size.height - 1) 1539 { 1540 float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1541 normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x)) + 1542 weight_eps; 1543 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) || 1544 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x)) 1545 weight += bad_region_penalty_; 1546 bottom_(y, x) = saturate_cast<int>(weight * 255.f); 1547 } 1548 else 1549 bottom_(y, x) = 0; 1550 } 1551 } 1552} 1553 1554 1555void GraphCutSeamFinderGpu::setGraphWeightsColorGrad( 1556 const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, 1557 const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, 1558 Mat &terminals, Mat &leftT, Mat &rightT, Mat &top, Mat &bottom) 1559{ 1560 const Size img_size = img1.size(); 1561 1562 terminals.create(img_size, CV_32S); 1563 leftT.create(Size(img_size.height, img_size.width), CV_32S); 1564 rightT.create(Size(img_size.height, img_size.width), CV_32S); 1565 top.create(img_size, CV_32S); 1566 bottom.create(img_size, CV_32S); 1567 1568 Mat_<int> terminals_(terminals); 1569 Mat_<int> leftT_(leftT); 1570 Mat_<int> rightT_(rightT); 1571 Mat_<int> top_(top); 1572 Mat_<int> bottom_(bottom); 1573 1574 // Set terminal weights 1575 for (int y = 0; y < img_size.height; ++y) 1576 { 1577 for (int x = 0; x < img_size.width; ++x) 1578 { 1579 float source = mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f; 1580 float sink = mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f; 1581 terminals_(y, x) = saturate_cast<int>((source - sink) * 255.f); 1582 } 1583 } 1584 1585 // Set regular edge weights 1586 const float weight_eps = 1.f; 1587 for (int y = 0; y < img_size.height; ++y) 1588 { 1589 for (int x = 0; x < img_size.width; ++x) 1590 { 1591 if (x > 0) 1592 { 1593 float grad = dx1.at<float>(y, x - 1) + dx1.at<float>(y, x) + 1594 dx2.at<float>(y, x - 1) + dx2.at<float>(y, x) + weight_eps; 1595 float weight = (normL2(img1.at<Point3f>(y, x - 1), img2.at<Point3f>(y, x - 1)) + 1596 normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x))) / grad + 1597 weight_eps; 1598 if (!mask1.at<uchar>(y, x - 1) || !mask1.at<uchar>(y, x) || 1599 !mask2.at<uchar>(y, x - 1) || !mask2.at<uchar>(y, x)) 1600 weight += bad_region_penalty_; 1601 leftT_(x, y) = saturate_cast<int>(weight * 255.f); 1602 } 1603 else 1604 leftT_(x, y) = 0; 1605 1606 if (x < img_size.width - 1) 1607 { 1608 float grad = dx1.at<float>(y, x) + dx1.at<float>(y, x + 1) + 1609 dx2.at<float>(y, x) + dx2.at<float>(y, x + 1) + weight_eps; 1610 float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1611 normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1))) / grad + 1612 weight_eps; 1613 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) || 1614 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1)) 1615 weight += bad_region_penalty_; 1616 rightT_(x, y) = saturate_cast<int>(weight * 255.f); 1617 } 1618 else 1619 rightT_(x, y) = 0; 1620 1621 if (y > 0) 1622 { 1623 float grad = dy1.at<float>(y - 1, x) + dy1.at<float>(y, x) + 1624 dy2.at<float>(y - 1, x) + dy2.at<float>(y, x) + weight_eps; 1625 float weight = (normL2(img1.at<Point3f>(y - 1, x), img2.at<Point3f>(y - 1, x)) + 1626 normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x))) / grad + 1627 weight_eps; 1628 if (!mask1.at<uchar>(y - 1, x) || !mask1.at<uchar>(y, x) || 1629 !mask2.at<uchar>(y - 1, x) || !mask2.at<uchar>(y, x)) 1630 weight += bad_region_penalty_; 1631 top_(y, x) = saturate_cast<int>(weight * 255.f); 1632 } 1633 else 1634 top_(y, x) = 0; 1635 1636 if (y < img_size.height - 1) 1637 { 1638 float grad = dy1.at<float>(y, x) + dy1.at<float>(y + 1, x) + 1639 dy2.at<float>(y, x) + dy2.at<float>(y + 1, x) + weight_eps; 1640 float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) + 1641 normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x))) / grad + 1642 weight_eps; 1643 if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) || 1644 !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x)) 1645 weight += bad_region_penalty_; 1646 bottom_(y, x) = saturate_cast<int>(weight * 255.f); 1647 } 1648 else 1649 bottom_(y, x) = 0; 1650 } 1651 } 1652} 1653#endif 1654 1655} // namespace detail 1656} // namespace cv 1657