Getting the samples out of vtkProbeFilter ...

Support
1

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

vtk_to_numpy(volume.GetPointData().GetArray("Temperature")).

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()
sampler.SetResolution(30)
sampler.SetPoint1(start)
sampler.SetPoint2(end)

probe = vtkProbeFilter()
probe.SetInputConnection(sampler.GetOutputPort())
probe.SetSourceData(volume)
probe.Update()

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?

2

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

3

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?

4

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())
W.SetName("W")

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

sg = vtkStructuredGrid()
sg.SetPoints(pts)
sg.SetDimensions(51, 21, 21)
sg.GetPointData().AddArray(W)

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

probe = vtkProbeFilter()
probe.SetInputConnection(sampler.GetOutputPort())
probe.SetSourceData(sg)
probe.Update()

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

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
    indexing='ij')

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
6

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

1 Like
7

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.

8

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] 
sg.GetDimensions(dims)
pts = sg.GetPoints()

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

print('xpoints:')
for i in range(dims[0]):
    print(pts.GetPoint(i))
print('ypoints:')
for j in range(dims[1]):
    print(pts.GetPoint(j*dims[0]))
print('zpoints:')
for k in range(dims[2]):
    print(pts.GetPoint(k*dims[0]*dims[1]))

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.

9

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.