python - 过滤数组，存储邻接信息

Question

假设我有一个 (N, N) 形状的二维 array:

import numpy as np
my_array = np.random.random((N, N))

现在我只想对该数组的某些“单元格”进行一些计算，例如数组中央部分的单元格。为了避免对我不感兴趣的单元格进行计算，我通常在这里做的是创建一个 bool 掩码，本着这种精神:

my_mask = np.zeros_like(my_array, bool)
my_mask[40:61,40:61] = True
my_array[my_mask] = some_twisted_computations(my_array[my_mask])

但是如果 some_twisted_computations() 涉及相邻单元格的值(如果它们在 mask 内)怎么办？在性能方面，创建一个形状为 (len(my_mask), 4) 的“邻接数组”是否是个好主意，将 4 个相邻单元格的索引存储在平面 my_array[mask] 我将在 some_twisted_computations() 中使用的数组？如果是，计算这种邻接数组的有效选项是什么？我应该切换到较低级别的语言/其他数据结构吗？

我的真实世界数组形状在 (1000,1000,1000) 左右，掩码仅涉及这些值的一小部分 (~100000)，并且具有相当复杂的几何形状。我希望我的问题是有道理的...

编辑:我制定的非常肮脏和缓慢的解决方案:

wall = mask

i = 0

top_neighbors = []
down_neighbors = []
left_neighbors = []
right_neighbors = []
indices = []

for index, val in np.ndenumerate(wall):
    if not val:
        continue
    indices += [index]
    if wall[index[0] + 1, index[1]]:
        down_neighbors += [(index[0] + 1, index[1])]
    else:
        down_neighbors += [i]
    if wall[index[0] - 1, index[1]]:
        top_neighbors += [(index[0] - 1, index[1])]
    else:
        top_neighbors += [i]
    if wall[index[0], index[1] - 1]:
        left_neighbors += [(index[0], index[1] - 1)]
    else:
        left_neighbors += [i]
    if wall[index[0], index[1] + 1]:
        right_neighbors += [(index[0], index[1] + 1)]
    else:
        right_neighbors += [i]
    i += 1


top_neighbors = [i if type(i) is int else indices.index(i) for i in top_neighbors]
down_neighbors = [i if type(i) is int else indices.index(i) for i in down_neighbors]
left_neighbors = [i if type(i) is int else indices.index(i) for i in left_neighbors]
right_neighbors = [i if type(i) is int else indices.index(i) for i in right_neighbors]

Answer 1

最佳答案可能取决于您要进行的计算的性质。例如，如果它们可以表示为相邻像素的总和，那么像 np.convolve 或 scipy.signal.fftconvolve 这样的东西可能是一个非常好的解决方案。

对于有效生成相邻索引数组的特定问题，您可以尝试这样的事情:

x = np.random.rand(100, 100)
mask = x > 0.9

i, j = np.where(mask)

i_neighbors = i[:, np.newaxis] + [0, 0, -1, 1]
j_neighbors = j[:, np.newaxis] + [-1, 1, 0, 0]

# need to do something with the edge cases
# the best choice will depend on your application
# here we'll change out-of-bounds neighbors to the
# central point itself.
i_neighbors = np.clip(i_neighbors, 0, 99)
j_neighbors = np.clip(j_neighbors, 0, 99)

# compute some vectorized result over the neighbors
# as a concrete example, here we'll do a standard deviation
result = x[i_neighbors, j_neighbors].std(axis=1)

结果是对应于屏蔽区域的值数组，包含相邻值的标准差。希望该方法适用于您想到的任何具体问题!

---

编辑:考虑到上面编辑过的问题，以下是如何调整我的回答以以矢量化方式生成索引数组:

x = np.random.rand(100, 100)
mask = x > -0.9

i, j = np.where(mask)

i_neighbors = i[:, np.newaxis] + [0, 0, -1, 1]
j_neighbors = j[:, np.newaxis] + [-1, 1, 0, 0]
i_neighbors = np.clip(i_neighbors, 0, 99)
j_neighbors = np.clip(j_neighbors, 0, 99)

indices = np.zeros(x.shape, dtype=int)
indices[mask] = np.arange(len(i))

neighbor_in_mask = mask[i_neighbors, j_neighbors]

neighbors = np.where(neighbor_in_mask,
                     indices[i_neighbors, j_neighbors],
                     np.arange(len(i))[:, None])

left_indices, right_indices, top_indices, bottom_indices = neighbors.T