vtkContourTriangulator (Triangulation failed)

[marcomusy]

I found an example where the filter seems to fail (but not clear to me why) - I wonder if there is a fix…:

from vedo import *

pts = np.array([
    [1.96, 2.48, -1.34],
    [1.89, 2.34, -1.33],
    [1.74, 2.41, -1.35],
    [1.68, 2.26, -1.35],
    [1.52, 2.32, -1.34],
    [1.59, 2.48, -1.34],
    [1.66, 2.63, -1.33],
    [1.81, 2.56, -1.34],
])
# pts[:,2] = -1.33 # this fixes it
# pts[:, [1,2]] = pts[:, [2,1]] # no effect in any case

line = Line(pts)
# line.write("line.vtk", binary=False)
mesh = line.clone().triangulate().lw(1).c("blue4")
labels = line.labels2d(c="yellow4", scale=2)
show(line, mesh, labels, bg="gray1", axes=1).close()

works fine if:
pts[:,2] = -1.33

fails without it (non-planar):

with error message:

2023-04-04 14:07:57.010 (   0.124s) [        C1F89180]vtkContourTriangulator.:91     ERR| vtkContourTriangulator (0x26bc1a0): Triangulation failed, output might have holes.

Pure VTK to reproduce the issue:

import vtk

reader = vtk.vtkPolyDataReader()
reader.SetFileName("line_bad.vtk")
reader.Update()
polydata = reader.GetOutput()

tf = vtk.vtkContourTriangulator()
tf.TriangulationErrorDisplayOn()
tf.SetInputData(polydata)
tf.Update()
tfpolydata = tf.GetOutput()

mapper = vtk.vtkPolyDataMapper()
mapper.SetInputData(tfpolydata)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().EdgeVisibilityOn()

ren = vtk.vtkRenderer()
ren.AddActor(actor)
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.Start()

lines.zip (699 Bytes)

@dgobbi

Related to:

create separated regions from polydata

I can reproduce what the vedo code from @marcomusy has proposed with the VTK API. But the split function from the vedo code gives anyway an uncorrect number of regions. I have prepared 3 vtk files to play with. https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_11.vtk [image] https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_07.vtk [image] https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_03.vtk [image] In t…

[dgobbi (David Gobbi)]

I have confirmed this locally. As you noted in your fix, your points do not lie in a plane. The vtkContourTriangulator was not designed for inputs that are not flat. The description of this filter is: Fill all 2D contours to create polygons.

The failure seems to be at the very last stage. One thing that I could do in order to make the filter more robust, is to use the normal to project the points onto a plane before doing the final triangulation. It would increase memory usage and the execution time, so it would have to be an option, not the default.

[PBrockmann (Patrick Brockmann)]

It seems that vtkContourTriangulator() does not handle 3D coordinates when non planar.

The big picture of this problem is how to get a anti-counter clockwise coordinates of polydata made of lines sometimes folded.

See: create separated regions from polydata

[marcomusy]

I’ve used this class a lot and never encountered this kind of problem, yet would be nice to have such an option to make it more robust as you suggest! Maybe it would be also a good idea to have the TriangulationErrorDisplayOn() on by default.
Thanks!

[marcomusy]

indeed projecting on the fitting plane makes it work fine:

from vedo import *

def triangulate_non_planar(line):
    data = line.points()
    dmean = data.mean(axis=0)
    pts = data - dmean
    vv = np.linalg.svd(pts)[2]
    n = np.cross(vv[0], vv[1])
    d = np.dot(pts, n)
    pts -= n * d[:, np.newaxis]
    faces = Line(pts).triangulate().faces()
    return Mesh([data, faces])

pts = np.array([
    [1.96, 2.48, -1.34],
    [1.89, 2.34, -1.33],
    [1.74, 2.41, -1.35],
    [1.68, 2.26, -1.35],
    [1.52, 2.32, -1.34],
    [1.59, 2.48, -1.34],
    [1.66, 2.63, -1.33],
    [1.81, 2.56, -1.34],
])
line = Line(pts)
mesh = triangulate_non_planar(line).lw(1).c("blue4")
labels = line.labels2d(c="yellow4", scale=2)
show(line, mesh, labels, bg="gray1", axes=1).close()

[PBrockmann (Patrick Brockmann)]

Ok. Indeed it works with this set of points.
But not sure to be able to find a fitting plane with case like captured below.

My purpose is to find all contigous lines that form convex polylines and describe their coordinates anticlockwise.

[PBrockmann (Patrick Brockmann)]

Here is another very basic example.
Described in Sort points on line in PolyData - #4 by PBrockmann