Subsampling 3D array using the neighbourhood sum

Question

The title is probably confusing. I have a reasonably large 3D numpy array. I'd like to cut it's size by 2^3 by binning blocks of size (2,2,2). Each element in the new 3D array should then contain the sum of the elements in it's respective block in the original array.

As an example, consider a 4x4x4 array:

input = [[[1, 1, 2, 2],
          [1, 1, 2, 2],
          [3, 3, 4, 4],
          [3, 3, 4, 4]],
         [[1, 1, 2, 2],
          [1, 1, 2, 2],
          [3, 3, 4, 4],
          [3, 3, 4, 4]],
              ...    ]]]

(I'm only representing half of it to save space). Notice that all the elements with the same value constitute a (2x2x2) block. The output should be a 2x2x2 array such that each element is the sum of a block:

output = [[[8, 16],
          [24, 32]],
             ... ]]]

So 8 is the sum of all 1's, 16 is the sum of the 2's, and so on.

Answer 1

There's a builtin to do those block-wise reductions - skimage.measure.block_reduce-

In [36]: a
Out[36]: 
array([[[1, 1, 2, 2],
        [1, 1, 2, 2],
        [3, 3, 4, 4],
        [3, 3, 4, 4]],

       [[1, 1, 2, 2],
        [1, 1, 2, 2],
        [3, 3, 4, 4],
        [3, 3, 4, 4]]])

In [37]: from skimage.measure import block_reduce

In [39]: block_reduce(a, block_size=(2,2,2), func=np.sum)
Out[39]: 
array([[[ 8, 16],
        [24, 32]]])

Use other reduction ufuncs, say max-reduction -

In [40]: block_reduce(a, block_size=(2,2,2), func=np.max)
Out[40]: 
array([[[1, 2],
        [3, 4]]])

Implementing such a function isn't that difficult with NumPy tools and could be done like so -

def block_reduce_numpy(a, block_size, func):
    shp = a.shape
    new_shp = np.hstack([(i//j,j) for (i,j) in zip(shp,block_size)])
    select_axes = tuple(np.arange(a.ndim)*2+1)
    return func(a.reshape(new_shp),axis=select_axes)