论文复现：Poisson Image Editing

一、论文介绍

Poisson Image Editing这篇论文使用泊松方程，实现了一种图像融合算法，将一张图像的某个区域无缝地复制到另一张图像中。如下图所示。

可以看到cloning是经过简单复制得到的图像，从源图像中抠出的区域并不能与目标图像很好地融合。右图是使用本文中提出的算法，经过无缝复制后，从源图像的抠出的区域可以无缝地导入目标图像中。

二、一些先验知识

1.泊松方程

泊松方程是泊松在物理学领域提出的一个偏微分方程，其形式为：
$$△f=\Omega$$
其中$△$表示二阶微分即拉普拉斯算子，$\Omega$为已知量。
对于泊松方程的求解，论文中提到了给定狄利克雷边界条件的泊松方程解法，即给定函数$f$在边界处的实际值。

2.Laplace算子

在二维图像中，图像函数$f$是离散函数，拉普拉斯算子为如下卷积核形式：
$$\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]$$
比如要计算$(i,j)$处的像素点的二阶微分：
$$△f(i,j)=\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right].\left[\begin{array}{ccc}f(i-1,j-1)& f(i-1,j)& f(i-1,j+1)\\ f(i,j-1)& f(i,j)& f(i,j+1)\\ f(i+1,j-1)& f(i+1,j)& f(i+1,j+1)\end{array}\right]$$

三、算法原理

算法的目的是将一张图像的某一个区域与另一张图像更好地融合，当两张图像梯度（图像灰度的变化率）越接近，则说明两张图像融合得越好。即最小化如下公式：
$$min\iint {\left|▽f-v\right|}^2$$
其中$v$是源图像的梯度，$▽f$是目标图像的梯度。
要使得该公式最小，即$$△f=▽v$$
如上图所示，将源图像中脸部区域复制到目标图像中，经过简单的复制，得到右图所示的cloing图像，这个过程可以抽象为如下图所示的模型。

图中，$g$是源图像圈出区域的图像函数，$v$是设置的与源图像有关的梯度向量场，
$$v=▽g$$
则$▽v=△f$.
$\Omega$是从源图像中抠出的部分区域，$f$是该区域像素的像素值；$S$是目标图像中除了$\Omega$区域外的部分，$f\ast$是该区域像素的像素值。$\partial \Omega$是$\Omega$区域的边界，边界像素值是对应的目标图像该处的像素值：
$$f_{\mid \partial \Omega }=f_{\mid \partial \Omega }^{\ast }$$

这样就构造了一个具有狄利克雷边界条件的泊松方程：$△f=△g$ $with$ $f_{\mid \partial \Omega }=f_{\mid \partial \Omega }^{\ast }$

四、一个例子

如下图所示，左图表示$4\ast 4$的源图像，右图是$4\ast 4$的目标图像，将左图中矩形框框住的区域复制到右图白色区域中，源图像$g_1,...,g_{16}$已知，目标图像中$f_1,...,f_{16}$已知，要求的是$x_1,x_2,x_3,x_4$

根据以上列出的泊松方程列出如下式子：
$-4x_1+x_2+x_3+f_2+f_5=-4g_6+g_2+g_5+g_7+g_{10}$
$-4x_2+x_1+x_4+f_3+f_8=-4g_7+g_3+g_6+g_8+g_{11}$
$-4x_3+x_1+x_4+f_9+f_{14}=-4g_{10}+g_6+g_9+g_{11}+g_{14}$
$-4x_4+x_2+x_3+f_{12}+f_{15}=-4g_{11}+g_7+g_{10}+g_{12}+g_{15}$

经过简单移项可写成以下矩阵形式：
$$\left[\begin{array}{cccc}-4& 1& 1& 0\\ 1& -4& 0& 1\\ 1& 0& -4& 1\\ 0& 1& 1& -4\end{array}\right]\left[\begin{array}{l}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right]=\left[\begin{array}{c}-4g_6+g_2+g_5+g_7+g_{10}-f_2-f_5\\ -4g_7+g_3+g_6+g_8+g_{11}-f_3-f_8\\ -4g_{10}+g_6+g_9+g_{11}+g_{14}-f_9-f_{14}\\ -4g_{11}+g_7+g_{10}+g_{12}+g_{15}-f_{12}-f_{15}\end{array}\right]$$
4个方程式解4个未知数，可得到唯一解。

五、具体实现

代码核心在于构建线性系统$Ax=b$

int o_start_r = 32, o_start_c = 14, d_start_r = 76, d_start_c = 550;
int h = 64, w = 200;

void BuildSystem() {
	SparseMatrix<double> A(h * w, h * w);
	for (int i = 0; i < h; i++) {
		for (int j = 0; j < w; j++) {
			int index = i * w + j;
			A.coeffRef(index, index) = -4;
			if (i != 0)		A.coeffRef(index, (i - 1) * w + j) = 1;
			if (i != h - 1) A.coeffRef(index, (i + 1) * w + j) = 1;
			if (j != 0)		A.coeffRef(index, i * w + j - 1) = 1;
			if (j != w - 1)	A.coeffRef(index, i * w + j + 1) = 1;
		}
	}

	Image o_img("./data/original.jpg");
	Image d_img("./data/destination.png");

	VectorXd R(h * w); VectorXd G(h * w); VectorXd B(h * w);
	R.setZero(); G.setZero(); B.setZero();

	for (int i = 0; i < h; i++) {
		for (int j = 0; j < w; j++) {
			int o_r = o_start_r + i, o_c = o_start_c + j;
			int d_r = d_start_r + i, d_c = d_start_c + j;

			Vec3d div_RGB = o_img.GetDivAt(o_r, o_c);
			if (j == 0)
				div_RGB -= d_img.colors[d_r][d_c - 1];
			if (j == w - 1)
				div_RGB -= d_img.colors[d_r][d_c + 1];
			if (i == 0)
				div_RGB -= d_img.colors[d_r - 1][d_c];
			if (i == h - 1)
				div_RGB -= d_img.colors[d_r + 1][d_c];

			R(i * w + j) = div_RGB[0];
			G(i * w + j) = div_RGB[1];
			B(i * w + j) = div_RGB[2];
		}
	}

	SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > solver;
	solver.compute(A);
	VectorXd r = solver.solve(R);
	VectorXd g = solver.solve(G);
	VectorXd b = solver.solve(B);

	for (int i = 0; i < h * w; i++) {
		int r_ = d_start_r + i / w;
		int c_ = d_start_c + i % w;
		//保证RGB值在0~255
		double new_r = min(max(r(i), 0.0), 255.0);
		double new_g = min(max(g(i), 0.0), 255.0);
		double new_b = min(max(b(i), 0.0), 255.0);
		d_img.colors[r_][c_] = Vec3d(new_r, new_g, new_b);
	}

	d_img.Save("./data/result.jpg");
	waitKey(0);
}

int main() {
	BuildSystem();
	return 0;
}

六、效果图

七、总结