灰度直方图是对图像上单个象素具有某个灰度进行统计的结果，而灰度共生矩阵是对图像上保持某距离的两象素分别具有某灰度的状况进行统计得到的。

取图像（NxN)中任意一点（x,y)及偏离它的另一点 (x+a,y+b),设该点对的灰度值为（g1,g2).令点（x,y)在整个画面上移动，则会得到各种（g1,g2).值，设灰度值的级数为0-7，则（g1,g2).的组合共有64种。对于整个画面，统计出每一种（g1,g2)值出现的次数，然后排列成一个方阵，在用（g1,g2)出现的总次数将它们归一化为出现的概率P(g1,g2)，这样的方阵称为灰度共生矩阵。

距离差分值（a,b)取不同的数值组合，可以得到不同情况下的联合概率矩阵。（a,b) 取值要根据纹理周期分布的特性来选择，对于较细的纹理，选取（1.0)、（1.1)、（2.0)等小的差分值。当a=1,b=0时，像素对是水平的，即0度扫描；当a=0,b=1时，像素对是垂直的，即
90度扫描；当a=1,b=1时，像素对是右对角线的，即45度扫描；当a=-1,b=-1时，像素对是左对角线，即135度扫描。

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

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import cv2
from PIL import Image
from skimage import data
from math import floor, ceil
from skimage.feature import greycomatrix, greycoprops

def main():
    pass

def image_patch(img2, slide_window, h, w):
    image = img2
    window_size = slide_window
    patch = np.zeros((slide_window, slide_window, h, w), dtype=np.uint8)

    for i in range(patch.shape[2]):
        for j in range(patch.shape[3]):
            patch[:, :, i, j] = img2[i: i + slide_window, j: j + slide_window]

    return patch

def calcu_glcm(img, vmin=0, vmax=255, nbit=64, slide_window=5, step=[2], angle=[0]):
    mi, ma = vmin, vmax
    h, w = img.shape[0],img.shape[1]

    # Compressed gray range：vmin: 0-->0, vmax: 256-1 -->nbit-1
    bins = np.linspace(mi, ma + 1, nbit + 1)
    img1 = np.digitize(img, bins) - 1

    # (512, 512) --> (slide_window, slide_window, 512, 512)
    img2 = cv2.copyMakeBorder(img1, floor(slide_window / 2), floor(slide_window / 2)
                              , floor(slide_window / 2), floor(slide_window / 2), cv2.BORDER_REPLICATE)  # 图像扩充

    patch = np.zeros((slide_window, slide_window, h, w), dtype=np.uint8)
    patch = image_patch(img2, slide_window, h, w)

    # Calculate GLCM (5, 5, 512, 512) --> (64, 64, 512, 512)
    # greycomatrix(image, distances, angles, levels=None, symmetric=False, normed=False)
    glcm = np.zeros((nbit, nbit, len(step), len(angle), h, w), dtype=np.uint8)
    for i in range(patch.shape[2]):
        for j in range(patch.shape[3]):
            glcm[:, :, :, :, i, j] = greycomatrix(patch[:, :, i, j], step, angle, levels=nbit)

    return glcm

def calcu_glcm_mean(glcm, nbit=64):
    '''
    calc glcm mean
    '''
    mean = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            mean += glcm[i, j] * i / (nbit) ** 2

    return mean

def calcu_glcm_variance(glcm, nbit=64):
    '''
    calc glcm variance
    '''
    mean = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            mean += glcm[i, j] * i / (nbit) ** 2

    variance = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            variance += glcm[i, j] * (i - mean) ** 2

    return variance


def calcu_glcm_homogeneity(glcm, nbit=64):
    '''
    calc glcm Homogeneity
    '''
    Homogeneity = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            Homogeneity += glcm[i, j] / (1. + (i - j) ** 2)

    return Homogeneity


def calcu_glcm_contrast(glcm, nbit=64):
    '''
    calc glcm contrast
    '''
    contrast = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            contrast += glcm[i, j] * (i - j) ** 2

    return contrast


def calcu_glcm_dissimilarity(glcm, nbit=64):
    '''
    calc glcm dissimilarity
    '''
    dissimilarity = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            dissimilarity += glcm[i, j] * np.abs(i - j)

    return dissimilarity


def calcu_glcm_entropy(glcm, nbit=64):
    '''
    calc glcm entropy
    '''
    eps = 0.00001
    entropy = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            entropy -= glcm[i, j] * np.log10(glcm[i, j] + eps)

    return entropy


def calcu_glcm_energy(glcm, nbit=64):
    '''
    calc glcm energy or second moment
    '''
    energy = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            energy += glcm[i, j] ** 2

    return energy


def calcu_glcm_correlation(glcm, nbit=64):
    '''
    calc glcm correlation (Unverified result)
    '''

    mean = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            mean += glcm[i, j] * i / (nbit) ** 2

    variance = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            variance += glcm[i, j] * (i - mean) ** 2

    correlation = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            correlation += ((i - mean) * (j - mean) * (glcm[i, j] ** 2)) / variance

    return correlation


def calcu_glcm_Auto_correlation(glcm, nbit=64):
    '''
    calc glcm auto correlation
    '''
    Auto_correlation = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
    for i in range(nbit):
        for j in range(nbit):
            Auto_correlation += glcm[i, j] * i * j

    return Auto_correlation

def MatrixToImage(data):
    data = data*255
    new_im = Image.fromarray(data.astype(np.uint8))
    return new_im

if __name__ == '__main__':
    main()

    nbit = 64  # gray levels
    mi, ma = 0, 255  # max gray and min gray
    slide_window = 7  # sliding window
    # step = [2, 4, 8, 16] # step
    # angle = [0, np.pi/4, np.pi/2, np.pi*3/4] # angle or direction
    step = [2]
    angle = [0]
    image = "t2.jpg"                                         #测试图像文件路径
    img = np.array(Image.open(image))
    h, w = img.shape[0],img.shape[1]
    img=img[:,:,0]
    glcm = calcu_glcm(img, mi, ma, nbit, slide_window, step, angle)
    a=glcm.shape[2]
    b=glcm.shape[3]
    for i in range(glcm.shape[2]):
        for j in range(glcm.shape[3]):
            glcm_cut = np.zeros((nbit, nbit, h, w), dtype=np.float32)
            glcm_cut = glcm[:, :, i, j, :, :]
            mean = calcu_glcm_mean(glcm_cut, nbit)   #在灰度共生矩阵的基础上计算特征（均值，对比度，相关度.....）

            #new_im = MatrixToImage(mean)
            #plt.imshow(mean)

            #new_im.show()

            #new_im.save('mean2.jpg')
            variance = calcu_glcm_variance(glcm_cut, nbit)

            #homogeneity = calcu_glcm_homogeneity(glcm_cut, nbit)
            #contrast = calcu_glcm_contrast(glcm_cut, nbit)
            #dissimilarity = calcu_glcm_dissimilarity(glcm_cut, nbit)
            #entropy = calcu_glcm_entropy(glcm_cut, nbit)
            energy = calcu_glcm_energy(glcm_cut, nbit)
            new_im = MatrixToImage(energy)
            plt.imshow(energy)
            new_im.show()
            new_im.save('energy2.jpg')
            #correlation = calcu_glcm_correlation(glcm_cut, nbit)
            #Auto_correlation = calcu_glcm_Auto_correlation(glcm_cut, nbit)

原始图像704*576

均值特征矩阵

方差

能量