6、尋找接縫線

6.1 原理

拼接影象的另一個重要的步驟是找到影象重疊部分內的一條接縫線，該接縫是重疊部分最相似的畫素的連線。當確定了接縫線後，在重疊部分，線的一側只選擇該側的影象部分，線的另一側只選擇這一側的影象部分，而不是把重疊部分的兩幅影象簡單融合起來。這麼做的目的可以避免影象的模糊及偽像。

目前，常用的尋找接縫線的方法有三種：逐點法、動態規劃法和圖割法。

逐點法比較簡單，它的原理就是重疊部分內的畫素距離哪個影象更近，就用哪個影象中的相應畫素值。這種方法並不是基於畫素值，只是基於距離準則。所以這種方法並不能真正得到最佳的接縫線，尋找的結果比較粗糙原始，僅僅是能夠分割重疊部分而已。

動態規劃法可以快速的找到最佳的接縫線，並且所需的記憶體較少，因此該方法很適合在移動終端裝置上進行高清晰度的全景拼接。

設影象1和影象2之間有重疊部分，我們需要得到它們之間的最佳接縫，首先定義重疊部分的誤差表面函式：

（88）

式中，I₁，I₂表示兩幅影象各自的重疊部分，││,||表示範數。

一般來說，接縫線有三個限制條件：一是如果重疊區域是寬大於高，則接縫是橫向走向的，而如果重疊區域是寬小於高，則接縫是縱向走向的，也就是要保證接縫線要有一定的長度；二是如果是橫向的接縫，則不允許有絕對垂直的接縫線，而如果是縱向的接縫，則不允許有絕對水平的接縫線；三是重疊區域是矩形，接縫線從矩形的一邊出發，必須到達與該邊平行的另一邊結束。

假設重疊區域的寬小於高，則接縫線是垂直的，我們需要橫向遍歷e值，並且計算所有可能到達當前畫素(i

, j)的路徑的累積最小誤差E：

（89）

式中，min中的三項內容分別表示當前畫素與其左上側、上側和右上側的E。在E中最後一行的最小值表明已經到達了最小垂直路徑的盡頭，那麼我們就可以追溯，從而得到最佳路徑，即接縫線。

同理，重疊區域的寬大於高也能夠做出類似的計算過程。

我們可以看出，動態規劃演算法非常適用於對式89的求解，這也是該尋找最佳接縫線方法名稱的由來。

對於彩色影象，計算式88可以有兩種方法：直接法和梯度法。設p和q是兩個相鄰畫素，直接法的計算公式為：

（90）

梯度法的計算公式為：

（91）

式中，G表示灰度影象的p到q方向上的梯度。

最後一種方法是圖割法。它的主要思想是將影象看成是一個有向圖G

={V, E}，V表示節點，E表示邊。節點分為兩種：普通節點和終端節點。影象的畫素可以作為普通節點，影象的聚類結果可以看成是終端節點，即某一聚類中的普通節點會匯聚於一個終端節點，一般圖只有兩個終端節點：s和t，也就是說一般圖割法只能夠把影象分割成兩類S和T。邊指的是節點之間的連線，對於普通節點，只有相鄰節點可以連線，該邊被稱為n-links，而所有普通節點都可與屬於它們的終端節點之間有連線，該邊被稱為t-links。用最大流圖割法可以得到哪些普通節點屬於終端節點s，那麼剩下的普通節點就會屬於t。

Opencv中所實現的最大流圖割演算法是基於Boykov和Kolmogorov所提出的演算法。該演算法是以s和t為根，以得到的兩棵無不相交的搜尋樹S和T為目標，它迭代的重複以下三個階段：

生長階段：搜尋樹S和T的生長，直到兩棵樹生長到了一起，形成了從s到t的通路（即為流）；

增廣階段：該通路被增廣，使相連的兩棵樹又重新分開；

收養階段：搜尋樹S和T收養孤立的節點（稱為孤兒）。

在生長階段，如果找不到任何一條從s到t的通路，則迭代結束。在增廣階段，開始孤兒集合是空集，而該階段結束後，可能會出現一些孤兒。在該階段，我們需要得到通路上最小的容量值（也稱為邊的權值），然後所有的容量值都減去該值（即通過流量更新容量值），則該通路上至少會有一個權值為0的邊（稱為飽和邊），假設飽和邊上的兩個節點為p和q，並且在正向通路中是由p流向q，如果p和q都屬於樹S，則q為孤兒，如果p和q都屬於樹T，則p為孤兒。在收養階段，檢查孤兒集合中的所有孤兒的鄰域畫素，連線那些它們的邊的容量值不為0的畫素，如果最終可以從該孤兒連通到s或t，則該孤兒就屬於樹S或樹T，否則由該孤兒組成的樹就屬於自由節點集，在下次迭代時的生長階段，樹S和樹T就通過連線自由節點集中的節點來使兩棵樹生長髮育的。

如何用最大流圖割法尋找接縫線呢？接縫線把重疊區域分割成了兩個區域，那麼我們通過最大流圖割法得到樹S和樹T，即得到了被分割的兩個區域，則區域之間的接縫線自然就可以得到。

普通節點之間的初始容量值（即邊的權值）也可以通過直接法和梯度法計算得到，它們的計算公式分別與式90和式91相同。終端節點與普通節點之間的初始容量值可以定義為一個很大的值，當重疊區域中的節點只屬於影象1或影象2時，節點與相應的終端節點的邊賦予該值，其他情況都置為0，表示它們與終端節點無直接連線。

6.2 原始碼

SeamFinder類是尋找縫的基類，該類內有一個最重要的虛擬函式是find，它的作用是尋找縫：

virtual void find(const std::vector<Mat> &src, const std::vector<Point> &corners,
                      std::vector<Mat> &masks) = 0;

它的三個輸入引數分別為：src表示待拼接的影象，corners表示待拼接影象在最終的全景影象座標系內的左上角的座標位置，masks輸入時表示上一步得到的對映變換掩碼，輸出時表示接縫掩碼，也就是該引數是要更新輸出的，輸出時如果畫素值更新為0了，則表示該畫素是在接縫的另一側，拼接時是不需要該畫素的。

基於SeamFinder類有四個子類，分別對應於四種不同的尋找縫的演算法——NoSeamFinder（無需尋找縫）、PairwiseSeamFinder（逐點尋找縫演算法）、DpSeamFinder（動態規劃尋找縫演算法）和GraphCutSeamFinder（圖割尋找縫演算法）。

我們先來介紹PairwiseSeamFinder類。逐點尋找縫演算法又可分為不同的方法，Opencv只實現了Voronoi圖方法，實現該方法的類為VoronoiSeamFinder類，它是PairwiseSeamFinder類的子類。VoronoiSeamFinder類的find函式為：

void VoronoiSeamFinder::find(const vector<Size> &sizes, const vector<Point> &corners,
                             vector<Mat> &masks)
{
    LOGLN("Finding seams...");
    if (sizes.size() == 0)    //沒有影象，則返回
        return;

#if ENABLE_LOG
    int64 t = getTickCount();    //計時
#endif
    //引數賦值
    sizes_ = sizes;
    corners_ = corners;
    masks_ = masks;
    run();    //呼叫run函式
    //計時輸出
    LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
}

VoronoiSeamFinder類沒有實現run函式，因此這裡的find函式呼叫的是它的父類PairwiseSeamFinder類的run函式：

void PairwiseSeamFinder::run()
{
    //遍歷影象對
    for (size_t i = 0; i < sizes_.size() - 1; ++i)
    {
        for (size_t j = i + 1; j < sizes_.size(); ++j)
        {
            Rect roi;    //表示影象間的重疊部分
            //如果影象i和圖j之間有重疊的區域，則執行findInPair函式
            if (overlapRoi(corners_[i], corners_[j], sizes_[i], sizes_[j], roi))
                findInPair(i, j, roi);
        }
    }
}

findInPair函式表示尋找兩幅影象之間的接縫：

void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi)
//first和second表示兩幅影象的尺寸大小，roi表示影象間重疊部分
{
    const int gap = 10;
    //分別表示影象first和影象second的重疊部分的掩碼，它們要比roi四周多gap個距離
    Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
    Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);

    Size img1 = sizes_[first], img2 = sizes_[second];    //得到兩幅影象的尺寸
    Mat mask1 = masks_[first], mask2 = masks_[second];    //得到兩幅影象的掩碼
    Point tl1 = corners_[first], tl2 = corners_[second];    //得到兩幅影象的左上角座標

    // Cut submasks with some gap
    //遍歷重疊部分，得到submask1和submask2
    for (int y = -gap; y < roi.height + gap; ++y)
    {
        for (int x = -gap; x < roi.width + gap; ++x)
        {
            //y1和x1表示重疊部分相對於影象first的座標
            int y1 = roi.y - tl1.y + y;
            int x1 = roi.x - tl1.x + x;
            //把影象first的掩碼賦值給submask1
            //表示當前座標在影象first內
            if (y1 >= 0 && x1 >= 0 && y1 < img1.height && x1 < img1.width) 
                submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1);
            else    //表示當前座標不在影象first內
                submask1.at<uchar>(y + gap, x + gap) = 0;
            //y2和x2表示重疊部分相對於影象second的座標
            int y2 = roi.y - tl2.y + y;
            int x2 = roi.x - tl2.x + x;
            //把影象second的掩碼賦值給submask2
            //表示當前座標在影象second內
            if (y2 >= 0 && x2 >= 0 && y2 < img2.height && x2 < img2.width)
                submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2);
            else    //表示當前座標不在影象second內
                submask2.at<uchar>(y + gap, x + gap) = 0;
        }
    }
    //collision表示submask1和submask2交集，roi區域可能有些區域是無效的部分（由於對映的原因），而collision（畫素值為1的部分）則是影象first和影象second的真正有意義部分的重疊區域，我們是需要在collision區域內尋找接縫的
    Mat collision = (submask1 != 0) & (submask2 != 0);
    //複製submask1和submask2為unique1和unique2，並且把collision部分的畫素清零
    Mat unique1 = submask1.clone(); unique1.setTo(0, collision);
    Mat unique2 = submask2.clone(); unique2.setTo(0, collision);

    Mat dist1, dist2;    //表示距離矩陣
    //分別得到unique1和unique2內為0的畫素值與最近的非0值之間的L1距離，顯然為0的區域肯定要比collision大
    distanceTransform(unique1 == 0, dist1, CV_DIST_L1, 3);
    distanceTransform(unique2 == 0, dist2, CV_DIST_L1, 3);
    //得到縫矩陣，dist1 < dist2成立，說明離first影象近，則seam內相應的畫素為1，否則為0
    Mat seam = dist1 < dist2; 
    //遍歷roi內的所有畫素
    for (int y = 0; y < roi.height; ++y)
    {
        for (int x = 0; x < roi.width; ++x)
        {
            //如果seam內的畫素值為1，說明collision區域內的相應畫素離影象first近，所以該點的值選擇影象first的畫素值，因此需把mask2清零，反之把mask1清零。也就是在重疊部分，哪些區域應該屬於影象first，哪些區域屬於影象second。
            if (seam.at<uchar>(y + gap, x + gap))
                mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0;
            else
                mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0;
        }
    }
}

DpSeamFinder類利用動態規劃演算法尋找縫，它的建構函式為：

DpSeamFinder::DpSeamFinder(CostFunction costFunc) : costFunc_(costFunc) {}

其中CostFunction為：

enum CostFunction { COLOR, COLOR_GRAD };

它表示動態規劃尋找接縫的兩種不同方法，直接法和梯度法，預設是直接法。

DpSeamFinder類的find函式為：

void DpSeamFinder::find(const vector<Mat> &src, const vector<Point> &corners, vector<Mat> &masks)
{
    LOGLN("Finding seams...");
#if ENABLE_LOG
    int64 t = getTickCount();    //用於計時
#endif

    if (src.size() == 0)    //沒有影象，則返回
        return;

    vector<pair<size_t, size_t> > pairs;    //表示影象對

    for (size_t i = 0; i+1 < src.size(); ++i)
        for (size_t j = i+1; j < src.size(); ++j)
            pairs.push_back(make_pair(i, j));    //每兩幅影象組成一個影象對
    //根據影象位置（基於左上角座標），對影象對進行排序
    sort(pairs.begin(), pairs.end(), ImagePairLess(src, corners));
    reverse(pairs.begin(), pairs.end());    //逆序

    for (size_t i = 0; i < pairs.size(); ++i)    //遍歷影象對，
    {
        size_t i0 = pairs[i].first, i1 = pairs[i].second;
        //尋找i0和i1這兩幅影象的接縫，該函式在後面有介紹
        process(src[i0], src[i1], corners[i0], corners[i1], masks[i0], masks[i1]);
    }

    LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
}

尋找兩幅影象間的接縫：

void DpSeamFinder::process(
        const Mat &image1, const Mat &image2, Point tl1, Point tl2,
        Mat &mask1, Mat &mask2)
//image1和image2分別表示待處理的兩幅影象
//tl1和tl2分別表示這兩幅影象的左上角座標
//mask1和mask2分別表示這兩幅影象的掩碼
{
    CV_Assert(image1.size() == mask1.size());    //確保影象與掩碼尺寸相同
    CV_Assert(image2.size() == mask2.size());
    //intersectTl和intersectBr分別表示image1和image2的交集區域（即重疊部分）的左上角和右下角座標
    Point intersectTl(std::max(tl1.x, tl2.x), std::max(tl1.y, tl2.y));

    Point intersectBr(std::min(tl1.x + image1.cols, tl2.x + image2.cols),
                      std::min(tl1.y + image1.rows, tl2.y + image2.rows));
    //如果if條件成立，則說明image1和image2沒有重疊區域，因此退出該函式
    if (intersectTl.x >= intersectBr.x || intersectTl.y >= intersectBr.y)
        return; // there are no conflicts
    //unionTl_和unionBr_分別表示image1和image2的並集區域（即外圍最大部分的矩形）的左上角和右下角座標
    unionTl_ = Point(std::min(tl1.x, tl2.x), std::min(tl1.y, tl2.y));

    unionBr_ = Point(std::max(tl1.x + image1.cols, tl2.x + image2.cols),
                     std::max(tl1.y + image1.rows, tl2.y + image2.rows));
    //得到image1和image2的並集區域的尺寸大小
    unionSize_ = Size(unionBr_.x - unionTl_.x, unionBr_.y - unionTl_.y);
    //定義兩幅影象的掩碼大小，先清零
    mask1_ = Mat::zeros(unionSize_, CV_8U);
    mask2_ = Mat::zeros(unionSize_, CV_8U);
    //tmp為mask1_中image1的區域
    Mat tmp = mask1_(Rect(tl1.x - unionTl_.x, tl1.y - unionTl_.y, mask1.cols, mask1.rows));
    mask1.copyTo(tmp);    //把image1的掩碼mask1賦值給mask1_
    //tmp為mask2_中image2的區域
    tmp = mask2_(Rect(tl2.x - unionTl_.x, tl2.y - unionTl_.y, mask2.cols, mask2.rows));
    mask2.copyTo(tmp);    //把image2的掩碼mask2賦值給mask2_

    // find both images contour masks
    //在並集中分別得到邊界掩碼contour1mask_和contour2mask_
    contour1mask_ = Mat::zeros(unionSize_, CV_8U);    //先清零
    contour2mask_ = Mat::zeros(unionSize_, CV_8U);
    //遍歷並集區域
    for (int y = 0; y < unionSize_.height; ++y)
    {
        for (int x = 0; x < unionSize_.width; ++x)
        {
            if (mask1_(y, x) &&
                ((x == 0 || !mask1_(y, x-1)) || (x == unionSize_.width-1 || !mask1_(y, x+1)) ||
                 (y == 0 || !mask1_(y-1, x)) || (y == unionSize_.height-1 || !mask1_(y+1, x))))
            {
                contour1mask_(y, x) = 255;    //得到image1的邊界掩碼
            }

            if (mask2_(y, x) &&
                ((x == 0 || !mask2_(y, x-1)) || (x == unionSize_.width-1 || !mask2_(y, x+1)) ||
                 (y == 0 || !mask2_(y-1, x)) || (y == unionSize_.height-1 || !mask2_(y+1, x))))
            {
                contour2mask_(y, x) = 255;    //得到image2的邊界掩碼
            }
        }
    }
    //以下三個函式在後面都有介紹
    findComponents();

    findEdges();

    resolveConflicts(image1, image2, tl1, tl2, mask1, mask2);
}

利用漫水填充法分割影象尋找元件，不同的元件對應於影象的不同區域：

void DpSeamFinder::findComponents()
{
    // label all connected components and get information about them

    ncomps_ = 0;    //表示元件的數量
    labels_.create(unionSize_);    //定義標籤的大小
    states_.clear();    //狀態清空
    tls_.clear();    //tls表示各種元件的左上角座標
    brs_.clear();    //brs表示各種元件的右下角座標
    contours_.clear();    //元件區域的邊界
    //遍歷兩幅影象的並集區域，為labels_賦值
    for (int y = 0; y < unionSize_.height; ++y)
    {
        for (int x = 0; x < unionSize_.width; ++x)
        {
            if (mask1_(y, x) && mask2_(y, x))    //表示兩幅影象的公共交集的部分
                labels_(y, x) = numeric_limits<int>::max();
            else if (mask1_(y, x))    //表示只有image1的部分
                labels_(y, x) = numeric_limits<int>::max()-1;
            else if (mask2_(y, x))    //表示只有image2的部分
                labels_(y, x) = numeric_limits<int>::max()-2;
            else    //表示沒有任何影象資訊
                labels_(y, x) = 0;
        }
    }
    //再次遍歷兩幅影象的並集區域，
    for (int y = 0; y < unionSize_.height; ++y)
    {
        for (int x = 0; x < unionSize_.width; ++x)
        {
            //表示當前畫素是並集區域內有影象資訊的
            if (labels_(y, x) >= numeric_limits<int>::max()-2)
            {
                if (labels_(y, x) == numeric_limits<int>::max())    //交集部分
                    states_.push_back(INTERS);    //狀態賦值
                else if (labels_(y, x) == numeric_limits<int>::max()-1)    //只有image1部分
                    states_.push_back(FIRST);    //狀態賦值
                else if (labels_(y, x) == numeric_limits<int>::max()-2)    //只有image2部分
                    states_.push_back(SECOND);    //狀態賦值
                //利用漫水填充演算法，從當前畫素點出發，填充++ncomps_值
                floodFill(labels_, Point(x, y), ++ncomps_);
                //存入當前畫素點，表示當前元件的左上角座標
                tls_.push_back(Point(x, y)); 
                //表示當前元件的右下角座標，+1的含義是用於區分左上角座標
                brs_.push_back(Point(x+1, y+1)); 
                //因為已經得到了一個元件，所以這時需要在contours_內預留一段空間，用於儲存當前元件的邊界畫素
                contours_.push_back(vector<Point>()); 
            }

            if (labels_(y, x))    //表示當前畫素屬於某幅影象
            {
                int l = labels_(y, x);    //當前標籤
                int ci = l-1;    //標籤值減1，表示當前元件

                //確定當前元件的左上角和右下角座標，在雙重for迴圈體內，不斷進入該if語句，所以實現了當前元件的左上角和右下角座標的不斷擴充套件，最終得到準確的值
                tls_[ci].x = std::min(tls_[ci].x, x);
                tls_[ci].y = std::min(tls_[ci].y, y);
                brs_[ci].x = std::max(brs_[ci].x, x+1);
                brs_[ci].y = std::max(brs_[ci].y, y+1);
                //if條件成立，說明當前畫素為當前元件的邊界
                if ((x == 0 || labels_(y, x-1) != l) || (x == unionSize_.width-1 || labels_(y, x+1) != l) ||
                    (y == 0 || labels_(y-1, x) != l) || (y == unionSize_.height-1 || labels_(y+1, x) != l))
                {
                    contours_[ci].push_back(Point(x, y));    //存入當前元件的邊界畫素
                }
            }
        }
    }
}

得到元件間的邊緣：

void DpSeamFinder::findEdges()
{
    // find edges between components
    //wedges表示邊緣權值，即邊緣畫素的數量，pair中的兩個引數表示兩個元件的索引
    map<pair<int, int>, int> wedges; // weighted edges
    //遍歷元件，初始化wedges為0
    for (int ci = 0; ci < ncomps_-1; ++ci)
    {
        for (int cj = ci+1; cj < ncomps_; ++cj)
        {
            wedges[make_pair(ci, cj)] = 0;
            wedges[make_pair(cj, ci)] = 0;
        }
    }

    for (int ci = 0; ci < ncomps_; ++ci)    //遍歷元件
    {
        for (size_t i = 0; i < contours_[ci].size(); ++i)    //遍歷當前元件的邊界
        {
            int x = contours_[ci][i].x;    //邊緣畫素座標
            int y = contours_[ci][i].y;
            int l = ci + 1;    //當前元件的標籤
            //if條件成立，表示當前畫素為左側邊緣，此時當前元件與其左側元件間的wedges權值加1
            if (x > 0 && labels_(y, x-1) && labels_(y, x-1) != l)
            {
                wedges[make_pair(ci, labels_(y, x-1)-1)]++; 
                wedges[make_pair(labels_(y, x-1)-1, ci)]++; 
            }
            //if條件成立，表示當前畫素為上側邊緣，此時當前元件與其上側元件間的wedges權值加1
            if (y > 0 && labels_(y-1, x) && labels_(y-1, x) != l)
            {
                wedges[make_pair(ci, labels_(y-1, x)-1)]++; 
                wedges[make_pair(labels_(y-1, x)-1, ci)]++; 
            }
            //if條件成立，表示當前畫素為右側邊緣，此時當前元件與其右側元件間的wedges權值加1
            if (x < unionSize_.width-1 && labels_(y, x+1) && labels_(y, x+1) != l)
            {
                wedges[make_pair(ci, labels_(y, x+1)-1)]++;
                wedges[make_pair(labels_(y, x+1)-1, ci)]++;
            }
            //if條件成立，表示當前畫素為下側邊緣，此時當前元件與其下側元件間的wedges權值加1
            if (y < unionSize_.height-1 && labels_(y+1, x) && labels_(y+1, x) != l)
            {
                wedges[make_pair(ci, labels_(y+1, x)-1)]++;
                wedges[make_pair(labels_(y+1, x)-1, ci)]++;
            }
        }
    }
    //edges_表示兩個元件是相鄰的，因為它們之間有邊界
    edges_.clear();    //清空
    //遍歷wedges
    for (int ci = 0; ci < ncomps_-1; ++ci)
    {
        for (int cj = ci+1; cj < ncomps_; ++cj)
        {
            //得到元件ci和cj之間的wedges
            map<pair<int, int>, int>::iterator itr = wedges.find(make_pair(ci, cj));
            if (itr != wedges.end() && itr->second > 0)    //表示ci和cj是相鄰的
                edges_.insert(itr->first);    //ci和cj存入edges中
            //得到元件cj和ci之間的wedges
            itr = wedges.find(make_pair(cj, ci));
            if (itr != wedges.end() && itr->second > 0)    //表示cj和ci是相鄰的
                edges_.insert(itr->first);    //cj和ci存入edges中
        }
    }
}

具體計算接縫的部分：

void DpSeamFinder::resolveConflicts(
        const Mat &image1, const Mat &image2, Point tl1, Point tl2, Mat &mask1, Mat &mask2)
{
    if (costFunc_ == COLOR_GRAD)    //梯度法
        computeGradients(image1, image2);    //呼叫computeGradients，計算梯度

    // resolve conflicts between components

    bool hasConflict = true;    //標識變數
    while (hasConflict)
    {
        int c1 = 0, c2 = 0;
        hasConflict = false;    //變換標識
        //遍歷edges_，得到相連元件
        for (set<pair<int, int> >::iterator itr = edges_.begin(); itr != edges_.end(); ++itr)
        {
            c1 = itr->first;    //得到連線邊緣的一個元件
            c2 = itr->second;    //得到連線邊緣的另一個元件
            //c1有交集部分，並且它不是與c2元件的交集
            if ((states_[c1] & INTERS) && (states_[c1] & (~INTERS)) != states_[c2])
            {
                hasConflict = true;    //賦值為true
                break;    //退出for迴圈
            }
        }

        if (hasConflict)
        {
            int l1 = c1+1, l2 = c2+1;    //表示c1和c2的標籤
            //表示c1只有一個鄰近元件，說明c1不是交集，不需要在該元件尋找接縫線
            if (hasOnlyOneNeighbor(c1)) 
            {
                // if the first components has only one adjacent component
                //遍歷c1區域，把標籤為l1的畫素改為l2，即合併兩個元件
                for (int y = tls_[c1].y; y < brs_[c1].y; ++y)
                    for (int x = tls_[c1].x; x < brs_[c1].x; ++x)
                        if (labels_(y, x) == l1)
                            labels_(y, x) = l2;
                //改變元件c1的狀態
                states_[c1] = states_[c2] == FIRST ? SECOND : FIRST;
            }
            //表示c1有不止一個鄰近元件，說明c1是交集部分，要在該元件尋找接縫線
            else
            {
                // if the first component has more than one adjacent component

                Point p1, p2;    //表示接縫的起點和終點
                //利用getSeamTips函式得到p1和p2，該函式在後面給出介紹
                if (getSeamTips(c1, c2, p1, p2))
                {
                    vector<Point> seam;    //表示接縫
                    bool isHorizontalSeam;    //表示接縫的基本走向
                    //利用estimateSeam函式得到接縫seam，該函式在後面給出介紹
                    if (estimateSeam(image1, image2, tl1, tl2, c1, p1, p2, seam, isHorizontalSeam))
                        //更新標籤，該函式在後面給出介紹
                        updateLabelsUsingSeam(c1, c2, seam, isHorizontalSeam);
                }
                //重新為c1的狀態賦值
                states_[c1] = states_[c2] == FIRST ? INTERS_SECOND : INTERS_FIRST;
            }

            const int c[] = {c1, c2};    //表示元件
            const int l[] = {l1, l2};    //表示標籤

            for (int i = 0; i < 2; ++i)    //遍歷c1和c2兩個元件
            {
                // update information about the (i+1)-th component
                //得到當前元件的左上角和右下角座標
                int x0 = tls_[c[i]].x, x1 = brs_[c[i]].x; 
                int y0 = tls_[c[i]].y, y1 = brs_[c[i]].y;
                //重新賦值
                tls_[c[i]] = Point(numeric_limits<int>::max(), numeric_limits<int>::max());
                brs_[c[i]] = Point(numeric_limits<int>::min(), numeric_limits<int>::min());
                contours_[c[i]].clear();    //當前元件的邊緣畫素清零

                for (int y = y0; y < y1; ++y)    //遍歷當前元件
                {
                    for (int x = x0; x < x1; ++x)
                    {
                        if (labels_(y, x) == l[i])
                        {
                            //更新當前元件的區域面積
                            tls_[c[i]].x = std::min(tls_[c[i]].x, x);
                            tls_[c[i]].y = std::min(tls_[c[i]].y, y);
                            brs_[c[i]].x = std::max(brs_[c[i]].x, x+1);
                            brs_[c[i]].y = std::max(brs_[c[i]].y, y+1);

                            if ((x == 0 || labels_(y, x-1) != l[i]) || (x == unionSize_.width-1 || labels_(y, x+1) != l[i]) ||
                                (y == 0 || labels_(y-1, x) != l[i]) || (y == unionSize_.height-1 || labels_(y+1, x) != l[i]))
                            {
                                //給出當前元件的邊緣畫素
                                contours_[c[i]].push_back(Point(x, y)); 
                            }
                        }
                    }
                }
            }

            // remove edges
            //刪除c1和c2所形成的邊緣
            edges_.erase(make_pair(c1, c2));
            edges_.erase(make_pair(c2, c1));
        }
    }

    // update masks
    //更新掩碼
    int dx1 = unionTl_.x - tl1.x, dy1 = unionTl_.y - tl1.y;    //影象1的左上角座標
    int dx2 = unionTl_.x - tl2.x, dy2 = unionTl_.y - tl2.y;    //影象2的左上角座標
    //遍歷mask2，更新掩碼2
    for (int y = 0; y < mask2.rows; ++y)
    {
        for (int x = 0; x < mask2.cols; ++x)
        {
             int l = labels_(y - dy2, x - dx2);    //當前畫素的標籤
             //if成立，說明當前畫素應該用影象1的內容，所以掩碼2清零
             if (l > 0 && (states_[l-1] & FIRST) && mask1.at<uchar>(y - dy2 + dy1, x - dx2 + dx1))
                mask2.at<uchar>(y, x) = 0;
        }
    }
    //遍歷mask1，更新掩碼1
    for (int y = 0; y < mask1.rows; ++y)
    {
        for (int x = 0; x < mask1.cols; ++x)
        {
             int l = labels_(y - dy1, x - dx1);    //當前畫素的標籤
             //if成立，說明當前畫素應該用影象2的內容，所以掩碼1清零
             if (l > 0 && (states_[l-1] & SECOND) && mask2.at<uchar>(y - dy1 + dy2, x - dx1 + dx2))
                mask1.at<uchar>(y, x) = 0;
        }
    }
}

計算彩色影象image1和image2的水平和垂直梯度：

void DpSeamFinder::computeGradients(const Mat &image1, const Mat &image2)
{
    //確保image1和image2是彩色影象
    CV_Assert(image1.channels() == 3 || image1.channels() == 4);
    CV_Assert(image2.channels() == 3 || image2.channels() == 4);
    CV_Assert(costFunction() == COLOR_GRAD);    //確保是梯度法

    Mat gray;
    //把彩色影象image1轉換為灰度影象
    if (image1.channels() == 3)
        cvtColor(image1, gray, CV_BGR2GRAY);
    else if (image1.channels() == 4)
        cvtColor(image1, gray, CV_BGRA2GRAY);

    Sobel(gray, gradx1_, CV_32F, 1, 0);    //得到image1的水平梯度gradx1_
    Sobel(gray, grady1_, CV_32F, 0, 1);    //得到image1的垂直梯度grady1_
    //把彩色影象image2轉換為灰度影象
    if (image2.channels() == 3)
        cvtColor(image2, gray, CV_BGR2GRAY);
    else if (image2.channels() == 4)
        cvtColor(image2, gray, CV_BGRA2GRAY);

    Sobel(gray, gradx2_, CV_32F, 1, 0);    //得到image2的水平梯度gradx2_
    Sobel(gray, grady2_, CV_32F, 0, 1);    //得到image2的垂直梯度grady2_
}

計算接縫的起點和終點：

bool DpSeamFinder::getSeamTips(int comp1, int comp2, Point &p1, Point &p2)
//comp1和comp2分別表示相鄰的元件
//p1和p2分別表示得到返回comp1中，聚類最遠的兩個聚類中的畫素點，即接縫的起點和終點
//返回為true，則說明得到了p1和p2，否則返回false
{
    CV_Assert(states_[comp1] & INTERS);    //確保comp1為影象間的重疊部分

    // find special points

    vector<Point> specialPoints;    //表示一些特殊的點，即靠近其他邊界的點
    int l2 = comp2+1;    //表示comp2的標籤

    for (size_t i = 0; i < contours_[comp1].size(); ++i)    //遍歷comp1的邊緣畫素
    {
        int x = contours_[comp1][i].x;    //得到當前邊緣畫素的座標
        int y = contours_[comp1][i].y;
        //if條件成立，說明當前畫素臨近並集的邊界，並且靠近comp2
        if (closeToContour(y, x, contour1mask_) &&
            closeToContour(y, x, contour2mask_) &&
            ((x > 0 && labels_(y, x-1) == l2) ||
             (y > 0 && labels_(y-1, x) == l2) ||
             (x < unionSize_.width-1 && labels_(y, x+1) == l2) ||
             (y < unionSize_.height-1 && labels_(y+1, x) == l2)))
        {
            specialPoints.push_back(Point(x, y));    //儲存這些特殊的點
        }
    }
    //特殊點數量太少，不能形成聚類，則直接返回
    if (specialPoints.size() < 2)
        return false;

    // find clusters

    vector<int> labels;
    //把specialPoints中的相距不大於10個單位的點聚成一類，並賦予相應的標籤labels
    cv::partition(specialPoints, labels, ClosePoints(10));
    //得到最大聚類的畫素數量
    int nlabels = *max_element(labels.begin(), labels.end()) + 1;
    if (nlabels < 2)    //聚類中的畫素數量太少，無法得到中心畫素，則直接返回
        return false;

    vector<Point> sum(nlabels);    //表示各個聚類中畫素點座標值之和
    vector<vector<Point> > points(nlabels);    //用於儲存各個聚類的畫素

    for (size_t i = 0; i < specialPoints.size(); ++i)    //遍歷specialPoints
    {
        sum[labels[i]] += specialPoints[i];
        points[labels[i]].push_back(specialPoints[i]);
    }

    // select two most distant clusters

    int idx[2] = {-1,-1};    //表示距離最大的兩個聚類的索引
    //表示兩個聚類間的最大距離，先初始化為最小值
    double maxDist = -numeric_limits<double>::max();
    //兩兩聚類遍歷，得到距離最大的兩個聚類
    for (int i = 0; i < nlabels-1; ++i)
    {
        for (int j = i+1; j < nlabels; ++j)
        {
            //分別表示聚類1和聚類2的特殊點的數量
            double size1 = static_cast<double>(points[i].size()), size2 = static_cast<double>(points[j].size());
            //分別表示聚類1和聚類2的橫、縱座標的均值，代表聚類的中心
            double cx1 = cvRound(sum[i].x / size1), cy1 = cvRound(sum[i].y / size1);
            double cx2 = cvRound(sum[j].x / size2), cy2 = cvRound(sum[j].y / size1);
            //聚類i和聚類j間的距離
            double dist = (cx1 - cx2) * (cx1 - cx2) + (cy1 - cy2) * (cy1 - cy2);
            if (dist > maxDist)
            {
                maxDist = dist;    //更新最大距離
                idx[0] = i;    //更新兩個聚類
                idx[1] = j;
            }
        }
    }

    // select two points closest to the clusters' centers

    Point p[2];
    //遍歷距離最大的那兩個聚類，得到聚類內離中心最近的點
    for (int i = 0; i < 2; ++i)
    {
        double size = static_cast<double>(points[idx[i]].size());    //聚類內的特殊點的數量
        double cx = cvRound(sum[idx[i]].x / size);    //聚類的橫座標的均值
        double cy = cvRound(sum[idx[i]].y / size);    //聚類的縱座標的均值

        size_t closest = points[idx[i]].size();
        double minDist = numeric_limits<double>::max();    //先設定一個最大值

        for (size_t j = 0; j < points[idx[i]].size(); ++j)    //遍歷聚類內的特殊點
        {
            //得到當前特殊點與該聚類中心的距離
            double dist = (points[idx[i]][j].x - cx) * (points[idx[i]][j].x - cx) +
                          (points[idx[i]][j].y - cy) * (points[idx[i]][j].y - cy);
            if (dist < minDist) 
            {
                minDist = dist;    //更新最小距離
                closest = j;    //更新聚類索引
            }
        }

        p[i] = points[idx[i]][closest];    //得到該點
    }

    p1 = p[0];    //得到這兩個畫素點
    p2 = p[1];
    return true;
}

計算得到接縫：

bool DpSeamFinder::estimateSeam(
        const Mat &image1, const Mat &image2, Point tl1, Point tl2, int comp,
        Point p1, Point p2, vector<Point> &seam, bool &isHorizontal)
//image1和image2分別表示重疊的兩幅影象
//tl1和tl2分別表示這兩幅影象在全景影象中的左上角的座標
//comp表示這兩幅影象的重疊部分，即交集元件
//p1和p2分別表示接縫的起點和終點
//seam表示得到的接縫
//isHorizontal表示接縫的基本走向是否為水平
//如果得到了接縫，該函式返回true，否則返回false
{
    CV_Assert(states_[comp] & INTERS);    //確保comp為交集元件

    Mat_<float> costV, costH;    //表示垂直和水平的代價
    //計算兩幅影象間重疊區域的水平和垂直的代價值，該函式在後面給出介紹
    computeCosts(image1, image2, tl1, tl2, comp, costV, costH);

    Rect roi(tls_[comp], brs_[comp]);    //定義重疊區域的矩形形式
    Point src = p1 - roi.tl();    //src表示p1點在roi的絕對位置
    Point dst = p2 - roi.tl();    //dst表示p2點在roi的絕對位置
    int l = comp+1;    //表示comp1的標籤

    // estimate seam direction

    bool swapped = false;    //定義一個標籤，表示src和dst交換
    //由src和dst構成的矩形如果是寬大於高，則isHorizontal為true，否則為false
    isHorizontal = std::abs(dst.x - src.x) > std::abs(dst.y - src.y);

    if (isHorizontal)    //寬大於高
    {
        if (src.x > dst.x)    //比較橫座標
        {
            std::swap(src, dst);    //交換，使src在dst的左邊
            swapped = true;
        }
    }
    else if (src.y > dst.y)    //高大於寬，則比較縱座標
    {
        swapped = true;
        std::swap(src, dst);    //交換，使src在dst的上邊
    }

    // find optimal control

    Mat_<uchar> control = Mat::zeros(roi.size(), CV_8U);    //表示控制量，即接縫的走向
    Mat_<uchar> reachable = Mat::zeros(roi.size(), CV_8U);    //表示已遍歷到該畫素
    Mat_<float> cost = Mat::zeros(roi.size(), CV_32F);    //表示代價值，即累積最小誤差

    reachable(src) = 1;    //從src開始遍歷
    cost(src) = 0.f;

    int nsteps;    //表示步長
    pair<float, int> steps[3];    //表示接縫的走向，並計算式89中min中的三個元素

    if (isHorizontal)    //寬大於高，接縫的走向基本上是從左向右
    {
        for (int x = src.x+1; x <= dst.x; ++x)    //遍歷重疊部分
        {
            for (int y = 0; y < roi.height; ++y)
            {
                // seam follows along upper side of pixels

                nsteps = 0;    //每次迴圈，步長都要清零

                if (labels_(y + roi.y, x + roi.x) == l)    //當前畫素的標籤為l
                {
                    if (reachable(y, x-1))
                        steps[nsteps++] = make_pair(cost(y, x-1) + costH(y, x-1), 1);
                    if (y > 0 && reachable(y-1, x-1))
                        steps[nsteps++] = make_pair(cost(y-1, x-1) + costH(y-1, x-1) + costV(y-1, x), 2);
                    if (y < roi.height-1 && reachable(y+1, x-1))
                        steps[nsteps++] = make_pair(cost(y+1, x-1) + costH(y+1, x-1) + costV(y, x), 3);
                }

                if (nsteps)    //表示該次迴圈，已經得到了steps變數
                {
                    //在steps陣列中，得到最小值，即式89
                    pair<float, int> opt = *min_element(steps, steps + nsteps);
                    cost(y, x) = opt.first;    //當前畫素的代價值
                    control(y, x) = (uchar)opt.second;    //當前畫素的控制量
                    reachable(y, x) = 255;    //表示已遍歷到該點
                }
            }
        }
    }
    else    //高大於寬，接縫的走向基本上是從上向下
    {
        for (int y = src.y+1; y <= dst.y; ++y)    //遍歷重疊部分
        {
            for (int x = 0; x < roi.width; ++x)
            {
                // seam follows along left side of pixels

                nsteps = 0;    //每次迴圈，步長都要清零

                if (labels_(y + roi.y, x + roi.x) == l)    //當前畫素的標籤為l
                {
                    if (reachable(y-1, x))
                        steps[nsteps++] = make_pair(cost(y-1, x) + costV(y-1, x), 1);
                    if (x > 0 && reachable(y-1, x-1))
                        steps[nsteps++] = make_pair(cost(y-1, x-1) + costV(y-1, x-1) + costH(y, x-1), 2);
                    if (x < roi.width-1 && reachable(y-1, x+1))
                        steps[nsteps++] = make_pair(cost(y-1, x+1) + costV(y-1, x+1) + costH(y, x), 3);
                }

                if (nsteps)    //表示該次迴圈，已經得到了steps變數
                {
                    //在steps陣列中，得到最小值，即式89
                    pair<float, int> opt = *min_element(steps, steps + nsteps);
                    cost(y, x) = opt.first;    //當前畫素的代價值
                    control(y, x) = (uchar)opt.second;    //當前畫素的控制量
                    reachable(y, x) = 255;    //表示已遍歷到該點
                }
            }
        }
    }

    if (!reachable(dst))    //沒有遍歷到dst點（終點），則直接退出
        return false;

    // restore seam

    Point p = dst;    //接縫的終點
    seam.clear();    //接縫清零
    seam.push_back(p + roi.tl());    //終點放入seam中

    if (isHorizontal)    //寬大於高
    {
        for (; p.x != src.x; seam.push_back(p + roi.tl()))
        {
            //根據控制量，得到接縫的座標，並把它存放入seam內
            if (control(p) == 2) p.y--;
            else if (control(p) == 3) p.y++;
            p.x--;
        }
    }
    else    //高大於寬
    {
        for (; p.y != src.y; seam.push_back(p + roi.tl()))
        {
            //根據控制量，得到接縫的座標，並把它存放入seam內
            if (control(p) == 2) p.x--;
            else if (control(p) == 3) p.x++;
            p.y--;
        }
    }

    if (!swapped)    //確保起點為p1，終點為p2
        reverse(seam.begin(), seam.end());

    CV_Assert(seam.front() == p1);    //確保起點為p1，終點為p2
    CV_Assert(seam.back() == p2);
    return true;
}

計算代價的函式，程式並不是直接計算重疊區域相對畫素之間的差值平方和，而是根據接縫線的走向，計算兩者之間水平方向或垂直方向上相差一個畫素之間的差值平方和，即式90，因此代價值分為水平值和垂直值兩類的：

void DpSeamFinder::computeCosts(
        const Mat &image1, const Mat &image2, Point tl1, Point tl2,
        int comp, Mat_<float> &costV, Mat_<float> &costH)
{
    CV_Assert(states_[comp] & INTERS);    //確保comp為交集部分

    // compute costs

    float (*diff)(const Mat&, int, int, const Mat&, int, int) = 0;    //定義一個函式指標
    //根據影象型別的不同，diff函式指標指向不同的函式，函式的本質都是比較兩幅影象中相對應元素之間的差值，即紅、綠、藍三通道強度差值平方之和
    if (image1.type() == CV_32FC3 && image2.type() == CV_32FC3)
        diff = diffL2Square3<float>;    //三通道
    else if (image1.type() == CV_8UC3 && image2.type() == CV_8UC3)
        diff = diffL2Square3<uchar>;    //三通道
    else if (image1.type() == CV_32FC4 && image2.type() == CV_32FC4)
        diff = diffL2Square4<float>;    //四通道
    else if (image1.type() == CV_8UC4 && image2.type() == CV_8UC4)
        diff = diffL2Square4<uchar>;    //四通道
    else
        CV_Error(CV_StsBadArg, "both images must have CV_32FC3(4) or CV_8UC3(4) type");

    int l = comp+1;    //表示標籤
    Rect roi(tls_[comp], brs_[comp]);    //得到comp的矩形形式
    //得到image1和image2左上角的相對座標
    int dx1 = unionTl_.x - tl1.x, dy1 = unionTl_.y - tl1.y;
    int dx2 = unionTl_.x - tl2.x, dy2 = unionTl_.y - tl2.y;
    //定義一個常數，用來表示無效的值
    const float badRegionCost = normL2(Point3f(255.f, 255.f, 255.f),
                                       Point3f(0.f, 0.f, 0.f));

    costV.create(roi.height, roi.width+1);    //建立costV的大小

    for (int y = roi.y; y < roi.br().y; ++y)    //遍歷roi區域
    {
        for (int x = roi.x; x < roi.br().x+1; ++x)
        {
            if (labels_(y, x) == l && x > 0 && labels_(y, x-1) == l)    //當前畫素的標籤為l
            {
                //在水平方向，計算式90
                float costColor = (diff(image1, y + dy1, x + dx1 - 1, image2, y + dy2, x + dx2) +
                                   diff(image1, y + dy1, x + dx1, image2, y + dy2, x + dx2 - 1)) / 2;
                if (costFunc_ == COLOR)    //直接法
                    costV(y - roi.y, x - roi.x) = costColor;
                else if (costFunc_ == COLOR_GRAD)    //梯度法
                {
                    //得到式91的分母部分
                    float costGrad = std::abs(gradx1_(y + dy1, x + dx1)) + std::abs(gradx1_(y + dy1, x + dx1 - 1)) +
                                     std::abs(gradx2_(y + dy2, x + dx2)) + std::abs(gradx2_(y + dy2, x + dx2 - 1)) + 1.f;
                    costV(y - roi.y, x - roi.x) = costColor / costGrad;    //式91
                }
            }
            else
                costV(y - roi.y, x - roi.x) = badRegionCost;    //無效的值
        }
    }

    costH.create(roi.height+1, roi.width);    //建立costV的大小

    for (int y = roi.y; y < roi.br().y+1; ++y)    //遍歷roi區域
    {
        for (int x = roi.x; x < roi.br().x; ++x)
        {
            if (labels_(y, x) == l && y > 0 && labels_(y-1, x) == l)    //當前畫素的標籤為l
            {
                //在垂直方向，計算式90
                float costColor = (diff(image1, y + dy1 - 1, x + dx1, image2, y + dy2, x + dx2) +
                                   diff(image1, y + dy1, x + dx1, image2, y + dy2 - 1, x + dx2)) / 2;
                if (costFunc_ == COLOR)    //直接法
                    costH(y - roi.y, x - roi.x) = costColor;
                else if (costFunc_ == COLOR_GRAD)        //梯度法
                {
                    //得到式91的分母部分
                    float costGrad = std::abs(grady1_(y + dy1, x + dx1)) + std::abs(grady1_(y + dy1 - 1, x + dx1)) +
                                     std::abs(grady2_(y + dy2, x + dx2)) + std::abs(grady2_(y + dy2 - 1, x + dx2)) + 1.f;
                    costH(y - roi.y, x - roi.x) = costColor / costGrad;    //式91
                }
            }
            else
                costH(y - roi.y, x - roi.x) = badRegionCost;    //無效的值
        }
    }
}

利用接縫來進行更新標籤：

void DpSeamFinder::updateLabelsUsingSeam(
        int comp1, int comp2, const vector<Point> &seam, bool isHorizontalSeam)
{
    //表示comp1的掩碼，先清零
    Mat_<int> mask = Mat::zeros(brs_[comp1].y - tls_[comp1].y,
                                brs_[comp1].x - tls_[comp1].x, CV_32S);
    //把comp1中的邊緣的掩碼賦值為255
    for (size_t i = 0; i < contours_[comp1].size(); ++i)
        mask(contours_[comp1][i] - tls_[comp1]) = 255;
    //把comp1中的接縫的掩碼賦值為255
    for (size_t i = 0; i < seam.size(); ++i)
        mask(seam[i] - tls_[comp1]) = 255;

    // find connected components after seam carving

    int l1 = comp1+1, l2 = comp2+1;    //得到兩個元件的標籤

    int ncomps = 0;    //表示元件數量
    //在mask1中，對除了邊緣和接縫以為的區域利用漫水填充演算法進行分割
    for (int y = 0; y < mask.rows; ++y)
        for (int x = 0; x < mask.cols; ++x)
            if (!mask(y, x) && labels_(y + tls_[comp1].y, x + tls_[comp1].x) == l1)
                floodFill(mask, Point(x, y), ++ncomps);

    for (size_t i = 0; i < contours_[comp1].size(); ++i)    //遍歷comp1元件邊緣畫素
    {
        int x = contours_[comp1][i].x - tls_[comp1].x;    //相對座標
        int y = contours_[comp1][i].y - tls_[comp1].y;

        bool ok = false;    //標識變數
        static const int dx[] = {-1, +1, 0, 0, -1, +1, -1, +1};    //表示8鄰域
        static const int dy[] = {0, 0, -1, +1, -1, -1, +1, +1};

        for (int j = 0; j < 8; ++j)    //遍歷當前畫素的8鄰域
        {
            int c = x + dx[j];
            int r = y + dy[j];
            //表示當前畫素的8鄰域屬於由上一步漫水填充演算法得到的分割區域
            if (c >= 0 && c < mask.cols && r >= 0 && r < mask.rows &&
                mask(r, c) && mask(r, c) != 255)
            {
                ok = true;    //賦值
                mask(y, x) = mask(r, c);    //屬於同一個分割區域
            }
        }

        if (!ok)
            mask(y, x) = 0;    //表示當前畫素為掩碼畫素
    }

    if (isHorizontalSeam)    //寬大於高
    {
        for (size_t i = 0; i < seam.size(); ++i)    //遍歷接縫
        {
            int x = seam[i].x - tls_[comp1].x;    //相對座標
            int y = seam[i].y - tls_[comp1].y;
            //當前畫素的下一行畫素的掩碼不為0或255
            if (y < mask.rows-1 && mask(y+1, x) && mask(y+1, x) != 255)
                mask(y, x) = mask(y+1, x);    //賦值
            else
                mask(y, x) = 0;
        }
    }
    else    //高大於寬
    {
        for (size_t i = 0; i < seam.size(); ++i)    //遍歷接縫
        {
            int x = seam[i].x - tls_[comp1].x;    //相對座標
            int y = seam[i].y - tls_[comp1].y;
            //當前畫素的右邊畫素的掩碼不為0或255
            if (x < mask.cols-1 && mask(y, x+1) && mask(y, x+1) != 255)
                mask(y, x) = mask(y, x+1);    //賦值
            else
                mask(y, x) = 0;
        }
    }

    // find new components connected with the second component and
    // with other components except the ones we are working with

    map<int, int> connect2;    //表示與第2個元件相鄰
    map<int, int> connectOther;    //表示與其他元件相鄰

    for (int i = 1; i <= ncomps; ++i)    //遍歷分割區域
    {
        connect2.insert(make_pair(i, 0));    //初始化
        connectOther.insert(make_pair(i, 0));
    }

    for (size_t i = 0; i < contours_[comp1].size(); ++i)    //遍歷comp1元件的邊緣畫素
    {
        int x = contours_[comp1][i].x;    //當前畫素的相對座標
        int y = contours_[comp1][i].y;
        //當前畫素的4鄰域範圍內只要有一個畫素的標籤為l2
        if ((x > 0 && labels_(y, x-1) == l2) ||
            (y > 0 && labels_(y-1, x) == l2) ||
            (x < unionSize_.width-1 && labels_(y, x+1) == l2) ||
            (y < unionSize_.height-1 && labels_(y+1, x) == l2))
        {
            //記錄該畫素，並且數量累計
            connect2[mask(y - tls_[comp1].y, x - tls_[comp1].x)]++;
        }
        //當前畫素的4鄰域範圍內只要有一個畫素的標籤為其他值
        if ((x > 0 && labels_(y, x-1) != l1 && labels_(y, x-1) != l2) ||
            (y > 0 && labels_(y-1, x) != l1 && labels_(y-1, x) != l2) ||
            (x < unionSize_.width-1 && labels_(y, x+1) != l1 && labels_(y, x+1) != l2) ||
            (y < unionSize_.height-1 && labels_(y+1, x) != l1 && labels_(y+1, x) != l2))
        {
            //記錄該畫素，並且數量累計
            connectOther[mask(y - tls_[comp1].y, x - tls_[comp1].x)]++;
        }
    }

    vector<int> isAdjComp(ncomps + 1, 0);    //表示鄰近元件
    //遍歷connect2內的畫素
    for (map<int, int>::iterator itr = connect2.begin(); itr != connect2.end(); ++itr)
    {
        //comp1元件內邊緣畫素的數量
        double len = static_cast<double>(contours_[comp1].size());
        //通過比較connect2和connectOther的數量來為isAdjComp賦值
        isAdjComp[itr->first] = itr->second / len > 0.05 && connectOther.find(itr->first)->second / len < 0.1;
    }

    // update labels
    //更新標籤
    for (int y = 0; y < mask.rows; ++y)
        for (int x = 0; x < mask.cols; ++x)
            if (mask(y, x) && isAdjComp[mask(y, x)])
                labels_(y + tls_[comp1].y, x + tls_[comp1].x) = l2;
}

GraphCutSeamFinder類為圖割法實現類，它的建構函式為：

GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float bad_region_penalty)
    : impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {}
//cost_type表示圖割法的型別，直接法還是梯度法，預設為梯度法
//terminal_cost表示終端節點與普通節點之間的初始邊權值，預設值為10000
//bad_region_penalty表示圖中無效邊的權值，預設值為1000

它主要例項化了impl_類為Impl，而GraphCutSeamFinder類中的find函式為：

void GraphCutSeamFinder::find(const vector<Mat> &src, const vector<Point> &corners,
                              vector<Mat> &masks)
{
    impl_->find(src, corners, masks);
}

所以對於GraphCutSeamFinder類，重點應該瞭解它內部的Impl類。

Impl類的find函式：

void GraphCutSeamFinder::Impl::find(const vector<Mat> &src, const vector<Point> &corners,
                                    vector<Mat> &masks)
{
    // Compute gradients
    //計算梯度
    dx_.resize(src.size());    //表示水平梯度
    dy_.resize(src.size());    //表示垂直梯度
    Mat dx, dy;    //表示當前影象的水平梯度和垂直梯度
    for (size_t i = 0; i < src.size(); ++i)    //遍歷所有待拼接的影象
    {
        CV_Assert(src[i].channels() == 3);    //確保影象是三通道的彩色影象
        Sobel(src[i], dx, CV_32F, 1, 0);    //得到當前影象的紅綠藍各個通道的水平梯度dx
        Sobel(src[i], dy, CV_32F, 0, 1);    //得到當前影象的紅綠藍各個通道的垂直梯度dy
        dx_[i].create(src[i].size(), CV_32F);    //定義大小
        dy_[i].create(src[i].size(), CV_32F);    //定義大小
        for (int y = 0; y < src[i].rows; ++y)    //遍歷當前影象的行
        {
            const Point3f* dx_row = dx.ptr<Point3f>(y);    //當前影象水平梯度行指標
            const Point3f* dy_row = dy.ptr<Point3f>(y);    //當前影象垂直梯度行指標
            float* dx_row_ = dx_[i].ptr<float>(y);    //所有影象的水平梯度行指標
            float* dy_row_ = dy_[i].ptr<float>(y);    //所有影象的垂直梯度行指標
            for (int x = 0; x < src[i].cols; ++x)    //遍歷當前行的每個元素
            {
                dx_row_[x] = normL2(dx_row[x]);    //三通道水平梯度平方和
                dy_row_[x] = normL2(dy_row[x]);    //三通道垂直梯度平方和
            }
        }
    }
    PairwiseSeamFinder::find(src, corners, masks);    //呼叫find函式
}

PairwiseSeamFinder類中的find函式主要是呼叫了run函式，而run函式在前面已給出介紹，它實際又呼叫了findInPair函式：

void GraphCutSeamFinder::Impl::findInPair(size_t first, size_t second, Rect roi)
{
    Mat img1 = images_[first], img2 = images_[second];    //表示相交的兩幅影象
    Mat dx1 = dx_[first], dx2 = dx_[second];    //表示這兩幅影象的水平梯度
    Mat dy1 = dy_[first], dy2 = dy_[second];    //表示這兩幅影象的垂直梯度
    Mat mask1 = masks_[first], mask2 = masks_[second];    //表示這兩幅影象的掩碼
    Point tl1 = corners_[first], tl2 = corners_[second];    //表示這兩幅影象的左上角座標

    const int gap = 10;
    //分別定義這兩幅影象相應的子影象，它們的面積要大一些
    Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3);
    Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3);
    Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
    Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
    Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
    Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
    Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
    Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);

    // Cut subimages and submasks with some gap
    for (int y = -gap; y < roi.height + gap; ++y)    //遍歷重疊區域
    {
        for (int x = -gap; x < roi.width + gap; ++x)
        {
            int y1 = roi.y - tl1.y + y;    //當前畫素在影象1中的座標
            int x1 = roi.x - tl1.x + x;
            if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols)
            //屬於影象1，則賦值相應的值
            {
                subimg1.at<Point3f>(y + gap, x + gap) = img1.at<Point3f>(y1, x1);
                submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1);
                subdx1.at<float>(y + gap, x + gap) = dx1.at<float>(y1, x1);
                subdy1.at<float>(y + gap, x + gap) = dy1.at<float>(y1, x1);
            }
            else    //不屬於影象1，則賦值為0
            {
                subimg1.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0);
                submask1.at<uchar>(y + gap, x + gap) = 0;
                subdx1.at<float>(y + gap, x + gap) = 0.f;
                subdy1.at<float>(y + gap, x + gap) = 0.f;
            }

            int y2 = roi.y - tl2.y + y;    //當前畫素在影象2中的座標
            int x2 = roi.x - tl2.x + x;
            if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols)
            //屬於影象2，則賦值相應的值
            {
                subimg2.at<Point3f>(y + gap, x + gap) = img2.at<Point3f>(y2, x2);
                submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2);
                subdx2.at<float>(y + gap, x + gap) = dx2.at<float>(y2, x2);
                subdy2.at<float>(y + gap, x + gap) = dy2.at<float>(y2, x2);
            }
            else    //不屬於影象1，則賦值為0
            {
                subimg2.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0);
                submask2.at<uchar>(y + gap, x + gap) = 0;
                subdx2.at<float>(y + gap, x + gap) = 0.f;
                subdy2.at<float>(y + gap, x + gap) = 0.f;
            }
        }
    }
    //表示普通節點數量
    const int vertex_count = (roi.height + 2 * gap) * (roi.width + 2 * gap);
    //表示邊數量
    const int edge_count = (roi.height - 1 + 2 * gap) * (roi.width + 2 * gap) +
                           (roi.width - 1 + 2 * gap) * (roi.height + 2 * gap);
    GCGraph<float> graph(vertex_count, edge_count);    //定義圖割變數

    switch (cost_type_)    //根據不同的方法，呼叫不同的函式
    {
    case GraphCutSeamFinder::COST_COLOR:    //直接法
        setGraphWeightsColor(subimg1, subimg2, submask1, submask2, graph);
        break;
    case GraphCutSeamFinder::COST_COLOR_GRAD:    //梯度法
        setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2,
                                 submask1, submask2, graph);
        break;
    default:    //目前只實現了上述那兩種演算法
        CV_Error(CV_StsBadArg, "unsupported pixel similarity measure");
    }

    graph.maxFlow();    //得到最大流

    for (int y = 0; y < roi.height; ++y)    //遍歷重疊區域，更新掩碼
    {
        for (int x = 0; x < roi.width; ++x)
        {
            //當前畫素屬於終端節點s
            if (graph.inSourceSegment((y + gap) * (roi.width + 2 * gap) + x + gap))
            {
                if (mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x))
                    mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0;    //掩碼影象2
            }
            else    //當前畫素屬於終端節點t
            {
                if (mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x))
                    mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0;    //掩碼影象1
            }
        }
    }
}

圖割中的直接法：

void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &img2,
                                                    const Mat &mask1, const Mat &mask2, GCGraph<float> &graph)
{
    const Size img_size = img1.size();    //得到重疊區域的尺寸

    // Set terminal weights
    for (int y = 0; y < img_size.height; ++y)    //遍歷重疊區域
    {
        for (int x = 0; x < img_size.width; ++x)
        {
            int v = graph.addVtx();    //為圖新增普通節點
            //為普通節點與兩個終端節點新增邊，即定義其初始權值
            graph.addTermWeights(v, mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f,
                                    mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f);
        }
    }

    // Set regular edge weights
    const float weight_eps = 1.f;    //設定一個規則邊權值
    for (int y = 0; y < img_size.height; ++y)    //遍歷重疊區域
    {
        for (int x = 0; x < img_size.width; ++x)
        {
            int v = y * img_size.width + x;    //表示當前普通節點，即當前畫素
            if (x < img_size.width - 1)
            {
                //計算當前畫素與其右側畫素的邊的權值，式90的分子部分
                float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
                               normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1)) +
                               weight_eps;
                //如果當前畫素和其右側畫素是被mask1和mask2掩碼掉的，則該邊無效
                if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) ||
                    !mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1))
                    weight += bad_region_penalty_;    //賦於一個較大的值
                //為當前畫素與其右側畫素之間的邊賦予初始權值，權值是雙向的
                graph.addEdges(v, v + 1, weight, weight);
            }
            if (y < img_size.height