点击上方"蓝色小字"关注我呀
本文不是原创哦，是经过知乎作者【蒸糕】授权转载。
原文链接：
https://zhuanlan.zhihu.com/p/96777721
正文
本篇文章主要为讲解图像处理的泊松融合的原理及实现。
泊松融合原理来源于这篇文章：《Poisson Image Editing》
本人为图像处理的小白，在机缘巧合下，看到了泊松融合的图像处理，觉得很强大也十分有趣，对其中的数学原理也十分感兴趣。因此便查找了很多资料，包括原理加实现。这篇文章主要基于我自己对泊松融合的理解，再进行解释一遍，并将它进行实现，一个是帮助自己加深了解，第二个也算是作为自己的笔记。
首先，我们来了解一些基本的概念。
微分和卷积
先从一维数组的微分开始介绍
当 h 趋向0时，上式的微分算式的结果会逐渐逼近真实的微分值。对于图像，图像的每个像素都是离散非连续的，因此这里的 h 放在实际的图像处理当中可看作为像素的间距，可视为1。
将1代入上面的二阶微分计算式，我们可以将二阶微分的计算结果看作是一个 1×3 的卷积核 [1,-2,1] 在一维度数组上进行卷积计算的结果。
当数组维度变为二维数组时，也就是图像处理的拉普拉斯算子：
此时卷积核尺寸应该是3×3 ,即
称为拉普拉斯卷积核。
泊松方程的求解
已知图像每点的二阶微分值（即散度 div）,求解各个图像点的像素值。
举个例子，假设有一张 4×4 的图像
Xi表示各个位置上的图像像素值，共16个未知参数需要被求解。
应用拉普拉斯卷积核后，得到4个方程式：
但是通过4个方程求解16个未知量是不可行的。但是如果边界的元素是已知的呢，即如果我们认为
【x1,x2,x3,x4,x8,x12,x16,x15,x14,x13,x9,x5】
这几个边界元素是已知的，那么就剩下4个未知量，就可以求解了。
矩阵化该方程，得此式 Ax = b
一维的融合
现在我们回到图像融合这个话题，同样，我们还是先从一维看起。
我们希望将左边图的红色小块块插入到右边的？？？当中，但是又希望插入后，边界处能够衔接的自然，那这和图像融合是一个道理。下面 fi 代表上图横轴上 i 的方块的高值， f1=6,f6=1 。
可以转化为下式：
然后我们可以得到下面这个式子：
转化为矩阵，即Ax = b的形式表示为：
解得 f2=6,f3=4,f4=5,f5=3
插入进去的效果图如下：
现在我们来看二维的问题：
有标号的像素为图像融合要插入的内容像素，红色的像素代表内容的边界。
大家现在可以回顾下上面的泊松方程求解的 [公式] 的图像的例子。
我们先以标号1的像素举个例子，方便大家理解下面的式子是怎么构造出来的
像素1up指的是像素1上面的像素，其他的类似。然后我们前面说了求解泊松方程时边界像素是已知的即，像素1的up、left、down的像素是已知的。
现在我们来创建求解矩阵的 Ax=b中的A,x,b。
上面的伪代码意思是，对角线的元素都为-4，如果行与列的俩元素是相邻元素则为1，比如5和2是相邻元素，在(5,2)和(2,5)位置都为1。
代码具体实现
代码语言：javascript复制#include "pch.h"
#include
#include
#include
#include
#include "vector"
#include "time.h"
#define elif else if
#define ATD at
#define vector vector
using namespace cv;
using namespace std;
//calculate horizontal gradient, img(i,j+1) - img(i,j)
Mat getGradientXp(Mat &img)
{
int height = img.rows;
int width = img.cols;
Mat cat = repeat(img, 1, 2);
Rect roi = Rect(1, 0, width, height);
Mat roimat = cat(roi);
return roimat - img;
}
//calculate vertical gradient, img(i+1,j) - img(i,j)
Mat getGradientYp(Mat &img)
{
int height = img.rows;
int width = img.cols;
Mat cat = repeat(img, 2, 1);
Rect roi = Rect(0, 1, width, height);
Mat roimat = cat(roi);
return roimat - img;
}
//calculate horizontal gradient, img(i,j-1) - img(i,j)
Mat getGradientXn(Mat &img)
{
int height = img.rows;
int width = img.cols;
Mat cat = repeat(img, 1, 2);
Rect roi = Rect(width-1, 0, width, height);
Mat roimat = cat(roi);
return roimat - img;
}
//calculate vertical gradient, img(i-1,j) - img(i,j)
Mat getGradientYn(Mat &img)
{
int height = img.rows;
int width = img.cols;
Mat cat = repeat(img, 2, 1);
Rect roi = Rect(0, height-1, width, height);
Mat roimat = cat(roi);
return roimat - img;
}
int getLabel(int i, int j, int height, int width)
{
return i * width + j;
}
//get Matrix A.
Mat getA(int height, int width)
{
Mat A = Mat::eye(height*width, height*width, CV_64FC1);
A *= -4;
Mat M = Mat::zeros(height, width, CV_64FC1);
Mat temp = Mat::ones(height, width - 2, CV_64FC1);
Rect roi = Rect(1, 0, width - 2, height);
Mat roimat = M(roi);
temp.copyTo(roimat);
temp = Mat::ones(height - 2, width, CV_64FC1);
roi = Rect(0, 1, width, height - 2);
roimat = M(roi);
temp.copyTo(roimat);
temp = Mat::ones(height - 2, width - 2, CV_64FC1);
temp *= 2;
roi = Rect(1, 1, width - 2, height - 2);
roimat = M(roi);
temp.copyTo(roimat);
for(int i=0; i for(int j=0; j int label = getLabel(i, j, height, width); if(M.ATD(i, j) == 0){ if(i == 0) A.ATD(getLabel(i + 1, j, height, width), label) = 1; elif(i == height - 1) A.ATD(getLabel(i - 1, j, height, width), label) = 1; if(j == 0) A.ATD(getLabel(i, j + 1, height, width), label) = 1; elif(j == width - 1) A.ATD(getLabel(i, j - 1, height, width), label) = 1; }elif(M.ATD(i, j) == 1){ if(i == 0){ A.ATD(getLabel(i + 1, j, height, width), label) = 1; A.ATD(getLabel(i, j - 1, height, width), label) = 1; A.ATD(getLabel(i, j + 1, height, width), label) = 1; }elif(i == height - 1){ A.ATD(getLabel(i - 1, j, height, width), label) = 1; A.ATD(getLabel(i, j - 1, height, width), label) = 1; A.ATD(getLabel(i, j + 1, height, width), label) = 1; } if(j == 0){ A.ATD(getLabel(i, j + 1, height, width), label) = 1; A.ATD(getLabel(i - 1, j, height, width), label) = 1; A.ATD(getLabel(i + 1, j, height, width), label) = 1; }elif(j == width - 1){ A.ATD(getLabel(i, j - 1, height, width), label) = 1; A.ATD(getLabel(i - 1, j, height, width), label) = 1; A.ATD(getLabel(i + 1, j, height, width), label) = 1; } }else{ A.ATD(getLabel(i, j - 1, height, width), label) = 1; A.ATD(getLabel(i, j + 1, height, width), label) = 1; A.ATD(getLabel(i - 1, j, height, width), label) = 1; A.ATD(getLabel(i + 1, j, height, width), label) = 1; } } } return A; } // Get the following Laplacian matrix // 0 1 0 // 1 -4 1 // 0 1 0 Mat getLaplacian(){ Mat laplacian = Mat::zeros(3, 3, CV_64FC1); laplacian.ATD(0, 1) = 1.0; laplacian.ATD(1, 0) = 1.0; laplacian.ATD(1, 2) = 1.0; laplacian.ATD(2, 1) = 1.0; laplacian.ATD(1, 1) = -4.0; return laplacian; } // Calculate b // using convolution. Mat getB1(Mat &img1, Mat &img2, int posX, int posY, Rect ROI){ Mat Lap; filter2D(img1, Lap, -1, getLaplacian()); int roiheight = ROI.height; int roiwidth = ROI.width; Mat B = Mat::zeros(roiheight * roiwidth, 1, CV_64FC1); for(int i=0; i for(int j=0; j double temp = 0.0; temp += Lap.ATD(i + ROI.y, j + ROI.x); if(i == 0) temp -= img2.ATD(i - 1 + posY, j + posX); if(i == roiheight - 1) temp -= img2.ATD(i + 1 + posY, j + posX); if(j == 0) temp -= img2.ATD(i + posY, j - 1 + posX); if(j == roiwidth - 1) temp -= img2.ATD(i + posY, j + 1 + posX); B.ATD(getLabel(i, j, roiheight, roiwidth), 0) = temp; } } return B; } // Calculate b // using getGradient functions. Mat getB2(Mat &img1, Mat &img2, int posX, int posY, Rect ROI){ Mat grad = getGradientXp(img1) + getGradientYp(img1) + getGradientXn(img1) + getGradientYn(img1); int roiheight = ROI.height; int roiwidth = ROI.width; Mat B = Mat::zeros(roiheight * roiwidth, 1, CV_64FC1); for(int i=0; i for(int j=0; j double temp = 0.0; temp += grad.ATD(i + ROI.y, j + ROI.x); if(i == 0) temp -= img2.ATD(i - 1 + posY, j + posX); if(i == roiheight - 1) temp -= img2.ATD(i + 1 + posY, j + posX); if(j == 0) temp -= img2.ATD(i + posY, j - 1 + posX); if(j == roiwidth - 1) temp -= img2.ATD(i + posY, j + 1 + posX); B.ATD(getLabel(i, j, roiheight, roiwidth), 0) = temp; } } return B; } // Solve equation and reshape it back to the right height and width. Mat getResult(Mat &A, Mat &B, Rect &ROI){ Mat result; solve(A, B, result); result = result.reshape(0, ROI.height); return result; } // img1: 3-channel image, we wanna move something in it into img2. // img2: 3-channel image, dst image. // ROI: the position and size of the block we want to move in img1. // posX, posY: where we want to move the block to in img2 Mat poisson_blending(Mat &img1, Mat &img2, Rect ROI, int posX, int posY){ int roiheight = ROI.height; int roiwidth = ROI.width; Mat A = getA(roiheight, roiwidth); // we must do the poisson blending to each channel. vector rgb1; split(img1, rgb1); vector rgb2; split(img2, rgb2); vector result; Mat merged, res, Br, Bg, Bb; // For calculating B, you can use either getB1() or getB2() Br = getB2(rgb1[0], rgb2[0], posX, posY, ROI); //Br = getB2(rgb1[0], rgb2[0], posX, posY, ROI); res = getResult(A, Br, ROI); result.push_back(res); cout<<"R channel finished..."< Bg = getB2(rgb1[1], rgb2[1], posX, posY, ROI); //Bg = getB2(rgb1[1], rgb2[1], posX, posY, ROI); res = getResult(A, Bg, ROI); result.push_back(res); cout<<"G channel finished..."< Bb = getB2(rgb1[2], rgb2[2], posX, posY, ROI); //Bb = getB2(rgb1[2], rgb2[2], posX, posY, ROI); res = getResult(A, Bb, ROI); result.push_back(res); cout<<"B channel finished..."< // merge the 3 gray images into a 3-channel image merge(result,merged); return merged; } int main(int argc, char** argv) { long start, end; start = clock(); Mat img1, img2; Mat in1 = imread("pic/airplane2.png"); Mat in2 = imread("pic/background1.png"); imshow("src", in1); imshow("dst", in2); in1.convertTo(img1, CV_64FC3); in2.convertTo(img2, CV_64FC3); int posXinPic2 = 350; int posYinPic2 = 50; Rect rc = Rect(0, 0, in1.cols, in1.rows); Mat result = poisson_blending(img1, img2, rc, posXinPic2, posYinPic2); result.convertTo(result, CV_8UC1); Rect rc2 = Rect(posXinPic2, posYinPic2, in1.cols, in1.rows); Mat roimat = in2(rc2); result.copyTo(roimat); end = clock(); cout<<"used time: "<<((double)(end - start)) / CLOCKS_PER_SEC<<" second"< imwrite("result/result.png", in2); imshow("roi", result); imshow("result", in2); waitKey(0); return 0; }效果图 上面的为结果图，由于用的是最简单的方法去求解泊松方程，所以速度会比较慢，建议大家用比较小的图片去尝试，泊松方程求解的加速方法有Jacobi, SOR, Conjugate Gradients, 和 FFT。 推荐阅读： 代码语言：javascript复制https://zhuanlan.zhihu.com/p/68349210THE END 我们有一个视觉算法交流群哦，加我微信备注【加群】，一起来学习。