1. 概述
我在之前的文章《基於均值坐標(Mean-Value Coordinates)的圖像融合算法的具體實現》中，根據《Coordinates for Instant Image Cloning》這篇論文，詳細論述了圖像融合中泊松融合算法的優化算法——均值坐標(Mean-Value Coordinates)融合算法的具體實現。其實在這篇論文中，還提出了兩種優化實現，能夠進一步提升效率，這里就論述一下其優化算法的具體實現。
2. 實現
2.1. 原理
均值坐標融合算法的核心思想是算出ROI中每個點的MVC(Mean-Value Coordinates)，如果ROI中存在n個點，ROI邊界像素為m個點，那么該算法的時間復雜度至少為O(nm)。根據《Coordinates for Instant Image Cloning》的描述，MVC融合算法修正值其實是一個線性區間，只要得到其中一些關鍵點的融合修正值，其他點的融合修正值就可以根據周圍的融合修正值線性插值出來。
因此，可以通過CGAL來對ROI多邊形邊界構建一個自適應三角網，以邊界上每個柵格點作為約束構網，為了滿足Delaunay特性，就會在ROI內部新添加一些點，這樣就會出現邊界小而密集，邊界大而稀疏的自適應三角網(可參看這篇文章《通過CGAL將一個多邊形剖分成Delaunay三角網》)：
這樣，n個點就會將為常數級別個數c，也就是時間復雜度至少為O(cm)。當然從三角面上插值也會有時間消耗，但時間復雜度會明顯小於計算均值坐標。
2.2. 核心代碼
具體核心代碼實現如下（完整的代碼實現地址見文章末尾）：
//三角網優化
void QImageShowWidget::MVCBlendOptimize(int posX, int posY)
{
QTime startTime = QTime::currentTime();
//Step1:找到邊界上所有的像素點
vector<Vector2d> ROIBoundPointList;
CalBoundPoint(ROIBoundPointList);
//
CDT cdt;
vector<Vertex_handle> vertexList;
for(int i = 0; i<ROIBoundPointList.size(); i++)
{
//cout<<ROIBoundPointList[i].x<<','<<ROIBoundPointList[i].y<<'\t';
//vertexList.push_back(cdt.insert(Point(pointList[i].x(), pointList[i].y() )));
vertexList.push_back(cdt.insert(CDTPoint(ROIBoundPointList[i].x, ROIBoundPointList[i].y )));
}
for(unsigned int i =0;i<vertexList.size()-1;i++)
{
cdt.insert_constraint(vertexList[i],vertexList[i+1]);
}
std::cout << "Number of vertices: " << cdt.number_of_vertices() <<std::endl;
std::cout << "Meshing the triangulation..." << std::endl;
CGAL::refine_Delaunay_mesh_2(cdt, Criteria());
std::cout << "Number of vertices: " << cdt.number_of_vertices() <<std::endl;
vector<Vector2d> vertex_list;
map<Vector2d, size_t> vertex_map;
for(CDT::Vertex_iterator vit = cdt.vertices_begin(); vit!= cdt.vertices_end(); ++vit)
{
vertex_map.insert(make_pair(Vector2d(vit->point().x(), vit->point().y()), vertex_list.size()));
vertex_list.push_back(Vector2d(vit->point().x(), vit->point().y()));
}
//計算邊界的像素差值
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);
//rMat.data[n] = d;
}
//clipMap[l] = false; //在多邊形邊上的點沒法計算MVC
}
cout<<"開始計算 mean-value coordinates..." << endl;
vector<Vec3d> tri_mesh_vertex_R(vertex_list.size());
#pragma omp parallel for //開啟OpenMP並行加速
for (int vi = 0; vi < vertex_list.size(); ++vi)
{
//逐點計算MVC
vector<double> alphaAngle(ROIBoundPointList.size());
for(size_t pi = 1; pi < ROIBoundPointList.size(); pi++)
{
alphaAngle[pi] = threePointCalAngle(ROIBoundPointList[pi-1], vertex_list[vi], ROIBoundPointList[pi]);
}
alphaAngle[0] = alphaAngle[ROIBoundPointList.size()-1];
vector<double> MVC(ROIBoundPointList.size()-1, 0);
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] - vertex_list[vi]).Mod();
MVC[pi-1] = w_a / w_b;
if(_isnan(MVC[pi-1])==1)
{
MVC[pi-1] = 0;
}
}
double sum = 0;
for(size_t pi = 0; pi < MVC.size(); pi++)
{
sum = sum + MVC[pi];
}
for(size_t pi = 0; pi < MVC.size(); pi++)
{
MVC[pi] = MVC[pi] / sum;
}
Vec3d r(0.0,0.0,0.0);
for(size_t pi = 0; pi < MVC.size(); pi++)
{
for(int bi = 0; bi < winBandNum; bi++)
{
r[bi] = r[bi] + MVC[pi] * diff[pi * winBandNum + bi];
}
}
tri_mesh_vertex_R[vi] = r;
}
cout<<"計算完成！" << endl;
//遍歷每一個三角面
vector<vector<size_t>> face_vertex_index;
CDT::Face_iterator fit;
for (fit = cdt.faces_begin(); fit!= cdt.faces_end(); ++fit)
{
vector<size_t> index(3);
for(int i = 0; i<3; i++)
{
auto iter = vertex_map.find(Vector2d(fit->vertex(i)->point().x(), fit->vertex(i)->point().y()));
if(iter == vertex_map.end())
{
continue;
}
index[i] = iter->second;
}
face_vertex_index.push_back(index);
}
size_t srcImgBufNum = static_cast<size_t>(srcImg.cols) * static_cast<size_t>(srcImg.rows);
vector<size_t> clipMap(srcImgBufNum, 0); //標識范圍內的點: 0標識初始不能寫入，1以上標識在那個三角形
#pragma omp parallel for //開啟OpenMP並行加速
for(int fi = 0; fi < face_vertex_index.size(); fi++)
{
Vector2d v0 = vertex_list[face_vertex_index[fi][0]];
Vector2d v1 = vertex_list[face_vertex_index[fi][1]];
Vector2d v2 = vertex_list[face_vertex_index[fi][2]];
double minX = std::min(std::min(v0.x, v1.x), v2.x);
double minY = std::min(std::min(v0.y, v1.y), v2.y);
double maxX = std::max(std::max(v0.x, v1.x), v2.x);
double maxY = std::max(std::max(v0.y, v1.y), v2.y);
int sX = std::max(int(floor(minX)), 0);
int sY = std::max(int(floor(minY)), 0);
int eX = std::max(int(ceil(maxX)), srcImg.cols - 1);
int eY = std::max(int(ceil(maxY)), srcImg.rows - 1);
for(int yi = sY; yi <= eY; yi++)
{
for(int xi = sX; xi <= eX; xi++)
{
if(PointinTriangle(Vector3d(v0), Vector3d(v1), Vector3d(v2), Vector3d(xi, yi, 0)))
{
size_t m = static_cast<size_t>(srcImg.cols) * static_cast<size_t>(yi) + xi;
clipMap[m] = fi+1;
}
}
}
}
cout<<"開始插值計算..." << endl;
//Mat result(srcImg.rows, srcImg.cols, CV_8UC1);
#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] == 0)
{
continue;
}
if(!Point_In_Polygon_2D(ci, ri, ROIBoundPointList))
{
continue;
}
size_t fi = clipMap[l]-1;
size_t index0 = face_vertex_index[fi][0];
size_t index1 = face_vertex_index[fi][1];
size_t index2 = face_vertex_index[fi][2];
vector<double> r(winBandNum, 0);
for(int bi = 0; bi < winBandNum; bi++)
{
Vector3d p0(vertex_list[index0].x, vertex_list[index0].y, tri_mesh_vertex_R[index0][bi]);
Vector3d p1(vertex_list[index1].x, vertex_list[index1].y, tri_mesh_vertex_R[index1][bi]);
Vector3d p2(vertex_list[index2].x, vertex_list[index2].y, tri_mesh_vertex_R[index2][bi]);
Vector3d vp(ci, ri, 0);
CalPlanePointZ(p0, p1, p2, vp);
r[bi] = vp.z;
}
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);
}
}
}
//imwrite("D:/result.tif", result);
cout<<"插值完成！" << endl;
QTime stopTime = QTime::currentTime();
int elapsed = startTime.msecsTo(stopTime);
cout<<"總結完成用時"<<elapsed<<"毫秒";
}
主要思路還是通過ROI多邊形柵格建立三角網，計算網格點的MVC，繼而計算融合修正值；而其他點的融合修正值則通過所在三角形頂點的融合修正值插值得到。
注意這里麻煩的地方是還得計算每個點是在那個三角形內，我這里是采取索引數組的辦法。如果直接采取遍歷每個點與每個三角形的辦法，那么時間復雜度可能會超過計算MVC的復雜度。而插值的算法可以參考這篇文章《已知空間三點組成的面求該面上某點的Z值》。
2.3. 第二種優化
在《Coordinates for Instant Image Cloning》這篇論文中，還介紹了第二種優化算法。算法思想是除了減少ROI內的點，還可以減少ROI邊界上的點：每個點的MVC值其實可以不用到邊界上所有的點，可以通過一種規則算法來指定需要的點。可惜這個規則算法我也沒看明白，有機會再進一步研究。
3. 結果
融合的源影像：
融合的目標影像：
融合的結果：
運行的時間：
這里可以看到，優化后的融合效率為501毫秒，而優化之前的效率為1秒，效率提升了50%。