# Complex Concave Polygons

**URL:** https://discourse.vtk.org/t/complex-concave-polygons/2857
**Category:** Support
**Created:** [March 20, 2020, 9:26am UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857 "2020-03-20T09:26:09Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![Gerald](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/g/ea666f/32.png) [@Gerald](https://discourse.vtk.org/u/Gerald)
#### Post date: [March 20, 2020, 9:26am UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/1 "2020-03-20T09:26:09Z")

</div>

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):

1. Just take all points
2. Do a Delaunay triangulation on (1)
3. Check wether the center point of each created triangle is inside the original polygon (skip if not)
4. 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

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [March 20, 2020, 12:27pm UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/2 "2020-03-20T12:27:49Z")

</div>

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).

---

<div class="post-metadata">

### Author: ![Gerald](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/g/ea666f/32.png) [@Gerald](https://discourse.vtk.org/u/Gerald)
#### Post date: [March 21, 2020, 8:45am UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/3 "2020-03-21T08:45:54Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Gerald](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/g/ea666f/32.png) [@Gerald](https://discourse.vtk.org/u/Gerald)
#### Post date: [March 21, 2020, 8:52am UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/4 "2020-03-21T08:52:11Z")

</div>

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()

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [March 21, 2020, 2:11pm UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/5 "2020-03-21T14:11:39Z")

</div>

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`  
`````

---

<div class="post-metadata">

### Author: ![Gerald](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/g/ea666f/32.png) [@Gerald](https://discourse.vtk.org/u/Gerald)
#### Post date: [March 21, 2020, 2:58pm UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/6 "2020-03-21T14:58:08Z")

</div>

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:

```python
#!/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):

```auto
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/](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

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [March 21, 2020, 3:09pm UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/7 "2020-03-21T15:09:32Z")

</div>

I won’t be able to review your code, but I can give you a working [concave curve mesher code in 3D Slicer](https://github.com/Slicer/Slicer/blob/1817652c93e6e88df61defca41a87490d1376275/Modules/Loadable/Markups/Logic/vtkSlicerMarkupsLogic.cxx#L1410-L1533). 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.

---

<div class="post-metadata">

### Author: ![Gerald](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/g/ea666f/32.png) [@Gerald](https://discourse.vtk.org/u/Gerald)
#### Post date: [March 21, 2020, 4:33pm UTC](https://discourse.vtk.org/t/complex-concave-polygons/2857/8 "2020-03-21T16:33:54Z")

</div>

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