基于偏微分方程的图像修复数值模型

Numerical Model of Image Inpainting Based on Partial Differential Equation

摘要: 图像修复是图像处理的一个分支，是利用周围区域的信息填充图像的缺失或损坏区域。本文为今后能够更好地对破损图像进行修复与技术处理，进而针对图像修复中常见的边缘问题总结了几类基于偏微分方程进行图像修复的数值模型，这几种数值模型能够在保证修复效果较好且稳定的同时减少计算量，并提高修复效率，进一步完善修复效果，能够更好地为图像修复发展研究提供理论与数值模型基础。

Abstract: Image inpainting is a branch of image processing, which uses the information of the surrounding area to fill in the missing or damaged area of the image. In order to better repair and process the damaged image in the future, this paper summarizes several kinds of numerical models based on partial differential equation for image restoration, aiming at the common edge problems in image restoration. These numerical models can ensure the good and stable repair effect, reduce the amount of calculation, improve the repair efficiency, and further improve the repair effect, it can provide a theoretical and numerical model foundation for the development of image inpainting.

1. 引言

图像修复技术是当前图像处理和计算机视觉研究的热点之一，在实际生活和科学研究中有着十分重要的应用。比如被广泛应用于文物修复、图像或视频后期处理中，在计算机、天文学、生物学等方面也都有很高的应用价值 [1] [2] [3]。图像修复实际是根据图像受损区域周围收集到的信息来填充图像的缺失部分。

近年来，图像修复技术的应用范围十分广泛，发展也十分迅速。许多学者对图像修复及其数值模型进行了研究。长久以来发展和优化了一系列用于图像修复的数值方法。包括基于纹理的图像修复算法 [4] 以及非纹理图像修复技术即基于偏微分方程(partial differential equation, PDE)的修复算法 [5]，全变分(Total Variation, TV)修复算法 [6]，快速图像修复算法 [7]，凸分裂 [8] 以及应用于二值图像修复的改进Cahn-Hilliard [9] 模型等，对于图像修复处理都有大量研究。但是，目前存在的部分图像修复技术仍存在一些不足之处。例如TV模型在图像平滑区域容易产生阶梯效应。本文将总结几类现有的高效稳定、且修复效果较好的图像修复数值模型。

2. 几类图像修复的数值模型

2.1. 改进的TV模型图像修复方法

CHAN和SHEN建立了基于能量最小化原则的统一修复模型，该模型能量函数的解是利用变分原理转化的，因此被称为整体变分(TV)模型 [10]，被应用于图像修复领域，并取得了良好效果。随后Chan等人在此基础上，建立起基于TV模型的图像修复模型。因此，TV模型是图像修复领域中最经典、应用最广泛的模型之一。首先引入一个以梯度模值作为参数的扩散函数，此类函数在小梯度区域，扩散函数取值较大，增强图像的平滑度；在大梯度区域，扩散函数取值较小，降低图像的平滑度，保护了图像的边缘信息，也提高了图像的修复速度。基于TV模型的代价函数为

$R\left(u\right)={\displaystyle \underset{\Omega }{\int }r\left|\nabla u\left(x,y,t\right)\right|\text{d}x\text{d}y},$ (1)

其中D为图像的待修补区域； $\Omega$ 为图像的全区域； $u^0=u\left(x,y,0\right)$ 为待修复图像的初始值； $u\left(x,y,t\right)$ 为t时刻修复后图像的像素值； $\nabla$ 是梯度算子；r为非负实数。需要选择一个恰当的r的函数，使得待修复区域以及修补边缘尽可能平滑。该式在修复过程中为保证修补区域及其边缘尽可能平滑并对噪声有良好处理效果，应满足下式：

$\frac{1}{\left|\Omega \backslash D\right|}{\displaystyle \underset{\Omega \backslash D}{\int }{\left|u-u^0\right|}^2\text{d}x\text{d}y}={\sigma }^2,$ (2)

其中 $\left|\Omega \backslash D\right|$ 表示有效信息区域， $\sigma$ 表示噪声的标准差。运用Lagrange乘子法将以上两个式子转化为无约束条件的极值问题，其新代价函数为

$E\left(u\right)={\displaystyle \underset{\Omega }{\int }\left|\nabla u\right|\text{d}x\text{d}y}+\frac{\lambda }{2}{\displaystyle \underset{\Omega }{\int }{\left|u-u^0\right|}^2\text{d}x\text{d}y},$ (3)

其中 $\lambda$ 是拉格朗日乘子。第一项保证去噪后图像尽可能的平滑，第二项则尽可能地维持去噪后图像的相似性。用梯度下降法求解 $E\left(u\right)$ 取得极小值所满足的Euler-Lagrange方程，即

$\frac{\partial u}{\partial t}=\text{div}\left(\frac{\nabla u}{\left|\nabla u\right|}\right)-\lambda \left(u-u^0\right).$ (4)

解得的u即为所求的最终修复图像。

最初的TV图像修复模型由去噪模型发展而来，TV模型能在修复的过程中去除噪声的同时又保持图像边缘，且数值PDE实现方便，能较好的实现图像划痕的修复，但原始算法的修复结果易产生明显的阶梯效应。针对此，提出一个结合各向异性和各向同性的非线性扩散修复方程，即

$\{\begin{array}{l}\frac{\partial u}{\partial t}=g\left(\left|G_{\nabla }\right|\right)\left(\Delta u\right)=0,\left|\nabla I_P<0.25\right|,\\ \frac{\partial u}{\partial t}=g\left(\left|G_{\nabla }\right|\right)\left(\text{div}\left(\frac{\nabla u}{\left|\nabla u\right|}\right)-{\lambda }_e\left(u-u^0\right)\right)=0,\left|\nabla I_P\ge 0.25\right|,\end{array}$ (5)

其中 $\left(x,y\right)\in D$ 时 ${\lambda }_e=\lambda$， $\left(x,y\right)\in \Omega \backslash D$ 时 ${\lambda }_e=0$ ； $g\left(\left|G_{\nabla }\right|\right)=1/\sqrt{G_x^2+G_y^2}$ 用于加强迭代的速度； $\sqrt{G_x^2+G_y^2}$

表示图像在点 $\left(x,y\right)$ 梯度向量的模值。式中的 $\nabla I_P$ 表示图像待修复点的梯度值，即其像素值变化的程度。该文献先将图像进行归一化处理，再用Sobel算子计算 $\nabla I_P$，因此可以通过统一的梯度阈值划分使用的修复迭代方程，进而提高选择迭代方程的准确性。接下来，在计算方法过程中，采用半点差分法，对扩散算子 $\text{div}\left(\nabla u/\left|\nabla u\right|\right)$ 以及拉普拉斯算子 $\Delta u$ 进行离散化，使用Sobel算子计算点的梯度向量模值，结合Gauss-Jacobi迭代算法对图像进行扩散迭代，根据给定阙值结合新旧图像的像素值判断迭代状态，停止迭代后，用新像素值代替对应像素点，进而完成图像修复(参见文献 [11] )。

而基于TV模型的图像修补算法，是由Tony Chan等人根据Rudin等人提出的图像去噪模型推广来的。主要是通过建立图像模型，根据整体变分原理，将图像建立的模型转化为约束最优化问题，进而利用Lagrange乘子法求解。本文介绍的算法在其基础上进行了改进，改进后的TV模型有效避免了TV模型算法在平滑区域造成的阶梯效应。同时，通过提高扩散速度，改进了原算法修复过程中只对图像边缘进行扩散造成的效率低的问题，提高了修复效率，减少了计算量。

2.2. 基于偏微分方程的图像修复方法

偏微分方程的主要思想可以追溯到Bertalmio等人最初对于数字图像修复的研究 [12]，他们的模型基于非线性偏微分方程，旨在模仿专门从事修复艺术家的技术。他们认为一个好的修复算法应该将周围区域的尖锐边缘扩散到需要填充的受损部分。并且，利用损伤区域的边缘信息来确定扩散信息和扩散方向。

由于TV算法输出的修复图像不满足视觉连通性，且TV图像修复算法在一些破损面积较大的图像修复中具有修复局限性。因此文献 [13] 在传统TV模型的基础上，结合非线性扩散的思想，提出了一种新的数值算法。并通过数值实验结果表明，这种针对参数的自适应迭代函数图像修复算法在修复过程中，结合扩散系数和自适应迭代函数的共同作用，消除了阶梯效应，完成了修复的图像自然过渡，将图像的边缘信息保存良好。文献 [14] 在其基础上对扩散函数进行修改，得到新的扩散函数如下：

$g\left(\left|\nabla u\right|\right)=\{\begin{array}{l}1,\left|\nabla u\right|=0,\\ 1/{\left(1+\frac{K}{\nabla u}\right)}^2,其他,\end{array}$ (6)

其中，正常数K是调整保护梯度的范围。该扩散函数可以增强小梯度区域的平滑，减少大梯度区域的平滑。使图像边缘清晰，过渡自然。并引入随迭代时间变化的参数 $\beta$ 的自适应函数 $\beta \left(t\right)=1/\left(\epsilon +t/K\right)$。其中 $\epsilon$ 为保证 $\beta \left(t\right)$ 分母不为零的任意正实数。随着迭代次数的增加，与扩散系数函数作用保证了图像的边缘特性。进一步得到修复方程为

$-g\left(\left|\nabla u\right|\right)\cdot \text{div}\left(\frac{\nabla u}{\sqrt{\beta {\left(t\right)}^2+{\left|\nabla u\right|}^2}}\right)+{\lambda }_e\left(u-u^0\right)=0.$ (7)

采用偏导数近似算法对散度算子项进行离散化处理，在目标像素点O周围增加4个点，构造更髙阶的差分格式。结合水平分量 $v_1$ 与垂直分量 $v_2$ 定义单位向量 $v=\left(v_1,v_2\right)=\left(\nabla u/\sqrt{\beta {\left(t\right)}^2+{\left|\nabla u\right|}^2}\right)$。选定目标像素点O为中心点和待修复点， $Q=\left\{N,E,S,W\right\}$ 为四邻域的参考点，以及步长为h，并设步长 $h=1$，并引入 $q=\left\{n,e,s,w\right\}$ 为半像素点以保证稳定性，即可得到

$\text{div}\left(v\right)=\frac{\partial v_1}{\partial x}+\frac{\partial v_2}{\partial y}\cong \left(v_{1e}-v_{1w}\right)+\left(v_{2n}-v_{2s}\right),$ (8)

其中 $U_i\left(i\in Q\right)$ 为各点的像素值； $U_j\left(j\in q\right)$ 为各中间点像素值； $U_o$ 为目标像素点像素值。估算半像素点的梯度值，以e点为例，即 $v_{1e}=\frac{1}{{\left|\nabla U_e\right|}_{\beta }}\left[\frac{\partial U}{\partial x}\right]=\frac{1}{{\left|\nabla U_e\right|}_{\beta }}\frac{U_E-U_O}{h}=\frac{U_E-U_O}{{\left|\nabla U_e\right|}_{\beta }}$，梯度模值为 ${\left|\nabla U_e\right|}_{\beta }=\frac{1}{h}\sqrt{{\left(U_E-U_O\right)}^2+{\left[\left(U_{NE}+U_N-U_S-U_{SE}\right)/4\right]}^2}$。依次算出其余半像素点梯度模值，代入(8)式，得到

$\begin{array}{c}\text{div}\left(v\right)=\left(v_{1e}-v_{1w}\right)+\left(v_{2n}-v_{2s}\right)\\ =\frac{U_E-U_O}{{\left|\nabla U_e\right|}_{\beta }}-\frac{U_O-U_W}{{\left|\nabla U_w\right|}_{\beta }}+\frac{U_N-U_O}{{\left|\nabla U_n\right|}_{\beta }}-\frac{U_O-U_S}{{\left|\nabla U_s\right|}_{\beta }}\\ ={\displaystyle \underset{\begin{array}{l}i\in Q\\ j\in q\end{array}}{\sum }\frac{U_i-U_o}{{\left|\nabla U_j\right|}_{\beta }}}\end{array}$ (9)

将上式代入对应的定常方程，则得到离散化的修复方程为

${\displaystyle \underset{\begin{array}{l}i\in Q\\ j\in q\end{array}}{\sum }\frac{g\left(\left|\nabla U_j\right|\right)}{{\left|\nabla U_j\right|}_{\beta }}}\left(U_o-U_i\right)+{\lambda }_o\left(0\right)\left(U_o-U_o^0\right)=0,$ (10)

其中， $U_o^0$ 表示待修复图像在O点的像素值。采用Gauss-Jacobi迭代法将(10)式简化为

$U_o^n={\left(\frac{{\displaystyle \underset{\begin{array}{l}i\in Q\\ j\in q\end{array}}{\sum }\frac{g\left(\left|\nabla U_j\right|\right)}{{\left|\nabla U_j\right|}_{\beta }}}}{{\displaystyle \underset{\begin{array}{l}i\in Q\\ j\in q\end{array}}{\sum }\frac{g\left(\left|\nabla U_j\right|\right)}{{\left|\nabla U_j\right|}_{\beta }}+{\lambda }_o\left(0\right)}}\right)}^{\left(n-1\right)}{U}_i^{\left(n-1\right)}+{\left(\frac{{\lambda }_o\left(0\right)}{{\displaystyle \underset{\begin{array}{l}i\in Q\\ j\in q\end{array}}{\sum }\frac{g\left(\left|\nabla U_j\right|\right)}{{\left|\nabla U_j\right|}_{\beta }}+{\lambda }_o\left(0\right)}}\right)}^{\left(n-1\right)}{U}_o^0,$ (11)

(11)式即为最终得到离散方程的差分格式，取 ${\lambda }_o\left(0\right)=0$，将各像素值统一带入(10)式并进行迭代，模型的稳定性分析请参考文献 [15]。

该种数值方法是一种对待修复区域实现不同权值的办法，即在TV模型中引入一个以梯度模值作为参数的扩散函数，此类函数在小梯度区域，扩散函数取值较大，增强图像的平滑度。能够有效避免TV模型算法在平滑区域造成的阶梯效应。不仅对小区域破损图像有明显修复效果，对大区域破损同样有效。并且，随着迭代次数不断增加，则与扩散系数函数进一步作用保证图像的边缘特性。同时，通过对峰值信噪比、迭代时间及迭代次数进行比较发现，该算法修复效果更好，提高了修复效率，减少了计算量。

2.3. 基于偏微分方程的CDD图像修复方法

文献 [16] 通过对TV模型的分析后可知，该算法的优点是数值实现简单，对破损图像的边缘有很好的保护作用，特别是对有划痕或者破损区域小的图像有很好的修复效果。但在图像梯度较小的区域易产生阶梯效应，修复效果不好。Chan等人为了解决TV模型的不足，提出了曲率驱动扩散模型。曲率驱动扩散模型在TV模型扩散项的基础上，将曲率项 $g\left(k\right)$ 引入扩散项中 [17]。 $g\left(k\right)$ 是一个关于等照度线的曲率函数，假设图像的待修复区域为D，图像的已知像素区域为E，此时的修复方程为

$\frac{\partial u}{\partial t}=-\nabla \left[\frac{g\left(\left|k\right|\right)}{\left|\nabla u\right|}\nabla u\right]+\lambda \left(u-u_0\right).$ (12)

当 $\left(x,y\right)\in E$ 时， $\lambda =\lambda ,g\left(\left|k\right|\right)=1$ ；当 $\left(x,y\right)\in D$ 时， $\lambda =0,g\left(\left|k\right|\right)=g\left(\left|k\right|\right)$。 $g\left(\left|k\right|\right)$ 是关于k的单调递增函数，在区域E中取值为1，此时模型退化为TV模型。对于函数 $g\left(\left|k\right|\right)$ 中的k定义为

$k=\nabla \cdot \left[\frac{\nabla u}{\left|\nabla u\right|}\right]=\frac{u_y^2{u}_{xx}+u_x^2{u}_{yy}-2u_x{u}_y{u}_{xy}}{{\left(u_{xx}^2+u_{yy}^2\right)}^{\frac{3}{2}}}.$ (13)

在区域D中， $g\left(\left|k\right|\right)$ 通常表示为

$g\left(\left|k\right|\right)={\left|k\right|}^p=k^p,k>0,p\ge 1,$ (14)

该模型中的扩散系数为 $g\left(\left|k\right|\right)/\left|\nabla u\right|$ ；在TV模型扩散项的基础上引入了 $g\left(\left|k\right|\right)$，使得在图像修复的过程中，扩散强度在和图像待修复点的梯度建立关系的同时，与图像像素点的特征信息也有关 [18]。在图像的平坦区域中，图像的信息同时和梯度以及平坦区域中的曲率面有关。根据函数的性质可知引入的 $g\left(\left|k\right|\right)$ 是关于梯度的增函数，在待修复图像的破损区域边缘处梯度较大，因此边缘处的曲率大，进而保证了图像边缘信息能够得到充分的扩散；而在图像破损区域的平坦区域，梯度的模值很小，此时为了保证像素信息的准确扩散，应该在这个方向上获得较弱的扩散。从扩散系数中也可以看出，在梯度变化较大的某些区域，梯度和曲率对扩散系数的影响形成反比。CDD模型的数值实现和BSCB模型一样，求出相应的 $u_x,u_y,u_{xx},u_{yy},u_{xy}$，带入到(12)式中，通常先对(14)式中的P进行取值，再结合(12)式对有划痕或有小破损的图像进行修复。CDD模型在TV模型的基础上，在扩散项中引入了曲率函数，使得图像的扩散和待修复点的梯度以及曲率建立联系，保证在像素信息的特征区域能够准确的扩散。为了方便计算对(14)式中P取值为0.5后，对有划痕的图像进行修复。从文献 [16] 的模型实验结果中可以看出，该模型能够对具有小破损区域的图像进行很好的修复。

该图像修复方法可以有效地对有划痕或者破损区域小的图像进行修复。与其他数值方法相比，该数值模型采用的技术对具有小破损区域的图像可以进行专业的修复，弥补了其他数值模型的不足，提高图像处理水平，一定程度上有助于完善图像修复效果，也是进行数字图像修复问题的一类有效的方法。

2.4. 基于Allen-Cahn方程算子分裂法的图像修复方法

继Li等人 [19] 提出一种快速局部图像修复方法，用算子分裂法求解Allen-Cahn模型后，文献 [20] 提出一种基于Allen-Cahn方程的算子分裂法得图像修复方法，并通过对合成图像和真实图像的数值实验，验证了该方法的准确性、稳定性和有效性。其核心思想是利用算子分裂法将原问题分解为线性方程和非线性方程。分别用有限差分Crank-Nicolson格式和解析法求解线性方程和非线性方程。

考虑具有如下边值条件的Allen-Cahn方程：

$\frac{\partial u\left(x,y,t\right)}{\partial t}=\{\begin{array}{l}\Delta u\left(x,y,t\right)-\frac{1}{{\epsilon }^2}{F}^{\prime }\left(u\left(x,y,t\right)\right),\left(x,y\right)\in {\Omega }_D\\ 0,\left(x,y\right)\notin {\Omega }_D\end{array}$ (15)

$\frac{\partial u\left(x,y,0\right)}{\partial t}=\{\begin{array}{l}0.5,\left(x,y\right)\in {\Omega }_D\\ \frac{g\left(x,y\right)-g_{\min }}{g_{\max }-g_{\min }},\left(x,y\right)\notin {\Omega }_D\end{array}$ (16)

$\nabla u\left(x,y,t\right)\cdot n_{\Omega }=0,\left(x,y\right)\in \partial {\Omega }_D$ (17)

其中 $\epsilon$ 为给定图像， $g_{\max }$ 和 $g_{\min }$ 分别为图像的最大值和最小值； $\partial {\Omega }_D$ 为图像的边界；n表示有界域和域边界上的单位外法向量；参数 $\epsilon$ 为与界面能量相关的梯度能量系数。 $F\left(u\right)=\left(u^2{\left(u-1\right)}^2\right)/4$ 为图像修复中常用函数；选取二维领域 $\Omega =\left[1,N_x\right]\times \left[1,N_y\right]$ 和网格大小为h，对图像区域进行剖分，即有

${\Omega }_h=\left\{\left(x_i,y_j\right)|x_i=\left(i-0.5\right)h,y_j=\left(j-0.5\right)h,1\le i\le N_x,1\le j\le N_y\right\}.$ (18)

令 $u_{ij}^n=u\left(x_i,y_j,n\tau \right)$，其中，时间步长为 $\tau =T/M_t$。Allen-Cahn方程的数值格式如下：

$\{\begin{array}{l}{u}_t=-\frac{F^{\prime }\left(u\right)}{{\epsilon }^2},t\in \left[t^n,t^{n+1/2}\right],\\ u_t=\Delta u,t\in \left[t^n,t^{n+1}\right],\\ u_t=-\frac{F^{\prime }\left(u\right)}{{\epsilon }^2},t\in \left[t^{n+1/2},t^{n+1}\right].\end{array}$ (19)

根据初始条件 $u^n$ 和Crank-Nicolson方法以及五点差分格式可得到第一个子问题 $u^{\ast }$ 以及第二个子问题 $u^{\ast }{}^{\ast }$，最后，用解决第一个子问题的方法来解决第三个子问题，综上即为

$\{\begin{array}{l}{u}^{\ast }=\frac{1}{2}+\left(u^n-0.5\right)/\sqrt{{\text{e}}^{\frac{-\tau }{2{\epsilon }^2}}+{\left(2u^n-1\right)}^2\left(1-{\text{e}}^{\frac{-\tau }{2{\epsilon }^2}}\right)},\\ \frac{u_{ij}^{\ast \ast }-u_{ij}^{\ast }}{\tau }=\frac{u_{i+1,j}^{\ast \ast }+u_{i-1,j}^{\ast \ast }-4u_{i,j}^{\ast \ast }+u_{i,j+1}^{\ast \ast }+u_{i,j-1}^{\ast \ast }}{2h^2}+\frac{u_{i+1,j}^{\ast }+u_{i-1,j}^{\ast }-4u_{i,j}^{\ast }+u_{i,j+1}^{\ast }+u_{i,j-1}^{\ast }}{2h^2},\\ u^{n+1}=\frac{1}{2}+\left(u^{\ast \ast }-0.5\right)/\sqrt{{\text{e}}^{\frac{-\tau }{2{\epsilon }^2}}+{\left(2u^{\ast \ast }-1\right)}^2\left(1-{\text{e}}^{\frac{-\tau }{2{\epsilon }^2}}\right)}.\end{array}$ (20)

该算法首先确定图像需要修复的区域为 ${\Omega }_D$。同时，令最接近修复区域的边界 ${\Omega }_D$ 作为修复区域计算时的边界条件。接下来，引入控制函数f，当 $\left(x_i,y_j\right)\in {\Omega }_D$ 时，将其定义为 $f_{ij}=1$，对于其他情况，则定义 $f_{ij}=0$。最后通过(19)式当且仅当 $f_{ij}=1$ 时更新 $u^{n+1}$。

采用基于Allen-Cahn方程算子分裂法用于图像修复，不仅保证修复过程的有效性、准确性以及稳定性，而且相对于含有四阶微分算子的Cahn-Hilliard方程计算得更加精简，其时间和空间精度都可以达到二阶。并且，在计算过程中选择五点差分法，减少了计算量。由于该方法只应用于修复域，其余区域的像素值与原始输入图像的像素值保持一致，因此大大提高了计算效率。

3. 总结

图像修复技术近年来发展十分迅速，并被广泛应用于计算机科学、天文学、考古、人工智能等领域。本文针对图像修复中常见的问题介绍了几类效果较好且高效的图像修复数值模型，限于篇幅，关于模型稳定性以及数值实验本文没有展开讨论，有兴趣的读者可以阅读文后的参考文献。现阶段现存的几类数值模型已经允许处理较大的数据集，且计算速度较快，能够提高图像处理效率，是进行数字图像去噪与修复问题有效的方法。