使用 NumPy 将 3D 函数映射到网格网格的有效方法

Question

我有一组标量 3D 函数的数据值，它们排列为形状为 (n,3) 的数组中的输入 x,y,z 和形状为 (n,) 的数组中的函数值 f(x,y,z)。

编辑：例如，考虑以下简单函数

data = np.array([np.arange(n)]*3).T
F = np.linalg.norm(data,axis=1)**2

我想将此函数与球核进行卷积，以执行 3D 平滑。我发现执行此操作的最简单方法是将函数值映射到 3D 空间网格中，然后使用我想要的内核应用 3D 卷积。

这很好用，但是将 3D 函数映射到 3D 网格的部分非常慢，因为我没有找到仅使用 NumPy 的方法。下面的代码是我的实际实现，其中data 是(n,3) 数组，其中包含函数被评估的3D 位置，F 是(n,) 数组，其中包含函数的相应值，M 是包含 3D 空间网格的 (N,N,N) 数组。

step = 0.1

# Create meshgrid
xmin = data[:,0].min()
xmax = data[:,0].max()
ymin = data[:,1].min()
ymax = data[:,1].max()
zmin = data[:,2].min()
zmax = data[:,2].max()

x = np.linspace(xmin,xmax,int((xmax-xmin)/step)+1)
y = np.linspace(ymin,ymax,int((ymax-ymin)/step)+1)
z = np.linspace(zmin,zmax,int((zmax-zmin)/step)+1)


# Build image
M = np.zeros((len(x),len(y),len(z)))

for l in range(len(data)):
    for i in range(len(x)-1):
        if x[i] < data[l,0] < x[i+1]:
            for j in range(len(y)-1):
                if y[j] < data[l,1] < y[j+1]:
                    for k in range(len(z)-1):
                        if z[k] < data[l,2] < z[k+1]:
                            M[i,j,k] = F[l]

有没有更有效的方法来用 3D 函数的值填充 3D 空间网格？

Answer 1

对于每一项数据，您都在扫描长方体的像素以检查它是否在里面。有一个选项可以跳过此扫描。您可以自己计算这些像素的对应索引，例如：

data = np.array([[1, 2, 3], #14 (corner1)
                 [4, 5, 6], #77 (corner2)
                 [2.5, 3.5, 4.5], #38.75 (duplicated pixel)
                 [2.9, 3.9, 4.9], #47.63 (duplicated pixel)
                 [1.5, 2, 3]]) #15.25 (one step up from [1, 2, 3])
step = 0.5                            
data_idx = ((data - data.min(axis=0))//step).astype(int)
M = np.zeros(np.max(data_idx, axis=0) + 1)
x, y, z = data_idx.T
M[x, y, z] = F

请注意，只有一个重复像素值被映射到 M。

Answer 2

您只需将F[:, 3]（仅 f(x, y, z)）重塑为网格。没有样本数据很难更精确：

如果数据没有排序，则需要排序：

F_sorted = F[np.lexsort((F[:,0], F[:,1], F[:,2]))]  # sort by x, then y, then z

只选择 f(x, y, z)

F_values = F_sorted[:, 3]

最后，将数据重塑为网格：

M = F_sorted.reshape(N, N, N)

Answer 3

这种方法比原来的方法更快（大约快 20 倍）：

step = 0.1
mins = np.min(data, axis=0)
maxs = np.max(data, axis=0)
ranges = np.floor((maxs - mins) / step + 1).astype(int)
indx = np.zeros(data.shape,dtype=int)
for i in range(3):
    x = np.linspace(mins[i], maxs[i], ranges[i])
    indx[:,i] = np.argmax(data[:,i,np.newaxis] <= (x[np.newaxis,:]), axis=1) -1
    
M = np.zeros(ranges)
M[indx[:,0],indx[:,1],indx[:,2]] = F

第一部分设置所需的网格变量。 argmax 函数提供了一种简单（快速）的方法来查找广播数组的第一个真值。这将为每个函数值生成一组 x、y 和 z 方向的索引。

结果数组M 与原始代码生成的数组不同，因为原始代码会丢失数据。 y[j] < data[l,1] < y[j+1] 的逻辑，其中y 是使用linspace 生成的向量，意味着将丢失每个方向的最小值和最大值（data[l,1] 可能等于y[j] 或y[j+1]！）。使用包含两个值的数据集运行它，每个值都有自己的坐标，M 数组将全为零。