Почти уверен, что в matplotlib нет ничего нативного. Однако найти лица, которые подходят друг другу, не особенно сложно. Основная идея, реализованная ниже, заключается в том, что вы создаете граф, где каждый узел представляет собой треугольник. Затем вы соединяете треугольники, лежащие в одной плоскости и смежные. Наконец, вы находите связанные компоненты графа, чтобы определить, какие треугольники образуют грань.
![введите описание изображения здесь](https://i.stack.imgur.com/2b7eW.png)
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
, вероятно, нуждается в доработке. Почти уверен, что это сломается для странных/невыпуклых форм, и я даже не уверен на 100%, что это всегда будет работать для выпуклых форм. Хотя отдых должен быть в порядке.
05.03.2018
nx
и сдался после того, как заблудился в 5-м подцикле или около того. 05.03.2018is_coplanar
был бы намного короче, если бы не ошибки округления, а переупорядочивание узлов для получения выпуклого патча составляет почти половину оставшихся строк. Аккуратная часть с подключенными компонентами довольно короткая. 06.03.2018networkx
иsympy
. Я почти уверен, что это слишком тяжелые зависимости для команды matplotlib, чтобы рассмотреть вопрос о включении в основную ветку matplotlib. Заменитьis_coplanar
и, таким образом, избавиться от библиотеки sympy, вероятно, довольно просто (я просто не доверял своей математике в тот вечер); Однако реализация структуры графа и поиск связанных компонентов — это как минимум страница кода. 06.03.2018reorder
все еще нужно сделать надежным, и в данный момент у меня нет хорошей идеи. Любые предложения, @ImportanceOfBeingErnest (пока мы привлекаем ваше внимание)? 06.03.2018reorder
вообще работает так, как есть, поскольку это кажется довольно произвольным. Когда я думал об этом вчера, у меня возникла идея использовать еще одинconvexhull
, потому что каждая грань лежит в 2D-плоскости, а выпуклая оболочка дает упорядоченный набор точек в 2D. Но для этого требуется проецировать вперед и назад между 3D-точками и 2D-точками. Это то, что вы имели ввиду? 06.03.2018plot_trisurf
, где соединяются копланарные треугольники (это также решило бы этот случай здесь, но было бы применимо ко многим другим случаям использования). Я сомневаюсь, что кто-то приложит к этому какие-либо усилия, matplotlib 3D на данный момент довольно мертв, но если есть запрос на его реализацию, это может приветствоваться. 06.03.2018plot_trisurf
. 06.03.2018