I have a cloud of 8 nodes which form a prism. I need the Information which node belongs to which side of the prism. It might be, that a side of the prism is plane, but this is not always the case.
So I thought I can search for triangles with ConvexHull and then find coplanar triangles with this answer from . This worked pretty well for a symmetric prism (like a cuboid) but in my case the object has no rectangular edges at all. When transferring the code to my problem, unfortunately always "holes" remain in the surface.
Here is my code so far.
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()
As you can see in the screenshot below I have a hole at the top of my prism. Do you know how I can fix this and get all triangles or quadrangles of the hull of my prism?
Thank you!
Screenshot:
Edit: I made some changes to my code (see below), so now all triangles/quadrangles are plottet seperatly and all faces are found. Could it be a plotting bug?
# 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()
