python共生矩阵统计量 求共生矩阵

纹理描述的共生矩阵方法是基于在纹理中某一灰度级结构重复出现的情况；这个结构在精细纹理中随着距离而快速的变化，而在粗糙纹理中则缓慢的变化。

假设待分析在精细纹理图像的一部分是一个M*N的矩形窗口。某一灰度级结构的出现情况可以由相对的频率

的矩阵来描述，他描绘了具有灰度级a,b的两个像素，在方向Φ上间隔距离为d，以多大的频率出现在窗口中。

共生矩阵求取方法：

这里有两种定义，书上的是第一种，也就是对称矩阵的定义，如下：

 

 

  

 

我程序用的是非对称的矩阵的定义，如下：

 

 

 

 

哦，我想这一大堆公式是不那么容易看明白的，不过下面这个图就让人一目了然了。下图中的Φ=0,d=1.

共生矩阵产生完之后就需要使用其派生出的准则了。

主要有如下六种：

1.能量，或角度二阶矩（图像均匀性的测度——图像越均匀，其值越大）：

2.熵：

3.最大概率：

4.对比度（局部图像变化的测度；典型k=2,λ=1）:

5.倒数差分矩：

6.相关性（图像线性度的测度，在方向Φ上的线性方向的结构在这个方向上会产生大的相关值）：

其中μ_x,μ_y 是均值，以及σ_x，σ_y是标准差，

通用算法如下：

1.给定方向和距离构造共生矩阵。

2.用方向、距离以及六个特征值构成纹理特征向量，这就产生了相关的特征了。

下面是matlab源代码：

main.m

clear all;
close all;
clc;

fai=0;
d=1;

img=imread('lena.jpg');
[m n]=size(img);
k=1;
for i=1:10:m-10     %取10*10的图像，这里我随便写了，准确的是一个挨一个的取10*10
   for j=1:10:n-10 
        re{k}=comatrix(img(i:i+9,j:j+9),fai,d);
        k=k+1;
   end
end

for j=1:8
    s=[];    %存储特征的地方
    for i=1:k-1
        s=[s re{i}(j)];
    end
    figure;
    plot(s)
end

comatrix.m

function re=comatrix(D,fai,d)

    [m n]=size(D);
    ma=max(D(:));
    mi=min(D(:));
    D=D-mi+1;       %将最小值统一为1
    k=ma-mi+1;      %产生合适大小的共生矩阵
    P=zeros(k,k);
    if fai==0
        h=0; w=d;       
    elseif fai==45
        h=d; w=d;
    elseif fai==90
        h=d; w=0;
    else
        h=d; w=-d;
    end

    if fai~=135
        for i=1:m-h
            for j=1:n-w 
                P(D(i,j),D(i+h,j+w))=P(D(i,j),D(i+h,j+w))+1;  
            end
        end
    else
        for i=1:m-h
            for j=1-w:n
                P(D(i,j),D(i+h,j+w))=P(D(i,j),D(i+h,j+w))+1;                
            end
        end
    end
    
    if sum(P(:))~=0
         P=P/sum(P(:));  %归一化矩阵
    end
    
    %re为八个特征值
    re=zeros(1,8);
    re(1)=fai;              %方向特征
    re(2)=d;                %步长特征
    re(3)=sum(P(:).^2);     %能量特征   
    re(4)=max(P(:));        %最大概率特征
    
    miux=sum(double(1:k).*sum(P(1:k,:)));
    miuy=sum(double(1:k).*sum(P(:,1:k)));
    sigmax=sum((double(1:k)-miux).^2.*sum(P(1:k,:)));
    sigmay=sum((double(1:k)-miuy).^2.*sum(P(:,1:k)));
    
    for i=1:k
       for j=1:k
            if P(i,j)~=0
               re(5)=re(5)+P(i,j)*log2(P(i,j)); %熵特征
            end
            re(6)=re(6)+(i-j).^2*P(i,j);        %对比度特征
            if i~=j
                re(7)=re(7)+P(i,j)/abs(i-j);    %倒数差分矩特征
            end
            re(8)=re(8)+i*j*P(i,j);            
       end
    end
    if sigmax*sigmay~=0
        re(8)=(re(8)-miux*miuy)/(sigmax*sigmay);    %相关性特征
    end
end

comatirx.m中的re就是当前取的纹理的特征向量，八个描述子。

最后产生的全是曲线，就不贴图了。

参考：

http://matlab.izmiran.ru/help/toolbox/images/graycomatrix.html