NumPy数组中固定大小子矩阵的索引

Question

NumPy数组中固定大小子矩阵的索引

14

我正在实现一个算法，需要查看在一个（严格二维）numpy数组中不重叠的连续子矩阵。例如，对于12x12的数组：

>>> a = np.random.randint(20, size=(12, 12)); a
array([[ 4,  0, 12, 14,  3,  8, 14, 12, 11, 18,  6,  6],
       [15, 13,  2, 18, 15, 15, 16,  2,  9, 16,  6,  4],
       [18, 18,  3,  8,  1, 15, 14, 13, 13, 13,  7,  0],
       [ 1,  9,  3,  6,  0,  4,  3, 15,  0,  9, 11, 12],
       [ 5, 15,  5,  6,  4,  4, 18, 13, 10, 17, 11,  8],
       [13, 17,  8, 15, 17, 12,  7,  1, 13, 15,  0, 18],
       [ 2,  1, 11, 12,  3, 16, 11,  9, 10, 15,  4, 16],
       [19, 11, 10,  7, 10, 19,  7, 13, 11,  9, 17,  8],
       [14, 14, 17,  0,  0,  0, 11,  1, 10, 14,  2,  7],
       [ 6, 15,  6,  7, 15, 19,  2,  4,  6, 16,  0,  3],
       [ 5, 10,  7,  5,  0,  8,  5,  8,  9, 14,  4,  3],
       [17,  2,  0,  3, 15, 10, 14,  1,  0,  7, 16,  2]])

并且在查看3x3子矩阵时，我希望第一个3x3子矩阵来自左上角：

>>> a[0:3, 0:3]
array([[ 4,  0, 12],
       [15, 13,  2],
       [18, 18,  3]])

下一个将被给予的是a[0:3, 3:6]等等。如果每行或每列的最后一组索引超出了数组的末尾，也没关系 - numpy的行为只需提供在切片中存在的部分即可。

我想要一种能够以编程方式生成这些切片索引的方法，适用于任意大小的矩阵和子矩阵。我目前有这个：

size = 3
x_max = a.shape[0]
xcoords = range(0, x_max, size)
xcoords = zip(xcoords, xcoords[1:])

类似地，生成y_coords的方法也是相似的，从而可以通过itertools.product(xcoords, ycoords)得到索引系列。

我的问题是：是否有更直接的方法来完成这个任务，比如使用numpy.mgrid或其他一些numpy技巧？

- lvc

3个回答

5

我把这个答案添加到一个老问题中，因为编辑使这个问题重新浮现。这里有一种计算块的替代方法：

size = 3
lenr, lenc = int(a.shape[0]/size), int(a.shape[1]/size)

t = a.reshape(lenr,size,lenc,size).transpose(0, 2, 1, 3)

分析结果表明这是最快的。使用Python 3.5进行分析，将结果传递给array()以保持兼容性，因为在3.5中map返回一个迭代器。

reshape/transpose:   643 ns per loop
reshape/index:       45.8 µs per loop
Map/split:           10.3 µs per loop

有趣的是，使用迭代器版本的map更快。然而，无论如何，使用reshape和transpose是最快的。

- user2699

毫无疑问，这是解决此问题的正确方法；毫无疑问。 - Eelco Hoogendoorn

5

您可以使用以下一行代码:

r = 3
c = 3
lenr = a.shape[0]/r
lenc = a.shape[1]/c
np.array([a[i*r:(i+1)*r,j*c:(j+1)*c] for (i,j) in np.ndindex(lenr,lenc)]).reshape(lenr,lenc,r,c)

- Saullo G. P. Castro

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Geoff · Accepted Answer

获取索引

下面是一个快速获取特定 size x size 块的方法：

base = np.arange(size) # Just the base set of indexes
row = 1                # Which block you want
col = 0                
block = a[base[:, np.newaxis] + row * size, base + col * size]

如果你愿意，你可以像你的xcoords一样建立类似的矩阵：

y, x = np.mgrid[0:a.shape[0]/size, 0:a.shape[1]/size]
y_coords = y[..., np.newaxis] * size + base
x_coords = x[..., np.newaxis] * size + base

那么你可以像这样访问一个代码块：

block = a[y_coords[row, col][:, np.newaxis], x_coords[row, col]]

直接获取块

如果您只想获取块（而不是块条目的索引），建议使用np.split函数（两次）：

blocks = map(lambda x : np.split(x, a.shape[1]/size, 1), # Split the columns
                        np.split(a, a.shape[0]/size, 0)) # Split the rows

那么你会有一个大小为size x size的2D列表块：

>>> blocks[0][0]
array([[ 4,  0, 12],
       [15, 13,  2],
       [18, 18,  3]])

>>> blocks[1][0]
array([[ 1,  9,  3],
       [ 5, 15,  5],
       [13, 17,  8]])

你可以将此转换为numpy数组并使用与上文相同的索引方式：

>>> blocks = np.array(blocks)
>>> blocks.shape
(4, 4, 3, 3)