用于 GLCM 计算和 window 大小的 Mahotas 库

Question

我正在使用 mahotas 库对卫星图像（250 x 200 像素）进行纹理分析 (GLCM)。 GLCM 计算在 window 尺寸内进行。因此，对于滑动window的两个相邻位置，我们需要从头开始计算两个共生矩阵。我读过我也可以设置步长以避免在重叠区域计算 GLCM。我提供了以下代码。

#Compute haralick features
def haralick_feature(image):
    haralick = mahotas.features.haralick(image, True)
    return haralick


img = 'SAR_image.tif'
win_s=32 #window size
step=32 #step size

rows = img.shape[0]
cols = img.shape[1]
array = np.zeros((rows,cols), dtype= object)
harList = []

for i in range(0, rows-win_s-1, step):
        print 'Row number: ', r
    for j in range(0, cols-win_s-1, step):
        harList.append(haralick_feature(image))

harImages = np.array(harList)     
harImages_mean = harImages.mean(axis=1)

对于上面的代码，我已将 window 和步长设置为 32。当代码完成时，我得到一个尺寸为 6 x 8（而不是 250 x 200）的图像，这是有道理的因为步长已设置为 32.

所以，我的问题是：通过设置步长（以避免在重叠区域进行计算以及代码变得更快），我可以以某种方式导出尺寸为 250 x 200 的整个图像的 GLCM 结果，而不是让它的一个子集（6 x 8 维）？或者我别无选择，只能以正常方式循环播放图像（不设置步长）？

Answer 1

您无法使用 mahotas as the function which computes the co-occurrence map is not available in this library. An alternative approach to the extraction of texture features from the GLCM would consist in utilizing skimage.feature.graycoprops (check out 了解详细信息）。

但是如果你想坚持使用 mahotas，你应该尝试使用 skimage.util.view_as_windows 而不是滑动 window，因为它可以加快图像的扫描速度。请务必阅读文档末尾关于 memory 用法的警告。如果使用 view_as_windows 对您来说是一种负担得起的方法，则以下代码可以完成工作：

import numpy as np
from skimage import io, util
import mahotas.features.texture as mht

def haralick_features(img, win, d):
    win_sz = 2*win + 1
    window_shape = (win_sz, win_sz)
    arr = np.pad(img, win, mode='reflect')
    windows = util.view_as_windows(arr, window_shape)
    Nd = len(d)
    feats = np.zeros(shape=windows.shape[:2] + (Nd, 4, 13), dtype=np.float64)
    for m in xrange(windows.shape[0]):
        for n in xrange(windows.shape[1]):
            for i, di in enumerate(d):
                w = windows[m, n, :, :]
                feats[m, n, i, :, :] = mht.haralick(w, distance=di)
    return feats.reshape(feats.shape[:2] + (-1,))

演示版

对于下面的示例运行，我已将 win 设置为 19，对应于形状 (39, 39) 的 window。我考虑了两种不同的距离。请注意，mht.haralick 为四个方向产生了 13 个 GLCM 特征。总而言之，这会为每个像素产生一个 104 维的特征向量。当应用于来自 Landsat 图像的 (250, 200) 像素裁剪时，特征提取在大约 7 分钟内完成。

In [171]: img = io.imread('landsat_crop.tif')

In [172]: img.shape
Out[172]: (250L, 200L)

In [173]: win = 19

In [174]: d = (1, 2)

In [175]: %time feature_map = haralick_features(img, win, d)
Wall time: 7min 4s

In [176]: feature_map.shape
Out[176]: (250L, 200L, 104L)

In [177]: feature_map[0, 0, :]
Out[177]: 
array([  8.19278030e-03,   1.30863698e+01,   7.64234582e-01, ...,
         3.59561817e+00,  -1.35383606e-01,   8.32570045e-01])

In [178]: io.imshow(img)
Out[178]: <matplotlib.image.AxesImage at 0xc5a9b38>

用于 GLCM 计算和 window 大小的 Mahotas 库

Mahotas library for GLCM calulation and window size

python

computer-vision

scikit-image

mahotas

glcm