problem with interpolating vtkImageData on vtkPolyData with vtkPointInterpolator

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!

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

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)

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

Try using a different locator to see if it works:

interpolator.SetLocator(vtk.vtkMergePoints())

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

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.

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