1. 概述
2. 實作
- 2.1. 準備
- 2.2. 核心
  - 2.2.1. 均值坐標（Mean-Value Coordinates）
  - 2.2.2. ROI邊界柵格化
  - 2.2.3. 核心實作
  - 2.2.4. 實作中的問題
3. 效果
- 3.1. 使用程序
- 3.2. 效率
4. 參考

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），對于如下多邊形內部的點：

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

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

2.2.2. ROI邊界柵格化

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

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

2.2.3. 核心實作

論文給出的演算法偽代碼如下： MCV融合

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

假設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<<"毫秒";