基於均值坐標(Mean-Value Coordinates)的影像融合演算法的具體實現

目錄

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()：

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

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

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融合

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

  1. 假設ROI區域內有n個點，其邊界由m個點組成。
  2. 那麼可以求每個點的MVC(均值坐標)，每個點有m個坐標值，一共有n個點，MVC就是就是一個n*m的矩陣。
  3. 求ROI區域邊界的像素差diff，顯然其是一個m*1的矩陣。
  4. 那麼新影像ROI區域的插值為：r = MVC * diff，矩陣乘法後r為n*1矩陣。
  5. 將插值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<<"毫秒";

2.2.4. 實現中的問題

  1. ROI邊界上的點無法計算MVC值，需要予以剔除，否則ROI邊界上會出現一圈白色的點。
  2. 用到了OpenMP加速，可以大幅提高性能。如有必要的話，可以通過顯示卡加速。

3. 效果

3.1. 使用過程

程式源程式碼可參見文章最末的鏈接，是一個OpenCV結合QT的GUI的程式。編譯運行後，點擊"打開"按鈕，介面會顯示源影像：

MCV融合

點擊"繪製"按鈕，在源影像區域內繪製一個多邊形，確定一個ROI：

MCV融合

準備一張想要融合的目標影像：

MCV融合

點擊"融合"按鈕，會載入目標影像，並會根據設置的位置，將源影像的ROI融合到目標影像中：

MCV融合

3.2. 效率

在Debug模式，不使用OpenMP加速的情況下，這個演算法的效率大約需要50秒左右的時間。

在Debug模式，使用OpenMP加速，演算法的效率可以優化到10秒，也就是不使用OpenMP加速時的5倍左右。而我使用的機器CPU是i7-8750H標壓6核CPU，考慮到一些IO操作造成的性能損耗，這個優化效率是正常的。

最後在使用Release模式，使用OpenMP加速之後，演算法的效率可以優化到1秒左右，這說明編譯器優化對程式性能也是有很大影響的，尤其是對並行程式而言。

這個實現只是這個演算法的初始實現，效率就已經達到了1秒左右，看來論文說的可以達到實時融合確實不是虛言。有機會再嘗試一下論文中提到的一些性能優化實現。

4. 參考

[1] 泊松融合及其優化演算法
[2] Coordinates for Instant Image Cloning
[3] 影像處理（十二）影像融合(1)Seamless cloning泊松克隆-Siggraph 2004
[4] 多尺度並行坐標插值實時影像克隆演算法

實現程式碼

VirMach 便宜 VPS

QNews

QNews