楊興明，陳海燕，王 剛，王彬彬，趙銀平

1.合肥工業(yè)大學(xué) 計(jì)算機(jī)與信息學(xué)院，合肥 230009

2.合肥工業(yè)大學(xué) 電氣與自動(dòng)化工程學(xué)院，合肥 230009

基于UDCT系數(shù)的改進(jìn)HMT和在圖像去噪中應(yīng)用

楊興明1，陳海燕1，王 剛1，王彬彬1，趙銀平2

1.合肥工業(yè)大學(xué) 計(jì)算機(jī)與信息學(xué)院，合肥 230009

2.合肥工業(yè)大學(xué) 電氣與自動(dòng)化工程學(xué)院，合肥 230009

1 引言

圖像處理的應(yīng)用非常廣泛，例如去噪、融合、分割等[1-2]，圖像噪聲去除的基本方法有空間域的和變換域的，空間域?yàn)V波能夠有效濾除光滑區(qū)域的噪聲但容易模糊邊緣；變換域去噪的方法主要由傅里葉變換和小波變換，而小波變換由于其多分辨率和時(shí)頻局部等特性，廣泛用于圖像去噪。由兩個(gè)一維正交小波基張成的二維小波具有各向同性使其在表示圖像邊界時(shí)候不具有稀疏性，不是最優(yōu)基。

沿圖像的曲線邊緣表示是圖像表示的一個(gè)突破。自2000年以后常用的方法有：曲波變換（Curvelet），輪廓波變換（Contourlet）等。Curvelet變換在描述二次光滑空間的曲線時(shí)具有獨(dú)到的優(yōu)勢(shì)，其通過幅度窗口函數(shù)和角度窗口函數(shù)實(shí)現(xiàn)對(duì)圖像曲線的完美重構(gòu)，并且Curvelet變換的快速變換（Fast Discrete Curvelet Transform，F(xiàn)DCT）在頻率域使用FFT算法實(shí)現(xiàn)的[3]。Contourlet變換是在空間域通過濾波器組結(jié)構(gòu)實(shí)現(xiàn)的參數(shù)化Curvelet變換，故Contourlet基函數(shù)不是帶限的[4]。Truong T.Nguyen和Hervé Chauris提出均勻離散曲波變換（UDCT）結(jié)合了Curvelet變換和Contourlet變換的優(yōu)點(diǎn)：它是通過在頻域用類似Contourlet結(jié)構(gòu)的多分辨率濾波器組實(shí)現(xiàn)的；UDCT頻率域的濾波器組的構(gòu)造是滿足Curvelet基函數(shù)要求；故UDCT既具有Curvelet完整的理論基礎(chǔ)也有類似于Contourlet的易于實(shí)現(xiàn)濾波器組結(jié)構(gòu)[5]。UDCT和折疊FDCT一樣都是使用FFT算法實(shí)現(xiàn)的，但是二者的區(qū)別是：UDCT采用濾波器組結(jié)構(gòu)下的采樣來降低冗余率，而折疊FDCT是通過頻率域的對(duì)頻率進(jìn)行折疊來降低冗余率的；UDCT基函數(shù)是在均勻的整數(shù)網(wǎng)格中的，而FDCT基函數(shù)是在非均勻網(wǎng)格中；UDCT系數(shù)具有和小波系數(shù)一樣的樹結(jié)構(gòu)而折疊FDCT系數(shù)則不具有樹的結(jié)構(gòu)。

HMT能很好地描述這種具有樹結(jié)構(gòu)系數(shù)之間的關(guān)系。四叉樹建立每個(gè)父系數(shù)和其四個(gè)子系數(shù)之間的隱狀態(tài)聯(lián)系。關(guān)于Wavwlet HMT模型和Contourlet HMT模型已有大量研究并且都證明其具有良好的去噪效果[6-7]。本文將HMT模型應(yīng)用到UDCT系數(shù)上，實(shí)驗(yàn)證明：該模型具有良好的去噪效果。

2 均勻離散曲波變化（UDCT）的實(shí)現(xiàn)

均勻離散曲波變換（UDCT）是可逆多分辨率變換，是在頻率域通過類似Contourlet濾波器組結(jié)構(gòu)實(shí)現(xiàn)的。實(shí)現(xiàn)過程如圖1示。

圖1 UDCT的正變換和反變換

在頻域?qū)⒂玫筋愃艭ontourlet濾波器組的參數(shù)化窗口，可以表示為有1個(gè)低通和2N個(gè)方向高通的二維濾波器組。在不同的分辨率下通過級(jí)聯(lián)相同的濾波器組，即是UDCT頻域離散分解。定義2-D參數(shù)化窗口函數(shù)族構(gòu)成式（1）的單位分解；當(dāng)N=3時(shí)，7個(gè)窗函數(shù)定義的如式（2）7帶濾波器組。

其中u0(ω)為離散低通濾波器頻率響應(yīng)，ul(ω)為6個(gè)2-D方向?yàn)V波器組頻率響應(yīng)，l為方向個(gè)數(shù)。

在空間域，合成濾波器gl(n)和分解濾波器是一樣的，gl(n)=fl(n)，抽取比率依然是2I，在濾波器的輸出端，重構(gòu)圖像是用傅里葉逆變換得到UDCT系數(shù)的實(shí)部進(jìn)行重構(gòu)的，如圖2。

圖2 UDCT頻域完全重構(gòu)濾波器組

3 UDCT系數(shù)統(tǒng)計(jì)特性分析

研究UDCT系數(shù)的統(tǒng)計(jì)特性首先要研究系數(shù)實(shí)部和虛部邊緣分布和聯(lián)合分布。對(duì)Lena圖像進(jìn)行4層分解，取最細(xì)節(jié)的一個(gè)方向子帶，統(tǒng)計(jì)其系數(shù)的實(shí)部和虛部的概率統(tǒng)計(jì)和聯(lián)合分布如圖3示。

計(jì)算系數(shù)實(shí)部和虛部分布的峰值系數(shù)kurtosis分別為84.11和50.44，得到子帶內(nèi)的系數(shù)的實(shí)部和虛部都是呈非高斯分布的特性。UDCT系數(shù)的實(shí)部和虛部是零均值、不相關(guān)、同變性。從圖3（c）中可以看出系數(shù)的實(shí)部和虛部關(guān)于圓點(diǎn)呈中心對(duì)稱分布的，復(fù)系數(shù)的概率密度函數(shù)只和系數(shù)的大小有關(guān)。

圖3 Lena高頻子帶系數(shù)邊緣統(tǒng)計(jì)和聯(lián)合分布

同時(shí)，還需要對(duì)尺度間的系數(shù)的關(guān)系進(jìn)行研究，系數(shù)X的鄰系數(shù)集包括：其上尺度的父系數(shù)（parentPX）、同一尺度內(nèi)的鄰系數(shù)（neighborNX）和同一尺度的不同方向子帶的方向系數(shù)（cousinCX），如圖4所示。對(duì)Lena圖像進(jìn)行UDCT后得到的系數(shù)，選取一個(gè)細(xì)尺度系數(shù)X，在已知父系數(shù)情況下實(shí)部和虛部的Kurtosis分別為：4.12和3.95，所以，在已知父系數(shù)PX的條件下，系數(shù)X實(shí)部和虛部的分布可以用高斯混合分布進(jìn)行建模。

4 基于信息理論系數(shù)依賴性分析

考察系數(shù)X和其PX，NX，CX的依賴性強(qiáng)弱，可以通過計(jì)算系數(shù)之間的互信息量。X和Y的互信息表示Y從X那得到的信息的多少，當(dāng)X，Y獨(dú)立時(shí)，互信息量為0，當(dāng)X=Y時(shí)，互信息量最大為1。X，Y的互信息的定義如下：

其中，fX(x)，fY(y)分別為X，Y的邊緣概率密度函數(shù)，fXY(x，y)為X，Y聯(lián)合概率密度函數(shù)。實(shí)際中用直方圖結(jié)合熵和互信息中的方法計(jì)算I(X;Y)[8]?；バ畔⒘砍１挥糜谛〔ㄓ蚝颓ㄓ騺碓u(píng)價(jià)系數(shù)的相關(guān)性[9-10]。由于UDCT系數(shù)的實(shí)數(shù)部分和虛數(shù)部分具有相似的性質(zhì)，這里僅研究系數(shù)X的實(shí)數(shù)部分和其PX、NX、CX的互信息量，如表1所示。

表1 系數(shù)X和其鄰系數(shù)PX、NX、CX的互信息量

從表1中，可以看出：

子帶內(nèi)系數(shù)之間的相關(guān)性具有對(duì)其他相關(guān)性的主導(dǎo)性。即：X服從什么分布和NX的依賴性最大的，這和已得的小波系數(shù)和曲波系數(shù)的依賴性結(jié)果是一致的[9]。

5 UDCT系數(shù)的HMT模型

5.1 UDCT系數(shù)的HMT模型建立

由第二、第三部分討論得出的UDCT系數(shù)在已知父系數(shù)PX的條件下，系數(shù)X實(shí)部和虛部的分布可以用高斯混合分布進(jìn)行建模的條件高斯特性和系數(shù)間的依賴性關(guān)系，選擇HMT進(jìn)行建模。首先進(jìn)行UDCT進(jìn)行4層分解，每一層的方向分別為：6、6、12、12，對(duì)其建立如圖4（c）的多樹模型。

圖4（a）所示的是Wavelet系數(shù)的父子關(guān)系，每一個(gè)小波系數(shù)的子系數(shù)只能在一個(gè)子帶內(nèi)，圖4（b）中UDCT系數(shù)的子系數(shù)可以在兩個(gè)子帶中，在圖4（c）中樹的第三層中可以看出兩個(gè)子帶對(duì)應(yīng)于一個(gè)父系數(shù)。

對(duì)于在尺度j，方向k位置n的系數(shù)記為c(j，k，n)，其相應(yīng)的隱狀態(tài)為S(j，k，n)，值為L(zhǎng)，S分別代表具有此狀態(tài)的系數(shù)是服從大的方差和小的方差的分布。當(dāng)S(j，k，n)=L對(duì)應(yīng)的分布的均值和方差，當(dāng)值為S時(shí)，對(duì)應(yīng)分布的均值和方差分別為。則UDCT系數(shù)的概率密度函數(shù)為：

其中，UDCT系數(shù)為零均值的，則

P(S(j，k，n)=m)表示系數(shù)s(j，k，n)是狀態(tài)m(L,S)的概率，

5.2 模型參數(shù)的初始化和估計(jì)

圖4 Wavelet系數(shù)和UDCT系數(shù)父子關(guān)系

用EM算法對(duì)UDCT系數(shù)的HMT模型的參數(shù)估計(jì)，該算法估計(jì)的參數(shù)在圖像去噪中取得非常好的效果，但是該算法由于計(jì)算量龐大導(dǎo)致運(yùn)算時(shí)間較長(zhǎng)，不利于對(duì)時(shí)間要求比較高的應(yīng)用，針對(duì)這個(gè)問題，可以通過選擇合適的初始參數(shù)來減小運(yùn)算時(shí)間[11]。通過研究UDCT系數(shù)的延續(xù)性，發(fā)現(xiàn)其具有規(guī)律性。以圖4（c）為例，UDCT分解層數(shù)為4層，6個(gè)方向，分別建立了6棵樹，不同層，同層不同方向子帶的狀態(tài)轉(zhuǎn)移概率是不同的，即：若原來父系數(shù)是服從大方差分布的但是其子系數(shù)是符合大方差還是小方差分布的概率和其所在的分解層數(shù)有關(guān)的；UDCT系數(shù)在高頻上具有很強(qiáng)的方向性，隨著尺度的增加，小系數(shù)產(chǎn)生小系數(shù)的概率不斷增加，到一定的尺度之后，邊界被完全分開只剩下平滑的區(qū)域，則小系數(shù)產(chǎn)生小系數(shù)的概率趨向于1，邊界在任一高分辨率子空間可能存在也可能消失，不存在邊界的高分辨率子空間的UDCT系數(shù)值較小，故大系數(shù)產(chǎn)生大系數(shù)的概率逐漸趨向于1/2。在中間層，由于UDCT函數(shù)的方向性，邊界沒有完全分開，出現(xiàn)大小系數(shù)的概率是差不多的。狀態(tài)轉(zhuǎn)移概率是用來刻畫UDCT系數(shù)的尺度間的延續(xù)性的，考察的是子帶與其父子帶的狀態(tài)概率的關(guān)系。為了研究這種關(guān)系，可以用子帶和其父帶的各狀態(tài)的統(tǒng)計(jì)來揭示它們的狀態(tài)概率的關(guān)系。由于用高斯混合模型來模擬UDCT系數(shù)的“大”、“小”狀態(tài)，則根據(jù)UDCT變換的特點(diǎn)，用一個(gè)自適應(yīng)的閾值T[n1，n2]把每個(gè)子帶的UDCT系數(shù)分成大、小兩類來分別統(tǒng)計(jì)其狀態(tài)。

其中，＜N[n1，n2]＞為以n1，n2為中心的鄰域；c(j,k,a,b)為鄰域內(nèi)的UDCT系數(shù)。利用式（6）得到子帶和其父帶的二值圖像B[n1，n2]，B[n1/2，n2/2]：

式（6）可以得到每個(gè)子帶和其父帶的二值圖像即狀態(tài)分布圖，分別統(tǒng)計(jì)子帶和其父帶“大”、“小”狀態(tài)個(gè)數(shù)，利用式（7）得到狀態(tài)轉(zhuǎn)移關(guān)系。

其中，number()表示事件發(fā)生的次數(shù)；m=L，S；按子帶內(nèi)每塊是8×8，其父子帶是4×4統(tǒng)計(jì)，經(jīng)過多次實(shí)驗(yàn)得如圖5所示關(guān)系。

圖5 UDCT系數(shù)狀態(tài)轉(zhuǎn)移概率和尺度關(guān)系

其中，j為分解層數(shù)。用式（8）初始化HMT參數(shù)，足以表征狀態(tài)轉(zhuǎn)移的趨勢(shì)。UDCT系數(shù)隨著尺度的增加成指數(shù)遞減的規(guī)律，如圖6所示?；旌戏讲畹闹狄卜磻?yīng)了UDCT系數(shù)的指數(shù)衰減性，故可以利用這種規(guī)律性初始化方差為：

圖6 UDCT系數(shù)隨尺度衰減性

用UDCT系數(shù)方差和狀態(tài)轉(zhuǎn)移的規(guī)律對(duì)參數(shù)估計(jì)初始化后用EM算法做參數(shù)估計(jì)。EM算法過程如下（6通道同時(shí)進(jìn)行）：

（2）計(jì)算隱狀態(tài)變量的聯(lián)合概率密度函數(shù)f(c(j，k，n))；

（3）更新模型參數(shù)θl+1=argmax[lnf(c(j，k，n))|c(j，k，n)，θl]，用初始化狀態(tài)轉(zhuǎn)移矩陣減少迭代次數(shù)；

（4）l=l+1，若滿足收斂條件則停止，否則轉(zhuǎn)（2）。

6 圖像去噪

圖像去噪問題一般被描述為如式（9）的問題。

y為去噪后的圖像UDCT系數(shù)；x為含噪圖像的UDCT系數(shù)；w為噪聲的UDCT系數(shù)；x為已知的，w為通過隨機(jī)噪聲計(jì)算出來的。具體算法步驟為：

（1）對(duì)含噪圖像進(jìn)行4層分解，每層的方向分別為6，6，12，12，得到含噪系數(shù)x。

（2）對(duì)含噪系數(shù)進(jìn)行HMT建模，用EM算法結(jié)合提出的改進(jìn)參數(shù)估計(jì)方法得到參數(shù)集θv，在實(shí)驗(yàn)中取cs=2，cl=5，ps1=(0.5，0.5)，ui，m=0。

（3）利用式（10）和模型參數(shù)θv得到去噪之后的參數(shù)θu：

圖7 Lena加噪圖像不同去噪結(jié)果的對(duì)比

為Monte-Carlo方法得到的在j尺度，k方向，n位置的噪聲方差。

（4）對(duì)于給定狀態(tài)的系數(shù)的分布是符合高斯分布的，求解過程可以簡(jiǎn)化為：

（5）對(duì)已求得去噪圖像的UDCT系數(shù)進(jìn)行重構(gòu)就得到去噪之后圖像。

通過這種初始化的狀態(tài)矩陣在保證去噪精度的同時(shí)，訓(xùn)練序列的時(shí)間減少了2/5，證明該初始化是有意義的。表2所示為改進(jìn)后的HMT模型算法用到圖像序列訓(xùn)練和去噪的時(shí)間和其他算法時(shí)間的對(duì)比。實(shí)驗(yàn)圖片均采用lena512。當(dāng)噪聲是加性噪聲Wiener2函數(shù)進(jìn)行濾波效果較線性濾波器是最好的，算法簡(jiǎn)單，但不適用于非平穩(wěn)信號(hào)濾波。

表2 UDCT IHMT（優(yōu)化的）與各種去噪方法時(shí)間對(duì)比 s

對(duì)去噪效果的評(píng)價(jià)選擇峰值信噪比（PSNR）和圖像重構(gòu)以后的結(jié)構(gòu)相似性（SSIM）[12]；對(duì)lena512進(jìn)行4層分解，分別用Wiener2、Wavelet HMT[9]、Contourlet HMT[7]、UDCT IHMT四種不同的去噪方法進(jìn)行對(duì)比。表3對(duì)應(yīng)的PSNR的值，表4是對(duì)應(yīng)的結(jié)構(gòu)相似性（SSIM）。

圖7是lena512加標(biāo)準(zhǔn)差為40的噪聲和各種去噪以后的圖像。選擇加高斯白噪聲的標(biāo)準(zhǔn)差為40，從左到右依次分別為加噪、Wiener2、Wavelet HMT、Contourlet HMT、UDCT IHMT去噪圖像。

表3 UDCT IHMT與各種去噪方法的PSNR值 dB

表4 UDCT IHMT與各種去噪方法的SSIM對(duì)比

7 結(jié)論

通過研究UDCT系數(shù)X與其父系數(shù)PX，鄰系數(shù)NX，方向系數(shù)CX的條件分布，根據(jù)互信息量表征各系數(shù)之間的聯(lián)系，用HMT模型來模擬UDCT系數(shù)之間的關(guān)系，提出初始化方差和狀態(tài)轉(zhuǎn)移矩陣的方法，保證模型精確性同時(shí)節(jié)省數(shù)據(jù)訓(xùn)練時(shí)間。通過對(duì)比不同的去噪方法證明：雖然算法簡(jiǎn)單性不如Wiener2，但UDCT的IHMT模型在圖像去噪效果方面優(yōu)于其他方法。

[1]劉仁金，高遠(yuǎn)飊，郝祥根.文本圖像頁(yè)面分割算法研究[J].中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2010，40（5）：500-504.

[2]汪一休.一種交互式圖像分割的修正優(yōu)化算法[J].中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2010，40（2）：129-132.

[3]Candes E J，Donoho D L，Ying L.Fast discrete curvelet transform[J].Multiscale Modeling and Simulation，2006，5（3）：861-899.

[4]Do M N，Vetterli M.The contourlet transform：an effect directional multi-resolution image representation[J].IEEE Trans on Image Processing，2005，14：2091-2106.

[5]Nguyen T T，Chauris H.Uniform discrete curvelet transform[J]. IEEE Transactions on Signal Processing，2010，58（7）：3618-3633.

[6]Crouse M，Nowak R，Baraniuk R.Wavelet-based statistical signal processing using hidden Markov models[J].IEEE Trans on Signal Processing，1998，46：886-902.

[7]Po D D Y，Do M N.Directional multi-scale modeling of images using the contourlet transform[J].IEEE Trans on Image Processing，2006，15：1610-1620.

[8]Moddemeijer R.On estimation of entropy and mutual information of continuous distributions[J].Signal Processing，1989，16（3）：233-246.

[9]Liu J，Moulin P.Information-theoretic analysis of inter-scale dependencies between image wavelet coefficients[J].IEEE Trans on Image Processing，2001，10：1647-1658.

[10]Alecu A，Munteanu A，Pizurical A.Analysis of the statistical dependencies in the curvelet domain and applications in image compression[C]//Lecture Notes in Computer Science，2007，4678：1061-1071.

[11]汪西莉，劉芳，焦李成.一種小波域HMT模型參數(shù)初始化方法[J].計(jì)算機(jī)科學(xué)，2003，30（1）：85-87.

[12]Wang Z，Bovik A C，Sheikh H R，et al.Image quality assessment：from error visibility to structural similarity[J].Image Processing，2004，13：600-612.

YANG Xingming1,CHEN Haiyan1,WANG Gang1,WANG Binbin1,ZHAO Yinping2

1.School of Computer and Information Science,Hefei University of Technology,Hefei 230009,China
2.School of Electrical and Automation Engineering,Hefei University of Technology,Hefei 230009,China

Based on the statistical properties of coefficients of the Uniform Discrete Curvelet Transform（UDCT）,and the analysis of correlation metric mutual information about the coefficients,this paper chooses the Hidden Markov Tree to model the coefficients finally and trains the sequence with the EM algorithm.With amount of time consuming,an optimization EM algorithm based on HMT of UDCT coefficients is presented;it further optimizes the algorithm by defining the variance and state transition matrix based on the attenuation of coefficients and continuity between the scales.Experimental results show that,in the use of similarity and Peak Signal to Noise Ratio effect as the measurement of image de-noising,under the same conditions,the algorithm proposed has better real-time and de-noising effect than the Wavelet HMT,Contourlet HMT,UDCT HMT algorithm.

Uniform Discrete Curvelet Transform（UDCT）;mutual information;Hidden Markov Tree model（HMT）;Expectation-Maximization（EM）algorithm;image denoising

通過對(duì)均勻離散曲波變換（Uniform Discrete Curvelet Transform，UDCT）系數(shù)的統(tǒng)計(jì)特性研究，同時(shí)對(duì)系數(shù)相關(guān)性度量指標(biāo)互信息量的分析，最終選擇隱馬爾可夫樹模型對(duì)其系數(shù)建模，且用EM算法訓(xùn)練序列；針對(duì)訓(xùn)練時(shí)間過長(zhǎng)問題，通過分析系數(shù)的衰減性和尺度間系數(shù)延續(xù)性，提出一種新的對(duì)算法參數(shù)初值的方差和狀態(tài)轉(zhuǎn)移矩陣的優(yōu)化方法，實(shí)驗(yàn)結(jié)果證明，在采用峰值信噪比和相似度作為圖像去噪效果的度量時(shí)，同等條件下文中提出的算法比Wavelet HMT、Contourlet HMT、UDCT HMT算法有較好的實(shí)時(shí)性和去噪效果。

均勻離散曲波變換；互信息；隱馬爾可夫樹模型（HMT）；最大期望（EM）算法；圖像去噪

A

TP751.1

10.3778/j.issn.1002-8331.1112-0463

YANG Xingming,CHEN Haiyan,WANG Gang,et al.Improvement of HMT based on uniform discrete curvelet coefficients and application in image denoising.Computer Engineering and Applications,2013,49（18）：195-199.

安徽省2009年度自然科學(xué)基金資助（No.090412041）。

楊興明（1977—），男，博士，副教授，研究方向：小波分析處理，自動(dòng)控制。E-mail：xmyang@hfut.edu.cn

2011-12-23

2012-02-27

1002-8331（2013）18-0195-05

CNKI出版日期：2012-06-05 http://www.cnki.net/kcms/detail/11.2127.TP.20120605.1519.001.html