Image texture with skimage

I'm trying to get texture properties from a GLCM I created using greycomatrix from skimage.feature. My input data is an image with multiple bands and I want the texture properties for each pixel (resulting in an image with the dimensions cols x rows x (properties *bands)), as it can be achieved using ENVI. But I'm too new to this to come to grips with greycomatrix and greycoprops. This is what I tried:

import numpy as np
from skimage import io
from skimage.feature import greycomatrix, greycoprops

array = io.imread('MYFILE.tif')
array = array.astype(np.int64)
props = ['contrast', 'dissimilarity', 'homogeneity', 'energy', 'correlation', 'ASM']
textures = np.zeros((array.shape[0], array.shape[1], array.shape[2] * len(props)), np.float32)
angles = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4]
bands = array.shape[2]
for b in range(bands):
    glcm = greycomatrix(array[:, :, b], [1], angles, np.nanmax(array) + 1,
                        symmetric=True, normed=True)
    for p, prop in enumerate(props):
        textures[:, :, b] = greycoprops(glcm, prop)

Unfortunately, this gives me a 1 x 4 matrix per prop, which I guess is one value per angle FOR THE WHOLE IMAGE, but this is not what I want. I need it per pixel, like contrast for each single pixel, computed from its respective surroundings. What am I missing?

This snippet should get the job done:

import numpy as np
from skimage import io, util
from skimage.feature.texture import greycomatrix, greycoprops

img = io.imread('fourbandimg.tif')

rows, cols, bands = img.shape

radius = 5
side = 2*radius + 1

distances = [1]
angles = [0, np.pi/2]
props = ['contrast', 'dissimilarity', 'homogeneity']
dim = len(distances)*len(angles)*len(props)*bands

padded = np.pad(img, radius, mode='reflect')
windows = [util.view_as_windows(padded[:, :, band].copy(), (side, side))
           for band in range(bands)]
feats = np.zeros(shape=(rows, cols, dim))

for row in range(rows):
    for col in range(cols):
        pixel_feats = []
        for band in range(bands):
            glcm = greycomatrix(windows[band][row, col, :, :],
                                distances=distances,
                                angles=angles)
            pixel_feats.extend([greycoprops(glcm, prop).ravel()
                                for prop in props])
        feats[row, col, :] = np.concatenate(pixel_feats)

The sample image has 128 rows, 128 columns and 4 bands (click here to download). At each image pixel a square local neighbourhood of size 11 is used to compute the grayscale matrices corresponding to the pixel to the right and the pixel above for each band. Then, contrast, dissimilarity and homogeneity are computed for those matrices. Thus we have 4 bands, 1 distance, 2 angles and 3 properties. Hence for each pixel the feature vector has 4 × 1 × 2 × 3 = 24 components.

Notice that in order to preserve the number of rows and columns the image has been padded using the image itself mirrored along the edge of the array. It this approach does not fit your needs you could simply ignore the outer frame of the image.

As a final caveat, the code could take a while to run.

Demo

In [193]: img.shape
Out[193]: (128, 128, 4)

In [194]: feats.shape
Out[194]: (128, 128, 24)

In [195]: feats[64, 64, :]
Out[195]: 
array([  1.51690000e+04,   9.50100000e+03,   1.02300000e+03,
         8.53000000e+02,   1.25203577e+01,   9.38930575e+00,
         2.54300000e+03,   1.47800000e+03,   3.89000000e+02,
         3.10000000e+02,   2.95064854e+01,   3.38267222e+01,
         2.18970000e+04,   1.71690000e+04,   1.21900000e+03,
         1.06700000e+03,   1.09729371e+01,   1.11741654e+01,
         2.54300000e+03,   1.47800000e+03,   3.89000000e+02,
         3.10000000e+02,   2.95064854e+01,   3.38267222e+01])

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

Edit

You could cast your data to the type required by greycomatrix through NumPy's uint8 or scikit-images's img_as_ubyte.