用于GLCM计算的Python中的滑动窗口 [英] Sliding window in Python for GLCM calculation

问题描述

我正在尝试使用GLCM算法在卫星图像中进行纹理分析。 scikit-image文档对此非常有帮助，但对于GLCM计算，我们需要在图像上循环一个窗口大小。这在Python中太慢了。我在stackoverflow上发现了许多关于滑动窗口的帖子，但计算需要永远。我有一个如下所示的例子，它可以工作，但需要永远。我想这一定是一种天真的做法

  image = np.pad（image，int（win / 2）， mode ='reflect'）
 row，cols = image.shape 
 feature_map = np.zeros（（M，N））
 
 for x in xrange（0，row） ：
 for x in xrange（0，cols）：
 window = image [m：m + win，n：n + win] 
 glcm = greycomatrix（window，d，theta，levels ）
 contrast = greycoprops（glcm，'contrast'）
 feature_map [m，n] = contrast 
  

我遇到了 skimage.util.view_as_windows 方法，这对我来说可能是一个很好的解决方案。我的问题是，当我尝试计算GLCM时，我得到一个错误，上面写着：

ValueError：参数 image 必须是二维数组

这是因为GLCM图像的结果有4d dimension和scikit-image view_as_windows 方法只接受2d数组。这是我的尝试

  win_w = 40 
 win_h = 40 
 
 features = np。零（image.shape，dtype ='uint8'）
 target = features [win_h // 2：-win_h // 2 + 1，win_w // 2：-win_w // 2 + 1] 
 windowed = view_as_windows（image，（win_h，win_w））
 
 GLCM = greycomatrix（窗口，[1]，[0，np.pi / 4，np.pi / 2,3 * np.pi / 4]，symmetric = True，normed = True）
 haralick = greycoprops（GLCM，'ASM'）
  

有没有人知道如何使用 skimage.util.view_as_windows 方法计算GLCM？

解决方案

您尝试执行的功能提取是一项计算机密集型任务。我通过仅为整个图像计算一次共现地图来加速你的方法，而不是在滑动窗口的重叠位置上反复计算共现图。

共现图是与原始图像大小相同的图像堆栈，其中 - 对于每个像素 - 强度级别由编码共现的整数代替两个强度，即该像素处的 Ii 和偏移像素处的 Ij 。共现图具有与我们考虑的偏移一样多的层（即所有可能的距离 - 角度对）。通过保留共生图，您不需要从头开始在滑动窗口的每个位置计算GLCM，因为您可以重复使用先前计算的共现图来获得每个距离的邻接矩阵（GLCM）角对。这种方法为您提供了显着的速度增益。

我想出的解决方案依赖于以下功能：

< pre class =lang-python prettyprint-override> 从skimage import导入numpy为np
来自scipy import stats的
来自skimage.feature的
导入greycoprops

def offset（长度，角度）：
返回给定长度和角度的偏移量
dv = length * np.sign（-np。 sin（角度））。astype（np.int32）
dh = length * np.sign（np.cos（angle））。astype（np.int32）
返回dv，dh

def crop（img，center，win）：
返回以中心为中心的方形裁剪（side = 2 * win + 1）
row，col =中心
side = 2 * win + 1
first_row = row - win
first_col = col - win
last_row = first_row + side
last_col = first_col + side
return img [first_row：last_row，first_col：last_col]

def cooc_maps（img，center，win， d = [1]，theta = [0]，levels = 256）：

为一个方形的
作物中心的不同d和theta返回一组共现图at center（side = 2 * w + 1）

shape =（2 * win + 1,2 * win + 1，len（d），len（theta））
cooc = np.zeros（shape = shape，dtype = np.int32）
row，col = center
Ii = crop（img，（row，col），win）
for d_index ，枚举的长度（d）：a_index的
，枚举的角度（theta）：
dv，dh = offset（长度，角度）
Ij = crop（img，center =（row） + dv，col + dh），win = win）
cooc [：，：，d_index，a_index] = encode_cooccurrence（Ii，Ij，levels）
返回cooc

def encode_cooccurrence（x，y，levels = 256）：
返回对应强度x和y共同出现的代码
返回x *级别+ y

def decode_cooccurrence（code，levels = 256）：
返回强度x，y对应代码
return co de // levels，np.mod（代码，级别）

def compute_glcms（cooccurrence_maps，levels = 256）：
计算共生图的同现频率
Nr，Na = cooccurrence_maps.shape [2：]
glcms = np.zeros（shape =（levels，levels，Nr，Na），dtype = np.float64）
表示r在范围内（Nr）：
的范围内（Na）：
table = stats.itemfreq（cooccurrence_maps [：，：，r，a]）
codes = table [：，0]
freqs = table [：，1] / float（table [：，1] .sum（））
i，j = decode_cooccurrence（codes，levels = levels）
glcms [i，j， r，a] = freqs
返回glcms

def compute_props（glcms，props =（'contrast'，））：
返回对应于集合的特征向量GLCM
Nr，Na = glcms.shape [2：]
features = np.zeros（shape =（Nr，Na，len（props）））
for index，枚举中的prop_name（props）：
features [：，：，index] = greycoprops（glcms，prop_name）
返回features.ravel（）

def haralick_features（img，win，d，theta，levels，props）：
返回Haralick功能的地图（一个功能）每像素矢量）
rows，cols = img.shape
margin = win + max（d）
arr = np.pad（img，margin，mode ='reflect'）
n_features = len（d）* len（theta）* len（道具）
feature_map = np.zeros（shape =（rows，cols，n_features），dtype = np.float64）
for x in xrange（rows）：
for n in xrange（cols）：
coocs = cooc_maps（arr，（m + margin，n + margin），win，d，theta，levels）
glcms = compute_glcms（coocs，levels）
feature_map [m，n，：] = compute_props（glcms，props）
return feature_map


---

DEMO

以下结果对应于Landsat图像中的（250,200）像素裁剪。我考虑了两个距离，四个角度和两个GLCM属性。这导致每个像素的16维特征向量。请注意，滑动窗口是平方的，其边是 2 * win + 1 像素（在此测试中，值 win = 19 被使用了）。此示例运行大约需要6分钟，这比永远短得多; - ）

 在[331]中：img.shape 
输出[331] :( 250L，200L）
 
在[332]中：img.dtype 
输出[332]：dtype（' uint8'）
 
在[333]中：d =（1,2）
 
在[334]中：theta =（0，np.pi / 4，np.pi / 2,3 * np.pi / 4）
 
在[335]中：props =（'contrast'，'homogeneity'）
 
在[336]中：levels = 256 
 
在[337]中：win = 19 
 
在[338]中：％time feature_map = haralick_features（img，win，d，theta，levels，props）
挂历时间：5分53秒
 
在[339]：feature_map.shape 
 Out [339] :( 250L，200L，16L）
 
在[340] ]：feature_map [0,0，：] 
 Out [340]：
 array（[10.3314,0.3477,25.1499,0.2738,25.1499,0.2738，
 25.1499,0.2738,23.5043,0.2755， 43.5523,0.1882，
 43.5523,0.1882,43.5523,0.1882]）
 
在[341]中：io.imshow（img）
输出[341]：< matpl otlib.image.AxesImage at 0xce4d160> 
  

I am trying to do texture analysis in a satellite imagery using GLCM algorithm. The scikit-image documentation is very helpful on that but for GLCM calculation we need a window size looping over the image. This is too slow in Python. I found many posts on stackoverflow about sliding windows but the computation takes for ever. I have an example shown below, it works but takes forever. I guess this must be a a naive way of doing it

image = np.pad(image, int(win/2), mode='reflect') 
row, cols = image.shape
feature_map = np.zeros((M, N))

for m in xrange(0, row):
    for n in xrange(0, cols):
        window = image[m:m+win, n:n+win]
        glcm = greycomatrix(window, d, theta, levels)
        contrast = greycoprops(glcm, 'contrast')
        feature_map[m,n] = contrast 

I came across with skimage.util.view_as_windows method which might be good solution for me. My problem is that, when I try to calculate the GLCM I get an error which says:

ValueError: The parameter image must be a 2-dimensional array

This is because the result of the GLCM image has 4d dimensions and scikit-image view_as_windows method accepts only 2d arrays. Here is my attempt

win_w=40
win_h=40

features = np.zeros(image.shape, dtype='uint8')
target = features[win_h//2:-win_h//2+1, win_w//2:-win_w//2+1]
windowed = view_as_windows(image, (win_h, win_w))

GLCM = greycomatrix(windowed, [1], [0, np.pi/4, np.pi/2, 3*np.pi/4], symmetric=True, normed=True)
haralick = greycoprops(GLCM, 'ASM')

Does anyone have an idea on how I can calculate the GLCM using skimage.util.view_as_windows method?

解决方案

The feature extraction you are trying to perform is a computer-intensive task. I have speeded up your method by computing the co-occurrence map only once for the whole image, rather than computing the co-occurrence map over and over on overlapping positions of the sliding window.

The co-occurrence map is a stack of images of the same size as the original image, in which - for each pixel - intensity levels are replaced by integer numbers that encode the co-occurrence of two intensities, namely Ii at that pixel and Ij at an offset pixel. The co-occurrence map has as many layers as we considered offsets (i.e. all the possible distance-angle pairs). By retaining the co-occurrence map you don't need to compute the GLCM at each position of the sliding window from the scratch, as you can reuse the previously computed co-occurrence maps to obtain the adjacency matrices (the GLCMs) for each distance-angle pair. This approach provides you with a significant speed gain.

The solution I came up with relies on the functions below:

import numpy as np
from skimage import io
from scipy import stats
from skimage.feature import greycoprops

def offset(length, angle):
    """Return the offset in pixels for a given length and angle"""
    dv = length * np.sign(-np.sin(angle)).astype(np.int32)
    dh = length * np.sign(np.cos(angle)).astype(np.int32)
    return dv, dh

def crop(img, center, win):
    """Return a square crop of img centered at center (side = 2*win + 1)"""
    row, col = center
    side = 2*win + 1
    first_row = row - win
    first_col = col - win
    last_row = first_row + side    
    last_col = first_col + side
    return img[first_row: last_row, first_col: last_col]

def cooc_maps(img, center, win, d=[1], theta=[0], levels=256):
    """
    Return a set of co-occurrence maps for different d and theta in a square 
    crop centered at center (side = 2*w + 1)
    """
    shape = (2*win + 1, 2*win + 1, len(d), len(theta))
    cooc = np.zeros(shape=shape, dtype=np.int32)
    row, col = center
    Ii = crop(img, (row, col), win)
    for d_index, length in enumerate(d):
        for a_index, angle in enumerate(theta):
            dv, dh = offset(length, angle)
            Ij = crop(img, center=(row + dv, col + dh), win=win)
            cooc[:, :, d_index, a_index] = encode_cooccurrence(Ii, Ij, levels)
    return cooc

def encode_cooccurrence(x, y, levels=256):
    """Return the code corresponding to co-occurrence of intensities x and y"""
    return x*levels + y

def decode_cooccurrence(code, levels=256):
    """Return the intensities x, y corresponding to code"""
    return code//levels, np.mod(code, levels)    

def compute_glcms(cooccurrence_maps, levels=256):
    """Compute the cooccurrence frequencies of the cooccurrence maps"""
    Nr, Na = cooccurrence_maps.shape[2:]
    glcms = np.zeros(shape=(levels, levels, Nr, Na), dtype=np.float64)
    for r in range(Nr):
        for a in range(Na):
            table = stats.itemfreq(cooccurrence_maps[:, :, r, a])
            codes = table[:, 0]
            freqs = table[:, 1]/float(table[:, 1].sum())
            i, j = decode_cooccurrence(codes, levels=levels)
            glcms[i, j, r, a] = freqs
    return glcms

def compute_props(glcms, props=('contrast',)):
    """Return a feature vector corresponding to a set of GLCM"""
    Nr, Na = glcms.shape[2:]
    features = np.zeros(shape=(Nr, Na, len(props)))
    for index, prop_name in enumerate(props):
        features[:, :, index] = greycoprops(glcms, prop_name)
    return features.ravel()

def haralick_features(img, win, d, theta, levels, props):
    """Return a map of Haralick features (one feature vector per pixel)"""
    rows, cols = img.shape
    margin = win + max(d)
    arr = np.pad(img, margin, mode='reflect')
    n_features = len(d) * len(theta) * len(props)
    feature_map = np.zeros(shape=(rows, cols, n_features), dtype=np.float64)
    for m in xrange(rows):
        for n in xrange(cols):
            coocs = cooc_maps(arr, (m + margin, n + margin), win, d, theta, levels)
            glcms = compute_glcms(coocs, levels)
            feature_map[m, n, :] = compute_props(glcms, props)
    return feature_map

---

DEMO

The following results correspond to a (250, 200) pixels crop from a Landsat image. I have considered two distances, four angles, and two GLCM properties. This results in a 16-dimensional feature vector for each pixel. Notice that the sliding window is squared and its side is 2*win + 1 pixels (in this test a value of win = 19 was used). This sample run took around 6 minutes, which is fairly shorter than "forever" ;-)

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

In [332]: img.dtype
Out[332]: dtype('uint8')

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

In [334]: theta = (0, np.pi/4, np.pi/2, 3*np.pi/4)

In [335]: props = ('contrast', 'homogeneity')

In [336]: levels = 256

In [337]: win = 19

In [338]: %time feature_map = haralick_features(img, win, d, theta, levels, props)
Wall time: 5min 53s    

In [339]: feature_map.shape
Out[339]: (250L, 200L, 16L)

In [340]: feature_map[0, 0, :]    
Out[340]:  
array([ 10.3314,   0.3477,  25.1499,   0.2738,  25.1499,   0.2738,
        25.1499,   0.2738,  23.5043,   0.2755,  43.5523,   0.1882,
        43.5523,   0.1882,  43.5523,   0.1882])

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

这篇关于用于GLCM计算的Python中的滑动窗口的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！