你在找这样的东西？

我已经采用了这个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:]给出了除了第一个值之外的所有值，所以相邻的值被求和并除以2。e.你得到中点。请注意，如果原始数组超过1d（在我们的例子中是3维），则它将获取（在我们的例子中是2d）子数组的“中点”。在第二步中，我们要做的是(x[:,:-1]+x[:,1:])/2。因为[:]给了你所有的值，midpoint函数会为每个维度做这些。