problem with interpolating vtkImageData on vtkPolyData with vtkPointInterpolator

Support
1

Hi all,
I’m trying to interpolate a .vti volume on a .vtp mesh using vtkPointInterpolator but I get “Segmentation fault (core dumped)” when I’m updating the vtkPointInterpolator filter. Here is my code:

import vtk
# Read vtkPolyData
reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName('surfacemesh.vtp')
reader.Update()
surface = reader.GetOutput()

# Read vtkImage
reader = vtk.vtkImageReader()
reader.SetFileName('volume.vti')
reader.Update()
image3D = reader.GetOutput()

# Gaussian kernel
gaussian_kernel = vtk.vtkGaussianKernel()
gaussian_kernel.SetSharpness(2)
gaussian_kernel.SetRadius(12)

interpolator = vtk.vtkPointInterpolator()
interpolator.SetInputData(surface)
interpolator.SetSourceData(image3D)
interpolator.SetKernel(gaussian_kernel)
interpolator.Update()

Any ideas to solve this issue?

Thanks!

2

Try the vtkXMLImageDataReader (the vtkImageReader is intended for raw files, not .vti files).

3

Hi David,
thanks for your advice but I got the same error. Here you can find the data: surfacemesh.vtp and volume.vti
surfacemesh.vtp (20.3 KB)
volume.vti (239.7 KB)

4

Hmm, the segfault occurs in vtkStaticPointLocator. Maybe a bug?

Try using a different locator to see if it works:

interpolator.SetLocator(vtk.vtkMergePoints())
1 Like
5

Thanks for your help @dgobbi , now it works ! Please find below the corrected code and the result.

CODE

import vtk
vtp_filename = 'surfacemesh.vtp'
vti_filename = 'volume.vti'
## Read vtkPolyData
reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(vtp_filename)
reader.Update()
surface = reader.GetOutput()
# Read vtkImage
reader = vtk.vtkXMLImageDataReader() 
reader.SetFileName(vti_filename)
reader.Update()
image = reader.GetOutput()
# Gaussian kernel
gaussian_kernel = vtk.vtkGaussianKernel()
gaussian_kernel.SetSharpness(2)
gaussian_kernel.SetRadius(12)
# Set vtkPointInterpolator
interpolator = vtk.vtkPointInterpolator()
interpolator.SetInputData(surface)
interpolator.SetSourceData(image)
interpolator.SetLocator(vtk.vtkMergePoints())
interpolator.SetKernel(gaussian_kernel)
interpolator.Update()
# Write result
writer = vtk.vtkXMLPolyDataWriter()
writer.SetInputData(interpolator.GetOutput())
writer.SetFileName('interpolated.vtp')
writer.Update()

RESULT

vtkPointInterpolator_result
vtkPointInterpolator_result939×503 40.5 KB

6

Good to hear that it works. It seems that the vtkStaticPointLocator cannot be used on vtkImageData, and that’s why a different locator must be used.

Note that if you are using VTK 9 and need extra speed, the vtkImageProbeFilter is specifically designed to efficiently interpolate points from an image (and it’s SMP, too, like vtkPointInterpolator).

interpolator = vtk.vtkImageProbeFilter()
interpolator.SetInputData(surface)
interpolator.SetSourceData(image3D)
interpolator.Update()

It uses linear interpolation by default, but it can also do higher order interpolation.

1 Like
7

Many thanks, I’tried the vtkImageProbeFilter and it’s really much faster than the vtkStaticPointLocator on vtkImageData :+1:

8

@dgobbi thank you, this helped me to solve similar problem.

I’m interpolating custom object inherited from vtkDataSet and I’m wondering is it possible to somehow tell the interpolator that it should look for the points in some order?

My dataset has some partly predictable structure and I would like to experiment with perfomance (the dataset may take pretty much of memory).

Is it vtkPointLocator I should try to look at and probably inherit from it? Or vtkDataSet::GetPoint()/GetCell() are the only methods I should improve?

9

Hi Kerim, if your data set has cells in addition to points, then try using vtkProbeFilter to do the interpolation. If it’s necessary to use vtkPointInterpolator instead of vtkProbeFilter, then a custom vtkPointLocator subclass is probably the best way to go (the locators have a lot of interface methods, so it might require a lot of work).

1 Like
10

Thank you for the hints,

I’m facing some problems when using vtkProbeFilter even though it works fine in most cases but sometimes there is an incorrect result.

I’m trying to debug and I set the breakpoint in MyDataset::FindCell() method but the program never goes there. I don’t understand why but according to the documentation it should.

What maybe the reason of that?

  vtkNew<vtkProbeFilter> interpolator;
  interpolator->SetFindCellStrategy(nullptr);
  interpolator->SetCellLocatorPrototype(nullptr);
  interpolator->SetSourceData(raisedGrid); // this is my dataset
  interpolator->SetInputData(imageData); // I want to interpolate my dataset `raisedGrid` on `vtkImageData`
  interpolator->Update();
11

I have looked to the source code of vtkProbeFilter and understood that when it somehow works with vtkImageData it doesn’t use vtkDataSet::FindCell() but instead it implements some more effective logic.
image