在 matplotlib 中绘制 3D 凸封闭区域
Plot 3D convex closed regions in matplotlib
我正在尝试在 3D 中绘制由一组不等式定义的多面体。本质上，我尝试在 matplotlib 中重现此 matlab plotregion 库的功能。
我的方法是获取相交顶点，构建它们的凸包，然后获取并绘制结果面（单纯形）。
问题是许多单纯形是共面的，它们无缘无故地使情节非常繁忙（请参见下图中的所有这些对角线边缘）。
有什么简单的方法可以只打印多面体的 "outside" 条边，而不必自己一个接一个地合并所有共面单形？
谢谢
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 中没有任何原生内容。不过，找到属于一起的面孔并不是特别困难。下面实现的基本思想是创建一个图形，其中每个节点都是一个三角形。然后连接 co-planar 和相邻的三角形。最后，您找到图形的连通分量以确定哪些三角形构成面。
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 可能需要做更多的工作。很确定这会因奇怪的/ non-convex 形状而中断，而且我什至不能 100% 确定它是否始终适用于凸形。不过休息应该没问题。
以下是我的解决方案版本。它类似于@Paul 的解决方案，因为它采用三角形，按它们所属的面对它们进行分组，并将它们连接到一个面。
不同之处主要在于此解决方案不使用 nx 或 simpy。许多必要的操作都是通过重建索引、广泛使用 unique 和一些线性代数来执行的。
最终面的顶点顺序由 ConvexHull 决定。我认为这不应该是一个限制，因为（我认为）任何一半 space 交叉点应该只产生凸形。但是，我还添加了另一种方法，如果形状不是凸的，可以使用它（基于 this question 的想法）。
from scipy.spatial import HalfspaceIntersection
from scipy.spatial import ConvexHull
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3
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)
simplices = hull.simplices
org_triangles = [verts[s] for s in simplices]
class Faces():
def __init__(self,tri, sig_dig=12, method="convexhull"):
self.method=method
self.tri = np.around(np.array(tri), sig_dig)
self.grpinx = list(range(len(tri)))
norms = np.around([self.norm(s) for s in self.tri], sig_dig)
_, self.inv = np.unique(norms,return_inverse=True, axis=0)
def norm(self,sq):
cr = np.cross(sq[2]-sq[0],sq[1]-sq[0])
return np.abs(cr/np.linalg.norm(cr))
def isneighbor(self, tr1,tr2):
a = np.concatenate((tr1,tr2), axis=0)
return len(a) == len(np.unique(a, axis=0))+2
def order(self, v):
if len(v) <= 3:
return v
v = np.unique(v, axis=0)
n = self.norm(v[:3])
y = np.cross(n,v[1]-v[0])
y = y/np.linalg.norm(y)
c = np.dot(v, np.c_[v[1]-v[0],y])
if self.method == "convexhull":
h = ConvexHull(c)
return v[h.vertices]
else:
mean = np.mean(c,axis=0)
d = c-mean
s = np.arctan2(d[:,0], d[:,1])
return v[np.argsort(s)]
def simplify(self):
for i, tri1 in enumerate(self.tri):
for j,tri2 in enumerate(self.tri):
if j > i:
if self.isneighbor(tri1,tri2) and \
self.inv[i]==self.inv[j]:
self.grpinx[j] = self.grpinx[i]
groups = []
for i in np.unique(self.grpinx):
u = self.tri[self.grpinx == i]
u = np.concatenate([d for d in u])
u = self.order(u)
groups.append(u)
return groups
f = Faces(org_triangles)
g = f.simplify()
ax = a3.Axes3D(plt.figure())
colors = list(map("C{}".format, range(len(g))))
pc = a3.art3d.Poly3DCollection(g, facecolor=colors,
edgecolor="k", alpha=0.9)
ax.add_collection3d(pc)
ax.dist=10
ax.azim=30
ax.elev=10
ax.set_xlim([0,5])
ax.set_ylim([0,5])
ax.set_zlim([0,5])
plt.show()
我正在尝试在 3D 中绘制由一组不等式定义的多面体。本质上，我尝试在 matplotlib 中重现此 matlab plotregion 库的功能。
我的方法是获取相交顶点，构建它们的凸包，然后获取并绘制结果面（单纯形）。
问题是许多单纯形是共面的，它们无缘无故地使情节非常繁忙（请参见下图中的所有这些对角线边缘）。
有什么简单的方法可以只打印多面体的 "outside" 条边，而不必自己一个接一个地合并所有共面单形？
谢谢
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 中没有任何原生内容。不过，找到属于一起的面孔并不是特别困难。下面实现的基本思想是创建一个图形，其中每个节点都是一个三角形。然后连接 co-planar 和相邻的三角形。最后，您找到图形的连通分量以确定哪些三角形构成面。
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 可能需要做更多的工作。很确定这会因奇怪的/ non-convex 形状而中断，而且我什至不能 100% 确定它是否始终适用于凸形。不过休息应该没问题。
以下是我的解决方案版本。它类似于@Paul 的解决方案，因为它采用三角形，按它们所属的面对它们进行分组，并将它们连接到一个面。
不同之处主要在于此解决方案不使用 nx 或 simpy。许多必要的操作都是通过重建索引、广泛使用 unique 和一些线性代数来执行的。
最终面的顶点顺序由 ConvexHull 决定。我认为这不应该是一个限制，因为（我认为）任何一半 space 交叉点应该只产生凸形。但是，我还添加了另一种方法，如果形状不是凸的，可以使用它（基于 this question 的想法）。
from scipy.spatial import HalfspaceIntersection
from scipy.spatial import ConvexHull
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3
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)
simplices = hull.simplices
org_triangles = [verts[s] for s in simplices]
class Faces():
def __init__(self,tri, sig_dig=12, method="convexhull"):
self.method=method
self.tri = np.around(np.array(tri), sig_dig)
self.grpinx = list(range(len(tri)))
norms = np.around([self.norm(s) for s in self.tri], sig_dig)
_, self.inv = np.unique(norms,return_inverse=True, axis=0)
def norm(self,sq):
cr = np.cross(sq[2]-sq[0],sq[1]-sq[0])
return np.abs(cr/np.linalg.norm(cr))
def isneighbor(self, tr1,tr2):
a = np.concatenate((tr1,tr2), axis=0)
return len(a) == len(np.unique(a, axis=0))+2
def order(self, v):
if len(v) <= 3:
return v
v = np.unique(v, axis=0)
n = self.norm(v[:3])
y = np.cross(n,v[1]-v[0])
y = y/np.linalg.norm(y)
c = np.dot(v, np.c_[v[1]-v[0],y])
if self.method == "convexhull":
h = ConvexHull(c)
return v[h.vertices]
else:
mean = np.mean(c,axis=0)
d = c-mean
s = np.arctan2(d[:,0], d[:,1])
return v[np.argsort(s)]
def simplify(self):
for i, tri1 in enumerate(self.tri):
for j,tri2 in enumerate(self.tri):
if j > i:
if self.isneighbor(tri1,tri2) and \
self.inv[i]==self.inv[j]:
self.grpinx[j] = self.grpinx[i]
groups = []
for i in np.unique(self.grpinx):
u = self.tri[self.grpinx == i]
u = np.concatenate([d for d in u])
u = self.order(u)
groups.append(u)
return groups
f = Faces(org_triangles)
g = f.simplify()
ax = a3.Axes3D(plt.figure())
colors = list(map("C{}".format, range(len(g))))
pc = a3.art3d.Poly3DCollection(g, facecolor=colors,
edgecolor="k", alpha=0.9)
ax.add_collection3d(pc)
ax.dist=10
ax.azim=30
ax.elev=10
ax.set_xlim([0,5])
ax.set_ylim([0,5])
ax.set_zlim([0,5])
plt.show()