How handle NaN to display a volume

[ivan-marroquin]

Hi,

Thanks for such a great package!

I am using Python 3.9.10 with vedo 2024.5.2 and vtk 9.0.3. I have a dataset where NaNs surround the numeric values (I attached the data for your review).

To display it, I followed the Slicer3DPlotter example. The data is loaded into a volume and then I generate a legosurface. The rendered window displays the volumetric legosurface and the 3D slices extracted from the volume. Although I replaced the NaNs with 0s, the color distribution is skewed by those 0s.

Any suggestions?

Thanks a lot for your help

Ivan

####
import numpy as np
from vedo import *
from vedo.addons import S
import os


def slider_function_x(widget, event):
    global current_i
    
    i= int(xslider.value)

    if (i == current_i):
        return
    
    current_i= i
    xslice= vol.xslice(i).lighting("", la, ld, 0)
    xslice.cmap(cmap_slicer, vmin= rmin, vmax= rmax)
    xslice.name= "XSlice"
    plt.remove("XSlice")
    
    if (0 < i < dims[0]):
        plt.add(xslice)
    
    plt.render()

def slider_function_y(widget, event):
    global current_j
    
    j= int(yslider.value)
            
    if (j == current_j):
        return
    
    current_j= j
    yslice= vol.yslice(j).lighting("", la, ld, 0)
    yslice.cmap(cmap_slicer, vmin= rmin, vmax= rmax)
    yslice.name= "YSlice"
    plt.remove("YSlice")
    
    if (0 < j < dims[1]):
        plt.add(yslice)
    
    plt.render()

def slider_function_z(widget, event):
    global current_k
    
    k= int(zslider.value)
    
    if (k == current_k):
        return
    
    current_k= k
    zslice= vol.zslice(k).lighting("", la, ld, 0)
    zslice.cmap(cmap_slicer, vmin= rmin, vmax= rmax)
    zslice.name= "ZSlice"
    plt.remove("ZSlice")
    
    if (0 < k < dims[2]):
        plt.add(zslice)
    
    plt.render()
    
    
if __name__ == "__main__":
    # get path of working directory
    working_dir= os.path.dirname(__file__)
    
    # set working directory
    os.chdir(working_dir)

    # load data set and convert into a DataFrame
    dataset= np.load('dataset.npy')
    
    # invert the 'z-axis' in a way the shallow depth values
    # are display on top in the rendered window
    dataset= np.flip(dataset, axis= 2)
    
    min_value= np.nanmin(dataset); max_value= np.nanmax(dataset)
    
    # replace NaNs with a value to mask them in
    # the rendered window
    nan_ind= np.isnan(dataset)
    dataset[nan_ind]= 0
    
    vol= Volume(dataset, spacing= [15, 15, 2]).mode(0)
    vol.cmap('turbo', alpha= 1.0, vmin= min_value, vmax= max_value)
    
    lego= vol.legosurface(vmin= min_value - 0.1, vmax= max_value + 0.1)
    lego.cmap('turbo', vmin= min_value, vmax= max_value)
    lego.linewidth(0).alpha(1.0)
    
    # add slider options to rendered window.
    # based on code at
    # https://github.com/marcomusy/vedo/blob/master/vedo/applications.py    
    plt= Plotter()
    
    plt += [lego]
        
    la, ld= 0.7, 0.3  # ambient, diffuse
    dims= vol.dimensions()
    rmin= np.nanquantile(S_wave_field, q= 0.20); rmax= np.nanquantile(S_wave_field, q= 0.90)
    
    cmap_slicer= 'turbo'
    current_i= None
    current_j= None
    current_k= int(dims[2] / 2)
   
    xslice= None
    yslice= None
    zslice= None
    
    zslice= vol.zslice(current_k).lighting("", la, ld, 0)
    zslice.name= "ZSlice"
    zslice.cmap(cmap_slicer, vmin= rmin, vmax= rmax)
    plt.add(zslice)
    
    cx, cy, cz, ch= "dr", "dg", "db", (0.3, 0.3, 0.3)
    
    xslider= plt.add_slider(slider_function_x, 0, dims[0], title= "", 
                            title_size=0.5, pos= [(0.8, 0.12), (0.95, 0.12)],
                            show_value= False, c= cx)
    
    yslider= plt.add_slider(slider_function_y, 0, dims[1], title= "",
                            title_size= 0.5, pos= [(0.8, 0.08), (0.95, 0.08)],
                            show_value= False, c= cy)
    
    zslider= plt.add_slider(slider_function_z, 0, dims[2], title= "",
                            title_size= 0.6, value= int(dims[2] / 2),
                            pos= [(0.8, 0.04), (0.95, 0.04)], show_value= False,
                            c= cz)
    
    plt.show(viewup= 'z', axes= 1).interactive().close()
###

dataset.zip

[ivan-marroquin]

sorry for the confusion, the line of code:
rmin= np.nanquantile(S_wave_field, q= 0.20); rmax= np.nanquantile(S_wave_field, q= 0.90)

should be replaced by:
rmin= np.nanquantile(dataset, q= 0.20); rmax= np.nanquantile(dataset, q= 0.90)

[marcomusy]

Hi @ivan-marroquin
thanks!

What do you mean by skewed? You mean to erase the blue part?

from vedo import *


def slider_function_x(widget, event):
    global current_i

    i = int(xslider.value)
    if i == current_i:
        return

    current_i = i
    xslice = vol.xslice(i).lighting("", la, ld, 0)
    xslice.cut_with_scalar(0.0, "input_scalars", invert=True)
    xslice.cmap(cmap_slicer, vmin=rmin, vmax=rmax)
    xslice.name = "XSlice"
    plt.remove("XSlice")

    if 0 < i < dims[0]:
        plt.add(xslice)


def slider_function_y(widget, event):
    global current_j

    j = int(yslider.value)

    if j == current_j:
        return

    current_j = j
    yslice = vol.yslice(j).lighting("", la, ld, 0)
    yslice.cut_with_scalar(0.0, "input_scalars", invert=True)
    yslice.cmap(cmap_slicer, vmin=rmin, vmax=rmax)
    yslice.name = "YSlice"

    if 0 < j < dims[1]:
        plt.remove("YSlice").add(yslice)


def slider_function_z(widget, event):
    global current_k

    k = int(zslider.value)
    if k == current_k:
        return

    current_k = k
    zslice = vol.zslice(k).lighting("", la, ld, 0)
    zslice.cut_with_scalar(rmin+.1, "input_scalars", invert=True)
    zslice.cmap(cmap_slicer, vmin=rmin, vmax=rmax)
    zslice.name = "ZSlice"

    if 0 < k < dims[2]:
        plt.remove("ZSlice").add(zslice)


if __name__ == "__main__":

    # load data set and convert into a DataFrame
    dataset = np.load("dataset.npy")

    # invert the 'z-axis' in a way the shallow depth values
    # are display on top in the rendered window
    dataset = np.flip(dataset, axis=2)

    min_value = np.nanmin(dataset)
    max_value = np.nanmax(dataset)

    # replace NaNs with a value to mask them in
    # the rendered window
    nan_ind = np.isnan(dataset)
    dataset[nan_ind] = 0

    settings.default_font = "Roboto"
    cmap_slicer = "turbo"

    vol = Volume(dataset, spacing=[15, 15, 2]).mode(0)
    vol.cmap(cmap_slicer, alpha=1.0, vmin=min_value, vmax=max_value)

    lego = vol.legosurface(vmin=min_value - 0.1, vmax=max_value + 0.1)
    lego.cmap(cmap_slicer, vmin=min_value, vmax=max_value)
    lego.add_scalarbar3d(c="black", title="scalar value")
    lego.scalarbar= lego.scalarbar.clone2d("center-right", size=0.2)
    lego.linewidth(0).alpha(0.2).lighting("off")

    # add slider options to rendered window.
    plt = Plotter(size=(1400,1200))

    la, ld = 0.7, 0.3  # ambient, diffuse
    dims = vol.dimensions()
    rmin = np.nanquantile(dataset, q=0.20)
    rmax = np.nanquantile(dataset, q=0.90)

    current_i = None
    current_j = None
    current_k = 0

    xslice = None
    yslice = None
    zslice = None

    zslice = vol.zslice(current_k).lighting("", la, ld, 0)
    zslice.name = "ZSlice"
    zslice.cmap(cmap_slicer, vmin=rmin, vmax=rmax)
    plt.add(zslice)

    cx, cy, cz, ch = "dr", "dg", "db", (0.3, 0.3, 0.3)

    xslider = plt.add_slider(
        slider_function_x,
        0,
        dims[0],
        title="",
        title_size=0.5,
        pos=[(0.8, 0.12), (0.95, 0.12)],
        show_value=False,
        c=cx,
    )

    yslider = plt.add_slider(
        slider_function_y,
        0,
        dims[1],
        title="",
        title_size=0.5,
        pos=[(0.8, 0.08), (0.95, 0.08)],
        show_value=False,
        c=cy,
    )

    zslider = plt.add_slider(
        slider_function_z,
        0,
        dims[2],
        title="",
        title_size=0.6,
        pos=[(0.8, 0.04), (0.95, 0.04)],
        show_value=False,
        c=cz,
    )

    plt.show(lego, viewup="z", axes=1)
    plt.interactive().close()

[ivan-marroquin]

The presence of 0s was skewing the colors, producing a big splash of dark blue.

Great tip to combine volume slicing with cutting a mesh. Thanks for your prompt reply!

[ivan-marroquin]

A quick question. Is it possible to display isolines on a slice plan?

I tried the following:

def slider_function_z(widget, event):
    global current_k
    
    k= int(zslider.value)
    
    if (k == current_k):
        return
    
    current_k= k
    zslice= vol.zslice(k).lighting("", la, ld, 0)
    zslice.cut_with_scalar(rmin + 0.1, "input_scalars", invert= True)
    zslice.isolines(n= 10).color("black")
    zslice.cmap(cmap_slicer, vmin= rmin, vmax= rmax)
    zslice.name= "ZSlice"
    
    if (0 < k < dims[2]):
        plt.remove("ZSlice").add(zslice)

Thanks!

[marcomusy]

You just forgot to add the isolines to the rendering:

def slider_function_z(widget, event):
    global current_k

    k = int(zslider.value)
    if k == current_k:
        return
    current_k = k

    zslice= vol.zslice(k)
    zslice.cut_with_scalar(rmin + 0.1, "input_scalars", invert= True)
    zslice.triangulate()
    zslice.name= "ZSlice"
    zslice.cmap(cmap_slicer, vmin= rmin, vmax= rmax).lighting("off")
    isos = zslice.isolines(vmin=3, vmax=4, n=12).c("black")
    isos.name = "ZSlice" # can give it the same name so it's simpler to remove

    if (0 < k < dims[2]):
        plt.remove("ZSlice").add(zslice, isos)

You can also remove the unnecessary call to .interactive() in the last line.
Using .lighting("off") will likely give you a better visualization.

[marcomusy]

i've been playing a bit with your script... maybe this is an improvement (?)

from vedo import *


def render_slice(vslice, name):
    vslice.cut_with_scalar(rmin, "input_scalars", invert=True)
    vslice.triangulate()
    vslice.cmap(cmap_slicer, vmin=rmin, vmax=rmax).lighting("off")
    isos = vslice.isolines(vmin=rmin, vmax=rmax, n=12).c("black")
    vslice.name = name
    isos.name = name
    plt.remove(name).add(vslice, isos)


def slider_function_x(widget, event):
    global previous_i
    i = int(widget.value)
    if i == previous_i:
        return
    previous_i = i
    render_slice(vol.xslice(i), "XSlice")

def slider_function_y(widget, event):
    global previous_j
    j = int(widget.value)
    if j == previous_j:
        return
    previous_j = j
    render_slice(vol.yslice(j), "YSlice")


def slider_function_z(widget, event):
    global previous_k
    k = int(widget.value)
    if k == previous_k:
        return
    previous_k = k
    render_slice(vol.zslice(k), "ZSlice")


if __name__ == "__main__":

    settings.default_font = "Roboto"
    cmap_slicer = "RdBu"

    # load data set and convert into a DataFrame
    dataset = np.load("dataset.npy")
    # dataset = np.flip(dataset, axis=2)

    min_value = np.nanmin(dataset)
    max_value = np.nanmax(dataset)

    # replace NaNs with a value to mask them in the rendered window
    nan_ind = np.isnan(dataset)
    dataset[nan_ind] = 0

    vol = Volume(dataset, spacing=[15, 15, 2])
    dims = vol.dimensions()
    rmin = np.nanquantile(dataset, q=0.55)
    rmax = np.nanquantile(dataset, q=0.95)
    print('quantile rmin, rmax', rmin, rmax)

    iso = vol.isosurface(rmin).smooth()
    iso.cmap(cmap_slicer, vmin=min_value, vmax=max_value)
    iso.add_scalarbar3d(c="black", title="scalar value")
    iso.scalarbar = iso.scalarbar.clone2d("center-right", size=0.2)
    iso.c("k5").alpha(0.1).lighting("off").wireframe().pickable(False)

    plt = Plotter(size=(1400, 1200))

    previous_i = 0
    previous_j = 0
    previous_k = 0

    xslider = plt.add_slider(
        slider_function_x,
        0, dims[0],
        pos=[(0.8, 0.12), (0.95, 0.12)],
        show_value=False,
        c="dr",
    )

    yslider = plt.add_slider(
        slider_function_y,
        0, dims[1],
        pos=[(0.8, 0.08), (0.95, 0.08)],
        show_value=False,
        c="dg",
    )

    zslider = plt.add_slider(
        slider_function_z,
        0, dims[2],
        pos=[(0.8, 0.04), (0.95, 0.04)],
        show_value=False,
        c="db",
    )

    plt.show(iso, viewup="z", axes=1).close()

[ivan-marroquin]

Wow! It is a great improvement to the volume visualization. Thanks for your input