基于Ginzburg-Landau泛函二值圖像復(fù)原的快速算法

孫寶宸，李 君，?；圪e

（天津師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院，天津 300387）

在圖像的采集和傳輸?shù)冗^程中，噪聲污染和信息缺失會(huì)導(dǎo)致圖像質(zhì)量降低.圖像噪聲有多種類型[1]，如高斯噪聲、泊松噪聲、椒鹽噪聲等.噪聲的多樣性使得圖像去噪成為一項(xiàng)富有挑戰(zhàn)且被廣泛關(guān)注的研究課題，相關(guān)學(xué)者提出了多種經(jīng)典的去噪方法[2-3].

1992年，Rudin等提出的一種以全變差作為正則項(xiàng)的凸ROF模型[2]，可以通過求解其歐拉-拉格朗日方程得到全局最優(yōu)解，且解是唯一的，之后，一系列快速求解的一階優(yōu)化算法被相繼提出[4].針對(duì)具體噪聲的統(tǒng)計(jì)性質(zhì)，可以通過保留全變差正則項(xiàng)改進(jìn)不同的保真項(xiàng)度量來提高復(fù)原圖像的質(zhì)量.ROF模型具有良好的性質(zhì)，但是全變差項(xiàng)的不可微性使算法設(shè)計(jì)變得困難，為解決這一問題，一系列光滑化技術(shù)被相繼提出.文獻(xiàn)[5]首次提出使用一類基于雙阱勢(shì)函數(shù)（doublewell potential function）的修正模型（簡(jiǎn)稱DP模型）處理圖像銳化問題，之后該模型也被應(yīng)用于圖像增強(qiáng)[6]、修復(fù)[7-8]和去噪[9]等問題.在文獻(xiàn)[10-11]中，基于非局部建模的DP模型首次用于圖像分割與數(shù)據(jù)分類問題.DP模型是逐點(diǎn)可微但非凸的，因此對(duì)模型求解具有雙重困難：一是直接求解歐拉-拉格朗日方程只能得到局部最優(yōu)解而無法保證達(dá)到全局最優(yōu)解；二是DP模型的一階變分是非線性的，很難設(shè)計(jì)有效的優(yōu)化算法和梯度流迭代格式.對(duì)于全隱梯度流格式，文獻(xiàn)[12]證明了若時(shí)間步長(zhǎng)滿足一定的約束條件，則可以保證能量穩(wěn)定；對(duì)于顯隱格式，在相場(chǎng)模型的求解中也得到了類似的結(jié)論[13].

近年來，一類保能量衰減的顯隱迭代格式被應(yīng)用于相場(chǎng)模型的計(jì)算中，如，乘子方法、IEQ格式和SAV格式[13].這類迭代格式通過引入輔助變量改寫目標(biāo)能量泛函，將非線性梯度流方程轉(zhuǎn)化為多變量線性系統(tǒng)，但在實(shí)際計(jì)算中沒有增加需要求解的變量，并能保證改寫的等價(jià)能量泛函隨迭代進(jìn)行是衰減的.此外，交替方向乘子算法在圖像處理的反問題求解中取得了巨大的成功，但目前未見其應(yīng)用于求解DP模型.

本文首先引入基于Ginzburg-Landau泛函松弛的雙阱勢(shì)模型和2種經(jīng)典算法，進(jìn)而通過構(gòu)造特殊的優(yōu)化問題，建立了Lie格式與MBO格式的聯(lián)系，對(duì)MBO格式提出了一種從算子分裂優(yōu)化技術(shù)角度的理解方式；然后通過引入輔助變量設(shè)計(jì)了一類新的算子分裂算法，其子優(yōu)化問題具有閉形式解；最后通過二值圖像復(fù)原實(shí)驗(yàn)評(píng)估各種算法的有效性.

1 雙阱勢(shì)模型和經(jīng)典算法

1.1 雙阱勢(shì)模型

ROF模型[2-3]的連續(xù)形式定義為

其中：λ為模型參數(shù)；Du為u的分布導(dǎo)數(shù)；f為污染圖像；u為待復(fù)原的圖像；Ω為二維歐氏空間的一個(gè)有界閉集.全變差函數(shù)空間為

顯然E1ROF（u）關(guān)于u是下凸的（若無特殊說明，本文所指的凸性均為下凸性）.對(duì)于上述優(yōu)化問題直接設(shè)計(jì)迭代格式求解歐拉-拉格朗日方程[4]，即式（1），獲得唯一的最優(yōu)解.

其中：ε為正常數(shù)；W（u）為雙阱勢(shì)函數(shù)（如W（u）=0.25（u2-1）2，它有2個(gè)全局極小值）.文獻(xiàn)[16]證明了當(dāng)ε→0+時(shí)，GL（u）會(huì)Γ收斂[3]到u的全變差.

基于上述結(jié)果將BV（Ω）松弛到H1（Ω）可建立雙阱勢(shì)模型（DP模型）如下

1.2 MBO格式（全隱格式）

文獻(xiàn)[17]提出了一種有效的算子分裂格式（記為MBO）.采用MBO格式進(jìn)行變量更新需要2步迭代：第1步利用向后Euler差分格式求解反應(yīng)擴(kuò)散方程ut=Δu；第2步使用1/2水平集進(jìn)行閾值化.由于DP模型中W（u）項(xiàng)含有2個(gè)駐點(diǎn)0和1，所以第2步迭代可表達(dá)為

將MBO應(yīng)用于求解DP模型得到一種修正的MBO格

1.3 凸分裂格式

對(duì)于非凸優(yōu)化問題，若通過能量恒等變換可將非凸的能量泛函轉(zhuǎn)化為凸能量與凹能量之和的形式，則凸分裂格式可視為一種簡(jiǎn)單有效的數(shù)值求解算法，但它通常只有一階收斂速度.文獻(xiàn)[7]將DP模型分解為EDP（u）=E1（u）+E2（u），其中E1（u）=‖Δu‖2+μ‖u（u-1）‖2，為構(gòu)造凸能量，在不改變目標(biāo)能量的前提下取E（2u）=λ‖f-u‖2，對(duì)E（iu）（i=1、2）設(shè)計(jì)凸分裂格式.將E（iu）分解為E（iu）=E（i1u）-E（i2u），其中：

文獻(xiàn)[7]證明了上述凸分裂格式（convex splitting with 2 mixtures，CS2M）是無條件能量穩(wěn)定的，即對(duì)任意的正整數(shù)n，有EDP（un+1）≤EDP（un）.

2 快速迭代求解方法

本節(jié)設(shè)計(jì)一些高效求解DP模型的迭代格式.首先將界面問題的新算法推廣到圖像復(fù)原問題的求解中，然后設(shè)計(jì)一類基于交替方向乘子方法（alternating direction method of multipliers，ADMM）的新算法.本文算法中均假設(shè)u滿足Neumann邊界條件，即0，其中：Γ為區(qū)域Ω的邊界，n為單位外法向量.

2.1 相場(chǎng)算法的推廣

利用梯度流方程求解能量泛函F的極小化問題，其一般框架為[13]

2.1.1 EQCN格式

文獻(xiàn)[21-22]對(duì)一類具有特殊結(jié)構(gòu)的目標(biāo)泛函設(shè)計(jì)了二階的顯隱差分格式，其核心是通過引進(jìn)輔助變量將能量泛函的非平方項(xiàng)轉(zhuǎn)化為二次形式，進(jìn)一步設(shè)計(jì)能量穩(wěn)定的時(shí)間、空間全離散數(shù)值格式，本文考慮其時(shí)間半離散格式.定義輔助變量V=u（u-1），將DP模型進(jìn)行變量分解，則式（3）轉(zhuǎn)化為

式（11）的梯度流AC方程組為

文獻(xiàn)[22]通過引入外插和Crank-Nicolson（CN）中心插商構(gòu)建了迭代格式較為簡(jiǎn)單且具有二階收斂速度的EQCN格式.記

利用式（14）進(jìn)行變量更新，最后使用式（4）進(jìn)行閾值化.EQCN格式可以保證能量穩(wěn)定且具有二階收斂速度[23].

2.1.2 SAVCN格式

文獻(xiàn)[13]中提出了具有一階收斂速度的向前Euler格式、高階CN和BDF時(shí)間半離散格式，本文主要考慮SAVCN格式，其變量更新滿足

結(jié)合式（15）和式（16）對(duì)線性系統(tǒng)進(jìn)行迭代求解，最后使用式（4）進(jìn)行閾值化.SAVCN格式可以保證能量穩(wěn)定且具有二階收斂速度[13].

2.2 交替方向乘子方法

交替方向乘子方法[24-26]可以處理非凸非光滑的復(fù)雜約束優(yōu)化問題.考慮如下優(yōu)化問題

上述框架可以自然地推廣到無窮維空間上.基于此，本文運(yùn)用算子分裂技術(shù)，通過引入輔助變量將無約束優(yōu)化問題轉(zhuǎn)化為約束優(yōu)化問題，設(shè)計(jì)了2種迭代格式，其子優(yōu)化問題具有閉形式解.

2.2.1 兩變量交替方向乘子算法

引入輔助變量v使得u=v，將W（u）進(jìn)行算子分裂，原無約束優(yōu)化問題可修改為約束優(yōu)化問題：

通過交替優(yōu)化變量和乘子可以得到求解DP問題的新算子分裂格式，其更新框架為

由于其目標(biāo)泛函是凸可微的，因此可由一階最優(yōu)性條件

同u子問題的求解過程一致，可得

乘子Λ依照式（18c）進(jìn)行更新.滿足終止條件后使用式（4）對(duì)u進(jìn)行閾值化.上述優(yōu)化算法簡(jiǎn)記為ADMM2V，其算法流程見表1.

表1 ADMM2V算法流程Tab.1 Algorithm flow of ADMM2V

2.2.2 三變量交替方向乘子算法

類比EQCN格式的解耦思想，本文設(shè)計(jì)了一種三變量交替方向乘子算法.引入2個(gè)輔助變量w和v將W（u）進(jìn)行算子分裂，無約束優(yōu)化問題可以轉(zhuǎn)化為約束優(yōu)化問題：

其中：r1和r2為正的增廣拉格朗日系數(shù)，Λ1和Λ2為乘子.通過交替優(yōu)化變量和乘子可以得到求解DP問題的新算子分裂格式，其更新框架為

對(duì)于乘子Λ1和Λ2按照式（21d）進(jìn)行更新，滿足終止條件后使用式（4）對(duì)u進(jìn)行閾值化.上述優(yōu)化算法簡(jiǎn)記為ADMM3V，其算法流程見表2.

表2 ADMM3V算法流程Tab.2 Algorithm flow of ADMM3V

3 數(shù)值實(shí)驗(yàn)

由于更新格式（12）和（15）考慮的是Neumann邊界條件，因此本文對(duì)離散圖像添加鏡像邊界，依據(jù)文獻(xiàn)[27]，對(duì)連續(xù)的梯度算子離散得Δu（i，j）=（Δux（i，j），Δuy（i，j）），令Nx和Ny分別為離散圖像矩陣的行數(shù)和列數(shù)，其中：

其中：p=（p1，p2），p1=Δux，p2=Δuy.外迭代終止條件設(shè)為迭代1 000次或變量u的相對(duì)殘差小于10-4，即‖un+1-un‖/‖un‖＜10-4.對(duì)子問題的線性方程組使用共軛梯度法迭代求解，并設(shè)內(nèi)迭代終止條件為迭代10次或相對(duì)誤差小于10-10.

3.1 圖像去噪

本文主要測(cè)試了3種加性噪聲（高斯噪聲、椒鹽噪聲和泊松噪聲）的復(fù)原問題.盡管針對(duì)不同種類的噪聲，可以設(shè)計(jì)不同的保真度量以獲得更好的去噪效果，但對(duì)于不可微的保真項(xiàng)會(huì)使基于梯度信息的SAVCN和EQCN格式失效，因此在圖像去噪實(shí)驗(yàn)中統(tǒng)一采用L2模數(shù)作為數(shù)據(jù)保真項(xiàng).采用信噪比（SNR）評(píng)價(jià)去噪效果，

其中：u為去噪后圖像，Img為無噪聲污染圖像.分別采用指紋、條形碼和二維碼圖像（Nx=Ny=256）進(jìn)行高斯去噪、椒鹽去噪和泊松去噪實(shí)驗(yàn)，計(jì)算mMBO、CS2M、EQCN、SAVCN、ADMM2V和ADMM3V等6種算法的SNR.基于2.2中對(duì)mMBO格式的解析，只討論其他5種算法的數(shù)值收斂性和能量衰減性.

對(duì)二值化的指紋圖像添加服從正態(tài)分布N（0，σ2）的高斯噪聲，取σ=0.2、0.5和0.8的噪聲圖像，見圖1，各算法復(fù)原圖像的SNR見表3，6種算法SNR的最大值用黑體標(biāo)出，下同.

圖1 指紋真實(shí)圖像和高斯噪聲圖像Fig.1 Real image and Gaussian noisy images of fingerprint

由表3可見，除mMBO外，其他5種算法的SNR基本相同，均高于mMBO.圖2給出了變量u的相對(duì)誤差和總能量EDP（u）的變化曲線.由圖2（a）—（c）可見，5種算法的收斂速度沒有顯著差異.由圖2（d）—（f）可見，各算法都保證了能量衰減性，但能量衰減的速度不一致.當(dāng)σ=0.2時(shí)，EQCN和SAVCN能量下降得最快；當(dāng)σ=0.5時(shí)，各算法的能量降速基本一致，當(dāng)σ=0.8時(shí)，ADMM能量下降得最快.

表3 高斯去噪實(shí)驗(yàn)的SNRTab.3 SNR of image denoising experiments of Gaussian noisydB

圖2 高斯去噪實(shí)驗(yàn)相對(duì)誤差與總能量的變化情況Fig.2 Changes of relative residuals and total energy in Gaussian denoising experiment

對(duì)二值化的條形碼圖像添加強(qiáng)度為ρ的椒鹽噪聲，取ρ=0.05、0.10和0.20的噪聲圖像，見圖3，各算法復(fù)原圖像的SNR見表4.由表4可見，除mMBO外，其他5種算法的SNR基本相同，均高于mMBO.圖4給出了變量u的相對(duì)誤差和總能量EDP（u）的變化曲線.由圖4（a）—（c）可見，5種算法的收斂速度沒有顯著差異.由圖4（d）—（f）可見，各算法都保證了能量衰減性，對(duì)于ρ=0.05、0.10和0.20，均為ADMM能量下降得最快.

圖3 條形碼真實(shí)圖像與椒鹽噪聲圖像Fig.3 Real image and salt and pepper noisy images of barcode

表4 椒鹽去噪實(shí)驗(yàn)的SNRTab.4 SNR of image denoising experiments of salt and pepper noisy dB

圖4 椒鹽去噪實(shí)驗(yàn)相對(duì)誤差與總能量的變化情況Fig.4 Changes of relative residuals and total energy in salt and pepper denoising experiment

對(duì)二值化的二維碼圖像逐點(diǎn)添加泊松噪聲，服從Poisson（ηf（i，j））/η分布，取η=4、8和16的噪聲圖像，見圖5，各算法復(fù)原圖像的SNR見表5.由表5可見，除mMBO外，其他5種算法的SNR均較高，ADMM2V的SNR最高.圖6給出了變量u的相對(duì)誤差和總能量EDP（u）的變化曲線.由圖6（a）—（c）可見，5種算法的收斂速度沒有顯著差異，但在優(yōu)先保證EDP（u）充分下降的前提下，2種ADMM算法的收斂速度較快.由圖6（d）—（f）可見，各算法都保證了能量衰減性，對(duì)于η=4、8和16，均為ADMM能量下降得最快.

圖6 泊松去噪實(shí)驗(yàn)相對(duì)誤差與總能量的變化情況Fig.6 Changes of relative residuals and total energy in Poisson denoising experiment

表5 泊松去噪實(shí)驗(yàn)的SNRTab.5 SNR of image denoising experiments of Poisson noisy dB

圖5 二維碼真實(shí)圖像與泊松噪聲圖像Fig.5 Real image and Poisson noisy images of QR code

3.2 圖像修復(fù)

對(duì)于圖像修復(fù)問題需要將DP模型中的λ∈R+修改為

其中：Ωb為待修復(fù)部分；λ0為一個(gè)正常數(shù).考慮簡(jiǎn)單結(jié)構(gòu)和復(fù)雜結(jié)構(gòu)2類圖像的修復(fù)問題，以修復(fù)后圖像與真實(shí)圖像的均方誤差（MSE）評(píng)價(jià)修復(fù)效果，

在簡(jiǎn)單結(jié)構(gòu)圖像修復(fù)實(shí)驗(yàn)中，對(duì)二值化的矩形拼接圖（Nx=Ny=128）設(shè)置不同比例的區(qū)域信息缺失，缺失比例分別設(shè)為30%、50%和70%，見圖7，各算法修復(fù)圖像的MSE見表6，6種算法MSE的最小值用黑體標(biāo)出，下同.

圖7 矩形拼接圖的真實(shí)圖像與待修復(fù)圖像Fig.7 Real image of rectangular and images to be inpainting

由表6可見，SAVCN算法的修復(fù)效果最優(yōu).圖8給出了變量u的相對(duì)誤差和總能量EDP（u）的變化曲線.由圖8（a）—（c）可見，在迭代初期，2種ADMM算法收斂較快，而隨著迭代次數(shù)增多，SAVCN格式的收斂速度超過其他算法.由圖8（d）—（f）可見，各算法都保證了能量衰減性，但CS2M的總能量隨迭代次數(shù)增多先上升后下降，最終達(dá)到穩(wěn)態(tài).

表6 矩形拼接圖像修復(fù)實(shí)驗(yàn)數(shù)值結(jié)果（MSE）Tab.6 Numerical results（MSE）of rectangular image inpainting experiments

圖8 矩形拼接圖像修復(fù)數(shù)值實(shí)驗(yàn)相對(duì)誤差與總能量的變化情況Fig.8 Changes of relative residuals and total energy in image inpainting numerical experiments for rectangular

在復(fù)雜結(jié)構(gòu)圖像修復(fù)實(shí)驗(yàn)中，對(duì)條形碼和二維碼圖像（Nx=Ny=256）分別設(shè)置不同的局部區(qū)域信息缺失，見圖9，各算法修復(fù)圖像的MSE見表7.

圖9 條形碼和二維碼的真實(shí)圖像與待修復(fù)圖像Fig.9 Real images of barcode and QR code and images to be inpainting

由表7可見，SAVCN算法的修復(fù)效果最優(yōu).圖10給出了變量u的相對(duì)誤差和總能量EDP（u）的變化曲線.由圖10（a）和（b）可見，對(duì)于條形碼圖像，CS2M收斂最快，對(duì)于二維碼圖像，2種ADMM算法和CS2M收斂速度基本一致，高于其他算法.由圖10（c）和（d）可見，各算法均可保證目標(biāo)總能量最終達(dá)到穩(wěn)態(tài).

圖10 條形碼和二維碼圖像修復(fù)數(shù)值實(shí)驗(yàn)相對(duì)誤差與總能量變化情況Fig.10 Changes of relative residuals and total energy in image inpainting numerical experiments for barcode and QR code

表7 條形碼和二維碼圖像修復(fù)實(shí)驗(yàn)數(shù)值結(jié)果（MSE）Tab.7 Numerical results（MSE）of barcode and QR code images inpainting experiments

對(duì)于以上圖像去噪和修復(fù)實(shí)驗(yàn)，圖11給出了ADMM和SAVCN算法中超參數(shù)對(duì)圖像復(fù)原效果的影響.由圖11（a）可見，ADMM2V的參數(shù)r2L對(duì)最終的結(jié)果數(shù)值影響是敏感的.由圖11（b）和（c）可見，SAVCN的參數(shù)τ和ADMM3V的參數(shù)r1、r2對(duì)最終的結(jié)果數(shù)值影響均不敏感.

圖11 算法參數(shù)對(duì)圖像復(fù)原結(jié)果的影響Fig.11 Influences of algorithm parameters on image restoration results

4 結(jié)論

本文利用DP模型處理二值圖像復(fù)原問題，主要考慮了三種成熟的半隱迭代格式和兩種新的ADMM算法.數(shù)值算例的結(jié)果表明，本文所提出的ADMM方法去噪實(shí)驗(yàn)中可以恢復(fù)到信噪比更高的結(jié)果且保持能量穩(wěn)定，在修復(fù)實(shí)驗(yàn)中，SAVCN格式最終修復(fù)所得結(jié)果更好，但其收斂速度較慢

鑒于ADMM算法在圖像去噪實(shí)驗(yàn)中表現(xiàn)出穩(wěn)定的優(yōu)勢(shì)，因此在后續(xù)的研究中，將對(duì)ADMM算法的收斂性進(jìn)行分析.除此之外，考慮到SAVCN格式在圖像修復(fù)中表現(xiàn)良好，從而將在后續(xù)研究中設(shè)計(jì)更為穩(wěn)定的能量衰減迭代格式.