Extract blocks or patches from NumPy Array
Using scikit-image:
import numpy as npfrom skimage.util import view_as_blocksa = np.array([[1,5,9,13], [2,6,10,14], [3,7,11,15], [4,8,12,16]])print(view_as_blocks(a, (2, 2)))
You can achieve it with a combination of np.reshape
and np.swapaxes
like so -
def extract_blocks(a, blocksize, keep_as_view=False): M,N = a.shape b0, b1 = blocksize if keep_as_view==0: return a.reshape(M//b0,b0,N//b1,b1).swapaxes(1,2).reshape(-1,b0,b1) else: return a.reshape(M//b0,b0,N//b1,b1).swapaxes(1,2)
As can be seen there are two ways to use it - With keep_as_view
flag turned off (default one) or on. With keep_as_view = False
, we are reshaping the swapped-axes to a final output of 3D
, while with keep_as_view = True
, we will keep it 4D and that will be a view into the input array and hence, virtually free on runtime. We will verify it with a sample case run later on.
Sample cases
Let's use a sample input array, like so -
In [94]: aOut[94]: array([[2, 2, 6, 1, 3, 6], [1, 0, 1, 0, 0, 3], [4, 0, 0, 4, 1, 7], [3, 2, 4, 7, 2, 4], [8, 0, 7, 3, 4, 6], [1, 5, 6, 2, 1, 8]])
Now, let's use some block-sizes for testing. Let's use a blocksize of (2,3)
with the view-flag turned off and on -
In [95]: extract_blocks(a, (2,3)) # Blocksize : (2,3)Out[95]: array([[[2, 2, 6], [1, 0, 1]], [[1, 3, 6], [0, 0, 3]], [[4, 0, 0], [3, 2, 4]], [[4, 1, 7], [7, 2, 4]], [[8, 0, 7], [1, 5, 6]], [[3, 4, 6], [2, 1, 8]]])In [48]: extract_blocks(a, (2,3), keep_as_view=True)Out[48]: array([[[[2, 2, 6], [1, 0, 1]], [[1, 3, 6], [0, 0, 3]]], [[[4, 0, 0], [3, 2, 4]], [[4, 1, 7], [7, 2, 4]]], [[[8, 0, 7], [1, 5, 6]], [[3, 4, 6], [2, 1, 8]]]])
Verify view
with keep_as_view=True
In [20]: np.shares_memory(a, extract_blocks(a, (2,3), keep_as_view=True))Out[20]: True
Let's check out performance on a large array and verify the virtually free runtime claim as discussed earlier -
In [42]: a = np.random.rand(2000,3000)In [43]: %timeit extract_blocks(a, (2,3), keep_as_view=True)1000000 loops, best of 3: 801 ns per loopIn [44]: %timeit extract_blocks(a, (2,3), keep_as_view=False)10 loops, best of 3: 29.1 ms per loop
Here's a rather cryptic numpy one-liner to generate your 3-d array, called result1
here:
In [60]: xOut[60]: array([[2, 1, 2, 2, 0, 2, 2, 1, 3, 2], [3, 1, 2, 1, 0, 1, 2, 3, 1, 0], [2, 0, 3, 1, 3, 2, 1, 0, 0, 0], [0, 1, 3, 3, 2, 0, 3, 2, 0, 3], [0, 1, 0, 3, 1, 3, 0, 0, 0, 2], [1, 1, 2, 2, 3, 2, 1, 0, 0, 3], [2, 1, 0, 3, 2, 2, 2, 2, 1, 2], [0, 3, 3, 3, 1, 0, 2, 0, 2, 1]])In [61]: result1 = x.reshape(x.shape[0]//2, 2, x.shape[1]//2, 2).swapaxes(1, 2).reshape(-1, 2, 2)
result1
is like a 1-d array of 2-d arrays:
In [68]: result1.shapeOut[68]: (20, 2, 2)In [69]: result1[0]Out[69]: array([[2, 1], [3, 1]])In [70]: result1[1]Out[70]: array([[2, 2], [2, 1]])In [71]: result1[5]Out[71]: array([[2, 0], [0, 1]])In [72]: result1[-1]Out[72]: array([[1, 2], [2, 1]])
(Sorry, I don't have time at the moment to give a detailed breakdown of how it works. Maybe later...)
Here's a less cryptic version that uses a nested list comprehension. In this case, result2
is a python list of 2-d numpy arrays:
In [73]: result2 = [x[2*j:2*j+2, 2*k:2*k+2] for j in range(x.shape[0]//2) for k in range(x.shape[1]//2)]In [74]: result2[5]Out[74]: array([[2, 0], [0, 1]])In [75]: result2[-1]Out[75]: array([[1, 2], [2, 1]])