<p><strong>问题陈述</strong></p>
<p>在一个大数组中有一个“实体”形状。你从上面刻出一个球。你的目标是找到球内固体表面的指数。曲面定义为具有6点连接的实体外部相邻的任何点。阵列的边也被视为曲面</p>
<p><strong>更快的环路解决方案</strong></p>
<p>您已经计算了表示实体和球的交点的遮罩。您可以更优雅地计算掩码,并将其转换为索引。我建议保持维度的顺序不变,而不是在不同的符号之间切换。例如,<code>RN</code>的顺序会受到影响,并且您会面临与轴限制不匹配的风险</p>
<pre><code>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)
</code></pre>
<p>通过执行以下操作,也可以在不重塑<code>RN</code>的情况下获得剪切</p>
<pre><code>cutout = ((np.rollaxis(np.indices(im.shape, sparse=False), 0, 4) - RN)**2).sum(axis=-1) <= rad**2
</code></pre>
<p>计算索引的好处在于,循环不再需要太大。通过使用<code>argwhere</code>,基本上去掉了外部的三个循环,只剩下<code>if</code>语句进行循环。您还可以将连接检查矢量化。这有一个很好的副作用,你可以为每个像素定义任意的连接</p>
<pre><code>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, :]
</code></pre>
<p>请注意<code>imtemp</code>实际上毫无意义。即使在原始循环中,也可以直接操作<code>mask</code>。当元素通过您的标准时,您可以将元素设置为<code>False</code>,而不是将它们设置为<code>-2</code>,如果它们没有通过标准,则可以将元素设置为<code>False</code></p>
<p>我在这里做类似的事情。我们检查实际选择的每个索引,并确定其中是否有任何索引位于实体内部。这些指数从遮罩中消除。然后根据掩码更新索引列表</p>
<p>检查<code>ind.all() and (ind < limit).all() and im[tuple((ind + connectivity).T)].all()</code>是对<code>or</code>条件执行操作的快捷方式,但相反(测试非曲面而非曲面)</p>
<ul>
<li><code>ind.all()</code>检查所有索引是否为零:即,不在顶部/前部/左侧曲面上</李>
<li><code>(ind < limit).all()</code>检查所有索引都不等于相应的图像大小减1</李>
<li><code>im[tuple((ind + connectivity).T)].all()</code>检查连接的像素是否为零<code>(ind + connectivity).T</code>是我们连接到的六个点的<code>(3, 6)</code>数组(当前由<code>(6, 3)</code>{<cd17>}数组在每个轴上定义为+/-1)。当你把它变成一个元组时,它就变成了一个奇特的索引,就像你做了类似<code>im[x + connectivity[:, 0], y + connectivity[:, 1], z + connectivity[:, 2]]</code>的事情一样。索引中的逗号只是将其组成一个元组。我展示的方式更适合任意数量的维度</李>
</ul>
<p>通过所有三个测试的像素位于实体内部,并被移除。当然,您可以通过另一种方式写入循环以进行检查,但随后您必须更改掩码:</p>
<pre><code>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, :]
</code></pre>
<p>手动循环无论如何都不理想。然而,它向您展示了如何缩短循环(可能是几个数量级),以及如何使用矢量化和广播定义任意连接,而不会陷入硬编码的困境</p>
<p><strong>完全矢量化解决方案</strong></p>
<p>上面的循环可以使用广播的魔力完全矢量化。我们可以将<code>connectivity</code>批量添加到<code>connectivity</code>中并批量过滤结果,而不是在<code>indices</code>中的每一行上循环。诀窍是添加足够的维度,以便将所有<code>connectivity</code>添加到<code>indices</code>的每个</em>元素中</p>
<p>您仍然希望忽略边缘处的像素:</p>
<pre><code>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, :]
</code></pre>
<p>我期望使用numba编译的正确编写的循环将比这个解决方案快得多,因为它将避免许多操作管道化的开销。它不需要大型临时缓冲或特殊处理</p>
<p><strong>TL;DR</strong></p>
<pre><code># 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, :]
</code></pre>