Numpy split cube into cubes

Question

There is a function np.split() which can split an array along 1 axis. I was wondering if there was a multi axis version where you can split along axes (0,1,2) for example.

Answer 1

Suppose the cube has shape (W, H, D) and you wish to break it up into N little cubes of shape (w, h, d). Since NumPy arrays have axes of fixed length, w must evenly divide W, and similarly for h and d.

Then there is a way to reshape the cube of shape (W, H, D) into a new array of shape (N, w, h, d).

For example, if arr = np.arange(4*4*4).reshape(4,4,4) (so (W,H,D) = (4,4,4)) and we wish to break it up into cubes of shape (2,2,2), then we could use

In [283]: arr.reshape(2,2,2,2,2,2).transpose(0,2,4,1,3,5).reshape(-1,2,2,2)
Out[283]: 
array([[[[ 0,  1],
         [ 4,  5]],

        [[16, 17],
         [20, 21]]],

...
       [[[42, 43],
         [46, 47]],

        [[58, 59],
         [62, 63]]]])

The idea here is to add extra axes to the array which sort of act as place markers:

 number of repeats act as placemarkers
 o---o---o
 |   |   |
 v   v   v
(2,2,2,2,2,2)
   ^   ^   ^
   |   |   |
   o---o---o
   newshape

We can then reorder the axes (using transpose) so that the number of repeats comes first, and the newshape comes at the end:

arr.reshape(2,2,2,2,2,2).transpose(0,2,4,1,3,5)

And finally, call reshape(-1, w, h, d) to squash all the placemarking axes into a single axis. This produces an array of shape (N, w, h, d) where N is the number of little cubes.

---

The idea used above is a generalization of this idea to 3 dimensions. It can be further generalized to ndarrays of any dimension:

import numpy as np
def cubify(arr, newshape):
    oldshape = np.array(arr.shape)
    repeats = (oldshape / newshape).astype(int)
    tmpshape = np.column_stack([repeats, newshape]).ravel()
    order = np.arange(len(tmpshape))
    order = np.concatenate([order[::2], order[1::2]])
    # newshape must divide oldshape evenly or else ValueError will be raised
    return arr.reshape(tmpshape).transpose(order).reshape(-1, *newshape)

print(cubify(np.arange(4*6*16).reshape(4,6,16), (2,3,4)).shape)
print(cubify(np.arange(8*8*8*8).reshape(8,8,8,8), (2,2,2,2)).shape)

yields new arrays of shapes

(16, 2, 3, 4)
(256, 2, 2, 2, 2)

---

To "uncubify" the arrays:

def uncubify(arr, oldshape):
    N, newshape = arr.shape[0], arr.shape[1:]
    oldshape = np.array(oldshape)    
    repeats = (oldshape / newshape).astype(int)
    tmpshape = np.concatenate([repeats, newshape])
    order = np.arange(len(tmpshape)).reshape(2, -1).ravel(order='F')
    return arr.reshape(tmpshape).transpose(order).reshape(oldshape)

---

Here is some test code to check that cubify and uncubify are inverses.

import numpy as np
def cubify(arr, newshape):
    oldshape = np.array(arr.shape)
    repeats = (oldshape / newshape).astype(int)
    tmpshape = np.column_stack([repeats, newshape]).ravel()
    order = np.arange(len(tmpshape))
    order = np.concatenate([order[::2], order[1::2]])
    # newshape must divide oldshape evenly or else ValueError will be raised
    return arr.reshape(tmpshape).transpose(order).reshape(-1, *newshape)

def uncubify(arr, oldshape):
    N, newshape = arr.shape[0], arr.shape[1:]
    oldshape = np.array(oldshape)    
    repeats = (oldshape / newshape).astype(int)
    tmpshape = np.concatenate([repeats, newshape])
    order = np.arange(len(tmpshape)).reshape(2, -1).ravel(order='F')
    return arr.reshape(tmpshape).transpose(order).reshape(oldshape)

tests = [[np.arange(4*6*16), (4,6,16), (2,3,4)],
         [np.arange(8*8*8*8), (8,8,8,8), (2,2,2,2)]]

for arr, oldshape, newshape in tests:
    arr = arr.reshape(oldshape)
    assert np.allclose(uncubify(cubify(arr, newshape), oldshape), arr)
    # cuber = Cubify(oldshape,newshape)
    # assert np.allclose(cuber.uncubify(cuber.cubify(arr)), arr)

Answer 2

In addition to my extra question to @unutbu's answer I think I have got the reverse to work (in case you want to split a cube into cubes, apply a function to each one and then combine them back).

import numpy as np
import pdb
np.set_printoptions(precision=3,linewidth=300)

class Cubify():
    def __init__(self,oldshape,newshape):
        self.newshape = np.array(newshape)
        self.oldshape = np.array(oldshape)
        self.repeats = (oldshape / newshape).astype(int)
        self.tmpshape = np.column_stack([self.repeats, newshape]).ravel()
        order = np.arange(len(self.tmpshape))
        self.order = np.concatenate([order[::2], order[1::2]])
        self.reverseOrder = self.order.copy()
        self.reverseOrder = np.arange(len(self.tmpshape)).reshape(2, -1).ravel(order='F')
        self.reverseReshape = np.concatenate([self.repeats,self.newshape])

    def cubify(self,arr):
        # newshape must divide oldshape evenly or else ValueError will be raised
        return arr.reshape(self.tmpshape).transpose(self.order).reshape(-1, *self.newshape)

    def uncubify(self,arr):
        return arr.reshape(self.reverseReshape).transpose(self.reverseOrder).reshape(self.oldshape)

if __name__ == "__main__":
    N = 9
    x = np.arange(N**3).reshape(N,N,N)
    oldshape = x.shape
    newshape = np.array([3,3,3])
    cuber = Cubify(oldshape,newshape)
    out = cuber.cubify(x)
    back = cuber.uncubify(out)