纹理——灰度共生矩阵公式及代码

灰度共生矩阵（Grey-Level Co-occurence Matrix，GLCM）被定义为从灰度为i的像素点出发，离开某个固定位置（相隔距离为d，方位为）的点上灰度值为的概率，即，所有估计的值可以表示成一个矩阵的形式，以此被称为灰度共生矩阵。

• 在图像中任意取一点(x, y)以及偏离它的另一点(x+a, y+b)，形成一个点对，设该点对的 灰度值为(i, j)，即点
(x, y)的灰度为i , 点(x+a, y+b)的灰度为j。
• 固定a 与b，令点(x, y)在整幅图像上移动，会得到各种(i, j)值。设灰度值的级数为L，则i 与j 的组合共有L2种。
• 在整幅图像中，统计出每一种(i, j)值出现的次数，再将它们归一化为出现的概率Pij ，则称方阵[Pij ]L× L 为灰度联合概率矩阵，也称为灰度共生矩阵。
• 灰度联合概率矩阵实际上就是两像素点的联合直方图距离差分值(a, b)取不同的数值组合，可以得到沿一定方向如 (0o，45o，90o，135o) 相隔一定距离d 的像元之间灰度联合概率矩阵。
• a 和b 的取值要根据纹理周期分布的特征来选择，对于较细的纹理，选取(1，0)，(0，1)，(1，1)，(-1，1)等小的差分值。
• 当a 与b 取值较小时，对应于变化较慢的纹理图像(粗纹理)，其灰度联合概率矩阵对角线上的数值较大，倾向于作对角线分布；若纹理的变化越快(细纹理)，则对角线上的数值越小，而对角线两侧上的元素值增大，倾向于均匀分布。

统计量：Haralick曾提出14种统计量：即：能量、熵、对比度、均匀性、相关性、方差、和平均、和方差、和熵、差方差、差平均、差熵、相关信息测度以及最大相关系数。（由于灰度共生矩阵的数据量较大，一般不直接作为区分纹理的特征，而是基于它构建的一些统计量作为纹理分类特征。）

C 语言代码块（原文链接：https://blog.csdn.net/guanyuqiu/article/details/53117507）

#define GLCM_DIS 3  //灰度共生矩阵的统计距离  
#define GLCM_CLASS 16 //计算灰度共生矩阵的图像灰度值等级化  
typedef enum	GLCM_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;  
} 

python代码 来源：https://blog.csdn.net/kmsj0x00/article/details/79463376

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import cv2
import math

#定义最大灰度级数
gray_level = 16

def maxGrayLevel(img):
    max_gray_level=0
    (height,width)=img.shape
    print height,width
    for y in range(height):
        for x in range(width):
            if img[y][x] > max_gray_level:
                max_gray_level = img[y][x]
    return max_gray_level+1

def getGlcm(input,d_x,d_y):
    srcdata=input.copy()
    ret=[[0.0 for i in range(gray_level)] for j in range(gray_level)]
    (height,width) = input.shape
    
    max_gray_level=maxGrayLevel(input)
    
    #若灰度级数大于gray_level，则将图像的灰度级缩小至gray_level，减小灰度共生矩阵的大小
    if max_gray_level > gray_level:
        for j in range(height):
            for i in range(width):
                srcdata[j][i] = srcdata[j][i]*gray_level / max_gray_level

    for j in range(height-d_y):
        for i in range(width-d_x):
             rows = srcdata[j][i]
             cols = srcdata[j + d_y][i+d_x]
             ret[rows][cols]+=1.0

    for i in range(gray_level):
        for j in range(gray_level):
            ret[i][j]/=float(height*width)

    return ret

def feature_computer(p):
    Con=0.0
    Eng=0.0
    Asm=0.0
    Idm=0.0
    for i in range(gray_level):
        for j in range(gray_level):
            Con+=(i-j)*(i-j)*p[i][j]
            Asm+=p[i][j]*p[i][j]
            Idm+=p[i][j]/(1+(i-j)*(i-j))
            if p[i][j]>0.0:
                Eng+=p[i][j]*math.log(p[i][j])
    return Asm,Con,-Eng,Idm

def test(image_name):
    img = cv2.imread(image_name)
    try:
        img_shape=img.shape
    except:
        print 'imread error'
        return

    img=cv2.resize(img,(img_shape[1]/2,img_shape[0]/2),interpolation=cv2.INTER_CUBIC)

    img_gray=cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    glcm_0=getGlcm(img_gray, 1,0)
    #glcm_1=getGlcm(src_gray, 0,1)
    #glcm_2=getGlcm(src_gray, 1,1)
    #glcm_3=getGlcm(src_gray, -1,1)

    asm,con,eng,idm=feature_computer(glcm_0)

    return [asm,con,eng,idm]

if __name__=='__main__':
    result = test("test.bmp")
    print(result)