平滑滤波方法研究

平滑滤波方法研究平滑滤波是低频增强的空间域滤波技术。

它的目的有两类：一类是模糊；另一类是消除噪音。

并且具有一定的处理要求，一是不能损坏图像的轮廓及边缘等重要信息；二是使图像清晰视觉效果好。

平滑滤波的方法有邻域平滑滤波，就是求邻近像元点的平均亮度值，双边滤波，中值滤波,以及非局部均值滤波等。

1、双边滤波法双边滤波是一种非线性滤波器，它可以达到保持边缘、降噪平滑的效果。

双边滤波的边缘保持特性主要是通过在卷积的过程中组合空域函数和值域核函数来实现的，典型的核函数为高斯分布函数，如下所示：其中：为归一化作用。

σs为空域高斯函数的标准差，σr为值域高斯函数的标准差，Ω表示卷积的定义域。

编写代码测试，当添加的噪声为0.05时，结果如下滤波后图像添加噪声为0.3时，结果如下滤波后图像由此可知，双边滤波具有去除噪音的作用2、邻域平均法邻域平滑滤波原理:邻域平均法就是对含噪声的原始图像f(x,y)的每一个像素点取一个邻域，计算S中所有像素灰度级的平均值，作为邻域平均处理后的图像g(x, y)的像素值。

即式中：x,y=0,1,…,N-1;S是以（x,y）为中心的邻域的集合，M是S 内的点数。

邻域平均法的思想是通过一点和邻域内像素点求平均来去除突变的像素点，从而滤掉一定噪声，其优点是算法简单，计算速度快，其代价会造成图像在一定程度上的模糊。

3、中值滤波法中值滤波就是用一个奇数点的移动窗口，将窗口的中心点的值用窗口内的各点中值代替。

假设窗口内有五点，其值为80、90、200、110和120，那么此窗口内各点的中值及为110。

设有一个一维序列f1,f2,…,fn,取窗口长度（点数）为m（m为奇数）,对其进行中值滤波，就是从输入序列中相继抽出m个数fi-v,…,fi-1,fi,fi+1,…,fi+v(其中fi为窗口中心值,v=(m-1)/2),再将这m个点按其数值大小顺序排序，取其序号的中心点的那个数作为滤波输出。

数学公式表示为：Yi=Med{fi-v,…,fi-1,fi,fi+1,…,fi+v} i∈N v=(m-1)/2 （式1-2）Yi称为序列fi-v,…,fi-1,fi,fi+1,…,fi+v的中值例如，有一序列{0，3，4，0，7}，重新排序后为{0，0，3，4，7}则Med{0，0，3，4，7}=3。

此列若用平滑滤波，窗口也取5，那么平滑滤波输出为（0+3+4+0+7）/5=2.8。

把一个点的特定长度或形状的邻域称作窗口。

在一维情况下，中值滤波器是一个含有奇数个像素的滑动窗口。

中值滤波很容易推广到二维，此时可以利用二维形式的窗口。

对于平面图像采用的二维中值滤波可以由下式表示：式中：A为窗口，{Xij}为二维数据序列,即数字图像各点的灰度值。

在对图像进行中值滤波时，如果窗口是关于中心点对称的，并且包含中心点在内，则中值滤波能保持任意方向的跳变边缘。

图像中的跳变边缘指图像中不同灰度区域之间的灰度突变边缘。

在实际使用窗口时，窗口的尺寸一般先取3，再取5，依次增大，直到滤波效果满意为止，对于有缓变的较长轮廓线物体的图像，采用方形或圆形窗口较合适，对于包含尖顶角物体的图像，采用十字形窗口较合适。

使用二维中值滤波值得注意的是要保持图像中有效线状体。

通过实验可得出中值滤波具有以下特性：（1）对于某些输入信号中值滤波具有不变性。

（2）中值滤波可以用来减弱随机干扰的脉冲干扰，具有较好的去噪声性能。

4、非局部均值（NLM）滤波法非局部均值滤波用于图像去噪，它就充分利用了图像的冗余性。

它的基本思想是图像包含了重复结构，取其平均值将降低噪声。

NLM 滤波是斯拉夫斯基的演变，它是从局部相似强度中取相似图像像素的平均值。

两种滤波主要不同之处是使用区域间比较得出像素间的相似性比像素间比较更具有鲁棒性。

况且匹配模式并未限制在局部区域。

也就是说，远离被过滤的像素不被惩罚。

给出一幅图像Y使用NLM方法在点处的过滤值可看成是计算邻近像素间的加权平均值Ni，公式如下:滤波前后图像对比平滑滤波即消除噪声是图像处理中重要的一个方面。

最基本的平滑方法在空域为均值滤波和中值滤波，频域有各种低通平滑滤波器。

平滑方法又可分为线性与非线性。

总之，不管我们使用什么样的平滑方法，其实就是满足去噪要求后尽最大可能地保持图像的原有特性与细节。

双边滤波clear all;close all;clc;filename = 'G:\蓝春\双边滤波\001.jpg';img=imread(filename);A = im2double(img);Info=imfinfo(filename);pixHeight=Info.Height;pixWidth=Info.Width;if Info.BitDepth>8img=rgb2gray(img);end% img = padarray(img,[3 3],'symmetric','both');% f=imnoise(f0,'salt & pepper',0.05);img=imnoise(img,'salt & pepper',0.1);img=imnoise(img,'gaussian',0.3);% img=mat2gray(img);[m n l]=size(img);% [m n 1] = size(img);dim = size(img);B = zeros(dim);figure(1);imshow(img);r=10;imgn=zeros(m+2*r+1,n+2*r+1);imgn(r+1:m+r,r+1:n+r)=img(:,:,1);imgn(1:r,r+1:n+r)=img(1:r,1:n);imgn(1:m+r,n+r+1:n+2*r+1)=imgn(1:m+r,n:n+r);imgn(m+r+1:m+2*r+1,r+1:n+2*r+1)=imgn(m:m+r,r+1:n+2*r+ 1);imgn(1:m+2*r+1,1:r)=imgn(1:m+2*r+1,r+1:2*r); sigma_d=3;sigma_r=0.1;% [x,y] = meshgrid(-r:r,-r:r);% w1=exp(-(x.^2+y.^2)/(2*sigma_d^2));[X,Y] = meshgrid(-r:r,-r:r);G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));% h = waitbar(0,'Applying bilateral filter...');% set(h,'Name','Bilateral Filter Progress');h=waitbar(0,'wait...');for i = 1:dim(1)for j = 1:dim(2)% Extract local region.iMin = max(i-r,1);iMax = min(i+r,dim(1));jMin = max(j-r,1);jMax = min(j+r,dim(2));I=A(iMin:iMax,jMin:jMax);% Compute Gaussian intensity weights.H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));% Calculate bilateral filter response.F = H.*G((iMin:iMax)-i+r+1,(jMin:jMax)-j+r+1); B(i,j) = sum(F(:).*I(:))/sum(F(:));endwaitbar(i/dim(1));endclose(h)Image_bf = B;imgout=uint8(Image_bf*255);% imgout=s(r+1:m+r,r+1:n+r);imwrite(imgout, 'G:\蓝春\双边滤波\001.jpg');figure(2);imshow(mat2gray(imgout));非局部均值（NLM）滤波% filename='original.bmp';inputpath ='G:\蓝春\双边滤波\';outputpath='G:\蓝春\双边滤波\';% file='128';% file='256';file='chu512';% file='1024';ex='.BMP';filename = [inputpath,file, ex];img=imread(filename);Info=imfinfo(filename);pixHeight=Info.Height;pixWidth=Info.Width;if Info.BitDepth>8img=rgb2gray(img);endfigure(1);imshow(img);tic;restored_img=NLmeansfilter(single(img),21,7,5); timeused = toc;figure(2);imshow(restored_img,[]);title('滤波后图像');imwrite(uint8(restored_img),'restored_img.bmp','bmp') ;function [output]=NLmeansfilter(input,t,f,h)[m n]=size(input);Output=zeros(m,n);input2 = padarray(input,[f f],'symmetric');kernel = make_kernel(f);kernel = kernel / sum(sum(kernel));for i=1:mfor j=1:ni1 = i+ f;j1 = j+ f;W1= input2(i1-f:i1+f , j1-f:j1+f);wmax=0;average=0;sweight=0;rmin = max(i1-t,f+1);rmax = min(i1+t,m+f);smin = max(j1-t,f+1);smax = min(j1+t,n+f);for r=rmin:1:rmaxfor s=smin:1:smaxif(r==i1 && s==j1)continue;end;W2= input2(r-f:r+f , s-f:s+f);d =sum(sum(kernel.*single((W1-W2).*(W1-W2))));w=exp(-d/h);if w>wmaxwmax=w;endsweight = sweight + w;W3 = input2(r,s);average = average + w*single(W3);endendaverage = average +wmax*single(input2(i1,j1)); ¨sweight = sweight + wmax;if sweight > 0output(i,j) = average / sweight;elseoutput(i,j) = input(i,j);endendwaitbar(i/m);endfunction [kernel] = make_kernel(f)kernel=zeros(2*f+1,2*f+1);for d=1:fvalue= 1 / (2*d+1)^2 ;for i=-d:dfor j=-d:dkernel(f+1-i,f+1-j)= kernel(f+1-i,f+1-j) + value ;endendendkernel = kernel ./f。