所以你正在寻找这样的东西?
我已采用此https://matplotlib.org/stable/gallery/mplot3d/voxels_rgb.html 来获取您的单工。
import itertools
import functools
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (18,9)
points = np.array([[1,2,2],[1,3,4],[4,1,4],[5,3,2]])
center = np.mean(points, axis=0)
def surface_normal_form(a,b,c):
v = b-a
w = c-b
n = np.cross(v,w)
#normal needs to point out
if (center-a)@n > 0:
n *= -1
return a, n
def midpoints(x):
sl = ()
for i in range(x.ndim):
x = (x[sl + np.index_exp[:-1]] + x[sl + np.index_exp[1:]]) / 2.0
sl += np.index_exp[:]
return x
x, y, z = (np.indices((60, 60, 60))-np.array([20,25,25]).reshape(-1,1,1,1))/8
mx = midpoints(x)
my = midpoints(y)
mz = midpoints(z)
conditions = []
for p1,p2,p3 in itertools.combinations(points, 3):
a, n = surface_normal_form(p1,p2,p3)
conditions.append((mx-a[0])*n[0]+(my-a[1])*n[1]+(mz-a[2])*n[2] <= 0)
simplex = conditions[0] & conditions[1] & conditions[2] & conditions[3]
ax = plt.figure().add_subplot(projection='3d')
ax.voxels(x, y, z, simplex, linewidth=0.5)
ax.set(xlabel='x', ylabel='y', zlabel='z')
ax.set_xlim(1.0,5.0)
ax.set_ylim(1.0,3.0)
ax.set_zlim(2.0,4.0)
你能解释一下你在 surface_normal_form 函数中做了什么吗?
当然,在 3d 中,您可以通过平面的一点和与平面正交的向量来描述平面。这特别有用,因为它很容易判断一个点是在平面的一侧还是另一侧。然后,如果您采用包含单纯形侧面的平面,则通过将一个体素立方体紧密贴合在一起并为每个平面移除其错误一侧的体素,您将获得单纯形体素。这就是我正在做的。
另外,能否请您解释一下中点函数?
首先要声明一点,这不是我自己写的。正如我所说,它来自 matplotlib 示例。但是如果你想计算一维数组的中点,你可以这样做。
(x[:-1]+x[1:])/2
首先x[:-1] 为您提供除最后一个之外的所有值,x[1:] 为您提供除第一个之外的所有值,因此彼此相邻的值相加并除以二,即您得到中点。请注意,如果原始数组大于 1d(在我们的示例中为 3 维),则它采用(在我们的示例中为 2d)子数组的“中点”。在第二步中,我们正在做的是(x[:,:-1]+x[:,1:])/2。
由于[:] 为您提供了所有值,midpoint 函数为每个维度执行此操作。