Scalar filed is not volume rendered with a colormap

Hi, all

I am a graduate student using mayavi for my research. I am trying to volome render a scalar field stored as a .csv file. The shape of the field is rectangular, halved into two parts: the upper part with a scalar of 0.25 and the lower part with a scalar of 1.0. Since I am using the cool colormap, I expected that the upper and lower parts should be rendered a translucent sky-blue volume and a translucent pink volume, respectively, like this example. However, as shown in the left figure attached below, the field is rendered mainly in translucent black, while the scalar field data itself is exactly what I expected (the upper and lower parts have a uniform scalar, respectively), as shown in the right plot with matplotlib. Actually, the surfaces of the parts are painted somewhat sky-blue and pink, but the inside of the parts is black. The shape is exactly what I expected. Does someone come up with why my scalar field is volume rendered like that? And, can anybody give me any advice to resolve it?

| |

Entries of the .csv file is an array of cartesian coordinates and a scalar at the location. The array's shape is (num_data_entries, (x, y, z, scalar)). The scalar is one of 0.0, 0.25, 1.0. The actual code to reproduce the rendering is as follows:

import numpy as np
from mayavi import mlab
from mayavi.core.lut_manager import LUTManager
from tvtk.util import ctf  # ColorTransferFunction

def _get_custom_colormap(name, num_colors=256, opacity=1., transparent_input=-1.):
    # Get a mayavi colormap as a numpy array
    lm = LUTManager(number_of_colors=num_colors, lut_mode=name, show_scalar_bar=True)
    rgbs = lm.lut.table.to_array()[:, :3]  # ~ (:, RGBA) -> (:, RGB)
    rgbs = rgbs.astype(float) / (num_colors - 1)  # [0.0, 1.0]**3
    # Setup a color transfer function with the above-generated colormap
    _ctf = ctf.ColorTransferFunction()
    for i, rgb in enumerate(rgbs):
        _ctf.add_rgb_point(i/(num_colors - 1), *rgb)

    # Setup a opacity transfer function
    _otf = ctf.PiecewiseFunction()
    _otf.add_point(transparent_input, 0.)  # transparent
    _otf.add_point(0., opacity)
    _otf.add_point(1., opacity)
    # Extend the input value range to allow transparent plot
    _ctf.range = [transparent_input, 1.]

    return _ctf, _otf

def _init_mayavi_volume_rendering_fig(x, y, z, scalars, _ctf=None, _otf=None):
    # Initial rendering of mass distr
    black = (0, 0, 0)
    white = (1, 1, 1)
    size = (720, 720)

    figure = mlab.figure(bgcolor=white, fgcolor=black, size=size)
    scalar_field = mlab.pipeline.scalar_field(x, y, z, scalars)
    volume = mlab.pipeline.volume(scalar_field, figure=figure)

    ctf_otf_update_necessary = False

    if ctf is not None:
        volume._ctf = _ctf
        volume._volume_property.set_color(_ctf)
        ctf_otf_update_necessary = True
    if _otf is not None:
        volume._otf = _otf
        volume._volume_property.set_scalar_opacity(_otf)
        ctf_otf_update_necessary = True

    if ctf_otf_update_necessary:
        # Apply the custom ctf
        volume.update_ctf = True

    return figure, volume

if __name__ == "__main__":

    # generate gt mass distr ===========================================
    aabb_scale = 0.64
    density_1 = 500
    density_2 = 2000
    # voxel grid spec
    level = 6
    res = 2**level  # resolution
    half_res = res// 2
    # scaled voxel spec
    vx_side = aabb_scale / half_res
    vx_volume = vx_side**3

    # coordinates and mass distribution
    ndc_coords = (np.arange(0, res) - half_res + .5) / half_res
    coords = aabb_scale * ndc_coords

    xyzs = np.array(np.meshgrid(coords, coords, coords, indexing="ij"))  # ~ (3, res, res, res)
    xyzs = xyzs.reshape(3, -1)  # (3, N = res*res*res)

    mass_distr = np.zeros_like(xyzs[0], dtype="float")  # ~ (N, )
    # upper part - - - - - - - - - - - - - - - - - - - - -
    zidx = np.where(.0 <= xyzs[2])
    xidx = np.where((-.25 * aabb_scale <= xyzs[0]) & (xyzs[0] <= .25 * aabb_scale))
    yidx = np.where((-.5 * aabb_scale <= xyzs[1]) & (xyzs[1] <= .5 * aabb_scale))
    xy_idx = np.intersect1d(xidx, yidx)
    idx1 = np.intersect1d(xy_idx, zidx)
    mass_distr[idx1] = density_1 * vx_volume
    # lower part - - - - - - - - - - - - - - - - - - - - -
    zidx = np.where(.0 >= xyzs[2])
    idx2 = np.intersect1d(xy_idx, zidx)
    mass_distr[idx2] = density_2 * vx_volume
    print(f"{mass_distr.sum()=}")

    # write the scalar field into a .csv ===============================
    csv_filename = "gt_mass_distr.csv"
    with open(csv_filename, 'w') as f:
        writer = csv.writer(f)
        for x, y, z, mass in zip(*xyzs, mass_distr):
            writer.writerow([x, y, z, mass])

    # load the .csv and plot w/ matplotlib =============================
    loaded = np.loadtxt(csv_filename, delimiter=',')  # ~ (N, (x, y, z, mass))
    _x, _y, _z, _mass_distr = loaded.T
    x = _x / aabb_scale
    y = _y / aabb_scale
    z = _z / aabb_scale

    # set plot data
    viridis = cm.get_cmap('viridis')
    c = _mass_distr / _mass_distr.max()
    s = np.zeros_like(_mass_distr)
    s[_mass_distr.nonzero()] = 1

    # plot
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.set_aspect('equal')
    size_const = 0.3
    ax.scatter(x, y, z, s=size_const*s, c=c, cmap=viridis)
    plt.show(block=False)

    # volume render w/ mayavi ========================================
    transparent_input = -1.
    meshgrid_shape = (res, res, res)
    scalars = transparent_input * np.ones(meshgrid_shape)
    zero_to_one_md = _mass_distr / _mass_distr.max()
    zero_to_one_md = zero_to_one_md.reshape(meshgrid_shape)
    scalars[zero_to_one_md.nonzero()] = zero_to_one_md[zero_to_one_md.nonzero()]

    print(f"{_mass_distr[_mass_distr.nonzero()].sum()=}")
    print(f"{scalars[scalars.nonzero()].sum()*_mass_distr.max()=}")

    _ctf, _otf = _get_custom_colormap("cool",
                                      num_colors=256,
                                      opacity=.7,
                                      transparent_input=transparent_input)

    figure, volume = _init_mayavi_volume_rendering_fig(x.reshape(meshgrid_shape),
                                                       y.reshape(meshgrid_shape),
                                                       z.reshape(meshgrid_shape),
                                                       scalars,
                                                       _ctf,
                                                       _otf)

    max_epochs = 121
    num_mayavi_camera_turn = 2
    azimuth, _, _, _ = mlab.view()
    azimuth_tick = num_mayavi_camera_turn * 360.0 / max_epochs

    mlab.show()

enthought / mayavi

Scalar filed is not volume rendered with a colormap #1305