参考:<<数字图像处理——使用MATLAB分析与实现>> 清华大学出版社 蔡利梅主编
1.图像滤波的分类
空间域滤波:利用模板运算,在原图像上进行处理
频率域滤波:利用正交变换(如小波变换,傅里叶变换),利用噪声点对应高频信息进行滤波
2.噪声的分类
加性噪声:
乘性噪声:
加性噪声一般指热噪声、散弹噪声等,它们与信号的关系是相加,不管有没有信号,噪声都存在。而乘性噪声一般由信道不理想引起,它们与信号的关系是相乘,信号在它在,信号不在他也就不在。
高斯噪声,椒盐噪声
通常我们希望在能够滤除噪声的同时尽可能的保留图象的细节信息
3.空间域滤波
a.均值滤波:即相邻像元取均值
代码:
%自己编写的函数可以将下面的zhongzhi_filter的median改为mean即可
%这里使用库函数
clear,clc;
image = imread("p2.JPG");
noiseI = imnoise(image,"gaussian");
subplot (221), imshow(image) ,title('原图');
subplot(222),imshow( noiseI ) ,title( '高斯噪声图像');
result1 = imfilter(noiseI,fspecial('average',3 )); %3×3均值滤波,提供了边界选项
result2 = imfilter(noiseI,fspecial('average',7)); %7×7均值滤波
subplot(223 ) , imshow(uint8(result1 ) ) , title( '3×3均值滤波');
subplot(224 ) , imshow(uint8(result2 ) ) , title( '7×7均值滤波');
b.高斯滤波:即相邻像元取加权均值
观察权矩阵,距离越小,权越小
代码:
function [result] = Gaussian_filter(image,sigma,r)
%编写高斯滤波函数实现
[ height , width] = size(image);
for x= -r:r
for y= -r:r
H(x+r+1, y+r+1) = 1/(2* pi * sigma ^2).*exp(( -x.^2-y.^2)/(2*sigma^2)); %即权重模板矩阵,距离越远权小
end
end
H = H/sum(H(:)); %H所有值的sum,前面除了,后面就不用除了
result3 = zeros( height, width) ; %最终的结果
midimg = zeros( height+2 * r , width + 2* r); %用来处理边界,边界默认填充0
midimg(r+1:height+r , r+1:width+r) = image;
for ai = r+ 1 :height +r %与权重矩阵相乘求值
for aj= r+ 1 : width + r
temp_row = ai-r; temp_col = aj- r; temp = 0;
for bi= 1:2*r+ 1
for bj=1:2*r+1
temp = temp + (midimg(temp_row + bi - 1, temp_col+ bj-1) *H( bi, bj) ) ;
end
end
result3( temp_row , temp_col) = temp;
end
end
result=uint8(result3);
% figure;imshow( uint8(result3)) ;title( 'myself高斯滤波');
end
Image = imread( 'p1.jpg') ;
Image = im2gray(Image);
sigma1 = 0.6;
sigma2= 10;
r= 3 ;
%高斯模板的参数
NoiseI = imnoise( Image, 'gaussian' ) ;
figure;subplot(221);imshow(NoiseI);title('噪声图');
%添加高斯噪声
gausFilter1 = fspecial( 'gaussian' , [2*r+1 2*r+ 1] , sigma1 ) ;
gausFilter2 = fspecial( 'gaussian' , [2*r+1 2*r+ 1] , sigma2 );
result1 = imfilter(NoiseI , gausFilter1, 'conv') ; %填充选项默认值为0
result2 = imfilter(NoiseI , gausFilter2, 'conv') ;
result3 = Gaussian_filter(Image,sigma1,r);
subplot(222);imshow(result1);title( 'sigma1 = 0.6高斯滤波');
subplot(223);imshow(result2) ; title( 'sigma2= 10高斯滤波');
subplot(224);imshow(result3) ;title( 'myself高斯滤波');
c.中值滤波:即相邻像元取中值
代码:
function [img] = zhongzhi_filter(image,m)
% image为原始图像,二维矩阵,unint8类型
% m为模板大小
n=m;
[height,width] = size(image);
x1 = double(image); %不损失精度,原始的
x2 = x1; %待改变的
for i = 1:height-n+1 %边缘的值无法改变
for j = 1:width-n+1
mb = x1(i:(i+n-1),j:(j+n-1)); %得到n*n的矩阵
mb = mb(:); %使二维数组变成一个列向量
mm = median(mb); %换成mean就是取中值
if abs(x2(i+(n-1)/2,j+(n-1)/2) - mm) >10
x2(i+(n-1)/2,j+(n-1)/2) = mm; %偶数没有中间值
end
end
end
img = uint8(x2);
end
f=imread('p2.JPG');
figure;subplot(221);imshow(f);title("原图像");
J = imnoise(f,"salt & pepper",0.05);
J = im2gray(J); %将三个rgb书改为一个灰度值,以一定的权重
subplot(222);imshow(J);title("添加噪声");
g = zhongzhi_filter(J,9);
subplot(223);imshow(g);title("平滑后图像1");
g = medfilt2(J,[9,9]); #库函数
subplot(224);imshow(g);title("平滑后图像2");
d.双边滤波:
代码:
function [result] = shuangbian_filter(image,sigma_s,sigma_r,w)
[X,Y] = meshgrid(-w:w,-w:w);%定义双边滤波窗口宽度,X按行复制,Y按列复制
Gs = exp(-(X.^2+Y.^2)/(2 * sigma_s ^2));%计算邻域内的空间权值
[hm , wn] = size(image) ;
result = zeros( hm, wn);
for i = 1 : hm
for j= 1 :wn
temp = image( max( i- w,1 ) : min(i+ w , hm ) , max(j- w,1 ) : min(j+ w , wn) );%滤波窗口
Gr = exp( - ( temp - image(i,j)).^2/(2 * sigma_r ^2) );%计算灰度邻近权值,每一点与其滤波窗口内的像素值作差,灰度相差越大,权越小
W = Gr.* Gs((max( i- w,1 ):min( i+ w , hm ) ) - i+w+1,...
(max(j- w, 1 ):min(j + w, wn)) -j+w+1);%w为Gs 和Gr的乘积,GS采用的是-r到r的邻域
result(i,j) = sum (W( : ).* temp( : )) / sum(W( : ) ) ;
end
end
% imshow(result),title( '双边滤波图像');
end
Image = im2double( im2gray(imread( 'p2.JPG')) );
Noise = Image + 0.05 * randn( size( Image)) ;
w = 15; sigma_s = 6; sigma_r= 0.1;
result = shuangbian_filter(Noise,sigma_s,sigma_r,w);
imshow(result),title( '双边滤波图像');
总结:
4.频率域滤波
a.理想低通滤波
类似于阈值法
Image = im2gray(imread( 'face2.jpg')) ;
FImage = fftshift(fft2( double( Image)) ); %里面的为低频
%傅里叶变换及频谱搬移
[N,M] = size(FImage);
g = zeros(N, M);
r1 = floor(M/2);
r2 = floor(N/2);
d0 = [5 11 45 68] ;
for i= 1:4
for x= 1:M %长,即为列
for y = 1 :N %高,即为行
d = sqrt((x-r1)^2+(y- r2)^2); %从中心出发指定半径
if d<= d0(i)
h =1;
else
h = 0;
end
g(y, x)= h * FImage( y,x) ;
end
end
c = real( ifft2( ifftshift(g)) ); %傅里叶逆转换后的函数为复数
figure,imshow(uint8(c)),title(['理想低通滤波DO = ' , num2str( d0(i))]);
end
出现振铃现象
b.巴特沃斯低通滤波器
Image = im2gray(imread( 'face2.jpg'));
Image = imnoise( Image, 'gaussian');
FImage = fftshift(fft2( double( Image))) ;
figure;imshow(Image),title("噪声图");
%傅里叶变换及频谱搬移
[N,M] = size(FImage);
g = zeros(N, M);
r1 = floor(M/2); r2= floor(N/2);d0= 30;
n =[1 2 6 7];
for i= 1:4
for x=1:M
for y= 1:N
d = sqrt( (x-r1)^2+(y-r2)^2); %直接处理
h=1/(1+ (d/d0)^ (2* n(i) ) );
g( y,x) = h * FImage( y,x);
end
end
g= ifftshift(g);
g = real(ifft2(g));
figure,imshow(uint8(g) ) ,title( ['Butterworth低通滤波n = ' , num2str(n(i))]);
end
c.指数低通滤波器
Image = im2gray(imread( 'face2.jpg'));
Image = imnoise( Image, 'gaussian' );
FImage = fftshift(fft2( double( Image) ) );
%傅里叶变换及频谱搬移
[N,M] = size(FImage);
g = zeros(N, M);
r1 = floor(M/2); r2= floor(N/2) ;
d0=[20 40 60];
n = 2;
for i= 1:length(d0)
for x = 1:M
for y= 1:N
d = sqrt((x -r1)^2+(y- r2)^2);
h = exp(-0.5 * (d/d0( i))^n);
g(y, x) = h * FImage(y, x);
end
end
g = ifftshift(g);g = real(ifft2(g));
figure,imshow( uint8(g)) ,title([ '指数低通滤波DO = ' , num2str(d0(i))]);
end
没有出现振铃现象
d.梯形滤波器
Image = im2gray(imread( 'face2.jpg'));
Image = imnoise( Image, 'gaussian' );
%加人噪声
FImage = fftshift( fft2 ( double( Image)));
%傅里叶变换及频谱搬移
[N,M] = size(FImage);
g = zeros(N, M);
r1 = floor(M/2) ; r2 = floor(N/2);
d0= [5 30] ; d1 = [45 70];
for i= 1:2
for x = 1:M
for y= 1:N
d = sqrt( (x- r1)^2+(y- r2)^2 ); %对每一个(x,y)
if d>d1
h= 0 ;
else
if d> d0
h= (d- d1 )/(d0 - d1 );
else
h = 1;
end
end
g(y,x) = h * FImage( y,x);
end
end
g = ifftshift(g);g = real( ifft2(g));
figure,imshow( uint8(g)),title(['梯度低通滤波D0= ', ...
num2str( d0(i)),' ,D1 = ' , num2str( d1 ( i))]);
end
总结:
噪声对应于高频,在傅里叶变换后,通过某一函数与频率域相乘,得到变换后的频率域图象,之后通过傅里叶逆变换复原原图像