24

I have a 2-d numpy array as follows:

a = np.array([[1,5,9,13],
              [2,6,10,14],
              [3,7,11,15],
              [4,8,12,16]]

I want to extract it into patches of 2 by 2 sizes with out repeating the elements.

The answer should exactly be the same. This can be 3-d array or list with the same order of elements as below:

[[[1,5],
 [2,6]],   

 [[3,7],
 [4,8]],

 [[9,13],
 [10,14]],

 [[11,15],
 [12,16]]]

How can do it easily?

In my real problem the size of a is (36, 72). I can not do it one by one. I want programmatic way of doing it.

15
  • I updated my answer at stackoverflow.com/questions/26871083/…. Given that question and stackoverflow.com/questions/31494190/…, I think we can close this one as a dupe. Commented Jul 21, 2015 at 1:13
  • @WarrenWeckesser Can you show me HERE how you would extract the patches as I extracted mannually in my question? Commented Jul 21, 2015 at 1:33
  • @WarrenWeckesser It's not about calculating mean as in oyur answer Commented Jul 21, 2015 at 1:42
  • I've already improved my explanation of how the answer at stackoverflow.com/questions/26871083/… works. Did you see the part that begins "To generalize..."? There are two steps: reshape the array to a 4-d array, and then average. The reshaping part is the same as what you are asking, so I'd rather not duplicate that here. Commented Jul 21, 2015 at 1:45
  • @WarrenWeckesser I think your answer is more than sufficient for him to generalize to an answer (you gave the exact formula lol). It definitely puts my little baby python coder attempt to shame. I'm glad I looked at it. Commented Jul 21, 2015 at 1:47

3 Answers 3

29

Using scikit-image:

import numpy as np
from skimage.util import view_as_blocks

a = 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)))
Sign up to request clarification or add additional context in comments.

Comments

15

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]: a
Out[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 loop

In [44]: %timeit extract_blocks(a, (2,3), keep_as_view=False)
10 loops, best of 3: 29.1 ms per loop

Comments

7

Here's a rather cryptic numpy one-liner to generate your 3-d array, called result1 here:

In [60]: x
Out[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.shape
Out[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]])

3 Comments

Thanks this is what I was looking for!
@Warren , you deserve a beer
How would you reverse this operation?

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.