緊框架小波和總廣義全變分聯(lián)合約束的醫(yī)學(xué)圖像復(fù)原算法

張 晶， 馬 瑾， 邵 晨， 桂志國， 張 權(quán)， 楊 婕,

(1. 山西財貿(mào)職業(yè)技術(shù)學(xué)院 信息工程系， 山西 太原 030031； 2. 山西中醫(yī)藥大學(xué) 醫(yī)藥管理學(xué)院， 山西 太原 030619； 3. 中北大學(xué) 信息與通訊工程學(xué)院， 山西 太原 030051)

0 引 言

在圖像的成像、 采集和傳輸過程中， 成像設(shè)備或其他物理因素會不可避免地造成圖像質(zhì)量退化， 為使圖像質(zhì)量滿足實際應(yīng)用(如醫(yī)療、 軍事、 交通等)要求， 研究圖像降噪、 圖像去模糊化、 圖像修復(fù)以及圖像重建等圖像復(fù)原技術(shù)[1-2]非常有必要. 通常情況下， 圖像復(fù)原可以看作一個典型的通過對退化圖像的近似數(shù)學(xué)估計模型進(jìn)行求逆的不適定問題

f=Aλ+η,(1)

式中：f表示觀測到的退化圖像；A是一個病態(tài)的線性算子；η為加性隨機(jī)噪聲. 由于設(shè)備和噪聲等因素的限制， 很難構(gòu)建模擬圖像退化的準(zhǔn)確數(shù)學(xué)模型； 即使能夠建立準(zhǔn)確的估計模型， 數(shù)據(jù)的大規(guī)模性及算法模型的非光滑性也會致使有效的優(yōu)化估計方法很難獲得. 因此， 圖像復(fù)原技術(shù)所需解決的兩個關(guān)鍵問題為圖像質(zhì)量退化模型的建立和相應(yīng)優(yōu)化估計算法的設(shè)計[3].

對于退化模型的構(gòu)建， 國內(nèi)外學(xué)者的研究焦點集中在由保真項和正則項構(gòu)成的能量泛函約束的正則化復(fù)原模型上. 就保真項而言， 由高斯、 泊松噪聲所致的圖像退化常由光滑的泛函L2范數(shù)來擬合； 由椒鹽噪聲引發(fā)的圖像退化則由非光滑的泛函L1范數(shù)來擬合[3，4]. 就正則項而言， 通過利用如全變分、 小波、 剪切波、 緊框架小波等變換的稀疏性， 圖像的先驗信息可以被很好地表征[5]. 對于退化模型的求解， 主要有直接方法和分解方法. 直接方法同時處理保真項和擬合項， 如果該兩項均為光滑泛函， 可通過梯度下降法和牛頓迭代法求解； 如果僅保真項是光滑的， 則當(dāng)正則項為有界全變分時可將模型轉(zhuǎn)換為變分偏微分方程， 進(jìn)而進(jìn)行數(shù)值求解； 此外， 根據(jù)保真項和擬合項的特性， 還可通過Fenchel變換將模型轉(zhuǎn)換為對偶問題求解等[4,6，7]. 從式(1)所示的線性不適定的逆問題中恢復(fù)高質(zhì)量圖像的關(guān)鍵是選取能夠合理表征圖像先驗的正則項， 近年比較成功的正則項約束的復(fù)原方法有Rudin-Osher-Fatemi模型[8]、 非局部變型[9，10]、 inf-卷積模型[11]、 TGV模型[12]、 聯(lián)合一階、 二階全變分模型[13]以及Gabor框架、 小波框架、 剪切波等多尺度應(yīng)用諧波分析方法[1,14，15]. 一般地， 基于泛函Lp(p=0,1,2)的范數(shù)正則項約束的復(fù)原模型等價于求解最小化問題

(2)

式中：Φ表示某種具有稀疏性的線性變換.

為了解決由全變分函數(shù)的“分片常數(shù)”效應(yīng)導(dǎo)致的階梯狀偽邊緣、 由稀疏逼近算法的精確性導(dǎo)致的尖銳邊緣及細(xì)節(jié)信息丟失等問題， 本文通過構(gòu)造能夠完備表征圖像特性的稀疏正則項及設(shè)計目標(biāo)函數(shù)迭代優(yōu)化估計方法來提高復(fù)原圖像的質(zhì)量. 僅采用一種小波基函數(shù)很難準(zhǔn)確描述圖像特征， 通過構(gòu)造能夠利用小波框架函數(shù)族對圖像進(jìn)行稀疏表示的緊框架小波來提高圖像描述的準(zhǔn)確性； 此外， TGV函數(shù)具有凸性、 旋轉(zhuǎn)不變性以及下半函數(shù)連續(xù)性， 通過逼近由多項式組成的任意階函數(shù)， 可以從原理上克服TV方法的“分片常數(shù)”效應(yīng)， 進(jìn)而解決階梯狀偽邊緣問題. 由此， 本文構(gòu)造了一種由緊框架小波和二階TGV組成的稀疏正則項， 并以所構(gòu)造的正則項為約束進(jìn)行圖像復(fù)原. 由于新的復(fù)原模型由光滑L2范數(shù)和非光滑L1范數(shù)泛函共同組成， 故直接求解難度較高. 本文首先通過輔助變量將具有混合泛函的目標(biāo)函數(shù)最小化問題分解為兩個子問題， 進(jìn)而分別采用均值增廣拉格朗日算法和Chambolle-Pock原始—對偶迭代方法獲得最優(yōu)解.

1 基本理論

1.1 緊框架小波

對于任意一個可數(shù)集χ?L2(Rd,d∈N)，L2(Rd) 上的緊框架[16-17]定義為

Rd),(3)

式中： 〈·,·〉表示L2(Rd)上的內(nèi)積， 且〈f,φ〉稱為f的典型系數(shù). 若給定Ψ={ψl∶l=1,…,r}?L2(Rd)， 則

χN∈N(Ψ)={ψl,n,k∶1≤l≤r,n∈Z,k∈Zd},(4)

則L2(Rd)上的緊框架小波χN(Ψ)的每一個緊框架ψ1,…,ψr函數(shù)由對稱緊支撐函數(shù)組成， 若

則基于多分辨率分析的緊框架Ψ構(gòu)造可通過尋找有限支撐的ql使其滿足

[k]φ(2x-k),

l=1,…,r.(7)

基于酉擴(kuò)張原理(Unitary Extension Principle， UEP)可得

Rd,(8)

例如，L2(Rd)上的分段線性B樣條(B-spline)函數(shù)為

(10)

二維離散快速緊框架變換的定義為

Wλ={Wl,αλ∶(l,α)∈Γ},(11)

式中： Γ=({0,…,L-1}×B)∪{(L-1,0)}， 緊支撐子帶B={0,…,r}2{0}， 則λ∈I2在第l層α子帶的緊支撐系數(shù)為

Wl,αλ=ql,α[-·](*)λ,(12)

此處， (*)表示具有一定邊界條件的離散卷積， 且

(13)

λ=WTWλ.(14)

1.2 總廣義全變分

近年來， 被廣泛采用的基于全變差(Total Variation， TV)正則項約束的復(fù)原模型雖然能夠保持圖像的尖銳邊緣， 但會產(chǎn)生塊狀或階梯狀偽邊緣. 針對這一不足， 一種廣義TV即TGV法被引入到圖像復(fù)原處理中. 與僅考慮一階導(dǎo)數(shù)的TV不同， TGV能夠考慮高階導(dǎo)數(shù)， 對于分段多項式函數(shù)或尖銳邊緣處理效果較好.k階TGV定義[18]為

‖divlφ‖∞≤αl},(15)

(17)

(18)

廣義變分有界空間定義為

(19)

2 算法模型及優(yōu)化估計

2.1 算法簡介

通過構(gòu)造能夠完備表征圖像特性的稀疏正則項， 本文提出了一種緊框架小波和總廣義全變分聯(lián)合約束的圖像復(fù)原算法. 算法模型為

s.t. ?=Wλ,(21)

式中： 第一項為保真項， 確保復(fù)原圖像盡可能保留原圖像的重要特征； 第二項和第三項為聯(lián)合正則項. 本文方法的設(shè)計主要是為了合理發(fā)揮緊框架小波和二階TGV表征待復(fù)原圖像先驗信息的優(yōu)越性. 一方面， 二階TGV能夠在平滑噪聲的同時克服“階梯偽影效應(yīng)”； 另一方面， 緊框架小波變換具有稀疏性， 能夠抑制重要紋理細(xì)節(jié)信息的丟失. 如式(21)所示， 本文方法所需解決的問題為條件約束的最小化問題， 采用雙重增廣拉格朗日算法[19](Doubly augmented Lagrangian， DAL)將該有約束的最小化問題變?yōu)闊o約束的問題， 即

(22)

其中，

(23)

由于該含光滑L2范數(shù)和非光滑L1范數(shù)泛函的最小化復(fù)原模型求解難度較高. 本文通過引入輔助變量μ來簡化求解過程， 即式(23)等價于

(24)

進(jìn)而可將式(24)分解為兩個子問題， 即

(25)

對于子問題(1)， 通過進(jìn)一步計算可得

(26)

進(jìn)而有

(27)

?j+1=Hβ1,κ,γ(Wλj+1+υj,?j),(28)

υj+1=υj+Wλj+1-?j+1.(29)

通過廣義硬閾值收縮求解式(28)， 得

(30)

令

采用一階原始對偶算法可以得出

其中

(33)

進(jìn)一步進(jìn)行迭代求解

(34)

其中，

2.2 算法流程

算法任務(wù)：

緊框架小波和總廣義全變分聯(lián)合約束的圖像復(fù)原算法.

參數(shù)初始化

While 不滿足停步準(zhǔn)則

第一步 求解子問題(1)， 給定μj，

1. 更新λj+1， 采用公式(27)；

2. 更新?j+1， 采用公式(28)；

3. 更新υj+1， 采用公式(29)；

第二步 求解子問題(2)， 采用求出的λj+1，

5. 更新μj+1， 采用公式(34)；

j=j+1

結(jié)束

輸出質(zhì)量改善的復(fù)原圖像

3 實驗與分析

本文實驗對象為以 “l(fā)ena”和“brain_tumor”命名的標(biāo)準(zhǔn)圖像(如圖 1 所示)， 對噪聲水平為σ=6時的退化圖像(如圖 2 所示)進(jìn)行復(fù)原實驗.

圖 1 本文實驗采用的標(biāo)準(zhǔn)圖像Fig.1 The standard images used in the experiments

圖 2 σ=6時的質(zhì)量退化圖像Fig.2 When σ=6, the degraded images used in the experiments

實驗中計算機(jī)硬件配置為： Intel(R) Core(TM) i5-6500 CPU @ 3.20 GHz， Windows7 32位操作系統(tǒng)； 軟件環(huán)境為： MATLAB R2014a. 復(fù)原圖像質(zhì)量量化評價指標(biāo)為峰值信噪比(Peak Signal to Noise Ratio， PSNR)和結(jié)構(gòu)相似度(Structural Similarity Index， SSIM)， 定義為

(36)

其中， 式(39)目標(biāo)函數(shù)的正則項為TV， 式(40)目標(biāo)函數(shù)的正則項為二階TGV. 本文算法迭代停步條件為

4.5×10-4.

圖 3 給出了分別采用TV、 TGV以及本文方法對退化lena進(jìn)行復(fù)原的實驗結(jié)果. 可以發(fā)現(xiàn)， TV的“分片常數(shù)”導(dǎo)致圖3(b)中有大量塊狀和階梯狀的偽邊緣出現(xiàn)； 圖3(c)很好地說明了二階TGV在抑制噪聲和階梯偽影方面的優(yōu)勢； 從圖3(d) 可以看出本文方法在細(xì)節(jié)和紋理信息方面復(fù)原情況更好. 圖4為圖3中各實驗結(jié)果圖中局部放大圖像， 直觀地， 如圖中箭頭所示， 本文方法可以很好地改善TV結(jié)果圖1中塊狀和階梯狀偽邊緣情況， 同時在細(xì)節(jié)和紋理結(jié)構(gòu)的復(fù)原方面又優(yōu)越于TGV方法.

圖 3 采用不同算法復(fù)原后的lena圖像Fig.3 The restored lena images obtained by different methods

為了進(jìn)一步說明本文所提方法的有效性， 對各復(fù)原圖的不同感興趣區(qū)(Region of Interest， ROI)進(jìn)行標(biāo)化， 如圖 5 所示， 并對其量化指標(biāo)表現(xiàn)情況進(jìn)行比較， 如表 1 所示. 結(jié)果標(biāo)明： 在ROI1和ROI2內(nèi)， 本文方法結(jié)果圖的PSNR值最大， 表明通過本文方法所獲復(fù)原圖失真較小， 與原始的標(biāo)準(zhǔn)圖像更加接近； 類似地， 本文方法結(jié)果圖的SSIM值也最大， 說明本文方法復(fù)原圖與原始標(biāo)準(zhǔn)圖像的結(jié)構(gòu)更加相似.

圖 4 采用不同算法復(fù)原后的lena局部圖像Fig.4 The partial images of the restored lena images in Fig.3

圖 5 不同感興趣區(qū)的lena和brain_tumor 圖Fig.5 The lena and brain_tumor images with region of interest (ROI) marked

算法PSNRSSIMROI1ROI2ROI1ROI2本文算法37．625925．46350．96110．8010TGV34．474923．67970．95980．7856TV34．314323．66260．93890．7776

圖 6 為分別采用TV、 TGV以及本文方法獲取的brain_tumor圖像復(fù)原圖， 圖7為圖6中各個結(jié)果圖像的局部放大圖. 表 2 給出了基于brain_tumor復(fù)原圖采用不同方法的量化指標(biāo)表示情況. 一方面， 觀察圖6和圖7可以發(fā)現(xiàn)， 與其他兩種方法相比較， 本文方法所獲結(jié)果圖像的細(xì)節(jié)和紋理信息復(fù)原情況明顯更好. 另一方面， 分析表2可以看出， 對于brain_tumor圖像， 在ROI3和ROI4內(nèi)， 本文方法結(jié)果圖的PSNR值和SSIM值也是最大的， 說明本文方法復(fù)原圖失真較小， 與原始標(biāo)準(zhǔn)圖像的結(jié)構(gòu)更加相似. 綜上， 無論從直觀視覺角度還是從客觀量化分析角度， 本文方法較傳統(tǒng)方法都表現(xiàn)出一定的優(yōu)越性， 是可行和有效的.

圖 6 采用不同算法復(fù)原后的brain_tumor圖像Fig.6 The restored brain_tumor images obtained by different methods

圖 7 采用不同算法復(fù)原后的brain_tumor局部圖像Fig.7 The partial images of the restored brain_tumor images in Fig.6

算法PSNRSSIMROI1ROI2ROI1ROI2本文算法27．573028．50380．82310．8378TGV24．655525．34440．80480．8157TV24．161424．34370．80290．8097

4 結(jié)束語

本文提出一種緊框架小波和總廣義全變分聯(lián)合約束的圖像復(fù)原算法. 首先， 構(gòu)造出一種由緊框架小波L1范數(shù)和二階TGVL2范數(shù)組成的聯(lián)合正則項約束的圖像復(fù)原模型； 其次， 采用交替方向迭代方法將所建模型的最小化問題分解為兩個子問題， 并分別采用均值增廣拉格朗日算法和Chambolle-Pock一階原始—對偶迭代方法獲得最優(yōu)解； 通過與其他復(fù)原算法進(jìn)行仿真實驗分析， 不管從視覺效果還是從量化指標(biāo)分析來看， 本文方法不僅可以有效地抑制噪聲， 而且在去模糊化后能夠很好地保留紋理和細(xì)節(jié)信息， 使得復(fù)原后的結(jié)果圖更接近原始圖像. 因此， 該算法的提出有助于圖像復(fù)原技術(shù)進(jìn)一步的發(fā)展.

[1] Dong B, Jiang Q, Shen Z. Image restoration: wavelet frame shrinkage, nonlinear evolution PDEs, and beyond[J]. Multiscale Modeling & Simulation, 2017, 15(1)： 606-660.

[2] Katsaggelos A K. Digital image restoration[M]. New York: Springer Publishing Company, 2012.

[3] 李旭超, 馬松巖, 邊素軒. 對偶算法在緊框架域TV-L1去模糊模型中的應(yīng)用[J]. 中國圖象圖形學(xué)報, 2015, 20(11)： 1434-1445.

Li Xuchao, Ma Sunyan, Bian Suxuan. Application of dual algorithm to TV-L1 deblurring model of frame domain[J]. Journal of Image and Graphics, 2015, 20(11)： 1434-1445. (in Chinese)

[4] 李旭超, 邊素軒, 李玉葉, 等. 圖像恢復(fù)中的凸能量泛函正則化模型綜述[J]. 中國圖象圖形學(xué)報, 2016, 21(4)： 405-415.

Li Xuchao, Bian Suxuan, Li Yuye. Survey on convex energy functional regularization model of image restoration[J]. Journal of Image and Graphics, 2016, 21(4)： 405-415. (in Chinese)

[5] 陳芳芳. 基于變分原理的圖像去噪研究[D]. 成都： 電子科技大學(xué), 2016.

[6] Bertsekas D P, Scientific A. Convex optimization algorithms[M]. Belmont： Athena Scientific, 2015.

[7] 李旭超, 宋博. 原始-對偶模型的牛頓迭代原理與圖像恢復(fù)[J]. 電子學(xué)報, 2015, 43(10)： 1984-1993.

Li Xuchao, Song Bo. Newton iterative principle of primal-dual model and image resroration[J]. Acta Electronica Sinica, 2015, 43(10)： 1984-1993. (in Chinese)

[8] 汪美玲, 周先春, 周林鋒, 等. 全變分耦合圖像去噪模型[J]. 通信學(xué)報, 2016, 37(4)： 182-191.

Wang Meiling, Zhou Xianchun, Zhou Linfeng. Coupling image denoising model based on total variation[J]. Journal on Communications, 2016, 37(4)： 182-191. (in Chinese)

[9] 張崢嶸, 黃麗麗, 費選, 等. 非局部TV正則化的圖像泊松去噪模型與算法[J]. 系統(tǒng)仿真學(xué)報, 2014, 26(9)： 2110-2115.

Zhang Zhengrong, Huang Lili, Fei Xuan. Image poisson denoising model and algorithm based on nonlocal TV regularization[J]. Journal of System Simulation, 2014, 26(9)： 2110-2115.

[10] Liu J, Zheng X. A block nonlocal TV method for image restoration[J]. UCLA CAM report, 2016： 16-25.

[11] Bergounioux M. Mathematical analysis of a inf-convolution model for image processing[J]. Journal of Optimization Theory and Applications, 2016, 168(1)： 1-21.

[12] Bredies K, Holler M. Regularization of linear inverse problems with total generalized variation[J]. Journal of Inverse and Ill-posed Problems, 2014, 22(6)： 871-913.

[13] Papafitsoros K, Sch?nlieb C B. A combined first and second order variational approach for image reconstruction[J]. Journal of mathematical imaging and vision, 2014, 48(2)： 308-338.

[14] Cai J F, Dong B, Shen Z. Image restoration： a wavelet frame based model for piecewise smooth functions and beyond[J]. Applied and Computational Harmonic Analysis, 2016, 41(1)： 94-138.

[15] 吳玉蓮, 馮象初. 利用平衡方法的非凸圖像修復(fù)[J]. 西安電子科技大學(xué)學(xué)報, 2014, 41(5)： 141-147.

Wu Yulian, Feng Xiangchu. Nonconvex image inpainting via balanced regularization approach [J]. Journal of Xidian University, 2014, 41(5)： 141-147. (in Chinese)

[16] Daubechies I, Han B, Ron A, et al. Framelets： MRA-based constructions of wavelet frames[J]. Applied and Computational Harmonic Analysis, 2003, 14(1)： 1-46.

[17] Skopina M A. Tight wavelet frames[C]∥Doklady Mathematics. MAIK Nauka/Interperiodica, 2008, 77(2)： 182-185.

[18] Bredies K, Kunisch K, Pock T. Total generalized variation[J]. SIAM Journal on Imaging Sciences, 2010, 3(3)： 492-526.

[19] Dong B, Zhang Y. An efficient algorithm for l0minimization in wavelet frame based image restoration[J]. Journal of Scientific Computing, 2013, 54(2-3)： 350-368.