参考论文中的文字：图像配准是图像处理的基本任务之一,用于将不同时间、不同传感器、不同视角及不同拍摄条件下获取的关于同一目标或场景的两幅或多幅图像进行主要是几何意义上的匹配套和的过程。在对图像配准的研究过程中,大量技术被应用于针对不同数据和问题的图像配准工作,产生了多种不同形式的图像配准技术。 图像配准的基本问题是找出一种图像转换方法,用以纠正图像的形变。造成图像形变的原因多种多样,例如对于遥感图像而言,传感器噪声、由传感器视点变化或平台不稳定造成的透视变化、被拍摄物体的移动、变形或生长等变化、闪电和大气变化,以及阴影和云层遮盖都使图像产生不同形式的形变。正是图像形变原因和形式的不同决定了多种多样的图像配准技术。 迄今已报道了多种图像配准方法,主要有互相关法、傅立叶变换法、点映射法口脚外和弹性模型法。其中傅立叶变换法基于傅立叶变换的相位匹配是利用傅立叶变换的性质而出现的一种图像配准方法。图像经过傅立叶变换,由空域变换到频率缘则两组数据在空何上的相关运算可以变为频谱的复数乘法运算,同时图像在变换域中还能获得在空域中很难获得的特征。 # 一，基于相位相关的图像配准方法 在时域中信号的平移运动可以通过在频域中相位的变化表现出来（这是傅里叶变换的特性）。同理,图像的旋转、平移和比例变化也能在傅里叶变换的频域中反映出来。而且使用频域方法的好处是计算简单,同时傅立叶变换可以采用方法提高执行的速度。因此,傅氏变换是图像配准中常用的方法之一。下面我们就具体分析当图像发生平移、旋转和缩放时,图像信号在频域中的表现。 ###1，原理阐述  通过求取互功率谱的傅立叶反变换,得到一个狄拉克函数（脉冲函数）,再寻找函数峰值点对应的坐标,即可得到我们所要求得的配准点。实际上,在计算机处理中,连续域要用离散域代替,这使得狄拉克函数转化为离散时间单位冲击函数序列的形式。在实际运算中,两幅图像互功率谱相位的反变换,总是含有一个相关峰值代表两幅图像的配准点,和一些非相关峰值,相关峰值直接反映两幅图像间的一致程度。更精确的讲,相关峰的能量对应重叠区域的所占百分比,非相关峰对应非重叠区域所占百分比。由此我们可以看出,当两幅图像重叠区域较小时,采用本方法就不能检测出两幅图像的平移量。 ###2，配准情况 当图像间仅存在平移时，正确的配准图像如图a所示（中心平移化了），最大峰的位置就是两图像的相对平移量，反之若不存在单纯的平移，则会出现如b所示的情况（多脉冲林立）  ### 3，算法流程  # 二，程序实现 ### 1，核心函数的实现： ~~~ void PhaseCorrelation2D(const BYTE *signal,//原信号 const BYTE *pattern,//带配准信号 int &height_offset,//高的偏移量 int &width_offset)//宽的偏移量 { fftw_complex *signal_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol); fftw_complex *pattern_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol); for (int i = 0; i < nRow*nCol; i++) { signal_img[i][0] = signal[i]; signal_img[i][1] = 0; } for (int j = 0; j < nRow*nCol; j++) { pattern_img[j][0] = pattern[j]; pattern_img[j][1] = 0; } // 对两幅图像傅里叶变换 fftw_plan signal_forward_plan = fftw_plan_dft_2d(nRow, nCol, signal_img, signal_img, FFTW_FORWARD, FFTW_ESTIMATE); fftw_plan pattern_forward_plan = fftw_plan_dft_2d(nRow, nCol, pattern_img, pattern_img, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(signal_forward_plan); fftw_execute(pattern_forward_plan); // 求互功率谱 fftw_complex *cross_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol); double temp; for (int i = 0; i < nRow*nCol; i++) { cross_img[i][0] = (signal_img[i][0] * pattern_img[i][0]) - (signal_img[i][1] * (-1.0*pattern_img[i][1])); cross_img[i][1] = (signal_img[i][0] * (-1.0*pattern_img[i][1])) + (signal_img[i][1] * pattern_img[i][0]); temp = sqrt(cross_img[i][0] * cross_img[i][0] + cross_img[i][1] * cross_img[i][1]) + 0.001; cross_img[i][0] /= temp; cross_img[i][1] /= temp; } //对互功率谱求反变换 fftw_plan cross_backward_plan = fftw_plan_dft_2d(nRow, nCol, cross_img, cross_img, FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute(cross_backward_plan); // 释放内存 fftw_destroy_plan(signal_forward_plan); fftw_destroy_plan(pattern_forward_plan); fftw_destroy_plan(cross_backward_plan); fftw_free(signal_img); fftw_free(pattern_img); double *cross_real=new double[nRow*nCol]; for (int i = 0; i < nRow*nCol; i++) cross_real[i] = cross_img[i][0]; int max_loc = 0;//准备存放最大值的位置坐标（注意，只有一个值） double max_vlaue = 0.0; for (int i = 0; i < nRow*nCol; i++) { if (max_vlaue 0.5*nRow) height_offset = height_offset - nRow; if (width_offset > 0.5*nCol) width_offset = width_offset - nCol; delete[] cross_real; cross_real = NULL; } ~~~ ### 2，配准模拟结果：  ### 3，实际的图像配准效果 注：下面这张图是单帧图像效果，有大量的椒盐噪声和散斑噪声  注：下面这张图像是八张单帧图像（图像内容可能一样，但是可能存在较小的位移）的叠加效果，散斑噪声和椒盐噪声得到明显去除  注：下面这张图像是八张图像先配准再叠加平均的效果，背景噪声及其散斑噪声得到消除，并且图像对比度得到一定程度提高  注：本博文为[EbowTang](http://my.csdn.net/EbowTang)原创，后续可能继续更新本文。如果转载，请务必复制本条信息！ 原文地址：http://blog.csdn.net/ebowtang/article/details/47816245 原作者博客：http://blog.csdn.net/ebowtang # 参考资源： 【1】西北工业大学，吕海霞，硕士论文，《自动图像配准技术研究》，2007.3 【2】MIT,FFTW库 【3】Joseph L. Horner and Peter D. Gianino. Phase-only matched ltering. Applied Optics, 23(6):812-816, 1984. 【4】复旦大学，宋智礼，博士论文，《图像配准技术及其应用的研究》，2010,4 【5】中国科学技术大学信号统计处理研究院，郑志彬，叶中付，《基于相位相关的图像配准算法》，2006,12 【6】(美)冈萨雷斯著; 阮秋琦译. 数字图像处理(MATLAB 版) [M]. 北京：电子工业出版社, 2005. 【7】(美)冈萨雷斯著; 阮秋琦等译. 数字图像处理(第二版) [M]. 北京：电子工业出版社, 2003.

# 一，经典滤波算法的基本原理 ###1，中值滤波和均值滤波的基本原理 参考以前转载的博客：http://blog.csdn.net/ebowtang/article/details/38960271 ###2，高斯平滑滤波基本原理 参考以前转载的博客：http://blog.csdn.net/ebowtang/article/details/38389747 # 二，噪声测试效果 ### 1，不同噪声效果 三幅图各噪声浓度分别是0.01 0.03,0.05（比如第一副图均是加入0.01的噪声浓度）    ### 2，实验代码 ~~~ %读入原始图像并显示 image_original=imread('dog.bmp'); figure(1) subplot(2,4,1); imshow(image_original); title('原输入图像'); axis square; %生成含高斯噪声图像并显示 pp=0.05; image_gaosi_noise=imnoise(image_original,'gaussian',0,pp); subplot(2,4,2); imshow(image_gaosi_noise); title('添加高斯噪声后图像'); axis square; %生成含椒盐噪声图像并显示 d=0.05; image_saltpepper_noise=imnoise(image_original,'salt & pepper',d); subplot(2,4,3); imshow(image_saltpepper_noise); title('添加椒盐噪声后图像'); axis square; %生成含乘性噪声图像并显示 var=0.05; image_speckle_noise=imnoise(image_original,'speckle',var); subplot(2,4,4); imshow(image_speckle_noise); title('添加乘性噪声后图像'); axis square; %原图像直方图 r=0:255; bb=image_original(:); pg=hist(bb,r); pgr1=pg/length(bb); subplot(245);bar(pgr1);title('源输入图像的直方图'); r=0:255; bl=image_gaosi_noise(:); pg=hist(bl,r); pgr2=pg/length(bl); subplot(246);bar(pgr2);title('高斯噪声污染后的直方图'); r=0:255; bh=image_saltpepper_noise(:); pu=hist(bh,r); pgr3=pu/length(bh); subplot(247);bar(pgr3);title('椒盐噪声污染后的直方图'); r=0:255; ba=image_speckle_noise(:); pa=hist(ba,r); pgr4=pa/length(ba); subplot(248);bar(pgr4);title('乘性噪声污染后直方图'); ~~~ # 三，椒盐噪声去除能力对比 ### 1，三大去噪效果 三幅图椒盐噪声浓度分别是0.01 0.03,0.05（比如第一副图均是加入0.01的椒盐噪声去噪对比）    ### 2，实现代码 ~~~

%读入原始图像并显示
image_original=imread('dog.bmp');
figure(1)
subplot(2,4,1);
imshow(image_original);
title('原输入图像');
axis square;

%生成含高斯噪声图像并显示
%pp=0.05;
%image_gaosi_noise=imnoise(image_original,'gaussian',0,pp);

%生成含椒盐噪声图像并显示
dd=0.05;
image_saltpepper_noise=imnoise(image_original,'salt & pepper',dd);

%生成含乘性噪声图像并显示
%var=0.05;
%image_speckle_noise=imnoise(image_original,'speckle',var);

image_saltpepper_noise_after1=medfilt2(image_saltpepper_noise,[3,3]);
subplot(2,4,2);
imshow(image_saltpepper_noise_after1);title('中值滤波去椒盐噪声效果图');
axis square;

h_gaosi1=fspecial('gaussian',3,1);
image_saltpepper_noise_after2=imfilter(image_saltpepper_noise,h_gaosi1);
subplot(2,4,3);
imshow(image_saltpepper_noise_after2);title('高斯平滑去椒盐噪声效果');
axis square;

image_saltpepper_noise_after3=wiener2(image_saltpepper_noise,[5 5]);
subplot(2,4,4);
imshow(image_saltpepper_noise_after3);title('维纳滤波去椒盐噪声效果');
axis square;

%原图像直方图
r=0:255;  
bb=image_original(:); 
pg=hist(bb,r);  
pgr1=pg/length(bb);  
subplot(245);bar(pgr1);title('源输入图像的直方图');

r=0:255;  
bl=image_saltpepper_noise_after1(:); 
pg=hist(bl,r);  
pgr2=pg/length(bl);  
subplot(246);bar(pgr2);title('中值滤波去椒盐噪声后的直方图');

r=0:255;  
bh=image_saltpepper_noise_after2(:); 
pu=hist(bh,r);  
pgr3=pu/length(bh);  
subplot(247);bar(pgr3);title('高斯平滑去椒盐噪声后的直方图');

r=0:255;  
ba=image_saltpepper_noise_after3(:); 
pa=hist(ba,r);  
pgr4=pa/length(ba);  
subplot(248);bar(pgr4);title('维纳滤波去除椒盐噪声后的直方图');
~~~

# 四，高斯噪声去除能力对比

### 1，去噪效果对比

![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526aa186a.jpg)
![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526ac1bee.jpg)
![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526ae8bee.jpg)

### 2，实验代码

~~~
%读入原始图像并显示
image_original=imread('dog.bmp');
figure(1)
subplot(2,4,1);
imshow(image_original);
title('原输入图像');
axis square;

%生成含高斯噪声图像并显示
pp=0.05;
image_gaosi_noise=imnoise(image_original,'gaussian',0,pp);

%生成含椒盐噪声图像并显示
%dd=0.01;
%image_saltpepper_noise=imnoise(image_original,'salt & pepper',dd);

%生成含乘性噪声图像并显示
%var=0.05;
%image_speckle_noise=imnoise(image_original,'speckle',var);

image_gaosi_noise_after1=medfilt2(image_gaosi_noise,[3,3]);
subplot(2,4,2);
imshow(image_gaosi_noise_after1);title('中值滤波去高斯噪声效果图');
axis square;

h_gaosi1=fspecial('gaussian',3,1);
image_gaosi_noise_after2=imfilter(image_gaosi_noise,h_gaosi1);
subplot(2,4,3);
imshow(image_gaosi_noise_after2);title('高斯平滑去高斯噪声效果');
axis square;

image_gaosi_noise_after3=wiener2(image_gaosi_noise,[5 5]);
subplot(2,4,4);
imshow(image_gaosi_noise_after3);title('维纳滤波去高斯噪声效果');
axis square;

%原图像直方图
r=0:255;  
bb=image_original(:); 
pg=hist(bb,r);  
pgr1=pg/length(bb);  
subplot(245);bar(pgr1);title('源输入图像的直方图');

r=0:255;  
bl=image_gaosi_noise_after1(:); 
pg=hist(bl,r);  
pgr2=pg/length(bl);  
subplot(246);bar(pgr2);title('中值滤波去高斯噪声后的直方图');

r=0:255;  
bh=image_gaosi_noise_after2(:); 
pu=hist(bh,r);  
pgr3=pu/length(bh);  
subplot(247);bar(pgr3);title('高斯平滑去高斯噪声后的直方图');

r=0:255;  
ba=image_gaosi_noise_after3(:); 
pa=hist(ba,r);  
pgr4=pa/length(ba);  
subplot(248);bar(pgr4);title('维纳滤波去除高斯噪声后的直方图');
~~~

# 五，乘性噪声去除能力对比

### 1，去噪效果对比

![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526b0fcaf.jpg)
![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526b2a4a6.jpg)
![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526b3f248.jpg)

### 2，实验代码

~~~
%读入原始图像并显示
image_original=imread('dog.bmp');
figure(1)
subplot(2,4,1);
imshow(image_original);
title('原输入图像');
axis square;

%生成含高斯噪声图像并显示
%pp=0.01;
%image_gaosi_noise=imnoise(image_original,'gaussian',0,pp);

%生成含椒盐噪声图像并显示
%dd=0.01;
%image_saltpepper_noise=imnoise(image_original,'salt & pepper',dd);

%生成含乘性噪声图像并显示
var=0.01;
image_speckle_noise=imnoise(image_original,'speckle',var);

image_speckle_noise_after1=medfilt2(image_speckle_noise,[3,3]);
subplot(2,4,2);
imshow(image_speckle_noise_after1);title('中值滤波去乘性噪声效果图');
axis square;

h_gaosi1=fspecial('gaussian',3,1);
image_speckle_noise_after2=imfilter(image_speckle_noise,h_gaosi1);
subplot(2,4,3);
imshow(image_speckle_noise_after2);title('高斯平滑去乘性噪声效果');
axis square;

image_speckle_noise_after3=wiener2(image_speckle_noise,[5 5]);
subplot(2,4,4);
imshow(image_speckle_noise_after3);title('维纳滤波去乘性噪声效果');
axis square;

%原图像直方图
r=0:255;  
bb=image_original(:); 
pg=hist(bb,r);  
pgr1=pg/length(bb);  
subplot(245);bar(pgr1);title('源输入图像的直方图');

r=0:255;  
bl=image_speckle_noise_after1(:); 
pg=hist(bl,r);  
pgr2=pg/length(bl);  
subplot(246);bar(pgr2);title('中值滤波去乘性噪声后的直方图');

r=0:255;  
bh=image_speckle_noise_after2(:); 
pu=hist(bh,r);  
pgr3=pu/length(bh);  
subplot(247);bar(pgr3);title('高斯平滑去乘性噪声后的直方图');

r=0:255;  
ba=image_speckle_noise_after3(:); 
pa=hist(ba,r);  
pgr4=pa/length(ba);  
subplot(248);bar(pgr4);title('维纳滤波去除乘性噪声后的直方图');
~~~
  

# 六，PNSR客观对比

（PNSR客观对比越高越好）

本对比也囊括了其他常见去噪方式的对比

![](https://docs.gechiui.com/gc-content/uploads/sites/kancloud/2016-03-02_56d6526b5f9bd.jpg)
  

  

  

  

#参考资源

【1】《百度百科》

【2】《维基百科》

【3】冈萨雷斯《数字图像处理》

【4】http://blog.csdn.net/ebowtang/article/details/38960271
                    

本文代码的实现严重依赖前面的两篇文章： [一维信号的小波阈值去噪](http://blog.csdn.net/ebowtang/article/details/40481393#t14) [小波变换一维Mallat算法的C++实现](http://blog.csdn.net/ebowtang/article/details/40433861#t0) 图像在获取或传输过程中会因各种噪声的干扰使质量下降,这将对后续图像的处理产生不利影响.所以必须对图像进行去噪处理,而去噪所要达到的目的就是在较好去除噪声的基础上,良好的保持图像的边缘等重要细节.在图像去噪领域得到广泛的应用.本博文根据小波的分解与重构原理，实现了基于硬阈值和软阈值函数的小波阈值去噪的C++版本，最终结果与matlab库函数运算结果完全一致。 注本文的大部分文字提取于参考论文 # 一，小波阈值去噪基本理论 ###小波阈值处理 小波阈值收缩法是Donoho和Johnstone提出的,其主要理论依据是,小波变换具有很强的去数据相关性,它能够使信号的能量在小波域集中在一些大的小波系数中;而噪声的能量却分布于整个小波域内.因此,经小波分解后,信号的小波系数幅值要大于噪声的系数幅值.可以认为,幅值比较大的小波系数一般以信号为主，而幅值比较小的系数在很大程度上是噪声.于是,采用阈值的办法可以把信号系数保留,而使大部分噪声系数减小至零.小波阈值收缩法去噪的具体处理过程为:将含噪信号在各尺度上进行小波分解,设定一个阈值，幅值低于该阈值的小波系数置为0,高于该阈值的小波系数或者完全保留,或者做相应的“收缩(shrinkage)”处理.最后将处理后获得的小波系数用逆小波变换进行重构,得到去噪后的图像. ###阈值函数的选取   阈值去噪中,阈值函数体现了对超过和低于阈值的小波系数不同处理策略，是阈值去噪中关键的一步。设w表示小波系数,T为给定阈值,sign(*)为符号函数,常见的阈值函数有： 硬阈值函数：     （小波系数的绝对值低于阈值的置零，高于的保留不变）       软阈值函数：           值得注意的是： 1) 硬阈值函数在阈值点是不连续的，在下图中已经用黑线标出。不连续会带来振铃，伪吉布斯效应等。 2) 软阈值函数，原系数和分解得到的小波系数总存在着恒定的偏差,这将影响重构的精度使得重构图像的边缘模糊等现象. 同时这两种函数不能表达出分解后系数的能量分布。见下图：                                                                                 图1 于是不少文章出现了折衷方案，一种新的阈值函数（非本文重点）：  参考代码： ~~~ %Generate signal and set threshold. y = linspace(-1,1,100); subplot(311); plot(y);title('原始线'); thr = 0.5; % Perform hard thresholding. ythard = wthresh(y,'h',thr); subplot(312); plot(ythard);title('硬阈值线'); % Perform soft thresholding. ytsoft = wthresh(y,'s',thr); subplot(313); plot(ytsoft);title('软阈值线'); ~~~ ### 阈值的确定 选取的阈值最好刚好大于噪声的最大水平，可以证明的是噪声的最大限度以非常高的概率低于 （此阈值是Donoho提出的），其中根号右边的这个参数（叫做sigma）就是估计出来的噪声标准偏差（第一级分解出的小波细节系数，即整个HH系数的绝对值的中值），本文将用此阈值去处理各尺度上的细节系数。 ### 阈值策略 曾经做的ppt挪用过来的。  ###二维信号分解与重构 小波的分解  小波的重构  # 二，Matlab库函数实现 ### 1，核心库函数说明 1）wavedec2 图像的多级小波分解，将返回分解出来的小波系数以及小波系数的各级长度 2）waverec2 多级小波系数的重构，重构出原信号 3）wthresh函数 对系数进行指定类型（全局阈值或者分层阈值）的阈值去噪处理 更具体的函数说明可以在，matlab里键入“doc 函数名”将得到很详细的说明，当然也可以百度 ### 2，软，硬阈值处理效果：     局部放大图像： 四幅图象均放大两倍，便于查看区别     ### 3，完整的代码实现 说明：代码实现对图像一层分解后的细节系数进行全局阈值处理（即LHL,LH,HH均用同一阈值处理），并且阈值是自定义的。 ~~~ clc; clear all; close all; %%%%%%%%%%%%%%%%%%%%%%%%通过matlab的函数来实现阈值去噪%%%%%%%%%%%%%%%%%%%%%%%%%% % X=imread('sample1.bmp');%有噪声的图像 X=double(X); figure(1); imshow(X,[]);title('原始图像'); %注意，原始图像带有噪声 % 获取输入参数 w = 'db3';%小波类型 n = 1;%分解层数 thr = 23.91;%自定义阈值 sorh1 = 'h';%硬阈值 sorh2 = 's';%软阈值 % 对图像进行小波分解 [c,l] = wavedec2(X,n,w); % 对小波系数全局阈值处理 cxch = c;% 保留近似系数 cxcs = c;% 保留近似系数 justdet = prod(l(1,:))+1:length(c);%截取细节系数（不处理近似系数） % 阈值处理细节系数 cxch(justdet) = wthresh(cxch(justdet),sorh1,thr);%硬阈值去噪 cxcs(justdet) = wthresh(cxcs(justdet),sorh2,thr);%软阈值去噪 %小波重建 xch = waverec2(cxch,l,w); xcs = waverec2(cxcs,l,w); figure(2); imshow(uint8(xch));title('硬阈值去噪图像'); figure(3); imshow(uint8(xcs));title('软阈值去噪图像'); ~~~ # 三，C加加实现 说明：如同一维的阈值去噪一样，在执行自己编写的wavedec2函数时必须先初始化，初始化的目的是为了获取信号的长度，选择的是什么小波，以及分解的等级等信息，然后计算出未来的各种信息，比如每个等级的系数的size，其中共有变量m_msgCL2D记录了这些信息。二维小波分解的初始化函数如下： ~~~ //初始化二维图像的分解信息，保存未来需要的信息 bool CWavelet::InitDecInfo2D( const int height,//预分解的图像的高度 const int width,//预分解的图像的宽度 const int Scale,//分解尺度 const int dbn//db滤波器编号，有默认值 ) { if (dbn != 3) SetFilter(dbn); if (height < m_dbFilter.filterLen - 1 || width < m_dbFilter.filterLen - 1) { cerr << "错误信息：滤波器长度大于信号的高度或者宽度!" << endl; return false; } int srcHeight = height; int srcWidth = width; m_msgCL2D.dbn = dbn; m_msgCL2D.Scale = Scale; m_msgCL2D.msgHeight.resize(Scale + 2); m_msgCL2D.msgWidth.resize(Scale + 2); //源图像的尺寸 m_msgCL2D.msgHeight[0] = height; m_msgCL2D.msgWidth[0] = width; //每一尺度上的尺寸 for (int i = 1; i <= Scale; i++)//注意：每个尺度的四个分量的长宽是一样的 { int exHeight = (srcHeight + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度的一半 srcHeight = exHeight; m_msgCL2D.msgHeight[i] = srcHeight; int exWidth = (srcWidth + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度一半 srcWidth = exWidth; m_msgCL2D.msgWidth[i] = srcWidth; } m_msgCL2D.msgHeight[Scale + 1] = srcHeight; m_msgCL2D.msgWidth[Scale + 1] = srcWidth; //计算总的数据个数 int tmpAllSize = 0; int curPartSize = 0; int prePartSize = 0; for (int i = 1; i <= Scale; i++) { curPartSize = m_msgCL2D.msgHeight[i] * m_msgCL2D.msgWidth[i]; tmpAllSize += curPartSize * 4 - prePartSize; prePartSize = curPartSize; } m_msgCL2D.allSize = tmpAllSize; m_bInitFlag2D = true; return true; } ~~~ ### 2，核心函数的实现 ### 1）二维信号的单次分解 说明：本函数建立在一维的小波分解函数基础上（DWT） ~~~ // 二维数据的小波分解 void CWavelet::DWT2( double *pSrcImage,//源图像数据(存储成一维数据，行优先存储) int height,//图像的高 int width,//图像的宽 double *pDstCeof//分解出来的图像系数 ) { if (!m_bInitFlag2D) { cerr << "错误信息：未初始化，无法对信号进行分解!" << endl; return; } if (pSrcImage == NULL || pDstCeof == NULL) { cerr << "错误信息：dwt2数据无内存" << endl; Sleep(3000); exit(1); } int exwidth = (width + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的宽度 int exheight = (height + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的高度 double *tempImage = new double[exwidth*height]; // 对每一行进行行变换 double *tempAhang = new double[width]; double *tempExhang = new double[exwidth]; // 临时存放每一行的处理数据 for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) tempAhang[j] = pSrcImage[i*width + j];//提取每一行的数据 DWT(tempAhang, width, tempExhang); for (int j = 0; j < exwidth; j++) tempImage[i*exwidth + j] = tempExhang[j]; } // 对每一列进行列变换 double *tempAlie = new double[height]; // 临时存放每一列的转置数据 double *tempexlie = new double[exheight]; // 临时存放每一列的处理数据 for (int i = 0; i < exwidth; i++) { // 列转置 for (int j = 0; j < height; j++) tempAlie[j] = tempImage[j*exwidth + i];//提取每一列数据 //执行变换 DWT(tempAlie, height, tempexlie); // 反转置 for (int j = 0; j < exheight; j++) pDstCeof[j*exwidth + i] = tempexlie[j]; } AdjustData(pDstCeof, exheight, exwidth);//调整数据 delete[] tempAlie; tempAlie = NULL; delete[] tempexlie; tempexlie = NULL; delete[] tempAhang; tempAhang = NULL; delete[] tempExhang; tempExhang = NULL; delete[] tempImage; tempImage = NULL; } ~~~ ### 2）二维信号的单次重构 说明： ~~~ //二维小波反变换 void CWavelet::IDWT2( double *pSrcCeof, //二维源图像系数数据 int dstHeight,//重构出来后数据的高度 int dstWidth,//重构出来后数据的宽度 double *pDstImage//重构出来的图像 ) { int srcHeight = (dstHeight + m_dbFilter.filterLen - 1) / 2 * 2; int srcWidth = (dstWidth + m_dbFilter.filterLen - 1) / 2 * 2;//pSrcCeof的高度 IAdjustData(pSrcCeof, srcHeight, srcWidth);//调整成LL,HL,LH,HH double *tempAline = new double[srcHeight]; // 临时存放每一列的数据 double *tempdstline = new double[dstHeight]; // 临时存放每一列的重构结果 double *pTmpImage = new double[srcWidth*dstHeight]; // 列重构 for (int i = 0; i < srcWidth; i++)//每一列 { // 列转置 for (int j = 0; j= 1; i--) { int nextheight = m_msgCL2D.msgHeight[i - 1];//重构出来的高度 int nextwidth = m_msgCL2D.msgWidth[i - 1];//重构出来的宽度 IDWT2(pTmpImage, nextheight, nextwidth, pDstData); if (i > 1)//i==1已经重构出来了，不再需要提取系数 { for (int j = 0; j < nextheight*nextwidth; j++) pTmpImage[j] = pDstData[j]; for (int j = 0; j < 3 * nextheight*nextwidth; j++) pTmpImage[nextheight*nextwidth + j] = pSrcCoef[gap + j]; gap += 3 * nextheight*nextwidth; } } delete[] pTmpImage; pTmpImage = NULL; return true; } ~~~ ### 3，函数正确性验证 ### 1）二维单次分解与重构测试  附带本次测试代码： ~~~ //小波函数的二维分解与重构 int main() { system("color 0A"); CWavelet cw; double s[36] = { 1, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27, 3, 4, 15, 6, 72, 8, 41, 5, 6, 7, 8, 9, 5, 64, 7, 8, 9, 14, 6, 27, 8, 9, 40, 31 }; int height = 6; int width = 6; for (int j = 0; j < 36; j++) { cout << s[j] << " "; if ((j + 1) % 6 == 0) cout << endl; } cout << endl; cout << endl; cw.InitDecInfo2D(height, width, 1, 3); double *dst = new double[cw.m_msgCL2D.allSize]; cw.DWT2(s, height, width, dst); for (int j = 0; j < cw.m_msgCL2D.allSize; j++) { cout << dst[j] << " "; if ((j + 1) % 10 == 0) cout << endl; } cout << endl; double *dsts = new double[36]; cw.IDWT2(dst, height, width, dsts); for (int j = 0; j < 36; j++) { cout << dsts[j] << " "; if ((j + 1) % 6 == 0) cout << endl; } delete[] dst; dst = NULL; delete[] dsts; dsts = NULL; system("pasue"); return 0; } ~~~ ### 2）二维多级分解与重构测试 说明：对二维数据进行了5层分解，选取的是小波族db3  附带本次测试代码： ~~~ //小波函数二维多级分解与重构测试 int main() { system("color 0A"); CWavelet cw; double s[48] = { 1, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27, 3, 4, 15, 6, 72, 8, 41, 5, 6, 7, 8, 9, 5, 64, 7, 8, 9, 14, 6, 27, 8, 9, 40, 31 , 1, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27 }; int height = 6; int width = 8; for (int j = 0; j < 48; j++) { cout << s[j] << " "; if ((j + 1) % 8 == 0) cout << endl; } cout << endl; int Scale = 5; int dbn = 2; cw.InitDecInfo2D(height, width, Scale, dbn); double *dstcoef = new double[cw.m_msgCL2D.allSize]; cw.WaveDec2(s,dstcoef); for (int i = 0; i < cw.m_msgCL2D.allSize; i++) { cout << dstcoef[i] << " "; if ((i + 1) % 10 == 0) cout << endl; } double *dst = new double[48]; for (int i = 0; i < 48; i++) dst[i] = 0.0; cw.WaveRec2(dstcoef, dst); cout << endl; cout << endl; for (int i = 0; i < 48; i++) { cout << dst[i] << " "; if ((i + 1) % 8 == 0) cout << endl; } cout << endl; delete[] dst; dst = NULL; delete[] dstcoef; dstcoef = NULL; system("pause"); return 0; } ~~~ ### 3，阈值去噪结果： 说明：本测试只是模拟测试，对图像的处理也是一样的（完全一致） ### 硬阈值去噪结果  ### 软阈值去噪结果  实际的图像处理结果为： 源噪声图像为：  注意以下是采用：db6，3层分解，软阈值去噪，阈值是在前文提及的阈值基础上缩小2.5倍得到的效果：  注意以下是采用：db6，3层分解，硬阈值去噪，阈值是在前文提及的阈值基础上缩小2.5倍得到的效果：  附带上述C++测试验证主函数 ~~~ //二维阈值去噪测试 int main() { system("color 0A"); CWavelet cw; double s[48] = { 10, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27, 3, 4, 15, 6, 72, 8, 41, 5, 6, 7, 8, 9, 5, 64, 7, 8, 9, 14, 6, 27, 8, 9, 40, 31, 10, 12, 30, 4, 50, 61, 2, 3, 41, 5, 6, 27}; int height = 6; int width = 8; for (int j = 0; j < 48; j++) { cout << s[j] << " "; if ((j + 1) % 8 == 0) cout << endl; } cout << endl; int Scale = 3; int dbn = 3; cw.InitDecInfo2D(height, width, Scale, dbn); double *dstcoef = new double[48]; if (!cw.thrDenoise2D(s, dstcoef)) cerr << "Error" << endl; for (int j = 0; j < 48; j++) { cout << dstcoef[j] << " "; if ((j + 1) % 8 == 0) cout << endl; } delete[] dstcoef; dstcoef = NULL; system("pause"); return 0; } ~~~ 附带上述matlab验证程序 ~~~ clc; clear all; close all; % %%%%%%%%%%%%%%%%%%%%%%%%通过matlab的函数来实现阈值去噪%%%%%%%%%%%%%%%%%%%%%%%%%% % X=[ 10, 12, 30, 4, 5, 61, 2, 3; 41, 5, 6, 27, 3, 4, 15, 6; 72, 8, 41, 5, 6, 7, 8, 9; 5, 64, 7, 8, 9, 14, 6, 27; 8, 9, 40, 31,10, 12, 30, 4; 50, 61, 2, 3, 41, 5, 6, 27]; X=double(X); % 获取输入参数 wname = 'db3';%小波类型 n = 3;%分解层数 sorh1 = 'h';%硬阈值 sorh2 = 's';%软阈值 % 对图像进行小波分解 [c,l] = wavedec2(X,n,wname); %求取阈值 N = numel(X); [chd1,cvd1,cdd1] = detcoef2('all',c,l,1); cvd1=cvd1(:)'; sigma = median(abs(cvd1))/0.6745;%提取细节系数求中值并除以0.6745 thr = sigma*sqrt(2*log(N)); % 对小波系数全局阈值处理 cxch = c;% 保留近似系数 cxcs = c;% 保留近似系数 justdet = prod(l(1,:))+1:length(c);%截取细节系数（不处理近似系数） % 阈值处理细节系数 cxch(justdet) = wthresh(cxch(justdet),sorh1,thr);%硬阈值去噪 cxcs(justdet) = wthresh(cxcs(justdet),sorh2,thr);%软阈值去噪 %小波重建 xch = waverec2(cxch,l,wname); xcs = waverec2(cxcs,l,wname); ~~~ 注：本博文为[EbowTang](http://my.csdn.net/EbowTang)原创，后续可能继续更新本文。如果转载，请务必复制本条信息！ 原文地址：http://blog.csdn.net/ebowtang/article/details/40481539 原作者博客：http://blog.csdn.net/ebowtang # 参考资源： 【1】《数字图像处理》（冈萨雷斯matlab第二版） 【2】http://ivm.sjtu.edu.cn/files/wavelet/第3章wavelet_original.pdf 【3】http://media.cs.tsinghua.edu.cn/~ahz/digitalimageprocess/chapter12/chapt12_ahz.htm#a1 【4】http://wenku.baidu.com/link?url=OYRL2n-cYkZ2J10zaMscZQ-lhR05kysQ_CaB1YM1e_aqr3DakexZRm8rtBYOHlDmxC0cNAtiCopjyog_yOIH1zliUmyz2fKfOzFTAQ1wWj3 【5】《维基百科》 【6】来自中国知网若干论文 【7】小波分析及其应用__孙延奎 【8】杨建国.小波分析及其工程应用[M].北京：机械工业出版社.2005 【9】毛艳辉.小波去噪在语音识别预处理中的应用.上海交通大学硕士学位论文.2010 【10】matlab各种函数说明，及其内部函数实现

本文代码的实现严重依赖前面的一篇文章： [小波变换Mallat算法的C++实现](http://blog.csdn.net/ebowtang/article/details/40433861) # 基本理论 小波阈值收缩法是Donoho和Johnstone提出的,其主要理论依据是,小波变换具有很强的去数据相关性,它能够使信号的能量在小波域集中在一些大的小波系数中;而噪声的能量却分布于整个小波域内.因此,经小波分解后,信号的小波系数幅值要大于噪声的系数幅值.可以认为,幅值比较大的小波系数一般以信号为主，而幅值比较小的系数在很大程度上是噪声.于是,采用阈值的办法可以把信号系数保留,而使大部分噪声系数减小至零.小波阈值收缩法去噪的具体处理过程为:将含噪信号在各尺度上进行小波分解,设定一个阈值，幅值低于该阈值的小波系数置为0,高于该阈值的小波系数或者完全保留,或者做相应的“收缩(shrinkage)”处理.最后将处理后获得的小波系数用逆小波变换进行重构,得到去噪后的信号. ### 1，阈值函数的选取   阈值去噪中,阈值函数体现了对超过和低于阈值的小波系数不同处理策略，是阈值去噪中关键的一步。设w表示小波系数,T为给定阈值,sign(*)为符号函数,常见的阈值函数有： 硬阈值函数：     （小波系数的绝对值低于阈值的置零，高于的保留不变）       软阈值函数：           值得注意的是： 1) 硬阈值函数在阈值点是不连续的，在下图中已经用黑线标出。不连续会带来振铃，伪吉布斯效应等。 2) 软阈值函数，原系数和分解得到的小波系数总存在着恒定的偏差,这将影响重构的精度使得重构图像的边缘模糊等现象. 同时这两种函数不能表达出分解后系数的能量分布。见下图：                                                                                 图1 于是不少文章出现了折衷方案，一种新的阈值函数（非本文重点）：  阈值线参考代码： ~~~ %Generate signal and set threshold.    y = linspace(-1,1,100);    subplot(311);   plot(y);title('原始线');   thr = 0.5;   % Perform hard thresholding.    ythard = wthresh(y,'h',thr);   subplot(312);   plot(ythard);title('硬阈值线');   % Perform soft thresholding.    ytsoft = wthresh(y,'s',thr);   subplot(313);   plot(ytsoft);title('软阈值线');   ~~~ ###2，阈值的确定 选取的阈值最好刚好大于噪声的最大水平，可以证明的是噪声的最大限度以非常高的概率低于 （此阈值是Donoho提出的），其中根号右边的这个参数（叫做sigma）就是估计出来的噪声标准偏差（根据第一级分解出的小波细节系数，即整个det1绝对值系数中间位置的值），本文将用此阈值去处理各尺度上的细节系数，注意所谓全局阈值就是近似系数不做任何阈值处理外，其他均阈值处理。 最后吐槽一下这个“绝对值系数中间位置的值” 1）如果det1的长度为偶数那么，这个“中值”便是中间位置的两个数之和的平均值，比如【2,2,3,5】，中值即是2.5而不是3 2）如果det1的长度为奇数那么，这个中值就是中间位置的那个数，比如【2,3,5】，中值即3 ###3，阈值策略 以前写的ppt挪用过来：  ### 4，一维信号的分解与重构 以下算法如果用简单的文字描述，可是：先将信号对称拓延（matlab的默认方式），然后再分别与低通分解滤波器和高通分解滤波器卷积，最后下采样，最后可以看出最终卷积采样的长度为floor(n-1)/2+n,如果想继续分解下去则继续对低频系数CA采取同样的方式进行分解。  # 一，matlab库函数实现 ### 1，核心库函数说明 1）wnoisest函数 作用：估计一维小波高频系数中的噪声偏差 这个估计值使用的是绝对值中间位置的值（估计的噪声偏差值）除以0.6745（Median Absolute Deviation / 0.6745），适合0均值的高斯白噪声 2）wavedec函数 一维信号的多尺度分解，将返回诸多细节系数和每个系数的长度，在matlab中键入“doc wavedec”具体功能一目了然 3）waverec函数 一维信号小波分解系数的重构，将返回重构后的信号在matlab中键入“doc waverec”具体功能一目了然，也可以键入“open waverec”查看matlab具体是怎么做的。 4）wdencmp函数 这个函数用于对一维或二维信号的压缩或者去噪，使用方法： 1 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('gbl',X,'wname',N,THR,SORH,KEEPAPP) 2 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',X,'wname',N,THR,SORH) 3 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',C,L,'wname',N,THR,SORH) wname是所用的小波函数， gbl(global的缩写)表示每层都采用同一个阈值进行处理， lvd表示每层用不同的阈值进行处理， N表示小波分解的层数， THR为阈值向量， 对于格式(2)(3)每层都要求有一个阈值，因此阈值向量THR的长度为N， SORH表示选择软阈值还是硬阈值(分别取为’s’和’h’)， 参数KEEPAPP取值为1时，则低频系数不进行阈值量化处理，反之，则低频系数进行阈值量化。 XC是消噪或压缩后的信号，[CXC，LXC]是XC的小波分解结构， PERF0和PERFL2是恢复和压缩L^2的范数百分比, 是用百分制表明降噪或压缩所保留的能量成分。 ###2，阈值去噪效果 从图中可以查出，除燥效果还是比较理想的，对于噪声比较重的地方软阈值去噪能力更加明显（因为没有无噪的信号参考，这并不能代表他比硬阈值更优秀）。  放大其中的细节部分，便于查看细节  ### 3，完整的matlab代码 ~~~ clc; clear; % 获取噪声信号 load leleccum; indx = 1:3450; noisez = leleccum(indx); %信号的分解 wname = 'db3'; lev = 3; [c,l] = wavedec(noisez,lev,wname); %求取阈值 sigma = wnoisest(c,l,1);%使用库函数wnoisest提取第一层的细节系数来估算噪声的标准偏差 N = numel(noisez);%整个信号的长度 thr = sigma*sqrt(2*log(N));%最终阈值 %全局阈值处理 keepapp = 1;%近似系数不作处理 denoisexs = wdencmp('gbl',c,l,wname,lev,thr,'s',keepapp); denoisexh = wdencmp('gbl',c,l,wname,lev,thr,'h',keepapp); % 作图 subplot(311), plot(noisez), title('原始噪声信号'); subplot(312), plot(denoisexs), title('matlab软阈值去噪信号') ; subplot(313), plot(denoisexh), title('matlab硬阈值去噪信号') ; ~~~ # 二，C加加实现 说明，一维信号的单尺度分解在前一篇文章中已经提及，这里不再累述，这里主要再次基础上的多尺度分解与重构 并且在执行自己编写的wavedec函数时必须先初始化，初始化的目的是为了获取信号的长度，选择的是什么小波，以及分解的等级等信息，然后计算出未来的各种信息，比如每个等级的系数的size，为了进行一维小波分解的初始化函数如下： ~~~ bool CWavelet::InitDecInfo( const int signalLen,//源信号长度 const int decScale,//分解尺度 const int decdbn//db滤波器的编号 ) { if (decdbn != 3) SetFilter(decdbn); if (signalLen < m_dbFilter.filterLen - 1) { cerr << "错误信息：滤波器长度大于信号!" << endl; return false; } int srcLen = signalLen; m_msgCL1D.dbn = decdbn; m_msgCL1D.Scale = decScale; m_msgCL1D.msgLen.resize(decScale + 2); m_msgCL1D.msgLen[0] = srcLen; for (int i = 1; i <= decScale; i++) { int exLen = (srcLen + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度 srcLen = exLen; m_msgCL1D.msgLen[i] = srcLen; } m_msgCL1D.msgLen[decScale + 1] = srcLen; for (int i = 1; i < decScale + 2; i++) m_msgCL1D.allSize += m_msgCL1D.msgLen[i]; m_bInitFlag1D = true;//设置为已经初始化 return true; } ~~~ ### 1，核心函数的实现 ### 1），信号的多级分解 注：本函数实现了对信号的任意级数分解，分解的全部系数与matlab的结果完全一致 ~~~ // 一维多尺度小波分解,必须先初始化 //分解的尺度等信息已经在初始化函数获取 bool CWavelet::WaveDec( double *pSrcData,//要分解的信号 double *pDstCeof//分解出来的系数 ) { if (pSrcData == NULL || pDstCeof == NULL) return false; if (!m_bInitFlag1D) { cerr << "错误信息：未初始化，无法对信号进行分解!" << endl; return false; } int signalLen = m_msgCL1D.msgLen[0]; int decLevel = m_msgCL1D.Scale; double *pTmpSrc = new double[signalLen]; double *pTmpDst = new double[m_msgCL1D.msgLen[1] * 2]; for (int i = 0; i < signalLen; i++) pTmpSrc[i] = pSrcData[i]; int gap = m_msgCL1D.msgLen[1] * 2; for (int i = 1; i <= decLevel; i++) { int curSignalLen = m_msgCL1D.msgLen[i - 1]; DWT(pTmpSrc, curSignalLen, pTmpDst); for (int j = 0; j < m_msgCL1D.msgLen[i] * 2; j++) pDstCeof[m_msgCL1D.allSize - gap + j] = pTmpDst[j]; for (int k = 0; k < m_msgCL1D.msgLen[i]; k++) pTmpSrc[k] = pTmpDst[k]; gap -= m_msgCL1D.msgLen[i]; gap += m_msgCL1D.msgLen[i + 1] * 2; } delete[] pTmpDst; pTmpDst = NULL; delete[] pTmpSrc; pTmpSrc = NULL; return true; } ~~~ ### 2），多级分解系数的重构 注：本函数只能还原成原始信号，不能还原到分解级数的中间某个位置 ~~~ // 重构出源信号 bool CWavelet::WaveRec( double *pSrcCoef,//源被分解系数 double *pDstData//重构出来的信号，两者的长度是一样的 ) { if (pSrcCoef == NULL || pDstData == NULL)//错误：无内存 return false; //从m_msgCL1D中获取分解信息 int signalLen = m_msgCL1D.msgLen[0];//信号长度 int decLevel = m_msgCL1D.Scale;//分解级数 int det1Len = m_msgCL1D.msgLen[1]; double *pTmpSrcCoef = new double[det1Len * 2]; for (int i = 0; i < m_msgCL1D.msgLen[decLevel] * 2; i++) pTmpSrcCoef[i] = pSrcCoef[i]; int gap = m_msgCL1D.msgLen[decLevel] * 2; for (int i = decLevel; i >= 1; i--) { int curDstLen = m_msgCL1D.msgLen[i - 1]; IDWT(pTmpSrcCoef, curDstLen, pDstData); if (i != 1) { for (int j = 0; j < curDstLen; j++) pTmpSrcCoef[j] = pDstData[j]; for (int k = 0; k < curDstLen; k++) pTmpSrcCoef[k + curDstLen] = pSrcCoef[k + gap]; gap += m_msgCL1D.msgLen[i - 1]; } } delete[] pTmpSrcCoef; pTmpSrcCoef = NULL; return true; } ~~~ ###3），阈值的获取 注：严格依照Donoho的阈值写的代码 ~~~ // 根据细节系数，以及信号长度计算阈值 double CWavelet::getThr( double *pDetCoef,//细节系数（应该是第一级的细节系数） int detLen,//此段细节系数的长度 bool is2D//当前细节系数是否来自是二维图像信号的 ) { double thr = 0.0; double sigma = 0.0; for (int i = 0; i < detLen; i++) pDetCoef[i] = abs(pDetCoef[i]); std::sort(pDetCoef, pDetCoef + detLen); if (detLen % 2 == 0 && detLen >= 2) sigma = (pDetCoef[detLen / 2-1] + pDetCoef[detLen / 2]) / 2 / 0.6745; else sigma = pDetCoef[detLen / 2] / 0.6745; if (!is2D)//一维信号 { double N = m_msgCL1D.msgLen[0]; thr = sigma *sqrt(2.0*log(N)); } else{//二维信号 double size = m_msgCL2D.msgHeight[0]*m_msgCL2D.msgWidth[0]; thr = sigma *sqrt(2.0*log(size)); } return thr; } ~~~ ### 4），高频系数阈值处理 注：本函数只对高频系数做处理，不对近似系数处理 ~~~ // 将系数阈值处理，一维二维均适用 void CWavelet::Wthresh( double *pDstCoef,//细节系数（应该是除近似系数外的所有的细节系数） double thr,//阈值 const int allsize,//分解出来的系数的总长度（非） const int gap,//跳过最后一层的近似系数 SORH ish//阈值函数的选取 ) { // if (ish)//硬阈值 { for (int i = gap; i < allsize; i++) { if (abs(pDstCoef[i]) < thr)//小于阈值的置零，大于的不变 pDstCoef[i] = 0.0; } } else//软阈值 { for (int i = gap; i < allsize; i++) { if (abs(pDstCoef[i]) < thr)//小于阈值的置零，大于的收缩 { pDstCoef[i] = 0.0; } else { if (pDstCoef[i] < 0.0) pDstCoef[i] = thr - abs(pDstCoef[i]); else pDstCoef[i] = abs(pDstCoef[i]) - thr; } } } } ~~~ ### 5），阈值去噪函数 注：本函数涉及到上面提及的多个函数，此函数是核心的对外接口 ~~~ bool CWavelet::thrDenoise( double *pSrcNoise,//源一维噪声信号 double *pDstData,//去噪后的信号 bool isHard//阈值函数的选取，有默认值 ) { if (pSrcNoise == NULL || pDstData == NULL) exit(1); if (!m_bInitFlag1D)//错误：未初始化 return false; double *pDstCoef = new double[m_msgCL1D.allSize]; WaveDec(pSrcNoise, pDstCoef);//分解出系数 int Det1Len = m_msgCL1D.msgLen[1]; int gapDet = m_msgCL1D.allSize - Det1Len; double *pDet1 = new double[Det1Len]; for (int i = gapDet, j = 0; i < m_msgCL1D.allSize; i++, j++) pDet1[j] = pDstCoef[i]; int gapApp = m_msgCL1D.msgLen[m_msgCL1D.Scale];//跳过最后一层的近似系数 double thr = getThr(pDet1, Det1Len);//获取阈值 Wthresh(pDstCoef, thr, m_msgCL1D.allSize, gapApp, isHard);//将细节系数阈值 WaveRec(pDstCoef, pDstData);//重构信号 delete[] pDstCoef; pDstCoef = NULL; delete[] pDet1; pDet1 = NULL; return true; } ~~~ ### 2，函数正确性验证 注：本次测试实现了10级分解，并与matlab进行了对比，结果完全一致（对比未完全展示）  以下数据就是上述中的第二行  附函数测试主函数： ~~~ //测试一维数据的任意级数的分解与还原（重构，只能重构为源信号） int _tmain(int argc, _TCHAR* argv[]) { system("color 0A"); double signal[23] = { 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92, 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92, 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79 }; int signalLen = sizeof(signal) / sizeof(double); cout << "原始信号：" << endl; for (int i = 0; i < signalLen; i++) cout << signal[i] << " "; CWavelet cw; int decdbn = 3; int declevel = 10; if (!cw.InitDecInfo(signalLen, declevel, decdbn)) { Sleep(5000); exit(1); } cout<> signal[i]; double *pDen = new double[signal_size]; CWavelet cw; int scale = 3; int dbn = 3; cw.InitDecInfo(signal_size,scale,dbn); cw.thrDenoise(signal, pDen, true); for (int i = 0; i < signal_size; i++) waveOut << pDen[i] << endl; delete[] pDen; pDen = NULL; delete[] signal; signal = NULL; system("pause"); return 0; } ~~~ 注：本博文为[EbowTang](http://my.csdn.net/EbowTang)原创，后续可能继续更新本文。本着共享、开源精神可以转载，但请务必复制本条信息！ 原文地址：http://blog.csdn.net/ebowtang/article/details/40481393 原作者博客：http://blog.csdn.net/ebowtang # 参考资源： 【1】网友，邹宇华，博客地址，http://blog.csdn.net/chenyusiyuan/article/details/2862768 【2】《维基百科----小波变换》 【3】乔世杰.小波图像编码中的对称边界延拓法[ J].中国图像图形学报,2000,5(2):725-729. 【4】MALLAT S.A theory formulti-resolution signal decompo-sition: the wavelet representation[ J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 1989, 11(4):674-693. 【5】《小波十讲》 【6】《小波与傅里叶分析基础》 【7】冈萨雷斯《数字图像处理》 【8】matlab小波算法说明文档 【9】阈值去噪鼻祖论文，Donoho, D.L. (1995), "De-noising by soft-thresholding," IEEE Trans. on Inf. Theory, 41, 3, pp. 613–627. **声明：** **本文部分文字学习并整理自网络，部分代码参考于网络资源.** **如果侵犯了您的版权，请联系本人tangyibiao520@163.com，本人将及时编辑掉！** [](http://info.flagcounter.com/4ycf) 

# 一，一维信号的拓延 在Mallat算法的推导中,假定输入序列是无限长的,而实际应用中常常是分时采样,即输入序列为有限长.此时,滤波器系数与输入序列卷积时就会出现轮空的现象.因此有必要对原始信号进行边界延拓,减小边界误差.解决的方法通常有补零法和周期延拓法. 1)补零法是在输入序列的末尾补零.补零法的缺点是可能人为造成输入序列边界的不连续,从而使得较高频率的小波系数很大. 2)周期延拓法是将原来有限长的输入序列拓展成周期的序列.周期延拓可适用于任何小波变换,但可能导致输入序列边缘的不连续,使得高频系数较大.这种方式的拓延卷积后与源信号的长度一致。 3)对称延拓（matlab默认采取这种方式）可避免输入序列边缘的不连续,但只适用于对称小波变换.本文根据Mallat算法公式,编写了对称延拓方式的小波变换的一般实现方法. 注：笔者采用的编译器为VS2013，当前系统为win10。 ### 1，常见的3种边缘拓延方法 设输入信号为f(n),长度为srcLen,滤波器长度为filterLen.下面给出信号边界处理几种方法的具体表达式如下： 1）周期拓延：             f(n+ srcLen)   ,  -(filterLen -1)≤ n< 0   f′(n)=  f(n)                 ,  0≤ n≤ srcLen -1             f(n -srcLen)    ,  srcLen -1< n≤srcLen+filterLen -2 举例说明：以“1 2 3 4 5 6 7 8”这个长度为8的信号为例，当滤波器的长度为4时，其具体的拓延长度为6（单边为3）： 6 7 8 （1 2 3 4 5 6 7 8）1 2 3 2）对称延拓（本文重点）             f(-n -1)                     ,   -(filterLen-1)≤ n< 0 f′(n)=    f(n)                          ,   0≤ n≤srcLen-1             f(2srcLen -n -1)       ,   srcLen-1< n≤ srcLen+ filterLen-2 举例说明：以“1 2 3 4 5 6 7 8”这个长度为8的信号为例，当滤波器的长度为4时，其具体的拓延长度为6（单边为3）： 3 2 1 （1 2 3 4 5 6 7 8）8 7 6 3）零值填补              0                ,   -(filterLen -1)≤ n< 0 f′(n)=     f(n)            ,  0≤ n≤srcLen -1              0                ,  srcLen -1< n≤srcLen+filterLen -2 举例说明：以“1 2 3 4 5 6 7 8”这个长度为8的信号为例，当滤波器的长度为4时，其具体的拓延长度为6（单边为3）： 0 0 0 （1 2 3 4 5 6 7 8）0 0 0  #二，Mallat算法分解过程 在Mallat算法中,信号分解过程按照系数形式,  其中:hd(n)和gd(n)分别为分解的低通和高通滤波器系数,长度为filterLen . 吐槽：信号滤波的过程其实就是滤波器与信号卷积的过程，也就是滤波器的频谱和信号的频谱相乘的过程。 那么低通滤波器频谱是高频低，低频高，与信号相乘时高频的会被过滤掉；高通滤波器同理！ ### 1，详细的分解过程 以db2小波为例，通过对称拓延后的详细计算过程如下：  从上图可知，小波的Mallat算法分解后的序列长度由原序列长srcLen和滤波器长filterLen决定。从Mallat算法（采用对称拓延）的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进行隔点抽取而来。即分解抽取的结果长度为(srcLen+filterLen-1)/2。 上述方法也是matlab默认采取的方法，在MATLAB中输入代码，查看结果便知： ~~~ [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db2');%获取小波系数 xx=[420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92]; wname='db2'; level=1; [c,l]=wavedec(xx,level,wname);%进行一层分解 ~~~ 结果为（被存储在系数矩阵c中，l记录的是每段系数的长度）： 595.429871699852  597.374655846484598.449909371632 604.108998389031607.245626013478 -2.03920021086643   0.509501871423169-1.289053089583672.33989974581161  -1.14513645475084 与上图一致 ### 2，C++实现分解算法 获取db2小波系数： ~~~ //系数精度很高，采自matlab extern double db2_Lo_D[4] = { -0.129409522550921, 0.224143868041857, 0.836516303737469, 0.482962913144690 }; extern double db2_Hi_D[4] = { -0.482962913144690, 0.836516303737469, -0.224143868041857, -0.129409522550921 }; extern double db2_Lo_R[4] = { 0.482962913144690, 0.836516303737469, 0.224143868041857, -0.129409522550921 }; extern double db2_Hi_R[4] = { -0.129409522550921, -0.224143868041857, 0.836516303737469, -0.482962913144690 }; ~~~ C++实现如下： ~~~ // 一维信号的小波分解 int CWavelet::DWT( double *pSrcData,//分解的源信号 int srcLen,//源信号的长度 double *pDstCeof//分解出来的,本函数将返回此长度 ) { //本程序禁止出现这种情况，否则数据出错（对称拓延长度为filterLen-1，如果大于了signalLen将越界） if (srcLen < m_dbFilter.filterLen - 1) { //实际上信号的长度可以是任意的（matlab顺序：信号拓延-》卷积-》下采样）， //但是本程序为了算法速度，写法上不允许 cerr << "错误信息：滤波器长度大于信号!" << endl; Sleep(1000); exit(1); } int exLen = (srcLen + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度 int k = 0; double tmp = 0.0; for (int i = 0; i < exLen; i++) { pDstCeof[i] = 0.0; pDstCeof[i + exLen] = 0.0; for (int j = 0; j < m_dbFilter.filterLen; j++) { k = 2 * i - j + 1; //信号边沿对称延拓 if ((k<0) && (k >= -m_dbFilter.filterLen + 1))//左边沿拓延 tmp = pSrcData[-k - 1]; else if ((k >= 0) && (k <= srcLen - 1))//保持不变 tmp = pSrcData[k]; else if ((k>srcLen - 1) && (k <= (srcLen + m_dbFilter.filterLen - 2)))//右边沿拓延 tmp = pSrcData[2 * srcLen - k - 1]; else tmp = 0.0; pDstCeof[i] += m_dbFilter.Lo_D[j] * tmp; pDstCeof[i + exLen] += m_dbFilter.Hi_D[j] * tmp; } } return 2 * exLen; } ~~~ 处理结果为：（与matlab一致）  附带测试主函数： ~~~ int main() { system("color 0A"); double s[8] = { 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92}; int signalLen = sizeof(s) / sizeof(double); cout << "原始信号：" << endl; for (int i = 0; i < signalLen; i++) cout << s[i] << " "; cout << endl; CWavelet cw(2); double *dst = new double[10]; cw.DWT(s, signalLen,dst); cout << "分解后的系数：" << endl; for (int i = 0; i < 10; i++) cout << dst[i] << " "; cout << endl; delete[] dst; dst = NULL; system("pause"); return 0; } ~~~ # 三，Mallat算法重构过程  ### 1，详细的重构过程 以db2小波为例，通过对称拓延后的详细重构过程如下：  在matlab中我们可以通过以下代码实现重构（完全实现原信号） ~~~ [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db2');%获取小波系数 xx=[420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92]; wname='db2'; level=1; [c,l]=wavedec(xx,level,wname);%进行一层分解 cc=waverec(c,l,wname); ~~~ ###2，C++实现重构算法 ~~~ // 一维小波反变换，重构出源信号 void CWavelet::IDWT( double *pSrcCoef,//源分解系数 int dstLen,//重构出来的系数的长度 double *pDstData//重构出来的系数 ) { int p = 0; int caLen = (dstLen + m_dbFilter.filterLen - 1) / 2; for (int i = 0; i < dstLen; i++) { pDstData[i] = 0.0; for (int j = 0; j < caLen; j++) { p = i - 2 * j + m_dbFilter.filterLen - 2; //信号重构 if ((p >= 0) && (p ';

之前一段我们提到的算法都是和平滑有关， 经过平滑算法之后， 图像锐度降低， 降低到一定程度， 就变成了模糊。 今天我们反其道行之， 我们看看锐化是怎么做的。 这里的锐化， 还是的从平滑谈开去。我们先来观察原来的图像和平滑图像的区别： 原图 raw：   模糊图 blur：   _________________________________________________________ 源图像减去模糊图像等于 锐化后的边缘mask效果：            这个时候， 我们发现，减法做完的这个图赫然勾勒出了原图的边缘！！ 这样给我们一个启示就是， 如果我们把这个mask加到原图上那岂不就是锐化了？ （不明白？ 锐化的意思就是边缘的色差比较大， 产生的图片貌似清晰的效果） 说干就干， 马上我们来做个新的算式: 老图 raw:    加上mask    _______________________________________________________ 等于锐化图  sharpen    怎么样， 是不是有了锐化的效果了？？所以我们实际上的锐化效果就是从这么简单的想法衍生出来的。 所以锐化的公式可以简单的表述为 sharp = raw + ( raw-blur ); 再来看看我们原来的高斯模版的话就是这样：  这样的话， 我们的锐化算法，也变得和之前的高斯平滑差不多了， 就是像素的加权平均值的计算就可以得到了。可以想见的事情是代码肯定也会出奇的一致！ 这是那个template改掉了： ~~~ void sharpenImage (unsigned char* gray, unsigned char* smooth, int width, int height) { int templates[25] = { -1, -4, -7, -4, -1, -4, -16, -26, -16, -4, -7, -26, 505, -26, -7, -4, -16, -26, -16, -4, -1, -4, -7, -4, -1 }; memcpy ( smooth, gray, width*height*sizeof(unsigned char) ); for (int j=2;j 255) sum = 255; if (sum <0) sum = 0; smooth [ j*width+i ] = sum; } } } ~~~ 当然， 这个锐化算法或者说锐化的模板只是我根据前面的算式自己计算的来的，其实还是有非常主流的锐化模版可以供使用的， 比如说著名的拉普拉斯算子。 参考资源： 【1】http://blog.csdn.net/hhygcy/article/details/4330939

#中值滤波 中值滤波法是一种非线性平滑技术，它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值.中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。方法是用某种结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升（或下降）的为二维数据序列。二维中值滤波输出为g（x,y）=med{f(x-k,y-l),(k,l∈W)} ，其中，f(x,y)，g(x,y)分别为原始图像和处理后图像。W为二维模板，通常为3*3，5*5区域，也可以是不同的的形状，如线状，圆形，十字形，圆环形等。 **以3*3的滤波窗口为例：** 中值滤波就是在3*3中的像素中寻找中值。 来看这样一个描述图（不然无图无真相）  这把可以清晰地看到， 这里有6，2，0，3，97，4，19，3，10这些像素， 然后中间的这些像素值就被这些像素的中位数也就是中值取代了。为了满足和前面一篇文章的格式相对应， 我们马上进入下一个单元， 来看看在平滑和降噪方面的功效！ 原图1                                                                           中值滤波之后    噪声图(5%)                                                                中值滤波后：    很明显中值滤波不仅使得图像变得平滑（也可以理解为图像变得模糊了），但是其对去除了椒盐噪声特别有效，因为中值滤波就是专门针对去除孤立点噪声孕育为生的。从这里中值的逻辑来看， 我们做中值操作的时候， 那么白色(接近255的像素点)和黑色(接近0的像素)因为在窗口排序时是最大最小值， 除非周围的颜色都是黑色或者白色，不然一般都会被剔除掉，所以，一般来说这个中值滤波是去除椒盐噪声的非常理想的选择。 最后还是贴一段我运行的代码： ~~~ void medianFilter(BYTE* pcorrPic, BYTE* psmoothPic, int width, int height) { memcpy(psmoothPic, pcorrPic, width*height*sizeof(BYTE)); for (int j = 1; j < height - 1; j++)//排除邻边 { for (int i = 1; i < width - 1; i++) { int k = 0; BYTE window[9]; for (int jj = j - 1; jj < j + 2; ++jj) for (int ii = i - 1; ii < i + 2; ++ii) window[k++] = pcorrPic[jj * width + ii]; // 元素排序 (一半即可) for (int m = 0; m < 5; ++m) { int min = m; for (int n = m + 1; n < 9; ++n) if (window[n] < window[min]) min = n; // 找到最小值应该放置的位置 BYTE temp = window[m]; window[m] = window[min]; window[min] = temp; } psmoothPic[j*width + i] = window[4]; } } } ~~~ #均值滤波 “**把每个像素都用周围的9个像素来做均值操作** ”， 比如说这里有一个例子：    非常明显的， 这个3*3区域像素的颜色值分别是5,3,6,2,1,9,8,4,7那么中间的1这个像素的过滤后的值就是这些值的平均值， 也就是前面的计算方法： (5+3+6+2+1+9+8+4+7)/9=5 一目了然。那么这个均值滤波器有什么用处呢？ 主要还是平滑图像的用处， 有的图像的锐度很高，用这样的均值算法，可以把锐度降低。使得图像看上去更加自然，下面就有几幅图我们可以看出一些端倪： 原图：                                                                          原图均值滤波处理之后：    这里还是可以明显的感觉到不同的， 没有好坏之分，就是第二幅图片看上去更为平滑。 继续我们的问题， 那这里均值平滑是否具有去除噪声的功能呢？ 我们搞来了椒盐噪声（就是随机的白点，黑点）来试试手： 噪声图（5%椒盐噪声）：                         均值滤波平滑处理之后：    首先这里的噪声还是比较小的， 只有5%，从均值的效果来看的话， 我可以说几乎没有用，其实直观的想也可以判断， 因为这里的处理并没有剔除这些噪声点， 而只是微弱地降低了噪声，所以效果可以想见的。。 好吧， 最后的时候还是贴上一段处理的C++均值滤波算法代码： ~~~ void meanFilter (unsigned char* corrupted, unsigned char* smooth, int width, int height) { memcpy ( smooth, corrupted, width*height*sizeof(unsigned char) ); //复制数据 for (int j=1;j ';

以另外一个滤波器而言----均值滤波器, 就是说某像素的颜色, 由以其为中心的九宫格的像素平均值来决定. 在这个基础上又发展成了带权的“平均”滤波器, 这里的高斯平滑或者说滤波器就是这样一种带权（通常我们认为距离要代替的点像素的作用大一些）的“平均”滤波器. 那么这些权重如何分布呢? 我们先来看几个经典的模板例子:   尝试了使用这些滤波器对我们原来的图进行操作, 得到了这样的一组结果: 原图:  3x3 高斯滤波处理后:  5x5 高斯处理后:    单纯从效果来看, 两个模板都起到了平滑的作用, 只是程度有深浅的区分. 那么从理论上来说为什么能起到平滑的作用呢? 很显然, 像素的颜色不仅由自身决定了, 同时有其周围的像素加权决定, 客观上减小了和周围像素的差异.同时这些权重的设定满足了越近权重越大的规律. 从理论来讲, 这些权重的分布满足了著名的所谓高斯分布: 这就是1维的计算公式：      这就是2维的计算公式：    x, y表示的就是当前点到对应点的距离, 而那些具体的模板就是由这里公式中的一些特例计算而来. 需要说明的是不只有这么一些特例, 从wikipedia可以方便地找到那些复杂的模板比如像: ### Sample Gaussian matrix This is a sample matrix, produced by sampling the Gaussian filter kernel (with σ = 0.84089642) at the midpoints of each pixel and then normalising. Note that the center element (at [4, 4]) has the largest value, decreasing symmetrically as distance from the center increases. | 0.00000067 | 0.00002292 | **0.00019117** | 0.00038771 | **0.00019117** | 0.00002292 | 0.00000067 | |-----|-----|-----|-----|-----|-----|-----| | 0.00002292 | 0.00078633 | 0.00655965 | 0.01330373 | 0.00655965 | 0.00078633 | 0.00002292 | | **0.00019117** | 0.00655965 | 0.05472157 | 0.11098164 | 0.05472157 | 0.00655965 | **0.00019117** | | 0.00038771 | 0.01330373 | 0.11098164 | **0.22508352** | 0.11098164 | 0.01330373 | 0.00038771 | | **0.00019117** | 0.00655965 | 0.05472157 | 0.11098164 | 0.05472157 | 0.00655965 | **0.00019117** | | 0.00002292 | 0.00078633 | 0.00655965 | 0.01330373 | 0.00655965 | 0.00078633 | 0.00002292 | | 0.00000067 | 0.00002292 | **0.00019117** | 0.00038771 | **0.00019117** | 0.00002292 | 0.00000067 | 是不是看到就头大了？不过没关系, 对于一般的应用来说, 前面的例子已经可以完成任务了.  代码的话我们还是给一份5x5的example: ~~~ void gaussianFilter2 (unsigned char* corrupted, unsigned char* smooth, int width, int height) { int templates[25] = { 1, 4, 7, 4, 1, 4, 16, 26, 16, 4, 7, 26, 41, 26, 7, 4, 16, 26, 16, 4, 1, 4, 7, 4, 1 }; //滤波器模板 memcpy ( smooth, corrupted, width*height*sizeof(unsigned char) ); //复制像素 for (int j=2;j 255) sum = 255; smooth [ j*width+i ] = sum; } } } ~~~ 参考资源： 【1】http://blog.csdn.net/hhygcy/article/details/4329056

> 原文出处：[图像处理算法](http://blog.csdn.net/column/details/imageproebow.html) 作者：[ebowtang](http://blog.csdn.net/ebowtang) **本系列文章经作者授权在看云整理发布，未经作者允许，请勿转载！** # 图像处理算法 > 记录各种图像算法以及相关的代码实现。