Python计算灰度共生矩阵和常用的纹理因子

前言

如何使用skimage计算GLCM和纹理因子可以参考利用python的skimage计算灰度共生矩阵。但是Skimage库只提供了常用的8中纹理因子（均值、方差、同质性、对比度、差异性、熵、角二阶矩、相关性、）中的5种，缺少均值、方差和熵的计算。这里讲解如何修改源码添加这三个因子的计算。添加后调用方式和其他因子的计算相同。

修改/添加部分

找到 skimage库greycoprops源码。Spyder可以直接右键该函数，转到定义处。
需要添加/修改的部分，保持了与源码相同的编码风格，全是numpy向量化计算。

...............
    elif prop in ['ASM', 'energy', 'correlation', 'mean', 'variance', 'entroy']:
...............
    elif prop == 'mean':
        results = np.apply_over_axes(np.mean, P, axes=(0, 1))[0, 0]
    elif prop == 'variance':
        results = np.apply_over_axes(np.sum,
                             (P - np.apply_over_axes(np.mean, P, axes=(0, 1)))**2,
                             axes=(0, 1))[0, 0]
    elif prop == 'entroy':
        with np.errstate(divide='ignore'):
            condition_list = [P != 0, P == 0]
            choice_list = [-np.log10(P), 0]
            P = np.select(condition_list, choice_list)
        results = np.apply_over_axes(np.sum, P, axes=(0, 1))[0, 0]
..............

完整源码

def greycoprops(P, prop='contrast'):
    """Calculate texture properties of a GLCM.

    Compute a feature of a grey level co-occurrence matrix to serve as
    a compact summary of the matrix. The properties are computed as
    follows:

    - 'contrast': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}(i-j)^2`
    - 'dissimilarity': :math:`\\sum_{i,j=0}^{levels-1}P_{i,j}|i-j|`
    - 'homogeneity': :math:`\\sum_{i,j=0}^{levels-1}\\frac{P_{i,j}}{1+(i-j)^2}`
    - 'ASM': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}^2`
    - 'energy': :math:`\\sqrt{ASM}`
    - 'correlation':
        .. math:: \\sum_{i,j=0}^{levels-1} P_{i,j}\\left[\\frac{(i-\\mu_i) \\
                  (j-\\mu_j)}{\\sqrt{(\\sigma_i^2)(\\sigma_j^2)}}\\right]

    Each GLCM is normalized to have a sum of 1 before the computation of texture
    properties.

    Parameters
    ----------
    P : ndarray
        Input array. `P` is the grey-level co-occurrence histogram
        for which to compute the specified property. The value
        `P[i,j,d,theta]` is the number of times that grey-level j
        occurs at a distance d and at an angle theta from
        grey-level i.
    prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \
            'correlation', 'ASM'}, optional
        The property of the GLCM to compute. The default is 'contrast'.

    Returns
    -------
    results : 2-D ndarray
        2-dimensional array. `results[d, a]` is the property 'prop' for
        the d'th distance and the a'th angle.

    References
    ----------
    .. [1] The GLCM Tutorial Home Page,
           http://www.fp.ucalgary.ca/mhallbey/tutorial.htm

    Examples
    --------
    Compute the contrast for GLCMs with distances [1, 2] and angles
    [0 degrees, 90 degrees]

    >>> image = np.array([[0, 0, 1, 1],
    ...                   [0, 0, 1, 1],
    ...                   [0, 2, 2, 2],
    ...                   [2, 2, 3, 3]], dtype=np.uint8)
    >>> g = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4,
    ...                  normed=True, symmetric=True)
    >>> contrast = greycoprops(g, 'contrast')
    >>> contrast
    array([[0.58333333, 1.        ],
           [1.25      , 2.75      ]])

    """
    check_nD(P, 4, 'P')

    (num_level, num_level2, num_dist, num_angle) = P.shape
    if num_level != num_level2:
        raise ValueError('num_level and num_level2 must be equal.')
    if num_dist <= 0:
        raise ValueError('num_dist must be positive.')
    if num_angle <= 0:
        raise ValueError('num_angle must be positive.')

    # normalize each GLCM
    P = P.astype(np.float64)
    glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1))
    glcm_sums[glcm_sums == 0] = 1
    P /= glcm_sums

    # create weights for specified property
    I, J = np.ogrid[0:num_level, 0:num_level]
    if prop == 'contrast':
        weights = (I - J) ** 2
    elif prop == 'dissimilarity':
        weights = np.abs(I - J)
    elif prop == 'homogeneity':
        weights = 1. / (1. + (I - J) ** 2)
    elif prop in ['ASM', 'energy', 'correlation', 'mean', 'variance', 'entroy']:
        pass
    else:
        raise ValueError('%s is an invalid property' % (prop))

    # compute property for each GLCM
    if prop == 'energy':
        asm = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
        results = np.sqrt(asm)
    elif prop == 'ASM':
        results = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
    elif prop == 'correlation':
        results = np.zeros((num_dist, num_angle), dtype=np.float64)
        I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
        J = np.array(range(num_level)).reshape((1, num_level, 1, 1))
        diff_i = I - np.apply_over_axes(np.sum, (I * P), axes=(0, 1))[0, 0]
        diff_j = J - np.apply_over_axes(np.sum, (J * P), axes=(0, 1))[0, 0]

        std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i) ** 2),
                                           axes=(0, 1))[0, 0])
        std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j) ** 2),
                                           axes=(0, 1))[0, 0])
        cov = np.apply_over_axes(np.sum, (P * (diff_i * diff_j)),
                                 axes=(0, 1))[0, 0]

        # handle the special case of standard deviations near zero
        mask_0 = std_i < 1e-15
        mask_0[std_j < 1e-15] = True
        results[mask_0] = 1

        # handle the standard case
        mask_1 = mask_0 == False
        results[mask_1] = cov[mask_1] / (std_i[mask_1] * std_j[mask_1])
    elif prop in ['contrast', 'dissimilarity', 'homogeneity']:
        weights = weights.reshape((num_level, num_level, 1, 1))
        results = np.apply_over_axes(np.sum, (P * weights), axes=(0, 1))[0, 0]
        
    elif prop == 'mean':
        results = np.apply_over_axes(np.mean, P, axes=(0, 1))[0, 0]
    elif prop == 'variance':
        results = np.apply_over_axes(np.sum, 
                             (P - np.apply_over_axes(np.mean, P, axes=(0, 1)))**2, 
                             axes=(0, 1))[0, 0]
    elif prop == 'entroy':
        with np.errstate(divide='ignore'):
            condition_list = [P != 0, P == 0]
            choice_list = [-np.log10(P), 0]
            P = np.select(condition_list, choice_list)
        results = np.apply_over_axes(np.sum, P, axes=(0, 1))[0, 0]
    
    return results

使用示例

分别为方差、均值、熵的计算。

greycoprops(glcm, 'variance')
greycoprops(glcm, 'mean')
greycoprops(glcm, 'entroy')