python – numpy数组上的边框/边操作

2023年12月13日 284次阅读

假设我有一个非零值的非numpy数组,并且“background”= 0.作为一个例子,我将采用随机值的范围：

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)功能吗？

编辑

我搞砸了一件事：我想用非零和非边界邻居的平均值替换所有边界体素.为此,我试图从ali_m的代码中清理邻居(2D情况)：

#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会变回空？

此外,纠正我,如果我错了,但不是tree.query_ball_point与半径1只加上6个下一个邻居(欧几里德距离1)？我应该将sqrt(3)(3D case)设置为半径以获得26邻域吗？

最佳答案 我认为最好先从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

要找到边缘像素,您可以对蒙版执行二进制侵蚀,然后使用蒙版对结果进行异或

# 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()

《python – numpy数组上的边框/边操作》

这种方法也将推广到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]

# ...
点赞