Error with vtkThreshold

[loisvazquezind]

Hello all.

I'm having problems with the colourmap of a vtkThreshold. I have a GUI in which I can change the field to be displayed as well as change the threshold range. Usually, I can interact perfectly with the field to be displayed and the threshold range. However, at some point, multiple errors occur:

1. The threshold does not work:
   

In this case, if I move the lower limit of the threshold to the middle and set it back to the extreme, I get this:

2. The colourmap does not work (the geometry is coloured with a solid colour):
   

The problem is that I can't quite find a pattern in the behaviour of the error that allows me to figure out where the error is coming from. Am I updating the colormap wrong or could it be a problem with the VTU file? In case the problem is with the VTU, what exactly is failing? How can I validate it? In my real code, I am already validating that there are no NaN values, cell types, etc.

My VTU file contains three fields, both in point data and cell data, as I show in the following log:

Available cell arrays (threshold_grid) (3): [
	[id: 0, name: 'Displacement', range: (2.5438041953596195e-12, 0.010237611981167797)],
	[id: 1, name: 'ElasticStrain'], range: (2.5075287430565503e-14, 0.00017659147922402104),
	[id: 2, name: 'VonMises'], range: (6.730424335011106e-05, 12361.403545681473),
]
Available point arrays (threshold_grid) (3): [
	[id: 0, name: 'Displacement', range: (0.0, 0.010264021022249822)],
	[id: 1, name: 'ElasticStrain'], range: (1.481919113327266e-14, 0.00017659147922402104),
	[id: 2, name: 'VonMises'], range: (5.494566415891106e-05, 12361.403545681473),
]

whose content is (I have written '...' instead of copying the entire arrays for simplicity):

<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkLZMADataCompressor">
  <UnstructuredGrid>
    <Piece NumberOfPoints="478248" NumberOfCells="2032442">
      <PointData>
        <DataArray type="Float64" Name="Displacement" format="binary" RangeMin="2.5438041953596195e-12" RangeMax="0.010237611981167797">
          ...
        </DataArray>
        <DataArray type="Float64" Name="ElasticStrain" format="binary" RangeMin="2.5075287430565503e-14" RangeMax="0.00017659147922402104">
          ...
        </DataArray>
        <DataArray type="Float64" Name="VonMises" format="binary" RangeMin="0.00006730424335011106" RangeMax="12361.403545681473">
          ...
        </DataArray>
      </PointData>
      <CellData>
        <DataArray type="Float64" Name="Displacement" format="binary" RangeMin="0" RangeMax="0.010264021022249822">
          ...
        </DataArray>
        <DataArray type="Float64" Name="ElasticStrain" format="binary" RangeMin="1.481919113327266e-14" RangeMax="0.00017659147922402104">
          ...
        </DataArray>
        <DataArray type="Float64" Name="VonMises" format="binary" RangeMin="0.00005494566415891106" RangeMax="12361.403545681473">
          ...
        </DataArray>
      </CellData>
      <Points>
        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="binary" RangeMin="0.019616892102132258" RangeMax="0.07953547271489829">
          ...
          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
            <Value index="0">
              0.019616892102
            </Value>
            <Value index="1">
              0.079535472715
            </Value>
          </InformationKey>
          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
            <Value index="0">
              0.019616892102
            </Value>
            <Value index="1">
              0.079535472715
            </Value>
          </InformationKey>
        </DataArray>
      </Points>
      <Cells>
        <DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="478247">
          ...
        </DataArray>
        <DataArray type="Int64" Name="offsets" format="binary" RangeMin="4" RangeMax="8129768">
          ...
        </DataArray>
        <DataArray type="UInt8" Name="types" format="binary" RangeMin="10" RangeMax="10">
          ...
        </DataArray>
      </Cells>
    </Piece>
  </UnstructuredGrid>
</VTKFile>

I have tested with a VTU that I have generated from the FiniteElementAnalysis example in which I have duplicated the field to have a VTU with two fields, both in point data and in cell data. With this file it doesn't seem to fail no matter how much I change the field to be displayed and the threshold range (I can provide it to you if necessary).

My (simplified) code is:

import os
import logging
# Trame
from trame.app import get_server
from trame.ui.vuetify import SinglePageWithDrawerLayout
from trame.widgets import vtk, vuetify, trame, html
# VTK
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonDataModel import vtkDataObject, vtkUnstructuredGrid
from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader
from vtkmodules.vtkRenderingAnnotation import vtkScalarBarActor
from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
from vtkmodules.vtkFiltersCore import vtkThreshold
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkDataSetMapper,
    vtkRenderer,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
)
from vtkmodules.vtkCommonCore import (
    vtkLookupTable,
)
# Required for interactor initialization
from vtkmodules.vtkInteractionStyle import vtkInteractorStyleSwitch  # noqa
# Required for rendering initialization, not necessary for
# local rendering, but doesn't hurt to include it
import vtkmodules.vtkRenderingOpenGL2  # noqa

# -----------------------------------------------------------------------------
# Trame initialization
# -----------------------------------------------------------------------------

server = get_server()
server.client_type = "vue2"
state, ctrl, ui = server.state, server.controller, server.ui
state.setdefault("active_ui", "threshold")  # set the default active UI

state.update({
    # Common parameters
    "dataset_arrays": [
        {
            "text": "Displacement",
            "value": 0,
        },
        {
            "text": "ElasticStrain",
            "value": 1,
        },
        {
            "text": "VonMises",
            "value": 2,
        },
    ],
    "default_color_array_idx": None,
    "color_array_idx": None,
    # Thresholding
    "threshold_enabled": False,
    "threshold_range": [0, 1],
    "threshold_full_min": 0,
    "threshold_full_max": 1,
    "threshold_step": 0.01,
    "threshold_step_factor": 0.001,
})

# -----------------------------------------------------------------------------
# VTK code
# -----------------------------------------------------------------------------
# Set colors (in RGBA)
colors = vtkNamedColors()
colors.SetColor("BkgColor", [82, 87, 110, 255])  # paraview background color

renderer = vtkRenderer()
renderer.GradientBackgroundOn()
renderer.SetBackground(colors.GetColor3d("BkgColor"))  # Other: SlateGray
render_window = vtkRenderWindow()
render_window.AddRenderer(renderer)
render_window_interactor = vtkRenderWindowInteractor()
interactor_style = vtkInteractorStyleTrackballCamera()
render_window_interactor.SetInteractorStyle(interactor_style)
render_window_interactor.SetRenderWindow(render_window)

# Declare variables
threshold_grid = vtkUnstructuredGrid()
threshold_filter = vtkThreshold()
threshold_filter_mapper = vtkDataSetMapper()
threshold_filter_actor = vtkActor()
threshold_scalar_bar_actor = vtkScalarBarActor()
threshold_lut = vtkLookupTable()

input_file_path = "/path/to/file.vtu"


def read_source(path: str):
    if not os.path.isfile(path):
        raise Exception(f"File not found at '{path}'")
    reader = vtkXMLUnstructuredGridReader()
    reader.SetFileName(path)
    reader.Update()
    # Extract array information from source
    data_arrays = []
    fields = [
        (reader.GetOutput().GetPointData(), vtkDataObject.FIELD_ASSOCIATION_POINTS),
        (reader.GetOutput().GetCellData(), vtkDataObject.FIELD_ASSOCIATION_CELLS),
    ]
    for field in fields:
        field_arrays, association = field
        for i in range(field_arrays.GetNumberOfArrays()):
            array = field_arrays.GetArray(i)
            _index = i
            data_arrays.append({
                "text": str(array.GetName()),
                "value": _index,
                "range": array.GetRange(),
                "type": association,
            })
    # Sort ascending by field 'text'
    data_arrays = sorted(data_arrays, key=lambda d: d["value"])
    return reader, data_arrays


# Read file and get default array (here, vtk_source is the vtkXMLUnstructuredGridReader())
vtk_source, dataset_arrays = read_source(path=input_file_path)
# Get default array (Displacement)
default_array_idx = 0
default_array = dataset_arrays[default_array_idx]

# Update state
state.dataset_arrays = dataset_arrays
state.default_color_array_idx = default_array_idx
state.color_array_idx = default_array_idx

# Initialize
threshold_filter.SetInputData(threshold_grid)
threshold_filter_mapper.SetInputConnection(threshold_filter.GetOutputPort())
threshold_filter_mapper.SetScalarRange(vtk_source.GetOutput().GetScalarRange())
# Set the lookup table as rainbow
threshold_lut = threshold_filter_mapper.GetLookupTable()
threshold_lut.SetHueRange(0.666, 0.0)
threshold_lut.SetSaturationRange(1.0, 1.0)
threshold_lut.SetValueRange(1.0, 1.0)
threshold_filter_mapper.SetLookupTable(threshold_lut)
threshold_filter_actor.SetMapper(threshold_filter_mapper)
# Scalar bar
threshold_scalar_bar_actor.SetLookupTable(threshold_filter_mapper.GetLookupTable())
threshold_scalar_bar_actor.SetTitle(default_array.get("text"))

threshold_grid.ShallowCopy(vtk_source.GetOutput())
cell_array = vtk_source.GetOutput().GetCellData().GetArray(default_array.get("text"))
threshold_filter.SetInputArrayToProcess(0, 0, 0, vtkDataObject.FIELD_ASSOCIATION_CELLS,
                                        default_array.get("text"))
threshold_filter.Update()
state.threshold_full_min = cell_array.GetRange()[0]
state.threshold_full_max = cell_array.GetRange()[1]
state.threshold_range = cell_array.GetRange()
state.threshold_step = state.threshold_step_factor * (state.threshold_full_max - state.threshold_full_min)

# Add actors to renderer
renderer.AddActor(threshold_filter_actor)
renderer.AddActor2D(threshold_scalar_bar_actor)

renderer.ResetCamera()


def update_colormap_threshold(arr_idx):
    # Print available point/cell arrays
    logging.debug(
        f"Available cell arrays (threshold_grid) ({threshold_grid.GetCellData().GetNumberOfArrays()}): [")
    for i in range(threshold_grid.GetCellData().GetNumberOfArrays()):
        logging.debug(f"\t[id: {i}, name: '{threshold_grid.GetCellData().GetArray(i).GetName()}'],")
    logging.debug(f"]")
    logging.debug(
        f"Available point arrays (threshold_grid) ({threshold_grid.GetPointData().GetNumberOfArrays()}): [")
    for i in range(threshold_grid.GetPointData().GetNumberOfArrays()):
        logging.debug(f"\t[id: {i}, name: '{threshold_grid.GetPointData().GetArray(i).GetName()}'],")
    logging.debug(f"]")

    # Get array from 'dataset_arrays' and from 'threshold_grid'
    logging.debug(f"Getting cell data array for field with id={arr_idx}")
    celldata_array = list(
        filter(lambda x: x["value"] == arr_idx and x["type"] == vtkDataObject.FIELD_ASSOCIATION_CELLS,
               state.dataset_arrays))[0]
    celldata_arr = threshold_grid.GetCellData().GetArray(arr_idx)
    if celldata_array.get("text") != celldata_arr.GetName():
        raise Exception(f"bad name. Current: {celldata_arr.GetName()}. "
                        f"Expected: {celldata_array.get('text')}")
    if celldata_array.get("range") != celldata_arr.GetRange():
        raise Exception(f"bad range. Current: {celldata_arr.GetRange()}. "
                        f"Expected: {celldata_array.get('range')}")
    logging.debug(f"Found field: {celldata_array}")

    # Associate the array to the threshold filter
    threshold_filter.SetInputArrayToProcess(0, 0, 0,
                                            vtkDataObject.FIELD_ASSOCIATION_CELLS,
                                            celldata_arr.GetName())
    # Update state
    state.threshold_full_min = celldata_arr.GetRange()[0]
    state.threshold_full_max = celldata_arr.GetRange()[1]
    state.threshold_range = celldata_arr.GetRange()
    state.threshold_step = state.threshold_step_factor * (state.threshold_full_max - state.threshold_full_min)
    # Update LUT
    threshold_lut.SetRange(celldata_arr.GetRange())
    threshold_lut.Build()
    threshold_filter_actor.GetMapper().SetLookupTable(threshold_lut)
    # Updating colormap
    update_colormap(actor=threshold_filter_actor, arr=celldata_array)
    # Update scalar bar
    threshold_scalar_bar_actor.SetTitle(f"{celldata_arr.GetName()}")
    threshold_scalar_bar_actor.SetLookupTable(threshold_filter_actor.GetMapper().GetLookupTable())

    # Update View
    ctrl.view_update()


def update_colormap(actor, arr):
    _min, _max = arr.get("range")
    mapper = actor.GetMapper()
    if not mapper:
        return
    mapper.SelectColorArray(arr.get("text"))
    mapper.SetScalarRange([_min, _max])
    if arr.get("type") == vtkDataObject.FIELD_ASSOCIATION_POINTS:
        mapper.SetScalarModeToUsePointFieldData()
    else:
        mapper.SetScalarModeToUseCellFieldData()
    actor.SetMapper(mapper)


@state.change("color_array_idx")
def update_color_by_name(color_array_idx, **kwargs):
    update_colormap_threshold(color_array_idx)
    state.color_array_idx = color_array_idx
    ctrl.view_update()


@state.change("threshold_range")
def update_threshold_range(threshold_range, **kwargs):
    if threshold_range:
        # Color the threshold geometry
        threshold_filter_mapper.SetScalarRange(threshold_range)
        threshold_filter.SetLowerThreshold(float(threshold_range[0]))
        threshold_filter.SetUpperThreshold(float(threshold_range[1]))
        threshold_filter.Update()
        # Update the scalar bar according to new threshold range
        threshold_lut.SetRange(*threshold_range)
        threshold_lut.Build()
        # Update view
        ctrl.view_update()


# -----------------------------------------------------------------------------
# GUI
# -----------------------------------------------------------------------------
def ui_card(title, ui_name):
    with vuetify.VCard(v_show=f"active_ui == '{ui_name}'"):
        vuetify.VCardTitle(
            title,
            classes="grey lighten-1 py-1 grey--text text--darken-3",
            style="user-select: none; cursor: pointer",
            hide_details=True,
            dense=True,
        )
        content = vuetify.VCardText(classes="py-2")
    return content


def threshold_card():
    with ui_card(title="Threshold", ui_name="threshold"):
        # Field and colormap selector
        with vuetify.VRow(classes="pt-2", dense=True):
            with vuetify.VCol(cols="12"):
                vuetify.VSelect(
                    label="Field",
                    v_model=("color_array_idx", state.default_color_array_idx),
                    items=("array_list", state.dataset_arrays),
                    hide_details=True,
                    dense=True,
                    outlined=True,
                    classes="pt-1",
                    disabled=(f"trame__busy",),
                )
        with html.Div(classes="text-center mt-6"):
            vuetify.VDivider(classes="mb-2")
            vuetify.VSpacer()

        with vuetify.VRow(classes="pt-2", dense=True):
            # Slider for threshold
            vuetify.VRangeSlider(
                label="Threshold",
                v_model=("threshold_range", state.threshold_range),
                min=("threshold_full_min", state.threshold_full_min),
                max=("threshold_full_max", state.threshold_full_max),
                step=("threshold_step", state.threshold_step),
                dense=True,
                hide_details=True,
                style="max-width: 400px",
                disabled=(f"trame__busy",),
            )


def actives_change(ids):
    # Do nothing
    pass


def visibility_change(event):
    # Do nothing
    pass


def left_menu():
    trame.GitTree(
        sources=("pipeline", [{"id": "1", "parent": "0", "visible": 1, "name": "Threshold"}, ],),
        actives_change=(actives_change, "[$event]"),
        visibility_change=(visibility_change, "[$event]"),
    )
    vuetify.VDivider(classes="mb-2")
    threshold_card()


with SinglePageWithDrawerLayout(server) as layout:
    layout.icon.click = ctrl.view_reset_camera

    with layout.toolbar:
        vuetify.VSpacer()
        vuetify.VDivider(vertical=True, classes="mx-2")
        vuetify.VProgressLinear(
            indeterminate=True,
            absolute=True,
            bottom=True,
            active=("trame__busy",),
        )

    # Drawer components
    with layout.drawer as drawer:
        drawer.width = 325
        left_menu()

    with layout.content:
        with vuetify.VContainer(fluid=True, classes="pa-0 fill-height"):
            view = vtk.VtkLocalView(
                render_window, namespace="view", mode="local", interactive_ratio=1,
            )
            ctrl.view_update = view.update
            ctrl.view_reset_camera = view.reset_camera
            ctrl.on_server_ready.add(view.update)

# -----------------------------------------------------------------------------
# Main
# -----------------------------------------------------------------------------

if __name__ == "__main__":
    server.start()

Thank you very much in advance.

[jourdain]

Did you try with VtkRemoteView to make sure your VTK code is valid? At that point the issue could be related to the VtkLocalView.

[loisvazquezind]

Indeed, with VtkRemoteView it works perfectly. Is it planned to make the threshold work in local mode in the near future?

[jourdain]

The local view should work but currently has some limitation. That is why it is good to validate your code with the remote one.
It will be good to investigate why the field get removed when local rendered. But in general threshold with colors work in local rendering since that app is working.