基于文章“用于图像处理的自适应中值滤波”的matlab代码如下:
%commonfilt2_1.m
%一种自适应调整窗口,一种自适应滤波的算法
function [y]=commonfilt2_1(x)
TD = 9;%判断噪声点所用阈值
N1 = 0; %子块1中噪声点的个数
N2 = 0; %子块2中噪声点的个数
N3 = 0; %子块3中噪声点的个数
N4 = 0; %子块4中噪声点的个数
W1 = 0; %子块1需要滤波的窗口的大小
W2 = 0; %子块2需要滤波的窗口的大小
W3 = 0; %子块3需要滤波的窗口的大小
W4 = 0; %子块4需要滤波的窗口的大小
y=x;
[m,n]=size(x);
F=zeros(m,n);%检查噪声点的矩阵。元素值为1则为噪声点,信号点则为0.初值设为0
%%%%%%%%%%%%%%第一步,检测子块中的噪声,按照文献,我们分图像为4块
%第一步,检测子块1中的噪声点:使用3*3窗口检测
for i=2:255
for j=2:255
sum1 = 0; %小窗口内非极值的像素之和
M1 = 0; %小窗口内非极值的像素个数
z1=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
z1max=max(z1(:));
z1min=min(z1(:));
%找出小窗口内与最值不相等的点的均值
for i1=i-1:i+1
for j1=j-1:j+1
if x(i1,j1) ~= z1max & x(i1,j1) ~= z1min
sum1 = sum1+x(i1,j1);
M1 = M1+1;
end
end
end
T1 = sum1/M1;%窗口内非极值的像素的均值
%下面开始判断噪声点
if( x(i,j) == z1max | x(i,j) == z1min ) & ( abs( x(i,j) - T1) >= TD )%满足该条件的像素是噪声点
F(i,j) = 1; %标记相应的点为噪声点
N1 = N1+1;
end
end
end
%第一步,检测子块2中的噪声点:使用3*3窗口检测
for i=2:255
for j=258:511
sum2 = 0; %小窗口内非极值的像素之和
M2 = 0; %小窗口内非极值的像素个数
z2=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
z2max=max(z2(:));
z2min=min(z2(:));
%找出小窗口内与最值不相等的点的均值
for i1=i-1:i+1
for j1=j-1:j+1
if x(i1,j1) ~= z2max & x(i1,j1) ~= z2min %非最值
sum2 = sum2+x(i1,j1);
M2 = M2+1;
end
end
end
T2 = sum2/M2;%窗口内非极值的像素的均值
%下面开始判断噪声点
if( x(i,j) == z2max | x(i,j) == z2min ) & ( abs( x(i,j) - T2) >= TD )%满足该条件的像素是噪声点
F(i,j) = 1; %标记相应的点为噪声点
N2 = N2+1;
end
end
end
%第一步,检测子块3中的噪声点:使用3*3窗口检测
for i=258:511
for j=2:255
sum3 = 0; %小窗口内非极值的像素之和
M3 = 0; %小窗口内非极值的像素个数
z3=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
z3max=max(z3(:));
z3min=min(z3(:));
%找出小窗口内与最值不相等的点的均值
for i1=i-1:i+1
for j1=j-1:j+1
if x(i1,j1) ~= z3max & x(i1,j1) ~= z3min
sum3 = sum3+x(i1,j1);
M3 = M3+1;
end
end
end
T3 = sum3/M3;%窗口内非极值的像素的均值
%下面开始判断噪声点
if( x(i,j) == z3max | x(i,j) == z3min ) & ( abs( x(i,j) - T3) >= TD )%满足该条件的像素是噪声点
F(i,j) = 1; %标记相应的点为噪声点
N3 = N3+1;
end
end
end
%第一步,检测子块4中的噪声点:使用3*3窗口检测
for i=258:511
for j=258:511
sum4 = 0; %小窗口内非极值的像素之和
M4 = 0; %小窗口内非极值的像素个数
z4=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
z4max=max(z4(:));
z4min=min(z4(:));
%找出小窗口内与最值不相等的点的均值
for i1=i-1:i+1
for j1=j-1:j+1
if x(i1,j1) ~= z4max & x(i1,j1) ~= z4min
sum4 = sum4+x(i1,j1);
M4 = M4+1;
end
end
end
T4 = sum4/M4;%窗口内非极值的像素的均值
%下面开始判断噪声点
if( x(i,j) == z4max | x(i,j) == z4min ) & ( abs( x(i,j) - T4) >= TD )%满足该条件的像素是噪声点
F(i,j) = 1; %标记相应的点为噪声点
N4 = N4+1;
end
end
end
%%%%%%%%%第二步,确定各个子块中滤波器窗口的大小
%%每一子块的大小都是相同的,都是256*256
W=[W1,W2,W3,W4];%存放窗口大小的矩阵,初试值为0
p1 = double(N1/(256*256));
p2 = double(N2/(256*256));
p3 = double(N3/(256*256));
p4 = double(N4/(256*256));
p=[p1,p2,p3,p4];
for i =1:4
if p(1,i) >= 0.002 & p(1,i) < 0.25
W(1,i) = 3;
end
if p(1,i) >= 0.25 & p(1,i) < 0.45
W(1,i) = 5;
end
if p(1,i) >= 0.45
W(1,i) = 7;
end
end
%%%%%%%%%%%第三步,噪声滤波
%对第一个子块滤波
k1=floor(W1/2)+1;
for i=k1:257-k1
for j=k1:257-k1
if F(i,j) == 1 %若是噪声点
z1 = x(i-k1+1:i+k1-1,j);
z11 = median(z1(:));%0°方向上像素的中值
z2 = x(i,j-k1+1:j+k1-1);
z22 = median(z2(:));%90°方向上像素的中值
%处理45°和135°要麻烦一点
k=1;z3=zeros(1,2*k1-1);
for l = 1-k1:k1-1
z3(1,k) = x(i+l,j+l);
k = k+1;
end
z33 =median(z3(:));%45°方向上像素的中值
k=1;z4=zeros(1,2*k1-1);
for l = 1-k1:k1-1
z4(1,k) = x(i+l,j-l);
k = k+1;
end
z44 =median(z4(:));%135°方向上像素的中值
%根据文献,求出4个方向的中值,再乘以相关的系数累加求和赋给中心点
y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
end
end
end
%对第二个子块滤波
k2=floor(W2/2)+1;
for i=k2:257-k2
for j=256+k2:513-k2
if F(i,j) == 1
z1 = x(i-k2+1:i+k2-1,j);
z11 = median(z1(:));%0°方向上像素的中值
z2 = x(i,j-k2+1:j+k2-1);
z22 = median(z2(:));%90°方向上像素的中值
%处理45°和135°要麻烦一点
k=1;z3=zeros(1,2*k2-1);
for l = 1-k2:k2-1
z3(1,k) = x(i+l,j+l);
k = k+1;
end
z33 =median(z3(:));%45°方向上像素的中值
k=1;z4=zeros(1,2*k2-1);
for l = 1-k2:k2-1
z4(1,k) = x(i+l,j-l);
k = k+1;
end
z44 =median(z4(:));%135°方向上像素的中值
%根据文献,求出4个方向的中值,再乘以相关的系数累加求和赋给中心点
y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
end
end
end
%对第三个子块滤波
k3=floor(W3/2)+1;
for i=256+k3:513-k3
for j=k3:257-k3
if F(i,j) == 1
z1 = x(i-k3+1:i+k3-1,j);
z11 = median(z1(:));%0°方向上像素的中值
z2 = x(i,j-k3+1:j+k3-1);
z22 = median(z2(:));%90°方向上像素的中值
%处理45°和135°要麻烦一点
k=1;z3=zeros(1,2*k3-1);
for l = 1-k3:k3-1
z3(1,k) = x(i+l,j+l);
k = k+1;
end
z33 =median(z3(:));%45°方向上像素的中值
k=1;z4=zeros(1,2*k3-1);
for l = 1-k3:k3-1
z4(1,k) = x(i+l,j-l);
k = k+1;
end
z44 =median(z4(:));%135°方向上像素的中值
%根据文献,求出4个方向的中值,再乘以相关的系数累加求和赋给中心点
y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
end
end
end
%对第四个子块滤波
k4=floor(W4/2)+1;
for i=256+k4:513-k4
for j=256+k4:513-k4
if F(i,j) == 1
z1 = x(i-k4+1:i+k4-1,j);
z11 = median(z1(:));%0°方向上像素的中值
z2 = x(i,j-k4+1:j+k4-1);
z22 = median(z2(:));%90°方向上像素的中值
%处理45°和135°要麻烦一点
k=1;z3=zeros(1,2*k4-1);
for l = 1-k4:k4-1
z3(1,k) = x(i+l,j+l);
k = k+1;
end
z33 =median(z3(:));%45°方向上像素的中值
k=1;z4=zeros(1,2*k4-1);
for l = 1-k4:k4-1
z4(1,k) = x(i+l,j-l);
k = k+1;
end
z44 =median(z4(:));%135°方向上像素的中值
%根据文献,求出4个方向的中值,再乘以相关的系数累加求和赋给中心点
y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
end
end
end
%commonfilt2_2.m
%库函数中值滤波
function [y]=commonfilt2_2(x,a,b)
y=x;
[m,n]=size(x);
k=floor(a/2);
for i=k+1:m-k
for j=k+1:n-k
z=x(i-k:i+k,j-k:j+k);
y(i,j)=median(z(:));
end
end
I=imread('1.png');
I=rgb2gray(I);
J=imnoise(I,'salt & pepper',0.4);
k1=commonfilt2_1(J);%一种自适应调整窗口,一种自适应滤波的算法
k2=commonfilt2_2(J,3,3);%库函数中值滤波方法
k3=commonfilt2_2(J,5,5);%库函数中值滤波方法
subplot(222),imshow( uint8(k2)),title('库函数中值滤波3*3');
subplot(224),imshow( uint8(k1) ),title('自适应调整窗口,自适应滤波');
subplot(223),imshow( uint8(k3)),title('库函数中值滤波5*5');
subplot(221),imshow( uint8(J)),title('噪声图像');
% 计算三种算法的峰值信噪比
[e f]=size(I);
B=8; %编码一个像素用多少二进制位
MAX=2^B-1; %图像有多少灰度级
I=double(I);
k1=double(k1);
k2=double(k2);
k3=double(k3);
%%%%%% psnr1=
MES1=sum(sum((I-k1).^2))/(e*f);
PSNR1=20*log10(MAX/sqrt(MES1));
%%%%%% psnr2=
MES2=sum(sum((I-k2).^2))/(e*f);
PSNR2=20*log10(MAX/sqrt(MES2));
%%%%%% psnr3=
MES3=sum(sum((I-k3).^2))/(e*f);
PSNR3=20*log10(MAX/sqrt(MES3));