MayaVi contour3d after coordinate transformation

Question

I have a 3D scalar field mesh given in non-cartesian coordinate system. After coordinate transformation back to conventional cartesian coordinates mlab.contour3d displays wrong result, while mlab.points3d works as expected. How can I view isosurfaces of given mesh in different coordinate systems?

This is my code

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

After watching 3D Visualization with Mayavi I managed to solve it myself. Problem is that mlab.contour3d works only with rectilinear grid data (grid generated with np.meshgrid or np.mgrid). In my case one can use tvtk.StructuredGrid object for structured grids which have same topology as rectilinear grid but nonuniform spacing and directions between points.

This is working code:

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