Iteration multi-dimensional array in Python

Question

I want to compute some operation (GLCM) over an image using a sliding window. The code to achieve it is:

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

from numpy.lib.stride_tricks import as_strided

image = np.arange(36).reshape((6,6))

window = 5

result = np.zeros(image.shape)
for i in xrange(window/2,image.shape[0]-window/2):
    for j in xrange(window/2,image.shape[1]-window/2):
        sample = image[i-(window/2):i+(window/2)+1, j - (window/2):j+(window/2)+1]
        glcm = greycomatrix(sample, [1], [0], 256, symmetric=False, normed=True)
        result[i,j] = greycoprops(glcm, 'contrast')[0, 0]

It works but two for loops is very expensive. I want to increase the speed and so, looking around on the web, I tried to use the as_stride trick:

from numpy.lib.stride_tricks import as_strided

image = np.arange(36).reshape((6,6))

window = 5

y = as_strided(image,shape=(image.shape[0] - window + 1,\
                            image.shape[1] - window + 1,) +\
               (window,window),  strides=image.strides * 2)

To compute, for e.g., the GLCM for the first window:

glcm = greycoprops(greycomatrix(y[0,0], [1], [0], 256, symmetric=False, normed=True))[0][0]

I tried to apply for all the sliding window as:

glcm[:,:] = greycoprops(greycomatrix(y[:,:], [1], [0], 256, symmetric=False, normed=True))[0][0]

But in this case, the y[:,:] doesn't have ndim==2 as y[0,0] but ndim==4, and so . I can't find a way to iterate in a smart way over all the subset keeping the ndim == 2 (required by the greycomatrix function).

Edit

I tried to use ravel and work over 1D vector and so, just 1 for-loop. This is the code:

a = y.ravel()
print a.shape
glcm=np.zeros(a.shape[0]/(window*window))
for i in np.arange(a.shape[0]/(window*window)):
    glcm[i] = greycoprops(greycomatrix(a[i*25:i*25+25].reshape(5,5), [1], [0], 256, symmetric=False, normed=True))[0][0]

result= glcm.reshape(y.shape[0],y.shape[1])

The processing time increase...

score 1 · Answer 1 · answered Jan 16 '15 at 13:36

Since you forgot to actually ask a question, I'm going to assume it was

How do I make this run fast?

Well, in this case the harsh truth is, as good as python is going to get, doing a python for loop with a lot of slicing is always going to be relatively expensive compared to individual pixel operations.

So if speed is your concern here, you should implement some function in a language that allows you to get python bindings (eg. C with cpython), and use that function.

score 0 · Answer 2 · answered Jan 16 '15 at 19:10

With the as_strided y you still need to access the sub-arrays individually, for example:

for i in range(y.shape[0]):
    for j in range(y.shape[1]):
        sample = y[i,j,...]
        print sample

or even

for row in y:
    for sample in row:
        print sample

except that you need to collect the results.

In iterations like this, as_strided has a benefit only if it makes access to the subarrays more efficient.

But its real benefit comes if you can rewrite the GCLM calculation to work with a 4d arrays. Some numpy operations are designed to operate on 1 or 2 of the axes, while the others just 'go along for the ride'. For example, if your calculation simply consisted of taking the mean of the image. With y that would simply be:

np.mean(y, axes=(-2,-1))

Iteration multi-dimensional array in Python

2 Answers2