灰度最小方差的均值滤波器

此为灰度最小方差的均值滤波器：在去噪能力上弱于传统的均值、中值滤波，但在保留图像边缘和细节能力方面要强于前者。

子函数：

%该种方法的基本思想是：先使用模板覆盖图像区域，若像素都属于一个区域，则区域内不含边界，则可以放心平滑去噪（高斯用均值，椒盐用中值）；否则像素属于不同部分，即区域含有边界，那我们不平滑去噪，保持边界即可。
%编程的思想是：9个模板找出方差最小的那个模板，这个模板覆盖像素同属于同一区域的概率最高，含有边界的可能性最低，找出其对其平滑去噪即可
function[y]=grayminvariance(x)
[m,n]=size(x);
x=double(x);
y=x;
for i=3:m-2
    for j=3:n-2
    z1=[x(i-1,j+2),x(i,j+2),x(i+1,j+2),x(i-1,j+1),x(i,j+1),x(i+1,j+1),x(i,j)];
    z2=[x(i-1,j-1),x(i,j-1),x(i+1,j-1),x(i-1,j-2),x(i,j-2),x(i+1,j-2),x(i,j)];
    z3=[x(i-2,j+1),x(i-1,j+1),x(i-2,j),x(i-1,j),x(i-2,j-1),x(i-1,j-1),x(i,j)];
    z4=[x(i+1,j+1),x(i+2,j+1),x(i+1,j),x(i+2,j),x(i+1,j-1),x(i+2,j-1),x(i,j)];
    z5=[x(i-2,j+2),x(i-1,j+2),x(i-2,j+1),x(i-1,j+1),x(i,j+1),x(i-1,j),x(i,j)];
    z6=[x(i+1,j+2),x(i+2,j+2),x(i,j+1),x(i+1,j+1),x(i+2,j+1),x(i+1,j),x(i,j)];
    z7=[x(i-1,j),x(i-2,j-1),x(i-1,j-1),x(i,j-1),x(i-2,j-2),x(i-1,j-2),x(i,j)];
    z8=[x(i+1,j),x(i,j-1),x(i+1,j-1),x(i+2,j-1),x(i+1,j-2),x(i+2,j-2),x(i,j)];
    z9=[x(i-1,j+1),x(i,j+1),x(i+1,j+1),x(i-1,j-1),x(i,j-1),x(i+1,j-1),x(i-1,j),x(i,j),x(i+1,j)];
    z=[z1,z2,z3,z4,z5,z6,z7,z8,z9];
    h=[var(z1),var(z2),var(z3),var(z4),var(z5),var(z6),var(z7),var(z8),var(z9)];
    %对9个方差进行冒泡排序法
        for i1=1:1:8
            for j1=1:1:9-i1
                if h(1,j1) > h(1,j1+1)
                    z0=z(1,j1);%z0为临时变量
                    z(1,j1)=z(1,j1+1);
                    z(1,j1+1)=z0;
                end
            end
        end
    %最小方差已经求出，即h(1,1)，对应z(1,1)这个向量，即对应的方差最小的模板
    y(i,j)=mean( z(1,1) );%中心像素获得最小方差模板对应的均值
    %y(i,j)=median( z(1,1) );%中心像素获得最小方差模板对应的中值
    end
end

I=imread(1.jpg);
I=rgb2gray(I);
[e f]=size(I);
%J=imnoise(I,salt & pepper,0.030);
J=imnoise(I,gaussian,0,0.0002);

%%第一种算法在保持边缘和细节较第二种算法要好很多，当然第一种算法在去噪方面不彻底，不如第二种
l=grayminvariance(J);%灰度最小方差的均值滤波器
m=commonfilt2(J,3,3);%自己编写的均或中值滤波方法

subplot(1,2,1),imshow( uint8(l) ),title(灰度最小方差的均值or中值滤波器);
subplot(1,2,2),imshow( uint8(m)),title(普通均值滤波);

%2 计算三种算法的峰值信噪比
B=8;                %编码一个像素用多少二进制位
MAX=2^B-1;          %图像有多少灰度级

I=double(I);
l=double(l);
m=double(m);


%%%%%% psnr2=42.0913
MES2=sum(sum((I-l).^2))/(e*f);     %自适应中值去噪的均方差
PSNR2=20*log10(MAX/sqrt(MES2));           %自适应中值算法去噪的峰值信噪比

%%%%%% psnr3=21.1225
MES3=sum(sum((I-m).^2))/(e*f);     %自己编写的中值去噪的均方差
PSNR3=20*log10(MAX/sqrt(MES3));           %自己编写的中值去噪的峰值信噪比