坐标变换后的MayaVi contour3d

Question

我有一个在非笛卡尔坐标系中给出的 3D 标量场网格。坐标变换回常规笛卡尔坐标后 mlab.contour3d 显示错误结果，而 mlab.points3d 按预期工作。如何在不同坐标系中查看给定网格的等值面？

这是我的代码

import os
import numpy as np
# fix incorrect order in foregroud objects
os.environ['ETS_TOOLKIT'] = 'qt4'
os.environ['QT_API'] = 'pyqt'
from mayavi import mlab


def plot_cell(cell, mlab):
    for nr, i in enumerate(cell):
        coord = np.zeros((4, 3), dtype=float)
        coord[1] = i
        for nr2, j in enumerate(cell):
            if nr == nr2:
                continue

            coord[2] = i + j
            for nr3, k in enumerate(cell):
                if nr3 == nr or nr3 == nr2:
                    continue
                coord[3] = i + j + k
                mlab.plot3d(*coord.T, color=(0, 0, 0), line_width=0.01)

# generate data in non-cartesian coordinate system
scaled_coord = [np.linspace(0, 1, 20, endpoint=False) for i in range(3)]
XYZ = np.array(np.meshgrid(*scaled_coord, indexing="ij"))
data = np.sin(2*np.pi*XYZ[0])*np.sin(2*np.pi*XYZ[1])*np.sin(2*np.pi*XYZ[2])

plot_cell(np.eye(3), mlab)
mlab.contour3d(*XYZ, data)
mlab.savefig("old_coord.png")
mlab.close()

# transform to cartesian coordinates
cell = np.array(
    [[ 1.  ,  0.  ,  0.  ],
     [-0.5 ,  0.87,  0.  ],
     [ 0.  ,  0.  ,  3.07]])
transformation_matrix = cell.T
XYZ2 = np.einsum('ij,jabc->iabc', transformation_matrix, XYZ)

# plot transformed result
plot_cell(cell, mlab)
mlab.contour3d(*XYZ2, data)
mlab.savefig("new_coord.png")
mlab.close()

# plot points
plot_cell(cell, mlab)
mlab.points3d(*XYZ2, data)
mlab.savefig("new_coord_expected.png")
mlab.close()

Answer 1

看完3D Visualization with Mayavi我自己解决了。问题是 mlab.contour3d 仅适用于直线网格数据（使用 np.meshgrid 或 np.mgrid 生成的网格）。在我的例子中，可以将 tvtk.StructuredGrid 对象用于具有相同拓扑结构的结构化网格直线网格但点之间的间距和方向不均匀。

这是工作代码：

import os
import numpy as np
from tvtk.api import tvtk
# fix incorrect order in foregroud objects
os.environ['ETS_TOOLKIT'] = 'qt4'
os.environ['QT_API'] = 'pyqt'
from mayavi import mlab


def plot_cell(cell, mlab):
    for nr, i in enumerate(cell):
        coord = np.zeros((4, 3), dtype=float)
        coord[1] = i
        for nr2, j in enumerate(cell):
            if nr == nr2:
                continue

            coord[2] = i + j
            for nr3, k in enumerate(cell):
                if nr3 == nr or nr3 == nr2:
                    continue
                coord[3] = i + j + k
                mlab.plot3d(*coord.T, color=(0, 0, 0), line_width=0.01)


scaled_coord = [np.linspace(0, 1, 20, endpoint=False) for i in range(3)]
ABC = np.array(np.meshgrid(*scaled_coord, indexing="ij"))
data = np.sin(2*np.pi*ABC[0])*np.sin(2*np.pi*ABC[1])*np.sin(2*np.pi*ABC[2])

# transform to cartesian coordinates
cell = np.array(
    [[ 1.  ,  0.  ,  0.  ],
     [-0.5 ,  0.87,  0.  ],
     [ 0.  ,  0.  ,  3.07]])
transformation_matrix = cell.T
x, y, z = np.einsum('ij,jabc->iabc', transformation_matrix, ABC)

def generate_structured_grid(x, y, z, scalars):
    pts = np.empty(z.shape + (3,), dtype=float)
    pts[..., 0] = x
    pts[..., 1] = y
    pts[..., 2] = z

    pts = pts.transpose(2, 1, 0, 3).copy()
    pts.shape = int(pts.size / 3), 3
    scalars = scalars.T.copy()

    sg = tvtk.StructuredGrid(dimensions=x.shape, points=pts)
    sg.point_data.scalars = scalars.ravel()
    sg.point_data.scalars.name = 'scalars'
    return sg

sgrid = generate_structured_grid(x, y, z, data)
src = mlab.pipeline.add_dataset(sgrid)
isosurface = mlab.pipeline.iso_surface(src)
plot_cell(cell, mlab)
mlab.show()

坐标变换后的MayaVi contour3d

MayaVi contour3d after coordinate transformation

python

plot

mayavi