在 matplotlib 中绘制 3D 凸封闭区域 [英] Plot 3D convex closed regions in matplotlib

问题描述

我试图在 3D 中绘制由一组不等式定义的多面体.本质上，我尝试重现这个 matlab

解决方案

很确定 matplotlib 中没有原生的东西.不过，找到属于一起的面孔并不是特别困难.下面实现的基本思想是创建一个图形，其中每个节点都是一个三角形.然后连接共面且相邻的三角形.最后，您找到图形的连通分量以确定哪些三角形形成一个面.

将 numpy 导入为 np从 sympy 导入平面，Point3D将 networkx 导入为 nxdef 简化(三角形):"""简化可迭代的三角形，使相邻和共面的三角形形成一个面.每个三角形是 3D 空间中的一组 3 个点."""# 创建一个图，其中节点代表三角形；# 如果对应的三角形相邻且共面，则节点连接G = nx.Graph()G.add_nodes_from(range(len(triangles)))对于 ii, a in enumerate(triangles):对于 jj, b 在 enumerate(triangles):if (ii < jj): # 仅以一种方式测试关系，因为邻接和共面是双射的如果 is_adjacent(a, b):如果 is_coplanar(a, b, np.pi/180.):G.add_edge(ii,jj)# 属于连通分量的三角形可以组合组件 = 列表(nx.connected_components(G))简化 = [set(flatten(triangles[index] for index in component)) for component in components]# 需要对节点重新排序，以便正确绘制补丁reordered = [reorder(face) for face in simple]退货重新排序def is_adjacent(a, b):return len(set(a) & set(b)) == 2 # 即三角形共享 2 个点，因此共享一个边def is_coplanar(a, b, tolerance_in_radians=0):a1, a2, a3 = ab1, b2, b3 = b平面_a = 平面(Point3D(a1), Point3D(a2), Point3D(a3))平面_b = 平面(Point3D(b1), Point3D(b2), Point3D(b3))如果不是tolerance_in_radians: # 只接受精确的结果返回plane_a.is_coplanar(plane_b)别的:角度 = plane_a.angle_between(plane_b).evalf()angle %= np.pi # 确保角度在 0 和 np.pi 之间返回(角度-tolerance_in_radians <= 0.)或((np.pi - 角度) - tolerance_in_radians <= 0.)flatten = lambda l: [l 中子列表的项目，子列表中的项目]定义重新排序(顶点):"""重新排序节点，使生成的路径对应于点集的外壳".笔记:-----未在边缘情况下进行测试，并且可能会损坏.可能只适用于凸面形状."""if len(vertices) <= 3: # 只是一个三角形返回顶点别的:# 取随机顶点(这里只是第一个)重新排序 = [vertices.pop()]# 获取下一个尚未重新排序的最近顶点# 重复直到只有一个顶点留在原始列表中顶点 = 列表(顶点)而 len(vertices) >1:idx = np.argmin(get_distance(reordered[-1], vertices))v = vertices.pop(idx)重新排序.append(v)# 将剩余的顶点添加到输出重新排序 += 顶点退货重新排序def get_distance(v1, v2):v2 = np.array(list(v2))差异 = v2 - v1ssd = np.sum(difference**2,axis=1)返回 np.sqrt(ssd)

应用于您的示例:

from scipy.spatial import HalfspaceIntersection从 scipy.spatial 导入 ConvexHull将 scipy 导入为 sp将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt将 mpl_toolkits.mplot3d 导入为 a3导入 matplotlib.colors 作为颜色w = np.array([1., 1., 1.])# ∑ᵢ hᵢ wᵢ qᵢ - ∑ᵢ gᵢ wᵢ <= 0# qᵢ - ubᵢ <= 0# -qᵢ + lbᵢ <= 0halfspaces = 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]])可行点 = np.array([0.1, 0.1, 0.1])hs = HalfspaceIntersection(半空间，可行点)verts = hs.intersections船体 = ConvexHull(verts)面孔 = hull.simplicesax = a3.Axes3D(plt.figure())ax.dist=10ax.azim=30ax.elev=10ax.set_xlim([0,5])ax.set_ylim([0,5])ax.set_zlim([0,5])三角形 = []对于 s 在脸上:平方 = [(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])]三角形.追加(平方)new_faces = 简化(三角形)对于 new_faces 中的 sq: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)# # 旋转轴并更新# 范围内的角度(0, 360):# ax.view_init(30, 角度)# plt.draw()# plt.pause(.001)

注意

经过反思，reordered 函数可能需要更多的工作.很确定这会打破奇怪/非凸形状，我什至不能 100% 确定它总是适用于凸形状.不过休息应该没问题.

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?

Thank you

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()

解决方案

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)

Applied to your example:

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)

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.

这篇关于在 matplotlib 中绘制 3D 凸封闭区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！