基于灰色系统理论的图像修复研究

曾卫波，邢永康，石 杨

(重庆大学计算机学院，重庆 400030)

1 概述

图像修复指根据图像已有信息来填充自然或人为缺损信息。目前图像修复技术一般分为三大类:基于局部信息扩散的结构修补，基于采样复制的纹理修补和基于图像分解的方法。基于局部信息扩散的传统结构修补法又可分为2 类:(1)基于偏微分方程(Partial Differential Equation，PDE)算法，典型代表有BSCB 方法［1］，曲率推动扩散方法(CDD)［2］和基于变分的TV 模型［3］等，这些方法通过全局反复迭代，求解偏微分方程，目标是将周围已知信息按等照度方向扩散至破损区域。这类方法适于小区域修补，由于其反复迭代求解偏微分方程，修补大块状区域速度过慢，当边缘区域为奇异值点时无法求解，不能保证产生强边缘;(2)邻域加权插值算法。这类算法结合了滤波算法思想，认为待修复值与邻近像素值相关。典型代表有TELA 的快速行进法(FMM)［4］采用邻域加权和，从外层向内层插值。这类算法与PDE 算法的反复迭代求解不同，它只需要作一次插值处理即可，因而速度较快。但是这类方法使用加权和并不能延续灰度值的变化趋势，并且需要严格控制修复顺序，否则修复区域过大时，拓扑结构延伸效果并不理想。基于采样复制的纹理修补在图中寻找与待修复点周围的纹理块最匹配的部分，并填充到待修复区［5］。它对大区域块修复效果较好，由于采用全局搜索匹配块，可能会产生错误匹配和填充速度过慢等问题。基于图像分解的方法将图像分解为结构图和纹理图，对纹理图采用纹理合成法，结构图采用局部信息扩散的结构修补，最后合成两者得到结果［6］。上述方法在修复时间上难以满足要求，在修复质量上也有所欠缺。

近年来，一些学者致力于推广灰色系统理论的应用，灰色系统理论［7］由我国控制论专家邓聚龙教授于20 世纪80 年代创立，在研究少数据、贫信息和不确定性问题上，有着较为不错的效果，因而受到普遍关注。国内对灰色理论在图像处理上的应用已有一定研究，如图像压缩、噪点去除和边界检测等［8］。而图像修复主要依靠灰色预测完成，灰色指信息不完全，介于已知与未知之间，灰色预测指通过部分已知信息来推导未知信息。文献［9］将灰色预测GM(1，1)模型引入图像修复中，通过一系列已知像素值建立模型来预测未知像素值，以此延续灰度值的走势。文献［10］将灰色预测用在UDP 传输丢包图像修复上，对这种丢失区域少的特定问题产生了快速且较好的修复效果。文献［9 -10］均用单条序列来预测像素值，并不能充分挖掘周围各方向传递过来的信息，其预测值的可信性并不高。当边缘点或者噪点出现时，预测所用的序列容易出现陡变情况，模型不能做很好的预测和反馈，此时等照度线并不能有效地保持。灰色预测法修复图像只需作一次插值，速度快，但其质量改善还有待进一步研究。

综上所述，针对图像中的结构区和由分解而得的结构图，难以寻找一种修复质量和修复时间兼顾的方法。本文在前人研究基础上，使用灰色理论，研究改善其修复图像的质量，并进行如下工作:在预测模型上选取SCGMmv(1，1)模型［11］，该模型不仅在精度和速度上要优于GM(1，1)模型，而且利用该模型的一些特性能很好地处理边界或者噪声出现在序列中的陡变情况。引用灰色关联度［12］区分噪点和正常点，并分别处理。在对每个未知像素点进行插值时，从它的多个邻域方向进行预测后并决策最终值，能更充分利用已知信息。在填充过程中省去了传统的优先度的计算问题，而是根据已知邻域序列条数和是否边缘点来控制填充顺序，该方法不需要复杂的计算，且起到了类似优先度的功能。这些改进使得修复结果不仅快速且质量较好。

2 图像局部性和SCGMmv(1，1)预测模型

2.1 图像的局部特性

数字图像由一系列离散像素点组成，由马尔科夫相关性可知，在图像局部区域内像素点之间可以近似归纳出一些特性，文献［13］对图像局部性进行了一定的阐述:图1 中的最后2 幅图像在局部区域不存在边缘时，像素值变化可分为平滑或者平坦变化。

图1 局部图像特征

图1 中前3 幅图像局部区域存在边缘时，边缘可以近似为一条像素值接近的直线，边缘两侧分别为平滑或者平坦区域。本文选用了SCGMmv(1，1)模型来挖掘这些关系，SCGMmv(1，1)模型是以系统云［14］为背景，按基于均值生成时序［11］和灰色趋势关联分析的灰色动态建模原理构造而成的时间序列预测模型。

2.2 均值生成时序

已知序列:

按积分生成原理，有α 权生成时序:

用矢量和矩阵形式表示:

令α=0.5，得均值生成时序:

均值生成运算能降低时序噪声影响，使观测值运算效果等同于“真值”效果。此时的均值生成时序可用于构建系统模型。

2.3 预测模型

设已知函数集为:

通过对x(0)和{f1，f2，…，fm}之间的趋势关联分析，得趋势关联度:

设θxfr=max{θxfj}，j=1，2，…，m，则fr是可以定义在时域上的一个确定函数，例如非齐次指数函数:

其中，a，b，c∈R，k=1，2，…，n，对θxfr满意就意味着认为fr隐含于x(0)，在此就可以直接用的数据拟合于求得估值和即依照系统云建模原理［14］，有:

当k=1，2，…，n 时，有估值:

其解(预测模型)为:

3 SCGMmv(1，1)模型在图像修复上的应用

如图2 所示D 为已知区域，Ω 是待修复未知区域，P 点是位于交界上的待修复点。选取P 点各条邻域方向上，5 个连续的已知像素值作为初始时序序列(图2 中为3 条满足5 个连续已知像素值的邻域序列，从左到右依次为X1、X2、X3，即:

图2 图像修复示意图

依照上节的建模步骤，式(11)中取k=6 时，分别得出3 个P 点的灰度估计值u1(i，j)u2(i，j)和u3(i，j)，最后通过某种决策机制来决定P 点最终的像素值up(i，j)。

3.1 单序列预测

文献［9 -10］均用GM(1，1)模型作为预测模型，这里对GM(1，1)模型不做赘述，而本文引入SCGMmv(1，1)预测模型。首先在挖掘图像数据变化特性上，本文的模型更为敏感:在平坦区域和平滑区域，本文模型比GM(1，1)模型更为可靠，而在陡变情况如噪点或者边缘阶跃变化出现时，本文模型中代表序列发展系数的ea计算结果相对较大，会预测出明显异常值，这个信息十分可贵，表明此时需对噪点或者边缘点出现情况时作特殊的预测处理，而GM(1，1)模型则无法挖掘这一点。其次在建模速度上，本文模型不需要像GM(1，1)模型一样做累加生成、累减还原等环节，计算量小，适于快速动态建模。

表1 针对几种典型的灰度值变化序列给出了2 种模型的预测案例。第1 行为平坦区域，2 种模型结果都正常(灰度值波动范围在15 以内)。第2 行～第5 行为从平坦区跳变到边缘区的情况，可以看到本文模型要么出现明显异常值，要么出现正常预测为边缘点。而GM(1，1)模型预测也正常(只针对序列变化)，但是无法判断出是否边缘点。第6 行以增量20 匀速递增，第7 行为加速递增，第8 行、第9 行为从平滑区过渡到边缘，本文模型预测精度要高得多。第10 行、第11 行显示噪点出现在序列中不同位置时的情况。从表中可以看到这些异常出现时受噪点或者边缘点的影响，序列出现陡变情况。下节将具体讨论这些异常处理。这里只列出了几种典型的递增或者随机平稳变化的例子，而递减变化预测也有类似的特性。

表1 GM(1，1)和SCGMmv(1，1)模型预测

3.2 异常检测和处理

假设序列x1，x2，x3，x4，x5，预测值为x6。选取一定阈值T，如果视为出现异常。T 为0 时则对每个预测值都要作异常检测处理，T 很大时虽然检测次数少，但可能会漏掉异常情况，实验中可以调节不同阈值T，本文实验选取阈值5，可以满足一般情况。异常出现时序列上的第5 个点可能出现2 种情况，噪点或者非噪点，如果是非噪点点如表1 中第3 行的边缘点，那么预测值将与之相同，同预测为边缘点。如果是噪点，则将预测值置为噪点周围正常点的平均值。问题关键便在如何判断已知序列上最后一个点是噪点还是正常点，特别是噪点和边缘点在某种程度上十分类似，难以区分。文献［15-16］介绍了基于灰色关联度的噪点检测法:其主要思想是噪点和正常点的灰色关联度不同，噪点的灰色关联度较小，正常点的灰色关联度较大，设定阈值可识别出正常点和噪点，并做不同处理。算法过程如下:

(1)P 点是待判断点，找出P 点周围与P 点距离小于2 的已知点，可以不论顺序，假设已知点灰度值序列为{x(0)，x(1)，…，x(n)}。

(2)求取序列灰度平均值arv。

(3)求差序列:

(4)求差序列中最大值max 和最小值min。

(5)求取各点的灰色关联度:

其中，包含了P 点的关联度ξ(p)，选取阈值T 为0.88，若ξ(p)＞T，则P 点可视为正常点，否则视为噪点，并求取上述序列中关联度大于T 的像素点的平均值。

在实用中还需要注意，模型中若序列值间隔相等时，式(5)出现分母为0 时，需要对某个像素值做加1 处理，若序列匀速增长(减小)，式(6)出现分母为0，则直接对最后一位值做加(减)变化量便可。出现正常预测情况下，预测值大于255 或者小于0，则预测值置为255 或0。

3.3 多序列决策

上面所述针对单条序列处理情况，而图2 中P点的最终值由多条序列预测值决定。此时需要做决策处理，决策约束如下:

(1)为了增加可信性，至少有4 个方向的序列参与预测，即至少有4 个预测值参与决策。

(2)选取预测值中的最大值和最小值之差，选取一定阈值，一般设为15 左右时，感官上会有明显灰度差，并且可以容许预测误差波动。若两者之差小于阈值，则视为平滑点，最终预测值取其均值，否则视为边界点。

(3)对于边界点处理，要结合下节的修复顺序，在修复过程中记录最近修复过的平滑区域的灰度值v 。然后对几个预测值分别以最大值和最小值为中心聚2 类，聚类直径为15 左右。每一类至少包含2 条序列。然后选取这2 类中与v 值差距大的一类平均值 作为边界点预测值，这样便完成了2 个区域的跳跃过程。

4 区域修复顺序

把图像修复过程看成蚕吃桑叶过程，叶子上的茎为边缘区，叶子区域即为平滑区域，蚕在吃桑叶过程中，先找叶子凸出的地方下口，然后慢慢蚕食，遇到茎或者凹进时停止吞噬，寻找另外的凸出，直到最后吃完整片叶子，蚕食过程中能完整保留叶子的茎，这正好对应图像中的边缘信息。

本文修复顺序模拟上述过程，凸出的地方代表已知邻域序列多，凹进的地方已知邻域序列少。创建3个队列:q1，q2 和q3，把已知邻域序列条数大于等于5的点，放入q1 中，已知邻域序列条数为4 的点放入队列q2 中，边缘点放入边缘队列q3 中，算法过程如下:

(1)在待修复区域边界上，搜索已知邻域序列条数最多的几个点，作为起始点按要求分别放入q1和q2 。

(2)如果q1，q2 和q3 都空，则修复完毕结束;否则转入步骤(3)。

(3)如果q1 不空，如下处理q1 中元素直至为空，按照上节方法判断q1 队列头像素点是否为边界点，是则插入边缘队列q3 滞后处理，否则修复该点，并搜索该点邻域未知点，按序列条数把邻域未知点点插入q1 和q2;否则转到步骤(4)。

(4)如果q2 不空，按照上节方法判断q2 队列头像素点是否为边界点，是则插入边缘队列q3 滞后处理，否则修复该点，并搜索该点邻域未知点，并按序列条数把邻域未知点插入q1 和q2，转到步骤(3)。如果q2 为空转到步骤(5)。

(5)如果q3 不空，如下处理q3 元素直至为空，按照上节修复边界点方法修复队列头像素点，并搜索该点邻域未知点，并按序列条数把邻域未知点插入q1 和q2 。如果q3 为空转到步骤(2)。

从图3 和图4 中可以看到处理过程，上述顺序控制可以保证，周围已知信息多的点优先修复，平滑区域优先处理。由于每个点至少要有4 条已知邻域序列，并且边缘滞后处理，实质上这种顺序控制使得修复是从八邻域方向的直线来连接边缘的。

蚕食顺序控制法使得图像的拓扑结构得到了比较好的延伸，且不需要复杂的优先度计算公式，其思想直观易懂。修复算法内存空间使用有如图5 所示的变化特性。

其空间占用主要由上述队列大小决定，为了方便分析，只有一个待处理队列作为修复缓存，存储就绪的待修复像素点，初始时队列缓存有a 个像素点，队列变化体现在修复一个像素点即出队列，搜索修复点邻域满足待修复点入队列。前期修复一个像素点平均要插2 个像素点入队列，中期每修复一个点，插入一个点，到了后期修复一个点，插入0 个点。队列最大值取决于某时刻修复区域内连接的待修复边界最大周长。当修复区域越大时，可能占用的内存也越多。在本次实验中，设置了队列容量400 个像素点，便可满足修复要求。

图3 lena 图鼻子修复过程

图4 lena 图镜框修复过程

图5 修复队列大小变化示意图

5 实验结果与分析

本文算法在英特尔双核2.5 GHz CPU，4 GB 内存微机中，在vs2010 平台下结合opencv 编程实现，分别对模拟图像和真实图像进行修复实现，并与基于曲率推动方法的CDD 算法和快进算法FMM 对比。

5.1 模拟图像修复实验

从图6、图7 中可以看到CDD 修复对边缘插值较为模糊，无法产生强边缘。在图7 圆盘修复实验中，可以看见上端的水平线一直往下延伸，说明CDD的曲率推动收敛性是一个问题。而对图8、图9，当修补区域较大时，需要迭代很多次，才初见效果，每次迭代计算量变大，耗时非常长，在实时修复应用中可以说是失败的。而FMM 算法虽然速度很快，基本上在1 s 内完成，但是修复质量却不尽人意。本文算法也在1 s 内完成，并且边缘保持的视觉效果明显优于其他算法。

图6 破损三角形修复

图7 破损圆盘修复

图8 椭圆修复

图9 T 形修复

5.2 真实图像修复实验

本文获取2 张结构性很强的真实图像，对图10地板砖破损图像，CDD 迭代5 000 次和10 000 次并无明显区别，对图11 的lena 破损图像CDD 迭代20 000次和40 000 次并无明显区别。

图10 星形图像修复

图10 实验中，CDD 修复失败，FMM 修复效果较好，是由于图像的对称性正符合了该算法的快进顺序。从图11 实验中可以看出，虽然在lena 背部平滑区CDD 修复效果较其他2 种好，但是对大块状区域和边缘连接上(背部头发断裂)并不理想，并且修复时间在小时级。FMM 算法和本文算法对这2 幅图像均在2 s 内完成修复，但是从图12 的修复细节来看，它和CDD 一样，对于复杂区域修复效果都不如本文算法好。由于图像修复是一个病态性问题，本身破损区域的真实像素值无从得知，因此本文没有做PSNR 对比，而图像修复更多是从主观视觉满足性来衡量，本文算法无论是时间上还是视觉质量上都取得了较优的效果。

图11 lena 图像修复

图12 lena 图像修复细节

6 结束语

本文提出一种简单快速的结构图像修复算法，使用预测模型从多个方向预测未知像素值，并引用灰色关联理论处理噪点。由于每个修复点只需在预测决策后作一次插值处理，因此速度较快。在修复顺序上模拟蚕食过程也保证了平滑区和边缘区的修复质量。但是对于纹理图像，虽然也有预测模型能对规则纹理做出较为准确的预测，但绝大部分纹理图像变化太过剧烈且无规律，所以预测模型对纹理图像修复作用不大。另外，当噪声密度过大时，其情形类似于纹理，预测模型和边界点检测也会出现失效，针对这种情况在修复前需进行去噪平滑，如使用TV 模型分解图像时便可得到具有较强局部特性的结构图。下一步研究方向是对图像进行分解后，用本文方法对结构图像快速修复，然后用纹理修补方法修复纹理图像，最后进行合成，完成图像修复。

［1］Bertalmio M，Sapiro G，Caselles V，et al.Image Inpainting［C］//Proc.of ACM SIGGRAPH Conference on Computer Graphics.New York，USA:ACM Press，2000:417-424.

［2］Chan T F，Shen J.Non-texture Inpainting by Curvaturedriven Diffusions (CDD)［J］.Journal of Visual Communication Image Representation，2001，12 (4):436-449.

［3］Chan T F，Shen J.Mathematical Models of Local Nontexture Inpaintings ［J］.SIAM Journal of Applied Mathematics，2002，62(3):1019-1043.

［4］Telen A.An Image Inpainting Technique Based on the Fast Marching Method［J］.Journal of Graphics Tools，2004，9(1):25-36.

［5］Criminisi A，Perez P，Toyama K.Object Removal by Exemplar-based Inpainting ［C］//Proc.of IEEE Computer Scoiety Conference on Computer Vision and Pattern Recognition.［S.1.］:IEEE Press，2003:18-20.

［6］Bertalmio M，Vese L，Sapiro G，et al.Simultaneous Structure and Texture Image Inpainting［J］.IEEE Transactions on Computer Vision，2003，12 (8):882-889.

［7］Deng Julong.Control Problems of Grey Systems［J］.Systems ＆ Control Letters，1982，1(5):288-294.

［8］马 苗，田红鹏，张艳宁.灰色理论在图像工程中的应用研究进展［J］.中国图象图形学报，2007，12(11):1943-1951.

［9］吴长勤，段汉根.基于GM(1，1)预测的残缺图像的修复算法［J］.计算机系统应用，2011，20(4):173-176.

［10］吴长勤，段汉根.基于灰色预测的残缺图像的修复算法［J］.计算机技术与发展，2010，20(5):124-127.

［11］陈为真，陈智洁，谢兆鸿，等.均值生成时序及SCGMmv(1，1)预测模型［J］.中南民族大学学报，2003，4(3):32-39.

［12］刘法贵，张愿章，李湘露.灰色数学及其应用［M］.开封:河南大学出版社，2002.

［13］屈 磊，韦 穗，梁 栋，等.快速自适应模板图像修复算法［J］.中国图象图形学报，2008，13(1):24-28.

［14］Chen Mianyun.Principle of Grey Dynamic Modeling［J］.SAMS，1996，26(3):69-79.

［15］黄春艳，张云鹏，黄红艳.基于灰色关联度的图像混合噪声的自适应滤波算法［J］.微电子学与计算机，2010，27(2):126-128.

［16］李俊峰，戴文战，潘海鹏，等.基于灰色系统理论的图像去噪算法研究［C］//第29 届中国控制会议论文集.北京:［出版者不详］，2010:2671-2676.