小波变换

《医学图像处理》实验报告实验十：小波变换日期： 2014年05月06日摘要本次实验的实验目的及主要内容是:一维小波变换和反变换二维小波变换和反变换二维小波细节置零、去噪一、技术讨论1.1实验原理小波变换的原理:是指一组衰减震动的波形，其振幅正负相间变化为零，是具有一定的带宽和中心频率波组。

小波变换是用伸缩和平移小波形成的小波基来分解（变换）或重构（反变换）时变信号的过程。

不同的小波具有不同带宽和中心频率，同一小波集中的带宽与中心频率的比是不变的，小波变换是一系列的带通滤波响应。

它的数学过程与傅立叶分析是相似的，只是在傅立叶分析中的基函数是单频的调和函数，而小波分析中的基函数是小波，是一可变带宽内调和函数的组合。

小波去噪的原理：利用小波变换把含噪信号分解到多尺度中，小波变换多采用二进型，然后在每一尺度下把属于噪声的小波系数去除，保留并增强属于信号的小波系数，最后重构出小波消噪后的信号。

其中关键是用什么准则来去除属于噪声的小波系数，增强属于信号的部分。

1.2实验方法1）dwt函数（实现1-D离散小波变换）[cA,cD]=dwt(X,’wname’)使用指定的小波基函数‘wname’对信号X进行分解，cA和cD分别是近似分量和细节分量；[cA,cD]=dwt(X,Lo_D,Hi_D)用指定的滤波器组Lo_D,Hi_D对信号进行分解2）idwt函数（实现1-D离散小波反变换）X=idwt(cA,cD,’wname’)X=idwt(cA,cD,Lo_R,Hi_R)X=idwt(cA,cD,’wname’,L)X=idwt(cA,cD,Lo_R,Hi_R,L)由近似分量cA和细节分量cD经过小波反变换，选择某小波函数或滤波器组，L为信号X中心附近的几个点3）dwt2函数（实现2-D离散小波变换）[cA,cH,cV,cD]=dwt2(X,’wname’)[cA,cH,cV,cD]=dwt2(X,’wname’)cA近似分量，cH水平细节分量，cV垂直细节分量，cD对角细节分量4）idwt2函数（实现2-D离散反小波变换）X=idwt2(cA,cH,cV,cD,’wname’)X=idwt2(cA,cH,cV,cD,Lo_R,Hi_R)X=idwt2(cA,cH,cV,cD,’wname’,S)X=idwt2(cA,cH,cV,cD,Lo_R,Hi_R,S)二、结果与讨论2.1实验结果一维小波变换的结果图一（a）（b）（a）此时size值取4，f为{1,4，-3,0,}（b）此时size值取6，f为{1,4，-3,0,8,3} 二维小波变换的结果:图二、正变换结果（a ） （b ）反变换结果（利用正变换结果进行反变换）（c ） （d ）置零处理结果：（三个细节分别置零处理）图三 右下角置0（a ）original pic （b ）zero（c ） wave result （d ）change图四 左下角置0（a ）original pic （b ）zero（c ）wave result （d ）change图五右上角置0（a）original pic （b）zero（c）wave result （d）change 去噪结果：图六去噪结果（a）original pic （b）zero（c ）wave result （d ）change2.2实验讨论1）一维小波变换原理是将信号分解为高频和低频两列信号，对于实验中的一维信号，实验结果图一验证了这一结果，即将信号分解为高频和低频。

2）对于二维小波变换，信号将图像信息分解为高频和低频信息。

正变换是将图像分解为高频、低频信息，图二（b ）是对（a ）进行分解的结果，（b ）的右上角、左上角以及左下角均为高频信息，左上角为低频信息。

而（d ）是根据（c ）的高频、低频信息进行重构的结果。

3）对于置零处理，分别对图像右下角、左下角以及左上角进行置零的时候，原图应该出现横向细节、垂直细节以及对角线细节被清理的后果，从图三到图五来看，可知由于图片选择问题，从当今结果不能清晰表达出来。

4）去噪是将小波变换结果的高频细节全部置零，得到的图像将比之前清晰，从图六可以简单得出此结论。

附录（实验代码）.pro程序如下：#-------------------------------------------------## Project created by QtCreator 2014-05-08T21:26:01##-------------------------------------------------QT += coreQT -= guiTARGET= shiyanshiCONFIG += consoleCONFIG -= app_bundleTEMPLATE = appSOURCES+= main.cppINCLUDEPATH+=d:\Qt\opencv2.2\include\opencv\d:\Qt\opencv2.2\include\opencv2\d:\Qt\opencv2.2\includeLIBS+=d:\Qt\opencv2.2\lib\libopencv_calib3d220.dll.a\d:\Qt\opencv2.2\lib\libopencv_contrib220.dll.a\d:\Qt\opencv2.2\lib\libopencv_core220.dll.a\d:\Qt\opencv2.2\lib\libopencv_features2d220.dll.a\d:\Qt\opencv2.2\lib\libopencv_flann220.dll.a\d:\Qt\opencv2.2\lib\libopencv_gpu220.dll.a\d:\Qt\opencv2.2\lib\libopencv_highgui220.dll.a\d:\Qt\opencv2.2\lib\libopencv_imgproc220.dll.a\d:\Qt\opencv2.2\lib\libopencv_legacy220.dll.a\d:\Qt\opencv2.2\lib\libopencv_ml220.dll.a\d:\Qt\opencv2.2\lib\libopencv_objdetect220.dll.a\d:\Qt\opencv2.2\lib\libopencv_video220.dll.aHEADERS += \I:/【快盘下载】实验十小波变换/头文件/targetver.h \I:/【快盘下载】实验十小波变换/头文件/stdafx.h \I:/【快盘下载】实验十小波变换/头文件/sdkddkver.h \.cpp程序一、一维小波变换#include "stdafx.h"#include "cv.h"#include "highgui.h"#include <stdio.h>#include <stdlib.h>#include <string.h>#include <assert.h>#include <math.h>#include <float.h>#include <limits.h>#include <time.h>#include <ctype.h>#include <math.h>#define SIZE 4#define layer 2void DWT_1_harr( double *in, double *out, int scale, int size ); void iDWT_1_harr2( double *in, double *out,int scale, int size ); int main(){double f[SIZE] = {1,4,-3,0};double q[SIZE]={0};DWT_1_harr( f, q, 2, SIZE );for (int i=0; i<SIZE; i++){printf("%lf\n",f[i]);}printf("\n");iDWT_1_harr2( f, q, 2, SIZE );for (int i=0; i<SIZE; i++){printf("%lf\n",f[i]);}char key_cmd = 1;while(1){key_cmd = cvWaitKey();if(key_cmd==27)break;}return 0;}//一维harr小波变换void DWT_1_harr( double *in, double*out,int scale, int size ){ //检查尺度是否吻合if( scale>(log(size)/log(2.0)) ){printf("Scale is too large\n");return;}int t = size;int i;while( scale>0 ){t = t/2;for (i=0; i<t; i++){out[i]=(in[2*i]+in[2*i+1])/1.4142;out[t+i]=(in[2*i]-in[2*i+1])/1.4142;}for (i=0; i<size; i++){in[i]=out[i];}scale--;}}//一维haar小波反变换（自动计算尺度）void iDWT_1_harr2( double *in, double *out, int scale, int size ){ //检查尺度是否吻合if( scale>(log(size)/log(2.0)) ){printf("Scale is too large\n");return;}int t = size/pow(2,scale);int i;while( scale>0 ){for( i=0; i<t; i++ ){out[2*i] = (in[i]+in[t+i])*0.7071;out[2*i+1] = (in[i]-in[t+i])*0.7071;}for(i=0;i<size; i++){in[i]=out[i];}t *=2;scale--;}}二、二维小波变换程序#include "stdafx.h"#include "cv.h"#include "highgui.h"#include <stdio.h>#include <stdlib.h>#include <string.h>#include <assert.h>#include <math.h>#include <float.h>#include <limits.h>#include <time.h>#include <ctype.h>#include <math.h>#define SIZE 452#define layer 1void DWT_1_harr( double *in, double *out, int scale, int size );void iDWT_1_harr( double *in, double *out, int size );void DWT_2_harr( IplImage *img_in, int scale );void iDWT_2_harr( IplImage *img_in, int scale );int main(){IplImage* img_in = cvLoadImage("d:/zero_screenshot_05.05.2014.png");CvSize size = cvSize( img_in->width, img_in->height );IplImage* img_gray = cvCreateImage(size,IPL_DEPTH_8U,1);cvCvtColor(img_in, img_gray, CV_BGR2GRAY );// IplImage* img_gray2 = cvCreateImage(size,IPL_DEPTH_64F,1);IplImage* img_gray3 = cvCreateImage(size,IPL_DEPTH_64F,1);double f[SIZE];double q[SIZE];for( int j = 0; j<size.height; j++ ) {uchar* ptr1 = (uchar*) ( img_gray->imageData + j * img_gray->widthStep );double* ptr2 = (double*) ( img_gray3->imageData + j * img_gray3->widthStep );for( int i = 0; i<size.width; i++ ) {ptr2[i] = (double)ptr1[i];}}// DWT_2_harr( img_gray3, layer );//小波变换iDWT_2_harr( img_gray3, layer );//反变换cvNormalize(img_gray3,img_gray3,1,0,CV_MINMAX);cvNamedWindow( "original pic", 1 );//创建窗口cvShowImage("originalpic", img_gray);//显示图像cvNamedWindow("wavelet pic",1);//创建窗口cvShowImage("wavelet pic", img_gray3);//显示图像char key_cmd = 1;while(1){key_cmd = cvWaitKey();if(key_cmd==27)break;}return 0;}//一维harr小波变换voidDWT_1_harr(double *in, double*out, int scale, int size){//检查尺度是否吻合if(scale>(log(size)/log(2.0))){printf("Scale is toolarge\n");return;}int t =size;int i;while(scale>0){t = t/2;for(i=0; i<t;i++){out[i]=(in[2*i]+in[2*i+1])/1.4142;out[t+i]=(in[2*i]-in[2*i+1])/1.4142;}for(i=0; i<size; i++){in[i]=out[i];}scale--;}}//一维haar小波反变换void iDWT_1_harr(double *in,double*out, int size ){//固定尺度为1int t= size/2;inti;for( i=0; i<t;i++){out[2*i]=(in[i]+in[t+i])*0.7071;out[2*i+1] = (in[i]-in[t+i])*0.7071;}for(i=0;i<size; i++){in[i]=out[i];}}//二维haar小波变换voidDWT_2_harr(IplImage*img_in, intscale){//检查尺度是否吻合CvSize size= cvSize(img_in->width,img_in->height);if(scale>(log(size.width)/log(2.0))){printf("Scaleistoolarge\n");return;}if( size.width> SIZE ){printf("Figureistoolarge,please change theSIZE parement\n");return;}//开辟变量及中间变量intstep_x =size.width;//这两个变量用于指向当前近似分量（LL）intstep_y=size.height;IplImage*img_temp =cvCreateImage(size,IPL_DEPTH_64F,1);doublef[SIZE];double q[SIZE];for( intk=1; k<=scale; k++,step_x/=2,step_y/=2){//水平分解for(intj = 0; j<step_y;j++){double* ptr1=(double*) ( img_in->imageData+j*img_in->widthStep );double*ptr2=(double*) ( img_temp->imageData+j * img_temp->widthStep );for(inti = 0; i<step_x;i++) {f[i] = ptr1[i];q[i] = 0;}DWT_1_harr( f, q, 1, step_x );for( int i = 0; i<step_x; i++ ) {ptr2[i] = f[i];}}//垂直分解cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);for( int j = 0; j<step_x; j++ ) {double* ptr1 = (double*) ( img_temp->imageData + j * img_temp->widthStep );double* ptr3 = (double*) ( img_in->imageData + j * img_in->widthStep );for( int i = 0; i<step_y; i++ ) {f[i] = ptr1[i];q[i] = 0;}DWT_1_harr( f, q, 1, step_y );for( int i = 0; i<step_y; i++ ) {ptr3[i] = f[i];}}cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);}//cvNormalize(img_temp,img_temp,1,0,CV_MINMAX);//cvNamedWindow( "c", 1 ); //创建窗口//cvShowImage( "c", img_temp ); //显示图像cvReleaseImage(&img_temp);return;}//二维haar小波反变换void iDWT_2_harr( IplImage *img_in, int scale ){//开辟变量及中间变量CvSize size = cvSize( img_in->width, img_in->height );int step_x = size.width/(pow(2,layer));int step_y = size.height/(pow(2,layer));IplImage* img_temp = cvCreateImage(size,IPL_DEPTH_64F,1);cvZero(img_temp);double f[SIZE];double q[SIZE];//cvNormalize(img_in,img_temp2,1,0,CV_MINMAX);//cvNamedWindow( "f", 1 ); //创建窗口//cvShowImage( "f", img_temp2 ); //显示图像for( int k=scale; k>=1; k--, step_x*=2, step_y*=2 ){//垂直反变换cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);for( int j = 0; j<step_x*2; j++ ) {double* ptr1 = (double*) ( img_in->imageData + j * img_in->widthStep );double* ptr2 = (double*) ( img_temp->imageData + j * img_temp->widthStep );for( int i = 0; i<step_y*2; i++ ) {f[i] = ptr1[i];q[i] = 0;}iDWT_1_harr( f, q, step_y*2 );for( int i = 0; i<step_y*2; i++ ) {ptr2[i] = f[i];}}cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);//水平反变换for( int j = 0; j<step_y*2; j++ ) {double* ptr1 = (double*) ( img_temp->imageData + j * img_temp->widthStep );double* ptr3 = (double*) ( img_in->imageData + j * img_in->widthStep );for( int i = 0; i<step_x*2; i++ ) {f[i] = ptr1[i];q[i] = 0;}iDWT_1_harr( f, q, step_x*2 );for( int i = 0; i<step_x*2; i++ ) {ptr3[i] = f[i];}}}//cvNormalize(img_temp,img_temp,1,0,CV_MINMAX);//cvNamedWindow( "d", 1 ); //创建窗口//cvShowImage( "d", img_temp ); //显示图像cvReleaseImage(&img_temp);return;}三、置零和去噪程序#include "stdafx.h"#include "cv.h"#include "highgui.h"#include <stdio.h>#include <stdlib.h>#include <string.h>#include <assert.h>#include <math.h>#include <float.h>#include <limits.h>#include <time.h>#include <ctype.h>#include <math.h>#define SIZE 300#define layer 1void DWT_1_harr( double *in, double *out, int scale, int size );void iDWT_1_harr( double *in,double*out, int size );void DWT_2_harr( IplImage *img_in, int scale );void iDWT_2_harr( IplImage *img_in, int scale );int main(){IplImage* img_in = cvLoadImage("d:/lena.tif");CvSize size = cvSize( img_in->width, img_in->height );IplImage* img_gray = cvCreateImage(size,IPL_DEPTH_8U,1);cvCvtColor(img_in, img_gray, CV_BGR2GRAY );IplImage* img_gray2 = cvCreateImage(size,IPL_DEPTH_64F,1);IplImage* img_gray3 = cvCreateImage(size,IPL_DEPTH_64F,1);double f[SIZE];double q[SIZE];for( int j = 0; j<size.height; j++ ) {uchar* ptr1 = (uchar*) ( img_gray->imageData + j * img_gray->widthStep );double* ptr2 = (double*) ( img_gray3->imageData + j * img_gray3->widthStep );for( int i = 0; i<size.width; i++ ) {ptr2[i] = (double)ptr1[i];}}//小波变换DWT_2_harr( img_gray3, layer );cvNormalize(img_gray3,img_gray2,1,0,CV_MINMAX);cvNamedWindow( "wavelet_result", 1 ); //创建窗口cvShowImage( "wavelet_result", img_gray2 ); //显示图像//置0for( int j = 128; j<size.width; j++ ) {//128<j<256 行数double* ptr2 = (double*) ( img_gray3->imageData + j * img_gray3->widthStep );for( int i = 128; i<size.width; i++ ) {//128<=i<256 列数ptr2[i] = 0;//把右下角细节置0}}for( int j = 128; j<size.width; j++ ) {//128<j<256double* ptr2 = (double*) ( img_gray3->imageData + j * img_gray3->widthStep );for( int i = 0; i<128; i++ ) {//0<=i<128ptr2[i] = 0;//把左下角细节置0}}for( int j = 0; j<128; j++ ) {//0<j<128double* ptr2 = (double*) ( img_gray3->imageData + j * img_gray3->widthStep );for( int i = 128; i<size.width; i++ ) {//128<=i<256ptr2[i] = 0;//把右上角细节置0}}cvNormalize(img_gray3,img_gray2,1,0,CV_MINMAX);cvNamedWindow( "change_coefficients", 1 ); //创建窗口cvShowImage( "change_coefficients", img_gray2 ); //显示图像iDWT_2_harr( img_gray3, layer );//反变换cvNormalize(img_gray3,img_gray3,1,0,CV_MINMAX);cvNamedWindow( "original pic", 1 ); //创建窗口cvShowImage( "original pic", img_gray ); //显示图像cvNamedWindow( "zero", 1 ); //创建窗口cvShowImage( "zero", img_gray3 ); //显示图像char key_cmd = 1;while(1){key_cmd = cvWaitKey();if(key_cmd==27)break;}return 0;}//一维harr小波变换void DWT_1_harr( double *in, double *out, int scale, int size ){ //检查尺度是否吻合if( scale>(log(size)/log(2.0)) ){printf("Scale is too large\n");return;}int t =size;int i;while( scale>0 ){t = t/2;for (i=0; i<t; i++){out[i]=(in[2*i]+in[2*i+1])/1.4142;out[t+i]=(in[2*i]-in[2*i+1])/1.4142;}for (i=0; i<size; i++){in[i]=out[i];}scale--;}}//一维haar小波反变换voidiDWT_1_harr( double *in, double*out,intsize){//固定尺度为1intt=size/2;int i;for(i=0;i<t; i++ ){out[2*i] =(in[i]+in[t+i])*0.7071;out[2*i+1] =(in[i]-in[t+i])*0.7071;}for(i=0; i<size; i++){in[i]=out[i];}}//二维haar小波变换voidDWT_2_harr(IplImage *img_in, int scale ){//检查尺度是否吻合CvSize size= cvSize( img_in->width,img_in->height );if( scale>(log(size.width)/log(2.0)) ){printf("Scaleis too large\n");return;}if( size.width > SIZE ){printf("Figureis too large, pleasechangetheSIZEparement\n");return;}//开辟变量及中间变量int step_x=size.width;//这两个变量用于指向当前近似分量（LL）int step_y= size.height;IplImage*img_temp=cvCreateImage(size,IPL_DEPTH_64F,1);doublef[SIZE];doubleq[SIZE];for(int k=1;k<=scale; k++, step_x/=2,step_y/=2){//水平分解for( int j =0;j<step_y;j++ ){double*ptr1 = (double*)( img_in->imageData+j*img_in->widthStep);double*ptr2= (double*) (img_temp->imageData+ j* img_temp->widthStep);for(inti= 0;i<step_x; i++){f[i]=ptr1[i];q[i]=0;}DWT_1_harr( f,q, 1,step_x );for(inti=0; i<step_x;i++){ptr2[i]=f[i];}}//垂直分解cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);for(intj=0; j<step_x;j++ ) {double*ptr1=(double*)( img_temp->imageData+j*img_temp->widthStep);double*ptr3=(double*)(img_in->imageData+j*img_in->widthStep );for(int i=0;i<step_y;i++){f[i]=ptr1[i];q[i] =0;}DWT_1_harr(f, q, 1,step_y );for( int i=0; i<step_y;i++){ptr3[i] = f[i];}}cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);}//cvNormalize(img_temp,img_temp,1,0,CV_MINMAX);//cvNamedWindow( "c", 1 ); //创建窗口//cvShowImage( "c", img_temp ); //显示图像cvReleaseImage(&img_temp);return;}//二维haar小波反变换void iDWT_2_harr( IplImage *img_in, int scale ){//开辟变量及中间变量CvSize size = cvSize( img_in->width, img_in->height );int step_x = size.width/(pow(2,layer));int step_y = size.height/(pow(2,layer));IplImage* img_temp = cvCreateImage(size,IPL_DEPTH_64F,1);cvZero(img_temp);double f[SIZE];double q[SIZE];//cvNormalize(img_in,img_temp2,1,0,CV_MINMAX);//cvNamedWindow( "f", 1 ); //创建窗口//cvShowImage( "f", img_temp2 ); //显示图像for( int k=scale; k>=1; k--, step_x*=2, step_y*=2 ){//垂直反变换cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);for( int j = 0; j<step_x*2; j++ ) {double* ptr1 = (double*) ( img_in->imageData + j * img_in->widthStep );double* ptr2 = (double*) ( img_temp->imageData + j * img_temp->widthStep );for( int i = 0; i<step_y*2; i++ ) {f[i] = ptr1[i];q[i] = 0;}iDWT_1_harr( f, q, step_y*2 );for( int i = 0; i<step_y*2; i++ ) {ptr2[i] = f[i];}}cvTranspose(img_temp,img_temp);cvTranspose(img_in,img_in);//水平反变换for( int j = 0; j<step_y*2; j++ ) {double* ptr1 = (double*) ( img_temp->imageData + j * img_temp->widthStep );double* ptr3 = (double*) ( img_in->imageData + j * img_in->widthStep );for( int i = 0; i<step_x*2; i++ ) {f[i] = ptr1[i];q[i] = 0;}iDWT_1_harr( f, q, step_x*2 );for( int i = 0; i<step_x*2; i++ ) {ptr3[i] = f[i];}}}//cvNormalize(img_temp,img_temp,1,0,CV_MINMAX);//cvNamedWindow( "d", 1 ); //创建窗口//cvShowImage( "d", img_temp ); //显示图像cvReleaseImage(&img_temp);return;}。