Point colour grid interpolation on a vtkPolyData

Julien · January 30, 2024, 2:27am

Hello,

I am trying to cut through a geometry which is the result of a finite element analysis.

Here is what I do (just a bit of context):

plane= vtkPlane()
clipper= vtkClipDataSet()
clipper.setClipFunction(plane)
clipper.setInputData(actorMapper.getInputAsDataSet())

clipperActorMapper= vtkDataSetMApper()
clipperActorMapper.setInputConnection(clipper.getOutputPort();
clippedActor= vtkActor()
clippedActor.SetMapper(clipperActorMapper);
ren.AddActor(clippedActor);

Then a I am using a vtkImplicitPlaneWidget2 to place my cutting plane and actually cut the geometry, it’s working well with its callback function

implicitPlaneWidget = vtkImplicitPlaneWidget2()
implicitPlaneWidget.AddObserver("InteractionEvent", cb, "planeUpdate")

Then in the callback function I am using a vtkFeatureEdges to get the boundary edges of the cut geometry, then I am using a vtkStripper() to get the polylines, then I am using a vtkContourTriangulator() (which works well even for “concave contour” on the flat surfaces) to get a vtkPolyData which is a meshed flat surface of the cut surface.

I can then use a mapper/actor and display it in a solid colour.

Now, the issue.
I am trying to interpolate the colours onto the cut surface from the stress values of the nearby Gauss points of the FEA results.

I am using the PointInterpolator example as a basis for this.
https://examples.vtk.org/site/Python/Meshes/PointInterpolator/
I will probably have to use another interpolator than the Gaussian one, but that is not the current issue here.

So I took my surface from my vtkContourTriangulator and adapted the PointInterpolator example.
I was confident, but as is, it’s actually quite poor.
To better show you the issue, I saved the surface to a *.vtk file, and wrote a minimal test case adapted from the PointInterpolator example, with a rectangular surface and only three points (and a vtkLookupTable from green to red)

After a while, I realized the issue is the colours are first interpolated from the Gauss points (here my three sample points) to be applied to the triangles vertex from the surface. Then, the colours for each triangle vertex are interpolated to fill the triangles area.
So, as the triangles are so flat and long with the vtkContourTriangulator(), the interpolations in each triangles are not what is expected at all.
It didn’t show on the PointInterpolator example because it’s a regular mesh with small triangles, so it’s good enough, but what I would need is some sort of direct mapping to the triangles from the Gauss points.
Or small nice triangles to make it look like it’s all right.

I tried to remesh the surface with a vtkDelaunay2D, and the mapping is already much better on the example (because there is single row of Gauss points) , but not so much on my real models where the mesh would need to be refined to get much smaller triangles, plus it totally fails on a complex concave surface anyway, so … Instead of remeshing with the same points, I would need to keep the contour and refine the triangles.

Anyway, here are my questions:

does it look like it’s more or less the way to go to do what I need? (interpolate the Gauss points onto my surface) or is there something better in vtk already?
is it possible to actually interpolate the colours directly to the triangles directly relative to the distance of the Gauss points?
if it’s not possible to interpolate the colours correctly directly to the triangles area, is there a refining class somewhere in vtk to get proper nice equilateral-ish triangles with a max size?

Thank you very much in advance!

Not really necessary, but if someone is interested, here is the full test case, with the content of the triangulatedSurface.vtk content attached to the post
triangulatedSurface.vtk (9.7 KB)

#!/usr/bin/env python
from typing import List

import numpy as np
# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonDataModel import vtkImageData, vtkPolyData, vtkPlane, vtkQuad, vtkPointData
from vtkmodules.vtkFiltersCore import vtkResampleWithDataSet, vtkPolyDataPlaneCutter
from vtkmodules.vtkFiltersGeneral import vtkTableToPolyData
from vtkmodules.vtkFiltersPoints import (
    vtkGaussianKernel,
    vtkPointInterpolator, vtkLinearKernel
)
from vtkmodules.vtkIOGeometry import vtkSTLReader
from vtkmodules.vtkIOInfovis import vtkDelimitedTextReader
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkPointGaussianMapper,
    vtkPolyDataMapper,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)

import vtk



def main():
    
    pdReader: vtk.vtkPolyDataReader= vtk.vtkPolyDataReader()
    pdReader.SetFileName("d:/triangulatedSurface.vtk")
    pdReader.Update()
    surfaceStl: vtkPolyData= pdReader.GetOutput()
    
    
    

    #
    remesh: bool= False
    if remesh:
        delaunay: vtk.vtkDelaunay2D= vtk.vtkDelaunay2D()
        delaunay.SetInputData(surfaceStl)
        delaunay.Update()
        surfaceStl= delaunay.GetOutput()
    #
    
    
    
    
    
    
    

    # Points
    points= vtk.vtkPoints()
    points.InsertNextPoint(2.0, 00.0, 0.0)
    points.InsertNextPoint(2.0, 10.0, 0.0)
    points.InsertNextPoint(2.0, 20.0, 0.0)


    # Scalars
    scalars = vtk.vtkFloatArray()
    scalars.SetNumberOfComponents(1)
    scalars.InsertNextValue(00.0)
    scalars.InsertNextValue(30.0)
    scalars.InsertNextValue(60.0)
    

    pd: vtkPolyData= vtkPolyData()
    pd.SetPoints(points)
    pd.GetPointData().SetScalars(scalars)

    
    
    valuesRange: List[float]= scalars.GetRange()









    # #########################################################################
    # below this, it's the original example (besides a vtkLookupTable addition)
    

    bounds = np.array(surfaceStl.GetBounds())

    dim= 11
    dims = np.array([dim, dim, dim])
    box: vtkImageData = vtkImageData()
    box.SetDimensions(dims)
    boxmax: np.ndarray= bounds[1::2]
    boxmin: np.ndarray= bounds[:-1:2]
    boxRange: np.ndarray= (boxmax - boxmin)
    box.SetSpacing(boxRange/ (dims - 1))
    boxOrigin: np.ndarray=bounds[::2]
    box.SetOrigin(boxOrigin)



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


    interpolator = vtkPointInterpolator()
    interpolator.SetInputData(box)
    interpolator.SetSourceData(pd)
    interpolator.SetKernel(gaussian_kernel)
    

    resample = vtkResampleWithDataSet()
    resample.SetInputData(surfaceStl)
    resample.SetSourceConnection(interpolator.GetOutputPort())
    
    

    mapperSurface: vtkPolyDataMapper= vtkPolyDataMapper()
    mapperSurface.SetInputConnection(resample.GetOutputPort())
    mapperSurface.SetScalarRange(valuesRange)
    
    actorSurface: vtkActor = vtkActor()
    actorSurface.SetMapper(mapperSurface)







    point_mapper: vtkPointGaussianMapper = vtkPointGaussianMapper()
    point_mapper.SetInputData(pd)
    point_mapper.SetScalarRange(valuesRange)
    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: vtkActor = vtkActor()
    point_actor.SetMapper(point_mapper)




    # lut mod begin
    lut= vtk.vtkLookupTable()
    lut.SetHueRange(0.25, 0.)
    lut.Build()
    mapperSurface.SetLookupTable(lut)
    point_mapper.SetLookupTable(lut)
    # lut mod end
    
    
    
    
    
    

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




    # actorSurface.GetProperty().EdgeVisibilityOn()


    renderer.AddActor(actorSurface)
    renderer.AddActor(point_actor)
    #renderer.AddActor(dactor)
    
    
    colors = vtkNamedColors()
    renderer.SetBackground(colors.GetColor3d('SlateGray'))

    renWin.SetSize(640, 480)
    renWin.SetWindowName('PointInterpolator')

    renderer.ResetCamera()
    #renderer.GetActiveCamera().Elevation(-45)

    iren.Initialize()

    renWin.Render()
    iren.Start()


if __name__ == '__main__':
    main()

Julien · January 30, 2024, 2:51am

Now that I posted my question, I just found VTK: vtkAdaptiveSubdivisionFilter Class Reference and it may look like it could be interesting to refine my mesh.
But, the question remains about the direct colour interpolation on the triangles rather than a mesh refinement anyway. And maybe it won’t work as I would need (we’ll see in a moment/tomorrow), and maybe there is a much better way, anyways!

So, the questions remain and thanks in advance to you guys!

…
…

Edit:
So, by using this refine block instead of the former remesh one

    refine: bool= True
    if refine:
        subdiv: vtk.vtkAdaptiveSubdivisionFilter= vtk.vtkAdaptiveSubdivisionFilter()
        subdiv.SetInputData(surfaceStl)
        subdiv.SetMaximumEdgeLength(1)
        subdiv.Update()
        surfaceStl= subdiv.GetOutput()

Here is what it now looks like. The triangles are still ugly, but as long as the edges are not shown, I guess it could be acceptable on the short term. But it would still be super nice if there were a (much) more elegant solution?

Julien · January 31, 2024, 12:16am

Note for later: maybe have a look at some sort of custom procedural textures?