沈阳师范大学，数学与系统科学学院，辽宁 沈阳

收稿日期：2021年2月8日；录用日期：2021年3月8日；发布日期：2021年3月16日

摘要

图像修复作为图像处理的一个分支，在计算机视觉、天文学、生物学等领域中有广泛的应用。本文使用修正Cahn-Hilliard方程进行二值图像修复，采用二阶有限差分方法将含有非线性项的方程在空间上进行离散，采用Crank-Nicolson方法将其在时间上进行离散，应用快速离散余弦变换结合不动点迭代法求解全离散格式下的方程组。基于该模型的图像修复数值方法具有参数少、存储量小、计算效率高等优点。最后，给出数值实验，数值结果验证该数值方法能有效地进行图像修复与去噪。

关键词

图像修复，Cahn-Hilliard方程，二值图像，有限差分法

Binary Image Inpainting Method Based on Cahn-Hilliard Equation

Junrong Huo, Rongpei Zhang, Hao Liu

School of Mathematics and Systematic Sciences, Shenyang Normal University, Shenyang Liaoning

Received: Feb. 8^th, 2021; accepted: Mar. 8^th, 2021; published: Mar. 16^th, 2021

ABSTRACT

As a branch of image processing, image inpainting is widely used in computer vision, astronomy, biology and other fields. In this paper, the modified Cahn-Hilliard equation is used for binary image inpainting. The second-order finite difference method is used to discretize the equation with nonlinear term in space, and the Crank-Nicolson method is used to discretize it in time. The fast discrete cosine transform combined with fixed point iteration method is used to solve the equations in the fully discrete scheme. The numerical method of image inpainting based on this model has the advantages of few parameters, small storage and high computational efficiency. Finally, numerical experiments are given to verify the effectiveness of the proposed method.

Keywords:Image Inpainting, Cahn-Hilliard Equation, Binary Image, Finite Difference Method

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

图像处理的一个重要任务是根据从图像受损区域周围收集到的信息来填充图像的缺失部分，称为图像修复。图像在传输和保存等过程中，可能会有部分信息丢失或者被污染，图像修复一般包括缺损信息的填补和去除噪声，通过未受破坏区域的像素信息，还原被破坏像素以及丢失像素。图像修复是图像退化过程的逆过程，可以恢复由于划伤、老化或自然因素等原因而发生部分破损的图像，也可以恢复被其他对象遮挡的图像。因此，被广泛用于修复古老字画 [1]，视频后期处理 [2]，减少医学成像伪影 [3] 以及卫星图像处理 [4] 等众多领域中。Bertalmio等人 [5] 首次提出数字图像修复的概念，并将偏微分方程的思想引入图像修复中。

近年来，许多学者对图像修复及其数值模型进行了研究。包括TV修复模型 [6] [7]，MS修复模型 [8]，快速图像修复算法 [9]，自适应图像修复算法 [10] 以及应用于二值图像修复的修正Cahn-Hilliard [11] 模型，由于Cahn-Hilliard模型的非凸性，一般采用凸分裂 [12] 的思想进行求解，但是需要定义多个参数。本文应用有限差分方法空间求解Cahn-Hilliard方程，利用Kronecker积写出二维拉普拉斯算子的微分矩阵并将其特征分解，结合快速离散余弦变换实现了快速求解。与原有的方法比较，本文所采用的方法选取的参数更少、存储量小、可以提高计算效率。

2. 数值方法

在二值图像情况下，修正的Cahn-Hilliard方程可更加快速高效地进行图像修复 [13]。本文考虑Cahn-Hilliard方程的修正方程(1)：

$\frac{\partial u}{\partial t}=\Delta \left(-\epsilon \Delta u+{\epsilon }^{-1}{W}^{\prime }\left(u\right)\right)+\lambda \left(f-u\right),\left(x,y,t\right)\in \Omega \times \left(0,T\right].$ (1)

其中 $\Omega =\left[a,b\right]\times \left[c,d\right]$ 为可操作区域； ${\epsilon }^{-1}{W}^{\prime }\left(u\right)$ 为自由能函数； $\Delta =\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$ 为拉普拉斯算子； $W\left(u\right)={\left(u^2-1\right)}^2/4$ 为Lyapunov泛函；当 $x,y\in D$ 时 $\lambda =0$，当 $x,y\in \Omega \backslash D$ 时 $\lambda ={\lambda }_0$ ； $f\left(x\right)$ 为给定区域 $\Omega$ 内的待处理图像， $D\subset \Omega$ 为信息缺失区域， $\Omega \backslash D$ 为有效信息区域； $\epsilon$ 为待定参数。图像修复，即为根据待修复图像在有效信息区域中 $u_0$ 的数值，对信息缺失区域中的像素值进行计算与匹配，进而得到完整图像的过程。设初始值为 $u_0$ 并且边界条件为齐次Neumann边界：

$\frac{\partial u}{\partial n}=\frac{\partial }{\partial n}\left(-\epsilon \Delta u+{\epsilon }^{-1}{W}^{\prime }\left(u\right)\right)=0,\left(x,y,t\right)\in \partial \Omega \times \left(0,T\right],$ (2)

本文考虑用直线 $x=x_i,y=y_j$ 在区域 $\Omega$ 上打网格，其中网格节点为 $\left(x_i,y_j\right)=\left(a+\left(i-1/2\right)\Delta x,c+\left(j-1/2\right)\Delta y\right),i=1,2,\cdots ,N_x;j=1,2,\cdots ,N_y$，网格步长为 $\Delta x=\left(a-b\right)/N_x$，$\Delta y=\left(c-d\right)/N_y$。定义离散解 $u_{i,j}$ 为 $N_x\times N_y$ 阶的矩阵 $U$，表示网格节点 $\left(x_i,y_j\right)$ 处的数值解。由齐次Neumann边界条件，以及网格内部点 $\left(x_i,y_j\right)$ 和边界处的二阶导数差分格式，得到该方程在 $x,y$ 方向的微分矩阵分别为 $B_{N_x}$ 和 $B_{N_y}$。其中矩阵 $B_{N_x}$ 为 $N_x\times N_x$ 阶的矩阵， ${\left(B_{N_x}\right)}_{i,i}=2,i=2,\cdots ,N_x-1$，${\left(B_{N_x}\right)}_{i,i-1}=-1,i=2,\cdots ,N_x$，${\left(B_{N_x}\right)}_{i,i+1}=-1,i=1,\cdots ,N_x-1$，${\left(B_{N_x}\right)}_{1,1}={\left(B_{N_x}\right)}_{N_x,N_x}=1$，矩阵 $B_{N_x}$ 中其他位置的元素均为0。同理有 $B_{N_y}$ 为 $N_y\times N_y$ 阶的矩阵。

将解矩阵 $U$ 按列向量化后得到 $\mathcal{U}$，记为 $\mathcal{U}=\text{vec}\left(U\right)$。设 $I_{N_x},I_{N_y}$ 分别为 $N_x,N_y$ 阶单位矩阵，利用微分矩阵 $B_{N_x}$ 和 $B_{N_y}$ 以及Kronecker积的定义，可将(1)式的中心差分格式写成如下形式：

$\frac{\text{d}\mathcal{U}}{\text{d}t}=-\epsilon K^2\mathcal{U}+{\epsilon }^{-1}K{\mathcal{W}}^{\prime }\left(\mathcal{U}\right)+\lambda \left(\mathcal{F}-\mathcal{U}\right),$ (3)

其中

$K=-\left(I_{N_y}\otimes \Delta x^{-2}{B}_{N_x}+\Delta y^{-2}{B}_{N_y}\otimes I_{N_x}\right),$

$K^2=\left(I_{N_y}\otimes \Delta x^{-4}{B}_{N_x}^2+2\Delta x^{-2}\Delta y^{-2}{B}_{N_y}\otimes B_{N_x}+\Delta y^{-4}{B}_{N_y}^2\otimes I_{N_x}\right).$

下面将矩阵 $B_{N_x}$ 特征分解得到

$B_{N_x}=C_{N_x}^{\text{T}}{\Lambda }_{N_x}{C}_{N_x}=C_{N_x}^{-\text{1}}{\Lambda }_{N_x}{C}_{N_x},$ (4)

其中 ${\Lambda }_{N_x}=\text{diag}\left({\lambda }_0^x,{\lambda }_1^x,\cdots ,{\lambda }_{{}_{N_x-1}}^x\right)$ 为对角矩阵，其对角元素 ${\lambda }_i^x=2-2\cos \left(i\pi /N_x\right),i=0,1,\cdots ,N_x-1$ 为矩阵 $B_{N_x}$ 的特征值，矩阵 $C_{N_x}$ 中 $C_{N_{x_{1,j}}}=\sqrt{1/N_x},j=1,\cdots ,N_x$，$C_{N_{x_{i,j}}}=\sqrt{2/N_x}\cos \left[\left(2i-1\right)\left(j-1\right)\pi /2N_x\right],$ $i=2,\cdots ,N_x;j=1,\cdots ,N_x$，同理，对 $B_{N_y}$ 有类似的结果。接下来，对矩阵 $K$ 做特征分解。由 $I_{N_x}=C_{N_x}^{\text{T}}{I}_{N_x}{C}_{N_x}$ 以及Kronecker积的性质 [13] 可以得到

$K=-{\left(C_{N_y}\otimes C_{N_x}\right)}^{\text{T}}\Lambda \left(C_{N_y}\otimes C_{N_x}\right),$ (5)

其中 $\Lambda =I_{N_y}\otimes \Delta x^{-2}{\Lambda }_{N_x}+\Delta y^{-2}{\Lambda }_{N_y}\otimes I_{N_x}$ 为 $N_x{N}_y\times N_x{N}_y$ 阶的对角矩阵，且对角元素均为正数。

本文应用Grank-Nicolson的差分方法对(3)式进行时间离散得到

(6)

其中 ${\mathcal{F}}_{n+1/2}=\left({\mathcal{F}}_{n+1}+{\mathcal{F}}_n\right)/2$，$K^2={\left(C_{N_y}\otimes C_{N_x}\right)}^{\text{T}}{\Lambda }^2\left(C_{N_y}\otimes C_{N_x}\right)$。上式中运算符“ ${\cdot }^2$ ”表示对向量中每个元素进行平方运算；运算符“ $\cdot *$ ”表示对两个向量中对应元素进行乘法运算。应用快速离散余弦变换结合不动点迭代法求解(6)式，迭代初值选取为 ${\mathcal{U}}_{n+1}^0={\mathcal{U}}_n$。

3. 数值实验

首先，对截断的条形图进行修复处理，如图1所示。其中，图1(b)中灰色区域表示修复区域，即信息缺失区域。采用修正的Cahn-Hilliard方程对截断的条形图分别进行修复10、20、35、55次，可以看出截断的条形图逐渐修复，到迭代55次时图像已经基本修补完成。

图1. 截断条形图像的修复

接下来，对带有涂鸦的图像进行修复处理，如图2所示。其中，图2(b)为带有涂鸦的熊猫图片，表示待修复图像，该图以涂鸦部分为修复区域。采用修正的Cahn-Hilliard方程对图片的缺失部分进行修复5、10、15、20次，可以看出涂鸦区域随着迭代次数的增加，变得越来越细小，直至消失。图2(f)为迭代20次后的修复图像，可以看出图像中破损的部分已经得到有效修补，该图像显示的即为修复后完整的熊猫图片。

图2. 图像中涂鸦的修复

最后，对带有水印的图像进行修复处理，如图3所示。其中，图3(b)为带有水印的向日葵图片，表示待修复图像，该图以水印部分为修复区域。采用修正的Cahn-Hilliard方程对图像分别进行修复5、10、20、30次，可以看出迭代30次时图像中的水印已消失，图像修复完成。

图3. 图像中水印的修复

4. 结论

本文利用四阶修正Cahn-Hilliard方程对截断的条形图以及带噪声的图像进行修复，在空间和时间上分别应用有限差分方法和Crank-Nicolson方法对修正的Cahn-Hilliard方程进行离散并进行快速求解。数值实验结果表明，该模型对于截断的条纹或有涂鸦、水印的二值图像均有良好的修复效果。本文提出的数值方法可以高效地对图像进行修复与去噪，使图像中受污染区域的像素信息得到合理填充，更加适合人眼视觉。同时，该方程可用的快速数值技术在处理较大的数据集上也更加有效，一定程度上有助于加快计算速度，提高图像处理效率。

基金项目

辽宁省自然科学基金(20180550996)资助的课题。

文章引用

霍俊蓉,张荣培,刘 昊. 基于Cahn-Hilliard方程的二值图像修复方法
Binary Image Inpainting Method Based on Cahn-Hilliard Equation[J]. 应用数学进展, 2021, 10(03): 674-679. https://doi.org/10.12677/AAM.2021.103073

参考文献