使用python对3D中的2D表面点进行Delaunay三角剖分？

Question

我有一组 3D 点。这些点以恒定水平采样 (z=0,1,...,7)。图片应该清楚：

这些点位于形状为(N, 3) 的numpy ndarray 中，称为X。上面的图是使用以下方法创建的：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

X = load('points.npy')
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(X[:,0], X[:,1], X[:,2])
ax.scatter(X[:,0], X[:,1], X[:,2])
plt.draw()

我想只对这个对象的表面进行三角测量，并绘制表面。但是，我不想要这个对象的凸包，因为这会丢失我希望能够检查的细微形状信息。

我试过ax.plot_trisurf(X[:,0], X[:,1], X[:,2])，但这会导致以下混乱：

有什么帮助吗？

示例数据

这是一个用于生成代表问题的 3D 数据的 sn-p：

import numpy as np
X = []
for i in range(8):
    t = np.linspace(0,2*np.pi,np.random.randint(30,50))
    for j in range(t.shape[0]):
        # random circular objects...
        X.append([
            (-0.05*(i-3.5)**2+1)*np.cos(t[j])+0.1*np.random.rand()-0.05,
            (-0.05*(i-3.5)**2+1)*np.sin(t[j])+0.1*np.random.rand()-0.05,
            i
        ])
X = np.array(X)

来自原始图像的示例数据

这是原始数据的粘贴箱：

http://pastebin.com/YBZhJcsV

这是沿常数 z 的切片：

Answer 1

我意识到您在问题中提到您不想使用凸包，因为您可能会丢失一些形状信息。我有一个简单的解决方案，它非常适合您的“抖动球形”示例数据，尽管它确实使用了scipy.spatial.ConvexHull。我想我还是会在这里分享它，以防万一它对其他人有用：

from matplotlib.tri import triangulation
from scipy.spatial import ConvexHull

# compute the convex hull of the points
cvx = ConvexHull(X)

x, y, z = X.T

# cvx.simplices contains an (nfacets, 3) array specifying the indices of
# the vertices for each simplical facet
tri = Triangulation(x, y, triangles=cvx.simplices)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.hold(True)
ax.plot_trisurf(tri, z)
ax.plot_wireframe(x, y, z, color='r')
ax.scatter(x, y, z, color='r')

plt.draw()

在这种情况下它做得很好，因为您的示例数据最终位于或多或少的凸面上。也许您可以制作一些更具挑战性的示例数据？环面是一个很好的测试用例，而凸包方法显然会失败。

从点云映射任意 3D 表面是一个真正棘手的问题。这是一个related question，其中包含一些可能有用的链接。

Answer 2

更新 3

这是我在更新 2 中描述的具体示例。如果您没有用于可视化的 mayavi，我建议使用 edm install mayavi pyqt matplotlib 安装它 via edm。

在 3D 中堆叠的玩具 2D 轮廓

轮廓 -> 3D 表面

生成图形的代码

from matplotlib import path as mpath
from mayavi import mlab
import numpy as np


def make_star(amplitude=1.0, rotation=0.0):
    """ Make a star shape
    """
    t = np.linspace(0, 2*np.pi, 6) + rotation
    star = np.zeros((12, 2))
    star[::2] = np.c_[np.cos(t), np.sin(t)]
    star[1::2] = 0.5*np.c_[np.cos(t + np.pi / 5), np.sin(t + np.pi / 5)]
    return amplitude * star

def make_stars(n_stars=51, z_diff=0.05):
    """ Make `2*n_stars-1` stars stacked in 3D
    """
    amps = np.linspace(0.25, 1, n_stars)
    amps = np.r_[amps, amps[:-1][::-1]]
    rots = np.linspace(0, 2*np.pi, len(amps))
    zamps = np.linspace
    stars = []
    for i, (amp, rot) in enumerate(zip(amps, rots)):
        star = make_star(amplitude=amp, rotation=rot)
        height = i*z_diff
        z = np.full(len(star), height)
        star3d = np.c_[star, z]
        stars.append(star3d)
    return stars

def polygon_to_boolean(points, xvals, yvals):
    """ Convert `points` to a boolean indicator mask
    over the specified domain
    """
    x, y = np.meshgrid(xvals, yvals)
    xy = np.c_[x.flatten(), y.flatten()]
    mask = mpath.Path(points).contains_points(xy).reshape(x.shape)
    return x, y, mask

def plot_contours(stars):
    """ Plot a list of stars in 3D
    """
    n = len(stars)

    for i, star in enumerate(stars):
        x, y, z = star.T
        mlab.plot3d(*star.T)
        #ax.plot3D(x, y, z, '-o', c=(0, 1-i/n, i/n))
        #ax.set_xlim(-1, 1)
        #ax.set_ylim(-1, 1)
    mlab.show()



if __name__ == '__main__':

    # Make and plot the 2D contours
    stars3d = make_stars()
    plot_contours(stars3d)

    xvals = np.linspace(-1, 1, 101)
    yvals = np.linspace(-1, 1, 101)

    volume = np.dstack([
        polygon_to_boolean(star[:,:2], xvals, yvals)[-1]
        for star in stars3d
    ]).astype(float)

    mlab.contour3d(volume, contours=[0.5])
    mlab.show()

更新 2

我现在这样做如下：

1. 我利用每个 z 切片中的路径是封闭且简单的这一事实，并使用matplotlib.path 来确定轮廓内部和外部的点。使用这个想法，我将每个切片中的轮廓转换为布尔值图像，然后组合成布尔值体积。
2. 接下来，我使用skimage's marching_cubes 方法获取表面的三角剖分以进行可视化。

这是该方法的一个示例。我认为数据略有不同，但您肯定可以看到结果更清晰，并且可以处理断开或有孔的表面。

原答案

好的，这是我想出的解决方案。我认为这在很大程度上取决于我的数据大致呈球形并在 z 中均匀采样。其他一些 cmets 提供了有关更强大的解决方案的更多信息。由于我的数据大致是球形的，因此我从数据点的球坐标变换中对方位角和天顶角进行三角测量。

import numpy as np
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri

X = np.load('./mydatars.npy')
# My data points are strictly positive. This doesn't work if I don't center about the origin.
X -= X.mean(axis=0)

rad = np.linalg.norm(X, axis=1)
zen = np.arccos(X[:,-1] / rad)
azi = np.arctan2(X[:,1], X[:,0])

tris = mtri.Triangulation(zen, azi)

fig = plt.figure()
ax  = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(X[:,0], X[:,1], X[:,2], triangles=tris.triangles, cmap=plt.cm.bone)
plt.show()

使用上面 pastebin 中的示例数据，得到：

---