Mapea una imagen en una esfera y traza trayectorias 3D

Lo que me gustaría hacer es definir una esfera en el centro de mi sistema de coordenadas 3D (con radio = 1), envolver un mapa planetario cilíndrico en la superficie de la esfera (es decir, realizar un mapeo de texturas en la esfera) y trazar trayectorias 3D alrededor del Objeto (como trayectorias de satélites). ¿Hay alguna manera de hacer esto usando matplotlib o mayavi?

Trazar las trayectorias es fácil usando mayavi.mlab.plot3d una vez que tenga su planeta, así que me concentraré en el mapeado de texturas de un planeta a una esfera usando mayavi. (En principio, podemos realizar la tarea utilizando matplotlib, pero el rendimiento y la calidad son mucho peores en comparación con mayavi; consulte el final de esta respuesta).

El buen escenario: esfera sobre una esfera.

Resulta que si desea mapear una imagen parametrizada de forma esférica en una esfera, tiene que ensuciarse un poco la mano y usar algo de vtk, pero en realidad hay muy poco trabajo por hacer y el resultado se ve excelente. Voy a usar una imagen de Blue Marble de la NASA para la demostración. Su readme dice que estas imágenes tienen

una proyección geográfica (Plate Carrée), que se basa en un espaciado de cuadrícula de latitud-longitud igual (¡no una proyección de área igual!)

Al buscarlo en wikipedia resulta que esto también se conoce como una proyección equirectangular . En otras palabras, los píxeles a lo largo de x corresponden directamente a la longitud y los píxeles a lo largo de y corresponden directamente a la latitud. Esto es lo que yo llamo “parametrizado esféricamente”.

Entonces, en este caso, podemos usar un TexturedSphereSource bajo nivel para generar una esfera en la que se pueda mapear la textura. Construir una malla de esferas nosotros mismos podría llevar a artefactos en el mapeo (más sobre esto más adelante).

Para el trabajo de vtk de bajo nivel reelaboré el ejemplo oficial que se encuentra aquí . Aquí está todo lo que se necesita:

 from mayavi import mlab from tvtk.api import tvtk # python wrappers for the C++ vtk ecosystem def auto_sphere(image_file): # create a figure window (and scene) fig = mlab.figure(size=(600, 600)) # load and map the texture img = tvtk.JPEGReader() img.file_name = image_file texture = tvtk.Texture(input_connection=img.output_port, interpolate=1) # (interpolate for a less raster appearance when zoomed in) # use a TexturedSphereSource, aka getting our hands dirty R = 1 Nrad = 180 # create the sphere source with a given radius and angular resolution sphere = tvtk.TexturedSphereSource(radius=R, theta_resolution=Nrad, phi_resolution=Nrad) # assemble rest of the pipeline, assign texture sphere_mapper = tvtk.PolyDataMapper(input_connection=sphere.output_port) sphere_actor = tvtk.Actor(mapper=sphere_mapper, texture=texture) fig.scene.add_actor(sphere_actor) if __name__ == "__main__": image_file = 'blue_marble_spherical.jpg' auto_sphere(image_file) mlab.show() 

El resultado es exactamente lo que esperamos:

mármol azul mapeado en una esfera

El escenario menos agradable: no una esfera.

Desafortunadamente, no pude averiguar cómo usar un mapeo no esférico con el método anterior. Además, puede suceder que no queramos mapear en una esfera perfecta, sino en un elipsoide o un objeto redondo similar. Para este caso, es posible que tengamos que construir la superficie nosotros mismos y probar el mapeado de texturas en esa superficie. Alerta de spoiler: no será tan bonita.

Partiendo de una esfera generada manualmente, podemos cargar la textura de la misma manera que antes y trabajar con objetos de nivel superior construidos por mlab.mesh :

 import numpy as np from mayavi import mlab from tvtk.api import tvtk import matplotlib.pyplot as plt # only for manipulating the input image def manual_sphere(image_file): # caveat 1: flip the input image along its first axis img = plt.imread(image_file) # shape (N,M,3), flip along first dim outfile = image_file.replace('.jpg', '_flipped.jpg') # flip output along first dim to get right chirality of the mapping img = img[::-1,...] plt.imsave(outfile, img) image_file = outfile # work with the flipped file from now on # parameters for the sphere R = 1 # radius of the sphere Nrad = 180 # points along theta and phi phi = np.linspace(0, 2 * np.pi, Nrad) # shape (Nrad,) theta = np.linspace(0, np.pi, Nrad) # shape (Nrad,) phigrid,thetagrid = np.meshgrid(phi, theta) # shapes (Nrad, Nrad) # compute actual points on the sphere x = R * np.sin(thetagrid) * np.cos(phigrid) y = R * np.sin(thetagrid) * np.sin(phigrid) z = R * np.cos(thetagrid) # create figure mlab.figure(size=(600, 600)) # create meshed sphere mesh = mlab.mesh(x,y,z) mesh.actor.actor.mapper.scalar_visibility = False mesh.actor.enable_texture = True # probably redundant assigning the texture later # load the (flipped) image for texturing img = tvtk.JPEGReader(file_name=image_file) texture = tvtk.Texture(input_connection=img.output_port, interpolate=0, repeat=0) mesh.actor.actor.texture = texture # tell mayavi that the mapping from points to pixels happens via a sphere mesh.actor.tcoord_generator_mode = 'sphere' # map is already given for a spherical mapping cylinder_mapper = mesh.actor.tcoord_generator # caveat 2: if prevent_seam is 1 (default), half the image is used to map half the sphere cylinder_mapper.prevent_seam = 0 # use 360 degrees, might cause seam but no fake data #cylinder_mapper.center = np.array([0,0,0]) # set non-trivial center for the mapping sphere if necessary 

Como puede ver los comentarios en el código, hay algunas advertencias. La primera es que el modo de mapeo esférico invierte la imagen de entrada por alguna razón (lo que lleva a una Tierra reflejada). Entonces, utilizando este método, primero tenemos que crear una versión invertida de la imagen de entrada. Esto solo debe hacerse una vez por imagen, pero dejé el bloque de código correspondiente en la parte superior de la función anterior.

La segunda advertencia es que si el atributo prevent_seam del prevent_seam de texturas se deja en el valor predeterminado 1, el mapeo pasa de 0 a 180 azimut y la otra mitad de la esfera obtiene un mapeo reflejado. Claramente no queremos esto: queremos mapear la esfera completa de 0 a 360 azimut. A medida que sucede, esta asignación podría implicar que veamos una costura (discontinuidad) en la asignación en phi=0 , es decir, en el borde del mapa. Esta es otra razón por la que se debe usar el primer método cuando sea posible. De todos modos, aquí está el resultado, que contiene el punto phi=0 (lo que demuestra que no hay costura):

versión en malla

Mapeo cilíndrico

La forma en que funcionan los mapeos esféricos anteriores es que cada punto en la superficie se proyecta en una esfera a través de un punto dado en el espacio. Para el primer ejemplo, este punto es el origen, para el segundo caso podemos establecer una matriz de 3 longitudes como el valor de cylinder_mapper.center para mapear en esferas no centradas en el origen.

Ahora, tu pregunta menciona un mapeo cilíndrico. En principio podemos hacerlo utilizando el segundo método:

 mesh.actor.tcoord_generator_mode = 'cylinder' cylinder_mapper = mesh.actor.tcoord_generator cylinder_mapper.automatic_cylinder_generation = 0 # use manual cylinder from points cylinder_mapper.point1 = np.array([0,0,-R]) cylinder_mapper.point2 = np.array([0,0,R]) cylinder_mapper.prevent_seam = 0 # use 360 degrees, causes seam but no fake data 

Esto cambiará la cartografía esférica a cilíndrica. Define la proyección en términos de los dos puntos ( [0,0,-R] y [0,0,R] ) que establecen el eje y la extensión del cilindro. Cada punto se asigna de acuerdo con sus coordenadas cilíndricas (phi,z) : el azimut de 0 a 360 grados y la proyección vertical de la coordenada. Los comentarios anteriores relativos a la costura todavía se aplican.

Sin embargo, si tuviera que hacer una asignación tan cilíndrica, definitivamente intentaría usar el primer método. En el peor de los casos, esto significa que tenemos que transformar el mapa parametrizado de forma cilíndrica en uno parametrizado esféricamente. Nuevamente, esto solo debe hacerse una vez por mapa, y se puede hacer fácilmente utilizando la interpolación 2d, por ejemplo, utilizando scipy.interpolate.RegularGridInterpolator . Para la transformación específica, debe conocer los detalles específicos de su proyección no esférica, pero no debería ser demasiado difícil transformarla en una proyección esférica, que luego puede usar de acuerdo con el primer caso de TexturedSphereSource .

Apéndice: matplotlib

Para completar, puedes hacer lo que quieras con matplotlib, pero necesitarás mucha más memoria y CPU (y debes tener en cuenta que debes usar mayavi o matplotlib, pero no puedes combinar ambos en una figura). La idea es definir una malla que corresponda a los píxeles del mapa de entrada y pasar la imagen como el argumento de palabras clave de Axes3D.plot_surface de Axes3D.plot_surface . La construcción es tal que la resolución de la esfera se acopla directamente a la resolución del mapeo. Solo podemos usar una pequeña cantidad de puntos para mantener la necesidad de memoria, pero luego el resultado se verá mal pixelado. De todas formas:

 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D def mpl_sphere(image_file): img = plt.imread(image_file) # define a grid matching the map size, subsample along with pixels theta = np.linspace(0, np.pi, img.shape[0]) phi = np.linspace(0, 2*np.pi, img.shape[1]) count = 180 # keep 180 points along theta and phi theta_inds = np.linspace(0, img.shape[0] - 1, count).round().astype(int) phi_inds = np.linspace(0, img.shape[1] - 1, count).round().astype(int) theta = theta[theta_inds] phi = phi[phi_inds] img = img[np.ix_(theta_inds, phi_inds)] theta,phi = np.meshgrid(theta, phi) R = 1 # sphere x = R * np.sin(theta) * np.cos(phi) y = R * np.sin(theta) * np.sin(phi) z = R * np.cos(theta) # create 3d Axes fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_surface(xT, yT, zT, facecolors=img/255, cstride=1, rstride=1) # we've already pruned ourselves # make the plot more spherical ax.axis('scaled') if __name__ == "__main__": image_file = 'blue_marble.jpg' mpl_sphere(image_file) plt.show() 

El parámetro de count en lo anterior es lo que define el submuestreo del mapa y el tamaño correspondiente de la esfera renderizada. Con los 180 ajustes anteriores obtenemos la siguiente figura:

versión matplotlib

Además, matplotlib usa un renderizador 2d, lo que implica que para la representación complicada de objetos 3d a menudo se producen artefactos extraños (en particular, los objetos extendidos pueden estar completamente uno frente al otro o uno detrás del otro, por lo que las geometrías entrelazadas generalmente parecen rotos). Teniendo en cuenta esto, definitivamente usaría mayavi para trazar una esfera texturizada. (Aunque el mapeo en el caso de matplotlib funciona cara a cara en la superficie, por lo que se puede aplicar directamente a superficies arbitrarias).