从NumPy数组中提取块或补丁

Question

我有一个二维numpy数组如下:

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

我想把它提取成2 * 2大小的块，并且不重复元素。

答案应该完全相同。这可以是3-d数组或列表，其元素顺序与以下相同:

[[[1,5],
 [2,6]],   
 [[3,7],
 [4,8]],
 [[9,13],
 [10,14]],
 [[11,15],
 [12,16]]]

怎么能轻松做到呢?

在我的实际问题中，a的大小是(36,72)。我不能一个一个地做。

Answer 1

使用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)))

Answer 2

您可以将np.reshape和np.swapaxes组合使用，如下-

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)

可以看到，有两种方法来使用它-与keep_as_view标志关闭(默认一个)或打开。对于keep_as_view = False，我们正在将交换轴重塑为3D的最终输出，而对于keep_as_view = True，我们将保持它为4D，这将是输入数组的视图，因此，在运行时几乎是免费的。我们将在稍后运行的示例案例中验证它。

样本情况下

让我们使用一个示例输入数组，像这样-

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]])

现在，让我们使用一些块大小进行测试。让我们使用块大小为(2,3)，并将视图标志关闭和打开-

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]]]])

用keep_as_view=True验证view

In [20]: np.shares_memory(a, extract_blocks(a, (2,3), keep_as_view=True))
Out[20]: True

让我们检查一个大数组上的性能，并验证前面讨论的虚拟自由运行时声明-

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

Answer 3

下面是一个相当神秘的numpy一行代码，用于生成3d数组，这里称为result1:

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就像一个二维数组的一维数组:

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]])

对不起，我现在没有时间详细说明它是如何工作的。也许以后…)

下面是一个使用嵌套列表推导式的不那么神秘的版本。在本例中，result2是一个二维numpy数组的python列表:

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]])