How to resize N-d numpy image?

I don't just want to subsample it, but to interpolate / average the pixels.

For example, if I start with

array([[[3, 1, 3, 1],
        [3, 1, 3, 1],
        [3, 1, 3, 1],
        [3, 1, 3, 1]],

       [[3, 1, 3, 1],
        [3, 1, 3, 1],
        [3, 1, 3, 1],
        [3, 1, 3, 1]]], dtype=uint8)

and shrink it by a factor of 2 in all dimensions, I want the output to be

array([[[2, 2],
        [2, 2]]], dtype=uint8)

Attempted solutions:

A. SciPy ndimage:

>>> scipy.ndimage.interpolation.zoom(x, .5, mode='nearest')

array([[[3, 1],
        [3, 1]]], dtype=uint8)

(The optional order parameter makes no difference)

B. Looping over 2**3 possible offsets: ugly, slow, works only for integer zoom factors, and needs extra steps to avoid overflows.

C. OpenCV and PIL work only with 2D images.

Solution

Reshape to split each axes into one more axis each of length 2, giving us a 6D array and then get the mean along those latter ones (axes : 1,3,5) -

m,n,r = a.shape
out = a.reshape(m//2,2,n//2,2,r//2,2).mean((1,3,5))

Extending to n-dim array, it would be -

def shrink(a, S=2): # S : shrink factor
    new_shp = np.vstack((np.array(a.shape)//S,[S]*a.ndim)).ravel('F')
    return a.reshape(new_shp).mean(tuple(1+2*np.arange(a.ndim)))

Sample run -

In [407]: a
Out[407]: 
array([[[1, 5, 8, 2],
        [5, 6, 4, 0],
        [8, 5, 5, 5],
        [1, 0, 0, 0]],

       [[0, 0, 7, 6],
        [3, 5, 4, 3],
        [4, 5, 1, 3],
        [6, 7, 4, 0]]])

In [408]: a[:2,:2,:2].mean()
Out[408]: 3.125

In [409]: a[:2,:2,2:4].mean()
Out[409]: 4.25

In [410]: a[:2,2:4,:2].mean()
Out[410]: 4.5

In [411]: a[:2,2:4,2:4].mean()
Out[411]: 2.25

In [412]: shrink(a, S=2)
Out[412]: 
array([[[ 3.125,  4.25 ],
        [ 4.5  ,  2.25 ]]])