正交小波变换去噪文档.doc

基于正交小波变换的图像去噪的研究与实现随着计算机科学和图像处理技术的迅速发展，数字图像在医学成像、模式识别等方面取得了广泛应用。但是，现场采集的数字图像一般都包含噪声，而且有些图像的噪声非常严重，因此，需要对数字图像进行去噪处理，便于更高层次的图像分析与理解。本系统利用正交小波变换和软阈值方法对数字图像的去噪进行了研究与实现。本系统主要包含正交小波变换、软阈值去噪与小波反变换部分，其中，正交小波反变换是指对包含噪声的数字图像进行正交小波变换，得到小波系数；软阈值处理是指对小波系数进行软阈值处理，去除噪声；正交小波反变换是指对去噪后的小波系数进行正交小波反变换，得到去噪图像。此外，为了减弱图像边缘失真的问题，最后进行了滤波处理。应用本系统方法对BMP灰度噪声图像进行了处理，由于正交小波变换去噪容易造成去噪后的图像边缘出现失真，在经过正交小波去噪后又进行了均值滤波处理，经过对比实验显示，算法实现了较好的去噪效果。1系统开发环境 硬件配置：硬盘：80GB 内存：512MB操作系统： Windows XP开发工具：Visual C+ 6.02系统设计2.1 构架概述本系统主要实现利用离散小波对数字图像进行去除噪声。首先导入待处理的256色位图，然后利用加噪算法为初始图像加入随机噪声，最后利用本系统的核心算法小波去噪算法实现图像噪声的去除。系统流程如图1所示。图1 程序流程图2.2读取图像模块读写BMP图像文件可以用DIB函数来实现，但是其得到的数据不便于后面小波去噪的相关处理，所以本系统重新编写函数读取和显示BMP图像数据。利用BMP的文件格式和CFile类的知识，建立一个打开BMP位图的函数，在函数中，先利用CFile类的成员函数Open打开文件，然后利用成员函数Read将BMP位图的各个部分的数据读入内存。再修改Viusual C+中的OnDraw函数，利用SetPixel函数把读到内存中的数据显示在屏幕上。需要注意的是，由于BMP数据是从下向上排列的，所以显示是从最低行开始的。2.3加噪模块从内存中循环取得BMP图像的各点像素值，然后利用随机函数产生随机数，为各点的像素值加上加入一个随机数，最后将改动后的数据复制到内存中。2.4去噪模块2.4.1小波滤波器选择小波变换具有很大的灵活性，在理论上可以有无数个小波基可供选择，同时这也为小波变换的应用提出了一个难题，那就是如何正确选择小波基。由于小波变换是将原始图像与小波基函数以及尺度函数进行内积运算，1989 年Daubechies 基于离散滤波器迭代的方法构造了紧支集的规范正交小波基，因而内积运算转换为信号和离散滤波器的卷积运算，小波变换中的小波基的选择转换为正交镜像滤波器组的选择。Haar小波是所有正交紧支撑小波中唯一具有对称性的小波，Haar小波的支撑极短，其高通和低通滤波器均只有两拍，可以节省计算量，所以本系统选择Haar小波基作为正交小波变换的小波基。Haar小波的滤波器系数可有资料查得为：Lo_D20=0.7071,0.7071,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 Hi_D20=-0.7071,0.7071,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,02.4.2Mallat算法Mallat以多分辨分析为基础提出了著名的快速小波算法Mallat算法。小波理论获得突破性进展，使得小波分析成为近年来迅速发展起来的新兴学科并得到广泛应用。设是中的正交小波基，则对于任意的，有如下展开: (1)其中：(2)由于式1是一无穷基数，其系数d，需要按照式子2来计算，但对于来说，其一般不具有初等解析表达式。在实际的图像处理过程中，输入信号一般以数值方式给出，所以通过式2直接求取比较复杂。因此需要一种离散算法，用于找出一种新的式子来代替式2。当时，充分逼近，因此，任取，可选到充分大的，使得在上的投影： (3)记的正交投影算子为，则上式可以表示为： (4)在数学上，为了方便的进行表示，可假定，并认为，因此，关于的分解，可以近似的认为是关于的分解。因为： (5)有： (6)其中，这样的分解是唯一的。事实上，因为，所以存在着，_使得成立。其中。显然：(7)其中为向及投影的正交投影算子。且，。记：(8)则有：(9)一般地：，(10)若记为的如下算子（）：(11)把分解为和的分解过程称为有限正交小波分解，对于数字图像处理来说，这一分解形式特别有用。我们可以把定义为待分解的数字信号，则分解过程完全是离散的。 同样，也可以从和出发来重构，因而通过模拟化可得到。若是数字信号，则这一模拟过程可以省略。记的共扼算子分别为 ,即有：(12)由于：(13)所以：(14)即为由和来重构的算法，重构过程也可由式15表示： (15)Mallat二维塔式快速小波变换的分解过程如图2所示，重构过程如图3所示。表示每两行中抽取一行表示每两列中抽取一列（，分别为低通和高通滤波器）图2 二维小波分解示意图表示每两行中插入一行表示每两列中插入一列（，分别为低通和高通滤波器）图3 二维小波重构示意图Mallat 算法通过一组分解滤波器H（低通滤波器LPF）和G（高通滤波器HPF）对信号进行滤波，然后对输出结果进行下二采样（指隔一取一）来实现小波分解，分解的结果是产生长度减半的两个部分，一个是经低通滤波器产生的原始信号的平滑部分，另一个则是经高通滤波器产生的原始信号细节部分。重构时使用一组h 和g 合成滤波器对小波分解的结果滤波，再进行上二采样（相邻两点间补零）来生成重构信号。多级小波分解通过级联的方式进行，每一级的小波变换都是在前一级分解产生的低频分量上的继续，重构是分解的逆运算。低频分量上的信息比较丰富，能量集中；高频分量上的信息分量多为零，细节信息丰富，能量较少。图像经过二维小波分解后，可以得到四幅子图像 , ，和)。它们分别表示在尺度上的水平低通-垂直低通子图像，水平带通-垂直低通子图像，水平低通-垂直带通子图像，水平带通-垂直带通子图像。可以对子图像再次小波分解，得到尺度上的四幅子带图像，类似的可以对图像再次分解，以此类推，可以得到图像的多级小波分解，得到不同分辨的子带图像。图4是图像的三级小波分解表示水平低通-垂直低通子图像，表示水平低通-垂直带通子图像，表示水平低通-垂直低通子图像，水平带通-垂直带通子图像(下标表示不同的分辨率)。LL3 HL3 HL2 LH3 HH3 HL1 LH2 HH2 水平高通-垂直低通LH1 HH1水平低通-垂直高通 水平高通-垂直高通图4 图像三级小波分解示意图图像经过小波变换后，能够获得良好的空间一频率多分辨率表示，小波变换具有以下主要特征：（1）不仅保持原图像的空间特性，而且很好的提取了图像的高频信息。在低频处有很好的频率特性，在高频处有很好的空间选择性；（2）小波分量有方向选择性，分为水平、垂直、斜向，这些特性都和人类的视觉特性相吻合；（3）能量主要集中在低频子带图像；（4）低通模糊子图具有很强的相关性，水平子带图像在水平方向相关系数大，而垂直方向小;垂直子带图像在垂直方向相关系数大，而水平方向小：斜子带图像在垂直方向和水平方向相关系数都小。2.4.3阈值处理函数选取Donoho8将阈值处理函数分为软阈值和硬阈值，设是小波系数的大小，是施加阈值后的小波系数大小，是阈值。 （1）硬阈值( hard thresholding) 当小波系数的绝对值小于给定阈值时，令其为0，而大于阈值时，保持其不变，即(16) （2）软阈值(soft thresholding)当小波系数的绝对值小于给定阈值时，令其为0，大于阈值时，令其都减去阈值，即(17)（a）信号 （b）硬阈值 （c）软阈值 图5 两种阈值处理函数示意图硬阈值函数在处是不连续的，容易造成去噪后的图像在奇异点附近出现明显的Pseudo-Gibbs现象。因此，本系统选用软阈值函数作为阈值处理函数。2.4.4阈值的选取阈值的选择是离散小波去噪中最关键的一步。在去噪过程中，小波阈值起到了决定性作用：如果阈值太小，则施加阈值后小波系数将包含过多的噪声分量，达不到去噪的效果；反之，如果阈值太大，则去除了有用的成分，造成失真.所以对阈值的估计非常重要。目前，所使用的阈值包括全局阈值和局部适应阈值，各种各样的阈值公式也曾出不穷，考虑到算法实现的复杂程度以及去噪的效果，本系统采用了Donoho和Johnstone统一阈值9。其中， 为噪声标准方差, 为信号的尺寸或长度。这是在正态高斯噪声模型下，针对多维独立正态变量联合分布,在维数趋向无穷时的研究得出的结论，即大于该阈值的系数含有噪声信号的概率趋于零。这个阈值由于同信号的尺寸对数的平方根成正比，所以当 较大时,阈值趋向于将所有小波系数置零。此时小波滤噪器退化为低通滤波器。然而，在实际环境中，图像中的噪声标准方差是不能知道的，因此在选取阈值时，要对用估计方法来确定噪声标准方差。其中较常用的估算方法多采用公式18：(18)其中是小波分解尺度，是MATLAB中求中值得运算命令。这种计算方法比较复杂，去噪效果也不能令人满意。本系统采用一种简单有效的阈值估计方法，在图像进行正交小波分解的第一级，取小波系数中的HH部分，以它的标准方差的估计值，再利用Donoho和Johnstone统一阈值计算出阈值。2.4.5均值滤波模板的选取均值滤波方法，是指对待处理的当前像素，选择一个模板，该模板为其近邻的若干像素组成，用模板中像素的均值来替代原图像素值。考虑到图像数据分布的平衡性，模板选用3*3模板。3 系统实现本系统分为的实现分为读取图像，加入图像噪声，去除噪声等部分，其界面如图6所示。图6图像去噪系统界面3.1读取图像模块 图像读取模块的显示界面如图7所示。图7读取图像模块界面要利用正交小波变换对图像进行处理，首先要将BMP图像文件读取出来并存入数组中，其中最关键的是读取各点的像素值。实现这一功能的主要函数是OpenBMP()。OpenBMP()函数的功能是通过文件操作将BMP文件的文件头、信息头以及像素信息一次读入。其代码如下：BYTE * OpenBMP(CString fileName,long *width,long *height)file.Read(&bmpFileHead,sizeof(BITMAPFILEHEADER);file.Read(&bmpInfo,sizeof(BITMAPINFOHEADER);BYTE *RGBQUAD=new BYTE1024;file.Read(RGBQUAD,1024); *width = bmpInfo.bmiHeader.biWidth; *height = bmpInfo.bmiHeader.biHeight;long lFileSize = bmpFileHead.bfSize - bmpFileHead.bfOffBits;pData = (BYTE*)new charlFileSize;if(!pData)return NULL; file.ReadHuge(pData,lFileSize);file.Close();return pData;函数先打开文件，打开成功后先读入BMP文件头，再读入BMP信息头，然后读入调色板，将文件头中的BMP文件的位图宽度和高度值分别赋给变量width和height，以便后面的处理。最后根据读入的文件结构大小为pData分配内存，利用ReadHuge（）函数将像素信息读入pData，已备下一步进行正交小波变换。3.2加噪模块加噪模块的显示界面如图8所示。图8加噪模块界面由于读入内存的BMP图像为无噪图像，所以进行去噪处理必须认为加入随机噪声，此功能由RandomNoise(BYTE *Image)函数实现。其代码如下：void CLiftingschemeView:RandomNoise(BYTE *Image)for( j=0;jm_orgHeight;j+)for( i=0;im_orgWidth;i+)noisepoint=rand()/1024;ImageAddNoisej*m_orgWidth+i=(ImageAddNoisej*m_orgWidth+i*224/256+(int)noisepoint);首先将原图像的像素值复制到ImageAddNoise中，然后通过行列循环加噪，再循环内部利用rand()生成随机数，将像素值与随机数的和赋给像素值。3.3去噪模块去噪模块的显示界面如图9所示。图9加噪模块界面3.3.1正交小波变换子模块本模块采用的小波算法是Mallat算法，将复杂深奥的小波变换函数用简单的滤波器系数表示出来，大大降低了正交小波变换应用的复杂程度，使通过程序实现小波变换变得容易，此处具体采用的是Haar小波函数。本模块主要由以下两个函数实现：void CLiftingschemeView:GeneralMallat()/处理多级变换，排列子带位置double* CLiftingschemeView:RealMallat()/进行卷积运算在GeneralMallat()中，首先完成初始变量的定义，然后通过循环来实现多级小波变换和子带位置排列。当进入循环过程后，首先判断现在进行的是否为第一级小波变换，如果是第一级小波变换则初始化中间变量并在其中存储原始图像信息，如果不是第一级小波变换则不执行初始化中间变量的工作。接下来，将中间变量中存储的图像数据传入小波变换函数部分进行处理，其返回值作为低频子带信息存入中间变量中。这里需要说明的是，每进行一次小波变换，中间变量的大小都被设为此次变换前中间变量大小的1/4，因为小波变换产生的低频子带信息始终是原始数据量的1/4。该函数的流程图如图10所示。图10 正交小波变换算法流程图具体代码为：void CLiftingschemeView:GeneralMallat (double h, int h0, int ln_h, double g, int g0, int ln_g, double *Image, long ww, long hh,int level)for(loop=0;looplevel;loop+)if(loop=0)tmpImage=new doubleww*hh;for(i=0;ihh;i+)for(j=0;jww;j+)tmpImagei*ww+j=Imagei*ww+j;HHHL=new doubleww*hh/2;LLLH=new doubleww*hh/2;LLLH=RealMallat( h, h0, ln_h, g, g0, ln_g, tmpImage, ww, hh,TRUE);HHHL=RealMallat( g, g0, ln_g, h, h0, ln_h, tmpImage, ww, hh,FALSE);for(i=(hhh-hh);ihhh;i+)for(j=0;jww/2;j+)ImageTransferedi*www+j=LLLH(i-(hhh-hh)*(ww/2)+j;ImageTransferedi*www+j+ww/2=HHHL(i-(hhh-hh)*(ww/2)+j;RealMallat()中主要实现了卷积运算，在开始部分，定义了一些初始变量，这里定义的初始变量包括保存临时计算结果的中间变量、返回计算结果的返回变量，图像的高与宽和滤波器系数等。接着，先对图像数据进行行卷积计算，依次读取图像中的像素点与所选择的滤波器系数进行卷积计算。计算时，对图像左侧和右侧的边界部分数据进行边界延拓。当图像数据全部计算完毕后，对行卷积的计算结果进行抽样处理，即隔二取一，取出原始图像数据的一半进行列卷积的计算。列卷积计算过程中，同样需要对图像上部和下部的边界部分数据进行边界延拓。当列卷积计算完毕后，取出参与列卷积计算的数据中的一半作为LL 低频子带数据。此小波去噪系统利用的只是每一级小波变换后产生的LL 低频子带数据，所以对于其他三个子带的数据并不需要，它们的相关计算也不用考虑。在离散卷积的实现过程中，程序先进行的是行卷积，然后在行卷积的基础上进行的是列卷积。经过行卷积处理后，图像数据大小缩小为原始图像数据大小的1/2。经过列卷积处理后，图像数据大小在行卷积处理结果的基础上又缩小1/2，即图像数据大小缩小为原始图像数据大小的1/4。这相当于原始图像数据1/4 的数据就是低频子带LL 的数据，将它重新排列后存入数组中返回给调用函数作为下一级小波变换的处理图像并参与下一级的小波变换。该函数的流程图如图11所示。图11 正交小波变换卷积算法流程图其代码如下：double* CLiftingschemeView:RealMallat(double h, int h0, int ln_h, double g, int g0, int ln_g, double *Image, long ww, long hh, bool Low)for(i=0;ihh;i+)for(j=0;jh0-1;j+)/图像的左侧边界需要对称延拓var=0;for(k=0;kln_h;k+)/卷积运算if(j+k-(h0-1)0)var=var+hk*Imagei*ww+(abs(j+k-(h0-1)-(ln_h-1)%2);elsevar=var+hk*Imagei*ww+(j+k-(h0-1);tmpi*ww+j=var;for(j=h0-1;jww-(ln_h-h0);j+)/卷积运算for(j=ww-(ln_h-h0);jww;j+)/图像的右侧边界需要对称延拓/卷积运算 /return tmp;for(i=0;ihh;i+)for(j=0;jww/2;j+)tmpi*ww+j=tmpi*ww+2*j;/抽样/行低通计算完成/列高通，LH HL3.3.2阈值处理子模块本模块采用此模块主要包含以下两个函数：double CLiftingschemeView:GenerateThreshold()/产生阈值void CLiftingschemeView:SoftThreshold()/软阈值处理GenerateThreshold()中，先定义了一些变量，然后进行一级正交小波变换，计算小波系数中的高频子带的标准偏差作为图像噪声的估计值，最后利用Donoho提出的计算全局阈值的公式求出阈值。其代码如下：double CLiftingschemeView:GenerateThreshold(double *Image)GeneralMallat(Lo_D,3,5,Hi_D,3,3,Image,m_orgWidth,m_orgHeight,1);value=sqrt(sum_square/k);value*=sqrt(log(2*m_orgWidth*m_orgHeight);return value;SoftThreshold()利用Donoho提出的操作简单快速的软阈值函数进行对小波系数进行处理。小于阈值的像素直接赋值为0，大于或等于阈值的像素与阈值作差。其代码如下：void CLiftingschemeView:SoftThreshold(double *Image,double value)for(i=0;im_orgHeight;i+)for(j=0;jm_orgWidth;j+)if(Imagei*m_orgWidth+j0;loop-)ww=www(loop-1);hh=hhh(loop-1);tmp=new doubleww*hh;tmp=RealReconstruct(Hi_R,Hi_R0,Hi_R_l,Lo_R,Lo_R0,Lo_R_l,Image,www,hhh,loop);for(i=hhh-hh;ihhh;i+)for(j=0;jww;j+)Imagei*www+j=tmp(i+hh-hhh)*ww+j;delete tmp;RealReconstruct ()中，首先完成初始变量的定义，然后分别求出LL、LH、HL、HH四个子带的卷积结果，最后将四个结果求和存入数组SubImage中。该函数的流程图如图13所示。图13 正交小波反变换重构算法流程图其代码如下：double* CLiftingschemeView:RealReconstruct (double Hi_R, int Hi_R0, int Hi_R_l, double Lo_R, int Lo_R0, int Lo_R_l, double *Image, long www, long hhh, int loop)tmp1=SubReconstruct(Lo_R,Lo_R0,Lo_R_l,Lo_R,Lo_R0,Lo_R_l,Image,ww,hh,LL);tmp2=SubReconstruct(Lo_R,Lo_R0,Lo_R_l,Hi_R,Hi_R0,Hi_R_l,Image,ww,hh,LH);tmp3=SubReconstruct(Hi_R,Hi_R0,Hi_R_l,Lo_R,Lo_R0,Lo_R_l,Image,ww,hh,HL);tmp4=SubReconstruct(Hi_R,Hi_R0,Hi_R_l,Hi_R,Hi_R0,Hi_R_l,Image,ww,hh,HH);for(i=0;i2*hh;i+) for(j=0;j2*ww;j+)subImagei*(2*ww)+j=tmp1i*(2*ww)+j+tmp4i*(2*ww)+j+tmp2i*(2*ww)+j+tmp3i*(2*ww)+j;在SubReconstruct ()中，主要实现了卷积运算，在开始部分，定义了一些初始变量，这里定义的初始变量包括保存临时计算结果的中间变量、返回计算结果的返回变量，图像的高与宽和滤波器系数等。接着，在对图像数据进列卷积前先进行上采样操作，即间隔插零，然后对图像数据进行列卷积计算，依次读取图像中的像素点与所选择的滤波器系数进行卷积计算。计算时，对图像上侧和下侧的边界部分数据进行边界延拓。当图像数据全部计算完毕后，对列卷积的计算结果上采样，对上采样后的图像数据进行行卷积的计算。行卷积计算过程中，同样需要对图像左侧和右侧的边界部分数据进行边界延拓。当行卷积计算完毕后，将结果保存并返回。从处的卷积运算同正变换时的卷积运算算法大体相同，需要注意的是，这里先进行列卷积后进行行卷积，此外就是卷积运算前的上采样操作。该函数的流程图如图14所示。图14 正交小波反变换卷积算法流程图其代码如下：double* CLiftingschemeView:SubReconstruct (double Lo_R, int Lo_R0, int Lo_R_l, double Hi_R, int Hi_R0, int Hi_R_l, double *Image, long ww, long hh, int p)int x,y,i,j,k;double var;if(p=LL)x=1;y=0;if(p=LH)x=0;y=0;if(p=HL)x=1;y=1;if(p=HH)x=0;y=1;double *tmp1=new d