梯度有关的加权滤波器

代码分两部分：

1 与梯度有关的加权滤波器，设计了2个：1）自适应加权滤波器 2）倒数加权滤波器

2 图像的梯度幅值图像，以及x方向的梯度幅值图像和y方向的梯度幅值图像

I=imread('1.png');
I=rgb2gray(I);
[M,N]=size(I);
I=double(I);
J=I;%权重因子处理后的图像
%k=4;%权重因子系数,数学函数表达式：exp( (-x*x)/(2*k) )，k值的选择也很关键
[Gx,Gy]=gradient(I);
W=sqrt(Gx.*Gx+Gy.*Gy); %求梯度幅值图像
%%%%%%%%%%%%%%%延拓图像来求合适的k值
W1=zeros(2*M,2*N);%新的延拓矩阵：（2M*2N）*（2M*2N）
for i=1:M
  for j=1:N
   W1(i,j) = W(i,j);
  end
end
for i=1:M
    for j=N+1:2*N
     W1(i,j) = -W(i,j-N); 
    end
end
%W1分成三部分：W  -W  0矩阵 
Wmean=mean(mean(W(:)));%W矩阵的均值
W2 =Wmean*ones(2*M,2*N);
W3 =W2-W1;
W3=W3.*W3;
%k^2 =e[ [未延拓图像的灰度梯度幅值均值-延拓后图像的灰度梯度幅值]^2 ]
k1 =sum(sum(W3(:)))/(4*M*N); 
k=sqrt(k1);
%%%%%%%%%%使用（功能类似于高通）滤波器，保留边缘与细节
F=zeros(M,N);%权重因子矩阵，存放计算得出的权重因子
for i=1:M
    for j=1:N
      F(i,j) = exp( double(-W(i,j)^2)/(2*k^2) );%自适应滤波器：效果一般
      %F(i,j) =double(1/W(i,j));%梯度倒数加权滤波器： 效果不错
    end
end
%权重因子矩阵F已经求出，联合原先图像，3*3窗口对中心像素进行处理
for i=2:M-1
    for j=2:N-1
      J(i,j)=I(i,j)*F(i,j)+I(i-1,j+1)*F(i-1,j+1)+I(i,j+1)*F(i,j+1)+I(i+1,j+1)*F(i+1,j+1)+I(i-1,j)*F(i-1,j)+I(i+1,j)*F(i+1,j)+I(i-1,j-1)*F(i-1,j-1)+I(i,j-1)*F(i,j-1)+I(i+1,j-1)*F(i+1,j-1); 
    end
end
%由滤波器的特性，我们发现：
%边缘处梯度值往往较大，F（i，j）=0，J（i，j）也基本等于0，所以处理后的图像边缘为黑
%非边缘处梯度值较小，F（i,j）等于1，J（i，j）基本等于九数之和，所以处理后的图像非边缘为白
subplot(221),imshow(Gx),title('x 方向的梯度图像');
subplot(222),imshow(Gy),title('y 方向的梯度图像');
subplot(223),imshow(uint16(J)),title('权重因子处理后的图像');
subplot(224),imshow(uint8(W)),title('梯度幅值图像');%可以凸出图像的边缘

%%计算一幅图像的平均梯度
%%平均梯度可敏感地反映图像对微小细节反差表达的能力，可用来评价图像的模糊程度
%%以此来衡量图像的清晰度，同时反映出图像中微小细节反差和纹理变换特征
diffX=zeros(M,N);% save the differential value of X orient
diffY=zeros(M,N);% save the differential value of Y orient

tempX=zeros(M,N);
tempY=zeros(M,N);
tempX(1:M,1:(N-1))=I(1:M,2:N);
tempY(1:(M-1),1:N)=I(2:M,1:N);

diffX=I-tempX;
diffY=I-tempY;
diffX(1:M,N)=0;       %%% the boundery set to 0，否则的话diffX最后一列保持与I一样；diffY最后一行与I一样
diffY(M,1:N)=0;
%计算平均梯度:((diffX.*diffX+diffY.*diffY)/2)/(m-1)*(n1)
diffX=diffX.*diffX;
diffY=diffY.*diffY;
diff=double((diffX+diffY))/2;
diff=sqrt(diff);
AVEGRAD=sum(sum(diff)); 
AVEGRAD=AVEGRAD/((M-1)*(N-1));