1.灰度共生矩陣生成原理

灰度共生矩陣(GLDM)的統計方法是20世紀70年代初由R.Haralick等人提出的，它是在假定影象中各畫素間的空間分佈關係包含了影象紋理資訊的前提下，提出的具有廣泛性的紋理分析方法。

度共生矩陣被定義為從灰度為i的畫素點出發，離開某個固定位置（相隔距離為d，方位為）的點上灰度值為的概率，即，所有估計的值可以表示成一個矩陣的形式，以此被稱為灰度共生矩陣。對於紋理變化緩慢的影象，其灰度共生矩陣對角線上的數值較大；而對於紋理變化較快的影象，其灰度共生矩陣對角線上的數值較小，對角線兩側的值較大。由於灰度共生矩陣的資料量較大，一般不直接作為區分紋理的特徵，而是基於它構建的一些統計量作為紋理分類特徵。Haralick曾提出了14種基於灰度共生矩陣計算出來的統計量：即：能量、熵、對比度、均勻性、相關性、方差、和平均、和方差、和熵、差方差、差平均、差熵、相關資訊測度以及最大相關係數。

本文是借用一篇文章的例子講解灰度共生矩陣，用文字說明感覺說不清，自己之前用該方法做過實驗，還是會忘，所以乾脆用例子的方式介紹，下一次再看也容易理解。

在影象中任意一點（x,y）及偏離它的一點(x+a,y+b)(其中a，b為整數，認為定義)構成點對。設該點對的灰度值為(f1,f2),假設影象的最大灰度級為L,則f1與f2的組合共有L*L種。對於整福影象，統計每一種(f1,f2)值出現的次數，然後排列成一個方陣，再用(f1,f2)出現的總次數將它們歸一化為出現的概率P(f1,f2)，由此產生的矩陣為灰度共生矩陣。

下圖為一個簡單的例子：

圖a為原影象，最大灰度級為16。為表示方便，這裡將灰度級數減小為4級，圖a變為圖b的形式。這樣(f1,f2)取值範圍便為[0,3]。取不同的間隔，將(f1,f2)各種組合出現的次數排列起來，就可得到圖e~g所示的灰度共生矩陣。

圖e表示圖b中(x,y)與偏離它的(x+1,y+0)構成點對時，(f1,f2)取值的情況(填充黃色部分為f1取0，f2取1時的情況，由圖b填充易知共10種)。同理，f,g分別表示圖c，d中(x,y)分別於點(x+1,y+1),(x+2,y+0)構成的點對(f1,f2)出現的情況(圖c填充黃色部分表示f1取0，f2取0時，對角線點對(0,0)出現的情況，共8種：圖d填充黃色部分表示f1取0，f2取2時水平點對(0,2)出現的情況，共9種)。例如，對於a=1，b=0，點對中(0,1)的組合共出現了10次。對比可以看出，(0,1),(1,2),(2,3)和(3,0)均有較高的出現頻數。圖b表明，影象中存在明顯的左上右下方向的紋理。

距離(a,b)的取值不同，灰度共生矩陣中的值不同。a和b的取值要根據紋理週期分佈的特徵來選擇，對於較細的紋理，選取(1,0),(1,1),(2,0)等這樣的值是有必要的。a,b取值較小對應於變化緩慢的紋理影象，其灰度共生矩陣對角線上的數值較大。紋理的變化越快，則對角線上的數值越小，而對角線兩側的值增大。

共生矩陣實際上是兩個畫素點的聯合直方圖，對於影象中細而規則的紋理，成對畫素點的二維直方圖傾向於均勻分佈；對於粗而規則的紋理，則傾向於最對角分佈。

2.灰度共生矩陣特徵量

2.1對比度

度量 矩陣的值是如何分佈和影象中區域性變化的多少，反應了影象的清晰度和紋理的溝紋深淺。紋理的溝紋越深，反差越大，效果越清晰；反之，對比值小，則溝紋淺，效果模糊。

2.2 能量

能量變換反映了影象灰度分佈均勻程度和紋理粗細度。若灰度共生矩陣的元素值相近，則能量較小，表示紋理細緻；若其中一些值大，而其它值小，則能量值較大。能量值大表明一種較均一和規則變化的紋理模式。

2.3 熵

影象包含資訊量的隨機性度量。當共生矩陣中所有值均相等或者畫素值表現出最大的隨機性時，熵最大；因此熵值表明了影象灰度分佈的複雜程度，熵值越大，影象越複雜。

2.4 逆方差

逆方差反映了影象紋理區域性變化的大小，若影象紋理的不同區域間較均勻，變化緩慢，逆方差會較大，反之較小。

2.5相關性

用來度量影象的灰度級在行或列方向上的相似程度，因此值得大小反應了局部灰度相關性，值越大，相關性也越大。

#define GLCM_DIS 3  //灰度共生矩陣的統計距離  
#define GLCM_CLASS 16 //計算灰度共生矩陣的影象灰度值等級化  
typedef enumGLCM_ANGLE
{
GLCM_ANGLE_HORIZATION,
GLCM_ANGLE_VERTICAL,
GLCM_ANGLE_DIGONAL_45,
GLCM_ANGLE_DIGONAL_135
}GLCM_ANGLE;

int CalGlCM(unsigned char* pImage,GLCM_ANGLE angleDirection,double* featureVector)
{  
    int i,j;  
  
    if(NULL == pImage)  
        return 1;  
  
    int * glcm = new int[GLCM_CLASS * GLCM_CLASS];  
    int * histImage = new int[width * height];  
  
    if(NULL == glcm || NULL == histImage)  
        return 2;  
  
    //灰度等級化---分GLCM_CLASS個等級  
    for(i = 0;i < height;i++)
    {  
        for(j = 0;j < width; j++)
	{  
            histImage[i * width + j] = (int)(data[width*i+j] * GLCM_CLASS/256);  
        }  
    }  
  
    //初始化共生矩陣  
    for (i = 0;i < GLCM_CLASS;i++)  
        for (j = 0;j < GLCM_CLASS;j++)  
            glcm[i * GLCM_CLASS + j] = 0;  
  
    //計算灰度共生矩陣  
    int w,k,l;  
    //水平方向  
    if(angleDirection == GLCM_ANGLE_HORIZATION)  
    {  
        for (i = 0;i < height;i++)  
        {  
            for (j = 0;j < width;j++)  
            {  
                l = histImage[i * width + j];  
                if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width)  
                {  
                    k = histImage[i * width + j + GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
                if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width)  
                {  
                    k = histImage[i * width + j - GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
            }  
        }  
    }  
    //垂直方向  
    else if(angleDirection == GLCM_ANGLE_VERTICAL)  
    {  
        for (i = 0;i < height;i++)  
        {  
            for (j = 0;j < width;j++)  
            {  
                l = histImage[i * width + j];  
                if(i + GLCM_DIS >= 0 && i + GLCM_DIS < height)   
                {  
                    k = histImage[(i + GLCM_DIS) * width + j];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
                if(i - GLCM_DIS >= 0 && i - GLCM_DIS < height)   
                {  
                    k = histImage[(i - GLCM_DIS) * width + j];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
            }  
        }  
    }  
    //對角方向  
    else if(angleDirection == GLCM_ANGLE_DIGONAL)  
    {  
        for (i = 0;i < height;i++)  
        {  
            for (j = 0;j < width;j++)  
            {  
                l = histImage[i * width + j];  
  
                if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width && i + GLCM_DIS >= 0 && i + GLCM_DIS < height)  
                {  
                    k = histImage[(i + GLCM_DIS) * width + j + GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
                if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width && i - GLCM_DIS >= 0 && i - GLCM_DIS < height)  
                {  
                    k = histImage[(i - GLCM_DIS) * width + j - GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
            }  
        }  
    }  
  
    //計算特徵值  
    double entropy = 0,energy = 0,contrast = 0,homogenity = 0;  
    for (i = 0;i < GLCM_CLASS;i++)  
    {  
        for (j = 0;j < GLCM_CLASS;j++)  
        {  
            //熵  
            if(glcm[i * GLCM_CLASS + j] > 0)  
                entropy -= glcm[i * GLCM_CLASS + j] * log10(double(glcm[i * GLCM_CLASS + j]));  
            //能量  
            energy += glcm[i * GLCM_CLASS + j] * glcm[i * GLCM_CLASS + j];  
            //對比度  
            contrast += (i - j) * (i - j) * glcm[i * GLCM_CLASS + j];  
            //一致性  
            homogenity += 1.0 / (1 + (i - j) * (i - j)) * glcm[i * GLCM_CLASS + j];  
        }  
    }  
    //返回特徵值  
    i = 0;  
    featureVector[i++] = entropy;  
    featureVector[i++] = energy;  
    featureVector[i++] = contrast;  
    featureVector[i++] = homogenity;  
  
    delete[] glcm;  
    delete[] histImage;  
    return 0;  
}