花式索引与元组答案

【问题标题】：Fancy indexing with tuples花式索引与元组
【发布时间】：2014-01-10 12:31:56
【问题描述】：

假设我在 numpy 中有一个 100x100 的数组，我想从这个数组中选择 10 个 (x*x) 的随机块像素并同时更改这些块的值。为每个块索引切片的最佳方法是什么？一个理想的解决方案是按照下面的思路，在元组对之间进行切片。

A = np.ones(100,100)
blockSize = 10

numBlocks = 15

blockCenter_Row = tuple(np.random.randint(blockSize,high=(100-blockSize),size=numBlocks))
blockCenter_Col = tuple(np.random.randint(blockSize,high=(100-blockSize),size=numBlocks))
rowLeft_Boundary = tuple((i-blockSize/2) for i in blockCenter_Row)
rowRight_Boundary = tuple((i+blockSize/2) for i in blockCenter_Row)
colLower_Boundary = tuple((i-blockSize/2) for i in blockCenter_Row)
colUpper_Boundary = tuple((i+blockSize/2) for i in blockCenter_Row)

for value in range(10):
    A[rowLeft_Boundary:rowRight_Boundary,colLower_Boundary:colUpper_Boundary] = value

【问题讨论】：

块的重叠怎么处理？
我不确定您所说的“最佳方式”是什么意思：写起来不那么麻烦，索引速度更快，创建索引速度更快，...？因为除了最后一行中*_Boundary 变量上缺少索引之外，这似乎工作正常。
我在这个例子中犯了一点错误，理想情况下应该计算块的中心，以便在最终实现中没有重叠。

标签： python arrays numpy tuples slice

【解决方案1】：

如果块可以重叠，我认为您可以使用as_strided() 来解决问题。

import pylab as pl
from numpy.lib.stride_tricks import as_strided

blockSize = 10
numBlocks = 15
n = 100
a = np.zeros((n, n))

itemsize = a.dtype.itemsize
new_shape = n-blockSize+1, n-blockSize+1, blockSize, blockSize
new_stride = itemsize*n, itemsize, itemsize*n, itemsize
b = as_strided(a, shape=new_shape, strides=new_stride)

idx0 = np.random.randint(0, b.shape[0], numBlocks)
idx1 = np.random.randint(0, b.shape[1], numBlocks)

b[idx0, idx1, :, :] = np.random.rand(numBlocks, blockSize, blockSize)*3 + np.arange(numBlocks).reshape(-1, 1, 1)

pl.imshow(a, cmap="gray", interpolation="nearest")

这是输出：

【讨论】：

+1 作为一个小问题，使用new_stride = a.strides * 2 更安全，无论原始数组的内存布局如何，它都可以工作。
这是一种适用于较小阵列的可能解决方案。然而，在最终的应用程序中，我的 Ny 和 Nx 分别为 1080 和 1920，这会产生大约 5GB 左右的数组，具体取决于块大小。
你试过 1080x1920 数组的代码吗？ as_strided() 创建一个视图而不是副本。这段代码使用的内存是：numBlocks x blockSize x blockSize x sizeof(dtype).
是的。任何超过 ~1000x1000 的东西在运行时都会崩溃。