matlab練習程序（共生矩陣）

紋理描述的共生矩陣方法是基於在紋理中某一灰度級結構重復出現的情況；這個結構在精細紋理中隨着距離而快速的變化，而在粗糙紋理中則緩慢的變化。

假設待分析在精細紋理圖像的一部分是一個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