Demystifying the vtkProbeFilter

vtkProbeFilter

Introduction

One of the most powerful and mystifying classes in the VTK arsenal is the vtkProbeFilter. This article attempts to explain a little bit about how it works and how to use it.

How to probe

From my understanding of the probe, in addition to the normal input, I would also need to set a 'Source' input. The probe would copy the input to the output dataset, setting output point data scalars to the interpolated point data scalars of the input. So the geometry of the output remains the same as the geometry of the input. Sounds good so far-- time to give it a try.

Elevation probe example

The following is an example of a python function to sample elevations in one dataset to another. For example, you might have a triangulated terrain surface and you might want to turn that surface into a grid. First I copy the Z values of the dataset I want to sample to the point data scalars. This can be done with a vtkElevationFilter if you know the Z-range, which I can get by flushing the input dataset and getting the bounds. Next, I flatten both datasets using a transform to ensure all input points are sampled. Now the datasets are ready for the probe. Since the probe only affect scalars data, I have to move that scalar to the Z values of the probe output and this can be done with the vtkWarpScalar class.

#!/usr/bin/env python
### fair warning, I haven't checked if this code runs or not ###
from vtk import *

def elevation_probe(input, source):
    source.Update()
    bounds = source.GetOutput().GetBounds()
    
    ### make a transform to set all Z values to zero ###
    flattener = vtkTransform()
    flattener.SetScale(1.0,1.0,0.0)

    ### flatten the input in case it's not already flat; it's Z values don't matter yet anyway ###
    i_flat = vtkTransformFilter()
    i_flat.SetInput(input)
    i_flat.SetTransform(flattener)

    ### transfer Elevation values to the source's point scalar data ###    
    s_elev = vtkElevationFilter()
    s_elev.SetInput(source)
    s_elev.SetHighPoint(0,0,bounds[5])
    s_elev.SetLowPoint(0,0,bounds[4])
    s_elev.SetScalarRange(bounds[4],bounds[5])

    ### flatten the source data; the Z elevations are already in the scalars data ###
    s_flat = vtkTransformFilter()
    s_flat.SetInput(s_elev.GetOutput())
    s_flat.SetTransform(flattener)

    ### Probe! Since both dataset are flat, every input probe point should get an elevation value
    ### if it falls within the source dataset somewhere. These values are stored in the output scalars
    probe = vtkProbeFilter()
    probe.SetInput(i_flat.GetOutput())
    probe.SetSource(s_flat.GetOutput())

    ### move the elevation scalars in the probe output to the Z values of the corresponding coordinates.
    warp = vtkWarpScalar()
    warp.SetInput(probe.GetOutput())
    warp.Update()
    
    return warp.GetOutput()

vtkProbeFilter Pitfalls

The first time I tried to probe a dataset, I used a triangulation that I read from a shapefile as the source dataset. I noticed in my output, that only one in every six points would actually sample a scalar value, despite the fact that my input clearly overlayed my source dataset. Worse still, the locations of the points missed by the probe seemed to be random. I was truly confounded by this data-processing cataclysm. I spent a fews days digging in VTK source code and think I figured out why this was happening.

Deep within the bowels of the vtkProbeFilter, the method vtkDataSet::FindAndGetCell is called, which returns the cell containing an input point of interest. Interestingly, VTK is exceptionally fast at retrieving cells because instead of iterating through every cell until it finds the one containing the input point, it finds the closest point to the input point using an internal locator and then finds all of the cells that include that closest point in their topologies. The result is that only a few cells in the dataset are iterated. This folks, is truly some brilliant O(log n) computer science.

So what's wrong with me? Well, a shapefile has no topology; points are copied for every cell. vtkDataSet::FindAndGetCell will find the closet point, but it will only use the first one it finds. So when the point locator picks a point, it will pick one of several copies of the same point, each of which connects to a different single cell. A closest point is picked and only one cell is checked against the input point. If it misses, the other cells are never iterated.

So the lesson here is that if you use a probe, make sure your topology is correct. For my example, I could use vtkDelaunay2D to retriangulate the dataset before I probe. I believe that vtkCleanPolyData and vtkConnectivetyFilter can be used together to remove the duplicate points and simplify topologies.

Checking Probe Output

You can check the output for missed points by using vtkProbeFilter::GetValidPoints. For example:

probe->Update(); assert(probe->GetOutput()->GetNumberOfPoints() == probe->GetValidPoints()->GetNumberOfTuples());

--Maik! 08:36, 4 Feb 2005 (EST)