灰度共生矩阵-python

本篇博客参考灰度共生矩阵（GLCM）附Python代码 ，对其进行了完善与修改。

非常感谢评论区 清慎Arthur、jy2lhz、 小小菜鸡升级ing 指出代码中的问题，目前已根据建议进行修改和优化。2019-01-07

灰度共生矩阵（Gray Level Co-occurrence Matrix，GLCM）统计了灰度图中像素间的灰度值分布规律以区分不同的纹理。灰度共生矩阵中每个元素的值可以定义为(x, y)点与(x + dx, y + dy)点的值对为(i, j)的概率。统计整幅图像中每一种灰度值组合出现的概率矩阵 P 即为灰度共生矩阵。

引自百度百科 灰度共生矩阵
灰度级量化
在实际应用中，一幅灰度影像的灰度级一般为256级，在计算由灰度共生矩阵推导出的纹理特征时，要求影像的灰度级远小于256，主要是因为影像共生矩阵的计算量由影像的灰度等级和影像的大小来确定。例如:假定影像G有L个灰度级，其大小为R行C列，则运算量大约是L² * R * C，按一般情况L=256, R=512, C=512来计算，其基本运算至少要1.7 * 10十次方次。以现行微机的运算速度每秒100万次为例，对上述一幅影像计算其灰度共生矩阵至少需要1.7 X 10³秒以上，约30分钟。由此可见，这样长的时间用来进行影像的识别是不太切合实际的。解决的办法是:在尽量保持影像原形的情况下大量削减影像灰度级的取值个数和影像分辨率。因此在计算空间灰度共生矩阵时，在不影响纹理特征的前提下往往先将原影像的灰度级压缩到较小的范围，一般取8级或16级，以便减小共生矩阵的尺寸。
因为灰度级由256变为了8级，所以影像颜色显得很暗淡。但是在进行纹理信息处理时对纹理特征的影响不大。

一张图像的灰度值范围一般为0-255，为减少计算量，在计算矩阵前，需要对灰度值进行量化，一般可设置为8, 16.由于灰度共生矩阵的维度较大，一般不直接作为区分纹理的特征，而是基于它构建的一些统计量作为纹理分类特征。例如Haralick曾提出了14种基于灰度共生矩阵计算出来的统计量：能量、熵、对比度、均匀性、相关性、方差、和平均、和方差、和熵、差方差、差平均、差熵、相关信息测度以及最大相关系数。

# -*- coding: utf-8 -*-
import cv2
import numpy as np
np.set_printoptions(suppress=True)

def glcm(arr, d_x, d_y, gray_level=16):
    '''计算并返回归一化后的灰度共生矩阵'''
    max_gray = arr.max()
    height, width = arr.shape
    arr = arr.astype(np.float64)  # 将uint8类型转换为float64，以免数据失真
    arr = arr * (gray_level - 1) // max_gray  # 若灰度级数大于gray_level，则将图像的灰度级缩小至gray_level，减小灰度共生矩阵的大小。量化后灰度值范围：0 ~ gray_level - 1
    ret = np.zeros([gray_level, gray_level])
    for j in range(height -  abs(d_y)):
        for i in range(width - abs(d_x)):  # range(width - d_x)  #注释为源代码，经评论指出错误后修改
            rows = arr[j][i].astype(int)
            cols = arr[j + d_y][i + d_x].astype(int)
            ret[rows][cols] += 1
    if d_x >= d_y:
        ret = ret / float(height * (width - 1))  # 归一化, 水平方向或垂直方向
    else:
        ret = ret / float((height - 1) * (width - 1))  # 归一化, 45度或135度方向
    return ret

if __name__=='__main__':
    '''归一化时分母值根据角度theta变化，0度或90度时为height * (width - 1), 45度或135度时为(height - 1) * (width - 1)'''
    fp = 'C:/Users/mamq/Pictures/test photos/house_1.jpg'
    img = cv2.imread(fp)
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转换为灰度图像，uint8
    glcm_0 = glcm(img_gray, 1, 0)  # 水平方向
    # glcm_1 = glcm(img_gray, 0, 1)  # 垂直方向
    # glcm_2 = glcm(img_gray, 1, 1)  # 45度方向
    # glcm_3 = glcm(img_gray, -1, 1)  # 135度方向

部分结果

[[0.0084385 0.00192513 0.00074866 0.00050802 0.00043316 0.00040107
0.00056684 0.00029947 0.00019786 0.00020856 0.00006417 0.00003209
0.00003743 0.00005348 0.00004813 0.00003209]
[0.0019893 0.00362032 0.00148128 0.00079679 0.00071658 0.00060428
0.00040107 0.00027807 0.00026738 0.00017112 0.00005882 0.00006417
0.00002139 0.00003743 0.00002139 0.00005882]
……
]

---

以上，欢迎批评指正