python - Matplotlib - 提取 3D 多边形图的 2D 轮廓

问题描述

我有一个由一组 Poly3DCollection 定义的 3D 图。集合的每个多边形都包含一个 3D 单纯形列表（单纯形 = 4 个点），如下所示。

[[[21096.4, 15902.1, 74.3],  
  [21098.5, 15904.3, 54.7],
  [21114.2, 15910.1, 63.0],
  [21096.4, 15902.1, 74.3]],
  ...
 [[21096.4, 15902.1, 74.3],
  [21114.8, 15909.9, 91.3],
  [21114.2, 15910.1, 63.0],
  [21096.4, 15902.1, 74.3]]]

从这些集合中，我绘制了一个 3D 网格给我这个结果

我想确定这个 3D 网格在投影到 2D 时的轮廓，以便在屏幕上绘制它，以便突出显示它。理想情况下，它会给我类似的东西

有没有办法做到这一点？

为了实现它，我正在考虑类似的东西

通过将每个点乘以 matplotlib 必须在内部进行最终渲染的投影矩阵来获取我的 3D 点的 2D 坐标，或者直接从 matplotlib 内部获取投影的 2D 坐标，我不知道它是否是可能的。
将某种 2D 轮廓检测算法应用于步骤 1 中的 2D 坐标
将在步骤 2 中找到的 2D 轮廓添加到现有的 3D 图

但是，我没有找到任何方法从我的 matplotlib Axes3D 对象公开的界面中实现这种轮廓检测。

只要我能实现绘制 2D 轮廓，如果我直接在我的原始 3D 数据集和投影上或从我的 matplotlib Axes3D 对象确定它，对我来说并不重要。

标签： pythonmatplotlib3dgeometrycontour

事实证明，这比我最初预期的要复杂得多。我解决它的方法是首先将对象旋转到正面视图（根据角度Axes3D elev和azim角度），将其投影到 yz 平面上，计算 2D 轮廓，重新添加第三维，然后旋转现在的 3D 轮廓回到当前视图。

旋转部分可以通过简单的矩阵运算来完成，只需注意 x、y 和 z 轴可能会被拉伸，并且需要在旋转前取消拉伸。

投影部分有点棘手，因为我不知道有什么聪明的方法可以找到如此大的点集合的外部点。因此，我通过分别计算每个单纯形的投影来解决它，计算它们的 2D 凸包（使用scipy），将它们转换为shapely多边形，最后计算所有这些多边形的并集。然后我添加了丢失的 x 坐标并将整个东西旋转回当前视图。

默认情况下，Axes3D对象使用透视，导致对象的实际轮廓与计算的投影不完全对齐。这可以通过使用正交视图（用设置ax.set_proj_type('ortho')）来避免。

最后，一旦图像旋转，轮廓/投影需要更新。因此，我按照这个示例将整个函数添加到事件队列中。

请询问是否还有其他问题。

from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
from matplotlib import pyplot as plt
import numpy as np

from shapely.geometry import Polygon
from scipy.spatial import ConvexHull

from scipy.spatial import Delaunay

##the figure
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))

##generating some random points:
points = np.random.rand(50,3)
xmin,xmax = 0,100
ymin,ymax = -10,10
zmin,zmax = -20,20
points[:,1] = (points[:,1]*(ymax-ymin)+ymin) * np.sin(points[:,0]*np.pi)
points[:,2] = (points[:,2]*(zmax-zmin)+zmin) * np.sin(points[:,0]*np.pi)
points[:,0] *= 100


##group them into simlices
tri =  Delaunay(points)
simplex_coords = np.array([tri.points[simplex] for simplex in tri.simplices])

##plotting the points
ax.scatter(points[:,0], points[:,1], points[:,2])

##visualizing simplices
line_coords = np.array(
    [[c[i],c[j]] for c in simplex_coords for i in range(len(c)) for j in range(i+1,len(c))]
)
simplex_lines = Line3DCollection(line_coords, colors='k', linewidths=1, zorder=10)
ax.add_collection3d(simplex_lines)    

##adjusting plot
ax.set_xlim([xmin,xmax])
ax.set_xlabel('x')
ax.set_ylim([2*ymin,2*ymax])
ax.set_ylabel('y')
ax.set_zlim([2*zmin,2*zmax])
ax.set_zlabel('z')


def compute_2D_outline():
    """
    Compute the outline of the 2D projection of the 3D mesh and display it as
    a Poly3DCollection or a Line3DCollection.
    """

    global collection
    global lines
    global elev
    global azim

    ##remove the previous projection (if it has been already created)
    try:
        collection.remove()
        lines.remove()
    except NameError as e:
        pass


    ##storing current axes orientation
    elev = ax.elev
    azim = ax.azim

    ##convert angles
    theta = -ax.elev*np.pi/180
    phi = -ax.azim*np.pi/180

    #the extend of each of the axes:
    diff = lambda t: t[1]-t[0]
    lx = diff(ax.get_xlim())
    ly = diff(ax.get_ylim())
    lz = diff(ax.get_zlim())

    ##to compute the projection, we 'unstretch' the axes and rotate them
    ##into the (elev=0, azmi=0) orientation
    stretch = np.diag([1/lx,1/ly,1/lz])
    rot_theta = np.array([
        [np.cos(theta), 0, -np.sin(theta)],
        [0, 1, 0],
        [np.sin(theta), 0,  np.cos(theta)],
    ])
    rot_phi = np.array([
        [np.cos(phi), -np.sin(phi), 0],
        [np.sin(phi),  np.cos(phi), 0],
        [0,0,1],
    ])
    rot_tot = np.dot(rot_theta,np.dot(rot_phi,stretch))

    ##after computing the outline, we will have to reverse this operation:
    bstretch = np.diag([lx,ly,lz])
    brot_theta = np.array([
        [ np.cos(theta), 0, np.sin(theta)],
        [0, 1, 0],
        [-np.sin(theta), 0, np.cos(theta)],
    ])
    brot_phi = np.array([
        [ np.cos(phi),  np.sin(phi), 0],
        [-np.sin(phi),  np.cos(phi), 0],
        [0,0,1],
    ])
    brot_tot = np.dot(np.dot(bstretch,brot_phi),brot_theta)

    ##To get the exact outline, we will have to compute the projection of each simplex
    ##separately and compute the convex hull of the projection. We then use shapely to
    ##compute the unity of all these convex hulls to get the projection (or shadow).
    poly = None
    for simplex in simplex_coords:
        simplex2D = np.dot(rot_tot,simplex.T)[1:].T
        hull = simplex2D[ConvexHull(simplex2D).vertices]
        if poly is None:
            poly = Polygon(hull)
        else:
            poly = poly.union(Polygon(hull))

    ##the 2D points of the final projection have to be made 3D and transformed back
    ##into the correct axes rotation
    outer_points2D = np.array(poly.exterior.coords.xy)
    outer_points3D = np.concatenate([[np.zeros(outer_points2D.shape[1])],outer_points2D])    
    outer_points3D_orig = np.dot(brot_tot, outer_points3D)

    ##adding the polygons
    collection = Poly3DCollection(
        [outer_points3D_orig.T], alpha=0.25, facecolor='b', zorder=-1
    )
    ax.add_collection3d(collection)

    ##adding the lines
    lines = Line3DCollection(
        [outer_points3D_orig.T], alpha=0.5, colors='r', linewidths=5, zorder=5
    )
    ax.add_collection3d(lines)    


def on_move(event):
    """
    For tracking rotations of the Axes3D object
    """

    if event.inaxes == ax and (elev != ax.elev or azim != ax.azim):
        compute_2D_outline()        
    fig.canvas.draw_idle()

##initial outline:
compute_2D_outline()

##the entire thing will only work correctly with an orthogonal view
ax.set_proj_type('ortho')

##saving ax.azim and ax.elev for on_move function
azim = ax.azim
elev = ax.elev

##adding on_move to the event queue
c1 = fig.canvas.mpl_connect('motion_notify_event', on_move)

plt.show()

最终结果（带有一些生成的随机数据）看起来像这样：

python - Matplotlib - 提取 3D 多边形图的 2D 轮廓

问题描述

解决方案

推荐阅读