多尺度輪廓波分解的群稀疏壁畫修復(fù)算法

陳 永，趙夢雪，陶美風(fēng)

(1.蘭州交通大學(xué) 電子與信息工程學(xué)院，甘肅 蘭州 730070；2.甘肅省人工智能與圖形圖像處理工程研究中心，甘肅 蘭州 730070)

敦煌莫高窟壁畫是享譽世界的文化瑰寶，具有極高的研究價值。然而，由于自然氣候腐蝕、人為破壞等因素導(dǎo)致大量壁畫出現(xiàn)空鼓、起甲、酥堿等亟待解決的病害問題。數(shù)字化修復(fù)可以克服人工物理修復(fù)風(fēng)險大且不可逆的問題，是目前的研究熱點[1]。

目前，傳統(tǒng)圖像修復(fù)方法可以分為：基于序列修復(fù)[2]和基于稀疏表示修復(fù)[3]。其中，基于序列修復(fù)方法又包括擴散法[4-5]和樣本合成法[6-7]。擴散法將完好信息向破損區(qū)域動力學(xué)擴散，適用于細(xì)小區(qū)域的修復(fù)。樣本合成法則是根據(jù)塊匹配原則，在完好區(qū)域選擇相似樣本塊對破損區(qū)域進(jìn)行復(fù)制填充，但易出現(xiàn)樣本塊填充錯誤等問題?；谙∈璞硎拘迯?fù)方法則是在過完備字典上計算破損塊的稀疏編碼，從而將稠密樣本轉(zhuǎn)化為模型復(fù)雜度較低的稀疏表達(dá)形式，借助稀疏編碼和過完備字典對受損區(qū)域進(jìn)行重構(gòu)[8-9]。ZHANG等[10]提出了組稀疏圖像修復(fù)方法，通過局部自相似性選擇匹配塊來構(gòu)建字典，但該方法未考慮全局圖像特征，易出現(xiàn)結(jié)構(gòu)斷裂問題。CHEN等[11]將自適應(yīng)分組稀疏表示引入到圖像修復(fù)中，但該方法采用歐式距離計算樣本塊之間的分組差異，未考慮圖像紋理和結(jié)構(gòu)信息的差異性。ZHA等[12]提出樣本塊聯(lián)合組稀疏重構(gòu)破損圖像的方法，但該方法建立聯(lián)合字典的過程也僅在局部區(qū)域進(jìn)行字典學(xué)習(xí)，導(dǎo)致修復(fù)結(jié)果存在紋理模糊的問題。馮象初等[13]利用博弈理論和稀疏先驗性提出了一種圖像修復(fù)模型，但該方法在完好區(qū)域與破損區(qū)域修復(fù)博弈時，未考慮樣本塊局部自相似性特點。

綜上所述，現(xiàn)有稀疏表示在壁畫圖像修復(fù)時，未考慮壁畫結(jié)構(gòu)信息與紋理信息的差異性，導(dǎo)致修復(fù)結(jié)果易出現(xiàn)紋理模糊和結(jié)構(gòu)線條斷裂等問題。針對上述問題，筆者提出了一種基于多尺度輪廓波分解的群稀疏壁畫修復(fù)算法。主要工作有：① 通過非下采樣輪廓波變換對待修復(fù)壁畫圖像進(jìn)行多尺度分解，將其分解為紋理重復(fù)低頻分量和結(jié)構(gòu)輪廓高頻分量，克服了現(xiàn)有稀疏表示壁畫修復(fù)時，未考慮壁畫結(jié)構(gòu)和紋理信息差異性的不足；② 對低頻分量采用提出的改進(jìn)群稀疏算法，完成對壁畫紋理分量的修復(fù)；③ 采用三次立方卷積插值算法對高頻分量壁畫結(jié)構(gòu)信息進(jìn)行插值修復(fù)，完成壁畫結(jié)構(gòu)信息的修復(fù)；④ 采用非下采樣輪廓波逆變換實現(xiàn)對修復(fù)后各尺度分量的融合重構(gòu)，從而完成壁畫修復(fù)。通過對敦煌壁畫修復(fù)實驗，結(jié)果證明所提方法與對比方法相比，在主觀及客觀定量指標(biāo)評價中均獲得了更好的修復(fù)效果。

1 文中算法

敦煌壁畫中不僅包含灰度值過渡平緩、較規(guī)則的紋理信息，還包含豐富的線條輪廓等像素值變化劇烈的結(jié)構(gòu)信息[14]。為了充分利用壁畫紋理和結(jié)構(gòu)信息，筆者提出了一種基于多尺度輪廓波分解的群稀疏壁畫修復(fù)算法，將壁畫分解為紋理重復(fù)變化緩慢的低頻分量和結(jié)構(gòu)輪廓變化劇烈的高頻分量，從而將壁畫修復(fù)視點從像素域轉(zhuǎn)移至頻域，利用頻率差異[15]對輸入圖像進(jìn)行多尺度分解，繼而對各分解分量分別進(jìn)行修復(fù)，并最終融合重構(gòu)，從而恢復(fù)圖像原始頻率分量。

1.1 基于非下采樣輪廓波變換的多尺度壁畫分解

非下采樣輪廓波變換(NonSubsampled Contourlet Transform，NSCT)使用非下采樣操作，使分解圖像具備位移不變、多尺度、多方向等特性，消除了因缺乏位移不變性出現(xiàn)在奇異點附近的偽吉布斯現(xiàn)象，其分解系數(shù)對曲線、奇異點等細(xì)節(jié)紋理信息具有更稀疏的表達(dá)[16]。首先，通過由雙通道二維濾波器組和采樣矩陣構(gòu)成的非下采樣金字塔濾波器(NonSubsampled Pyramid，NSP)對尺寸大小為M×M的壁畫圖像進(jìn)行a層多分辨率特征提取，將待修復(fù)壁畫圖像分解為高頻分量和低頻分量。對得到的高頻分量繼續(xù)進(jìn)行方向濾波，而低頻分量則通過金字塔NSP濾波器繼續(xù)進(jìn)行下一尺度迭代分解，直至達(dá)到設(shè)定分解層數(shù)。以2尺度2方向分解為例，其分解結(jié)構(gòu)框架如圖1所示。

圖1 壁畫分解結(jié)構(gòu)框架示意圖

與其他圖像單一分解方法不同，為了獲得豐富的結(jié)構(gòu)信息，筆者對圖1中分解得到的高頻分量，繼續(xù)采用由雙通道扇形濾波器組和梅花矩陣構(gòu)成的非下采樣方向濾波器組(NonSubsampled Directional Filter Bank，NSDFB)進(jìn)行兩級方向濾波疊加，其四通道濾波結(jié)構(gòu)圖以及濾波頻域分布如圖2所示。第1級為扇形濾波器組，首先對高頻分量進(jìn)行垂直方向U0(z)和水平方向U1(z)濾波，然后分別將濾波結(jié)果送入第2級象限濾波器中，與各象限濾波方向的重疊部分即為該方向子圖的頻域分布范圍，輸出重疊區(qū)域的方向子圖yk，完成方向濾波過程。經(jīng)過濾波各個尺度的高通分量分別被分解為2bi(i=1，2，…，a)個方向子帶圖像，i為分解層數(shù)，bi為當(dāng)前層高頻子帶濾波方向數(shù)，各方向頻域分布如圖2(b)所示。

(a) 四通道濾波結(jié)構(gòu)圖

該濾波結(jié)構(gòu)輸出等價公式可表示為

(1)

NSP和NSDFB對圖像的疊加過程中均無下采樣環(huán)節(jié)，使其獲得平移不變性，以圖3敦煌莫高窟第57窟被譽為唐代美人的“菩薩圖”的壁畫分解結(jié)果為例，可以看出，分解后的低頻分量圖3(b)中分布著大量像素變化緩慢的紋理特征；而高頻分量圖3(c)～圖3(h)則攜帶著明顯的圖像變化劇烈的結(jié)構(gòu)特征，并包含了壁畫圖像的主要輪廓結(jié)構(gòu)信息，其中圖3(c)和圖3(d)分別對應(yīng)第一尺度扇形濾波器組水平方向和垂直方向的濾波結(jié)果，圖3(e)～圖3(h)為二尺度四個方向的結(jié)構(gòu)特征，較細(xì)尺度分解子圖則依據(jù)濾波劃分頻帶提取四個方向的結(jié)構(gòu)特征，將分布在同方向的奇異點合成NSCT分解系數(shù)。從圖3可以發(fā)現(xiàn)，經(jīng)過雙通道扇形濾波器組和梅花矩陣對其進(jìn)行聯(lián)合作用后，可以在多個方向的高頻分量上提取到充分的結(jié)構(gòu)特征，可以在重構(gòu)壁畫時具有足夠完備的結(jié)構(gòu)特征信息選擇。

圖3 壁畫分解結(jié)果示意圖

1.2 改進(jìn)群稀疏的壁畫低頻分量修復(fù)

在采用NSCT對壁畫分解為低頻分量和高頻分量的基礎(chǔ)上，其低頻分量為大量重復(fù)和規(guī)則變化的紋理信息，根據(jù)其局部自相似性和全局不同的特點，采用改進(jìn)群稀疏方法對其進(jìn)行重構(gòu)修復(fù)。組稀疏算法以圖像結(jié)構(gòu)組作為稀疏表示的原子，在求解模型時首先考慮樣本塊之間的相關(guān)性，即構(gòu)造相似圖像塊集合，相似集合由結(jié)構(gòu)紋理相似并且尺寸一致的匹配塊矩陣構(gòu)成[3，10]。但傳統(tǒng)組稀疏算法構(gòu)建結(jié)構(gòu)組時使用歐式距離作為衡量匹配塊相似度的基準(zhǔn)，衡量標(biāo)準(zhǔn)單一，未考慮圖像紋理和結(jié)構(gòu)信息的差異性。

針對上述不足，根據(jù)分解后壁畫低頻分量為大量重復(fù)和局部相似的特征，以及紋理信息像素之間存在強相關(guān)性的特點，在構(gòu)造相似集合時，提出將結(jié)構(gòu)相似性作為匹配原則，將圖像結(jié)構(gòu)信息以及對比度、亮度等更多衡量因素考慮進(jìn)去，充分利用先驗信息來提高字典的構(gòu)建準(zhǔn)確性。如圖4所示，對低頻分量壁畫圖像群稀疏修復(fù)時，首先將壁畫的低頻分量劃分為n個N×N大小的重疊圖像塊Xi；然后將Xi表示為Bi×1維列矩陣，形成待匹配樣本塊候選集；接著計算匹配塊之間的結(jié)構(gòu)相似性，選擇匹配度最高的c個樣本塊組成集合SXi}。其結(jié)構(gòu)相似性的計算公式為

(2)

再將集合{SXi}內(nèi)的圖像塊橫向展開即轉(zhuǎn)化為大小為Bi×c的相似群XGi，其中XGi為第Xi個圖像塊對應(yīng)的相似群。

在構(gòu)造完字典后，進(jìn)行字典學(xué)習(xí)與系數(shù)更新。對于已經(jīng)構(gòu)建的群稀疏模型，每個樣本塊都由相似群域中的原子結(jié)合稀疏系數(shù)近似重構(gòu)，通過rGi奇異值分解估計各相似群字典，即

(3)

(4)

圖4 壁畫結(jié)構(gòu)組構(gòu)造示意圖

圖5 低頻分量與訓(xùn)練字典

其中，DGi為相似群XGi的對應(yīng)字典；dGi×k為第i個群字典中的第k個原子列，k=1，2，…，n；uGi、vGi分別為UGi、VGi的列向量。通過提出的多尺度輪廓波分解變換，將大量的紋理信息保留在低頻分量上，而將結(jié)構(gòu)信息分離在高頻分量上，從而避免了傳統(tǒng)組稀疏表示方法在字典構(gòu)造時，未考慮結(jié)構(gòu)和紋理信息差異的問題。對于分解后的低頻分量，采用式(4)構(gòu)建得到的訓(xùn)練字典，如圖5所示。

對每個樣本塊相似群XGi的訓(xùn)練字典DGi，采用分裂伯格曼算法求解其對應(yīng)的稀疏向量αGi=[αGi?1，αGi?2，…，αGi?c]，其數(shù)學(xué)定義為

(5)

(6)

1.3 壁畫高頻分量修復(fù)及融合重構(gòu)

對分解后的高頻結(jié)構(gòu)信息，采用三次立方卷積插值算法進(jìn)行插值修復(fù)。利用鄰域已知點像素值進(jìn)行卷積操作，增強其本身攜帶的點奇異特征信息，提高各尺度高頻分量的分辨率。三次立方卷積插值算法利用浮點坐標(biāo)(i+u，j+v)鄰域內(nèi)16個坐標(biāo)像素值，通過下式求得待修復(fù)像素值：

f(i+u，j+v)=[A][B][C] ，

(7)

其中，f(i+u，j+v)為待修復(fù)坐標(biāo)像素值，[A]=[S(u+1)S(u)S(u-1)S(u-2)]，[B]為目標(biāo)點16個鄰點像素值構(gòu)成的4×4維矩陣，[C]=[S(v+1)S(v)S(v-1)S(v-2)]T，S是對sin(xπ)/(xπ)的逼近函數(shù)。

最后將低頻分量的修復(fù)結(jié)果圖XL與各尺度各方向高頻分量的插值修復(fù)結(jié)果圖XH進(jìn)行逆變換得到最終修復(fù)圖像

X=XL+XH。

(8)

2 實驗結(jié)果與分析

為了驗證所提方法的有效性，進(jìn)行不同修復(fù)方法對比實驗。對比實驗均在相同的軟硬件環(huán)境下進(jìn)行：MATLAB R2016a編程實現(xiàn)，硬件環(huán)境為Intel(R)Core(TM)i5-1155G7 @ 2.50 GHz，16.0 GB RAM，NVIDIA GeForce MX450。采用評價準(zhǔn)則峰值信噪比PSNR和結(jié)構(gòu)相似性SSIM分別進(jìn)行人為添加破損和真實敦煌壁畫修復(fù)實驗對比及客觀性量化評價，并與曲率驅(qū)動擴散(Curvature Driven Diffusions，CDD)算法、Criminisi圖像修復(fù)算法、文獻(xiàn)[10]組稀疏和文獻(xiàn)[12]聯(lián)合樣本塊稀疏表示方法進(jìn)行對比分析。

2.1 人為添加破損壁畫修復(fù)

選擇4幅完好壁畫圖像經(jīng)人為添加破損后進(jìn)行修復(fù)比較，結(jié)果如圖6所示。圖6(a)為完好壁畫，圖6(b)為添加破損圖像。圖6(c)為CDD算法修復(fù)結(jié)果，因該算法采用流體擴散修復(fù)，其擴散能力有限，出現(xiàn)了修復(fù)不徹底的現(xiàn)象，如壁畫1右上菩薩頭光部分未完成修復(fù)。圖6(d)為 Criminisi算法結(jié)果，該算法采用塊匹配復(fù)制填充的方法，易出現(xiàn)填充錯誤的現(xiàn)象，如壁畫1左側(cè)矩形框內(nèi)結(jié)構(gòu)線條出現(xiàn)了錯誤匹配，壁畫4矩形框也出現(xiàn)了紋理紊亂的現(xiàn)象。圖6(e)為文獻(xiàn)[10]組稀疏修復(fù)結(jié)果，因該方法未考慮壁畫紋理和結(jié)構(gòu)的差異性，壁畫中線條信息修復(fù)后存在斷裂及不連貫的問題，如圖6(e)壁畫3“狩獵·野牛圖”的身體部分出現(xiàn)線條結(jié)構(gòu)斷裂。圖6(f)為文獻(xiàn)[12]聯(lián)合樣本塊稀疏修復(fù)結(jié)果，在壁畫1和壁畫4中存在結(jié)構(gòu)斷裂和塊狀模糊問題。圖6(g)為文中算法結(jié)果，與其他算法相比，修復(fù)后圖像線條連貫更自然，獲得了更好的視覺效果。

為了進(jìn)一步定量比較分析，采用PSNR和SSIM對圖6的修復(fù)結(jié)果定量比較，結(jié)果見表1。其中PSNR和SSIM分別表示圖像失真程度和結(jié)構(gòu)相似性對比，其值越大，說明失真越小，修復(fù)效果越好。從表1可以看出，文中方法在 PSNR 和 SSIM 上均優(yōu)于對比算法，表明筆者所提方法修復(fù)質(zhì)量更好。

表1 人為添加破損壁畫修復(fù)客觀比較

2.2 真實破損壁畫修復(fù)

下面繼續(xù)對4組真實破損壁畫圖像進(jìn)行修復(fù)對比實驗，修復(fù)結(jié)果如圖7所示。其中，圖7(c)為CDD算法修復(fù)結(jié)果，對于壁畫2和壁畫4，該算法均無法完成擴散修復(fù)。圖7(d)為Criminisi算法修復(fù)結(jié)果，其中壁畫1和壁畫3矩形框內(nèi)出現(xiàn)修復(fù)線條錯位以及塊匹配出錯現(xiàn)象，導(dǎo)致修復(fù)結(jié)果平滑度欠缺，不滿足人眼視覺效果。圖7(e)為文獻(xiàn)[10]組稀疏修復(fù)結(jié)果，該方法對于細(xì)小裂紋有較好的修復(fù)效果，但對于較大區(qū)域其修復(fù)結(jié)果中會出現(xiàn)結(jié)構(gòu)模糊的問題，如壁畫4矩形框內(nèi)存在模糊和暈染現(xiàn)象。圖7(f)為文獻(xiàn)[12]修復(fù)結(jié)果，該方法存在線條不連續(xù)的和塊效應(yīng)問題，如壁畫3和壁畫4矩形框內(nèi)出現(xiàn)了結(jié)構(gòu)斷裂和紋理模糊的現(xiàn)象。圖7(g)為文中算法修復(fù)結(jié)果，因采用結(jié)構(gòu)和紋理分解修復(fù)，其對真實壁畫的修復(fù)效果相較對比方法，取得了較好的視覺修復(fù)效果。

圖7 真實破損壁畫修復(fù)實驗對比圖

對于真實破損壁畫修復(fù)結(jié)果的評價，因缺少相應(yīng)的標(biāo)準(zhǔn)參考圖像，故采用無參考評價方法進(jìn)行評價。首先，采用無參考主觀評價指標(biāo)(Mean Opinion Scores，MOS)進(jìn)行評價[17]。MOS值代表主觀評分值，其值越大，表明修復(fù)效果越滿足人眼視覺效果，其與破損量對應(yīng)關(guān)系見表2。為了對主觀評價進(jìn)行量化分析，進(jìn)一步通過計算皮爾森線性相關(guān)系數(shù)(Pearson Linear Correlation Coefficient，PLCC)進(jìn)行量化評價。通過PLCC系數(shù)計算MOS值與信息熵的一致程度，該值越大，代表修復(fù)效果越好，比較結(jié)果見表3。從表3中可以看出，筆者所提方法修復(fù)效果優(yōu)于對比算法。

表2 真實破損壁畫修復(fù)主觀評價及量化表

表3 真實破損壁畫修復(fù)結(jié)果主觀評價對比

在完成真實壁畫主觀評價后，下面進(jìn)一步采用無參考客觀評價指標(biāo)(Sum of Modulus of gray Difference，SMD)和信息熵對圖7修復(fù)結(jié)果進(jìn)行量化比較，結(jié)果見表4。其中SMD和信息熵分別代表灰度方差函數(shù)和圖像信息豐富程度，二者值越大，表明修復(fù)后信息越豐富，修復(fù)結(jié)果更清晰。從表4看出，文中方法的SMD和信息熵評價值均優(yōu)于對比算法，進(jìn)一步驗證了文中方法對于真實破損壁畫修復(fù)的有效性。

表4 真實破損壁畫修復(fù)客觀比較

2.3 算法泛化性能分析

通過平均準(zhǔn)確率對不同方法進(jìn)行泛化能力對比。為了便于比較，采用自建數(shù)據(jù)集中有參照人為添加破損壁畫進(jìn)行量化比較。在平均準(zhǔn)確率度量時，使用均方誤差 (Mean Square Error，MSE)和平均峰值信噪比(Mean Peak Signal Noise Ratio，MPSNR)進(jìn)行評價，二者均表示修復(fù)后壁畫圖像的失真程度，其中MSE值越小，代表修復(fù)誤差越小，MPSNR值越大，表明修復(fù)效果更好，泛化能力更強。各算法比較結(jié)果如表5 所示?？梢钥闯觯?種方法中，文中方法均方根誤差值最小，且MPSNR值最大，說明文中方法修復(fù)后，壁畫失真最少，其泛化性能優(yōu)于其他比較算法。

表5 不同算法MSE和MPSNR比較

3 結(jié)束語

針對稀疏表示未考慮壁畫結(jié)構(gòu)信息與紋理信息的差異性，導(dǎo)致修復(fù)結(jié)果易出現(xiàn)紋理模糊和結(jié)構(gòu)線條斷裂等問題，提出了一種基于多尺度輪廓波分解的群稀疏壁畫修復(fù)算法。利用非下采樣輪廓波變換將待修復(fù)壁畫分解為紋理低頻分量和結(jié)構(gòu)高頻分量。由于低頻分量包含著大量待修復(fù)塊的可匹配塊，對于低頻分量采用改進(jìn)群稀疏算法進(jìn)行修復(fù)，以結(jié)構(gòu)相似性作為塊匹配準(zhǔn)則，同時曲線細(xì)節(jié)特征被提取分離至高頻分量上，因此避免了稀疏算法破壞壁畫輪廓特征導(dǎo)致線條斷裂的現(xiàn)象，并有效克服了稀疏算法修復(fù)結(jié)果中出現(xiàn)的塊效應(yīng)問題。采用三次立方卷積插值算法加強紋理特征，獲得高分辨率的輪廓線條。通過NSCT逆變換將各分量融合重構(gòu)得到最終修復(fù)效果。通過對比，筆者所提方法的修復(fù)結(jié)果在主客觀評價上均優(yōu)于比較算法，更加符合視覺一致性。