问题描述
我试图以3D方式绘制由一组不等式定义的多面体.本质上，我尝试重现此Matlab的功能 plotregion 库.
I am trying to plot in 3D a polytope defined by a set of inequalities. Essentially, I try to reproduce the functionality of this matlab plotregion library in matplotlib.
我的方法是获取相交顶点，构造它们的凸包，然后获取并绘制所得面(单形).
My approach is to get the intersection vertices, construct the convex hull of them, and then get and plot the resulting faces (simplices).
问题在于许多单纯形是共面的，它们无缘无故地使绘图变得非常繁忙(请参见下图中的所有这些对角线边缘).
The problem is that many simplices are coplanar, and they are making the plot very busy for no reason (see all of these diagonal edges in the plot below).
是否有任何简便的方法可以只打印多面体的外部"边缘，而不必由我自己一个接一个地合并所有共面的单纯形?
Is there any easy way to just print the "outside" edges of the polyhedron, without having to consolidate by my self, one by one, all of the coplanar simplices?
谢谢
from scipy.spatial import HalfspaceIntersection
from scipy.spatial import ConvexHull
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3
import matplotlib.colors as colors
w = np.array([1., 1., 1.])
# ∑ᵢ hᵢ wᵢ qᵢ - ∑ᵢ gᵢ wᵢ <= 0
# qᵢ - ubᵢ <= 0
# -qᵢ + lbᵢ <= 0
halfspaces = np.array([
[1.*w[0], 1.*w[1], 1.*w[2], -10 ],
[ 1., 0., 0., -4],
[ 0., 1., 0., -4],
[ 0., 0., 1., -4],
[-1., 0., 0., 0],
[ 0., -1., 0., 0],
[ 0., 0., -1., 0]
])
feasible_point = np.array([0.1, 0.1, 0.1])
hs = HalfspaceIntersection(halfspaces, feasible_point)
verts = hs.intersections
hull = ConvexHull(verts)
faces = hull.simplices
ax = a3.Axes3D(plt.figure())
ax.dist=10
ax.azim=30
ax.elev=10
ax.set_xlim([0,5])
ax.set_ylim([0,5])
ax.set_zlim([0,5])
for s in faces:
sq = [
[verts[s[0], 0], verts[s[0], 1], verts[s[0], 2]],
[verts[s[1], 0], verts[s[1], 1], verts[s[1], 2]],
[verts[s[2], 0], verts[s[2], 1], verts[s[2], 2]]
]
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()
推荐答案
很确定matplotlib中没有任何本机.不过，找到属于彼此的面孔并不是特别困难.下面实现的基本思想是创建一个图形，其中每个节点都是一个三角形.然后，您可以连接共面和相邻的三角形.最后，您找到图形的连接组件，以确定哪些三角形构成了一个面.
Pretty sure there is nothing native in matplotlib. Finding the faces that belong together is not particularly hard, though. The basic idea implemented below is that you create a graph, where each node is a triangle. You then connect triangles that are co-planar and adjacent. Finally, you find the connected components of the graph to determine which triangles form a face.
import numpy as np
from sympy import Plane, Point3D
import networkx as nx
def simplify(triangles):
"""
Simplify an iterable of triangles such that adjacent and coplanar triangles form a single face.
Each triangle is a set of 3 points in 3D space.
"""
# create a graph in which nodes represent triangles;
# nodes are connected if the corresponding triangles are adjacent and coplanar
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): # test relationships only in one way as adjacency and co-planarity are bijective
if is_adjacent(a, b):
if is_coplanar(a, b, np.pi / 180.):
G.add_edge(ii,jj)
# triangles that belong to a connected component can be combined
components = list(nx.connected_components(G))
simplified = [set(flatten(triangles[index] for index in component)) for component in components]
# need to reorder nodes so that patches are plotted correctly
reordered = [reorder(face) for face in simplified]
return reordered
def is_adjacent(a, b):
return len(set(a) & set(b)) == 2 # i.e. triangles share 2 points and hence a side
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):
"""
Reorder nodes such that the resulting path corresponds to the "hull" of the set of points.
Note:
-----
Not tested on edge cases, and likely to break.
Probably only works for convex shapes.
"""
if len(vertices) <= 3: # just a triangle
return vertices
else:
# take random vertex (here simply the first)
reordered = [vertices.pop()]
# get next closest vertex that is not yet reordered
# repeat until only one vertex remains in original list
vertices = list(vertices)
while len(vertices) > 1:
idx = np.argmin(get_distance(reordered[-1], vertices))
v = vertices.pop(idx)
reordered.append(v)
# add remaining vertex to output
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)
应用于您的示例:
from scipy.spatial import HalfspaceIntersection
from scipy.spatial import ConvexHull
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3
import matplotlib.colors as colors
w = np.array([1., 1., 1.])
# ∑ᵢ hᵢ wᵢ qᵢ - ∑ᵢ gᵢ wᵢ <= 0
# qᵢ - ubᵢ <= 0
# -qᵢ + lbᵢ <= 0
halfspaces = np.array([
[1.*w[0], 1.*w[1], 1.*w[2], -10 ],
[ 1., 0., 0., -4],
[ 0., 1., 0., -4],
[ 0., 0., 1., -4],
[-1., 0., 0., 0],
[ 0., -1., 0., 0],
[ 0., 0., -1., 0]
])
feasible_point = np.array([0.1, 0.1, 0.1])
hs = HalfspaceIntersection(halfspaces, feasible_point)
verts = hs.intersections
hull = ConvexHull(verts)
faces = hull.simplices
ax = a3.Axes3D(plt.figure())
ax.dist=10
ax.azim=30
ax.elev=10
ax.set_xlim([0,5])
ax.set_ylim([0,5])
ax.set_zlim([0,5])
triangles = []
for s in faces:
sq = [
(verts[s[0], 0], verts[s[0], 1], verts[s[0], 2]),
(verts[s[1], 0], verts[s[1], 1], verts[s[1], 2]),
(verts[s[2], 0], verts[s[2], 1], verts[s[2], 2])
]
triangles.append(sq)
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)
# # rotate the axes and update
# for angle in range(0, 360):
# ax.view_init(30, angle)
# plt.draw()
# plt.pause(.001)
注意
反射时，功能reordered可能需要做更多的工作.可以肯定的是，这对于怪异的/非凸的形状会破裂，而且我什至不能100％地确定它总是可以用于凸形.休息应该没问题.
Note
Upon reflection, the function reordered probably needs some more work. Pretty sure this will break for weird / non-convex shapes, and I am not even 100% sure that it will always work for convex shapes. Rest should be fine though.
这篇关于在Matplot lib中绘制3D凸封闭区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！