Search code examples
pythonvtkparaview

How to render a simple 2D vtkImageData object in Paraview?


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.


Solution

  • 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()