将图像映射到球体上并绘制3D轨迹 [英] Map an image onto a sphere and plot 3D trajectories

查看:400
本文介绍了将图像映射到球体上并绘制3D轨迹的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

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

解决方案

使用来自NASA的蓝色大理石图像进行演示.他们的自述说这些图片有

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

在Wikipedia中进行查找,结果发现这也称为 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(p 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])定义投影.每个点均根据其圆柱坐标(pz)进行映射:0到360度的方位角以及坐标的垂直投影.先前有关接缝的说明仍然适用.

但是,如果我必须进行这样的圆柱映射,那么我肯定会尝试使用第一种方法.在最坏的情况下,这意味着我们必须将圆柱参数化的贴图转换为球面参数化的贴图.同样,每个地图只需要做一次,并且可以使用2d插值轻松完成,例如,使用

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

此外,matplotlib使用2d渲染器,这意味着对于复杂的3d对象,渲染通常会以怪异的伪像结束(特别是,扩展的对象可以完全位于彼此之间,也可以彼此之间互斥,因此互锁的几何形状通常看起来很破碎). .考虑到这些,我肯定会使用mayavi绘制纹理球体. (尽管matplotlib情况下的映射在表面上并排工作,因此可以直接应用于任意表面.)

What I would like to do is to define a sphere in the center of my 3D coordinate system (with radius=1), wrap a cylindrical planet map onto the sphere's surface (i.e. perform texture mapping on the sphere) and plot 3D trajectories around the object (like satellite trajectories). Is there any way I can do this using matplotlib or mayavi?

解决方案

Plotting trajectories is easy using mayavi.mlab.plot3d once you have your planet, so I'm going to concentrate on texture mapping a planet to a sphere using mayavi. (In principle we can perform the task using matplotlib, but the performance and quality is much worse compared to mayavi, see the end of this answer.)

The nice scenario: sphere on a sphere

It turns out that if you want to map a spherically parametrized image onto a sphere you have to get your hand a little dirty and use some bare vtk, but there's actually very little work to be done and the result looks great. I'm going to use a Blue Marble image from NASA for demonstration. Their readme says that these images have

a geographic (Plate Carrée) projection, which is based on an equal latitude- longitude grid spacing (not an equal area projection!)

Looking it up in wikipedia it turns out that this is also known as an equirectangular projection. In other words, pixels along x directly correspond to longitude and pixels along y directly correspond to latitude. This is what I call "spherically parametrized".

So in this case we can use a low-level TexturedSphereSource in order to generate a sphere which the texture can be mapped onto. Constructing a sphere mesh ourselves could lead to artifacts in the mapping (more on this later).

For the low-level vtk work I reworked the official example found here. Here's all that it takes:

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

The result is exactly what we expect:

The less nice scenario: not a sphere

Unfortunately I couldn't figure out how to use a non-spherical mapping with the above method. Furthermore, it might happen that we don't want to map on a perfect sphere, but rather on an ellipsoid or similar round object. For this case we may have to construct the surface ourselves and try texture-mapping onto that. Spoiler alert: it won't be as pretty.

Starting from a manually generated sphere we can load the texture in the same way as before, and work with higher-level objects constructed by 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

As you can see the comments in the code, there are a few caveats. The first one is that the spherical mapping mode flips the input image for some reason (which leads to a reflected Earth). So using this method we first have to create a flipped version of the input image. This only has to be done once per image, but I left the corresponding code block at the top of the above function.

The second caveat is that if the prevent_seam attribute of the texture mapper is left on the default 1 value, the mapping happens from 0 to 180 azimuth and the other half of the sphere gets a reflected mapping. We clearly don't want this: we want to map the entire sphere from 0 to 360 azimuth. As it happens, this mapping might imply that we see a seam (discontinuity) in the mapping at phi=0, i.e. at the edge of the map. This is another reason why using the first method should be used when possible. Anyway, here is the result, containing the phi=0 point (demonstrating that there's no seam):

Cylindrical mapping

The way the above spherical mappings work is that each point on the surface gets projected onto a sphere via a given point in space. For the first example this point is the origin, for the second case we can set a 3-length array as the value of cylinder_mapper.center in order to map onto non-origin-centered spheres.

Now, your question mentions a cylindrical mapping. In principle we can do this using the second method:

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

This will change spherical mapping to cylindrical. It defines the projection in terms of the two points ([0,0,-R] and [0,0,R]) that set the axis and extent of the cylinder. Every point is mapped according to its cylindrical coordinates (phi,z): the azimuth from 0 to 360 degrees and the vertical projection of the coordinate. Earlier remarks concerning the seam still apply.

However, if I had to do such a cylindrical mapping, I'd definitely try to use the first method. At worst this means we have to transform the cylindrically parametrized map to a spherically parametrized one. Again this only has to be done once per map, and can be done easily using 2d interpolation, for instance by using scipy.interpolate.RegularGridInterpolator. For the specific transformation you have to know the specifics of your non-spherical projection, but it shouldn't be too hard to transform it into a spherical projection, which you can then use according to the first case with TexturedSphereSource.

Appendix: matplotlib

For the sake of completeness, you can do what you want using matplotlib but it will take a lot more memory and CPU (and note that you have to use either mayavi or matplotlib, but you can't mix both in a figure). The idea is to define a mesh that corresponds to pixels of the input map, and pass the image as the facecolors keyword argument of Axes3D.plot_surface. The construction is such that the resolution of the sphere is directly coupled to the resolution of the mapping. We can only use a small number of points to keep the memory need tractable, but then the result will look badly pixellated. Anyway:

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

The count parameter in the above is what defines the downsampling of the map and the corresponding size of the rendered sphere. With the above 180 setting we get the following figure:

Furthermore, matplotlib uses a 2d renderer, which implies that for complicated 3d objects rendering often ends up with weird artifacts (in particular, extended objects can either be completely in front of or behind one another, so interlocking geometries usually look broken). Considering these I would definitely use mayavi for plotting a textured sphere. (Although the mapping in the matplotlib case works face-by-face on the surface, so it can be straightforwardly applied to arbitrary surfaces.)

这篇关于将图像映射到球体上并绘制3D轨迹的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆