Это оказалось намного сложнее, чем я сначала ожидал.Я решил, что сначала нужно повернуть объект во фронтальный вид (с точки зрения углов Axes3D elev и azim), спроецировать его на плоскость yz, вычислить двухмерный контур, заново добавитьтретье измерение, а затем поворачивая теперь 3D-контур обратно в текущий вид.
Поворотная часть может быть выполнена с помощью простой матричной операции, нужно только позаботиться о том, чтобы оси x, y и z могли быть растянуты, и перед вращением их необходимо растянуть.
Часть проекции была немного хитрой, поскольку я не знаю ни одного умного способа найти внешние точки такого большого набора точек.Поэтому я решил ее, вычислив проекцию каждого симплекса отдельно, вычислил их двумерные выпуклые оболочки (используя 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()
Окончательный результат (с некоторымислучайные данные) выглядит примерно так:
