我正在使用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 。

因此,我的问题是:通过设置步长(以避免在重叠区域中进行计算以及代码变得更快),我可以以某种方式得出尺寸为250 x 200的整个图像的GLCM结果,而不是它的子集( 6 x 8尺寸)?或者我别无选择,只能以正常方式遍历图像(无需设置步长)?

最佳答案

您不能使用mahotas作为计算共现映射的函数来执行此操作,该库在此库中不可用。从GLCM中提取纹理特征的另一种方法是利用skimage.feature.graycoprops(有关详细信息,请参见this thread)。

但是,如果您希望坚持使用mahotas,则应尝试使用skimage.util.view_as_windows而不是滑动窗口,因为它可以加快整个图像的扫描速度。请确保阅读文档末尾有关内存使用的警告。如果使用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>


python - Mahotas库,用于GLCM计算和窗口大小-LMLPHP

关于python - Mahotas库,用于GLCM计算和窗口大小,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/42624729/

10-09 19:18