INTERPOLATE ON STL: Plotting scalar field on 50 points on An 70000-points STL surface

Hello,

I have an STL surface file. The data were sampled on 50 selected points only
(wind tunnel test).
I want to visualize the field as a contour on the STL file.

I know basic stuff on how to use VTK python library. But this task I want it
involves:

  • interpolating the field inside the probes-area
  • extrapolating outside the probes-area
  • then creating the final VTK file.
  • maybe triangulation is also needed…

OR, there might be a crazy filter in paraview that can do the whole thing!

Any help would be appreciated…

Can you share the file?

Perhaps the vtki Python package would be a goodd place to start:

import vtki
data = vtki.read(‘my_file.stl’)

tri = data.tri_filter()
contours = tri.contour()
contours.plot()

contours.save(‘my_new_output.vtk’)
1 Like

Hello Bane,

Thank you very much for your reply. I could not attach the files (new user) but here is a dropbox link:

You will find an STL file (mesh only) and a txt file that contains the coordinates of the points and the value of the field at each point: [x, y, z, val].

What I need is to interpolate this field on the STL and visualize it as a VTK on paraview.

The script you shared does not interpolate the field I have on the STL, does it?

Thanks,

There are a lot of different ways you could resample/interpolate your sparse points onto your surface. Here are a few using numpy, scipy and vtki but you probably want to check out libraries that implement kriging.

import vtki
import numpy as np
# Load STL file
surface = vtki.read('InterpolatingOnSTL_final.stl')
# Load numpy file
foo = np.loadtxt('points.txt', skiprows=1, delimiter='\t')
lookup = vtki.PolyData(foo[:, 0:3])
lookup.point_arrays['val'] = foo[:,3]
# Inspect the given data
p = vtki.Plotter()
p.add_mesh(surface)
p.add_mesh(lookup, render_points_as_spheres=True, point_size=10)
p.show()

# Use scipy to run interpolation
from scipy.interpolate import griddata

method = 'linear'
surface.point_arrays[method] = griddata(
    lookup.points, lookup.point_arrays['val'],
    surface.points, method=method)

method = 'nearest'
surface.point_arrays[method] = griddata(
    lookup.points, lookup.point_arrays['val'],
    surface.points, method=method)

# Assisted linear
assisted = surface.point_arrays['linear'].copy()
badi = np.argwhere(np.isnan(assisted))
assisted[badi] = surface.point_arrays['nearest'][badi]
surface.point_arrays['assisted linear'] = assisted
# Plot the nearest neighbor interpolated data
p = vtki.Plotter()
p.add_mesh(surface, scalars='nearest')
p.add_mesh(lookup, render_points_as_spheres=True, 
           point_size=10, color='w')
p.show()

# Plot the linear interpolated data
p = vtki.Plotter()
p.add_mesh(surface, scalars='linear')
p.add_mesh(lookup, render_points_as_spheres=True, 
           point_size=10, color='w')
p.show()

# Plot the assisted linear interpolated data
p = vtki.Plotter()
p.add_mesh(surface, scalars='assisted linear')
p.add_mesh(lookup, render_points_as_spheres=True, 
           point_size=10, color='w')
p.show()

1 Like

These toolsets looks promising:

1 Like

The vtkPointInterpolator with a Gaussian Kernel (or other kernel) is an alternative method to interpolate and extrapolate a bit more smoothly.

import vtk
import numpy as np

points_reader = vtk.vtkDelimitedTextReader()
points_reader.SetFileName('points.txt')
points_reader.DetectNumericColumnsOn()
points_reader.SetFieldDelimiterCharacters('\t')
points_reader.SetHaveHeaders(True)

table_points = vtk.vtkTableToPolyData()
table_points.SetInputConnection(points_reader.GetOutputPort())
table_points.SetXColumn('x')
table_points.SetYColumn('y')
table_points.SetZColumn('z')
table_points.Update()

points = table_points.GetOutput()
points.GetPointData().SetActiveScalars('val')
range = points.GetPointData().GetScalars().GetRange()

# Read a probe surface
stl_reader = vtk.vtkSTLReader()
stl_reader.SetFileName('InterpolatingOnSTL_final.stl')
stl_reader.Update()

surface = stl_reader.GetOutput()
bounds = np.array(surface.GetBounds())

dims = np.array([101, 101, 101])
box = vtk.vtkImageData()
box.SetDimensions(dims)
box.SetSpacing((bounds[1::2] - bounds[:-1:2])/(dims - 1))
box.SetOrigin(bounds[::2])

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

interpolator = vtk.vtkPointInterpolator()
interpolator.SetInputData(box)
interpolator.SetSourceData(points)
interpolator.SetKernel(gaussian_kernel)

resample = vtk.vtkResampleWithDataSet()
resample.SetInputData(surface)
resample.SetSourceConnection(interpolator.GetOutputPort())

mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(resample.GetOutputPort())
mapper.SetScalarRange(range)

actor = vtk.vtkActor()
actor.SetMapper(mapper)

point_mapper = vtk.vtkPointGaussianMapper()
point_mapper.SetInputData(points)
point_mapper.SetScalarRange(range)
point_mapper.SetScaleFactor(0.6)
point_mapper.EmissiveOff();
point_mapper.SetSplatShaderCode(
    "//VTK::Color::Impl\n"
    "float dist = dot(offsetVCVSOutput.xy,offsetVCVSOutput.xy);\n"
    "if (dist > 1.0) {\n"
    "  discard;\n"
    "} else {\n"
    "  float scale = (1.0 - dist);\n"
    "  ambientColor *= scale;\n"
    "  diffuseColor *= scale;\n"
    "}\n"
)

point_actor = vtk.vtkActor()
point_actor.SetMapper(point_mapper)

renderer = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(renderer)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

renderer.AddActor(actor)
renderer.AddActor(point_actor)

iren.Initialize()

renWin.Render()
iren.Start()

PointInterpolator_Gaussian

2 Likes

Bane & Kenichiro,

Wow, thank you very much for the quick amazing solutions to this. Exactly what I needed!

Thank you very much!
Hosam

@Kenichiro-Yoshimi

Would you mind if I include a C++ version of your excellent solution in the VTKExamples Project (https://lorensen.github.io/VTKExamples/site/).
@amaclean Perhaps you can add the python version if Kenichiro gives us permission.

@Hosam
Can I add your data to the VTKExamples Project?

Bill

Bill, I would love to add it.

@lorensen
I submitted a pull request to add C++ example based on this code.

Yes of course you can use my example.

Hello Kenichiro,

Thank you very much again for solving my issue.
I have a couple of questions about the code.

Most importantly: how to save the result as a vtk to open it later in paraview?
Here is an example of the VTK files I am used to:

Also:

  • Why you have used 101 in the following:
    dims = np.array([101, 101, 101])

  • Also 2, 12 in:
    gaussian_kernel.SetSharpness(2)
    gaussian_kernel.SetRadius(12)

Could you please refer me to the reference where I can choose the type of interpolating kernel e.g. cubic,…etc.

  • The choice of 0.6 in point_mapper.SetScaleFactor(0.6)

  • The choice of 1.0 in:
    “if (dist > 1.0) {\n”
    " discard;\n

Thank you very much in advance!

Hello Bane,

Thank you very much for introducing this solution. I am interested also in the nearest neighbor interpolation to show the areas influenced directly by each probe.
Could you please show me how to save this results as a VTK file. Here is an example of the VTK files I am used to:

Thank you!
Hosam

Hello Hosam,

To save the result in legacy vtk format(*.vtk), you can usually use the vtkDataSetWriter.

I previously interpolated point data on a surface through 3d image data, but it is really not necessary. So I rewrite the previous code exclusive of image data as below and you can ignore the corresponding part of vtkImageData.

Their values can be set as you prefer. Also there are really various interpolating kernels. See the links below and try:
https://vtk.org/gitweb?p=VTK.git;a=blob;f=Filters/Points/Testing/Python/TestEllipsoidalGaussianKernel.py
https://vtk.org/gitweb?p=VTK.git;a=blob;f=Filters/Points/Testing/Python/TestPointInterpolator.py

import vtk

points_reader = vtk.vtkDelimitedTextReader()
points_reader.SetFileName('points.txt')
points_reader.DetectNumericColumnsOn()
points_reader.SetFieldDelimiterCharacters('\t')
points_reader.SetHaveHeaders(True)

table_points = vtk.vtkTableToPolyData()
table_points.SetInputConnection(points_reader.GetOutputPort())
table_points.SetXColumn('x')
table_points.SetYColumn('y')
table_points.SetZColumn('z')
table_points.Update()

points = table_points.GetOutput()
points.GetPointData().SetActiveScalars('val')
range = points.GetPointData().GetScalars().GetRange()

# Read a probe surface
stl_reader = vtk.vtkSTLReader()
stl_reader.SetFileName('InterpolatingOnSTL_final.stl')
stl_reader.Update()

surface = stl_reader.GetOutput()

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

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

writer = vtk.vtkDataSetWriter()
writer.SetInputConnection(interpolator.GetOutputPort())
writer.SetFileName('interpolated.vtk')
writer.Write()
1 Like

@Hosam: vtki objects are instances of VTK objects, so you could save them in the same way you would save any VTK data object.

To make life a bit easier, we have added a .save() method on each of the objects:

...
# after the code I shared in my previous post
surface.save('my-data.vtk')

Or you could pass that object to a VTK pipeline if you’re more familiar with that:

...
# after the code I shared in my previous post
import vtk
writer = vtk.vtkDataSetWriter()
writer.SetInputDataObject(surface)
writer.SetFileName('my-data.vtk')
writer.Write()
1 Like

@Kenichiro-Yoshimi I have added your Python example to VTKExamples.
It will appear as PointInterpolator when the web pages are updated.

cc: @lorensen

@amaclean
Thank you very much!