在 3D numpy 数组中绘制和/或获取 2D 平面的索引

编程问答 2022-08-26

问题描述

我有一个 3D numpy 数组，例如(200,200,200) 形状。

我也有一个二维平面，我可以用下面的形式表达 ax + by + cz + d = 0 或 (x,y,z) 点和以三个 x,z 轴表示的方向向量。

我想在 3D numpy 数组上“绘制”这个平面，方法是将位于该平面上的所有点的值更改为一个常数值。

像这样，但只有一种颜色：

enter image description here

对此（我认为）最通用的解决方案是获取沿这条线的所有点的整数索引。然后我可以使用它来索引 numpy 数组。然而，如果一个库提供了直接绘制平面的能力，而不提供索引，那么这会很好。

我在这个站点上看到过一些表面上与我的问题相似的问题，例如使用 x,z 查找沿 scipy.ndimage.map_coordinates 向量的值的方法，但当然我正在处理一个平面不是一条线（这也只是返回值，而不是索引）。

我考虑过的另一种方法，但似乎很难（而且可能很慢）类似于 this question，其中回复显示了如何在 3D 空间中绘制 2D 三角形。我可以改为在 3D 空间中绘制一个正方形，但这表明填充这个正方形并非微不足道，而且如果至少有一个轴平行于 x、y 或 z（或“角”），则正方形只会填充整个平面将保持空置）。

有人知道我如何实现这一目标吗？

解决方法

如果平面接近水平，您可以给 x 和 y 赋值并计算 z：

a,b,c,d = 1,2,3,-600                     # Plane parameters
v = np.arange(200)
x,y = np.meshgrid(v,v)                       # All xy combinations
z = np.rint((a*x + b*y + d) / -c).astype(int)  # Rounded
plane_voxels = np.dstack([x,y,z]).reshape(-1,3)
print(plane_voxels)

结果：

[[  0   0 200]
 [  1   0 200]
 [  2   0 199]
 ...
 [197 199   2]
 [198 199   1]
 [199 199   1]]

对于一般情况，您需要找到可变性较小的维度，这将是计算变量，并为其他两个提供值：

a,-600
v = np.arange(200)
dim = np.argmax(np.abs((a,c)))
if dim == 0:
    y,z = np.meshgrid(v,v)
    x = np.rint((b*y + c*z + d) / -a).astype(int)
elif dim == 1:
    x,v)
    y = np.rint((a*x + c*z + d) / -b).astype(int)
elif dim == 2:
    x,v)
    z = np.rint((a*x + b*y + d) / -c).astype(int)
plane_voxels = np.dstack([x,3)
,

这是一个想法。首先，假设您有一个 (3,3)。然后定义平面函数等于 0。用 meshgrid 创建坐标并在平面函数中使用它。最后，如果您想要一些容差，请使用 isclose 到 0 作为数组的掩码来更改值

n = 3 # change to 200 for your real case or arr.shape[0]

# just for the example
arr =  np.ones([n]*3)

# plane function = 0
f = lambda x,z: 2*x + 3*y**2 + 4*z - 3

# create the mask
mask = np.isclose(f(*np.meshgrid(*[np.arange(n)]*3)),# because your plane function is equal to 0
                  atol=1) # this is if you want some tolerance to catch nearby points

# change the value to whatever
arr[mask] = 5

print(arr)
[[[1. 5. 1.]
  [5. 1. 1.]
  [5. 1. 1.]]

 [[5. 1. 1.]
  [1. 1. 1.]
  [1. 1. 1.]]

 [[1. 1. 1.]
  [1. 1. 1.]
  [1. 1. 1.]]]
numpypythonscipyscipy