numpy ndarray 中的高效邻域搜索而不是嵌套条件 for 循环
Efficient neighbourhood search in numpy ndarray instead of nested conditional for loops
尽管问题有很多实例：“嵌套 for 循环的 numpy 替代方法是什么”，但我无法为我的案例找到合适的答案。开始了：
我有一个 3D numpy 数组，背景为“0”，前景为其他整数。我想找到并存储落在预定义掩码（定义距参考节点给定距离的球体）内的前景体素。我已经使用嵌套 'for' 循环和一系列 'if' 条件成功完成了任务，如下所示。我正在寻找一种更高效、更紧凑的替代方法来避免这种邻域搜索算法的循环和长条件。
示例输入数据：
import numpy as np
im = np.array([[[ 60, 54, 47, 52, 57, 53, 46, 48]
, [ 60, 57, 53, 53, 54, 53, 50, 55]
, [ 60, 63, 56, 58, 59, 57, 50, 50]
, [ 70, 70, 64, 69, 74, 72, 64, 47]
, [ 73, 76, 77, 80, 82, 76, 58, 37]
, [ 85, 85, 86, 86, 78, 62, 38, 20]
, [ 94, 94, 92, 78, 54, 33, 16, 255]
, [ 94, 90, 72, 51, 32, 19, 255, 255]
, [ 65, 53, 29, 18, 255, 255, 255, 255]
, [ 29, 22, 255, 255, 255, 255, 255, 0]]
, [[ 66, 67, 70, 69, 75, 73, 72, 63]
, [ 68, 70, 73, 74, 78, 80, 74, 53]
, [ 75, 87, 87, 83, 89, 86, 61, 33]
, [ 81, 89, 88, 98, 99, 77, 41, 18]
, [ 84, 94, 100, 100, 82, 49, 21, 255]
, [ 99, 101, 92, 75, 48, 25, 255, 255]
, [ 93, 77, 52, 32, 255, 255, 255, 255]
, [ 52, 40, 25, 255, 255, 255, 255, 255]
, [ 23, 16, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]]
, [[ 81, 83, 92, 101, 101, 83, 49, 19]
, [ 86, 96, 103, 103, 95, 64, 28, 255]
, [ 94, 103, 107, 98, 79, 41, 255, 255]
, [101, 103, 98, 79, 51, 28, 255, 255]
, [102, 97, 76, 49, 27, 255, 255, 255]
, [ 79, 62, 35, 21, 255, 255, 255, 255]
, [ 33, 23, 15, 255, 255, 255, 255, 255]
, [ 16, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]]
, [[106, 107, 109, 94, 58, 26, 15, 255]
, [110, 104, 90, 66, 37, 19, 255, 255]
, [106, 89, 61, 35, 22, 255, 255, 255]
, [ 76, 56, 34, 19, 255, 255, 255, 255]
, [ 40, 27, 18, 255, 255, 255, 255, 255]
, [ 17, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]]
, [[ 68, 51, 33, 19, 255, 255, 255, 255]
, [ 45, 34, 20, 255, 255, 255, 255, 255]
, [ 28, 18, 255, 255, 255, 255, 255, 255]
, [ 17, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 0, 0, 0, 0, 0, 0, 0]]
, [[255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 255, 0, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]]
, [[255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 255, 0, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 255, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]]
, [[255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 255, 0, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 255, 0, 0, 0, 0, 0, 0]
, [255, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]]])
实现方法：
[Z,Y,X]=im.shape
RN = np.array([3,4,4])
################Loading Area search
rad = 3
a,b,c = RN
x,y,z = np.ogrid[-c:Z-c,-b:Y-b,-a:X-a]
neighborMask = x*x + y*y + z*z<= rad*rad
noNodeMask = im > 0
mask = np.logical_and(neighborMask, noNodeMask)
imtemp = im.copy()
imtemp[mask] = -1
for i in range (X):
for j in range (Y):
for k in range (Z):
if imtemp[i,j,k]==-1:
if i in (0, X-1) or j in (0, Y-1) or k in (0, Z-1):
imtemp[i,j,k]=-2
elif imtemp[i+1,j,k] == 0 or imtemp[i-1,j,k] == 0 or imtemp[i,j+1,k] == 0 or imtemp[i,j-1,k] == 0 or imtemp[i,j,k+1] == 0 or imtemp[i,j,k-1] == 0:
imtemp[i,j,k]=-2
LA = np.argwhere(imtemp==-2)
上述示例代码生成的 LA 为：
In [90]:LA
Out[90]:
array([[4, 4, 0],
[4, 4, 6],
[4, 5, 5],
[4, 6, 4],
[4, 6, 5],
[4, 7, 3],
[5, 3, 5],
[5, 4, 4],
[5, 4, 5],
[5, 5, 3],
[5, 5, 4],
[5, 6, 2],
[5, 6, 3],
[6, 2, 4],
[6, 3, 3],
[6, 3, 4],
[6, 4, 2],
[6, 4, 3],
[6, 5, 1],
[6, 5, 2]])
Z 方向的切片（XY 平面实例）显示不同的未触及、屏蔽 (-1) 和目标 (-2) 节点：
由于您的循环仅使用直接 Numpy 索引，因此您可以使用 Numba 的 @njit 以 更高效的方式执行此操作 方式。
@njit
def compute_imtemp(imtemp, X, Y, Z):
for i in range (Z):
for j in range (Y-1):
for k in range (X-1):
if imtemp[i,j,k]==-1:
if i==(Z-1):
imtemp[i,j,k]=-2
elif imtemp[i+1,j,k] == 0 or imtemp[i-1,j,k] == 0 or imtemp[i,j+1,k] == 0 or imtemp[i,j-1,k] == 0 or imtemp[i,j,k+1] == 0 or imtemp[i,j,k-1] == 0:
imtemp[i,j,k]=-2
[...]
imtemp = im.copy()
imtemp[mask] = -1
compute_imtemp(imtemp, X, Y, Z)
LA = np.argwhere(imtemp==-2)
以下是我机器上的性能结果：
281 µs ± 1.43 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
776 ns ± 16.4 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
Numba 实施速度 362 倍。
请注意，由于编译，第一次调用 compute_imtemp 会很慢。克服这个问题的一种方法是在一个空的 Numpy 数组上调用 compute_imtemp。另一种方法是使用 Numba API 手动编译函数并显式向 Numba 提供类型。
问题陈述
你有一个大阵列中的“实心”形状。你从中切出一个球。您的目标是找到球内固体表面的指数。表面定义为与实体外部相邻且具有 6 点连通性的任何点。阵列的边缘也被认为是表面。
更快的循环解决方案
您已经计算出代表实体和球相交的掩码。您可以更优雅地计算掩码并将其转换为索引。我建议保持维度的顺序不变，而不是在不同的符号之间切换。例如，RN 的顺序会受到影响，并且您 运行 存在轴限制不匹配的风险。
RN = np.array([4, 4, 3])
rad = 3
im = ...
cutout = ((np.indices(im.shape) - RN.reshape(-1, 1, 1, 1))**2).sum(axis=0) <= rad**2
solid = im > 0
mask = solid & cutout
indices = np.argwhere(mask)
您也可以在不重塑 RN 的情况下获得切口，方法是
cutout = ((np.rollaxis(np.indices(im.shape, sparse=False), 0, 4) - RN)**2).sum(axis=-1) <= rad**2
计算索引的好处是您的循环不再需要很大。通过使用 argwhere，您基本上剥离了外部三个循环，只留下 if 语句进行循环。您还可以矢量化连接检查。这有一个很好的副作用，您可以为每个像素定义任意连接。
limit = np.array(im.shape) - 1 # Edge of `im`
connectivity = np.array([[ 1, 0, 0], # Add rows to determine connectivity
[-1, 0, 0],
[ 0, 1, 0],
[ 0, -1, 0],
[ 0, 0, 1],
[ 0, 0, -1]], dtype=indices.dtype)
index_mask = np.ones(len(indices), dtype=bool)
for n, ind in enumerate(indices):
if ind.all() and (ind < limit).all() and im[tuple((ind + connectivity).T)].all():
index_mask[n] = False
LA = indices[index_mask, :]
请注意 imtemp 根本没有意义。即使在您原来的循环中，您也可以直接操作 mask 。如果元素没有通过您的标准，您可以将元素设置为 False，而不是将它们设置为 -2。
我在这里做类似的事情。我们检查实际选择的每个索引，并确定它们中是否有任何一个在实体内部。这些索引从掩码中消除。然后根据掩码更新索引列表。
检查 ind.all() and (ind < limit).all() and im[tuple((ind + connectivity).T)].all() 是您在 or 条件下执行的操作的快捷方式，但相反（测试 non-surface 而不是表面）。
ind.all() 检查 none 的索引是否为零：即不在 top/front/left 表面上。(ind < limit).all() 检查 none 个索引等于相应的图像大小减一。im[tuple((ind + connectivity).T)].all() 检查连接像素中的 none 是否为零。 (ind + connectivity).T 是我们连接到的六个点的 (3, 6) 数组（当前由 (6, 3) connectivity 数组在每个轴中定义为 +/-1）。当你把它变成一个元组时，它就变成了一个奇特的索引，就好像你做了类似 im[x + connectivity[:, 0], y + connectivity[:, 1], z + connectivity[:, 2]] 的事情一样。索引中的逗号只是使它成为一个元组。我展示的方式更适合任意数量的维度。
通过所有三个测试的像素都在实体内部，并被移除。你当然可以编写循环来检查另一种方式，但是你必须改变你的掩码：
index_mask = np.zeros(len(indices), dtype=bool)
for n, ind in enumerate(indices):
if (ind == 0).any() or (ind == limit).any() or (im[tuple((ind + connectivity).T)] == 0).any():
index_mask[n] = True
LA = indices[index_mask, :]
手动循环无论如何都不理想。但是，它向您展示了如何缩短循环（可能缩短几个数量级），以及如何使用矢量化和广播定义任意连接，而不会陷入 hard-coding 它的困境。
完全矢量化的解决方案
上面的循环可以使用广播的魔力完全矢量化。我们可以批量添加 connectivity 并批量过滤结果，而不是遍历 indices 中的每一行。诀窍是添加足够的维度，将所有 connectivity 添加到 indices.
的 each 元素
您仍然希望忽略边缘的像素：
edges = (indices == 0).any(axis=-1) | (indices == limit).any(axis=-1)
conn_index = indices[~edges, None, :] + connectivity[None, ...]
index_mask = np.empty(len(indices), dtype=bool)
index_mask[edges] = True
index_mask[~edges] = (im[tuple(conn_index.T)] == 0).any(axis=0)
LA = indices[index_mask, :]
我希望使用 numba 编译的正确编写的循环会比此解决方案快得多，因为它将避免流水线操作的大部分开销。它不需要大的临时缓冲区或特殊处理。
TL;DR
# Parameters
RN = np.array([4, 4, 3])
rad = 3
im = ...
# Find subset of interest
cutout = ((np.indices(im.shape) - RN.reshape(-1, 1, 1, 1))**2).sum(axis=0) <= rad**2
solid = im > 0
# Convert mask to indices
indices = np.argwhere(solid & cutout)
# Find image edges among indices
edges = (indices == 0).any(axis=-1) | (indices == limit).any(axis=-1)
# Connectivity elements for non-edge pixels
conn_index = indices[~edges, None, :] + connectivity[None, ...]
# Mask the valid surface pixels
index_mask = np.empty(len(indices), dtype=bool)
index_mask[edges] = True
index_mask[~edges] = (im[tuple(conn_index.T)] == 0).any(axis=0)
# Final result
LA = indices[index_mask, :]
尽管问题有很多实例：“嵌套 for 循环的 numpy 替代方法是什么”，但我无法为我的案例找到合适的答案。开始了：
我有一个 3D numpy 数组，背景为“0”，前景为其他整数。我想找到并存储落在预定义掩码（定义距参考节点给定距离的球体）内的前景体素。我已经使用嵌套 'for' 循环和一系列 'if' 条件成功完成了任务，如下所示。我正在寻找一种更高效、更紧凑的替代方法来避免这种邻域搜索算法的循环和长条件。
示例输入数据：
import numpy as np
im = np.array([[[ 60, 54, 47, 52, 57, 53, 46, 48]
, [ 60, 57, 53, 53, 54, 53, 50, 55]
, [ 60, 63, 56, 58, 59, 57, 50, 50]
, [ 70, 70, 64, 69, 74, 72, 64, 47]
, [ 73, 76, 77, 80, 82, 76, 58, 37]
, [ 85, 85, 86, 86, 78, 62, 38, 20]
, [ 94, 94, 92, 78, 54, 33, 16, 255]
, [ 94, 90, 72, 51, 32, 19, 255, 255]
, [ 65, 53, 29, 18, 255, 255, 255, 255]
, [ 29, 22, 255, 255, 255, 255, 255, 0]]
, [[ 66, 67, 70, 69, 75, 73, 72, 63]
, [ 68, 70, 73, 74, 78, 80, 74, 53]
, [ 75, 87, 87, 83, 89, 86, 61, 33]
, [ 81, 89, 88, 98, 99, 77, 41, 18]
, [ 84, 94, 100, 100, 82, 49, 21, 255]
, [ 99, 101, 92, 75, 48, 25, 255, 255]
, [ 93, 77, 52, 32, 255, 255, 255, 255]
, [ 52, 40, 25, 255, 255, 255, 255, 255]
, [ 23, 16, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]]
, [[ 81, 83, 92, 101, 101, 83, 49, 19]
, [ 86, 96, 103, 103, 95, 64, 28, 255]
, [ 94, 103, 107, 98, 79, 41, 255, 255]
, [101, 103, 98, 79, 51, 28, 255, 255]
, [102, 97, 76, 49, 27, 255, 255, 255]
, [ 79, 62, 35, 21, 255, 255, 255, 255]
, [ 33, 23, 15, 255, 255, 255, 255, 255]
, [ 16, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]]
, [[106, 107, 109, 94, 58, 26, 15, 255]
, [110, 104, 90, 66, 37, 19, 255, 255]
, [106, 89, 61, 35, 22, 255, 255, 255]
, [ 76, 56, 34, 19, 255, 255, 255, 255]
, [ 40, 27, 18, 255, 255, 255, 255, 255]
, [ 17, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]]
, [[ 68, 51, 33, 19, 255, 255, 255, 255]
, [ 45, 34, 20, 255, 255, 255, 255, 255]
, [ 28, 18, 255, 255, 255, 255, 255, 255]
, [ 17, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 0, 0, 0, 0, 0, 0, 0]]
, [[255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 255]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 255, 0, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]]
, [[255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 255, 0]
, [255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 255, 0, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 255, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]]
, [[255, 255, 255, 255, 255, 255, 0, 0]
, [255, 255, 255, 255, 255, 0, 0, 0]
, [255, 255, 255, 255, 0, 0, 0, 0]
, [255, 255, 255, 0, 0, 0, 0, 0]
, [255, 255, 0, 0, 0, 0, 0, 0]
, [255, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]
, [ 0, 0, 0, 0, 0, 0, 0, 0]]])
实现方法：
[Z,Y,X]=im.shape
RN = np.array([3,4,4])
################Loading Area search
rad = 3
a,b,c = RN
x,y,z = np.ogrid[-c:Z-c,-b:Y-b,-a:X-a]
neighborMask = x*x + y*y + z*z<= rad*rad
noNodeMask = im > 0
mask = np.logical_and(neighborMask, noNodeMask)
imtemp = im.copy()
imtemp[mask] = -1
for i in range (X):
for j in range (Y):
for k in range (Z):
if imtemp[i,j,k]==-1:
if i in (0, X-1) or j in (0, Y-1) or k in (0, Z-1):
imtemp[i,j,k]=-2
elif imtemp[i+1,j,k] == 0 or imtemp[i-1,j,k] == 0 or imtemp[i,j+1,k] == 0 or imtemp[i,j-1,k] == 0 or imtemp[i,j,k+1] == 0 or imtemp[i,j,k-1] == 0:
imtemp[i,j,k]=-2
LA = np.argwhere(imtemp==-2)
上述示例代码生成的 LA 为：
In [90]:LA
Out[90]:
array([[4, 4, 0],
[4, 4, 6],
[4, 5, 5],
[4, 6, 4],
[4, 6, 5],
[4, 7, 3],
[5, 3, 5],
[5, 4, 4],
[5, 4, 5],
[5, 5, 3],
[5, 5, 4],
[5, 6, 2],
[5, 6, 3],
[6, 2, 4],
[6, 3, 3],
[6, 3, 4],
[6, 4, 2],
[6, 4, 3],
[6, 5, 1],
[6, 5, 2]])
Z 方向的切片（XY 平面实例）显示不同的未触及、屏蔽 (-1) 和目标 (-2) 节点：
由于您的循环仅使用直接 Numpy 索引，因此您可以使用 Numba 的 @njit 以 更高效的方式执行此操作 方式。
@njit
def compute_imtemp(imtemp, X, Y, Z):
for i in range (Z):
for j in range (Y-1):
for k in range (X-1):
if imtemp[i,j,k]==-1:
if i==(Z-1):
imtemp[i,j,k]=-2
elif imtemp[i+1,j,k] == 0 or imtemp[i-1,j,k] == 0 or imtemp[i,j+1,k] == 0 or imtemp[i,j-1,k] == 0 or imtemp[i,j,k+1] == 0 or imtemp[i,j,k-1] == 0:
imtemp[i,j,k]=-2
[...]
imtemp = im.copy()
imtemp[mask] = -1
compute_imtemp(imtemp, X, Y, Z)
LA = np.argwhere(imtemp==-2)
以下是我机器上的性能结果：
281 µs ± 1.43 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
776 ns ± 16.4 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
Numba 实施速度 362 倍。
请注意，由于编译，第一次调用 compute_imtemp 会很慢。克服这个问题的一种方法是在一个空的 Numpy 数组上调用 compute_imtemp。另一种方法是使用 Numba API 手动编译函数并显式向 Numba 提供类型。
问题陈述
你有一个大阵列中的“实心”形状。你从中切出一个球。您的目标是找到球内固体表面的指数。表面定义为与实体外部相邻且具有 6 点连通性的任何点。阵列的边缘也被认为是表面。
更快的循环解决方案
您已经计算出代表实体和球相交的掩码。您可以更优雅地计算掩码并将其转换为索引。我建议保持维度的顺序不变，而不是在不同的符号之间切换。例如，RN 的顺序会受到影响，并且您 运行 存在轴限制不匹配的风险。
RN = np.array([4, 4, 3])
rad = 3
im = ...
cutout = ((np.indices(im.shape) - RN.reshape(-1, 1, 1, 1))**2).sum(axis=0) <= rad**2
solid = im > 0
mask = solid & cutout
indices = np.argwhere(mask)
您也可以在不重塑 RN 的情况下获得切口，方法是
cutout = ((np.rollaxis(np.indices(im.shape, sparse=False), 0, 4) - RN)**2).sum(axis=-1) <= rad**2
计算索引的好处是您的循环不再需要很大。通过使用 argwhere，您基本上剥离了外部三个循环，只留下 if 语句进行循环。您还可以矢量化连接检查。这有一个很好的副作用，您可以为每个像素定义任意连接。
limit = np.array(im.shape) - 1 # Edge of `im`
connectivity = np.array([[ 1, 0, 0], # Add rows to determine connectivity
[-1, 0, 0],
[ 0, 1, 0],
[ 0, -1, 0],
[ 0, 0, 1],
[ 0, 0, -1]], dtype=indices.dtype)
index_mask = np.ones(len(indices), dtype=bool)
for n, ind in enumerate(indices):
if ind.all() and (ind < limit).all() and im[tuple((ind + connectivity).T)].all():
index_mask[n] = False
LA = indices[index_mask, :]
请注意 imtemp 根本没有意义。即使在您原来的循环中，您也可以直接操作 mask 。如果元素没有通过您的标准，您可以将元素设置为 False，而不是将它们设置为 -2。
我在这里做类似的事情。我们检查实际选择的每个索引，并确定它们中是否有任何一个在实体内部。这些索引从掩码中消除。然后根据掩码更新索引列表。
检查 ind.all() and (ind < limit).all() and im[tuple((ind + connectivity).T)].all() 是您在 or 条件下执行的操作的快捷方式，但相反（测试 non-surface 而不是表面）。
ind.all()检查 none 的索引是否为零：即不在 top/front/left 表面上。(ind < limit).all()检查 none 个索引等于相应的图像大小减一。im[tuple((ind + connectivity).T)].all()检查连接像素中的 none 是否为零。(ind + connectivity).T是我们连接到的六个点的(3, 6)数组（当前由(6, 3)connectivity数组在每个轴中定义为 +/-1）。当你把它变成一个元组时，它就变成了一个奇特的索引，就好像你做了类似im[x + connectivity[:, 0], y + connectivity[:, 1], z + connectivity[:, 2]]的事情一样。索引中的逗号只是使它成为一个元组。我展示的方式更适合任意数量的维度。
通过所有三个测试的像素都在实体内部，并被移除。你当然可以编写循环来检查另一种方式，但是你必须改变你的掩码：
index_mask = np.zeros(len(indices), dtype=bool)
for n, ind in enumerate(indices):
if (ind == 0).any() or (ind == limit).any() or (im[tuple((ind + connectivity).T)] == 0).any():
index_mask[n] = True
LA = indices[index_mask, :]
手动循环无论如何都不理想。但是，它向您展示了如何缩短循环（可能缩短几个数量级），以及如何使用矢量化和广播定义任意连接，而不会陷入 hard-coding 它的困境。
完全矢量化的解决方案
上面的循环可以使用广播的魔力完全矢量化。我们可以批量添加 connectivity 并批量过滤结果，而不是遍历 indices 中的每一行。诀窍是添加足够的维度，将所有 connectivity 添加到 indices.
您仍然希望忽略边缘的像素：
edges = (indices == 0).any(axis=-1) | (indices == limit).any(axis=-1)
conn_index = indices[~edges, None, :] + connectivity[None, ...]
index_mask = np.empty(len(indices), dtype=bool)
index_mask[edges] = True
index_mask[~edges] = (im[tuple(conn_index.T)] == 0).any(axis=0)
LA = indices[index_mask, :]
我希望使用 numba 编译的正确编写的循环会比此解决方案快得多，因为它将避免流水线操作的大部分开销。它不需要大的临时缓冲区或特殊处理。
TL;DR
# Parameters
RN = np.array([4, 4, 3])
rad = 3
im = ...
# Find subset of interest
cutout = ((np.indices(im.shape) - RN.reshape(-1, 1, 1, 1))**2).sum(axis=0) <= rad**2
solid = im > 0
# Convert mask to indices
indices = np.argwhere(solid & cutout)
# Find image edges among indices
edges = (indices == 0).any(axis=-1) | (indices == limit).any(axis=-1)
# Connectivity elements for non-edge pixels
conn_index = indices[~edges, None, :] + connectivity[None, ...]
# Mask the valid surface pixels
index_mask = np.empty(len(indices), dtype=bool)
index_mask[edges] = True
index_mask[~edges] = (im[tuple(conn_index.T)] == 0).any(axis=0)
# Final result
LA = indices[index_mask, :]