哈尔滨工业大学学报  2021, Vol. 53 Issue (8): 72-80  DOI: 10.11918/202101038
0

引用本文 

陈永, 陶美风, 陈锦. 块核范数的RPCA分解与熵权类稀疏的壁画修复[J]. 哈尔滨工业大学学报, 2021, 53(8): 72-80. DOI: 10.11918/202101038.
CHEN Yong, TAO Meifeng, CHEN Jin. Mural inpainting based on RPCA decomposition of block nuclear norm and entropy weighted clustering sparse representation[J]. Journal of Harbin Institute of Technology, 2021, 53(8): 72-80. DOI: 10.11918/202101038.

基金项目

教育部人文社会科学研究青年基金(19YJC760012);兰州交通大学天佑创新团队(TY202003);甘肃省人文社会科学一般项目(20ZC11)

作者简介

陈永(1979—),男,教授,硕士生导师

通信作者

陈永,edukeylab@126.com

文章历史

收稿日期: 2021-01-11
块核范数的RPCA分解与熵权类稀疏的壁画修复
陈永1,2, 陶美风1, 陈锦1    
1. 兰州交通大学 电子与信息工程学院,兰州 730070;
2. 甘肃省人工智能与图形图像处理工程研究中心,兰州 730070
摘要: 针对图像修复过程中,颜色纹理光学属性分离不彻底,以及在稀疏表示图像修复时字典设计单一,导致壁画图像修复结果易出现结构不连贯和模糊效应等问题,提出了一种基于块核范数的鲁棒主成分分析(robust principal component analysis,RPCA)分解与熵权类稀疏的壁画修复方法。首先,采用提出的基于块核范数的RPCA图像分解算法,将壁画图像分解为结构层和纹理层,利用块核范数进行纹理矫正操作,克服了RPCA结构纹理分离不完全的问题。然后,提出熵加权k-means方法对结构层图像进行聚类,构建得到稀疏子类字典,并通过奇异值分解和分裂Bregman迭代优化的类稀疏修复方法,完成结构层图像的重构。最后,利用双三次插值算法实现对纹理层图像的修复,将修复后的结构层和纹理层进行融合,完成破损壁画的修复。通过对真实敦煌壁画数字化修复, 实验结果表明,该算法能够有效地保护壁画图像的边缘和纹理等重要特征信息,无论从视觉效果还是从峰值信噪比等定量评价方面,提出的方法修复效果均优于比较算法,且修复执行效率更高。
关键词: 图像重构    壁画修复    RPCA分解    块核范数    类稀疏表示    
Mural inpainting based on RPCA decomposition of block nuclear norm and entropy weighted clustering sparse representation
CHEN Yong1,2, TAO Meifeng1, CHEN Jin1    
1. School of Electronic and Information Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China;
2. Gansu Provincial Engineering Research Center for Artificial Intelligence and Graphics & Image Processing, Lanzhou 730070, China
Abstract: In the process of image restoration, in order to solve the problems of incomplete separation of color and texture optical properties and single dictionary design in image inpainting of sparse representation, which leads to the structural incoherence and blurring effect of mural image inpainting results, a mural inpainting method based on robust principal component analysis (RPCA) decomposition of block nuclear norm and entropy weighted clustering sparse representation was proposed. First, the proposed RPCA image decomposition algorithm based on block nuclear norm was used to decompose the mural image into a structure layer and a texture layer, and the block nuclear norm was used to perform texture correction operations, which could overcome the problem of incomplete separation of structure and texture by RPCA. Then, the entropy weighted k-means method was proposed to cluster the structure layer image and construct sparse sub-cluster dictionaries, and the reconstruction of the structure layer image was completed by sparse value decomposition and split Bregman iterative optimization method. Finally, the image of texture layer was repaired by using the bicubic interpolation algorithm, and the repaired structure layer and texture layer were fused to complete the repair of the damaged murals. Experimental results of digital restoration on the real Dunhuang murals show that the proposed method could effectively protect the edges, textures, and other important features in the mural image. In terms of visual quality and quantitative evaluation such as peak signal-to-noise ratio (PSNR), the proposed method had better performance than the comparison algorithms, and the restoration efficiency was higher.
Keywords: image reconstruction    mural inpainting    RPCA decomposition    block nuclear norm    clustering sparse representation    

敦煌莫高窟是世界上现存规模最宏大、内容最丰富的佛教石窟壁画宝库,其内所存的壁画、经卷等具有珍贵的研究价值。然而,由于自然风化的破坏及人为因素的影响,窟内壁画出现了地仗脱落、划痕、褪色、裂纹等严重的灾害,亟待保护[1]。因此,研究病态敦煌壁画的修复极其重要。但人工修复存在风险大、耗时长且修复不可逆等问题,将数字化虚拟修复应用于古代壁画的保护,是近年来国内外文化遗产保护的研究热点[2]

目前,传统图像修复方法主要分为3类:基于偏微分方程、基于样本块和基于稀疏表示的修复方法。基于偏微分方程的主要修复模型有TV(total variation)模型[3]和CDD(curvature-driven diffusion)模型[4],该类方法可扩散完成小区域破损图像的修复[5]。基于样本块的修复方法以Criminisi算法[6]为代表,该类算法可以完成较大区域破损图像的修复,但易出现填充错误的问题。陶兆胜等[7]在Criminisi算法的基础上,提出了基于边缘特征和像素结构相似度的图像修复算法,减少了像素误匹配率。

基于稀疏表示的修复方法利用字典中极少量原子及其系数的线性组合,完成对破损图像的重构[8]。Zhang等[9]提出了组稀疏图像修复方法,但该方法未考虑图像全局信息,会导致修复结果出现结构传播错误和线条断裂的现象。Zha等[10]提出将图像块与组稀疏联合稀疏表示的方法,但该方法也仅对局部固定大小区域进行联合块组字典建立。Ghorai等[11]提出一种考虑局部图像块统计信息和几何特征的多尺度金字塔稀疏表示修复方法。Zhao等[12]通过采用RPCA分解的方法,提出了一种结合阈值选择和中值滤波器的图像去噪算法。Qiang等[13]通过形态成分分析(morphological component analysis, MCA)方法,将图像分解为结构层和纹理层的修复方法,并应用于云南剑川石宝山壁画修复中。Hu等[14]提出一种基于图像分解的唐卡图像修复方法,后采用模糊集和局部二值模式完成修复。然而,上述稀疏表示算法字典设计单一,仅对局部固定大小区域构建字典,导致壁画图像修复结果易出现结构不连贯和模糊块效应等问题,此外在图像分解修复时,还存在颜色纹理光学属性分离不彻底的问题[15-16]

针对上述问题,本文提出了一种基于块核范数的RPCA分解与熵权类稀疏的壁画修复方法。主要工作有:1)首先,采用提出的基于块核范数的RPCA图像分解算法,通过引入凸先验信息块核范数,利用块核范数进行纹理矫正操作,将待修复图像分解为结构层和纹理层,更精确地完成结构和纹理信息的分离。2)然后,对结构层采用提出的熵权类稀疏表示进行修复,通过熵加权k-means算法对相似图像块进行聚类,构建稀疏子类字典,完成结构层图像的重构,克服了传统稀疏表示字典单一的问题。3)最后,对纹理层采用双三次插值进行修复,融合得到修复后的图像。通过对敦煌壁画的修复实验结果表明,本文方法较对比算法,获得了较好的主客观评价效果。

1 相关理论

稀疏表示就是用极少量原子和系数的线性组合对信号进行重构[17]。稀疏表示的关键在于寻找给定信号的最佳稀疏解,其数学定义为

$ \underset{\boldsymbol{\alpha}}{\arg \min }\|\boldsymbol{\alpha}\|_{0} \quad \text { s. t. } \boldsymbol{T}=\boldsymbol{D} \boldsymbol{\alpha} $ (1)

式中:α为稀疏系数,T为待观测信号,D为字典,‖ · ‖0l0范数,用来衡量α中非零元素个数。一般自然信号存在噪声或重构图像会出现一定的误差ε,因此,为了平衡图像的稀疏性和重构误差,式(1)可以改写为

$ \underset{\boldsymbol{\alpha}}{\arg \min } \frac{1}{2}\|\boldsymbol{T}-\boldsymbol{D} \boldsymbol{\alpha}\|_{2}+\lambda\|\boldsymbol{\alpha}\|_{0} \text { s.t. }\|\boldsymbol{T}-\boldsymbol{D} \boldsymbol{\alpha}\| \leqslant \varepsilon $ (2)

式中:$\frac{1}{2}{\left\| {\mathit{\boldsymbol{T}} - \mathit{\boldsymbol{D\alpha }}} \right\|_2}$为数据保真项,用来度量重构图像与原始图像的拟合程度;‖ α0为稀疏正则化项,用来约束图像的稀疏度;λ为平衡数据保真项和稀疏正则化项的参数。

在稀疏表示的图像修复中,图像退化的过程可定义为

$ \boldsymbol{Y}=\boldsymbol{J} \boldsymbol{X}+\boldsymbol{n} $ (3)

式中:Y为退化后的图像,X为未退化的原始图像,J为退化算子,n为高斯白噪声。因此,图像修复过程就是对退化图像Y进行推测估计得到原始图像X的过程。

根据式(2),在稀疏字典D给定的条件下,可得到稀疏表示图像修复的正则化模型为

$ \boldsymbol{\alpha}=\underset{\boldsymbol{\alpha}}{\arg \min } \frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{J} \cdot \boldsymbol{D} \boldsymbol{\alpha}\|_{2}^{2}+\lambda\|\boldsymbol{\alpha}\|_{1} $ (4)

通过式(4)求解稀疏系数α,最后通过Y = 得到修复后的图像。

2 本文算法 2.1 基于块核范数的RPCA图像分解

一般图像中既包含结构信息,也包含纹理信息,其中结构信息主要是图像边缘和轮廓等变化剧烈的信息部分;纹理信息主要是图像中亮度或灰度值变化缓慢或周期性变化区域,具有重复性、规则性和方向性,其反映了图像中同质现象的一种视觉特征[15, 18]。RPCA[19-20]是一种高维数据降维分析方法,通过RPCA可以将图像M (MR m×n)分解为两个矩阵的线性相加,一个是结构矩阵L (LR m×n,rank(L)≤m, n),主要用来表征图像的结构信息;另一个是纹理矩阵S (SR m×n),主要是图像的纹理细节信息,其图像分解模型定义为

$ \min \limits_{\boldsymbol{L}, \boldsymbol{S}} \operatorname{rank}(\boldsymbol{L})+\gamma\|\boldsymbol{S}\|_{0} \quad \text { s.t. } \boldsymbol{M}=\boldsymbol{L}+\boldsymbol{S} $ (5)

式中:‖ · ‖0表示l0范数,γ为平衡图像的结构和纹理分解的因子。

图 1敦煌莫高窟第92窟“衔花双鹿”的局部壁画为例进行RPCA分解说明。其中,图 1(a)为原始壁画,图 1(b)1(c)分别为RPCA分解后的结构层和纹理层,图 1(b)结构层中出现了模糊现象,图 1(c)纹理层也出现了分离不彻底的现象。从图 1实验可以发现:采用RPCA分解后纹理层中仍含有部分结构信息。这是因为采用RPCA空间变换分解后,破坏了纹理层的低秩性,导致结构信息仍残留在纹理层[21]

图 1 RPCA算法分解结果 Fig. 1 Decomposition results of RPCA algorithm

针对以上RPCA图像分解不彻底的问题,本文提出了一种基于块核范数的RPCA图像分解方法,以更好地分解结构与纹理信息。以图 2纹理图像为例,图像块2(a)和2(b)表现出全局不同但局部相似的特点,此时不满足全局重复、相似的低秩性要求。因此,为了将纹理层恢复到低秩性,通过对局部纹理图像块进行矫正操作,使其恢复为低秩矩阵。低秩矩阵的恢复是一个非凸连续函数优化问题,是一个NP难题,可以将求解最小秩的问题转化为核函数最小化的问题,即核函数为秩的凸包络,是矩阵秩的凸近似[16]

图 2 纹理的局部相似及全局不相似性示意 Fig. 2 Schematic of local similarity and global dissimilarity of texture

块核范数‖ HBNN是纹理部分的局部图像块,经过方向矫正和重叠取块操作后形成的矩阵奇异值之和,其定义为

$ \|\boldsymbol{H}\|_{\mathrm{BNN}}=\left\|\boldsymbol{P}_{m, \delta}\circ \boldsymbol{S}_{\theta}(\boldsymbol{H})\right\|_{\mathrm{m}^{*}}^{\mathrm{w}} $ (6)

式中:Pm, δ是周期扩展操作子,该操作子将大小为nv×nh的矩阵H进行周期扩展,扩展为$\frac{{{m_{\rm{v}}}}}{{{\delta _{\rm{v}}}}}{n_{\rm{v}}} \times \frac{{{m_{\rm{h}}}}}{{{\delta _{\rm{h}}}}}{n_{\rm{h}}}$的矩阵,m=(mv, mh)是块核范数的图像块大小,δ=(δv, δh)是垂直和水平的移位步长控制块重叠值;Sθ是矫正操作子,θ为矫正角度,θ∈[0, 45°],矫正方向为d={′t′, ′b′, ′l′, ′r′},分别代表上、下、左、右;‖·‖m*w为预核范数。

图 3为例对块核范数求解进行说明,其中,取nvnh为64,mvmh为16,δvδh为8,d=′r′,θ=45°。图中左侧矩阵H的原始坐标为(i, j),通过向右45°矫正操作后,变为秩为1的低秩矩阵,矫正操作后坐标为$\left( {\hat i, \hat j} \right)$,然后将Sθ(H)中重叠的矩阵块分离为不重叠的矩阵,即得到矩阵H的块核范数,如图 3中右侧用白色分割线分出的所有矩阵奇异值的和。

图 3 块核范数示意 Fig. 3 Schematic of block nuclear norm

将块核范数引入到RPCA图像分解模型,则得到如下的分解模型

$ \begin{gathered} \min \limits_{\boldsymbol{L}, \boldsymbol{t}_{i} \in \boldsymbol{Z}} \operatorname{rank}(\boldsymbol{L})+\gamma\left(\sum\limits_{i=1}^{K}\left\|t_{i}\right\|_{\mathrm{BNN}}+\right. \\ F_{z}\left(\varPhi\left(\boldsymbol{L}+\sum\limits_{i=1}^{K} t_{i}\right)\right) \quad \text { s.t. } \boldsymbol{M}=\boldsymbol{L}+\sum\limits_{i=1}^{K} t_{i} \end{gathered} $ (7)

式中:ti是由矫正操作子Sθ所确定的第i个子纹理,K为子纹理数,Z ={xRN|$\sum\limits_{i = 1}^N {{x_i} = 0} $}是所有零均值向量集,Fz是数据保真项,Fz(x)=tCz(x),tCz是闭凸集Cz={xRM|x=z的标示函数;Φ是线性度量操作子,为图像破损的掩膜标记。然后,采用交替乘子方向法(alternating direction method of multipliers,ADMM)[22]对式(7)进行迭代求解。

为了验证本文分解方法的有效性,以图 1所选取壁画为例进行实验,并与RPCA分解结果比较,见图 4。其中,图 4(a)(b)分别为RPCA分解和本文方法得到的结构层图像,比较发现:本文方法可以使结构层信息更加清晰;图 4(c)(d)为纹理层分解结果,发现通过本文方法分离的纹理和结构信息更加彻底。

图 4 引入块核范数的RPCA分解结果比较 Fig. 4 Comparison of RPCA decomposition results with block nuclear norm
2.2 熵加权k-means结构层图像聚类

对结构层图像进行类稀疏修复时,在k-means聚类算法的基础上,对欧氏距离进行熵加权改进作为新的聚类准则,提出了熵加权的k-means聚类算法,从而得到结构相似的子类图像。待分类图像块A的信息熵E(A)为

$ E(\boldsymbol{A})=-\sum\limits_{i=0}^{255} P_{\boldsymbol{A}}(i) \log _{2}\left[P_{\boldsymbol{A}}(i)\right] / 3 $ (8)

式中:PA (i)为图像块A中RGB颜色通道中第i个等级的像素占比率;i表示像素等级,范围为0~255。定义图像块A与聚类中心B的信息熵特征差异性系数为

$ q_{j}=|E(\boldsymbol{A})-E(\boldsymbol{B})| $ (9)

式中0≤qj≤1,qj越小,图像块AB越相似。然后,采用式(10)得到图像块A的熵权重vj

$ v_{j}=\frac{q_{j}}{\sum\limits_{j=1}^{k} q_{j}} $ (10)

式中k为质心个数,且满足$\sum\limits_{j = 1}^n {{v_j} = 1} $

从而得到图像块AB的熵权欧氏距离计算准则

$ d(\boldsymbol{A}, \boldsymbol{B})=\sqrt{v_{j}(\boldsymbol{A}-\boldsymbol{B})^{2}} $ (11)

通过式(11)将待分类图像块A划分到熵权欧氏距离最小的类中,从而构造得到结构层的子类图像。为了验证熵加权k-means聚类算法对壁画聚类的有效性,以图 4(b)结构图像为例进行实验,结果见图 5。其中,图 5(a)为原始图像,图 5(b)为聚类结果,图 5(c)(d)(e)为各子类图像,可以发现,具有相似结构的图像块能够较好地进行聚类,得到的结构层图像聚类更为准确。

图 5 熵加权的k-means算法聚类结果 Fig. 5 Clustering results of entropy weighted k-means algorithm
2.3 类稀疏表示的结构层图像修复

在结构层子类图像划分的基础上,对各子类进行类稀疏修复,定义类稀疏修复模型为

$ \boldsymbol{\alpha}_{k}=\underset{\boldsymbol{\alpha}_{k}}{\operatorname{argmin}} \frac{1}{2}\left\|\boldsymbol{y}_{k}-\boldsymbol{D}_{k} \boldsymbol{\alpha}_{k}\right\|_{2}^{2}+\lambda\left\|\boldsymbol{\alpha}_{k}\right\|_{1} $ (12)

式中:yk是第k类结构层图像块;Dkαk分别是对应的字典和稀疏系数。然后,采用奇异值分解和分裂Bregman迭代优化算法更新求解字典和系数[9]。将图像块yk的估计值rk进行奇异值分解,即

$ \boldsymbol{r}_{k}=\boldsymbol{U}_{k} \sum\nolimits_{k} \boldsymbol{V}_{k}^{\mathrm{T}} $ (13)

式中:UkVk分别为yk的左奇异正交矩阵和右奇异正交矩阵, ∑k为对角矩阵,则字典Dk中每个原子dki

$ \boldsymbol{d}_{k_{i}}=\boldsymbol{u}_{k_{i}} \boldsymbol{v}_{k_{i}}, i=1,2, \cdots, t $ (14)

式中ukivki分别为UkVk的列向量,t为第k类字典的原子数。所以,第k类自适应字典Dk

$ \boldsymbol{D}_{k}=\left\{\boldsymbol{d}_{k_{1}}, \boldsymbol{d}_{k_{2}}, \boldsymbol{d}_{k_{3}}, \cdots, \boldsymbol{d}_{k_{i}}\right\}, i=1,2, \cdots, t $ (15)

图 5结构层图像聚类结果为例,采用奇异值分解学习更新得到自适应子类字典,见图 6。其中,图 6(a)(b)(c)分别为图 5子类对应的自适应字典。可以发现,本文提出的类稀疏表示方法得到的各子类字典更能够体现图像的结构特征,且聚类字典训练更加高效。

图 6 自适应字典 Fig. 6 Adaptive dictionaries

在完成子类字典的更新后,接着对子类图像进行类稀疏系数求解

$ {\mathit{\boldsymbol{\alpha }}_k} = \mathop {\arg \min }\limits_{{\mathit{\boldsymbol{\alpha }}_k}} \frac{1}{2}\left\| {{\mathit{\boldsymbol{y}}_k} - \mathit{\boldsymbol{J}} \cdot {\mathit{\boldsymbol{D}}_k}{\mathit{\boldsymbol{\alpha }}_k}} \right\|_2^2 + \lambda {\left\| {{\mathit{\boldsymbol{\alpha }}_k}} \right\|_1} $ (16)

采用分裂Bregman迭代优化算法对式(16)求解,从而完成破损图像块的稀疏编码。最后,利用各子类图像的字典和系数完成结构层各子类图像的修复,并按照索引矩阵放回原始结构层图像,即得到结构层图像的修复结果YS

$ \boldsymbol{Y}_{\mathrm{S}}=\boldsymbol{D} \boldsymbol{\alpha}=\sum\limits_{i=1}^{k} \boldsymbol{D}_{i} \boldsymbol{\alpha}_{i} $ (17)
2.4 分层图像融合重构

采用双三次插值(bicubic interpolation)对纹理层破损区域进行插值修复,完成图像纹理细节信息的重建。双三次线性插值是通过对待修复点最近的16个采样点的加权平均完成像素点的估计,按下式计算

$ f(p, q)=\sum\limits_{i=0}^{3} \sum\limits_{j=0}^{3} f\left(p_{i}, q_{j}\right) W\left(p-p_{i}\right) W\left(q-q_{j}\right) $ (18)

式中:f(p, q)为待修复点坐标,W(p-pi)和W(q-qj)分别为待修复点的横坐标权重和纵坐标权重。权重的计算如下

$ W(u)=\left\{\begin{array}{lr} (a+2)|u|^{3}-(a+3)|u|^{2}+1, & |u| \leqslant 1 \\ a|u|^{3}-5 a|u|^{2}+8 a|u|-4 a, &1<|u|<2 \\ 0, & \text { otherwise } \end{array}\right. $ (19)

最后将结构层的修复结果YS与纹理层的修复结果YT融合,得到最终的修复图像Y

$ \boldsymbol{Y}=\boldsymbol{Y}_{\mathrm{S}}+\boldsymbol{Y}_{\mathrm{T}} $ (20)
2.5 算法步骤

Step1    输入待修复壁画图像,并获得其掩膜图像。

Step2    将待修复壁画图像利用块核范数的RPCA方法进行分解,得到结构层和纹理层。

Step3    对结构层图像采用熵加权的k-means算法进行聚类,对每个子类图像构造类稀疏字典,并采用奇异值分解和分裂Bregman进行类稀疏表示修复。

Step4    对纹理层图像进行双三次插值修复。

Step5    将结构层和纹理层的修复结果融合操作,得到修复后的壁画图像。

3 实验结果与分析

实验运行软件环境为Windows 10操作系统,采用MATLAB R2016a软件,硬件配置为Intel(R) Core i7-9700K CPU@3.6 GHz,16.0 GB RAM, NVIDIA GeForce GTX 1660。修复结果采用主观和客观两种方式进行评价,客观评价使用峰值信噪比(peak signal to noise ratio, PSNR)、结构相似性(structural similarity index measurement, SSIM)和算法时间复杂度进行评价。此外,为了验证本文方法的有效性,采用人为添加破损和真实破损的敦煌壁画进行修复实验比较,并与经典修复CDD算法[4]、Criminisi算法[6]、文献[7]和文献[9]的修复结果进行对比分析。

3.1 人为添加破损壁画修复

考虑到壁画修复是一个病态逆问题, 为了定量分析修复算法的有效性,一般对完好壁画图像人工添加破损后进行定量分析。首先对壁画进行人为添加划痕修复对比实验,如图 7所示,从上往下分别为“北壁弥勒经变”局部壁画、“双笛合奏”局部壁画和“观无量寿经变”局部壁画。其中,第1列为原始图像,第2列为掩膜图像,第3列为CDD算法修复结果,从图 7(c)可以发现,CDD算法存在修复不彻底的问题,如第二幅壁画左侧矩形框内出现了明显的修复痕迹。第4列为Criminisi算法修复结果,出现了块匹配错误和结构传播错误的现象,如图 7(d)第三幅壁画的左侧眉毛部分和右侧矩形框中均出现了误匹配现象。第5列为文献[7]修复结果,图 7(e)第三幅壁画的左侧眉毛的连续性较好,但第一幅壁画的修复结果中仍存在像素匹配错误的问题。第6列为文献[9]修复结果,其修复结果较符合人类的视觉感受,但因为该算法仅考虑了图像局部相关性,导致修复结果中出现了线条断裂和模糊的现象,如图 7(f)第三幅壁画的眉毛部分出现了线条不连贯的问题。最后一列图 7(g)为本文修复结果,与其他算法比较,本文算法修复后图像更加连贯自然,获得了更好的视觉效果。

图 7 不同算法对人为添加破损壁画的修复结果对比 Fig. 7 Comparison of mural inpainting results of artificial damages by different algorithms

为了进一步对图 7的修复结果做出定量评价,采用PSNR和SSIM进行比较,见表 1。PSNR和SSIM是衡量图像质量的重要定量评价指标,其中PSNR值越大,代表失真越少,即修复效果越好;SSIM是一种衡量两幅图像相似度的指标,该值越大,表明修复后图像的失真程度越小[18]。从表 1中可以发现,本文算法在PSNR和SSIM上均优于其他对比算法,表明本文方法修复质量更好,从而验证了本文方法对于人工破损壁画修复的有效性。

表 1 不同算法修复结果PSNR和SSIM对比 Tab. 1 Comparison of PSNR and SSIM results of different algorithms

下面对不同算法时间复杂度进行比较分析,以待修复壁画图像大小R×W为例。CDD算法[4]通过计算待修复点的梯度和曲率扩散完成修复,修复过程共迭代n次,则其时间复杂度为O(R×W×n),因在时间复杂度分析时一般将常数忽略,进一步简化后,其时间复杂度为O(n)。Criminisi算法一般将待修复图像块像素取值为9×9[6],然后在完好区域选取与待修复块大小相同的r个匹配块,通过计算欧氏距离得到最佳匹配块,从而完成修复,修复过程迭代n2次,则其时间复杂度为Ο(R×W×r×n2×92)=Ο(r×n2)。文献[7]利用样本相似度和信息熵确定最佳样本块集合,并采用果蝇优化算法寻找最优匹配块,果蝇算法在迭代j次后完成寻优,则文献[7]时间复杂度为Ο(R×W×r×n2×92×j)= Ο(r×j×n2)。文献[9]采用组稀疏修复方法,将图像划分为Bs×Bs大小的图像块,遍历n2次图像块,然后在L×L的邻域中得到相似结构组,并进行n次字典和系数迭代更新,则文献[9]时间复杂度为${\rm{O}}\left( {{n^2} \times \frac{R}{{{B_{\rm{s}}}}} \times \frac{W}{{{B_{\rm{s}}}}} \times \frac{L}{{{B_{\rm{s}}}}} \times \frac{L}{{{B_{\rm{s}}}}} \times \mathit{n}} \right) = {\rm{O}}\left( {{n^3}} \right)$。本文将结构图像划分为Bs×Bs大小的图像块,采用熵权k-means分类算法遍历n2次后,划分为k类完成修复,则本文算法时间复杂度为${\rm{O}}\left( {{n^2} \times \frac{R}{{{B_{\rm{s}}}}} \times \frac{W}{{{B_{\rm{s}}}}} \times \mathit{k}} \right) = {\rm{O}}\left( {k \times {n^2}} \right)$,综上不同算法时间复杂度见表 2

表 2 不同算法修复时间复杂度对比 Tab. 2 Comparison of inpainting time complexity of different algorithms

通过表 2可以发现:CDD算法时间复杂度最低,但CDD算法修复效果较差,无法彻底完成修复;文献[9]组稀疏时间复杂度最高,修复时间耗时最长;Criminisi算法和文献[7]时间复杂度均低于文献[9],但远高于本文算法,这是因为文献[7]中果蝇优化迭代j值以及Criminisi算法中匹配块搜索值r远大于本文划分类k值。本文算法时间复杂度仅高于CDD算法,远小于其他3种比较算法,说明本文在达到较好修复效果的同时,算法执行效率也较高。

3.2 真实破损壁画修复

为了进一步验证本文算法的有效性,采用4组真实破损壁画图像进行修复实验,见图 8。其中,图 8(a)为真实壁画图像,从左往右分别为“第263窟·人字坡·供养菩萨”、“第61窟·东壁维摩诘经变·骑队”、“第285窟·北壁四龛楣·供养菩萨”和“第220窟·飞天”的局部壁画,图 8(b)为掩膜图像。图 8(c)为CDD修复结果,从第一幅和第三幅可以看出该算法出现了结构模糊的问题。图 8(d)为Criminisi修复结果,第二幅和第三幅壁画右侧矩形框内的供养菩萨头光部分出现了结构传播错误问题,第四幅壁画出现了块匹配错误。图 8(e)为文献[7]修复结果,同样在第三幅和第四幅壁画中存在块匹配错误;图 8(f)为文献[9]组稀疏修复结果,从第一幅壁画图像和第四幅壁画矩形框内箭头所指区域存在修复痕迹,第二幅和第三幅壁画出现了线条断裂现象。图 8(g)为本文结果,可以看出,在线条连续性和清晰性方面本文修复效果均有所提高。

图 8 不同算法对真实破损壁画的修复效果对比 Fig. 8 Comparison of mural inpainting results of real damages by different algorithms
4 结论

1) 提出了一种基于块核范数的RPCA图像分解方法,将块核范数引入RPCA图像分解模型,对RPCA分解后的纹理层通过矫正操作,克服了图像结构层与纹理层分离不彻底的问题。

2) 提出了熵权类稀疏的修复方法,采用熵加权的k-means聚类算法得到结构相似的子类图像,并通过奇异值分解和分裂Bregman的类稀疏表示完成结构层图像的重构,克服了单一字典学习出现的结构传播错误和模糊现象。

3) 破损敦煌壁画修复实验表明,本文方法取得了较好的视觉效果和较高的定量评价,修复效率也进一步提高。

4) 尽管本文方法具有较好的修复效果,但未对壁画图像的历史、宗教文化等语义信息进行考虑,后续将采用语义分析等深度学习方法进一步研究。

参考文献
[1]
付心仪, 麻晓娟, 孙志军. 破损壁画的数字化复原研究—以敦煌壁画为例[J]. 装饰, 2019(1): 21.
FU Xinyi, MA Xiaojuan, SUN Zhijun. Digital restoration of damaged murals: based on Dunhuang murals[J]. Art & Design, 2019(1): 21. DOI:10.16272/j.cnki.cn11-1392/j.2019.01.007
[2]
WANG Huan, LI Qingquan, JIA Sen. A global and local feature weighted method for ancient murals inpainting[J]. International Journal of Machine Learning and Cybernetics, 2020, 11(6): 1197. DOI:10.1007/s13042-019-01032-2
[3]
SHEN J H, CHAN T F. Mathematical models for local nontexture inpaintings[J]. SIAM Journal on Applied Mathematics, 2002, 62(3): 1019. DOI:10.1137/S0036139900368844
[4]
CHAN T F, SHEN J H. Nontexture inpainting by curvature-driven diffusions[J]. Journal of Visual Communication and Image Representation, 2001, 12(4): 436. DOI:10.1006/jvci.2001.0487
[5]
CAO Xin, SUN Yi, KANG Fei, et al. A novel denoising framework for Cerenkov luminescence imaging based on spatial information improved clustering and curvature-driven diffusion[J]. Journal of Innovative Optical Health Sciences, 2018, 11(4): 1850017. DOI:10.1142/S1793545818500177
[6]
CRIMINISI A, PEREZ P, TOYAMA K. Region filling and object removal by exemplar-based image inpainting[J]. IEEE Transactions on Image Processing, 2004, 13(9): 1200. DOI:10.1109/TIP.2004.833105
[7]
陶兆胜, 张敬寒, 王磊, 等. 基于边缘特征和像素结构相似度的图像修复算法[J]. 计算机辅助设计与图形学学报, 2019, 31(10): 1768.
TAO Zhaosheng, ZHANG Jinghan, WANG Lei, et al. Image inpainting algorithm based on edge feature and pixel structure similarity[J]. Journal of Computer-Aided Design & Computer Graphics, 2019, 31(10): 1768. DOI:10.3724/SP.J.1089.2019.17671
[8]
张杰, 史小平, 张焕龙, 等. 高噪声遥感图像稀疏去噪重建[J]. 哈尔滨工业大学学报, 2019, 51(10): 47.
ZHANG Jie, SHI Xiaoping, ZHANG Huanlong, et al. High noise remote sensing image sparse denoising reconstruction[J]. Journal of Harbin Institute of Technology, 2019, 51(10): 47. DOI:10.11918/j.issn.0367-6234.201806161
[9]
ZHANG Jian, ZHAO Debin, GAO Wen. Group-based sparse representation for image restoration[J]. IEEE Transactions on Image Processing, 2014, 23(8): 3336. DOI:10.1109/TIP.2014.2323127
[10]
ZHA Zhiyuan, YUAN Xin, WEN Bihan, et al. Image restoration using joint patch-group-based sparse representation[J]. IEEE Transactions on Image Processing, 2020, 29: 7735. DOI:10.1109/TIP.2020.3005515
[11]
GHORAI M, SAMANTA S, MANDAL S, et al. Multiple pyramids based image inpainting using local patch statistics and steering kernel feature[J]. IEEE Transactions on Image Processing, 2019, 28(11): 5495. DOI:10.1109/TIP.2019.2920528
[12]
ZHAO Chenyi, YAN Yue, WANG Yanme, et al. White spots noise removal of neutron images using improved robust principal component analysis[J]. Fusion Engineering and Design, 2020(156): 111739. DOI:10.17632/ww4gzym4yj.1
[13]
QIANG Zhenpin, HE Libo, CHEN Yaqiong, et al. Research on mural inpainting method based on MCA image decomposition[C]//Proceedings of the 2nd International Conference on Intelligent Information Processing. New York: Association for Computing Machinery, 2017, 23: 1. DOI: 10.1145/3144789.3144816
[14]
HU Wenjin, LIU Zhongmin, YE Yuqi. A new quality assessment for Thangka image inpainting[J]. Concurrency Computation, 2018, 30(22): 1. DOI:10.1145/3144789.3144816
[15]
FAN Yaru, HUANG Tingzhu, MA Tianhui, et al. Cartoon-texture image decomposition via non-convex low-rank texture regularization[J]. Journal of the Franklin Institute, 2017, 354: 3170. DOI:10.1109/TIP.2014.2299067
[16]
ONO S, MIYATA T, YAMADA I. Cartoon-texture image decomposition using blockwise low-rank texture characterization[J]. IEEE Transactions on Image Processing, 2014, 23(3): 1128. DOI:10.1109/TIP.2014.2299067
[17]
GUO Kanghui, LABATE D, AYLLON J R P. Image inpainting using sparse multiscale representations: image recovery performance guarantees[J]. Applied and Computational Harmonic Analysis, 2020(49): 343. DOI:10.1016/j.acha.2020.05.001
[18]
陈永, 郭红光, 艾亚鹏. 基于双域分解的多尺度深度学习单幅图像去雾[J]. 光学学报, 2020, 40(2): 0210003.
CHEN Yong, GUO Hongguang, AI Yapeng. Single image dehazing of multi-scale deep learning based on dual domain decomposition[J]. Acta Optica Sinica, 2020, 40(2): 0210003. DOI:10.3788/AOS202040.0210003
[19]
CANDES E J, LI Xiaodong, MA Yi, et al. Robust principal component analysis?[J]. Journal of the ACM, 2011, 58(3): 1. DOI:10.1145/1970392.1970395
[20]
BOUWMANS T, JAVED S, ZHANG Hongyang, et al. On the applications of robust PCA in image and video processing[J]. Proceedings of the IEEE, 2018, 106(8): 1427. DOI:10.1109/JPROC.2018.2853589
[21]
XIE Ting, LI Shutao, SUN Bin. Hyperspectral images denoising via nonconvex regularized low-rank and sparse matrix decomposition[J]. IEEE Transactions on Image Processing, 2020, 29: 44. DOI:10.1109/TIP.2019.2926736
[22]
LI Peng, CHEN Wengu, NG M K. Compressive total variation for image reconstruction and restoration[J]. Computers & Mathematics with Applications, 2020, 80(5): 874. DOI:10.1016/j.camwa.2020.05.006