instalar - Delaunay ¿Triangulación de puntos de superficie 2D en 3D con python?

python matplotlib line plot

Question

Tengo una colección de puntos 3D. Estos puntos se muestrean en niveles constantes (z = 0,1, ..., 7). Una imagen debe dejarlo claro:

Estos puntos están en una gran cantidad de formas (N, 3) llamadas X La trama anterior se crea utilizando:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

X = load(''points.npy'')
fig = plt.figure()
ax = fig.gca(projection=''3d'')
ax.plot_wireframe(X[:,0], X[:,1], X[:,2])
ax.scatter(X[:,0], X[:,1], X[:,2])
plt.draw()

En cambio, me gustaría triangular solo la superficie de este objeto y trazar la superficie. Sin embargo, no quiero el casco convexo de este objeto porque pierde información de forma sutil que me gustaría poder inspeccionar.

He probado ax.plot_trisurf(X[:,0], X[:,1], X[:,2]) , pero esto resulta en el siguiente lío:

¿Alguna ayuda?

Ejemplo de datos

Aquí hay un fragmento para generar datos 3D que son representativos del problema:

import numpy as np
X = []
for i in range(8):
    t = np.linspace(0,2*np.pi,np.random.randint(30,50))
    for j in range(t.shape[0]):
        # random circular objects...
        X.append([
            (-0.05*(i-3.5)**2+1)*np.cos(t[j])+0.1*np.random.rand()-0.05,
            (-0.05*(i-3.5)**2+1)*np.sin(t[j])+0.1*np.random.rand()-0.05,
            i
        ])
X = np.array(X)

Ejemplo de datos de la imagen original

Aquí hay un pastebin a los datos originales:

http://pastebin.com/YBZhJcsV

Aquí están las rebanadas a lo largo de la constante z:

Answer 1

actualización 2

Ahora hago esto de la siguiente manera:

1. Utilizo el hecho de que las rutas en cada corte en z son cerradas y simples, y uso matplotlib.path para determinar los puntos dentro y fuera del contorno. Usando esta idea, convierto los contornos de cada sector en una imagen de valor booleano, que se combina en un volumen de valor booleano.
2. A continuación, utilizo marching_cubes método marching_cubes skimage para obtener una triangulación de la superficie para la visualización.

Aquí hay un ejemplo del método. Creo que los datos son ligeramente diferentes, pero definitivamente puedes ver que los resultados son mucho más limpios y pueden manejar superficies que están desconectadas o tienen agujeros.

actualización 1 (todavía mal)

Debería actualizar esto para las futuras personas que lo encuentren. Si bien el método anterior funciona la mayor parte del tiempo, asume (a través de la transformación de coordenadas esféricas) que no hay dos puntos a lo largo del mismo rayo. Si observa el artefacto en la parte central izquierda de la imagen de arriba, esta es la razón.

Un mejor enfoque es hacer "cirugía" en la superficie. Pensando en la superficie como una cáscara de naranja, la cortas por un lado y luego la abres, estirándola. Luego tienes un plano 2D que puedes triangular e interpolar. Simplemente debe realizar un seguimiento de cómo volver al lugar correspondiente en 3D. Para implementar esta idea se necesita bastante trabajo, y la implementación también requiere un cuidado especial exclusivo de cómo se representan mis datos.

De todos modos, esto es solo para dar una indicación de cómo se podría abordar esto de manera más sólida.

Respuesta original

Ok, aquí está la solución que se me ocurrió. Depende en gran medida de que mis datos sean aproximadamente esféricos y se muestreen de manera uniforme en z, creo. Algunos de los otros comentarios proporcionan más información sobre soluciones más robustas. Como mis datos son aproximadamente esféricos, triangulo los ángulos de azimut y cenit de la transformación de coordenadas esféricas de mis puntos de datos.

import numpy as np
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri

X = np.load(''./mydatars.npy'')
# My data points are strictly positive. This doesn''t work if I don''t center about the origin.
X -= X.mean(axis=0)

rad = np.linalg.norm(X, axis=1)
zen = np.arccos(X[:,-1] / rad)
azi = np.arctan2(X[:,1], X[:,0])

tris = mtri.Triangulation(zen, azi)

fig = plt.figure()
ax  = fig.add_subplot(111, projection=''3d'')
ax.plot_trisurf(X[:,0], X[:,1], X[:,2], triangles=tris.triangles, cmap=plt.cm.bone)
plt.show()

Usando los datos de muestra del pastebin anterior, esto produce:

Answer 2

Me doy cuenta de que mencionaste en tu pregunta que no querías usar el casco convexo porque podrías perder alguna información de forma. Tengo una solución simple que funciona bastante bien para sus datos de ejemplo ''esféricos tontas'', aunque usa scipy.spatial.ConvexHull . Pensé que lo compartiría aquí de todos modos, en caso de que sea útil para otros:

from matplotlib.tri import triangulation
from scipy.spatial import ConvexHull

# compute the convex hull of the points
cvx = ConvexHull(X)

x, y, z = X.T

# cvx.simplices contains an (nfacets, 3) array specifying the indices of
# the vertices for each simplical facet
tri = Triangulation(x, y, triangles=cvx.simplices)

fig = plt.figure()
ax = fig.gca(projection=''3d'')
ax.hold(True)
ax.plot_trisurf(tri, z)
ax.plot_wireframe(x, y, z, color=''r'')
ax.scatter(x, y, z, color=''r'')

plt.draw()

Lo hace bastante bien en este caso, ya que los datos de su ejemplo terminan en una superficie más o menos convexa. ¿Quizás podrías hacer algunos ejemplos de datos más desafiantes? Una superficie toroidal sería un buen caso de prueba en el que el método del casco convexo obviamente fallaría.

Mapear una superficie 3D arbitraria desde una nube de puntos es un problema realmente difícil. Aquí hay una pregunta relacionada que contiene algunos enlaces que pueden ser útiles.