python - numpy 数组的边界/边缘操作

Question

假设我有一个非零值和 "background"= 0 的 3D numpy 数组。作为一个例子，我将采用一个随机值的球体:

array = np.random.randint(1, 5, size = (100,100,100))
z,y,x = np.ogrid[-50:50, -50:50, -50:50]
mask = x**2 + y**2 + z**2<= 20**2
array[np.invert(mask)] = 0

首先，我想找到“边界体素”(所有在其 3x3x3 邻域内具有零的非零值)。其次，我想用它们的非零邻居的平均值替换所有边界体素。到目前为止，我尝试通过以下方式使用 scipy 的通用过滤器:

应用于每个元素的函数:

def borderCheck(values):
    #check if the footprint center is on a nonzero value
    if values[13] != 0:
        #replace border voxels with the mean of nonzero neighbours
        if 0 in values:
            return np.sum(values)/np.count_nonzero(values)
        else:
            return values[13]
    else:
        return 0

通用过滤器:

from scipy import ndimage
result = ndimage.generic_filter(array, borderCheck, footprint = np.ones((3,3,3)))

这是处理这个问题的正确方法吗？我觉得我正在尝试在这里重新发明轮子，并且必须有一种更短、更好的方法来实现结果。我可以使用其他合适的(numpy、scipy)函数吗？

编辑

我搞砸了一件事:我想用它们的非零 AND 非边界 邻居的平均值替换所有边界体素。为此，我尝试从 ali_m 的代码(2D 案例)中清理 neighbors:

#for each neighbour voxel, check whether it also appears in the border/edges
non_border_neighbours = []
for each in neighbours:
    non_border_neighbours.append([i for i in each if nonzero_idx[i] not in edge_idx])

现在我不明白为什么 non_border_neighbours 返回时是空的？

此外，如果我错了，请纠正我，但半径为 1 的 tree.query_ball_point 不会仅寻址 6 个下一个邻居(欧氏距离 1)吗？我是否应该将 sqrt(3)(3D 案例)设置为半径以获得 26 邻域？

Answer 1

我认为最好先从 2D 案例开始，因为它可以更容易地可视化:

import numpy as np
from matplotlib import pyplot as plt

A = np.random.randint(1, 5, size=(100, 100)).astype(np.double)
y, x = np.ogrid[-50:50, -50:50]
mask = x**2 + y**2 <= 30**2
A[~mask] = 0

要找到边缘像素，您可以对蒙版执行二元腐 eclipse ，然后将结果与蒙版异或

# rank 2 structure with full connectivity
struct = ndimage.generate_binary_structure(2, 2)
erode = ndimage.binary_erosion(mask, struct)
edges = mask ^ erode

找到每个边缘像素最近的非零邻居的一种方法是使用 scipy.spatial.cKDTree :

from scipy.spatial import cKDTree

# the indices of the non-zero locations and their corresponding values
nonzero_idx = np.vstack(np.where(mask)).T
nonzero_vals = A[mask]

# build a k-D tree
tree = cKDTree(nonzero_idx)

# use it to find the indices of all non-zero values that are at most 1 pixel
# away from each edge pixel
edge_idx = np.vstack(np.where(edges)).T
neighbours = tree.query_ball_point(edge_idx, r=1, p=np.inf)

# take the average value for each set of neighbours
new_vals = np.hstack(np.mean(nonzero_vals[n]) for n in neighbours)

# use these to replace the values of the edge pixels
A_new = A.astype(np.double, copy=True)
A_new[edges] = new_vals

一些可视化:

fig, ax = plt.subplots(1, 3, figsize=(10, 4), sharex=True, sharey=True)
norm = plt.Normalize(0, A.max())
ax[0].imshow(A, norm=norm)
ax[0].set_title('Original', fontsize='x-large')
ax[1].imshow(edges)
ax[1].set_title('Edges', fontsize='x-large')
ax[2].imshow(A_new, norm=norm)
ax[2].set_title('Averaged', fontsize='x-large')
for aa in ax:
    aa.set_axis_off()
ax[0].set_xlim(20, 50)
ax[0].set_ylim(50, 80)
fig.tight_layout()
plt.show()

这种方法也将推广到 3D 情况:

B = np.random.randint(1, 5, size=(100, 100, 100)).astype(np.double)
z, y, x = np.ogrid[-50:50, -50:50, -50:50]
mask = x**2 + y**2 + z**2 <= 20**2
B[~mask] = 0

struct = ndimage.generate_binary_structure(3, 3)
erode = ndimage.binary_erosion(mask, struct)
edges = mask ^ erode

nonzero_idx = np.vstack(np.where(mask)).T
nonzero_vals = B[mask]

tree = cKDTree(nonzero_idx)

edge_idx = np.vstack(np.where(edges)).T
neighbours = tree.query_ball_point(edge_idx, r=1, p=np.inf)

new_vals = np.hstack(np.mean(nonzero_vals[n]) for n in neighbours)

B_new = B.astype(np.double, copy=True)
B_new[edges] = new_vals

针对您的版本进行测试:

def borderCheck(values):
    #check if the footprint center is on a nonzero value
    if values[13] != 0:
        #replace border voxels with the mean of nonzero neighbours
        if 0 in values:
            return np.sum(values)/np.count_nonzero(values)
        else:
            return values[13]
    else:
        return 0

result = ndimage.generic_filter(B, borderCheck, footprint=np.ones((3, 3, 3)))

print(np.allclose(B_new, result))
# True

我敢肯定这不是最有效的方法，但它仍然比使用 generic_filter 快得多。

---

更新

通过减少被视为边缘像素/体素候选邻域的点的数量，可以进一步提高性能:

# ...

# the edge pixels/voxels plus their immediate non-zero neighbours
erode2 = ndimage.binary_erosion(erode, struct)
candidate_neighbours = mask ^ erode2

nonzero_idx = np.vstack(np.where(candidate_neighbours)).T
nonzero_vals = B[candidate_neighbours]

# ...