1. 概述

泊松融合是圖像融合處理效果最好的算法，其來自於2004年Siggraph的經典paper：《Poisson Image Editing》。以這篇文章為發端，很多大神提出了一系列的優化算法。2009年, Zeev Farbman 在的SIGGRAPH上面提出的基於Mean-Value Coordinates方法的泊松融合加速算法《Coordinates for Instant Image Cloning》（文獻二）。在這篇文章中，泊松方程被轉換成拉普拉斯方程，並且提出了用均值坐標Mean-Value Coordinates來近似求解這個方程，從而達到實時運算的效果。

初步了解了一下原生的泊松融合算法和均值坐標融合算法，其原理包含的內涵十分豐富，包含一些諸如列散度、拉普拉斯算子、梯度場、泊松方程等等數學知識，要完全弄明白確實需要一定的基礎。這里就重點關注一下根據《Coordinates for Instant Image Cloning》（文獻二）實現圖像融合的過程，有機會的話再詳細推導一下其原理。

2. 實現

2.1. 准備

在OpenCV中，已經收錄了泊松融合算法，也就是函數seamlessClone()：

這個算法要求輸入一個源圖像，一個目標圖像，源圖像希望融合到目標圖像的位置，以及一個mask圖像。這個mask圖像也就是一張二值化圖像，用來標識圖像的ROI(region of interest感興趣區域)。均值坐標融合算法的輸入參數也是一樣的，不過mask圖像很難以處理，OpenCV自帶的GUI難以滿足需求。所以我這里通過QT來做GUI，通過OpenCV將圖像顯示到QT窗體上，然后再QT窗體的圖像區域內繪制多邊形，多邊形內部即為ROI。可以參考我的這兩篇文章：
《使用QT顯示OpenCV讀取的圖片》
《使用QT繪制一個多邊形》

2.2. 核心

2.2.1. 均值坐標（Mean-Value Coordinates）

在論文中提出了一個很重要的概念也就是均值坐標（Mean-Value Coordinates）。對於如下多邊形內部的點：

都有一系列與多邊形邊界相關的坐標值：

也就是說，只要確定了ROI，也就確定了ROI區域內每個點的均值坐標（Mean-Value Coordinates），每個點會有m個值（m為ROI邊界多邊形的頂點）。

2.2.2. ROI邊界柵格化

論文中是以ROI邊界多邊形為例的，實際用到圖像處理中是不會只用幾個多邊形的節點來計算的，而應該是ROI邊界上連續的點。實際上不用想也知道，圖像融合最關鍵的部分就是ROI邊界部分的像素值。必須要用到ROI邊界上所有的像素值來計算。

也就是說這里還需要一個工作，就是將ROI邊界多邊形柵格化，取得其上連續的像素位置，得到准確的柵格化多邊形邊界。這里可以參看我的這篇文章《矢量線的一種柵格化算法》。按照順序逐條將多邊形的邊柵格化，即可以得到ROI的柵格化多邊形邊界。

2.2.3. 核心實現

論文給出的算法偽代碼如下：

這段算法描述並不復雜，轉換成自然語言如下：

假設ROI區域內有n個點，其邊界由m個點組成。
那么可以求每個點的MVC(均值坐標)，每個點有m個坐標值，一共有n個點，MVC就是就是一個n*m的矩陣。
求ROI區域邊界的像素差diff，顯然其是一個m*1的矩陣。
那么新圖像ROI區域的插值為：r = MVC * diff，矩陣乘法后r為n*1矩陣。
將插值r與原圖像g矩陣相加：f = g + r，替換目標圖像相應位置的值。

核心部分具體的實現代碼如下：

    QTime startTime = QTime::currentTime();

    //Step1:找到邊界上所有的像素點
    vector<Vector2d> ROIBoundPointList;
    CalBoundPoint(ROIBoundPointList);

    //Step2:計算范圍內每個點的 mean-value coordinates
    size_t srcImgBufNum = static_cast<size_t>(srcImg.cols) * static_cast<size_t>(srcImg.rows);
    vector<vector<double>> MVC(srcImgBufNum);
    for(size_t i = 0; i < srcImgBufNum; i++)
    {
        MVC[i].resize(ROIBoundPointList.size()-1, 0);
    }
    vector<bool> clipMap(srcImgBufNum, true);           //標識范圍內的點

    cout<<"開始計算 mean-value coordinates..." << endl;
    #pragma omp parallel for        //開啟OpenMP並行加速
    for (int ri = 0; ri < srcImg.rows; ++ri)
    {
        for (int ci = 0; ci < srcImg.cols; ++ci)
        {
            //點是否在多邊形內
            size_t m = static_cast<size_t>(srcImg.cols) * ri + ci;
            if(!Point_In_Polygon_2D(ci, ri, ROIBoundPointList))
            {
                clipMap[m] = false;
                continue;
            }

            //逐點計算MVC
            Vector2d P(ci, ri);
            vector<double> alphaAngle(ROIBoundPointList.size());
            for(size_t pi = 1; pi < ROIBoundPointList.size(); pi++)
            {
                alphaAngle[pi] = threePointCalAngle(ROIBoundPointList[pi-1], P, ROIBoundPointList[pi]);
            }
            alphaAngle[0] = alphaAngle[ROIBoundPointList.size()-1];


            for(size_t pi = 1; pi < ROIBoundPointList.size(); pi++)
            {
                double w_a = tan(alphaAngle[pi-1]/2) + tan(alphaAngle[pi]/2);
                double w_b = (ROIBoundPointList[pi-1] - P).Mod();
                MVC[m][pi-1] = w_a / w_b;
                if(_isnan(MVC[m][pi-1])==1)
                {
                    MVC[m][pi-1] = 0;
                }
            }

            double sum = 0;
            for(size_t pi = 0; pi < MVC[m].size(); pi++)
            {
                sum = sum + MVC[m][pi];
            }

            for(size_t pi = 0; pi < MVC[m].size(); pi++)
            {
                MVC[m][pi] = MVC[m][pi] / sum;
            }
        }
    }
    cout<<"計算完成！" << endl;

    //Step3:計算邊界的像素插值
    vector<int> diff;
    for(size_t i = 0; i < ROIBoundPointList.size()-1; i++)
    {
        size_t l = (size_t) srcImg.cols * ROIBoundPointList[i].y + ROIBoundPointList[i].x;
        for(int bi = 0; bi < winBandNum; bi++)
        {
            size_t m = (size_t) dstImg.cols * winBandNum * (ROIBoundPointList[i].y + posY)+ winBandNum * (ROIBoundPointList[i].x + posX) + bi;
            size_t n = (size_t) srcImg.cols * winBandNum * ROIBoundPointList[i].y + winBandNum * ROIBoundPointList[i].x + bi;
            int d = (int)(dstImg.data[m]) - (int)(srcImg.data[n]);
            diff.push_back(d);
        }
        clipMap[l] = false;         //在多邊形邊上的點沒法計算MVC
    }

    //Step4:插值計算
    cout<<"開始插值計算..." << endl;
    //Mat rMat(srcImg.rows, srcImg.cols, CV_64FC3);
    #pragma omp parallel for
    for (int ri = 0; ri < srcImg.rows; ++ri)
    {
        for (int ci = 0; ci < srcImg.cols; ++ci)
        {
            size_t l = (size_t) srcImg.cols * ri + ci;
            if(!clipMap[l])
            {
                continue;
            }

            vector<double> r(winBandNum, 0);

            for(size_t pi = 0; pi < MVC[l].size(); pi++)
            {
                for(int bi = 0; bi < winBandNum; bi++)
                {
                    r[bi] = r[bi] + MVC[l][pi] * diff[pi * winBandNum + bi];
                }
            }

            for(int bi = 0; bi < winBandNum; bi++)
            {
                size_t n = (size_t) srcImg.cols * winBandNum * ri + winBandNum * ci + bi;
                size_t m = (size_t) dstImg.cols * winBandNum * (ri + posY)+ winBandNum * (ci + posX) + bi;

                dstImg.data[m] = min(max(srcImg.data[n] + r[bi], 0.0), 255.0);
            }
        }
    }
    cout<<"插值完成！" << endl;

    QTime stopTime = QTime::currentTime();
    int elapsed = startTime.msecsTo(stopTime);
    cout<<"總結完成用時"<<elapsed<<"毫秒";