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

Question

我有一组数据值，用于标量3D函数，按照输入 x,y,z 的顺序排列在形如 (n,3) 的数组中，并且函数值 f(x,y,z) 排列在形如 (n,) 的数组中。

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

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

我希望使用球形卷积核来对此函数进行卷积，以执行3D平滑。我发现最简单的方法是将函数值映射到一个3D空间网格中，然后使用我想要的卷积核进行3D卷积。

这种方法可以正常工作，但是将3D函数映射到3D网格的部分非常缓慢，因为我没有找到仅使用NumPy完成它的方法。下面的代码是我的实际实现，其中data是包含函数在其中评估的3D位置的(n,3)数组，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空间网格中？

- Toool

---

我对问题不是很清楚，没有一些数据样本很难进行一些测试。但是你尝试过使用np.meshgrid来创建3D网格吗？ - Ulises Bussi

“3D功能”在哪里？所有这些if语句将使执行“整个数组”操作变得困难。它们本质上是标量。 - hpaulj

我已经更新了我的问题，希望现在清楚了。所谓的“3D函数”，是指以三维向量为输入的标量函数。 - Toool

请注意，由于条件：“具有多个元素的数组的真值是模棱两可的。请使用a.any()或a.all()”，当前代码无法正常工作。 - Jérôme Richard

你能解释一下循环中那些 if 条件的动机吗？ - amzon-ex

@JérômeRichard 这段代码在我的电脑上执行没有问题（刚刚测试过）。 - Toool

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。

- mathfux

---

你能详细解释一下它正在做什么吗？我测试过了，但是无法得到与我上面的代码相同的结果。 - Toool

我认为与原始代码的区别在于比较。该方法使用截断以高效地设置值。这相当于检查 x[i] <= data[l,0] < x[i+1]（其他维度也是如此）。原始代码测试 x[i] < data[l,0] < x[i+1]。请注意，下限缺少的地方有一个“=”。我不确定这是否是有意的，但我发现原始代码排除了一些值，因为输入只包含整数。因此，在计算原始代码后，M 总是为零（即循环什么也不做）... - Jérôme Richard

@user14900935的回答基本上是一样的，但帮助我澄清了你的解决方案，谢谢你们两个！ - Toool

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)

- dankal444

---

你的回答提供了一种方法，可以将形状为(N,)的函数值转换为形状为(N,3)的3D向量网格上的值，并将其转换为排序后的3D网格(n,n,n)，其中n=N**(1/3)。虽然我的问题有点不同，但这对其他人可能会有用。 - Toool

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与原始代码生成的数组不同，因为原始代码会丢失数据。 使用使用linspace生成的向量y 的逻辑'y [j]

- Mat

---

谢谢，您的解决方案实际上与 @mathfux 的相同，但更清晰。由于您们两个都提供了有效的解决方案，我将给您奖励并接受另一个方案。 - Toool

谢谢。我猜他们的代码会更快一些，因为它去掉了循环和“linspace”命令，但很高兴能帮到你。 - Mat