5

I'm trying to implement a function that would sum (or eventually average) Each 3x3 window of a given matrix, and create a matrix 9x smaller with the result of each window.

I can't figure out an efficient and concise way to do this with numpy.

Any ideas ?

Thank you!

mtourne
  • 440
  • 1
  • 6
  • 13

2 Answers2

8

The simplest numpy only approach, which does much less work than convolution and will therefore be likely faster than filter based methods, is to resize your original array to one with extra dimensions, then reduce it back to normal by summing over the new dimensions:

>>> arr = np.arange(108).reshape(9, 12)
>>> rows, cols = arr.shape
>>> arr.reshape(rows//3, 3, cols//3, 3).sum(axis=(1, 3))
array([[117, 144, 171, 198],
       [441, 468, 495, 522],
       [765, 792, 819, 846]])

If you wanted the mean, you would simply divide the resulting array by the number of elements:

>>> arr.reshape(rows//3, 3, cols//3, 3).sum(axis=(1, 3)) / 9
array([[ 13.,  16.,  19.,  22.],
       [ 49.,  52.,  55.,  58.],
       [ 85.,  88.,  91.,  94.]])

This method only works if your array has a shape which is itself a multiple of 3.

ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • 1
    I spent a long time wondering why this wasn't working, turned out sum() was updated with numpy 1.7 and axis=(1,3) wasn't working in my version (1.6) – mtourne Nov 18 '13 at 05:47
  • Yes, you are right. You could get the same result in earlier versions doing `.sum(axis=3).sum(axis=1)` instead of `.sum(axis=(1, 3))`. Note that the order is important, i.e. `.sum(axis=1).sum(axis=3)` will raise an error because after the first summation the array is only 3D, and hence `axis=3` is out of bounds. – Jaime Nov 18 '13 at 06:55
4

To achieve exactly what you are asking for I would apply a [3x3] box-filter on the image and than I would resize the matrix using nearest neighbor interpolation.

# Pseudo code
kernel = np.array([[1/9, 1/9, 1/9],
                   [1/9, 1/9, 1/9],
                   [1/9, 1/9, 1/9]])
avg_data= ndimage.convolve(data, kernel)

smaller_data = scipy.misc.imresize(avg_data, org_size/3, interp='nearest', mode=None)

In case you want something more efficient - as pointed by @Jaime - you can do something like this How can I efficiently process a numpy array in blocks similar to Matlab's blkproc (blockproc) function :

from numpy.lib.stride_tricks import as_strided as ast

def block_view(A, block= (3, 3)):
    """Provide a 2D block view to 2D array. No error checking made.
    Therefore meaningful (as implemented) only for blocks strictly
    compatible with the shape of A."""
    # simple shape and strides computations may seem at first strange
    # unless one is able to recognize the 'tuple additions' involved ;-)
    shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
    strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
    return ast(A, shape= shape, strides= strides)

if __name__ == '__main__':
    B = block_view(A).sum(axis=(2,3))

When you try to understand what's going on, remember that strides represent the amount of bytes we need to offset in memory to traverse to the next cell in each dimension. So if you are dealing with int32 data type it would be a multiply of 4.

Community
  • 1
  • 1
zenpoy
  • 19,490
  • 9
  • 60
  • 87
  • You could do the last step by slicing, i.e. `avg_data[1::3, 1::3]` gives you the right result, almost certainly much faster. And it kind of points out that you are having to throw away 89% of the values you calculate, so using a filter here is probably far from optimal. – Jaime Nov 17 '13 at 13:40
  • It works, but I think the correct sum is block_view(A).sum(axis=(2,3)). Btw, is there a good explanations of "as_strided" for humans somewhere ? – mtourne Nov 18 '13 at 05:48