Getting the samples out of vtkProbeFilter ...

I have a structured grid with a temperature array. I have verified that it has nonzero data stored in this array, checked with:


I am trying to probe a ray through this volume to get the temperatures along this ray. By construction, I know that the endpoints of the ray are on the surface of the volume (and I have visually verified it in ParaView). The code I am using for the probe is:

sampler = vtkLineSource()

probe = vtkProbeFilter()

temps = vtk_to_numpy(probe.GetOutputDataObject(0).GetPointData().GetArray("Temperature"))

However when I inspect the temps array, I get all zeros!
Any ideas what I am doing wrong?

Apparently I have been down this road before, and nobody seemed to know then either:

A VTK wiki that I found that seemed to have some useful information.
I guess another hint for me might be this check (no valid points returned):

(Pdb) p vtk_to_numpy(probe.GetValidPoints())
array([], dtype=int64)

I can’t understand why the geometry would be ill-defined - it’s a vtkStructuredGrid? What’s to mess up? If I can see the file in ParaView, how can it be that wrong?

To add insult to injury - the “demo code” I wrote works :frowning: … I guess that means that there is something wrong with reading the data in from the volume?

EDIT: If I changed volume to be closer to what I’m working with.

#!/usr/bin/env python3
import numpy as np
from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
from vtkmodules.vtkCommonCore import vtkPoints
from vtkmodules.vtkCommonDataModel import vtkStructuredGrid
from vtkmodules.vtkFiltersCore import vtkProbeFilter
from vtkmodules.vtkFiltersSources import vtkLineSource

X, Y, Z = np.meshgrid(
    np.linspace( 0.120, 0.225, 51),
    np.linspace(-0.006, 0.027, 21),
    np.linspace(-0.012, 0.036, 21))
R = np.hypot(X, np.hypot(Y, Z))
W = np.sin(R) / R
W = numpy_to_vtk(W.ravel())

pts = vtkPoints()
for x, y, z in zip(X.ravel(), Y.ravel(), Z.ravel()):
    pts.InsertNextPoint(x, y, z)

sg = vtkStructuredGrid()
sg.SetDimensions(51, 21, 21)

sampler = vtkLineSource()
sampler.SetPoint1(0.120, -0.006, -0.012)
sampler.SetPoint2(0.225,  0.027,  0.036)

probe = vtkProbeFilter()

print(f"valid points: {vtk_to_numpy(probe.GetValidPoints())}")
print(f"probed values: {vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('W'))}")

The issue is the ordering of numpy multidimensional arrays, and in particular the use of meshgrid.

Try switching your meshgrid around like this:

Z, Y, X = np.meshgrid(
    np.linspace(-0.012, 0.036, 21),  # Z space
    np.linspace(-0.006, 0.027, 21),  # Y space
    np.linspace( 0.120, 0.225, 51),  # X space

If your meshgrid is built correctly, then when you print out the points, you should see the X values changing the most rapidly, followed by Y, and then followed by Z:

# Points:
#  0      1      2    ... 22488  22489  22490
[0.12   0.1221 0.1242 ... 0.2208 0.2229 0.225]  # X.ravel()
[-0.006 -0.006 -0.006 ... 0.027  0.027  0.027]  # Y.ravel()
[-0.012 -0.012 -0.012 ... 0.036  0.036  0.036]  # Z.ravel()

On the other hand, if Point 0 and Point 1 have the same X value, then you know that meshgrid has produced a grid that’s not ordered correctly for VTK.

Also be sure to switch around the numpy array dimensions before you use them in VTK:

z_size, y_size, x_size = X.shape
sg.SetDimensions(x_size, y_size, z_size)

Some people prefer to transpose the numpy arrays before using them with VTK, but I prefer to leave all the arrays in their original ordering, and just re-order the indices instead.

1 Like

One minor edit above: use indexing='ij' with meshgrid to avoid swapping the first two dimensions.

1 Like

EDIT: there is something most definetly wrong with my point locations now … I’m still investigating where it went wrong. Viewing my volume in ParaView shows something is terribly wrong.

My suspicion is still that there is something wrong with the structure of your vtkStructuredGrid. Check the points along some of the edges of the grid to make sure they’re correct:

dims = [1,1,1] 
pts = sg.GetPoints()

if pts.GetNumberOfPoints() != dims[0]*dims[1]*dims[2]:
    print('Either points or dimensions are wrong')

for i in range(dims[0]):
for j in range(dims[1]):
for k in range(dims[2]):

For example, if your grid is aligned with the x, y, z axes, then in the xpoints only x will vary, in ypoints only y will vary, and in zpoints only z will vary.

After I corrected the order (X then Y then Z) - there was a transpose operation in there. After that the shape of the array was incorrectly reported to vtkStructuredGrid. In particular, the data was (21, 21, 51) but I told VTK (21, 51, 21) … so that was a disaster.

Thanks for your help.