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.uint8input imagevmin: intminimum value of input imagevmax: intmaximum value of input imagelevels: intnumber of grey-levels of GLCMkernel_size: intPatch size to calculate GLCM around the target pixeldistance: floatpixel pair distance offsets [pixel] (1.0, 2.0, and etc.)angle: floatpixel pair angles [degree] (0.0, 30.0, 45.0, 90.0, and etc.)Returns-------Grey-level co-occurrence matrix for each pixelsshape = (levels, levels, h, w)'''# digitizebins = np.linspace(vmin, vmax+1, levels+1)gl1 = np.digitize(img, bins) - 1# make shifted imagedx = 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.shapegl2 = cv2.warpAffine(gl1, mat, (w,h), flags=cv2.INTER_NEAREST,borderMode=cv2.BORDER_REPLICATE)# make glcmglcm = 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] = 1kernel = 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**2ent  = 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)**2std = 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)