2

我正在使用 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 尺寸）？或者我别无选择，只能以正常方式循环图像（不设置步长）？

4

1 回答 1

1

You cannot do that using 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 this thread for details).

But if you wish to stick to mahotas, you should try using skimage.util.view_as_windows rather than a sliding window as it could speed up the scanning across the image. Please, be sure to read the caveat about memory usage at the end of the documentation. If using view_as_windows is an affordable approach for you, the following code gets the job done:

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,))

DEMO

For the sample run below I have set win to 19 which corresponds to a window of shape (39, 39). I have considered two different distances. Notice that mht.haralick yields 13 GLCM features for four directions. Altogether this results in a 104-dimensional feature vector for each pixel. When applied to a (250, 200) pixels crop from a Landsat image, feature extraction finishes in around 7 minutes.

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>

Landsat image

于 2017-03-15T13:04:32.307 回答