3

I would like to use Paraview to plot simple 2D meshes with either different colors per cell or different colors per vertices. As far as I can tell, the Paraview documentation does not explain how to Show() a user-defined VTK object.

I read from the Paraview guide how the VTK data model works and from the VTK User's Guide how to generate a vtkImageData object.

From what I could gather, the following code should yield a vtkImageData object of a 10x5 2D mesh spanning [0.;10.]x[0.;5.] with 50 blue elements.

But now I don't know how to actually plot it in Paraview.

from paraview import vtk
import paraview.simple as ps
import numpy as np
from paraview.vtk.util.numpy_support import numpy_to_vtk

def main():
    # create the vtkImageData object
    myMesh = vtk.vtkImageData()
    myMesh.SetOrigin(0.,0.,0.)
    myMesh.SetExtent(0,10,0,5,0,0)
    myMesh.SetSpacing(1.,1.,0.)

    # create the numpy colors for each cell
    blue = np.array([15, 82, 186], dtype=np.ubyte)  # 8 bits each [0, 255]
    npColors = np.tile(blue, (myMesh.GetNumberOfCells(), 1))

    # transform them to a vtkUnsignedCharArray
    # organized as 50 tuples of 3
    vtkColors = numpy_to_vtk(npColors, deep=1, array_type=vtk.VTK_UNSIGNED_CHAR)

    # allocate the sets of 3 scalars to each cell
    myMesh.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 3)  # set 3 scalars per cell
    myMesh.GetCellData().SetScalars(vtkColors)  # returns 0, the index to which
                                                # vtkColors is assigned
    # Something... generate a proxy using servermanager??

    ps.Show(myMesh?)
    ps.Interact()  # or ps.Render()

if __name__ == "__main__":
    main()

From what I could gather is that I have to apply a geometric filter first, such as vtkImageDataGeometryFilter(). But this does not exist in paraview.vtk, only by importing the vtk module directly.

Another option, according to this VTK C++ SO question is using vtkMarchingSquares.

Either way, apparently paraview.simple.Show only accepts a proxy object as input. Which then begs the question of how to create a proxy out of the filtered vtkImageData object? To be honest, I have not quite grasped how the visualization pipeline works, despite reading the docs.

Thus far I've only found ways of visualizing VTK objects using vtk directly via the Kitware examples in GitHub, without using the higher level features of Paraview.

Breno
  • 748
  • 8
  • 18

2 Answers2

2

ProgrammableSource is what you want to use. See this example or this

Utkarsh
  • 1,492
  • 1
  • 11
  • 19
1

I managed to render it using a TrivialProducer object and the method .GetClientSideObject(). This interfaces ParaView to serverside objects.

Sources: the source code and the tip given by Mathieu Westphal from ParaView support.

from paraview import simple as ps
from paraview import vtk
from paraview.vtk.util.numpy_support import numpy_to_vtk
import numpy as np

def main():
    # Create an image (this is a data object)
    myMesh = vtk.vtkImageData()
    myMesh.SetOrigin(0., 0., 0.)
    myMesh.SetSpacing(0.1, 0.1, 0.)
    myMesh.SetExtent(0, 10, 0, 5, 0, 0)

    # coloring
    blue = np.array([15, 82, 186], dtype=np.ubyte)
    # numpy colors
    scalarsnp = np.tile(blue, (myMesh.GetNumberOfCells(), 1))
    scalarsnp[[9, 49]] = np.array([255, 255, 0], dtype=np.ubyte)  # yellow

    # vtk array colors. Organized as 50 tuples of 3
    scalarsvtk = numpy_to_vtk(scalarsnp, deep=1, array_type=vtk.VTK_UNSIGNED_CHAR)
    scalarsvtk.SetName("colorsArray")
    # allocate the scalars to the vtkImageData object
    # myMesh.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 3)  # set 3 scalars per cell
    # myMesh.GetCellData().SetScalars(scalarsvtk)  # do not use this in ParaView!!
    colorArrayID = myMesh.GetCellData().AddArray(scalarsvtk)
    myMesh.GetCellData().SetActiveScalars(scalarsvtk.GetName())

    # TrivialProducer to interface ParaView to serverside objects
    tp_mesh = ps.TrivialProducer(registrationName="tp_mesh")
    myMeshClient = tp_mesh.GetClientSideObject()
    # link the vtkImageData object to the proxy manager
    myMeshClient.SetOutput(myMesh)
    tp_mesh.UpdatePipeline()

    # Filter for showing the ImageData to a plane
    mapTexture2Plane = ps.TextureMaptoPlane(registrationName="TM2P_mesh", Input=tp_mesh)

    renderViewMesh = ps.CreateView("RenderView")
    renderViewMesh.Background = [1, 1, 1]
    renderViewMesh.OrientationAxesVisibility = 0
    display = ps.Show(proxy=mapTexture2Plane, view=renderViewMesh)
    display.SetRepresentationType("Surface")
    display.MapScalars = 0  # necessary so as to not generate a colormap
    ps.Interact()  # or just ps.Render()


if __name__ == "__main__":
    main()
Breno
  • 748
  • 8
  • 18