import numpy as np
import cv2

"""计算灰度共生矩阵"""
def fast_glcm(img, vmin=0, vmax=255, levels=8, kernel_size=5, distance=1.0, angle=0.0):
    '''
    Parameters
    ----------
    img: array_like, shape=(h,w), dtype=np.uint8
        input image
    vmin: int
        minimum value of input image
    vmax: int
        maximum value of input image
    levels: int
        number of grey-levels of GLCM
    kernel_size: int
        Patch size to calculate GLCM around the target pixel
    distance: float
        pixel pair distance offsets [pixel] (1.0, 2.0, and etc.)
    angle: float
        pixel pair angles [degree] (0.0, 30.0, 45.0, 90.0, and etc.)
    Returns
    -------
    Grey-level co-occurrence matrix for each pixels
    shape = (levels, levels, h, w)
    '''
    
    # digitize
    bins = np.linspace(vmin, vmax+1, levels+1)
    gl1 = np.digitize(img, bins) - 1

    # make shifted image
    dx = distance*np.cos(np.deg2rad(angle))
    dy = distance*np.sin(np.deg2rad(-angle))
    mat = np.array([[1.0,0.0,-dx], [0.0,1.0,-dy]], dtype=np.float32)
    h,w = img.shape
    gl2 = cv2.warpAffine(gl1, mat, (w,h), flags=cv2.INTER_NEAREST,
                         borderMode=cv2.BORDER_REPLICATE)

    # make glcm
    glcm = np.zeros((levels, levels, h, w), dtype=np.uint8)
    for i in range(levels):
        for j in range(levels):
            mask = ((gl1==i) & (gl2==j))
            glcm[i,j, mask] = 1

    kernel = np.ones((kernel_size, kernel_size), dtype=np.uint8)
    for i in range(levels):
        for j in range(levels):
            glcm[i,j] = cv2.filter2D(glcm[i,j], -1, kernel)

    glcm = glcm.astype(np.float32)
    return glcm
 
"""计算纹理特征"""
def fast_textural_features(glcm, vmin=0, vmax=255, levels=8, kernel_size=5):
    _,_,h,w = glcm.shape
    # 计算最大值
    max_  = np.max(glcm, axis=(0,1))
    # 计算熵
    pnorm = glcm / np.sum(glcm, axis=(0,1)) + 1./kernel_size**2
    ent  = np.sum(-pnorm * np.log(pnorm), axis=(0,1))
    # 均值/标准差/对比度/相异性/同质性/角二阶距
    mean,std2,cont,diss,homo,asm = [np.zeros((h,w), dtype=np.float32) for x in range(6)]
    for i in range(levels):
        for j in range(levels):
            mean += glcm[i,j] * i / (levels)**2 # 计算均值
            cont += glcm[i,j] * (i-j)**2 # 计算对比度
            diss += glcm[i,j] * np.abs(i-j) # 计算相异性
            homo += glcm[i,j] / (1.+(i-j)**2) # 计算同质性
            asm  += glcm[i,j]**2 # 计算角二阶距
    # 计算能量
    ene = np.sqrt(asm)
    # 计算标准差
    for i in range(levels):
        for j in range(levels):
            std2 += (glcm[i,j] * i - mean)**2
    std = np.sqrt(std2)   
    
    return [max_,ent,mean,std,cont,diss,homo,asm,ene]


if __name__ == '__main__':
    
    img = cv2.imread("input.png",0) # 以灰度读取
    glcm = fast_glcm(img, vmin=0, vmax=255, levels=8, kernel_size=5, distance=1.0, angle=0.0) # 计算共生矩阵
    textural = fast_textural_features(glcm, vmin=0, vmax=255, levels=8, kernel_size=5)