Complex Concave Polygons

Gerald · March 20, 2020, 9:26am

Hi all,

I am dealing with complex concave polygons. My polygons are:

concave
self-touching (either point/point or point/edge, edge/edge) at least once, most of the time multiple
never self-intersecting
aligned on a plane (xy, z=CONST)

I need to create a 3D structure by doing a linear extrude. I can do a extrusion for the polygon outline by adding all the points to a line and pass all the data as a polydata object to the vtkLinearExtrusionFilter.

Using either CappingOn() or using SetPolys() end up in a corrupted top/bottom surface.

I tried vtkTriangleFilter and vtkDelaunay2D without success.

I can try to create a minimal data set (the original polygons have about >1000 edges) to reproduce the behavior (+sample code), if needed.

I did the following test outside of VTK (using octave):

Just take all points
Do a Delaunay triangulation on (1)
Check wether the center point of each created triangle is inside the original polygon (skip if not)
Plot resulting surface
This works for all my test data.

I also tried to mesh the original polygon using gmsh and netgen. While gmsh failed in all cases, netgen sometimes could create a mesh (polygon has only one intersection in this case)

Loading the polygon into FreeCAD and do an extrude from sketch succeeded for all of my testdata.

Any idea? Best regards,

Gerald

lassoan · March 20, 2020, 12:27pm

vtkDelaunay2D can robustly mesh 2D concave polygons. Make sure you don’t reuse any point IDs (duplicate points at self-touching positions and use each point ID only once).

Gerald · March 21, 2020, 8:45am

Hi Andras,

thanks for the quick reply. I am still facing the same problem. To be more clear about the initial issue:

The meshing (any algo) idea comes from extrusion problems with complex polygons. I tried the ‘standard’ way of using vtkTriangleFilter, but ended up with malformed/empty results. I created a minimal test setup with 3 example polygons working with vtkTriangleFilter and 1 failing example (wrong results).

Goal would be to create a linear extrusion of an arbitrary complex input polygon.

I need to inline the test program (no uploads allowed for new users)

from vtk import *

define test polygon ---------------------------------------------------------

x=[0,6,6,4,4,5,5,2,2,0] # concave no touch/intersect

y=[0,0,6,6,4,4,2,2,6,6]

x=[0,6,6,2,4,5,5,2,2,0] # concave touch at 1 point/no intersect

y=[0,0,6,6,4,4,2,2,6,6]

x=[0,6,6,2,2,5,5,2,2,0] # concave touch at one edge/no intersect

y=[0,0,6,6,4,4,2,2,6,6]

20 4 3 concave polygon CCW

|17|6 | Self-touching two edges

| |_________| | 4,19 at same coords

|16 | 7 | 5,18 at same coords

| 13| 9 |10 | 9,14 at same coords

| 12|_____|11 |

|_____________|

1 2

x=[0,6,6,3,3,5,5,3,3,4,4,2,2,3,3,1,1,3,3,0]
y=[0,0,6,6,5,5,4,4,3,3,1,1,3,3,4,4,5,5,6,6]

npts = len(x)

creation of the source ------------------------------------------------------

points = vtkPoints()
polygon = vtkPolygon()
polygon.GetPointIds().SetNumberOfIds( npts )

for i in range( npts ):
points.InsertNextPoint( x[i], y[i], 0 ) # xy plane aligned
polygon.GetPointIds().SetId(i, i)

polys = vtkCellArray()
polys.InsertNextCell(polygon)

polyData = vtkPolyData()
polyData.SetPoints(points)
polyData.SetPolys(polys)
result = polyData # set inital result

filter the polydata ---------------------------------------------------------

#pdflt = vtkDelaunay2D()
#pdflt.SetInputData( polyData )
#pdflt.Update()
#result = pdflt.GetOutput()

#pdflt = vtkTriangleFilter()
#pdflt.SetInputData( polyData )
#pdflt.Update()
#result = pdflt.GetOutput()

map it ----------------------------------------------------------------------

mapper = vtkPolyDataMapper()
mapper.SetInputData( result )

render it -------------------------------------------------------------------

actor = vtkActor()
actor.SetMapper(mapper)

render = vtkRenderer()
render.AddActor(actor)

renwin = vtkRenderWindow()
renwin.AddRenderer(render)

rwiact = vtkRenderWindowInteractor()
rwiact.SetRenderWindow(renwin)

rwiact.SetInteractorStyle(vtkInteractorStyleTrackballCamera())
rwiact.Initialize()
rwiact.Start()

Best regards
Gerald

Gerald · March 21, 2020, 8:52am

I removed all comments from the code, as this creates some strange formating:

from vtk import *

x=[0,6,6,4,4,5,5,2,2,0]
y=[0,0,6,6,4,4,2,2,6,6]

x=[0,6,6,2,4,5,5,2,2,0]
y=[0,0,6,6,4,4,2,2,6,6]

x=[0,6,6,2,2,5,5,2,2,0]
y=[0,0,6,6,4,4,2,2,6,6]

x=[0,6,6,3,3,5,5,3,3,4,4,2,2,3,3,1,1,3,3,0]
y=[0,0,6,6,5,5,4,4,3,3,1,1,3,3,4,4,5,5,6,6]

npts = len(x)

points = vtkPoints()
polygon = vtkPolygon()
polygon.GetPointIds().SetNumberOfIds( npts )

for i in range( npts ):
points.InsertNextPoint( x[i], y[i], 0 )
polygon.GetPointIds().SetId(i, i)

polys = vtkCellArray()
polys.InsertNextCell(polygon)

polyData = vtkPolyData()
polyData.SetPoints(points)
polyData.SetPolys(polys)
result = polyData

pdflt = vtkTriangleFilter()
pdflt.SetInputData( polyData )
pdflt.Update()
result = pdflt.GetOutput()

mapper = vtkPolyDataMapper()
mapper.SetInputData( result )

actor = vtkActor()
actor.SetMapper(mapper)

render = vtkRenderer()
render.AddActor(actor)

renwin = vtkRenderWindow()
renwin.AddRenderer(render)

rwiact = vtkRenderWindowInteractor()
rwiact.SetRenderWindow(renwin)

rwiact.SetInteractorStyle(vtkInteractorStyleTrackballCamera())
rwiact.Initialize()
rwiact.Start()

lassoan · March 21, 2020, 2:11pm

You can define a triangle strip cell between the two polygons, fill the two polygons using vtkDelaunay2D, and merge them into a watertight mesh using vtkAppendPolyData.

You need to put three backtick characters (```) before and after source code (otherwise text is interpreted as markdown). You may also indicate programming language to optimize syntax highlighting:

```python
some source code
some more source code
```

Gerald · March 21, 2020, 2:58pm

Thanks again for the fast reply. Good to know about the backtick characters. So the test code may got rendered somehow misleading:

I only have one polygon as the input.

x=[0,6,6,3,3,5,5,3,3,4,4,2,2,3,3,1,1,3,3,0]
y=[0,0,6,6,5,5,4,4,3,3,1,1,3,3,4,4,5,5,6,6]

The triangulation is just (needed?) because of the extrusion done in the 2nd step. If I look at the result of vtkDelaunay2D I find some missing triangles at the outside boundary (hull is not convex).

I also tried to SetSourceData() to create a constrained delaunay. This results in:

Warning: In /build/vtk7-w4DzBd/vtk7-7.1.1+dfsg1/Filters/Core/vtkDelaunay2D.cxx, line 1340
vtkDelaunay2D (0x263da50): Edge not recovered, polygon fill suspect

This is the (now readable ) test code:

#!/usr/bin/env python

from vtk import *

# define test polygon ---------------------------------------------------------

#  20 ______4______ 3             concave polygon CCW
#    |17____|____6 |              Self-touching two edges
#    | |_________| |              4,19 at same coords
#    |16  __|__  7 |              5,18 at same coords
#    | 13|  9  |10 |              9,14 at same coords
#    | 12|_____|11 |
#    |_____________| 
#   1              2

x=[0,6,6,3,3,5,5,3,3,4,4,2,2,3,3,1,1,3,3,0]
y=[0,0,6,6,5,5,4,4,3,3,1,1,3,3,4,4,5,5,6,6]

npts = len(x)

# creation of the source ------------------------------------------------------

points  = vtkPoints()
polygon = vtkPolygon()
polygon.GetPointIds().SetNumberOfIds( npts )

for i in range( npts ):
    points.InsertNextPoint( x[i], y[i], 0 )     # xy plane aligned
    polygon.GetPointIds().SetId(i, i)

polys = vtkCellArray()
polys.InsertNextCell(polygon) 
    
polyData = vtkPolyData()
polyData.SetPoints(points)
polyData.SetPolys(polys)
result = polyData                               # set inital result

# filter the polydata ---------------------------------------------------------

pdflt = vtkDelaunay2D()
pdflt.SetSourceData( polyData )
pdflt.SetInputData( polyData )
pdflt.Update()
result = pdflt.GetOutput()

# info ------------------------------------------------------------------------

numpoints = result.GetNumberOfPoints()
print( '#points = ', numpoints )
numpolys  = result.GetNumberOfPolys()
print( '#polys  = ', numpolys )

# map it ----------------------------------------------------------------------

mapper = vtkPolyDataMapper()
mapper.SetInputData( result )

# render it -------------------------------------------------------------------

actor =  vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetRepresentationToWireframe()

render = vtkRenderer()
render.AddActor(actor)

renwin = vtkRenderWindow()
renwin.AddRenderer(render)

rwiact = vtkRenderWindowInteractor()
rwiact.SetRenderWindow(renwin)

rwiact.SetInteractorStyle(vtkInteractorStyleTrackballCamera())
rwiact.Initialize()
rwiact.Start()

Here is a Octave version (for comparison):

vx=[0,6,6,3,3,5,5,3,3,4,4,2,2,3,3,1,1,3,3,0];
vy=[0,0,6,6,5,5,4,4,3,3,1,1,3,3,4,4,5,5,6,6];

tri = delaunay (vx, vy);

cx = [ vx(tri(:,1)) + vx(tri(:,2)) + vx(tri(:,3)) ] ./ 3;
cy = [ vy(tri(:,1)) + vy(tri(:,2)) + vy(tri(:,3)) ] ./ 3;

[in,on] = inpolygon (cx, cy, vx, vy);

fin = tri(find(in),:);

clf;
triplot (tri, vx, vy);
hold on;
triplot( fin, vx, vy, "c" ); 
plot (vx, vy, "r*");
plot (cx, cy, "g*");

axis( 'equal' );

The octave version is internaly based on http://www.qhull.org/ library. As a second step I simply check if all center points of the generated triangles are within the original polygon. Any others are dropped. I was assuming the vtkDelaunay2D works somehow similar.

Gerald

lassoan · March 21, 2020, 3:09pm

I won’t be able to review your code, but I can give you a working concave curve mesher code in 3D Slicer. The code is more complex than you need, as it allows creation of arbitrary space curves (it does not have to be planar and may be oriented arbitrarily). As you can see in the code, there are a few things to pay attention to (project the curve points to XY plane, ensure correct polygon winding, etc.). If these are taken care of then triangulation works perfectly.

Gerald · March 21, 2020, 4:33pm

Thanks again, and yes it is ‘slighly’ more complex than I need. I was hoping for some kind of single filter usage to get this done. For now I will stuck with some preprocessing (as done by the matlab/octave code) and put the precalculated triangles into the polydata object rather than the inital polygon.

I still think that somethings is wrong with the vtkDelaunay2D (for the case of complex polygons).

Gerald