14

I am apply an operation on a moving window of constant size across a 2D array. Is there an efficient vectorize-like operation I can implement to do this without looping in Python? My current structure looks something like this

 for i in range(1,xmax-1):
     for j in range(1,ymax-1):
        out[i][j] = f(in[i][j],in[i+1][j],in[i-1][j],in[i][j+1],in[i][j-1],...)

The comments that eat left in this question allude to the possibility of vectorizing this operation this, but without further details vectorized indexing/slicing in numpy/scipy?

Community
  • 1
  • 1
Rich
  • 12,068
  • 9
  • 62
  • 94
  • For applying a generic NumPy ufunc, you can put every block into a column, similar to MATLAB has with [`im2col`](http://www.mathworks.com/help/images/ref/im2col.html). A vectorized implementation of the same in NumPy/Python is listed in [`Implement Matlab's im2col 'sliding' in Python`](http://stackoverflow.com/questions/30109068/implement-matlabs-im2col-sliding-in-python). Also, you can look [`here`](http://stackoverflow.com/a/34327984/3293881) to see some examples. – Divakar Jan 02 '16 at 08:39

2 Answers2

15

You can use the rolling window technique as explained here, here and here, but for 2D array.

The source code for 2D rolling window in NumPy:

# Rolling window for 2D arrays in NumPy
import numpy as np

def rolling_window(a, shape):  # rolling window for 2D array
    s = (a.shape[0] - shape[0] + 1,) + (a.shape[1] - shape[1] + 1,) + shape
    strides = a.strides + a.strides
    return np.lib.stride_tricks.as_strided(a, shape=s, strides=strides)

a = np.array([[0,  1,  2,  3,  4,  5],
              [6,  7,  8,  9, 10,  11],
              [12, 13, 14, 15, 7,   8],
              [18, 19, 20, 21, 13, 14],
              [24, 25, 26, 27, 19, 20],
              [30, 31, 32, 33, 34, 35]], dtype=np.int)
b = np.arange(36, dtype=np.float).reshape(6,6)
present = np.array([[7,8],[13,14],[19,20]], dtype=np.int)
absent  = np.array([[7,8],[42,14],[19,20]], dtype=np.int)

found = np.all(np.all(rolling_window(a, present.shape) == present, axis=2), axis=2)
print(np.transpose(found.nonzero()))
found = np.all(np.all(rolling_window(b, present.shape) == present, axis=2), axis=2)
print(np.transpose(found.nonzero()))
found = np.all(np.all(rolling_window(a, absent.shape) == absent, axis=2), axis=2)
print(np.transpose(found.nonzero()))

Array present is occurred in array a two times on [1,1] and [2,4].

More examples in my CoLab notebook "Rolling window on NumPy arrays without for loops".

FooBar167
  • 2,721
  • 1
  • 26
  • 37
  • 1
    It looks like the `rolling_window` function is implemented at [numpy.lib.stride_tricks.sliding_window_view](https://numpy.org/devdocs/reference/generated/numpy.lib.stride_tricks.sliding_window_view.html) as of numpy 1.20.0 – sappjw May 04 '22 at 19:48
7

If you can express the function 

f(in[i][j],in[i+1][j],in[i-1][j],in[i][j+1],in[i][j-1],…)

as a linear operator, you could use scipy's signal.convolve2d function to do exactly that. For instance, say you have an 50x50 array, A, and you want to calculate a second array B where each of its element b[ij] is the average over a[i,j], a[(i-1),j], a[i,(j-1)], a[(i-1),(j-1)] from the array A. You could do that simply doing :

A = # your first array
B = numpy.ones((2,2))/4
C = scipy.signal.convolve2d(A,B, 'valid')

When the convolution is performed, the array B "slides" across A, multiplying the corresponding elements and summing up the result. Because of border effects, you must be careful when using the resulting array C. Here, C is of shape 49x49, because of the 'valid' argument in convolve2d, to discard the first row and column since they contain border effects. If you wanted to have a 50x50 array, without discarding, you would swap that argument for 'same'

EDIT: Perhaps if you could tell me more about that function you need, I could help you more specifically in turning it into an array that would be used to do the 2D convolution.

Hope that helps!

matehat
  • 5,214
  • 2
  • 29
  • 40