Probing/sampling/interpolating VTK data using python TVTK or MayaVi

Question

I would like to visualise a VTK data file (OpenFOAM output) using python. The plot I would like to make is a 1-d line plot of a quantity between two endpoints. To do so, the unstructured data should be interpolated on the points which lie between the two endpoints.

I've used the package Mayavi to visualise the VTK data. At the mayavi webpage there is a description of probing a single value from a scalarfield. This function does not work on a VTK file.

Also I've found a delaunay3d method (mlab.pipeline.delaunay3d) at the mayavi webpage. I did not get this one to work either.

Could anyone advise me how to interpolate my data?

Astrict

Your question is extremely unclear. Can you explain better what you are trying to do? – Intrigue

Excuses for not stating the question clear enough. I've got a VTK data file (OpenFOAM output) and I would like to make a 1-d line plot of a quantity between two points. – Astrict

The phrase "making a 1-d line plot of a quantity between two points means" is almost a literal restatement of what you said in the question and I still don't know what it means. What does the data in your vtk file represent and what about it are you trying to measure? – Intrigue

I data in the VTK-file represents the pressure in a fluid domain. I want to visualise variation in pressure along a line in the 3D space. Then I would get a plot with on the x-axis the distance along the line and on the y-axis the pressure at that point on the line. In this picture you can see the 3D volume with the line. And here an interpolation over a line. For this example I am using paraView, which does not allow for scripting. Also it does not capture well discontinuities. – Astrict

Answer 1

Eventually, I found the answer to my own question. Just in case anyone visits this topic ans has the same problems, I will post my solution. I've used the vtkProbeFilter to interpolate my VTK-data. After the interpolation I've transformed the VTK-line to a numpy array for plotting convenience.

#!/usr/bin/env python
import numpy as np
from vtk.util import numpy_support as VN
from matplotlib import pyplot as plt
import vtk

def readVTK(filename):
    #read the vtk file with an unstructured grid
    reader = vtk.vtkUnstructuredGridReader()
    reader.SetFileName(filename)
    reader.ReadAllVectorsOn()
    reader.ReadAllScalarsOn()
    reader.Update()
    return reader

def createLine(p1,p2,numPoints):
    # Create the line along which you want to sample
    line = vtk.vtkLineSource()
    line.SetResolution(numPoints)
    line.SetPoint1(p1)
    line.SetPoint2(p2)
    line.Update()
    return line

def probeOverLine(line,reader):
    #Interpolate the data from the VTK-file on the created line.
    data = reader.GetOutput()
    # vtkProbeFilter, the probe line is the input, and the underlying dataset is the source.
    probe = vtk.vtkProbeFilter()
    probe.SetInputConnection(line.GetOutputPort())
    probe.SetSource(data)
    probe.Update()
    #get the data from the VTK-object (probe) to an numpy array
    q=VN.vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('U'))
    numPoints = probe.GetOutput().GetNumberOfPoints() # get the number of points on the line
    #intialise the points on the line    
    x = np.zeros(numPoints)
    y = np.zeros(numPoints)
    z = np.zeros(numPoints)
    points = np.zeros((numPoints , 3))
    #get the coordinates of the points on the line
    for i in range(numPoints):
        x[i],y[i],z[i] = probe.GetOutput().GetPoint(i)
        points[i,0]=x[i]
        points[i,1]=y[i]
        points[i,2]=z[i]
    return points,q

def setZeroToNaN(array):
    # In case zero-values in the data, these are set to NaN.
    array[array==0]=np.nan
    return array

#Define the filename of VTK file
filename='a-VTK-file.vtk'

#Set the points between which the line is constructed.
p1=[0.0,-0.1,0.0]
p2=[0.0,-0.1,1.0]

#Define the numer of interpolation points
numPoints=100

reader = readVTK(filename) # read the VTKfile
line=createLine(p1,p2,numPoints) # Create the line
points,U =  probeOverLine(line,reader) # interpolate the data over the line

U = setZeroToNaN(U) # Set the zero's to NaN's
plt.plot(points[:,2],U[:,0]) #plot the data
plt.show()

Astrict