尚翠娟 林 勇 周 青

（山東大學(xué)控制科學(xué)與工程學(xué)院生物醫(yī)學(xué)工程研究所，濟(jì)南 250061）

引言

眼底圖像是眼科疾病及其他疾病的有效診斷依據(jù)，對(duì)現(xiàn)代醫(yī)學(xué)有著很重要的價(jià)值。在臨床上，視網(wǎng)膜血管網(wǎng)絡(luò)對(duì)高血壓、糖尿病、動(dòng)脈硬化、腎炎等疾病的診斷治療有重要意義［1］。眼底圖像多由眼底相機(jī)獲得，由于眼睛結(jié)構(gòu)的特殊性以及客觀條件的限制，眼底圖像的質(zhì)量往往令人不是很滿意，主要表現(xiàn)為中央視神經(jīng)乳頭灰度較高，周邊較暗，圖像背景亮度不均且血管對(duì)比度較差，整體圖像不夠明快等。這樣就使得醫(yī)生很難從圖中看清眼底的各種信息，給診斷帶來(lái)了不便。為了改善眼底圖像的視覺(jué)效果，便于醫(yī)生從眼底圖像中獲得各種必要信息，同時(shí)也為了眼底圖像的后續(xù)處理，便于提取視網(wǎng)膜血管網(wǎng)絡(luò)，需要對(duì)眼底圖像進(jìn)行必要的增強(qiáng)處理。

近年來(lái)，小波變換在圖像增強(qiáng)領(lǐng)域得到廣泛應(yīng)用，相比于傳統(tǒng)的增強(qiáng)方法如直接調(diào)整圖像對(duì)比度、直方圖均衡化等，能更好地增強(qiáng)圖像的細(xì)節(jié)特征，并取得了不錯(cuò)的效果。但是由于小波變換不能很好地刻畫(huà)噪聲和邊緣，對(duì)于一些存在噪聲的圖像，基于小波變換的增強(qiáng)算法在增強(qiáng)邊緣的同時(shí)明顯地?cái)U(kuò)大了噪聲?；谛〔ǖ牟蛔?，多尺度幾何分析被提了出來(lái)，并迅速應(yīng)用于圖像處理的各個(gè)領(lǐng)域。多尺度幾何分析中的 Contourlet變換繼承了Curvelet變換的各項(xiàng)異性尺度關(guān)系，實(shí)現(xiàn)了“真正的圖像二維表示方法”，是圖像的一種稀疏表示方法，其在圖像處理方面表現(xiàn)出了很大的優(yōu)勢(shì)［2～4］。在圖像去噪聲增強(qiáng)領(lǐng)域，Contourlet變換有一個(gè)很大的缺點(diǎn)，即不具有平移不變性，而平移不變性在圖像去噪中是一個(gè)比較重要的參數(shù)［5］，因此 M.N.Do等人于2006年提出了非下采樣 Contourlet變換，并分析了其在增強(qiáng)去噪方面的應(yīng)用［6］。

考慮到眼底圖像具有背景亮度不均、血管對(duì)比度較差、微小細(xì)節(jié)繁多的特點(diǎn)［7］，而非下采樣Contourlet變換能夠很好地刻畫(huà)圖像邊緣的紋理信息且具有平移不變性，所以該變換適合用來(lái)分析視網(wǎng)膜圖像中的病變和血管信息。圖像去噪增強(qiáng)算法的關(guān)鍵在于如何區(qū)分噪聲和紋理信息，尤其是弱紋理信息。使用閾值方法進(jìn)行去噪，比較簡(jiǎn)潔有效，因此大多數(shù)去噪算法的重點(diǎn)都在于如何選取區(qū)分噪聲和弱紋理信息的閾值［8］。閾值的選取多依賴于對(duì)圖像噪聲方差的估計(jì)，而且多數(shù)的研究都以高斯模型作為噪聲估計(jì)模型。目前比較常見(jiàn)的是Donoho提出的全局閾值和最近研究較多的貝葉斯萎縮閾值，這兩種方法的閾值估計(jì)過(guò)程也用到了噪聲方差。實(shí)際上，通過(guò)固定模型對(duì)存在很大隨機(jī)性的噪聲進(jìn)行估計(jì)，其結(jié)果必然存在著偏差。為了避免對(duì)噪聲方差的估計(jì)，本研究提出了基于非下采樣Contourlet變換的主分量分析方法，并將其應(yīng)用于眼底圖像的增強(qiáng)處理上。

1 基本理論

1.1 非下采樣Contourlet變換（NSCT）

NSCT是去除了下采樣環(huán)節(jié)的Contourlet變換，由非下采樣金字塔分解和多方向子帶分解兩部分組成［6］。首先非下采樣金字塔分解基于二維兩通道濾波器組，將圖像分解成大小與原圖像大小一樣的低通子帶和高通子帶。高通子帶經(jīng)方向?yàn)V波器組被分解成大小不變的多個(gè)方向子帶。對(duì)低通子帶重復(fù)上述操作就可實(shí)現(xiàn)圖像的多尺度多方向分解。圖1是NSCT的濾波器組結(jié)構(gòu)和頻率分解圖。

圖1 NSCT。（a）濾波器組結(jié)構(gòu)圖；（b）頻率分解圖Fig.1 NSCT.（a）Construction of the filter banks；（b）Resulting of frequency division

圖像增強(qiáng)方法的關(guān)鍵在于如何區(qū)分噪聲和邊緣信息，尤其是弱邊緣信息。在頻域里，弱邊緣和噪聲都表現(xiàn)為幅度小的系數(shù)。噪聲有很大的隨機(jī)性因此不具有幾何結(jié)構(gòu)，而弱邊緣信息具有幾何結(jié)構(gòu)，這是 NSCT區(qū)分噪聲和弱邊緣信息的基礎(chǔ)。NSCT具有平移不變性，圖像經(jīng)NSCT變換后子帶的每一個(gè)像素和原圖像的像素一一對(duì)應(yīng)，因此可從每一個(gè)NSCT系數(shù)中收集像素的幾何信息。通過(guò)分析不同子帶的系數(shù)分布，可以把像素分成三大類：強(qiáng)邊緣，弱邊緣和噪聲。強(qiáng)邊緣對(duì)應(yīng)著同一尺度所有方向子帶中具有大幅值系數(shù)的像素，弱邊緣對(duì)應(yīng)于同一尺度下在一部分方向子帶中有大幅值但在另一部分方向子帶中有小幅值的像素，噪聲對(duì)應(yīng)于在所有方向子帶中有小幅值的像素。

一個(gè)簡(jiǎn)單的像素分類方法就是計(jì)算每一個(gè)像素同一尺度下方向子帶中系數(shù)幅度的均值和最大值，然后按式（1）分類：

式中，T是估計(jì)閾值，c是一個(gè)閾值調(diào)節(jié)參數(shù)。

圖像增強(qiáng)是為了抑制噪聲和增強(qiáng)弱邊緣，根據(jù)上面的像素分類，通過(guò)式（2）中的匹配函數(shù)對(duì)NSCT系數(shù)修正，達(dá)到抑制噪聲和增強(qiáng)弱邊緣的效果。

式中，x是方向子帶中的像素值，p為增強(qiáng)調(diào)節(jié)參數(shù)，范圍取在0到1之間。

1.2 主分量分析（PCA）

PCA是一種基于統(tǒng)計(jì)學(xué)的特征提取方法，近年來(lái)在圖像分析和模式識(shí)別領(lǐng)域應(yīng)用比較廣泛［10－11］。PCA方法的思想是在 K－L變換基礎(chǔ)上建立起來(lái)的，根據(jù)K－L變換在測(cè)量空間中找到一組正交向量，將模式矢量從原來(lái)的n維空間投影到這組正交矢量張成的m維子空間上（通常m?n），其投影系數(shù)構(gòu)成原模式的特征矢量，這組向量能最大的表示出數(shù)據(jù)的方差。PCA有多種不同的計(jì)算方法，本研究是通過(guò)計(jì)算數(shù)據(jù)矩陣X的協(xié)方差矩陣COV的特征值和特征向量，得到正交變換矩陣W進(jìn)行求解的。

根據(jù)矩陣分析理論，令M×N階的圖像矩陣表示為X=［x1，x2，…，xn］′，其協(xié)方差矩陣為

式中，

協(xié)方差矩陣的特征值分解為

式中，U=［μ1，μ2，…，μN(yùn)］是 COV 的特征向量構(gòu)成的正交矩陣，V是COV的特征值構(gòu)成的對(duì)角矩陣，diag（V）= ［λ1，λ2，…，λN］。當(dāng)特征值 λN按從大到小的順序排列時(shí)，那么所對(duì)應(yīng)U的各個(gè)基向量就是PCA的最優(yōu)投影方向，按該方向?qū)?shù)據(jù)進(jìn)行投影，得到的各主分量互不相關(guān)。

根據(jù)矩陣分析原理，特征值的和等于協(xié)方差矩陣的跡，而協(xié)方差矩陣對(duì)角線上的元素是原各變量的方差，考慮特征值中前d個(gè)較大的特征值，定義方差貢獻(xiàn)率為：

當(dāng)方差貢獻(xiàn)率η足夠大（比如達(dá)到70％、80％或 90％），那么前d個(gè)（d?N）基向量W= ［μ1，μ2，…，μd］可以表征xi的主要特征，這樣就可以應(yīng)用前d個(gè)主分量來(lái)估計(jì)x的重構(gòu)模型：

按式（7）完成對(duì)M個(gè)行向量的重建，得到原圖像矩陣X?的重建矩陣。

由于方向子帶中的Contourlet系數(shù)多數(shù)是噪聲引起的系數(shù)幅值，所以通過(guò)PCA方法重構(gòu)的矩陣反映了由噪聲引起的Contourlet系數(shù)的幅值特征。使用重構(gòu)矩陣絕對(duì)值的中值作為對(duì)該方向子帶中由噪聲引起的Contourlet系數(shù)幅值的估計(jì)。

2 算法的實(shí)現(xiàn)步驟

步驟1：對(duì)眼底圖像進(jìn)行NSCT分解。

步驟2：對(duì)于某一尺度下的一個(gè)方向子帶。

a.根據(jù)式（3）～式（7）對(duì)分解后各方向子帶進(jìn)行重構(gòu)，根據(jù)式（8）計(jì)算出每一個(gè)方向子帶的噪聲估計(jì)值。

b.根據(jù)式（1）對(duì)各個(gè)方向子帶的像素進(jìn)行分類，然后根據(jù)式（2）對(duì)系數(shù)進(jìn)行修正處理。

步驟3：按步驟2逐一處理所有尺度下的每一個(gè)方向子帶。

步驟4：對(duì)處理后的NSCT系數(shù)進(jìn)行重建。

3 算法的驗(yàn)證

為了驗(yàn)證算法的正確性和有效性，選擇多幅眼底圖像進(jìn)行測(cè)試。圖像來(lái)自公開(kāi)的DRIVE數(shù)據(jù)庫(kù)，圖像大小為768像素×584像素，RGB彩色圖像，以TIFF格式保存。因?yàn)閳D像的綠色分量對(duì)比度較高，所以測(cè)試時(shí)僅取用了綠色分量部分。

實(shí)驗(yàn)中對(duì)基于小波的Donoho軟閾值去噪（WSoft）、基于小波的 PCA去噪（W-PCA）、基于 NSCT的貝葉斯去噪（NSCT-B）、基于NSCT的PCA（NSCTPCA）去噪四種方法進(jìn)行了比較。本文使用峰值信噪比（PSNR）作為評(píng)論標(biāo)準(zhǔn)，定義為

式中，Iij是表示原始圖像I中位置在（i，j）的像素的灰度值，原始圖像I是一個(gè)理想的不含噪聲圖像；理想圖像經(jīng)噪聲干擾后即為實(shí)際獲取的圖像X，即待處理的圖像；Yij是經(jīng)處理后得到的圖像Y中相應(yīng)位置上的灰度值；M×N是圖像的大小，Imax是原始圖像中灰度值的最大值，灰度圖像中Imax一般取為255。

理論上，NSCT分解的層次越高、方向子帶分解的個(gè)數(shù)越多，對(duì)圖像的刻畫(huà)越細(xì)致，去噪增強(qiáng)效果越好，但同時(shí)也制約了算法的速度。分解層次每增加一層，用時(shí)增加一倍左右；在相同的分解層次下，方向子帶的個(gè)數(shù)每增加一倍，用時(shí)增加一倍或一倍以上。同時(shí)，圖像尺寸的大小對(duì)速度也有較大的影響，圖像尺寸縮小一倍，而用時(shí)可縮短四倍左右。因此本文對(duì)圖像做子塊劃分，每個(gè)子塊的大小為76×76，分別對(duì)子塊進(jìn)行去噪處理。為了保證本文算法的去噪效果，并達(dá)到可以和小波去噪算法相比擬的速度，實(shí)驗(yàn)對(duì)眼底圖像進(jìn)行NSCT分解的尺度為［2、3、4］。

PCA部分主分量個(gè)數(shù)的選取由方差貢獻(xiàn)率調(diào)整，NSCT是對(duì)圖像的稀疏表示，分解后各方向子帶中的系數(shù)主要為噪聲，故實(shí)驗(yàn)中方差貢獻(xiàn)率保持在98.5％ ～99.5％。經(jīng)過(guò)實(shí)驗(yàn)測(cè)試，式（2）中的參數(shù)p對(duì)PSNR的影響程度較小，為了得到最佳的 PSNR，最終選取為p=0.6。式（1）中的參數(shù)c對(duì)實(shí)驗(yàn)結(jié)果影響較大，c值的選取同各層次方向子帶的分解個(gè)數(shù)有較大的關(guān)系，因?yàn)榉较蜃訋У膫€(gè)數(shù)越多，方向子帶對(duì)圖像的表示越稀疏，子帶內(nèi)的細(xì)節(jié)分量越少，噪聲越多，c參數(shù)取較大值效果較好；反之，方向子帶的個(gè)數(shù)少，c參數(shù)取較小值效果較好?；谛〔ㄗ儞Q的去噪算法，對(duì)圖像進(jìn)行3層分解。

選取一幅質(zhì)量較高的眼底圖像作為理想圖像，對(duì)其加入高斯噪聲，然后對(duì)噪聲圖像進(jìn)行處理。

4 結(jié)果

圖2顯示了眼底圖像原圖和加入高斯噪聲后的圖像，圖3～圖5分別顯示了規(guī)一化噪聲方差sigma=0.04～0.08時(shí)，各種算法處理后的圖像，為了便于觀察，噪聲圖像和處理后的圖像僅顯示了原圖方框區(qū)域內(nèi)部分。

圖2 原圖及含噪圖像。（a）原圖；（b）含噪圖像：sigma=0.08；（c）含噪圖像：sigma=0.06；（d）含噪圖像：sigma=0.04Fig.2 Originalimage and noise images.（a）Original image（b）noise image：sigma=0.08；（c）noise image：sigma=0.06；（d）noise image：sigma=0.04

表1為不同噪聲情況下對(duì)噪聲圖像進(jìn)行各種方法處理后的PSNR值。通過(guò)比較可知：基于PCA方法估計(jì)噪聲能量進(jìn)行的去噪增強(qiáng)處理要優(yōu)于基于估計(jì)噪聲方差的去噪增強(qiáng)處理；基于NSCT的去噪增強(qiáng)處理在噪聲較大的情況下，要優(yōu)于基于小波變換的去噪增強(qiáng)處理；而本文提出的算法在各種噪聲情況下，PSNR值都高于其他方法，視覺(jué)效果也都優(yōu)于其他算法。

對(duì)DRIVE圖片庫(kù)中的20幅測(cè)試圖像隨意選取10幅進(jìn)行NSCT-PCA增強(qiáng)處理，由專業(yè)眼科醫(yī)師和圖像處理人員進(jìn)行主觀評(píng)價(jià)，得出結(jié)論：圖像質(zhì)量較好時(shí)，對(duì)圖像細(xì)節(jié)的增強(qiáng)效果比較明顯；圖像含噪較多時(shí)，對(duì)圖像的整體去噪效果明顯。本文提出的算法能夠增強(qiáng)眼底圖像的紋理細(xì)節(jié)，提高了血管對(duì)比度，尤其在圖像周邊較暗區(qū)域的微小血管也得到了增強(qiáng)。

圖3 sigma=0.08時(shí)的去噪結(jié)果。（a）W-Soft；（b）W-PCA；（c）NSCT-B；（d）NSCT-PCAFig.3 Images after ehancing by the methords with sigma=0.08.（a）W-Soft；（b）W-PCA；（c）NSCT-B；（d）NSCT-PCA

圖4 sigma=0.06時(shí)的去噪結(jié)果。（a）W-Soft；（b）W-PCA；（c）NSCT-B；（d）NSCT-PCAFig.4 Images after ehancing by the methords with sigma=0.06.（a）W-Soft；（b）W-PCA；（c）NSCT-B；（d）NSCT-PCA

圖5 sigma=0.04時(shí)的去噪結(jié)果。（a）W-Soft；（b）W-PCA；（c）NSCT-B；（d）NSCT-PCAFig.5 images after ehancing by the methords with sigma=0.04.（a）W-Soft；（b）W-PCA；（c）NSCT-B；（d）NSCT-PCA

表1 噪聲圖像經(jīng)各種方法增強(qiáng)后的PSNR/dBTab.1 PSNR of noise images after enhancing

5 討論和結(jié)論

通過(guò)客觀和主觀兩方面的評(píng)價(jià)，可以看出本文提出的算法在眼底圖像去噪增強(qiáng)方面可以取得不錯(cuò)的效果。NSCT多尺度、多方向的分析特性使得NSCT對(duì)圖像細(xì)節(jié)和邊緣的刻畫(huà)更為細(xì)致，同時(shí)以主分量分析來(lái)估計(jì)去噪閾值的方法，避免了對(duì)噪聲方差的估計(jì)，去噪效果比基于小波變換的算法更加優(yōu)良。但是NSCT算法的速度和小波變換相比還有一定的差距，因此進(jìn)行去噪處理時(shí)，對(duì)圖像進(jìn)行分解的層次數(shù)及方向子帶分解的個(gè)數(shù)不宜太大。本研究采用了對(duì)圖像進(jìn)行子塊劃分的方法來(lái)縮短算法的用時(shí)，這使得算法的速度可以和基于小波的算法相比擬，但若圖像尺寸過(guò)大，此方法的效果就不再明顯。

由于非下采樣Contourlet變換具有多尺度、多方向的分析特性，在對(duì)圖像細(xì)節(jié)和邊緣特征的刻畫(huà)上具有很大的優(yōu)勢(shì)。并且NSCT還具有圖像去噪增強(qiáng)領(lǐng)域里比較重要的特性——平移不變性，因此被認(rèn)為在圖像增強(qiáng)領(lǐng)域有很大的發(fā)展?jié)摿Α１狙芯刻岢隽嘶贜SCT的主分量分析去噪方法，該方法充分利用了NSCT的優(yōu)勢(shì)，并通過(guò)PCA方法避免了比較困難的噪聲方差估計(jì)，并首次應(yīng)用于眼底圖像的處理上。實(shí)驗(yàn)表明，該方法能夠很好地改善眼底圖像的視覺(jué)效果，有效地提高眼底圖像的對(duì)比度，圖像周邊較暗區(qū)域的微小血管也可以得到增強(qiáng)。

［1］Noronha K，Nayak J，Bhat SN.Enhancement of retinal fundus Image to highlight the features for detection of abnormal eyes［A］.In：Luk KM，eds.TENCON 2006.2006 IEEE Region 10 Conference［C］.Hong Kong：IEEE Hong Kong Section，2006.1－4.

［2］Do MN，Vetterli M.Contourlets：a directional multiresolution image representation［A］.In：Reibman A，Knox K，eds.Proceedings of IEEE International Conference on Image Processing［C］.Rochester：IEEE Rochester Section，2002.357－360.

［3］Do MN，Vetterli M.The contourlet transform：an efficient directional multiresolution image representation ［J］.IEEE Transactions on Image Processing，2005，14（2）：2091－2106.

［4］Po DDY，Do MN.Directional multiscale modeling of images using the contourlet transform［J］.IEEE Transactions on Image Processing，2006，15（6）：1610－1620.

［5］Cunha AL， ZhouJianping， DoMN.Thenonsubsampled contourlet transform：theory，design，and application ［J］.IEEE Transactions on Image Processing，2006，15（10）：3089－3101.

［6］Zhou Jianping，Cunha AL，Do MN.Nonsubsampled contourlet transform：construction and application in enhancement［A］.In：Marsh A，eds. Proceedings of IEEE International Conference on Image Processing［C］.Genoa：IEEE Genoa Section，2005.469－472.

［7］Feng Peng，Pan Yingjun，Wei Biao，et al.Enhancing retinal image by the Contourlet transform ［J］.Medical Engineering，2004，17（6）：452－456.

［8］Donoho DL. De-noising by soft-thresholding ［J］.IEEE Transactions on Information Theory，1995，41（3）：613－627.

［9］劉盛鵬，方勇.基于貝葉斯估計(jì)的 Contourlet域圖像降噪方法［J］.計(jì)算機(jī)工程，2007，33（18）：31－33.

［10］張久文，李建征，孟令鋒.基于Contourlet的圖像PCA去噪方法［J］.計(jì)算機(jī)工程與應(yīng)用，2007，43（21）：46－48.

［11］芮挺，王金巖，沈春林，等.基于 PCA的圖像小波去噪方法［J］.小型微型計(jì)算機(jī)系統(tǒng)，2006，27（1）：158－161.