volume rendering for volumetric data

[wangfudong]

Hi, sorry for my bothering you~ I have nowadays utilized this useful project to render my research datas, but I got some troubles with it. I have output some volumetric data of shape [H, W, D, C] consisting of [R,G,B,A] channels following the Optical Models for Direct Volume Rendering method https://courses.cs.duke.edu/cps296.8/spring03/papers/max95opticalModelsForDirectVolumeRendering.pdf, where channel A is the so-called 'extinction coefficient (= 1- opacity)' or 'density' in some literatures. In this data, each voxel in the grid is assigned with [r,g,b,a] values.

Next I used pyvista to render it as follows:

import pyvista as pv
import plotly.graph_objects as go
import mcubes
import numpy as np

inside = np.load('./data/inside_bool.npy') # shape = [198, 381, 60], bool. Only inside part of the grid data is useful, the other outside part is empty space.
rgba = np.load('./data/rgba_inside.npy') # of shape N x 3, float

dims = inside.shape 

scalars_r = np.zeros(dims);scalars_r[inside] = rgba[:,0]
scalars_g = np.zeros(dims);scalars_g[inside] = rgba[:,1]
scalars_b = np.zeros(dims);scalars_b[inside] = rgba[:,2]
scalars_a = np.zeros(dims);scalars_a[inside] = rgba[:,3]
scalars_a[inside ==0] = None

# scalars_rgb = np.zeros((dims[0],dims[1],dims[2],3))
# scalars_rgb[inside] = rgba[:,0:3]
# grid.point_arrays['scalars'] = scalars_rgb.reshape(-1,3)

## visulization as volume data
grid = pv.UniformGrid(dims)
grid.point_arrays['scalars'] = scalars_a.ravel()
grid.set_active_scalars('scalars')

cpos = [(-600, 220, -360), (99.0, 190.0, 30), (0,-1,0)]
p = pv.Plotter()
opacity = [0, 1, 0]
p.add_volume(grid, opacity=opacity)
p.camera_position = cpos
p.show(screenshot='vis_volume_rgba.png')

## visulization as point cloud/mesh with 'density' rgba[:,3]
Pxyz = np.where(inside==1)
PC = pv.PolyData(np.asarray(Pxyz).transpose())
PC['scalars'] = rgba[:,3] #.max() - rgba[:,3]
# PC.plot(eye_dome_lighting=True)
pc = pv.Plotter()
opacity = [0, 1, 0] 
# opacity = [0, 1]
pc.add_mesh(PC, opacity=opacity, lighting=True, smooth_shading=True)
pc.camera_position = cpos
pc.show(screenshot='vis_PC_rgba_opacity1.png')

##using marching-cubes to extracte the iso-surface of volumetric grid rgba[:,3]
Mcube = np.zeros(dims)
Mcube[inside == 1] = rgba[:,3]
# Mcube = np.pad(Mcube, 10, mode='constant')
vertices, triangles = mcubes.marching_cubes(Mcube, isovalue = 10)
mesh = pv.PolyData(vertices, np.concatenate((np.ones([triangles.shape[0],1])*3,triangles),axis=1).astype(np.int))

pc1 = pv.Plotter()
pc1.add_mesh(mesh, lighting=True, smooth_shading=True)
pc1.camera_position = cpos
pc1.show(screenshot='vis_isosurface_rgba.png')

print('done')

In the code pc.add_mesh(PC, opacity=opacity, lighting=True, smooth_shading=True), I used opacity = [0, 1, 0] or opacity = [0, 1], respectively. The output figures are as follows:

With these rendering results., I have some doubts,

1. Why the volume rendering result looks so strange? Did I make someything wrong with it? And how to render it with the [r,g,b] channels?
2. Comparing the rendering results of pointcloud data with opacity = [0, 1, 0] or opacity = [0, 1], the later looks similar to the isosurface mesh, which is closed to my desired. I read the python code 'add_volume()' and learn how the 'opacity' is calculated numerically and spanned to the the same size with 'n_colors'. However, why can not assign each voxel/point/cell an opacity value? Can you tell me the principle that pyvista used to implement the parameter 'opacity', or articles/books/links are ok.

Thank you very much :)

[adamgranthendry]

This also seems odd to me. I can confirm I get similar results.

@akaszynski Is there something we are doing wrong when volume rendering (upper left corner image)? (FYI, the majority of the volume rendering examples on the docs don't appear for me.)

[adamgranthendry]

@wangfudong Actually, I think I see your problem. You need to switch the line

grid.point_arrays['scalars'] = scalars_a.ravel()

to

grid.point_arrays['scalars'] = scalars_a.flatten(order="F")

FORTRAN ordering is required. This solved the problem for me. Can you confirm it does the same for you?

[adamgranthendry]

@wangfudong In regards to your second question, pyvista.add_volume uses 256 colors on the scale bar by default (controlled by the n_colors argument). Your data has 3 color channels and an alpha channel (R, G, B, and A), but assuming each color channel is 8-bit, this means you can represent (2**8)**3 = 2**24 colors (a big gamut!). The n_colors is made 256 (i.e. 2**8) so the lookup table (LUT) for opacity can be contained in one byte and thus provide faster lookups. The opacity argument is intended to be the same length as n_colors, but you can specify fewer arguments.

In essence, when you specify [0, 1, 0] for opacity, you are setting the opacity for the first 3 colors (out of the total 256 colors) to 0, 1, and 0, respectively. The remaining 253 colors preserve their default values. When you specify [0, 1], you set only the first two colors, and the remaining 254 retain their value.

I believe from your code that you think you are setting the opacity of the color channels, but that is not actually how opacity works in pyvista.add_volume. You can instead create an opacity array of length 256 that has varying degrees of opacity from 0 to 1. For example, opacity = np.linspace(0, 1, 256) would make a linearly increasing opacity from 0 at the low values of your data to 1 at the high values. Instead, pyvista already comes with a bunch of pre-made opacity functions (called opacity transfer functions in VTK and pyvista), so you can specify a string that maps to a pre-made opacity transfer function.

The opacity transfer functions are specified in pyvista.plotting.tools in the function opacity_transfer_function. You have the following to choose from:

 transfer_func = {
        'linear': np.linspace(0, 255, n_colors, dtype=np.uint8),
        'geom': np.geomspace(1e-6, 255, n_colors, dtype=np.uint8),
        'geom_r': np.geomspace(255, 1e-6, n_colors, dtype=np.uint8),
        'sigmoid': sigmoid(np.linspace(-10.,10., n_colors)),
        'sigmoid_3': sigmoid(np.linspace(-3.,3., n_colors)),
        'sigmoid_4': sigmoid(np.linspace(-4.,4., n_colors)),
        'sigmoid_5': sigmoid(np.linspace(-5.,5., n_colors)),
        'sigmoid_6': sigmoid(np.linspace(-6.,6., n_colors)),
        'sigmoid_7': sigmoid(np.linspace(-7.,7., n_colors)),
        'sigmoid_8': sigmoid(np.linspace(-8.,8., n_colors)),
        'sigmoid_9': sigmoid(np.linspace(-9.,9., n_colors)),
        'sigmoid_10': sigmoid(np.linspace(-10.,10., n_colors)),
    }

For example, sigmoid will give you a nice continuous variation in opacity. In addition, if you would prefer to reverse the color bar so that low values have more opacity, you can add _r to any of the strings above and it will reverse the colors for you.

I would recommend using sigmoid as that seems to work nice most of the time. In addition, if you want to "filter out" some values from your data set, I would use clim because I think None values won't work, though I could be wrong. For example, you could try clim=[np.nanmin(scalars_a[scalars_a>0]), np.nanmax(scalars_a)].

@MatthewFlamm @akaszynski This would be another great thing to add to our documentation! The volume rendering section could use some TLC.

[adamgranthendry]

@wangfudong Are you able to share your data set with us so we can analyze the issue on our end as well?

@akaszynski @adeak @MatthewFlamm I've also had scenarios like in vis_volume_rgba happen to me before even with Fortran ordering. I'd like to figure out why and investigate further. Has this every happened to you?