Trazar 3D regiones cerradas convexas en matplot lib

Estoy tratando de trazar en 3D un politopo definido por un conjunto de desigualdades. Esencialmente, trato de reproducir la funcionalidad de esta biblioteca de trazados gráficos de matlab en matplotlib.

Mi enfoque es obtener los vértices de intersección, construir el casco convexo de ellos, y luego obtener y trazar las caras resultantes (simplices).

El problema es que muchas simplicidades son coplanares y están haciendo que la ttwig esté muy ocupada sin ninguna razón (vea todos estos bordes diagonales en la ttwig a continuación).

¿Hay alguna manera fácil de imprimir los bordes “externos” del poliedro, sin tener que consolidarme uno a uno, todas las sencillas coplanares?

Gracias

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

Resultado del código anterior.

Bastante seguro de que no hay nada nativo en matplotlib. Sin embargo, encontrar las caras que pertenecen juntas no es particularmente difícil. La idea básica implementada a continuación es crear un gráfico, donde cada nodo es un triángulo. Luego conectas triangularjs que son coplanarios y adyacentes. Finalmente, encuentra los componentes conectados de la gráfica para determinar qué triangularjs forman una cara.

introduzca la descripción de la imagen aquí

 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 # ie 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) 

Aplicado a tu ejemplo:

 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) 

Nota

Después de reflexionar, la función reordered probablemente necesita un poco más de trabajo. Bastante seguro de que esto se romperá para formas extrañas / no convexas, y ni siquiera estoy 100% seguro de que siempre funcionará para formas convexas. El descanso debe estar bien sin embargo.

La siguiente sería mi versión de una solución. Es similar a la solución de @ Paul, ya que toma los triangularjs, los agrupa por cara a la que pertenecen y los une en una sola cara.

La diferencia sería principalmente que esta solución no utiliza nx o simpy . Muchas de las operaciones necesarias se realizan mediante reindexación, uso extensivo de álgebra unique y un poco lineal.
El orden de los vértices de las caras finales está determinado por ConvexHull . Creo que esto no debería ser una limitación, ya que (creo que) cualquier intersección de medio espacio debería dar como resultado formas convexas solamente. Sin embargo, también agregué otro método que se puede utilizar si las formas no son convexas (según la idea de esta pregunta ).

 from scipy.spatial import HalfspaceIntersection from scipy.spatial import ConvexHull import numpy as np import matplotlib.pyplot as plt import mpl_toolkits.mplot3d as a3 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) simplices = hull.simplices org_triangles = [verts[s] for s in simplices] class Faces(): def __init__(self,tri, sig_dig=12, method="convexhull"): self.method=method self.tri = np.around(np.array(tri), sig_dig) self.grpinx = list(range(len(tri))) norms = np.around([self.norm(s) for s in self.tri], sig_dig) _, self.inv = np.unique(norms,return_inverse=True, axis=0) def norm(self,sq): cr = np.cross(sq[2]-sq[0],sq[1]-sq[0]) return np.abs(cr/np.linalg.norm(cr)) def isneighbor(self, tr1,tr2): a = np.concatenate((tr1,tr2), axis=0) return len(a) == len(np.unique(a, axis=0))+2 def order(self, v): if len(v) <= 3: return v v = np.unique(v, axis=0) n = self.norm(v[:3]) y = np.cross(n,v[1]-v[0]) y = y/np.linalg.norm(y) c = np.dot(v, np.c_[v[1]-v[0],y]) if self.method == "convexhull": h = ConvexHull(c) return v[h.vertices] else: mean = np.mean(c,axis=0) d = c-mean s = np.arctan2(d[:,0], d[:,1]) return v[np.argsort(s)] def simplify(self): for i, tri1 in enumerate(self.tri): for j,tri2 in enumerate(self.tri): if j > i: if self.isneighbor(tri1,tri2) and \ self.inv[i]==self.inv[j]: self.grpinx[j] = self.grpinx[i] groups = [] for i in np.unique(self.grpinx): u = self.tri[self.grpinx == i] u = np.concatenate([d for d in u]) u = self.order(u) groups.append(u) return groups f = Faces(org_triangles) g = f.simplify() ax = a3.Axes3D(plt.figure()) colors = list(map("C{}".format, range(len(g)))) pc = a3.art3d.Poly3DCollection(g, facecolor=colors, edgecolor="k", alpha=0.9) ax.add_collection3d(pc) ax.dist=10 ax.azim=30 ax.elev=10 ax.set_xlim([0,5]) ax.set_ylim([0,5]) ax.set_zlim([0,5]) plt.show() 

introduzca la descripción de la imagen aquí