塊核范數的RPCA分解與熵權類稀疏的壁畫修復

2021-08-10 12:59:34陶美風

哈爾濱工業大學學報 2021年8期

陳永，陶美風，陳錦

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

敦煌莫高窟是世界上現存規模最宏大、內容最豐富的佛教石窟壁畫寶庫，其內所存的壁畫、經卷等具有珍貴的研究價值。然而，由于自然風化的破壞及人為因素的影響，窟內壁畫出現了地仗脫落、劃痕、褪色、裂紋等嚴重的災害，亟待保護[1]。因此，研究病態敦煌壁畫的修復極其重要。但人工修復存在風險大、耗時長且修復不可逆等問題，將數字化虛擬修復應用于古代壁畫的保護，是近年來國內外文化遺產保護的研究熱點[2]。

目前，傳統圖像修復方法主要分為3類：基于偏微分方程、基于樣本塊和基于稀疏表示的修復方法。基于偏微分方程的主要修復模型有TV(total variation)模型[3]和CDD(curvature-driven diffusion)模型[4]，該類方法可擴散完成小區域破損圖像的修復[5]?；跇颖緣K的修復方法以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]。稀疏表示的關鍵在于尋找給定信號的最佳稀疏解，其數學定義為

(1)

式中：α為稀疏系數，T為待觀測信號，D為字典，‖·‖0為l0范數，用來衡量α中非零元素個數。一般自然信號存在噪聲或重構圖像會出現一定的誤差ε，因此，為了平衡圖像的稀疏性和重構誤差，式(1)可以改寫為

(2)

在稀疏表示的圖像修復中，圖像退化的過程可定義為

Y=JX+n

(3)

式中：Y為退化后的圖像，X為未退化的原始圖像，J為退化算子，n為高斯白噪聲。因此，圖像修復過程就是對退化圖像Y進行推測估計得到原始圖像X的過程。

根據式(2)，在稀疏字典D給定的條件下，可得到稀疏表示圖像修復的正則化模型為

(4)

通過式(4)求解稀疏系數α，最后通過Y=Dα得到修復后的圖像。

2 本文算法

2.1 基于塊核范數的RPCA圖像分解

一般圖像中既包含結構信息，也包含紋理信息，其中結構信息主要是圖像邊緣和輪廓等變化劇烈的信息部分；紋理信息主要是圖像中亮度或灰度值變化緩慢或周期性變化區域，具有重復性、規則性和方向性，其反映了圖像中同質現象的一種視覺特征[15,18]。RPCA[19-20]是一種高維數據降維分析方法，通過RPCA可以將圖像M(M∈Rm×n)分解為兩個矩陣的線性相加，一個是結構矩陣L(L∈Rm×n，rank(L)≤m,n)，主要用來表征圖像的結構信息；另一個是紋理矩陣S(S∈Rm×n)，主要是圖像的紋理細節信息，其圖像分解模型定義為

(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

塊核范數‖H‖BNN是紋理部分的局部圖像塊，經過方向矯正和重疊取塊操作后形成的矩陣奇異值之和，其定義為

(6)

圖3 塊核范數示意Fig.3 Schematic of block nuclear norm

將塊核范數引入到RPCA圖像分解模型，則得到如下的分解模型

(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)為

(8)

式中：PA(i)為圖像塊A中RGB顏色通道中第i個等級的像素占比率；i表示像素等級，范圍為0～255。定義圖像塊A與聚類中心B的信息熵特征差異性系數為

qj=|E(A)-E(B)|

(9)

式中0≤qj≤1，qj越小，圖像塊A與B越相似。然后，采用式(10)得到圖像塊A的熵權重vj

(10)

從而得到圖像塊A與B的熵權歐氏距離計算準則

(11)

通過式(11)將待分類圖像塊A劃分到熵權歐氏距離最小的類中，從而構造得到結構層的子類圖像。為了驗證熵加權k-means聚類算法對壁畫聚類的有效性，以圖4(b)結構圖像為例進行實驗，結果見圖5。其中，圖5(a)為原始圖像，圖5(b)為聚類結果，圖5(c)、(d)、(e)為各子類圖像，可以發現，具有相似結構的圖像塊能夠較好地進行聚類，得到的結構層圖像聚類更為準確。

2.3 類稀疏表示的結構層圖像修復

在結構層子類圖像劃分的基礎上，對各子類進行類稀疏修復，定義類稀疏修復模型為

(12)

式中：yk是第k類結構層圖像塊；Dk和αk分別是對應的字典和稀疏系數。然后，采用奇異值分解和分裂Bregman迭代優化算法更新求解字典和系數[9]。將圖像塊yk的估計值rk進行奇異值分解，即

(13)

式中：Uk、Vk分別為yk的左奇異正交矩陣和右奇異正交矩陣,∑k為對角矩陣，則字典Dk中每個原子dki為

dki=ukivki,i=1,2,…,t

(14)

式中uki、vki分別為Uk、Vk的列向量，t為第k類字典的原子數。所以，第k類自適應字典Dk為

Dk={dk1,dk2,dk3,…,dki},i=1,2,…,t

(15)

以圖5結構層圖像聚類結果為例，采用奇異值分解學習更新得到自適應子類字典，見圖6。其中，圖6(a)、(b)、(c)分別為圖5子類對應的自適應字典?？梢园l現，本文提出的類稀疏表示方法得到的各子類字典更能夠體現圖像的結構特征，且聚類字典訓練更加高效。

圖5 熵加權的k-means算法聚類結果Fig.5 Clustering results of entropy weighted k-means algorithm

圖6 自適應字典Fig.6 Adaptive dictionaries

在完成子類字典的更新后，接著對子類圖像進行類稀疏系數求解

(16)

采用分裂Bregman迭代優化算法對式(16)求解，從而完成破損圖像塊的稀疏編碼。最后，利用各子類圖像的字典和系數完成結構層各子類圖像的修復，并按照索引矩陣放回原始結構層圖像，即得到結構層圖像的修復結果YS

(17)

2.4 分層圖像融合重構

采用雙三次插值(bicubic interpolation)對紋理層破損區域進行插值修復，完成圖像紋理細節信息的重建。雙三次線性插值是通過對待修復點最近的16個采樣點的加權平均完成像素點的估計，按下式計算

(18)

式中：f(p,q)為待修復點坐標，W(p-pi)和W(q-qj)分別為待修復點的橫坐標權重和縱坐標權重。權重的計算如下

(19)

最后將結構層的修復結果YS與紋理層的修復結果YT融合，得到最終的修復圖像Y

Y=YS+YT

(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的修復結果做出定量評價，采用PSNR和SSIM進行比較，見表1。PSNR和SSIM是衡量圖像質量的重要定量評價指標，其中PSNR值越大，代表失真越少，即修復效果越好；SSIM是一種衡量兩幅圖像相似度的指標，該值越大，表明修復后圖像的失真程度越小[18]。從表1中可以發現，本文算法在PSNR和SSIM上均優于其他對比算法，表明本文方法修復質量更好，從而驗證了本文方法對于人工破損壁畫修復的有效性。

圖7 不同算法對人為添加破損壁畫的修復結果對比Fig.7 Comparison of mural inpainting results of artificial damages by different algorithms

表1 不同算法修復結果PSNR和SSIM對比Tab.1 Comparison of PSNR and SSIM results of different algorithms

表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)為本文結果，可以看出，在線條連續性和清晰性方面本文修復效果均有所提高。