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=40features = 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 greycopropsdef 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, dhdef crop(img, center, win):"""Return a square crop of img centered at center (side = 2*win + 1)"""row, col = centerside = 2*win + 1first_row = row - winfirst_col = col - winlast_row = first_row + side    last_col = first_col + sidereturn 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 = centerIi = 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 coocdef encode_cooccurrence(x, y, levels=256):"""Return the code corresponding to co-occurrence of intensities x and y"""return x*levels + ydef 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] = freqsreturn glcmsdef 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.shapemargin = 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 = 256In [337]: win = 19In [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>