Andras Deak 给出了一个非常全面的answer 讨论了 matplotlibs 3D 绘图功能在手头任务中的问题/限制。在他的回答结束时——以肯定的结尾——他给出了一个使用替代库的解决方案。
我开始尝试在 matplotlib 中找到一个 hacky/专业的解决方案。让我先说明原因。我想在 2D 表面上绘制轨迹，我开始使用 matplotlib。我将它用于我所有的 2D 绘图，并希望为这个特定的 3D 绘图应用程序找到解决方案。 matplotlibs 3D 图的好处在于它们是矢量化的，因为它们基本上只是通过将 3D 元素投影到相机平面上并覆盖它们（根据它们与相机的距离按顺序绘制它们）生成的 2D 图。可以为绘图的每个元素单独控制光栅化，而不会影响轴、标签等。使用光线追踪的“真实”3D 绘图库通常无法生成完全矢量化的绘图。我认为mayavi 就是一个例子，我知道 Mathematica 在这方面也非常有限。
提出我的解决方案：我查看了 plot_surface 的代码，该代码最终基于 Poly3DCollection，以了解 matplotlib 如何决定首先绘制表面上的哪些多边形/元素。的方法_do_3d_projection
Poly3DCollection 将投影到 2d 相机平面上的多边形按（原始 3D 对象的）到相机的距离排序。首先绘制远离相机的元素，然后绘制靠近相机的元素。这对于大多数绘图都可以很好地创建正确的视角（但该方法有局限性，例如，请参阅mplot3d FAQ。但是，这种排序是我解决方案的关键。给定一组点 pts和一个表面surf（必须使用show 或savefig 绘制才能设置其相机/投影变量）：
-
surf 中所有 3D 多边形的 2D 投影 segments_2d 到相机平面的计算包括它们基于到相机的距离的排序（存储在 segments_idxs ）。
- 所有点都与 3D 表面上的元素/多边形相关联。
- 计算 3D 点到相机平面的 2D 投影。
- 为了确定一个点是否可见，我们检查它是否被一个多边形覆盖在它所关联的那个之后（从第 2 步开始）。为此，我们使用来自
matplotlib.path 的contains_points 方法，另请参阅相关问题What's the fastest way of checking if a point is inside a polygon in python。
- 我包含了一个动态更新（改编自How to obscure a line behind a surface plot in matplotlib?）。警告：具有大量多边形的表面的代码/绘图可能会变得非常缓慢。
这里是必要的代码/最小工作示例，表面由 OP 给出，样本点位于单位圆。
import matplotlib.pyplot as plt
import numpy as np
import copy
import matplotlib.path as mpltPath
from mpl_toolkits.mplot3d import proj3d
from matplotlib import cm
def clip_on_surface(surf,pts):
## Get projection of 3d surface onto 2d camera plane
## [Code form [mpl_toolkits/mplot3d/art3d.py - Poly3DCollection._do_3d_projection(self, renderer=None)] to ]
txs, tys, tzs = proj3d._proj_transform_vec(surf._vec, surf.axes.M)
xyzlist = [(txs[sl], tys[sl], tzs[sl]) for sl in surf._segslices]
cface = surf._facecolor3d
cedge = surf._edgecolor3d
if len(cface) != len(xyzlist):
cface = cface.repeat(len(xyzlist), axis=0)
if len(cedge) != len(xyzlist):
if len(cedge) == 0:
cedge = cface
else:
cedge = cedge.repeat(len(xyzlist), axis=0)
if xyzlist:
# sort by depth (furthest drawn first)
z_segments_2d = sorted(
((surf._zsortfunc(zs), np.column_stack([xs, ys]), fc, ec, idx)
for idx, ((xs, ys, zs), fc, ec)
in enumerate(zip(xyzlist, cface, cedge))),
key=lambda x: x[0], reverse=True)
# z_segments_2d = sorted(z_segments_2d,key=lambda x:x[4])
segments_zorder, segments_2d, facecolors2d, edgecolors2d, segments_idxs = zip(*z_segments_2d)
segments_paths = [mpltPath.Path(x) for x in segments_2d]
## Get polygons in 3d space
xs, ys, zs = surf._vec[0:3,:]
xyzlist = [(xs[sl], ys[sl], zs[sl]) for sl in surf._segslices]
segments_3d=[]
segments_3d_centroid=[]
for q in xyzlist:
vertices = np.transpose(np.array([q[0],q[1],q[2]]))
segments_3d.append( vertices )
segments_3d_centroid.append( sum(list(vertices))/len(list(vertices)) ) # centroid of polygon (mean of vertices)
## Process points
pts_info = [[0,0,True] for x in range(len(pts))]
# 0: index of closest 3d polygon
# 1: index of closest 3d polygon in segments_idxs: drawing order
# 2: True if visible (not overlapped by polygons drawn after associated polygon), False else
pts_visible = copy.copy(pts) # visible points (invisible set to np.nan)
pts_invisible = copy.copy(pts) # invisible points (visible set to np.nan)
# compute pts_info[:,0] and pts_info[:,1] -- index of closest 3d polygon and its position in segments_idxs
for i in range(len(pts)):
# Associate by distance
dist = np.inf
x=[pts[i][0],pts[i][1],pts[i][2]]
for j in range(len(segments_3d_centroid)):
yc=segments_3d_centroid[j]
dist_tmp = np.sqrt( (x[0]-yc[0])**2 + (x[1]-yc[1])**2 + (x[2]-yc[2])**2 )
if dist_tmp<dist:
dist=dist_tmp
pts_info[i][0]=j
pts_info[i][1] = segments_idxs.index( pts_info[i][0] )
# compute projection of 3d points into 2d camera plane
pts_2d_x, pts_2d_y, pts_2d_z = proj3d._proj_transform_vec(np.transpose(np.array([[x[0],x[1],x[2],1.0] for x in pts])), surf.axes.M)
# decide visibility
for i in range(len(pts_info)):
for j in range(pts_info[i][1]+1,len(segments_paths)):
b=segments_paths[j].contains_points( [[pts_2d_x[i],pts_2d_y[i]]] )
if b==True:
pts_info[i][2]=False
break
if pts_info[i][2]:
pts_invisible[i][0]=np.nan
pts_invisible[i][1]=np.nan
pts_invisible[i][2]=np.nan
else:
pts_visible[i][0]=np.nan
pts_visible[i][1]=np.nan
pts_visible[i][2]=np.nan
return { 'pts_visible': pts_visible, 'pts_invisible':pts_invisible, 'pts_info':pts_info }
def f(x, y):
return np.sin(2*x) * np.cos(2*y)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev=30., azim=55.)
# Generate surface plot (surf)
xs = np.linspace(-2, 2, 25)
ys = np.linspace(-2, 2, 25)
Xs, Ys = np.meshgrid(xs, ys)
zs = np.array(f(np.ravel(Xs), np.ravel(Ys)))
Zs = zs.reshape(Xs.shape)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
surf = ax.plot_surface(Xs, Ys, Zs, rstride=1, cstride=1,
cmap=cm.get_cmap('viridis'),linewidth=0.0,edgecolor='black',
antialiased=True,rasterized=False)
# Generate pts on surf
t = np.linspace(0, 1, 200)
xp = np.sin(t*2*np.pi)
yp = np.cos(t*2*np.pi)
zp = f(xp,yp)
pts=np.transpose(np.array([xp,yp,zp]))
def rotate(event):
if event.inaxes == ax:
surf_pts=clip_on_surface(surf,pts)
ax.plot(surf_pts['pts_visible'][:,0],surf_pts['pts_visible'][:,1],surf_pts['pts_visible'][:,2],'.', zorder=10,c='red',markersize=2)
ax.plot(surf_pts['pts_invisible'][:,0],surf_pts['pts_invisible'][:,1],surf_pts['pts_invisible'][:,2],'.', zorder=10,c='green',markersize=2)
fig.canvas.draw_idle()
c1 = fig.canvas.mpl_connect('motion_notify_event', rotate)
plt.show()
代码仍然有点混乱，它仍然不能完美运行，但这里有一些结果，表面上有 25*25=625 个四边形，单位圆上有 200 个点。
红色点是可见点，绿色点是不可见点（此处为说明目的而绘制，但为了最初的问题/问题，人们会诅咒从图中省略它们）。有些点应该清晰可见，但被检测为不可见。我还不确定那里出了什么问题，但对我来说，这种有限的未检测到并没有太大问题，因为我最终想绘制很多（任意密集）点的线/轨迹。如果未命中的不聚集，我可以忍受一些丢失的。
另一个固有的问题/限制是，当前的方法没有真正的概念，即点是在表面之上还是之下，这意味着从表面下方看时，表面之上/之上的点是可见的。这是此行为的示例：
这与 Andras Deak 已经提出的观点相关，即当前的问题在没有额外限制的情况下有些不明确或至少模棱两可。例如，可以要求将所有点放置在指向相机的表面上。在目前的方法中实现这一点是困难的。在几何方面，当前的实现将有限大小的球放置在无穷小的多边形上，使它们从两侧都可见（这在某些用例中实际上可能是可行的/理想的）。
代码仍在进行中，如果我发现重大改进，我可能会更新此答案。非常欢迎对一般方法和/或实施发表评论。我绝不是 python 专家（我几乎只将它用于绘图和相关的非常轻量级的数据处理），因此它们在代码性能和范围方面可能有很大的改进空间。