基于Navier-Stokes方程的圖像修復(fù)算法*

李率杰,李 鵬,馮兆永，姚正安

(中山大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院，廣東 廣州 510275)

圖像修復(fù)技術(shù)是指針對圖像中遺失或者損壞的部分,利用周圍未被損壞的圖像信息,按照一定的規(guī)則進(jìn)行填補(bǔ),使修復(fù)后的圖像接近或達(dá)到原圖的視覺效果，在本質(zhì)上是一種典型的插值技術(shù)。圖像修復(fù)具有廣泛的應(yīng)用前景，如填補(bǔ)美術(shù)作品上出現(xiàn)的裂痕，以使其清晰并恢復(fù)它的完整性；恢復(fù)舊照片中的劃痕；去除圖像中的文字或物體等特效制作；視頻通信中錯誤隱匿；圖像遮擋物體的去除等。

近來的數(shù)字圖像修復(fù)技術(shù)大多是利用偏微分方程(partial differential equation，PDE)和變分法。而這種基于PDE的數(shù)字圖像修復(fù)技術(shù)最早是由Bertalmio等[1]首先引入到圖像處理領(lǐng)域中的。在[1] 中，他們利用待修復(fù)區(qū)域的周圍信息，采用一種由粗到精的方法來估計等照度線的方向，并采用傳播機(jī)制將信息傳播到待修復(fù)區(qū)域內(nèi)。之后Balleater等[2]通過提出一個能量泛函并利用變分法極小化該能量得到非線性的偏微分方程計算出修復(fù)后的圖像。在文獻(xiàn)[3-4]中作者考慮圖像的等照度線方向和不可壓Navier-Stokes方程之間的聯(lián)系將圖像灰度作為二維不可壓流體的流函數(shù)，將Laplacian算子作用后的待修復(fù)圖像作為流的渦流函數(shù)。基于由流函數(shù)定義的向量場將渦流函數(shù)傳播到待修復(fù)區(qū)域，從而得到修復(fù)后的圖像。由此Tai等[5]提出了基于TV-Stokes 方程的圖像修復(fù)兩步算法,首先通過求解一個非線性TV-Stokes方程將等照度線延伸到待修復(fù)區(qū)域,然后圖像就以此等照度線方向進(jìn)行修復(fù)。

不同于以上方法，Chan等[6]引入了Rudin等[7]提出的全變差(total variation，TV)圖像去噪模型的思想，得到基于TV模型的圖像修復(fù)模型，該方法通過求解最小化能量泛函而得到非線性方程，無論待修復(fù)圖像是否帶有噪聲都能夠有效地對圖像進(jìn)行修復(fù)，得到了較好的效果。但是TV模型在圖像平滑區(qū)域容易產(chǎn)生階梯效應(yīng)。在后面的工作如文獻(xiàn)[8-10]中，Chan和Shen提出了基于圖像等照度線的曲率的能量泛函，試圖在保證光滑條件下連接修復(fù)區(qū)域周圍的等照度線，由此得到的方程是非線性的四階方程，該方法可以較好的修復(fù)面積較大的非紋理區(qū)域，但是非線性的高階方程給求解帶來了困難?；谝陨螩han和Shen的工作，Esedoglu等[11]基于Mumford-Shah圖像分割模型提出了Mumford-Shah-Euler圖像修復(fù)模型[12]，利用Ambrosio和Tortorelli提出的Mumford-Shah模型橢圓逼近方法得到相應(yīng)的拋物方程。方法的優(yōu)點(diǎn)是方程的最高階偏導(dǎo)是線性的可以快速的計算。但是卻有難以修復(fù)較大的區(qū)域的缺點(diǎn)。Grossauer等[13]利用復(fù)Ginzburg-Landau方程解決圖像修復(fù)問題，近來Bornemann等[14]提出基于一階傳輸方程的圖像修復(fù)的非迭代方法。但是這些方法計算量大且非常復(fù)雜。

另一方面紋理圖像的修復(fù)也受到廣泛的研究，在文[15]中，作者首先將待修復(fù)區(qū)域分解為紋理和結(jié)構(gòu)，然后對紋理采用紋理合成的方法，對結(jié)構(gòu)采用一般的非紋理圖像修復(fù)模型進(jìn)行分別修復(fù)。另外最近紋理合成的方法也同樣被廣泛的運(yùn)用在紋理圖像修復(fù)中[16-20]。

本文基于文獻(xiàn)[3-5]中的思想，我們在等照度線方向上向待修復(fù)區(qū)域傳輸圖像信息。在文獻(xiàn)[3]通過求解不可壓流體的Navier-Stokes方程，得到修復(fù)后的圖像函數(shù)。假設(shè)D為待修復(fù)區(qū)域，?D為待修復(fù)區(qū)域的邊界。一旦給定待修復(fù)區(qū)域D此模型就可以將信息沿等照度方向連續(xù)的從邊界?D傳輸?shù)酱迯?fù)區(qū)域D。如此圖像修復(fù)依據(jù)連接和待修復(fù)區(qū)域邊界相交的等照度線得以完成。但是該修復(fù)算法只是基于待修復(fù)區(qū)域的邊界?D上的信息，因此難以得到與待修復(fù)區(qū)域相交的等照度線的準(zhǔn)確方向。另外如果待修復(fù)區(qū)域D是規(guī)則的矩形區(qū)域，算法容易實(shí)現(xiàn)，但是實(shí)際上待修復(fù)區(qū)域往往是形狀不規(guī)則的區(qū)域，實(shí)現(xiàn)此算法將面臨執(zhí)行有限差分格式和邊界條件的困難，最后該修復(fù)算法對噪聲圖像的修復(fù)無能為力。

與上面算法不同的是，本文修復(fù)算法的思想是能夠在待修復(fù)區(qū)域進(jìn)行信息填充，而在待修復(fù)區(qū)域外異性光滑圖像去除噪聲(如果存在)，且保持邊緣信息。在本文中，基于不可壓流的Navier-Stokes方程，我們主要提出一種新的圖像修復(fù)算法，目的是任給一個帶有噪聲的且信息丟失或被破壞的待修復(fù)圖像，我們可以同時的進(jìn)行圖像的修復(fù)和去噪，得到清晰、完整的圖像。

1 基于Navier-Stokes方程的圖像修復(fù)

不可壓的流體一般遵循如下所謂的Navier-Stokes方程

υt+υ·▽υ=-▽p+νΔυ,▽·υ=0

(1)

其中υ為流體的速度場，p為流體的壓強(qiáng)，ν為黏性系數(shù)。對于二維空間，我們引進(jìn)流函數(shù)Ψ，滿足

▽⊥Ψ=υ

(2)

同時設(shè)渦量ω=▽×υ，且由(1)可得滿足渦流方程

ωt+υ·▽ω=νΔω

(3)

在二維空間里，渦量ω是一個標(biāo)量，并且滿足ω=ΔΨ。當(dāng)黏性系數(shù)ν=0時，可以得到無黏性流的歐拉方程。在二維時黏性流或無黏性流，方程都是適定的，對于任意光滑初始條件下解都是存在的，并且連續(xù)的依賴初始條件和邊界條件[21]。

根據(jù)流函數(shù)的定義，黏性流方程(3)的穩(wěn)定態(tài)一定滿足下列條件

▽⊥Ψ·▽ΔΨ≈0

(4)

此式說明Laplacian算子作用后流函數(shù)ΔΨ，其等照度線的方向和流體的流速方向是平行的。即流體是沿著ΔΨ等照度線方向流動，如此其和在圖像修復(fù)中將周圍信息沿等照度線方向傳播到待修復(fù)區(qū)域的原理相同。

將二維不可壓流體力學(xué)中的流函數(shù)Ψ和圖像修復(fù)中待修復(fù)圖像的灰度函數(shù)I0對應(yīng)起來，Bertalmi等[3]提出了基于Navier-Stokes方程的圖像修復(fù)模型。設(shè)Ω為整個圖像區(qū)域，D為待修復(fù)區(qū)域，?D為待修復(fù)區(qū)域的邊界，待修復(fù)圖像I0在ΩD是足夠光滑，并且I0和ΔI0在ΩD是已知的。用圖像灰度函數(shù)I對應(yīng)流函數(shù)Ψ，將二維不可壓流體理論和圖像修復(fù)理論對應(yīng)關(guān)系總結(jié)如下：

Bertalmio等[3]目的是在待修復(fù)區(qū)域上解決Navier-Stokes方程的一種形式，即渦流方程(3)。他們解決如下方程：

表1 流體力學(xué)和圖像修復(fù)物理量對應(yīng)關(guān)系

ωt+υ·▽ω=ν▽·(g(|▽ω|)▽ω)

(5)

其中渦量ω=ΔI，函數(shù)g為單調(diào)遞減函數(shù)，可以使函數(shù)ω異性擴(kuò)散，速度場υ=▽⊥I。由求解方程(5)所得的渦量ω,通過求解Poisson方程：

ΔI=ω，I?D=I0

(6)

從而得到修復(fù)后的圖像I。

2 圖像修復(fù)算法

本文中，我們的目的是在整個區(qū)域上求解Navier-Stokes方程的一種不同于(5)式的一種形式，可以使得在修復(fù)區(qū)域D內(nèi)進(jìn)行圖像修復(fù)，而在修復(fù)區(qū)域外ΩD進(jìn)行圖像去噪。我們求解下列渦流方程：

(7)

(8)

其中χD為修復(fù)區(qū)域D的特征函數(shù)，定義為χD=1,x∈D，χD=0,x∈ΩD。由方程求得的ω，我們同樣可以在區(qū)域Ω求與(6)式類似的Poisson方程，得到修復(fù)后的圖像I。

由提出的模型,可以看出它能夠在修復(fù)區(qū)域的內(nèi)部和外部起著不同的作用。其修復(fù)過程包括：修復(fù)區(qū)域內(nèi)，基于周圍信息對破壞或丟失的信息進(jìn)行填充；而在修復(fù)區(qū)域外部，對ω進(jìn)行異性光滑，從而可以達(dá)到去噪目的。

3 數(shù)值計算

對方程數(shù)值計算，我們假設(shè)m×n的圖像I是在連續(xù)區(qū)域[0,m]×[0,n]上，離散為m×n網(wǎng)格，定義時間步長Δt。我們用有限差分半隱格式逼近此問題。

Δ

Δtν▽·(g(|▽ωk|)▽ωk+1)i,j+

Δtλ(1-χD)(ω0-ωk+1)i,j

(9)

▽·(g(|▽ωk|)▽ωk+1)i,j=

(10)

(11)

其中K是取定的常數(shù)。

對于Poisson方程ΔIk+1=ωk+1，我們采用五點(diǎn)中心離散格式，將Ik+1和ωk+1重組成mn的向量，得到線性方程組AIk+1=ωk+1，這里A為一個稀疏的mn×mn的矩陣，求解這個線性方程組由ωk+1我們快速的得到修復(fù)后的圖像Ik+1。詳細(xì)的計算流程如下：

3)求解Poisson方程得到Ik；

其中k=1,2,…，直到穩(wěn)定狀態(tài)，得到的Ik既是修復(fù)后的圖像。

4 仿真實(shí)驗(yàn)

在本節(jié)將用數(shù)值實(shí)驗(yàn)結(jié)果展示本模型的有效性和實(shí)用性。我們選了不同類型的圖像，他們有若干區(qū)域丟失信息，劃痕，或者同時帶有噪聲，我們對其進(jìn)行修復(fù)。

圖1 圖像修復(fù)結(jié)果比較

實(shí)驗(yàn)2 本文算法針對兩幅合成圖像修復(fù)結(jié)果。如圖2所示， (a)和(d)分別為原始圖像，(b)和(e)為丟失信息的圖像，為了清晰顯示修復(fù)區(qū)域，我們同樣進(jìn)行了隨機(jī)像素填充。通過修復(fù)，我們可以分別得到它們修復(fù)后的圖像(c)和(f)，從結(jié)果可以看出修復(fù)效果非常明顯。由于二值合成圖像灰度對比較強(qiáng)，在邊緣處和原始圖像有所差別，但是圖像信息基本得到了恢復(fù)。本實(shí)驗(yàn)中參數(shù)選擇如表2所示。

表2 實(shí)驗(yàn)中參數(shù)選擇

圖2 圖像修復(fù)結(jié)果

實(shí)驗(yàn)3 本文算法針對兩幅真實(shí)圖像的修復(fù)結(jié)果，所用參數(shù)和迭代次數(shù)見表2。在圖3中，(a)為一幅原始圖像，(b)為其帶有劃痕的圖像，(c)為修復(fù)后的圖像。從效果上來看修復(fù)圖像和原始圖像非常相近，可以證明本文算法對這種劃痕是有效的。(d)為選擇的另一幅原始圖像，(e)為在其不同位置有信息缺失的圖像，其缺失區(qū)域大多分布在圖像邊緣上，(f)為修復(fù)后的圖像。我們可以看出修復(fù)結(jié)果仍然十分理想。

圖3 圖像修復(fù)結(jié)果

實(shí)驗(yàn)4 本文算法同時去除噪聲和修復(fù)圖像信息的仿真實(shí)驗(yàn)。本實(shí)驗(yàn)利用本章所提的修復(fù)算法同時去除噪聲和修補(bǔ)圖像中缺失的信息，結(jié)果如圖4所示。我們同樣選擇了兩幅圖像，其中(a)和(d)為原始圖像，(b)和(e)為破損且加噪圖像，(c)和(f)為本文算法修復(fù)的圖像。從圖中可以看出，本文算法可以有效地修復(fù)破損的噪聲圖像。

圖4 圖像修復(fù)結(jié)果

實(shí)驗(yàn)5 本文算法移除圖中不需要的景物，修復(fù)圖畫的仿真實(shí)驗(yàn)。結(jié)果如圖5所示，(a)和(d)為原始圖像，(b)和(e)是在圖中標(biāo)出修復(fù)區(qū)域的圖像，(c)和(f)為修復(fù)結(jié)果。從修復(fù)后的圖像可以看出，雖然移除景物面積較大但是本文算法仍然得到了很好的修復(fù)效果，使觀察者難以察覺圖像曾經(jīng)被修改過。但是修復(fù)后的圖像亮度發(fā)生了輕微的變化，相關(guān)計算過程值得我們研究改進(jìn)。

圖5 圖像修復(fù)結(jié)果

5 結(jié) 論

本文基于Navier-Stokes方程，提出了一個新的圖像修復(fù)的數(shù)學(xué)模型。該模型能夠在待修復(fù)區(qū)域?qū)D像丟失或破損的信息進(jìn)行修補(bǔ)，而在待修復(fù)區(qū)域外去除噪聲。

該模型的思想來自于不可壓流體力學(xué)的Navier-Stokes方程的物理背景，以及流體速度方向和圖像的等照度線方向的之間的一致性。對于圖像信息丟失或損壞，并且?guī)в性肼暤膱D像，修復(fù)將被自適應(yīng)方程完成，在修復(fù)區(qū)域內(nèi)表現(xiàn)為信息由邊界向內(nèi)部傳播，在修復(fù)區(qū)域外表現(xiàn)為異性光滑。實(shí)驗(yàn)結(jié)果顯示了模型在處理破損圖像(圖2，圖3和圖5)或噪聲圖像(圖4)時的有效性和實(shí)用性。但是模型也存在某些不足，如在保持修復(fù)區(qū)域圖像邊緣方面如圖2，圖4，如何避免這些缺點(diǎn)將是我們以后繼續(xù)研究的工作。

參考文獻(xiàn)：

[1]BERTALMIO M,SAPIRO G,BALLESTER C,et al.Image inpainting[C]∥Siggraph 2000,Computer Graphics proceedings,K.Akeley,Ed ACM Press/ ACM SIGGRAPH/ Addison Wesley Longman,2000：417-424.

[2]BALLESTER C,BERTALMIO M,CASELLES V.Filling in by joint interpolation of vector fields and gray levels[J].IEEE Trans Image Processing,2001,10(8): 1200-1211.

[3]BERTALMIO M,BERTOZZI A L,SAPIRO G.Navier-Stokes,fluid dynamics and image and video inpainting[J].IEEE Computer Vision and Pattern Recognition (CVPR),2001,1: 355-362.

[4]WILSON A,RYO T.Inpainting with the Navier-Stokes equations [R].http://www.math.ucla.edu/～rrtakei/gradProj/930project.pdf.

[5]TAI X C,OSHER S,HOLM R.Image inpainting using TV-Stokes equation[C]∥Image Processing Based on Partial Differential Equations,Springer,Heidelberg,2006:3-22.

[6]CHAN T F,SHEN J.Mathematical models for local non-texture inpaintings[J].SIAM J Appl Math,2002,62(3): 1019-1043.

[7]RUDIN L,OSHER S，F(xiàn)ATEMI E.Nonlinear total variation based noise removal algorithms[J].Physica D,1992,60:259-268.

[8]CHAN T F,KANG S H,SHEN J.Euler’s elastica and curvature-based inpainting[J].SIAM J Appl Math,2002,63(2): 564-592.

[9]CHAN T F,SHEN J,VESE L.Variational PDE models in image processing[J].Notices Am Math Soc,2003,50(1): 14-26.

[10]CHAN T,SHEN J.Non-texture inpainting by curvature-driven diffusions (CDD)[J].J Visual Comm.Image Rep,2001,12(4):436-449.

[11]ESEDOGLU S,SHEN J.Digital inpainting based on the Mumford-Shah-Euler image model[J].European Journal of Applied Mathematics,2002,13:353-370.

[12]MUMFORD D,SHAH J.Optimal approximations by piecewise smooth functions and associated variational problems[J].Communications on Pure and Applied Mathematics,1989,42: 577-685.

[13]GROSSAUER H,SCHERZER O.Using the complex Ginzburg-Landau equation for digital inpainting in 2d and 3d,scale space methods in computer vision[J].Lecture Notes in Computer Science 2695,2003: 225-236.

[14]BORNEMANN F,MRZ T.Fast image inpainting based on coherence transport[J].Journal of Mathematical Imaging and Vision,2007,28(3): 259-278.

[15]BERTALMIO M,VESE L,SAPIRO G,et al.Simultaneous texture and structure image inpainting[J].IEEE Tran Image Process.2003,12(8): 882-889.

[16]CRIMINISI A,PEREZ P,TOYAMA K.Region filling and object removal by exemplar-based image inpainting[J].IEEE Trans Image Processing,2004,13(9): 1200-1212.

[17]WEXLER Y,SHECHTMAN E,IRANI M.Space-time video completion[J].IEEE Trans Pattern Anal Mach Intell,2007,29(3): 1463-1476.

[18]FADILI M,STARCK J,MURTAGH F.Inpainting and zooming using sparse representations[J].Comput J,2009,52(1):64-79.

[19]KOMODAKIS N,TZIRITAS G.Image completion using efficient belief propagation via priority scheduling and dynamic pruning[J].IEEE Transactions on Image Processing,2007,16(11): 2649-2661.

[20]AUJOL J,LADJAL S,MASNOU S.Exemplar-based inpainting from a variational point of view[J].SIAM Journal on Mathematical Analysis,2010,42(3): 1246-1285.

[21]MAJDA A,BERTOZZI A.Vorticity and incompressible flow [M].Cambridge Univ Press,2001.

[22]TAI X.Global extrapolation with a parallel splitting method[J].Numerical Algorithm,1991,3:527-440.

[23]WEICKERT J,ROMENY B,VIERGEVER M.Efficient and reliable schemes for nonlinear diffusion filtering[J].IEEE Trans on Image Processing,1998,7(3):398-410.

[24]KUHNE G,WEICKERT J,VIERGEVER M.Fast implicit active contours models lecture[J].Notes on Computer Science,2002,2449: 133-140.