Summing each 3x3 window of a M*N matrix, into a M/3*N/3 matrix with numpy

Question

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!

Answer 1

The simplest numpy only approach, which does much less work thatn convolution and will therefore be likely faster than fileter 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.

Answer 2

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.