2

I have a 3D NumPy array a of shape (2, 9, 9) like this one:

a = np.array([
       [[4, 5, 1, 3, 8, 8, 0, 6, 6],
        [9, 2, 2, 1, 8, 2, 2, 4, 5],
        [2, 3, 2, 2, 5, 3, 1, 2, 4],
        [9, 6, 2, 9, 1, 0, 6, 2, 3],
        [4, 2, 7, 7, 9, 1, 3, 7, 2],
        [5, 8, 9, 4, 6, 3, 1, 6, 7],
        [3, 6, 4, 7, 2, 9, 8, 3, 4],
        [0, 4, 1, 2, 3, 7, 3, 7, 5],
        [6, 9, 2, 6, 0, 0, 5, 1, 4]],

       [[4, 2, 0, 1, 6, 7, 1, 0, 8],
        [1, 5, 3, 6, 4, 2, 4, 8, 3],
        [7, 4, 9, 9, 1, 9, 7, 3, 1],
        [3, 6, 1, 2, 5, 4, 1, 3, 0],
        [3, 3, 6, 6, 9, 8, 4, 2, 8],
        [7, 9, 1, 3, 0, 2, 0, 7, 4],
        [6, 7, 9, 3, 0, 2, 1, 9, 2],
        [1, 0, 3, 4, 7, 8, 1, 6, 5],
        [4, 4, 7, 8, 3, 7, 0, 4, 7]]])

I would like to get 3D chunks of shape 2 × 3 × 3 using a moving window along latter two dimensions (in this case 9 × 9). The size of the first dimension (I'd call it "depth") is arbitrary. The example of the first chunk would be:

>>> array([
       [[np.nan, np.nan, np.nan],
        [np.nan, 4, 5],
        [np.nan, 9, 2]],

        [[np.nan, np.nan, np.nan],
        [np.nan, 4, 2],
        [np.nan, 1, 5]]])

The second would be:

>>> array([
       [[np.nan, np.nan, np.nan],
        [4, 5, 1],
        [9, 2, 2]],

        [[np.nan, np.nan, np.nan],
        [4, 2, 0],
        [1, 5, 3]]])

And so on...

I later need to apply a more complicated function to these chunks, not a simple average or such, so I would appreciate a new array with them (I guess that is quite memory intensive, is there a different approach? Possibly vectorized? But it's not necessary)

I tried applying np.lib.stride_tricks.as_strided to my case as in #44305987 and played around with fancy indexing as in #15722324, but did not achieve the desired result.

Thanks!

janchytry
  • 151
  • 1
  • 9
  • 1
    I posted a general solution to something similar [here](https://stackoverflow.com/q/53263678/1782792) - you would just need to [`pad`](https://numpy.org/doc/stable/reference/generated/numpy.pad.html) with NaN beforehand. – jdehesa Jul 22 '20 at 09:30
  • Well served @jdehesa! I need to learn more about stride tricks. – janchytry Jul 22 '20 at 09:33

1 Answers1

4

You could use skimage.util.view_as_windows for this. Since it appears that you want a minimum size of 2 elements for those windowed views, you could assign the array to a larger array of np.nan, and take the strided view of the resulting array:

from skimage.util import view_as_windows

i,j,k= a.shape
a_exp = np.full((i,j+2,k+2), np.nan)
a_exp[:,1:j+1,1:k+1] = a

Or you could also do the same with np.pad:

a_exp = np.pad(a.astype('float'), 
               pad_width=((0,0),(1,1),(1,1)), 
               constant_values=np.nan)

And view as strided with:

out = view_as_windows(a_exp, (a.shape[0],3,3))

out
array([[[[[[nan, nan, nan],
           [nan,  4.,  5.],
           [nan,  9.,  2.]],

          [[nan, nan, nan],
           [nan,  4.,  2.],
           [nan,  1.,  5.]]],


         [[[nan, nan, nan],
           [ 4.,  5.,  1.],
           [ 9.,  2.,  2.]],

          [[nan, nan, nan],
           [ 4.,  2.,  0.],
           [ 1.,  5.,  3.]]],
         ...
yatu
  • 75,195
  • 11
  • 47
  • 89
  • Thanks for the `skimage.util.view_as_windows`, extremely useful! I've been cruising around scikit but never got to this method. – janchytry Jul 22 '20 at 09:14
  • 1
    To be agnostic, use `view_as_windows(a_exp, (1,3,3))[...,0,:,:]` and then maybe `transpose(1,2,0,3,4)[0,1]` if OP needs in that first, second chunk order. More info - https://stackoverflow.com/a/51890064/ – Divakar Jul 22 '20 at 09:15
  • @Divakar I'll look at that, thanks. I am just interested, do you think you could point to a NumPy solution to the problem? To my best knowledge, `np.convolve` sorts only 1D arrays... – janchytry Jul 22 '20 at 09:18
  • @janchytry Convolution is for summing. Are you summing those windows? If you are looking to sum, you might want multidim convolution - https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.convolve.html#scipy.ndimage.convolve – Divakar Jul 22 '20 at 09:23
  • @Divakar No, I do not, that's probably my misconduct of simply moving window vs. convolution kernel. The function I want to apply later is not expressible by a kernel. – janchytry Jul 22 '20 at 09:30
  • Added as alternative, though unsure of why the strided view could be somewhow dependant, or why do you suggest otherwise? @Divakar – yatu Jul 22 '20 at 09:38
  • @yatu I meant agnostic to make it work irrespective of array shape. I suppose you could altternatively do `view_as_windows(a_exp, (a.shape[0],3,3))`. – Divakar Jul 22 '20 at 10:38