skimage下的灰度共生矩阵计算源码

在scikit-image中，计算灰度共生矩阵的函数为`skimage.feature.greycomatrix`。以下是其代码实现：

import numpy as np
from skimage.util.dtype import dtype_range
from skimage.util import img_as_ubyte

def greycomatrix(image, distances, angles, levels=None, symmetric=False,
                 normed=False):
    if levels is None:
        if np.issubdtype(image.dtype, np.integer):
            # For integer images, the number of levels is determined by the
            # range of values.
            imin, imax = dtype_range[image.dtype.type]
            levels = np.arange(imin, imax + 1)
        else:
            # For floating point images, we use 256 levels by default.
            levels = np.linspace(0, 1, 256, endpoint=False)

    # Convert the image to uint8, so that we can use it as an index into the
    # levels array.
    image = img_as_ubyte(image)

    glcm = np.zeros((len(levels), len(levels), len(distances), len(angles)),
                    dtype=np.double)

    for d, distance in enumerate(distances):
        for a, angle in enumerate(angles):
            coords1, coords2 = _offsets_to_coords(image.shape, distance, angle)
            patch1 = image[coords1[:, 0], coords1[:, 1]]
            patch2 = image[coords2[:, 0], coords2[:, 1]]

            # Find the indices into the levels array for each pixel in the
            # patch. This allows us to construct a joint histogram of levels.
            index1 = np.searchsorted(levels, patch1)
            index2 = np.searchsorted(levels, patch2)

            # Clip the indices to ensure that they are within the range [0,
            # len(levels)-1].
            index1 = np.clip(index1, 0, len(levels) - 1)
            index2 = np.clip(index2, 0, len(levels) - 1)

            # Update the co-occurrence matrix by counting the number of times
            # each combination of levels occurs.
            for i, j in zip(index1, index2):
                glcm[i, j, d, a] += 1

            if symmetric:
                # Make the matrix symmetric by mirroring over the diagonal.
                glcm = glcm + np.transpose(glcm, (1, 0, 2, 3))

    if normed:
        glcm /= np.sum(glcm)

    return glcm

该函数接受以下参数：

- `image`：输入的灰度图像。
- `distances`：灰度共生矩阵中用于计算概率的像素对之间的距离列表。例如，`[1, 2, 3]`表示将在距离为1、2和3的像素对之间计算概率。
- `angles`：灰度共生矩阵中用于计算概率的像素对之间的角度列表。例如，`[0, np.pi/4, np.pi/2, 3*np.pi/4]`表示将在水平、45度、垂直和135度方向上计算概率。
- `levels`：灰度级别列表。默认情况下，对于整数图像，该列表将由图像数据类型的范围确定（例如在[0, 255]范围内），而对于浮点图像，该列表将由256个等间距的值组成。您也可以指定自定义级别列表。
- `symmetric`：如果为True，则将灰度共生矩阵设置为对称矩阵。这在某些应用中很有用，例如，当您需要比较两个灰度图像时。
- `normed`：如果为True，则将灰度共生矩阵归一化为概率矩阵。

函数返回一个4D数组，其中每个元素都是一个灰度共生矩阵，表示给定距离和角度的像素对之间的概率。数组的第一个和第二个维度表示像素级别，第三个维度表示距离，第四个维度表示角度。