Home | History | Annotate | Download | only in src
      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 
     46 namespace cv {
     47 namespace detail {
     48 
     49 void 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 
     72 void 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 
     85 void 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 
     91 void 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 
    111 void 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 
    165 DpSeamFinder::DpSeamFinder(CostFunction costFunc) : costFunc_(costFunc) {}
    166 
    167 
    168 void 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 
    202 void 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 
    266 void 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 
    332 void 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 
    399 void 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 
    526 void 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 
    552 bool 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 
    561 bool 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 
    584 bool 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 
    686 namespace
    687 {
    688 
    689 template <typename T>
    690 float 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 
    699 template <typename T>
    700 float 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 
    711 void 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 
    790 bool 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 
    933 void 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 
   1060 class GraphCutSeamFinder::Impl : public PairwiseSeamFinder
   1061 {
   1062 public:
   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 
   1071 private:
   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 
   1085 void 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 
   1116 void 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 
   1164 void 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 
   1218 void 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 
   1316 GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float bad_region_penalty)
   1317     : impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {}
   1318 
   1319 GraphCutSeamFinder::~GraphCutSeamFinder() {}
   1320 
   1321 
   1322 void 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
   1330 void 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 
   1361 void 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 
   1465 void 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 
   1555 void 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