scikit-image Python中greycomatrix的levels参数

Question

我正在把我的Matlab图像处理算法转到Python，使用scikit-image工具。我在计算灰度共生矩阵（GLCM），用的是greycomatrix这个函数。现在我遇到一个问题，就是当参数levels小于图像的最大强度值（image.max()）时，就会出错。例如：

import numpy as np
from skimage.feature import greycomatrix
image = np.array([[0, 0, 1, 1],[0, 0, 1, 1],[0, 2, 2, 2],[2, 2, 3, 3]], dtype=np.uint8)
result = greycomatrix(image, distances = [1], angles = [0], levels = 4, symmetric=True)

输出结果是：

glcm = result[:,:,0,0]

array([[4, 2, 1, 0],
   [2, 4, 0, 0],
   [1, 0, 6, 1],
   [0, 0, 1, 2]], dtype=uint32)

这个结果是对的，是一个4x4的矩阵。但是如果我把levels=3，我就无法计算GLCM，错误信息是：

result = greycomatrix(image, distances = [1], angles = [0], levels = 3, symmetric=True)

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/usr/lib/python3.4/site-packages/skimage/feature/texture.py", line 97, in greycomatrix
assert image.max() < levels
AssertionError

当然，我收到了错误提示，但我应该能计算出一个levels小于image.max()的GLCM（3x3的矩阵）。比如，对于：

result = greycomatrix(image, distances = [1], angles = [0], levels = 3, symmetric=True)

我应该得到以下的GLCM（我在Matlab里可以做到）：

4     3     0
3    10     1
0     1     2

当我处理很大的图像时，我会减少GLCM的级别，以此来缩短计算时间。请问greycomatrix有没有问题，还是我想错了？谢谢！

Answer 1

也许我来得有点晚，但希望我的回答能对将来其他人有所帮助...

根据Matlab的文档，'NumLevels'可以小于max(image(:))，这是因为灰度值被截断了：

小于或等于低值的灰度值会被调整为1，而大于或等于高值的灰度值会被调整为'NumLevels'。

你也可以在Python中使用类似的方法，通过NumPy的clip来实现：

In [11]: import numpy as np

In [12]: from skimage.feature import greycomatrix

In [13]: image = np.array([[0, 0, 1, 1], 
    ...:                   [0, 0, 1, 1], 
    ...:                   [0, 2, 2, 2], 
    ...:                   [2, 2, 3, 3]], dtype=np.uint8)
    ...: 

In [14]: clipped = np.clip(image, a_min=0, a_max=2)

In [15]: clipped
Out[15]: 
array([[0, 0, 1, 1],
       [0, 0, 1, 1],
       [0, 2, 2, 2],
       [2, 2, 2, 2]], dtype=uint8)

In [16]: result = greycomatrix(clipped, 
    ...:                       distances=[1], 
    ...:                       angles=[0], 
    ...:                       levels=3, 
    ...:                       symmetric=True)
    ...: 

In [17]: result[:, :, 0, 0]
Out[17]: 
array([[ 4,  2,  1],
       [ 2,  4,  0],
       [ 1,  0, 10]], dtype=uint32)

如果你想减少灰度共生矩阵的大小，我建议重新量化图像，而不是简单地截断灰度值。这个方法可以通过NumPy的digitize轻松实现：

In [59]: nlevels = 256

In [60]: nbins = 64

In [61]: arr = np.arange(nlevels).reshape((16, 16)).astype(np.uint8).T

In [62]: np.set_printoptions(threshold=300, linewidth=100)

In [63]: arr
Out[63]: 
array([[  0,  16,  32,  48,  64,  80,  96, 112, 128, 144, 160, 176, 192, 208, 224, 240],
       [  1,  17,  33,  49,  65,  81,  97, 113, 129, 145, 161, 177, 193, 209, 225, 241],
       [  2,  18,  34,  50,  66,  82,  98, 114, 130, 146, 162, 178, 194, 210, 226, 242],
       [  3,  19,  35,  51,  67,  83,  99, 115, 131, 147, 163, 179, 195, 211, 227, 243],
       [  4,  20,  36,  52,  68,  84, 100, 116, 132, 148, 164, 180, 196, 212, 228, 244],
       [  5,  21,  37,  53,  69,  85, 101, 117, 133, 149, 165, 181, 197, 213, 229, 245],
       [  6,  22,  38,  54,  70,  86, 102, 118, 134, 150, 166, 182, 198, 214, 230, 246],
       [  7,  23,  39,  55,  71,  87, 103, 119, 135, 151, 167, 183, 199, 215, 231, 247],
       [  8,  24,  40,  56,  72,  88, 104, 120, 136, 152, 168, 184, 200, 216, 232, 248],
       [  9,  25,  41,  57,  73,  89, 105, 121, 137, 153, 169, 185, 201, 217, 233, 249],
       [ 10,  26,  42,  58,  74,  90, 106, 122, 138, 154, 170, 186, 202, 218, 234, 250],
       [ 11,  27,  43,  59,  75,  91, 107, 123, 139, 155, 171, 187, 203, 219, 235, 251],
       [ 12,  28,  44,  60,  76,  92, 108, 124, 140, 156, 172, 188, 204, 220, 236, 252],
       [ 13,  29,  45,  61,  77,  93, 109, 125, 141, 157, 173, 189, 205, 221, 237, 253],
       [ 14,  30,  46,  62,  78,  94, 110, 126, 142, 158, 174, 190, 206, 222, 238, 254],
       [ 15,  31,  47,  63,  79,  95, 111, 127, 143, 159, 175, 191, 207, 223, 239, 255]], dtype=uint8)

In [64]: binned = np.uint8(np.digitize(arr, np.arange(0, nlevels, nbins))) - 1

In [65]: binned
Out[65]: 
array([[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]], dtype=uint8)

In [66]: glcm = greycomatrix(binned, 
    ...:                     distances=[1], 
    ...:                     angles=[0], 
    ...:                     levels=binned.max()+1, 
    ...:                     symmetric=True)
    ...: 

In [67]: glcm[:, :, 0, 0]    
Out[67]: 
array([[96, 16,  0,  0],
       [16, 96, 16,  0],
       [ 0, 16, 96, 16],
       [ 0,  0, 16, 96]], dtype=uint32)