用于 GLCM 计算和窗口大小的 Mahotas 库答案

【问题标题】：Mahotas library for GLCM calulation and window size用于 GLCM 计算和窗口大小的 Mahotas 库
【发布时间】：2017-07-26 06:28:46
【问题描述】：

我正在使用 mahotas 库对卫星图像（250 x 200 像素）进行纹理分析 (GLCM)。 GLCM 计算在窗口大小内进行。因此，对于滑动窗口的两个相邻位置，我们需要从头开始计算两个共现矩阵。我读过我也可以设置步长，以避免在重叠区域上计算 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)

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

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

【问题讨论】：

您的代码示例很可能不完整：您总是使用与haralick_feature 相同的参数
嗨 Luispedro，你可能是对的，我可能错过了在这里粘贴一些代码。我创建了这个脚本的几个不同版本（所有脚本都可以正常工作），试图弄清楚如何获得整个图像的 GLCM 结果（设置步长时）而不是它的一个小主题。如果我不设置步长，这个过程非常耗时。

标签： python computer-vision scikit-image mahotas glcm

【解决方案1】：

您不能使用mahotas 来执行此操作，因为此库中没有计算共现图 的函数。从 GLCM 中提取纹理特征的另一种方法是使用skimage.feature.graycoprops（详情请查看this thread）。

但是，如果您想坚持使用 mahotas，则应尝试使用 skimage.util.view_as_windows 而不是滑动窗口，因为它可以加快图像的扫描速度。请务必阅读文档末尾关于 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) 的窗口。我考虑了两种不同的距离。请注意，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>

【讨论】：

感谢您的详细解答。我很感激。