将图像映射到球体上并绘制3D轨迹

Question

我想要做的是在我的3D坐标系的中心定义一个球体(半径= 1),将圆柱形行星地图包裹在球体表面上(即在球体上执行纹理映射)并绘制围绕该球体的3D轨迹.物体(如卫星轨迹).有什么方法可以使用matplotlib或mayavi来做到这一点吗？

Answer 1

mayavi.mlab.plot3d一旦拥有了你的星球,绘制轨迹很容易使用,所以我将专注于使用mayavi将行星映射到球体的纹理.(原则上我们可以使用matplotlib执行任务,但与mayavi相比,性能和质量要差得多,请参阅本答案的结尾.)

漂亮的场景:球体上的球体

事实证明,如果你想将一个球形的参数化图像映射到一个球体上,你必须让你的手有点脏,并使用一些裸的vtk,但实际上很少有工作要做,结果看起来很棒.我将使用美国宇航局的蓝色大理石图像进行演示.他们的自述文件说这些图像有

地理(PlateCarrée)投影,基于相等的纬度 - 经度网格间距(不是等面积投影!)

在维基百科中查找它,结果证明这也称为equirectangular投影.换句话说,沿x直接对应于经度的像素和沿y直接对应于纬度的像素.这就是我所说的"球形参数化".

因此,在这种情况下,我们可以使用低级别TexturedSphereSource来生成可以将纹理映射到的球体.自己构造一个球体网格可能会导致映射中的工件(稍后会详细介绍).

对于低级别的vtk工作,我重写了这里发现的官方示例.这就是它所需要的一切:

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, a.k.a. 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()

结果正是我们所期望的:

不太好的场景:不是球体

不幸的是我无法弄清楚如何使用上述方法的非球面映射.此外,我们可能不希望在完美的球体上进行映射,而是在椭圆体或类似的圆形物体上进行映射.对于这种情况,我们可能必须自己构建表面并尝试在其上进行纹理映射.剧透警报:它不会那么漂亮.

从手动生成的球体开始,我们可以像以前一样加载纹理,并使用由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

正如您可以在代码中看到的那样,有一些注意事项.第一个是球形映射模式由于某种原因翻转输入图像(导致反射的地球).因此,使用此方法,我们首先必须创建输入图像的翻转版本.这只需要对每个图像执行一次,但是我将相应的代码块留在上面函数的顶部.

第二个警告是,如果prevent_seam纹理映射器的属性保留在默认的1值上,则映射从0到180方位角发生,而另一半的球体获得反射映射.我们显然不希望这样:我们想要将整个球体从0​​方位映射到360方位角.碰巧的是,这种映射可能意味着我们在映射中看到了一个接缝(不连续)phi=0,即在地图的边缘.这是在可能的情况下使用第一种方法的另一个原因.无论如何,这是结果,包含phi=0点(证明没有接缝):

圆柱映射

上述球形映射的工作方式是通过空间中的给定点将表面上的每个点投影到球体上.对于第一个示例,此点是原点,对于第二种情况,我们可以将3长度数组设置为值,cylinder_mapper.center以便映射到非原点居中的球体.

现在,您的问题提到了圆柱映射.原则上我们可以使用第二种方法:

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

这会将球面映射更改为圆柱形.它根据设置圆柱轴和范围的两个点([0,0,-R]和[0,0,R])来定义投影.每个点都根据其圆柱坐标进行映射(phi,z):0到360度的方位角和坐标的垂直投影.关于接缝的早期评论仍然适用.

但是,如果我不得不做这样的圆柱映射,我肯定会尝试使用第一种方法.在最坏的情况下,这意味着我们必须将圆柱参数化地图转换为球形参数化地图.同样,这只需要在每个地图上完成一次,并且可以使用2d插值轻松完成,例如通过使用scipy.interpolate.RegularGridInterpolator.对于特定的变换,您必须知道非球面投影的细节,但是将其转换为球形投影应该不会太难,您可以根据第一种情况使用它TexturedSphereSource.

附录:matplotlib

为了完整起见,你可以使用matplotlib做你想做的事情,但它需要更多的内存和CPU(并注意你必须使用mayavi或matplotlib,但你不能在图中混合两者).我们的想法是定义一个与输入映射的像素相对应的网格,并将图像作为facecolors关键字参数传递Axes3D.plot_surface.该结构使得球体的分辨率直接耦合到映射的分辨率.我们只能使用少量的点来保持内存需要易处理,但结果看起来很像素化.无论如何:

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(x.T, y.T, z.T, 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()

count上面的参数定义了地图的下采样和渲染球体的相应大小.使用上面的180设置,我们得到下图:

此外,matplotlib使用2D渲染器,这意味着对于复杂的3d对象渲染通常最终会产生奇怪的伪像(特别是,扩展对象可能完全在彼此的前面或后面,因此互锁几何通常看起来很破碎).考虑到这些,我肯定会使用mayavi来绘制纹理球体.(尽管matplotlib案例中的映射在表面上是面对面的,因此它可以直接应用于任意表面.)