几种常见空间滤波器MATLAB实现

%******************************************************%
%目的：比较几种不同均值滤波器和自适应中值滤波器的去噪声效果
%日期：2018.5.18
%******************************************************%

%读入图像
img = imread(‘img.tif‘);
figure;
subplot(1,3,1);
imshow(img);
title(‘原始图像‘);
%分别添加加入高斯噪声和椒盐噪声
f1 = imnoise(img,‘gaussian‘);
subplot(1,3,2)
imshow(f1);
title(‘添加高斯噪声‘);
f2 = imnoise(f1,‘salt & pepper‘);
subplot(1,3,3)
imshow(f2);
title(‘再添加椒盐噪声‘);
f2 = double(f2);
[M,N] = size(f2);

%%%%比较几种均值滤波器及中值滤波器的滤波结果%%%%%%
%算术均值滤波器
n = 5; %模板大小
f11 = f2;
for i = 1:M-n+1
for j = 1:N-n+1
g1 = f2(i:i+n-1,j:j+n-1);
s1 = sum(sum(g1));
%中心点的值用子图像的算术均值代替
f11(i+(n-1)/2,j+(n-1)/2) = s1/(n*n); 
end
end
img1 = uint8(f11);
figure;
subplot(1,2,1);
imshow(img1);
title(‘算术均值滤波结果‘);

%几何均值滤波器
n = 5; %模板大小
f21 = f2 + ones(M,N); %防止有个像素点为0而导致乘积为0
f12 = f21;
for i = 1:M-n+1
for j = 1:N-n+1
g2 = f21(i:i+n-1,j:j+n-1);
s2 = prod(g2(:));
%中心点的值用子图像的几何均值代替
f12(i+(n-1)/2,j+(n-1)/2) = s2.^(1/numel(g2)); 
end
end
img2 = uint8(f12);
subplot(1,2,2);
imshow(img2);
title(‘几何均值滤波结果‘);

%中值滤波器
img3 = uint8(medfilt2(f2,[5 5]));
figure;
subplot(1,2,1);
imshow(img3);
title(‘中值滤波结果‘);

%修正的阿尔法均值滤波器
n = 5; %模板大小
d = 5; %去掉最值灰度值的个数
f13 = f2;
for i = 1:M-n+1
for j = 1:N-n+1
g3 = f2(i:i+n-1,j:j+n-1);
g3 = sort(g3(:)); %对邻域内的像素点进行排序
min_num = ceil(d/2); %去掉最小灰度值的个数
max_num = floor(d/2); %去掉最大灰度值的个数
%去掉d个最值灰度级后求算术均值
g3(1:min_num) = zeros(min_num,1);
g3(n-max_num+1:n) = zeros(max_num,1);
s3 = sum(g3);
%中心点的值用子图像的算术均值代替
f13(i+(n-1)/2,j+(n-1)/2) = s3/(n*n-d); 
end
end
img4 = uint8(f13);
subplot(1,2,2);
imshow(img4);
title(‘修正的阿尔法均值滤波结果‘);


%%%%%%%比较中值滤波器及自适应中值滤波器的滤波结果%%%%%%%
f2 = imnoise(img,‘salt & pepper‘);
figure;
subplot(1,4,1);
imshow(img);
title(‘原图‘);
subplot(1,4,2);
imshow(f2);
title(‘添加椒盐噪声的图像‘);

%中值滤波器
img5 = uint8(medfilt2(f2,[5 5]));
subplot(1,4,3);
imshow(img5);
title(‘中值滤波结果‘);

%自适应中值滤波器
f14 = f2; 
alreadyPro = false(size(f2)); %看是否完成进程
Smax=5; %最大窗口尺寸
for s = 3:2:Smax %起始窗口尺寸设为3
%得到特定的灰度值
zmin = ordfilt2(f2, 1, ones(s, s), ‘symmetric‘); 
zmax = ordfilt2(f2, s * s, ones(s, s), ‘symmetric‘); 
zmed = medfilt2(f2, [s s], ‘symmetric‘); 
%进程B
processB = (zmed > zmin) & (zmax > zmed) & ~alreadyPro; 
%进程A
processA = (f2 > zmin) & (zmax > f2); 
outZxy = processB & processA; 
outZmed = processB & ~processA; 
f14(outZxy) = f2(outZxy); 
f14(outZmed) = zmed(outZmed); 
alreadyPro = alreadyPro | processB; 
if all(alreadyPro(:)) 
break; 
end 
end 
f14(~alreadyPro) = zmed(~alreadyPro); 
img6 = uint8(f14);
subplot(1,4,4);
imshow(img6);
title(‘自适应中值滤波结果‘);