本篇文章主要为讲解图像处理的泊松融合的原理及实现。
泊松融合原理来源于这篇文章：《Poisson Image Editing》
本人为图像处理的小白，在机缘巧合下，看到了泊松融合的图像处理，觉得很强大也十分有趣，对其中的数学原理也十分感兴趣。因此便查找了很多资料，包括原理加实现。这篇文章主要基于我自己对泊松融合的理解，再进行解释一遍，并将它进行实现，一个是帮助自己加深了解，第二个也算是作为自己的笔记。
首先，我们来了解一些基本的概念。
微分和卷积
先从一维数组的微分开始介绍
∇
\nabla
∇表示位置
x
x
x一阶微分计算（一阶中心导）:
d
f
(
x
)
d
x
=
f
(
x
+
h
)
−
f
(
x
−
h
)
2
h
\frac{df(x)}{dx}=\frac{f(x+h)-f(x-h)}{2h}
dxdf(x)=2hf(x+h)−f(x−h)
△ \triangle △表示位置 x x x二阶微分计算（二阶中心导）: d 2 f ( x ) d 2 x = f ( x + h ) − 2 f ( x ) + f ( x − h ) 2 h 2 \frac{d^2f(x)}{d^2x}=\frac{f(x+h)-2f(x)+f(x-h)}{2h^2} d2xd2f(x)=2h2f(x+h)−2f(x)+f(x−h)
当 h → 0 h\rightarrow0 h→0时，上式的微分算式的结果会逐渐逼近真实的微分值。对于图像，图像的每个像素都是离散非连续的，因此这里的 h h h放在实际的图像处理当中可看作为像素的间距，可视为1。将1代入上面的二阶微分计算式，我们可以将二阶微分的计算结果看作是一个 1 × 3 1\times3 1×3 的卷积核 [ 1 , − 2 , 1 ] [1,-2,1] [1,−2,1] 在一维度数组上进行卷积计算的结果。
当数组维度变为二维数组时，也就是图像处理的拉普拉斯算子：
△
=
∂
∂
x
2
+
∂
∂
y
2
\triangle=\frac{\partial}{\partial x^2}+\frac{\partial}{\partial y^2}
△=∂x2∂+∂y2∂
此时卷积核尺寸应该是
3
×
3
3\times3
3×3,即
[
0
0
0
1
−
2
1
0
0
0
]
+
[
0
1
0
0
−
2
0
0
1
0
]
=
[
0
1
0
1
−
4
1
0
1
0
]
\begin{bmatrix}0 & 0& 0\\1&-2&1\\0&0&0\end{bmatrix}+\begin{bmatrix}0 & 1& 0\\0&-2&0\\0&1&0\end{bmatrix}=\begin{bmatrix}0 & 1& 0\\1&-4&1\\0&1&0\end{bmatrix}
⎣⎡0100−20010⎦⎤+⎣⎡0001−21000⎦⎤=⎣⎡0101−41010⎦⎤， 称为拉普拉斯卷积核。
泊松方程的求解
已知图像每点的二阶微分值（即散度
d
i
v
div
div）,求解各个图像点的像素值。
举个例子，假设有一张
4
×
4
4\times4
4×4 的图像
X
:
[
x
1
x
2
x
3
x
4
x
5
x
6
x
7
x
8
x
9
x
10
x
11
x
12
x
13
x
14
x
15
x
16
]
X:\left[\begin{array}{llll}{x_{1}} & {x_{2}} & {x_{3}} & {x_{4}} \\ {x_{5}} & {x_{6}} & {x_{7}} & {x_{8}} \\ {x_{9}} & {x_{10}} & {x_{11}} & {x_{12}} \\ {x_{13}} & {x_{14}} & {x_{15}} & {x_{16}}\end{array}\right]
X:⎣⎢⎢⎡x1x5x9x13x2x6x10x14x3x7x11x15x4x8x12x16⎦⎥⎥⎤，
x
i
x_i
xi
表示各个位置上的图像像素值，共16个未知参数需要被求解。
应用拉普拉斯卷积核后，得到4个方程式：
{
x
2
+
x
5
+
x
7
+
x
10
−
4
x
6
=
div
x
6
x
3
+
x
6
+
x
8
+
x
11
−
4
x
7
=
div
x
7
x
6
+
x
9
+
x
11
+
x
14
−
4
x
10
=
div
x
10
x
7
+
x
10
+
x
12
+
x
15
−
4
x
11
=
div
x
11
\left\{\begin{array}{l}{x_{2}+x_{5}+x_{7}+x_{10}-4 x_{6}=\operatorname{div} x_{6}} \\ {x_{3}+x_{6}+x_{8}+x_{11}-4 x_{7}=\operatorname{div} x_{7}} \\ {x_{6}+x_{9}+x_{11}+x_{14}-4 x_{10}=\operatorname{div} x_{10}} \\ {x_{7}+x_{10}+x_{12}+x_{15}-4 x_{11}=\operatorname{div} x_{11}}\end{array}\right.
⎩⎪⎪⎨⎪⎪⎧x2+x5+x7+x10−4x6=divx6x3+x6+x8+x11−4x7=divx7x6+x9+x11+x14−4x10=divx10x7+x10+x12+x15−4x11=divx11
但是通过4个方程求解16个未知量是不可行的。但是如果边界的元素是已知的呢，即如果我们认为
x
1
,
x
2
,
x
3
,
x
4
,
x
8
,
x
12
,
x
16
,
x
15
,
x
14
,
x
13
,
x
9
,
x
5
x_1,x_2,x_3,x_4,x_8,x_{12},x_{16},x_{15},x_{14},x_{13},x_9,x_5
x1,x2,x3,x4,x8,x12,x16,x15,x14,x13,x9,x5这几个边界元素是已知的，那么就剩下4个未知量，就可以求解了。
矩阵化该方程，的此式
A
x
=
b
Ax=b
Ax=b。
一维的融合
现在我们回到图像融合这个话题，同样，我们还是先从一维看起。
我们希望将左边图的红色小块块插入到右边的？？？当中，但是右希望插入后，边界处能够衔接的自然，那这和图像融合是一个道理。下面
f
i
f_i
fi代表上图横轴上
i
i
i的方块的高度值，
f
1
=
6
f_1=6
f1=6,
f
6
=
1
f_6=1
f6=1。
{
M
i
n
(
(
f
2
−
f
1
)
−
1
)
2
M
i
n
(
(
f
3
−
f
2
)
−
(
−
1
)
)
2
M
i
n
(
(
f
4
−
f
3
)
−
2
)
2
M
i
n
(
(
f
5
−
f
4
)
−
(
−
1
)
)
2
\left\{ \begin{aligned} Min&\quad((f_2-f_1)-1)^2 \\ Min&\quad((f_3-f_2)-(-1))^2 \\ Min&\quad((f_4-f_3)-2)^2 \\ Min&\quad((f_5-f_4)-(-1))^2 \\ \end{aligned} \right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧MinMinMinMin((f2−f1)−1)2((f3−f2)−(−1))2((f4−f3)−2)2((f5−f4)−(−1))2
可以转化为下式：
Q
=
M
i
n
(
f
2
2
+
49
−
14
f
2
+
f
3
2
+
f
2
2
+
1
−
2
f
3
f
2
+
2
f
3
−
2
f
2
+
f
4
2
+
f
3
2
+
4
−
2
f
3
f
4
−
4
f
4
+
4
f
3
+
f
5
2
+
f
4
2
+
1
−
2
f
5
f
4
+
2
f
5
−
2
f
4
+
f
5
2
+
4
−
4
f
5
)
Q=Min(f^2_2+49-14f_2 \\ \quad\quad\quad\quad+f_3^2+f_2^2+1-2f_3f_2+2f_3-2f_2\\ \quad\quad\quad\quad+f_4^2+f_3^2+4-2f_3f_4-4f_4+4f_3\\ \quad\quad\quad\quad+f_5^2+f_4^2+1-2f_5f_4+2f_5-2f_4\\ \quad\quad\quad\quad+f_5^2+4-4f_5 )
Q=Min(f22+49−14f2+f32+f22+1−2f3f2+2f3−2f2+f42+f32+4−2f3f4−4f4+4f3+f52+f42+1−2f5f4+2f5−2f4+f52+4−4f5)
然后我们可以得到下面这个式子：
d
Q
d
f
2
=
2
f
2
+
2
f
2
−
2
f
3
−
16
d
Q
d
f
3
=
2
f
3
−
2
f
2
+
2
+
2
f
3
−
2
f
4
+
4
d
Q
d
f
4
=
2
f
4
−
2
f
3
−
4
+
2
f
4
−
2
f
5
−
2
d
Q
d
f
5
=
2
f
5
−
2
f
4
+
2
+
2
f
5
−
4
\begin{array}{l}{\frac{d Q}{d f_{2}}=2 f_{2}+2 f_{2}-2 f_{3}-16} \\ {\frac{d Q}{d f_{3}}=2 f_{3}-2 f_{2}+2+2 f_{3}-2 f_{4}+4} \\ {\frac{d Q}{d f_{4}}=2 f_{4}-2 f_{3}-4+2 f_{4}-2 f_{5}-2} \\ {\frac{d Q}{d f_{5}}=2 f_{5}-2 f_{4}+2+2 f_{5}-4}\end{array}
df2dQ=2f2+2f2−2f3−16df3dQ=2f3−2f2+2+2f3−2f4+4df4dQ=2f4−2f3−4+2f4−2f5−2df5dQ=2f5−2f4+2+2f5−4
转化为矩阵，即
A
x
=
b
Ax=b
Ax=b 的形式表示为：
(
4
−
2
0
0
−
2
4
−
2
0
0
−
2
4
−
2
0
0
−
2
4
)
(
f
2
f
3
f
4
f
5
)
=
(
16
−
6
6
2
)
\left(\begin{array}{cccc}{4} & {-2} & {0} & {0} \\ {-2} & {4} & {-2} & {0} \\ {0} & {-2} & {4} & {-2} \\ {0} & {0} & {-2} & {4}\end{array}\right)\left(\begin{array}{l}{f_{2}} \\ {f_{3}} \\ {f_{4}} \\ {f_{5}}\end{array}\right)=\left(\begin{array}{c}{16} \\ {-6} \\ {6} \\ {2}\end{array}\right)
⎝⎜⎜⎛4−200−24−200−24−200−24⎠⎟⎟⎞⎝⎜⎜⎛f2f3f4f5⎠⎟⎟⎞=⎝⎜⎜⎛16−662⎠⎟⎟⎞
解得
f
2
=
6
,
f
3
=
4
,
f
4
=
5
,
f
5
=
3
f_2=6, f_3=4, f_4=5, f_5=3
f2=6,f3=4,f4=5,f5=3
插入进去的效果图如下：
现在我们来看二维的问题：
有标号的像素为图像融合要插入的内容像素，红色的像素代表内容的边界。
大家现在可以回顾下上面的泊松方程求解的
4
×
4
4\times4
4×4 的图像的例子。
我们先以标号1的像素举个例子，方便大家理解下面的式子是怎么构造出来的
像 素 1 u p ∗ 1 + 像 素 1 l e f t ∗ 1 + 像 素 1 r i g h t ∗ 1 + 像 素 1 d o w n ∗ 1 − 4 ∗ 像 素 1 = d i v ( 像 素 1 ) 像素1_{up}*1+像素1_{left}*1+像素1_{right}*1+像素1_{down}*1-4*像素1=div(像素1) 像素1up∗1+像素1left∗1+像素1right∗1+像素1down∗1−4∗像素1=div(像素1)
像 素 1 u p 像素1_{up} 像素1up指的是像素1上面的像素，其他的类似。然后我们前面说了求解泊松方程时边界像素是已知的即，像素1的up、left、down的像素是已知的。
像 素 1 d o w n ∗ 1 − 4 ∗ 像 素 1 = d i v ( 像 素 1 ) − 像 素 1 u p ∗ 1 − 像 素 1 l e f t ∗ 1 − 像 素 1 r i g h t ∗ 1 像素1_{down}*1-4*像素1=div(像素1)-像素1_{up}*1-像素1_{left}*1-像素1_{right}*1 像素1down∗1−4∗像素1=div(像素1)−像素1up∗1−像素1left∗1−像素1right∗1
现在我们来创建求解矩阵的
A
x
=
b
Ax=b
Ax=b中的
A
A
A,
x
x
x,
b
b
b。
x
=
[
像
素
1
，
像
素
2
，
像
素
3
，
⋯
,
像
素
10
]
x=[像素1，像素2，像素3，\cdots,像素10]
x=[像素1，像素2，像素3，⋯,像素10]
A
为
10
×
10
的
矩
阵
，
创
建
方
式
如
下
：
A为10\times10的矩阵，创建方式如下：
A为10×10的矩阵，创建方式如下：
创建出的矩阵如下：
上面的伪代码意思是，对角线的元素都为-4，如果行与列的俩元素是相邻元素则为1，比如5和2是相邻元素，在(5,2)和(2,5)位置都为1。
b
b
b的值为，
b
[
i
]
=
d
i
v
(
像
素
i
)
−
（
像
素
i
相
邻
已
知
像
素
值
的
和
）
b[i]=div(像素i)-（像素i相邻已知像素值的和）
b[i]=div(像素i)−（像素i相邻已知像素值的和）
这样就创建好了
A
x
=
b
Ax=b
Ax=b，即
x
=
A
−
1
∗
b
x=A^{-1}*b
x=A−1∗b
代码具体实现
#include "pch.h"
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include "vector"
#include "time.h"
#define elif else if
#define ATD at<double>
#define vector vector<Mat>
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<height; i++){
for(int j=0; j<width; 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<roiheight; i++){
for(int j=0; j<roiwidth; 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<roiheight; i++){
for(int j=0; j<roiwidth; 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..."<<endl;
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..."<<endl;
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..."<<endl;
// 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"<<endl;
imwrite("result/result.png", in2);
imshow("roi", result);
imshow("result", in2);
waitKey(0);
return 0;
}效果图
上面的为结果图，由于用的是最简单的方法去求解泊松方程，所以速度会比较慢，建议大家用比较小的图片去尝试，泊松方程求解的加速方法有Jacobi, SOR, Conjugate Gradients, and FFT
以上内容主要参考如下俩个博客的内容，然后再进行了一定的整理。
http://eric-yuan.me/poisson-blending/
https://zhuanlan.zhihu.com/p/68349210
8216
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
