Java模拟泊松分布

一、论文介绍

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

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

二、一些先验知识

1.泊松方程

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

2.Laplace算子

在二维图像中，图像函数$f$是离散函数，拉普拉斯算子为如下卷积核形式：
$\begin{bmatrix} 0& 1& 0\\ 1& -4& 1\\ 0& 1& 0 \end{bmatrix}$
比如要计算$(i,j)$处的像素点的二阶微分：
$\bigtriangleup f(i,j)=\begin{bmatrix} 0& 1& 0\\ 1& -4& 1\\ 0& 1& 0 \end{bmatrix}.\begin{bmatrix} 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{bmatrix}$

三、算法原理

算法的目的是将一张图像的某一个区域与另一张图像更好地融合，当两张图像梯度（图像灰度的变化率）越接近，则说明两张图像融合得越好。即最小化如下公式：

$min\iint \left | \bigtriangledown f-v \right | ^{2}$

其中$v$是源图像的梯度，$\bigtriangledown f$是目标图像的梯度。

要使得该公式最小，即$\bigtriangleup f=\bigtriangledown v$

如上图所示，将源图像中脸部区域复制到目标图像中，经过简单的复制，得到右图所示的cloing图像，这个过程可以抽象为如下图所示的模型。

图中，$g$是源图像圈出区域的图像函数，$v$是设置的与源图像有关的梯度向量场，

$v=\bigtriangledown g$

则$\bigtriangledown v=\bigtriangleup f$.

$\Omega$是从源图像中抠出的部分区域，$f$是该区域像素的像素值；$S$是目标图像中除了$\Omega$区域外的部分，$f*$是该区域像素的像素值。$\partial \Omega$是$\Omega$区域的边界，边界像素值是对应的目标图像该处的像素值：

$f_{\mid \partial \Omega}=f_{\mid \partial \Omega}^{*}$

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

四、一个例子

如下图所示，左图表示$4*4$的源图像，右图是$4*4$的目标图像，将左图中矩形框框住的区域复制到右图白色区域中，源图像$g_{1},...,g_{16}$已知，目标图像中$f_{1},...,f_{16}$已知，要求的是$x_{1},x_{2},x_{3},x_{4}$

根据以上列出的泊松方程列出如下式子：

$-4 x_{1}+x_{2}+x_{3}+f_{2}+f_{5}=-4 g_{6}+g_{2}+g_{5}+g_{7}+g_{10}$

$-4 x_{2}+x_{1}+x_{4}+f_{3}+f_{8}=-4 g_{7}+g_{3}+g_{6}+g_{8}+g_{11}$

$-4 x_{3}+x_{1}+x_{4}+f_{9}+f_{14}=-4 g_{10}+g_{6}+g_{9}+g_{11}+g_{14}$

$-4 x_{4}+x_{2}+x_{3}+f_{12}+f_{15}=-4 g_{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} -4 g_{6}+g_{2}+g_{5}+g_{7}+g_{10}-f_{2}-f_{5} \\ -4 g_{7}+g_{3}+g_{6}+g_{8}+g_{11}-f_{3}-f_{8} \\ -4 g_{10}+g_{6}+g_{9}+g_{11}+g_{14}-f_{9}-f_{14} \\ -4 g_{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;
}

六、效果图

七、总结