Scalar filed is not volume rendered with a colormap
barikata1984 opened this issue · 1 comments
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()Okay, the black area was assumed to be generated by the shading feature of vtk. After I disabled it by adding the line volume._volume_property.shade = False, the black part disappeared