基于非下采样Contourlet变换的眼底图像PCA增强方法

尚翠娟 林 勇 周 青

（山东大学控制科学与工程学院生物医学工程研究所，济南 250061）

引言

眼底图像是眼科疾病及其他疾病的有效诊断依据，对现代医学有着很重要的价值。在临床上，视网膜血管网络对高血压、糖尿病、动脉硬化、肾炎等疾病的诊断治疗有重要意义［1］。眼底图像多由眼底相机获得，由于眼睛结构的特殊性以及客观条件的限制，眼底图像的质量往往令人不是很满意，主要表现为中央视神经乳头灰度较高，周边较暗，图像背景亮度不均且血管对比度较差，整体图像不够明快等。这样就使得医生很难从图中看清眼底的各种信息，给诊断带来了不便。为了改善眼底图像的视觉效果，便于医生从眼底图像中获得各种必要信息，同时也为了眼底图像的后续处理，便于提取视网膜血管网络，需要对眼底图像进行必要的增强处理。

近年来，小波变换在图像增强领域得到广泛应用，相比于传统的增强方法如直接调整图像对比度、直方图均衡化等，能更好地增强图像的细节特征，并取得了不错的效果。但是由于小波变换不能很好地刻画噪声和边缘，对于一些存在噪声的图像，基于小波变换的增强算法在增强边缘的同时明显地扩大了噪声。基于小波的不足，多尺度几何分析被提了出来，并迅速应用于图像处理的各个领域。多尺度几何分析中的 Contourlet变换继承了Curvelet变换的各项异性尺度关系，实现了“真正的图像二维表示方法”，是图像的一种稀疏表示方法，其在图像处理方面表现出了很大的优势［2～4］。在图像去噪声增强领域，Contourlet变换有一个很大的缺点，即不具有平移不变性，而平移不变性在图像去噪中是一个比较重要的参数［5］，因此 M.N.Do等人于2006年提出了非下采样 Contourlet变换，并分析了其在增强去噪方面的应用［6］。

考虑到眼底图像具有背景亮度不均、血管对比度较差、微小细节繁多的特点［7］，而非下采样Contourlet变换能够很好地刻画图像边缘的纹理信息且具有平移不变性，所以该变换适合用来分析视网膜图像中的病变和血管信息。图像去噪增强算法的关键在于如何区分噪声和纹理信息，尤其是弱纹理信息。使用阈值方法进行去噪，比较简洁有效，因此大多数去噪算法的重点都在于如何选取区分噪声和弱纹理信息的阈值［8］。阈值的选取多依赖于对图像噪声方差的估计，而且多数的研究都以高斯模型作为噪声估计模型。目前比较常见的是Donoho提出的全局阈值和最近研究较多的贝叶斯萎缩阈值，这两种方法的阈值估计过程也用到了噪声方差。实际上，通过固定模型对存在很大随机性的噪声进行估计，其结果必然存在着偏差。为了避免对噪声方差的估计，本研究提出了基于非下采样Contourlet变换的主分量分析方法，并将其应用于眼底图像的增强处理上。

1 基本理论

1.1 非下采样Contourlet变换（NSCT）

NSCT是去除了下采样环节的Contourlet变换，由非下采样金字塔分解和多方向子带分解两部分组成［6］。首先非下采样金字塔分解基于二维两通道滤波器组，将图像分解成大小与原图像大小一样的低通子带和高通子带。高通子带经方向滤波器组被分解成大小不变的多个方向子带。对低通子带重复上述操作就可实现图像的多尺度多方向分解。图1是NSCT的滤波器组结构和频率分解图。

图1 NSCT。（a）滤波器组结构图；（b）频率分解图Fig.1 NSCT.（a）Construction of the filter banks；（b）Resulting of frequency division

图像增强方法的关键在于如何区分噪声和边缘信息，尤其是弱边缘信息。在频域里，弱边缘和噪声都表现为幅度小的系数。噪声有很大的随机性因此不具有几何结构，而弱边缘信息具有几何结构，这是 NSCT区分噪声和弱边缘信息的基础。NSCT具有平移不变性，图像经NSCT变换后子带的每一个像素和原图像的像素一一对应，因此可从每一个NSCT系数中收集像素的几何信息。通过分析不同子带的系数分布，可以把像素分成三大类：强边缘，弱边缘和噪声。强边缘对应着同一尺度所有方向子带中具有大幅值系数的像素，弱边缘对应于同一尺度下在一部分方向子带中有大幅值但在另一部分方向子带中有小幅值的像素，噪声对应于在所有方向子带中有小幅值的像素。

一个简单的像素分类方法就是计算每一个像素同一尺度下方向子带中系数幅度的均值和最大值，然后按式（1）分类：

式中，T是估计阈值，c是一个阈值调节参数。

图像增强是为了抑制噪声和增强弱边缘，根据上面的像素分类，通过式（2）中的匹配函数对NSCT系数修正，达到抑制噪声和增强弱边缘的效果。

式中，x是方向子带中的像素值，p为增强调节参数，范围取在0到1之间。

1.2 主分量分析（PCA）

PCA是一种基于统计学的特征提取方法，近年来在图像分析和模式识别领域应用比较广泛［10－11］。PCA方法的思想是在 K－L变换基础上建立起来的，根据K－L变换在测量空间中找到一组正交向量，将模式矢量从原来的n维空间投影到这组正交矢量张成的m维子空间上（通常m≪n），其投影系数构成原模式的特征矢量，这组向量能最大的表示出数据的方差。PCA有多种不同的计算方法，本研究是通过计算数据矩阵X的协方差矩阵COV的特征值和特征向量，得到正交变换矩阵W进行求解的。

根据矩阵分析理论，令M×N阶的图像矩阵表示为X=［x1，x2，…，xn］′，其协方差矩阵为

式中，

协方差矩阵的特征值分解为

式中，U=［μ1，μ2，…，μN］是 COV 的特征向量构成的正交矩阵，V是COV的特征值构成的对角矩阵，diag（V）= ［λ1，λ2，…，λN］。当特征值 λN按从大到小的顺序排列时，那么所对应U的各个基向量就是PCA的最优投影方向，按该方向对数据进行投影，得到的各主分量互不相关。

根据矩阵分析原理，特征值的和等于协方差矩阵的迹，而协方差矩阵对角线上的元素是原各变量的方差，考虑特征值中前d个较大的特征值，定义方差贡献率为：

当方差贡献率η足够大（比如达到70％、80％或 90％），那么前d个（d≪N）基向量W= ［μ1，μ2，…，μd］可以表征xi的主要特征，这样就可以应用前d个主分量来估计x的重构模型：

按式（7）完成对M个行向量的重建，得到原图像矩阵X?的重建矩阵。

由于方向子带中的Contourlet系数多数是噪声引起的系数幅值，所以通过PCA方法重构的矩阵反映了由噪声引起的Contourlet系数的幅值特征。使用重构矩阵绝对值的中值作为对该方向子带中由噪声引起的Contourlet系数幅值的估计。

2 算法的实现步骤

步骤1：对眼底图像进行NSCT分解。

步骤2：对于某一尺度下的一个方向子带。

a.根据式（3）～式（7）对分解后各方向子带进行重构，根据式（8）计算出每一个方向子带的噪声估计值。

b.根据式（1）对各个方向子带的像素进行分类，然后根据式（2）对系数进行修正处理。

步骤3：按步骤2逐一处理所有尺度下的每一个方向子带。

步骤4：对处理后的NSCT系数进行重建。

3 算法的验证

为了验证算法的正确性和有效性，选择多幅眼底图像进行测试。图像来自公开的DRIVE数据库，图像大小为768像素×584像素，RGB彩色图像，以TIFF格式保存。因为图像的绿色分量对比度较高，所以测试时仅取用了绿色分量部分。

实验中对基于小波的Donoho软阈值去噪（WSoft）、基于小波的 PCA去噪（W-PCA）、基于 NSCT的贝叶斯去噪（NSCT-B）、基于NSCT的PCA（NSCTPCA）去噪四种方法进行了比较。本文使用峰值信噪比（PSNR）作为评论标准，定义为

式中，Iij是表示原始图像I中位置在（i，j）的像素的灰度值，原始图像I是一个理想的不含噪声图像；理想图像经噪声干扰后即为实际获取的图像X，即待处理的图像；Yij是经处理后得到的图像Y中相应位置上的灰度值；M×N是图像的大小，Imax是原始图像中灰度值的最大值，灰度图像中Imax一般取为255。

理论上，NSCT分解的层次越高、方向子带分解的个数越多，对图像的刻画越细致，去噪增强效果越好，但同时也制约了算法的速度。分解层次每增加一层，用时增加一倍左右；在相同的分解层次下，方向子带的个数每增加一倍，用时增加一倍或一倍以上。同时，图像尺寸的大小对速度也有较大的影响，图像尺寸缩小一倍，而用时可缩短四倍左右。因此本文对图像做子块划分，每个子块的大小为76×76，分别对子块进行去噪处理。为了保证本文算法的去噪效果，并达到可以和小波去噪算法相比拟的速度，实验对眼底图像进行NSCT分解的尺度为［2、3、4］。

PCA部分主分量个数的选取由方差贡献率调整，NSCT是对图像的稀疏表示，分解后各方向子带中的系数主要为噪声，故实验中方差贡献率保持在98.5％ ～99.5％。经过实验测试，式（2）中的参数p对PSNR的影响程度较小，为了得到最佳的 PSNR，最终选取为p=0.6。式（1）中的参数c对实验结果影响较大，c值的选取同各层次方向子带的分解个数有较大的关系，因为方向子带的个数越多，方向子带对图像的表示越稀疏，子带内的细节分量越少，噪声越多，c参数取较大值效果较好；反之，方向子带的个数少，c参数取较小值效果较好。基于小波变换的去噪算法，对图像进行3层分解。

选取一幅质量较高的眼底图像作为理想图像，对其加入高斯噪声，然后对噪声图像进行处理。

4 结果

图2显示了眼底图像原图和加入高斯噪声后的图像，图3～图5分别显示了规一化噪声方差sigma=0.04～0.08时，各种算法处理后的图像，为了便于观察，噪声图像和处理后的图像仅显示了原图方框区域内部分。

图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为不同噪声情况下对噪声图像进行各种方法处理后的PSNR值。通过比较可知：基于PCA方法估计噪声能量进行的去噪增强处理要优于基于估计噪声方差的去噪增强处理；基于NSCT的去噪增强处理在噪声较大的情况下，要优于基于小波变换的去噪增强处理；而本文提出的算法在各种噪声情况下，PSNR值都高于其他方法，视觉效果也都优于其他算法。

对DRIVE图片库中的20幅测试图像随意选取10幅进行NSCT-PCA增强处理，由专业眼科医师和图像处理人员进行主观评价，得出结论：图像质量较好时，对图像细节的增强效果比较明显；图像含噪较多时，对图像的整体去噪效果明显。本文提出的算法能够增强眼底图像的纹理细节，提高了血管对比度，尤其在图像周边较暗区域的微小血管也得到了增强。

图3 sigma=0.08时的去噪结果。（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时的去噪结果。（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时的去噪结果。（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 噪声图像经各种方法增强后的PSNR/dBTab.1 PSNR of noise images after enhancing

5 讨论和结论

通过客观和主观两方面的评价，可以看出本文提出的算法在眼底图像去噪增强方面可以取得不错的效果。NSCT多尺度、多方向的分析特性使得NSCT对图像细节和边缘的刻画更为细致，同时以主分量分析来估计去噪阈值的方法，避免了对噪声方差的估计，去噪效果比基于小波变换的算法更加优良。但是NSCT算法的速度和小波变换相比还有一定的差距，因此进行去噪处理时，对图像进行分解的层次数及方向子带分解的个数不宜太大。本研究采用了对图像进行子块划分的方法来缩短算法的用时，这使得算法的速度可以和基于小波的算法相比拟，但若图像尺寸过大，此方法的效果就不再明显。

由于非下采样Contourlet变换具有多尺度、多方向的分析特性，在对图像细节和边缘特征的刻画上具有很大的优势。并且NSCT还具有图像去噪增强领域里比较重要的特性——平移不变性，因此被认为在图像增强领域有很大的发展潜力。本研究提出了基于NSCT的主分量分析去噪方法，该方法充分利用了NSCT的优势，并通过PCA方法避免了比较困难的噪声方差估计，并首次应用于眼底图像的处理上。实验表明，该方法能够很好地改善眼底图像的视觉效果，有效地提高眼底图像的对比度，图像周边较暗区域的微小血管也可以得到增强。

［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］刘盛鹏，方勇.基于贝叶斯估计的 Contourlet域图像降噪方法［J］.计算机工程，2007，33（18）：31－33.

［10］张久文，李建征，孟令锋.基于Contourlet的图像PCA去噪方法［J］.计算机工程与应用，2007，43（21）：46－48.

［11］芮挺，王金岩，沈春林，等.基于 PCA的图像小波去噪方法［J］.小型微型计算机系统，2006，27（1）：158－161.