基于灰色系統理論的圖像修復研究

曾衛波，邢永康，石 楊

(重慶大學計算機學院，重慶 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.