如何找到没有孔的 3D 对象的整个外壳？

Question

我有一个由 8 个节点组成的云，它们形成一个棱镜。我需要哪个节点属于棱镜哪一侧的信息。棱镜的一侧可能是平面，但情况并非总是如此。

所以我想我可以用 ConvexHull 搜索三角形，然后用 this answer from 找到共面三角形。这对于对称棱镜（如长方体）非常有效，但在我的情况下，物体根本没有矩形边缘。将代码转移到我的问题时，不幸的是，表面上总是有“洞”。

到目前为止，这是我的代码。

import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import networkx as nx
from sympy import Plane, Point3D
import mpl_toolkits.mplot3d as a3
import matplotlib.colors as colors
import scipy as sp

nodes = np.array([[5.05563104322024e+04, 9.86881840921214e+04, -1.86894465727364e+03], [5.05703988843260e+04, 9.86813866643775e+04, -2.17336823588979e+03], [5.02542761707438e+04, 9.87087037212873e+04, -1.87234751337548e+03], [5.05535885714918e+04, 9.90000078125000e+04, -1.88856455463216e+03], [5.02534400234038e+04, 9.86956918140383e+04, -2.16834889960677e+03], [5.02542854985772e+04, 9.90000078125000e+04, -1.86873609902935e+03], [5.05494853841027e+04, 9.90000078125000e+04, -2.20374533059710e+03], [5.02533554687500e+04, 9.90000078125000e+04, -2.19881323242188e+03]])

fig = plt.figure()
ax = Axes3D(fig)
ax.dist=10
ax.azim=30
ax.elev=30
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

verts = zip(*nodes)

ax.plot(verts[0], verts[1], verts[2], 'bo')
hull = ConvexHull(nodes)
faces = hull.simplices

triangles = []
for i in hull.simplices:
    i = np.append(i, i[0])
    plt.plot(nodes[i,0], nodes[i,1], nodes[i,2])

for j in faces:
    tri = [(nodes[j[0], 0], nodes[j[0], 1], nodes[j[0], 2]), (nodes[j[1], 0], nodes[j[1], 1], nodes[j[1], 2]), (nodes[j[2], 0], nodes[j[2], 1], nodes[j[2], 2])]
    triangles.append(tri)

def simplify(triangles):
    G = nx.Graph()
    G.add_nodes_from(range(len(triangles)))
    for ii, a in enumerate(triangles):
        for jj, b in enumerate(triangles):
            if (ii < jj): 
                if is_adjacent(a, b):
                    if is_coplanar(a, b, np.pi / 180.):
                        G.add_edge(ii,jj)
    components = list(nx.connected_components(G))
    simplified = [set(flatten(triangles[index] for index in component)) for component in components]
    reordered = [reorder(face) for face in simplified]
    return reordered

def is_adjacent(a, b):
    return len(set(a) & set(b))

def is_coplanar(a, b, tolerance_in_radians=0):
    a1, a2, a3 = a
    b1, b2, b3 = b
    plane_a = Plane(Point3D(a1), Point3D(a2), Point3D(a3))
    plane_b = Plane(Point3D(b1), Point3D(b2), Point3D(b3))
    if not tolerance_in_radians: # only accept exact results
        return plane_a.is_coplanar(plane_b)
    else:
        angle = plane_a.angle_between(plane_b).evalf()
        angle %= np.pi # make sure that angle is between 0 and np.pi
        return (angle - tolerance_in_radians <= 0.) or \
            ((np.pi - angle) - tolerance_in_radians <= 0.)

flatten = lambda l: [item for sublist in l for item in sublist]

def reorder(vertices):
    if len(vertices) <= 3: 
        return vertices
    else:
        reordered = [vertices.pop()]
        vertices = list(vertices)
        while len(vertices) > 1:
            idx = np.argmin(get_distance(reordered[-1], vertices))
            v = vertices.pop(idx)
            reordered.append(v)
        reordered += vertices
        return reordered

def get_distance(v1, v2):
    v2 = np.array(list(v2))
    difference = v2 - v1
    ssd = np.sum(difference**2, axis=1)
    return np.sqrt(ssd)

new_faces = simplify(triangles)

for sq in new_faces:
    f = a3.art3d.Poly3DCollection([sq])
    f.set_color(colors.rgb2hex(sp.rand(3)))
    f.set_edgecolor('k')
    f.set_alpha(0.1)
    ax.add_collection3d(f)

plt.show()

---

正如您在下面的屏幕截图中所见，我的棱镜顶部有一个孔。你知道我该如何解决这个问题并获得棱镜外壳的所有三角形或四边形吗？

谢谢！

---

截图：

编辑： 我对我的代码进行了一些更改（见下文），所以现在所有三角形/四边形都是单独绘制的，并且所有面都找到了。会不会是绘图错误？

# verts = zip(*nodes)
verts = nodes.T

hull = ConvexHull(nodes, incremental=True)
faces = hull.simplices

triangles = []
for i in hull.simplices:
    i = np.append(i, i[0])
    plt.plot(nodes[i,0], nodes[i,1], nodes[i,2])

for j in faces:
    tri = [(nodes[j[0], 0], nodes[j[0], 1], nodes[j[0], 2]), (nodes[j[1], 0], nodes[j[1], 1], nodes[j[1], 2]), (nodes[j[2], 0], nodes[j[2], 1], nodes[j[2], 2])]
    triangles.append(tri)


new_faces = simplify(triangles)

for sq in new_faces:
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot(verts[0], verts[1], verts[2], 'bo')
    ax.dist=10
    ax.azim=110
    ax.elev=30
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    f = a3.art3d.Poly3DCollection([sq])
    f.set_color(colors.rgb2hex(sp.rand(3)))
    f.set_edgecolor('k')
    # f.set_alpha(0.1)
    ax.add_collection3d(f)

    plt.show()

Answer 1

您需要一些共面性的公差标准。

在您的情况下，只需将函数 is_coplanar(a, b, tolerance_in_radians=0) 中的 tolerance_in_radians 关键字设置为适当的值。

其他：

• 计算表面法线，也许你的凸包已经返回它们。 它应该返回具有固定顺序（顺时针，逆时针）的三角形

• 计算节点的质心

• 作为检查计算表面法线是否指向远离质心

• 计算曲面法线之间的角度

• 对于低于 0.5° 的角度，三角形对于您的情况来说足够“共面”