平滑算法1:
%roburstmedfilt1.m
function [y]=roburstmedfilt1(x,a,b)
x=double(x);
y=x;
[m,n]=size(x);
k=floor(a/2)+1;
for i=k:m+1-k
for j=k:n+1-k
z=zeros(1,4*k^2-4*k);
num=1;
for i1=i-k+1:i+k-1
for j1=j-k+1:j+k-1
if i1==i & j1==j
continue;
end
z(1,num)=x(i1,j1);
num=num+1;
end
end
%z即为不含中心像素的行向量
zmax=max(z);
zmin=min(z);
%中心像素为噪声点
if x(i,j) > zmax
y(i,j) = zmax;
end
%中心像素为噪声点
if x(i,j) < zmin
y(i,j) = zmin;
end
if x(i,j) < zmax & x(i,j) > zmin
y(i,j) = x(i,j);
end
end
end
平滑算法2:
%roburstmedfilt2.m
function [y]=roburstmedfilt2(x,a,b)
x=double(x);
y=x;
[m,n]=size(x);
k=floor(a/2)+1;
for i=k:m+1-k
for j=k:n+1-k
num=0;
num1=1;
%统计窗口范围内和中心灰度值不相等的像素个数
for i1=i-k+1:i+k-1
for j1=j-k+1:j+k-1
if x(i1,j1)~= x(i,j);
num=num+1;
end
end
end
z=zeros(1,num);%来存储与中心像素不等的行向量,共计num个
for i1=i-k+1:i+k-1
for j1=j-k+1:j+k-1
if x(i1,j1)~= x(i,j);
z(1,num1) = x(i1,j1);
num1=num1+1;
end
end
end
%z即为与中心像素值不等的像素构成的行向量
zmax=max(z);
zmin=min(z);
zmed=median(z);
if x(i,j) > zmax
if zmed < zmax & zmed > zmin%检验中值是否是脉冲
y(i,j) = zmed;
else
y(i,j) = zmax;
end
end
if x(i,j) < zmin
if zmed < zmax & zmed > zmin
y(i,j) = zmed;
else
y(i,j) = zmin;
end
end
if x(i,j) < zmax & x(i,j) > zmin
y(i,j) = x(i,j);
end
end
end
主函数:
I=imread('1.png');
I=rgb2gray(I);
[e f]=size(I);
J=imnoise(I,'salt & pepper',0.7);
%J=imnoise(I,'gaussian',0,0.0002);
%%第一种算法在保持边缘和细节较第二种算法要差一点
l=roburstmedfilt1(J,3,3);%第一种稳健平滑滤波器
l1=roburstmedfilt2(J,3,3);%第二种稳健平滑滤波器
l2=commonfilt2(J,3,3);%普通的中值滤波器
subplot(131),imshow( uint8(l) ),title('第一种稳健平滑滤波器');
subplot(132),imshow( uint8(l1) ),title('第二种稳健平滑滤波器');
subplot(133),imshow( uint8(l2) ),title('普通的中值滤波器');
% 计算三种算法的峰值信噪比
B=8; %编码一个像素用多少二进制位
MAX=2^B-1; %图像有多少灰度级
I=double(I);
l=double(l);
l1=double(l1);
l2=double(l2);
%%%%%% psnr1=
MES1=sum(sum((I-l).^2))/(e*f); %自适应中值去噪的均方差
PSNR1=20*log10(MAX/sqrt(MES1)); %自适应中值算法去噪的峰值信噪比
%%%%%% psnr2=
MES2=sum(sum((I-l1).^2))/(e*f); %自己编写的中值去噪的均方差
PSNR2=20*log10(MAX/sqrt(MES2)); %自己编写的中值去噪的峰值信噪比
%%%%%% psnr3=
MES3=sum(sum((I-l2).^2))/(e*f); %自己编写的中值去噪的均方差
PSNR3=20*log10(MAX/sqrt(MES3)); %自己编写的中值去噪的峰值信噪比