【问题标题】:Slice a submatrix with center element indices切片具有中心元素索引的子矩阵
【发布时间】:2019-11-06 03:33:48
【问题描述】:

给定一个矩阵 A、一个行索引列表和一个列索引列表,如何有效地提取以行索引和列索引为中心的大小为 k 的平方子矩阵?

例如:

A = array([[12,  6, 14,  8,  4,  1],
       [18, 13,  8, 10,  9, 19],
       [ 8, 15,  6,  5,  6, 18],
       [ 3,  0,  2, 14, 13, 12],
       [ 4,  4,  5, 19,  0, 14],
       [16,  8,  7,  7, 11,  0],
       [ 3, 11,  2, 19, 11,  5],
       [ 4,  2,  1,  9, 12, 12]])
r = np.array([2, 5])
c = np.array([3, 2])
k = 3

输出应该是A[1:4, 2:5] 和A[4:7, 1:4]。所以基本上,输出是大小为 kxk 的平方子矩阵,以 [r,c] 元素为中心(在这种情况下为 A[2,3] 和 A[5,2])

如何高效而优雅地做到这一点?谢谢

【问题讨论】:

  • 我认为没有什么特别的技巧。对于r 和c 中的每一对值,确定相关的切片(只是一些基本的数学运算),然后进行切片。
  • 是的,当然可以。但是如果 r 和 c 的长度非常大,那么逐个循环遍历每个 case 可能会很慢。
  • 这是viewing 数组作为一堆(可能重叠)窗口的一种方式。您可以从中选择一个子集。它使用as_strided 函数。它很有效,但有点难以理解和正确操作。
  • 那么,所有子矩阵的形状都一样吗?

标签: python numpy numpy-ndarray array-broadcasting numpy-slicing


【解决方案1】:

你的意思是这样的?

for x,y in zip(r,c):
    s = k // 2
    print("position:",[x - s,x + s + 1], [y - s,y + s + 1])
    print(A[x - s:x + s + 1,y - s:y + s + 1])
    print()

输出:

position: [1, 4] [2, 5]
[[ 8 10  9]
 [ 6  5  6]
 [ 2 14 13]]

position: [4, 7] [1, 4]
[[ 4  5 19]
 [ 8  7  7]
 [11  2 19]]

注意k在这里应该是奇数

【讨论】:

  • 这将完成这项工作,但效率不高,因为它需要一个 for 循环来遍历 r 和 c 中的每个元素。当 r 和 c 的长度很大时,它会变慢。我想知道是否有一种很好的切片方法可以有效地获得结果,即使 r 和 c 中的元素数量很大。
【解决方案2】:

对于子矩阵具有相同形状的情况,我们可以获得滑动窗口,然后沿着行和列对那些具有起始索引的窗口进行索引,以获得所需的输出。要获得这些窗口,我们可以利用基于scikit-image's view_as_windows 的np.lib.stride_tricks.as_strided。 More info on use of as_strided based view_as_windows -

from skimage.util.shape import view_as_windows

# Get all sliding windows
w = view_as_windows(A,(k,k))

# Select relevant ones for final o/p
out = w[r-k//2,c-k//2]

【讨论】:

  • 谢谢。这是一个非常好的解决方案。还有一个问题,因为view_as_windows 的返回值返回了二维数组中的所有切片窗口。这可能会消耗额外的内存。有没有办法只在给定的 r 和 c 上获得切片窗口?
  • @Tao 这些是滑动窗口,作为输入数组的视图,因此没有使用额外的内存。所以,我们在那里很好。您可以阅读帖子中的最后第三个链接以获取更多信息。
  • 谢谢。但是当我打印 w 的形状时。它说 (6, 4, 3, 3)。 A 矩阵的形状为 (8, 6)。所以我认为确实生成了所有可能的滑动窗口。
  • @Tao 使用np.shares_memory(A,w)。它将打印True,这意味着w 是A 的视图。有关np.shares_memory 的更多信息,您可以谷歌搜索其文档。
  • 我明白了。谢谢。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2022-07-10
  • 2013-01-03
  • 1970-01-01
  • 1970-01-01
  • 2012-09-23
相关资源
最近更新 更多