(1.河海大學(xué) 土木工程學(xué)院， 南京 210098;2.南京信息工程大學(xué) 遙感學(xué)院， 南京 210044)

摘 要：基于偏微分方程的圖像修復(fù)算法運(yùn)行緩慢，且無(wú)法恢復(fù)紋理細(xì)節(jié)，實(shí)用性較差。基于地統(tǒng)計(jì)學(xué)思想，提出一種簡(jiǎn)單有效的基于各向異性插值模型的圖像修復(fù)方法。實(shí)驗(yàn)表明該方法具有計(jì)算復(fù)雜度低和能夠恢復(fù)圖像紋理細(xì)節(jié)的優(yōu)點(diǎn)，對(duì)于圖像小區(qū)域劃痕具有很好的修復(fù)效果，具有較高的實(shí)用價(jià)值。

關(guān)鍵詞：圖像修復(fù)；插值；各向異性

中圖分類(lèi)號(hào)：TP391文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-3695(2009)04-1554-03

Fast image inpainting algorithm based on anisotropic interpolation model

CHEN Ren-xi1，LI Xin-hui2

(1.College of Civil Engineering， Hohai University， Nanjing 210098， China;2.School of Remote Sensing， Nanjing University of Information Science Technology， Nanjing 210044， China)Abstract:Partial differential equation based image inpainting algorithms have disadvantage of high time consuming and are not practical. Based on the thinking of spatial statistic，this paper presented an inpainting algorithm using an anisotropic interpolation model. Experiment results show that the algorithm has the advantages of low computation cost and texture preserved ability.The algorithm is suitable for restoring missing information in small regions such as scribed gaps and is practicable in applications.

Key words:image inpainting;interpolation;anisotropic

0 引言

數(shù)字圖像修復(fù)是目前新穎而前沿的研究熱點(diǎn)，在照片潤(rùn)飾、視頻文字去除、圖像目標(biāo)隱藏、圖像壓縮、印刷出版等方面具有廣闊的應(yīng)用前景。圖像修復(fù)的目的是恢復(fù)圖像中缺失的信息。根據(jù)待處理區(qū)域的大小，圖像修復(fù)總體上分成兩大類(lèi)：基于圖像幾何模型的修復(fù)方法(inpainting)，主要適用于劃痕等小區(qū)域的處理；基于紋理的圖像補(bǔ)全技術(shù)(completion)，適合于大面積缺失區(qū)域的恢復(fù)[1]。這兩類(lèi)方法各有優(yōu)缺點(diǎn)，適合于不同的情況。其中，inpainting最早由Bertalmio等人[2]提出，他們采用基于偏微分方程(partial differential equation，PDE)的算法，利用擴(kuò)散方程，將待修復(fù)區(qū)域外圍的信息沿著等照度線(xiàn)(isophote)方向向內(nèi)傳播，達(dá)到恢復(fù)丟失信息的目的；Chan等人[3]提出基于總變分(total variation)的修復(fù)方法;張紅英等人[4]提出p-Laplace各向異性擴(kuò)散算法。這些基于PDE的方法取得了一定的成果，但在實(shí)際應(yīng)用中還存在一些缺陷：

a)需要大量迭代運(yùn)算，處理速度非常慢，實(shí)用性不夠強(qiáng)；

b)不能恢復(fù)某些紋理細(xì)節(jié)，造成修復(fù)區(qū)域的模糊現(xiàn)象。

為了克服這些缺陷，本文仍針對(duì)圖像劃痕等小面積修復(fù)的問(wèn)題，將圖像看做三維曲面，模擬三維局部曲面重建，力圖找出一種簡(jiǎn)單而有效的解決方法。通過(guò)將地統(tǒng)計(jì)學(xué)(geostatistic)的思想引入其中，提出一種簡(jiǎn)單快速的各向異性插值修復(fù)模型，并進(jìn)行了算法實(shí)驗(yàn)，取得了較好的結(jié)果。

1 修復(fù)模型

1.1 圖像的各向異性分析

在自然界，某些自然現(xiàn)象具有各向異性的空間分布特性，如礦產(chǎn)分布、空氣污染、地貌等。對(duì)于圖像的灰度值，通過(guò)觀(guān)察發(fā)現(xiàn)也具有相似的分布特性，特別在圖像邊緣或方向性紋理區(qū)域，灰度值并不是在各個(gè)方向上具有相同的分布規(guī)律，而具有各向異性。在地學(xué)領(lǐng)域，分析這種各向異性現(xiàn)象的強(qiáng)大工具是地統(tǒng)計(jì)學(xué)。地統(tǒng)計(jì)學(xué)通過(guò)變異函數(shù)來(lái)定量地描述空間變異性，變異函數(shù)中三個(gè)重要的參數(shù)是變程、基臺(tái)值和塊金值［5］。

將地統(tǒng)計(jì)學(xué)工具應(yīng)用于圖像分析，通過(guò)對(duì)圖像局部區(qū)域計(jì)算各個(gè)方向的變異函數(shù)曲線(xiàn)發(fā)現(xiàn)：沿著紋理方向和邊緣方向，變異函數(shù)曲線(xiàn)具有較小的變程和基臺(tái)值，根據(jù)地統(tǒng)計(jì)學(xué)原理，說(shuō)明該方向上像素的變異程度最小，像素間的相關(guān)性最大。這與對(duì)圖像的觀(guān)察結(jié)果一致，為各向異性插值模型提供了理論支持。

1.2 插值模型

通過(guò)以上分析，認(rèn)為在圖像局部區(qū)域，灰度值具有各向異性分布特征，據(jù)此提出一種各向異性插值模型（圖1）。設(shè)圖像正常區(qū)域記為Φ，缺損區(qū)域記為Ω，缺損區(qū)域的邊界記為Ω。假設(shè)當(dāng)前需要修補(bǔ)的像素位于邊界Ω上的p點(diǎn)，取p點(diǎn)的一個(gè)橢圓鄰域Na，b(p)。其中：a表示橢圓的長(zhǎng)半軸；b表示橢圓的短半軸。橢圓的長(zhǎng)半軸方向與該局部區(qū)域的紋理走向(或邊緣方向)相同。本文認(rèn)為p的像素灰度值取決于橢圓鄰域Na，b(p)中的像素值，理論上可以采用基于地統(tǒng)計(jì)學(xué)原理的克里金(Kriging)插值方法來(lái)估計(jì)p的灰度值。

本文采用橢圓鄰域窗口來(lái)取代傳統(tǒng)的圓形鄰域窗口，并且橢圓長(zhǎng)軸方向自適應(yīng)性地保持與圖像局部邊緣走向一致。沿著紋理方向，像素具有更大的相關(guān)性；而垂直于紋理方向，像素具有較小的相關(guān)性。因此，橢圓鄰域窗口充分考慮到數(shù)據(jù)空間分布的特點(diǎn)，顧及到圖像的邊緣和紋理方向，具有更大的合理性。

基于地統(tǒng)計(jì)學(xué)原理的克里金插值本質(zhì)上也是一種加權(quán)平均插值法，不過(guò)它能夠充分考慮數(shù)據(jù)空間場(chǎng)的分布規(guī)律，在插值過(guò)程中能夠反映空間場(chǎng)的各向異性，并能充分利用數(shù)據(jù)點(diǎn)之間的相關(guān)性。若直接采用克里金插值來(lái)修復(fù)圖像，則計(jì)算量比較大。首先需要統(tǒng)計(jì)局部區(qū)域的半變異函數(shù)值，然后選擇合適的半變異函數(shù)模型進(jìn)行參數(shù)擬合，涉及到比較復(fù)雜的計(jì)算過(guò)程。另外，對(duì)每個(gè)未知點(diǎn)進(jìn)行插值時(shí)，需要進(jìn)行線(xiàn)性方程組的求解，當(dāng)插值的未知像素比較多時(shí)，需要大量的計(jì)算時(shí)間。為了簡(jiǎn)化插值過(guò)程，基于地統(tǒng)計(jì)學(xué)的思想，采用一種簡(jiǎn)易的各向異性插值模型(圖2)。

假設(shè)點(diǎn)p(x0，y0)是當(dāng)前需要插值的點(diǎn)，且該點(diǎn)的方向矢量(即邊緣走向)E=(Ex，Ey)。以p(x0，y0)為中心點(diǎn)，確定一個(gè)矩形窗口(橢圓窗口的簡(jiǎn)化)，窗口的長(zhǎng)軸方向與方向矢量E相同，矩形大小為a×b(且a>b)。采用矩形窗口內(nèi)部的已知像素點(diǎn)對(duì)未知點(diǎn)p(x0，y0)進(jìn)行插值，其灰度值G(x0，y0)由窗口內(nèi)已知像素的灰度值加權(quán)平均來(lái)估計(jì)：

G(x0，y0)=Ni=1wiG(xi，yi)/Ni=1wi(1)

其中：N為窗口內(nèi)已知像素的總數(shù)；權(quán)值wi的計(jì)算主要考慮已知點(diǎn)與p(x0，y0)之間的距離和方向兩個(gè)因素。另設(shè)p(x1，y1)為窗口內(nèi)的某個(gè)已知點(diǎn)，p(x0，y0)與p(x1，y1)間的歐式距離為

dist=(x0-x1)2+(y0-y1)2(2)

點(diǎn)p(x0，y0)與p(x1，y1)之間的向量為

V=(Vx，Vy)=(x1-x0，y1-y0)(3)

向量V與方向矢量E之間的夾角記為θ：

cos(θ)=‖V#8226;E‖/‖V‖#8226;‖E‖(4)

則點(diǎn)p(x1，y1)的權(quán)值計(jì)算如下：

w1=1/(1+dist）(1+|cos(θ)|)(5)

從式(5)可知，這種權(quán)值計(jì)算方法同時(shí)考慮到數(shù)據(jù)點(diǎn)的空間距離和分布方向兩個(gè)因素。若該點(diǎn)距離待插點(diǎn)越遠(yuǎn)，則貢獻(xiàn)越小，權(quán)值較??；若該點(diǎn)的方向越靠近主軸方向(方向矢量的方向)，則權(quán)值越大。根據(jù)前面的空間統(tǒng)計(jì)分析，位于主軸方向的像素與待插點(diǎn)具有更大的相似性，賦予較大的權(quán)值是合理的。

根據(jù)實(shí)驗(yàn)表明，矩形窗口的長(zhǎng)軸a和短軸b一般不能太大，否則影響插值效果。一般情況下，取a=5~7像素，b=3~5像素時(shí)能夠較好地保持紋理細(xì)節(jié)，若再增大窗口的尺寸會(huì)帶來(lái)一些模糊效應(yīng)。

2 算法實(shí)現(xiàn)

2.1 預(yù)處理

本文算法需要計(jì)算圖像邊緣方向場(chǎng)，宏觀(guān)反映圖像邊緣和紋理走向。然而圖像噪聲往往影響計(jì)算的穩(wěn)定性，為此，先對(duì)圖像進(jìn)行平滑去噪。通常的平滑算子易造成圖像邊緣模糊現(xiàn)象，考慮采用能夠保持邊緣的平滑方法，如反梯度平均的卷積方法[6]。在計(jì)算卷積模板時(shí)，每個(gè)像素處的權(quán)值根據(jù)反梯度值來(lái)計(jì)算，其思想是區(qū)域內(nèi)部的亮度變化一般比相鄰區(qū)域間要小。設(shè)卷積模板(模板一般為奇數(shù)大小，如3×3、5×5)中心像素坐標(biāo)為(p， q)，則模板中的點(diǎn)(i， j)相對(duì)于中心點(diǎn)(p， q)的反梯度值為：

δ(i， j)=1/|g(p，q)-g(i， j)|當(dāng)g(p，q)≠g(i， j)

2當(dāng)g(p，q)=g(i， j)(6)

對(duì)于N×N大小的卷積模板h，需要進(jìn)行歸一化處理，其中每個(gè)點(diǎn)(中心位置除外)的數(shù)值計(jì)算如下：

h(i， j)=0.5×δ(i， j)/Ni=1

Nj=1δ(i， j)(7)

另外，模板中心位置的數(shù)值取為0.5。當(dāng)采用模板h對(duì)圖像進(jìn)行平滑時(shí)，區(qū)域中孤立噪聲具有小的權(quán)值，采用鄰域點(diǎn)平均時(shí)噪聲得到消除，且能夠保持圖像邊緣不被模糊。

2.2 方向場(chǎng)的計(jì)算

把圖像中每個(gè)像素所在的邊緣走向(即邊緣的切線(xiàn)方向)定義為方向矢量，則全圖形成一個(gè)方向矢量場(chǎng)，反映圖像紋理和邊緣的方向。切線(xiàn)方向往往可以根據(jù)梯度(gradient)來(lái)計(jì)算，因?yàn)樘荻扰c切線(xiàn)的方向具有互相垂直的關(guān)系。

對(duì)于每個(gè)像素點(diǎn)(i， j)，采用如圖3所示的梯度算子來(lái)計(jì)算該點(diǎn)的梯度值G(i， j)，然后將G(i， j)旋轉(zhuǎn)90°或者-90°，作為邊緣矢量E(i， j)。

采用上述方法直接計(jì)算得到的圖像方向場(chǎng)，結(jié)果不太理想，矢量方向不穩(wěn)定，存在局部擺動(dòng) (圖4(a))。實(shí)際上，某點(diǎn)的方向矢量應(yīng)該由該點(diǎn)周?chē)徲虻乃邢袼氐姆较蚴噶抗餐瑳Q定，因此，采用鄰域加權(quán)平均的方法來(lái)重新計(jì)算每個(gè)像素點(diǎn)的方向矢量。假設(shè)一個(gè)以點(diǎn)(p， q)為中心的(2n+1)×(2n+1)大小的一個(gè)鄰域窗口，對(duì)于窗口中的每個(gè)像素位置(i， j)，采用“距離+1”的倒數(shù)作為該點(diǎn)的權(quán)值：

w(i， j)=1/[1+(p-i)2+(q-j)2](8)

實(shí)驗(yàn)表明，采用3×3、5×5大小的鄰域窗口已經(jīng)足夠，若增大窗口，效果不再改變，反而增加計(jì)算量。圖4（b）是采用3×3窗口對(duì)方向場(chǎng)進(jìn)行平滑后的結(jié)果。

2.3 插值修復(fù)

完成前面的步驟后，最后進(jìn)行圖像像素的修復(fù)，分為兩步進(jìn)行：

a）進(jìn)行方向場(chǎng)的恢復(fù)，以便后面的各向異性灰度插值。方向場(chǎng)插值可以直接采用加權(quán)平均插值方法：以待插點(diǎn)為中心選擇7×7大小的窗口，利用窗口內(nèi)已知點(diǎn)上的方向矢量來(lái)估計(jì)待插點(diǎn)處的方向矢量，權(quán)值可采用距離倒數(shù)法。插值時(shí)從缺損區(qū)域邊界逐層向內(nèi)進(jìn)行。圖5是對(duì)缺失區(qū)域方向場(chǎng)重建的示例。

b)進(jìn)行灰度值的恢復(fù)，利用前述的各向異性插值模型完成缺失區(qū)域的灰度值重建。處理順序與方向場(chǎng)插值一樣，由邊界向內(nèi)逐層進(jìn)行。

3 實(shí)驗(yàn)

根據(jù)前面提出的各向異性圖像修復(fù)模型，采用C++實(shí)現(xiàn)了本文的算法。為了與其他的圖像修復(fù)算法進(jìn)行比較，同時(shí)也實(shí)現(xiàn)了Bertalmio[2]和Oliveira算法[7]。另外，為了比較基于地統(tǒng)計(jì)學(xué)思想的各向異性插值方法和一般的插值方法的優(yōu)劣，同時(shí)也采用了反距離加權(quán)平均(inverse distance weight，IDW)的插值方法對(duì)圖像進(jìn)行修復(fù)。本文采用不同的方法對(duì)不同類(lèi)別的圖像進(jìn)行了大量的修復(fù)實(shí)驗(yàn)。實(shí)驗(yàn)中，在進(jìn)行方向場(chǎng)計(jì)算時(shí)，采用5×5的鄰域窗口；在進(jìn)行灰度值各向異性插值時(shí)，矩形窗口的大小選為a=7，b=3。

圖6是采用不同的方法對(duì)遙感影像上紋理區(qū)域進(jìn)行修復(fù)的結(jié)果，其劃痕寬度在7~8個(gè)像素。顯而易見(jiàn)，采用各向異性插值方法能夠恢復(fù)精細(xì)的紋理細(xì)節(jié)，而其他方法則存在不同程度的模糊現(xiàn)象。

圖7是對(duì)受損的遙感影像的修復(fù)實(shí)驗(yàn)，劃痕寬度在8~12個(gè)像素。從實(shí)驗(yàn)結(jié)果可以看出，Bertalmio、Oliveira修復(fù)方法對(duì)于邊緣的修復(fù)效果不太好，存在比較嚴(yán)重的模糊現(xiàn)象；IDW插值方法也不能很好地恢復(fù)斷裂邊緣；本文的方法對(duì)斷裂邊緣的修補(bǔ)效果較好。

表1統(tǒng)計(jì)了不同修復(fù)方法的運(yùn)行時(shí)間。從表中可以看出，Bertalmio算法非常耗時(shí)，遠(yuǎn)遠(yuǎn)超過(guò)其他方法的運(yùn)行時(shí)間，因?yàn)樵撍惴ㄐ枰M(jìn)行反復(fù)的迭代(本實(shí)驗(yàn)中迭代1 000次)。本文的方法運(yùn)行速度很快，與Oliveira快速圖像修復(fù)方法和IDW插值方法相當(dāng)。

表1 圖像修復(fù)時(shí)間統(tǒng)計(jì)

圖像Bertalmio/minOliveira/sIDW/s本文方法/s

圖67233

圖72122

為了進(jìn)一步定量地評(píng)價(jià)該算法的效果，采用三種不同類(lèi)型的彩色圖像（自然風(fēng)景圖像rock.bmp中國(guó)畫(huà)paint.bmp、遙感圖像RS.bmp）進(jìn)行模擬實(shí)驗(yàn)。先對(duì)正常圖像添加人工劃痕，模擬圖像受損，然后采用不同的修復(fù)方法對(duì)缺損區(qū)域進(jìn)行修復(fù)，再比較修復(fù)后的圖像與原始正常圖像之間的誤差。這里采用均方誤差(MSE)來(lái)進(jìn)行衡量，對(duì)于彩色圖像分R、G、B三個(gè)通道分別計(jì)算。 表2列出了相應(yīng)的MSE統(tǒng)計(jì)數(shù)據(jù)(修復(fù)結(jié)果圖略)，從表中可見(jiàn)本文方法的MSE最小。

表2 圖像修復(fù)均方差(MSE)統(tǒng)計(jì)結(jié)果

圖像通道BertalmioOliveiraIDW本文方法

rock.bmpR295.92360.23406.34250.79

G291.05340.14375.89241.90

B288.57300.92332.27229.68

paint.bmp

R650.56526.85572.95437.34

G490.16474.50492.52387.20

B402.32437.26438.15346.08

RS.bmp

R2 118.201 067.841 197.05802.42

G1 981.72943.871 189.40816.01

B1 619.74734.78954.24614.96

通過(guò)以上實(shí)驗(yàn)表明，從主觀(guān)視覺(jué)效果來(lái)看，本文方法能夠較好地連接圖像中的斷裂邊緣，恢復(fù)紋理結(jié)構(gòu)；從定量統(tǒng)計(jì)結(jié)果來(lái)看，本文方法對(duì)于圖像缺損區(qū)域重建的誤差較小。

4 結(jié)束語(yǔ)

本文針對(duì)圖像上小區(qū)域缺失信息的修復(fù)問(wèn)題，提出一種簡(jiǎn)單高效的處理算法。該方法基于圖像像素的各向異性分布特征，采用一種各向異性插值模型，充分顧及到像素的空間分布規(guī)律，能夠較好地恢復(fù)圖像邊緣和紋理。該方法克服了PDE修復(fù)方法計(jì)算量大的缺陷，也避免了傳統(tǒng)IDW插值方法和Oliveira快速修復(fù)方法帶來(lái)的模糊問(wèn)題，不失為一種簡(jiǎn)單高效的解決方案。

參考文獻(xiàn)：

［1］張紅英，彭啟琮.數(shù)字圖像修復(fù)技術(shù)綜述[J].中國(guó)圖象圖形學(xué)報(bào)，2007，12(1):1-10.

[2]BERTALMIO M，SAPIRO G，CASELLES V，et a1.Image inpainting[C]//Proc of the ACM SIGGRAPH Conference on Computer Gra-phics.New York：ACM Press，2000：417-424.

[3]CHAN T F，SHEN J H.Mathematical models for local nontexture inpainting[J].SIAM Journal of Applied Mathematics，2001，62(3):1019-1043.

[4]張紅英，彭啟琮，吳亞?wèn)|.數(shù)字破損圖像的非線(xiàn)性各向異性擴(kuò)散修補(bǔ)算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2001，18(10):1541-1546.

［5］侯景儒，黃競(jìng)先.地質(zhì)統(tǒng)計(jì)學(xué)及其在礦產(chǎn)儲(chǔ)量計(jì)算中的應(yīng)用[M].北京:地質(zhì)出版社，1982.

[6]SONKA M，HLAVAC V，BOYLE R.圖像處理、分析與機(jī)器視覺(jué)[M].艾海舟，武勃，等譯.北京: 人民郵電出版社， 2003.

[7]OLIVEIRA M M，BOWEN B，MCKENNA R，et al.Fast digital image inpainting[C]//Proc of International Conference on Visualization， Imaging， and Image Processing.2001:261-266.