Binning multidimensional array in numpy

user_5558281 0
Viewed 301

I have a 4d numpy array (these are stacks of imaging data) and would like to perform mean binning along all but one of the axes.

starting with say

x=np.random.random((3,100,100,100))

I want to apply binning to axes 1,2,3 with bin size 10 and average the values in each bin.

expected result would be an array of shape (3,10,10,10)

I have looked into np.reshape like so:

result=x.reshape(3,-1,10,100,100).mean(axis=1)
result=result.reshape(3,10,-1,10,100).mean(axis=2)

and so on, but this messes up the structure of the image arrays

is there a more straightforward way to do this?

numpy arrays python
4 Answers

How about this:

import numpy as np
import skimage.measure

a = np.arange(36).reshape(6, 6)

b = skimage.measure.block_reduce(a, (2,2), np.mean)

output:

a = 
[[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]
 [12 13 14 15 16 17]
 [18 19 20 21 22 23]
 [24 25 26 27 28 29]
 [30 31 32 33 34 35]]
b = 
[[ 3.5  5.5  7.5]
 [15.5 17.5 19.5]
 [27.5 29.5 31.5]]

But instead of my 2d example, you can do that for a block size of (1, 10, 10, 10) of your data.

user_153066060 
#block size 
bs = (10,10,10)
s = 1
shape = [3,
         x.shape[s+0]//bs[0], bs[0],
         x.shape[s+1]//bs[1], bs[1]
         x.shape[s+2]//bs[2], bs[2]]
result = x.reshape(*shape).mean(axis = (2,4,6))

Possibly a redundant answer at this point, but if you prefer not to use skimage then this should do the same thing.

user_118953990

In my experience cooltools (https://cooltools.readthedocs.io/en/latest/index.html) works terrific for rescaling (either up or down) images of any shape. Can be 1D, 2D, 3D etc.

from cooltools import numutils

a = np.random.random((3,100,100,100))

#down-sample
rescale1 = numutils.zoom_array(a, (3,10,10,10))
#up-sample
rescale2 = numutils.zoom_array(a, (3,200,200,200))

rescale1.shape, rescale2.shape

Out[1]: ((3, 10, 10, 10), (3, 200, 200, 200))
user_115456940

Try this:

x=np.random.random((3,100,100,100))

x_resized=np.zeros((3,10,10,10))

for i in range(len(x_resized[0])):
    for j in range(len(x_resized[0][0])):
        for k in range(len(x_resized[0][0][0])):
            
            x_resized[0,i,j,k]=np.average(x[0,i*10:i*10+10,j*10:j*10+10,k*10:k*10+10])
            x_resized[1,i,j,k]=np.average(x[1,i*10:i*10+10,j*10:j*10+10,k*10:k*10+10])
            x_resized[2,i,j,k]=np.average(x[2,i*10:i*10+10,j*10:j*10+10,k*10:k*10+10])

which performs the averaging blockwise.

user_101466260
Related