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