新算法对应的子函数:
主函数:%一种高效快速的高密度椒盐噪声滤除算法 %quickdelhighnoise.m function [R] = quickdelhighnoise(J) [e,f]=size(J); R=zeros(e,f); Nmax=3; %确定最大的滤波半径,那么滤波窗口最大为7 %下面是边界扩展,图像上下左右各增加Nmax像素 z=zeros(e+2*Nmax,f+2*Nmax);%引入一个新矩阵,用以贴上图片 z(Nmax+1:e+Nmax,Nmax+1:f+Nmax)=J;%将图片贴在新矩阵的中心位置 %扩展图像边界 z(1:Nmax,Nmax+1:f+Nmax)=J(1:Nmax,1:f); %扩展上边界 z(1:e+Nmax,f+Nmax+1:f+2*Nmax)=z(1:e+Nmax,f:f+Nmax-1); %扩展右边界 z(e+Nmax+1:e+2*Nmax,Nmax+1:f+2*Nmax)=z(e:e+Nmax-1,Nmax+1:f+2*Nmax); %扩展下边界 z(1:e+2*Nmax,1:Nmax)=z(1:e+2*Nmax,Nmax+1:2*Nmax); %扩展左边界 %z即为一个扩展后的图像 y=z;%y复制z,y为最后输出的图像 %%%%%1 噪声检测阶段 F=zeros(e+2*Nmax,f+2*Nmax);%噪声标记矩阵,初始所有像素默认为信号点,标记为0 for i=Nmax+1:e+Nmax for j=Nmax+1:f+Nmax %检测窗口开始用3*3 if z(i,j) == 0 | z(i,j) == 255 %疑似噪声点,需要进一步判断 F(i,j) = 1; end end end %%%%%% 2 噪声滤除阶段 for i=Nmax+1:e+Nmax for j=Nmax+1:f+Nmax %初始滤波窗口为3*3 if F(i,j) == 1 %噪声点,需要进行滤波 N=sum(sum(F(i-1:i+1,j-1:j+1) == 0));%统计高密度区域内信号点的个数 %N=sum(sum( F(i-1:i+1,j-1:j+1) == 0 & z(i-1:i+1,j-1:j+1) ~= max(z1) & z(i-1:i+1,j-1:j+1) ~= min(z1) )) if N > 0 z6 = zeros(1,N);%放置信号点的矩阵 %存放信号点矩阵 for i1 = i-1:i+1 for j1 = j-1:j+1 if F(i1,j1) == 0 %if F(i1,j1) == 0 & z(i1,j1) ~= max(z1) & z(i1,j1) ~= min(z1) z6(1,N) = z(i1,j1); N = N-1; end end end %z6内存放着信号点灰度值 %for i1 = 1:sum(sum(F(i-1:i+1,j-1:j+1) == 0)) %for i1 = 1:sum(sum( F(i-1:i+1,j-1:j+1) == 0 & z(i-1:i+1,j-1:j+1) ~= max(z1) & z(i-1:i+1,j-1:j+1) ~= min(z1) )) % sum1 = sum1 + z6(1,i1)/( 1 + ( z6(1,i1) - median(z6) )^2 ); % sum2 = sum2 + 1.0/( 1 + ( z6(1,i1) - median(z6) )^2 ); %end y(i,j) = mean(z6); %medy4 =median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] );%已经处理过的4领域系数 %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - medy4 )^2 ) + y(i-1,j)/( 1 + ( y(i-1,j) - medy4 )^2 ) + y(i-1,j+1)/( 1 + ( y(i-1,j+1) - medy4 )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - medy4 )^2 ) ; %y(i,j) =( y(i-1,j-1) + y(i-1,j) + y(i-1,j+1) + y(i,j-1) ) / 4; %y(i,j) = median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] ); %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - z(i,j) )^2 ) + y(i-1,j)/( 1 + ( y(i-1,j) - z(i,j) )^2 ) + y(i-1,j+1)/( 1 + ( y(i-1,j+1) - z(i,j) )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - z(i,j) )^2 ) ; end %F(i,j) = 0;%处理完了噪声点 %z(i,j) = y(i,j); end %处理完该噪声点后,将处理完的新像素用来迭代 end end %第一步滤波结束,第二波开始在预处理的图像上进行滤波 z=y; for i=Nmax+1:e+Nmax for j=Nmax+1:f+Nmax %初始滤波窗口为3*3 if F(i,j) == 1 %噪声点,需要进行滤波 N =0; for i1= i-1:i+1 for j1 =j-1:j+1 if i1 == i & j1 == j%中心像素被排除在外,不进行统计 continue; end if z(i1,j1) ~= 0 & z(i1,j1) ~= 250 N = N + 1;%统计信号点的个数 end end end if N > 0 z6 = zeros(1,N);%放置信号点的矩阵 %存放信号点矩阵 for i1 = i-1:i+1 for j1 = j-1:j+1 if i1 == i & j1 == j%中心像素被排除在外,不进行统计 continue; end if z(i1,j1) ~= 0 & z(i1,j1) ~= 250 %if F(i1,j1) == 0 & z(i1,j1) ~= max(z1) & z(i1,j1) ~= min(z1) z6(1,N) = z(i1,j1); N = N-1; end end end y(i,j) = mean(z6); %medy4 =median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] );%已经处理过的4领域系数 %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - medy4 )^2 ) + y(i-1,j)/( 1 + ( y(i-1,j) - medy4 )^2 ) + y(i-1,j+1)/( 1 + ( y(i-1,j+1) - medy4 )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - medy4 )^2 ) ; %y(i,j) =( y(i-1,j-1) + y(i-1,j) + y(i-1,j+1) + y(i,j-1) ) / 4; %y(i,j) = median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] ); %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - z(i,j) )^2 ) + y(i-1,j)/( 1 + ( y(i-1,j) - z(i,j) )^2 ) + y(i-1,j+1)/( 1 + ( y(i-1,j+1) - z(i,j) )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - z(i,j) )^2 ) ; end %F(i,j) = 0;%处理完了噪声点 %z(i,j) = y(i,j); end %处理完该噪声点后,将处理完的新像素用来迭代 end end R(1:e,1:f) = y( Nmax+1:e+Nmax,Nmax+1:f+Nmax );%拷出滤波后的图像
I=imread('2.png'); I=rgb2gray(I); [e,f]=size(I); J=imnoise(I,'salt & pepper',0.7); %J=imnoise(J,'gaussian',0.002); R=quickdelhighnoise(J); tic;R=quickdelhighnoise(J);toc; subplot(1,2,1),imshow( uint8(I) ),title('原图像'); subplot(1,2,2),imshow( uint8(R) ),title('新算法滤波后的图像'); % 计算三种算法的峰值信噪比 B=8; %编码一个像素用多少二进制位 MAX=2^B-1; %图像有多少灰度级 I=double(I); R=double(R); %%%%%% psnr1= MES1=sum(sum((I-R).^2))/(e*f); %自适应中值去噪的均方差 PSNR1=20*log10(MAX/sqrt(MES1)); %自适应中值算法去噪的峰值信噪比