delaunay_2d fails for (some) polygons with holes

[TysonRayJones]

Describe the bug, what's wrong, and what you expected.

delaunay_2d outputs a warning (Edge not recovered, polygon fill suspect) and fails to correctly filter PolyData encoding a basic polygon with two holes.

More precisely, I create a PolyData from a set of counterclockwise outer coordinates and clockwise inner coordinates, which together encode a polygon with two holes. I call delaunay_2d upon this data, but receive VTK warning

2023-03-07 13:59:08.217 (   0.288s) [                ]      
vtkDelaunay2D.cxx:1419  WARN| vtkDelaunay2D (0000016829DFAE40): 
Edge not recovered, polygon fill suspect

Viewing the resulting polydata (with .plot()) reveals an incorrect result:

The expected and correct result (rendered with matplotlib) resembles:

This error does not occur with all holed polygons. Indeed, PyVista filters the above geometry correctly if one of the holes is removed. In that case, it respectively yields:

and

Several other holed polygons (with both convex, complex, and general polygon holes) filter and plot totally fine, as per this SO thread.

Steps to reproduce the bug.

Simply run this and behold the VTK warning and erroneous plot:

import pyvista

outer_xy = [(14, 0), (14, 7), (0, 7), (0, 0)]
holes_xy = [
    [(6, 1), (6, 2), (7, 2), (7,3), (8,3), (8,2), (9,2), (9,1)],
    [(3, 1), (3, 5), (5, 5), (5, 1)]
]

def get_polydata(xy, z):
    xyz = [(x,y,z) for x,y in xy]
    faces = [len(xyz), *range(len(xyz))]
    data = pyvista.PolyData(xyz, faces=faces)
    return data

data = get_polydata(outer_xy, z=0)
for hole_xy in holes_xy:
    data += get_polydata(hole_xy, z=0)

polygon = data.delaunay_2d(edge_source=data) # error

polygon.plot(cpos=[(6.5, 4, 35), (7, 3.5, 0), (0, 1, 0)], show_edges=True)

System Information

Date: Tue Mar 07 14:34:25 2023 GMT Standard Time

                OS : Windows
            CPU(s) : 20
           Machine : AMD64
      Architecture : 64bit
               RAM : 15.7 GiB
       Environment : Python
        GPU Vendor : Intel
      GPU Renderer : Intel(R) Iris(R) Xe Graphics
       GPU Version : 4.5.0 - Build 31.0.101.3301

  Python 3.10.10 (tags/v3.10.10:aad5f6a, Feb  7 2023, 17:20:36) [MSC v.1929 64
  bit (AMD64)]

           pyvista : 0.38.3
               vtk : 9.2.6
             numpy : 1.23.4
           imageio : 2.26.0
            scooby : 0.7.1
             pooch : v1.7.0
        matplotlib : 3.6.2
           IPython : 8.11.0
        ipywidgets : 8.0.4
             scipy : 1.10.1
             trame : 2.3.2
      trame_client : 2.7.3
      trame_server : 2.9.1
         trame_vtk : 2.2.1
      nest_asyncio : 1.5.6

Screenshots

(screenshots embedded above)

[TysonRayJones]

As pointed out by Andras on this SO question, the error can be avoided by introducing additional interpolating vertices:

def interpolate_polygon(xy):
    new_xy = []
    for i in range(len(xy)):
        a, b = xy[i], xy[(i+1)%len(xy)]
        for coord in numpy.linspace(a, b, num=5):
            new_xy.append(coord)

    return new_xy

outer_xy = interpolate_polygon(outer_xy)
for i in range(len(holes_xy)):
    holes_xy[i] = interpolate_polygon(holes_xy[i])

For the MWE above, this yields correct result:

This is a somewhat unsatisfying solution however, since I do not a priori know how many interpolating points (above, 5 assumed) are necessary, and my application is a user-facing package. Further, always interpolating will slow down an interactive visualisation in my application.

Finally, this solution actually worsens the problem for my more complicated non-MWE geometries, causing filtering to fail when it previously worked for non-interpolated polygons!

[darikg]

Is using triangle feasible? I've found it to be very robust

import pyvista as pv
from triangle import triangulate
import numpy as np

outer_xy = [(14, 0), (14, 7), (0, 7), (0, 0)]

holes_xy = [
    [(6, 1), (6, 2), (7, 2), (7,3), (8,3), (8,2), (9,2), (9,1)],
    [(3, 1), (3, 5), (5, 5), (5, 1)]
]

def edge_idxs(nv):
    i = np.append(np.arange(nv), 0)
    return np.stack([i[:-1], i[1:]], axis=1)

nv = 0
verts, edges = [], []
for loop in (outer_xy, *holes_xy):
    verts.append(loop)
    edges.append(nv + edge_idxs(len(loop)))
    nv += len(loop)

verts, edges = np.concatenate(verts), np.concatenate(edges)

# Triangulate needs to know a single interior point for each hole
# Using the centroid works here, but for very non-convex holes may need a more sophisticated method,
# e.g. shapely's `polylabel`
holes = np.array([np.mean(h, axis=0) for h in holes_xy])

# Because triangulate is a wrapper around a C library the syntax is a little weird, 'p' here means planar straight line graph
d = triangulate(dict(vertices=verts, segments=edges, holes=holes), opts='p')

# Convert back to pyvista
v, f = d['vertices'], d['triangles']
nv, nf = len(v), len(f)
points = np.concatenate([v, np.zeros((nv, 1))], axis=1)
faces = np.concatenate([np.full((nf, 1), 3), f], axis=1).reshape(-1)

pv.PolyData(points, faces=faces).plot(show_edges=True)

[adeak]

# Triangulate needs to know a single interior point for each hole
# Using the centroid works here, but for very non-convex holes may need a more sophisticated method,

I think simply connected polygons work fine with VTK, in which case one could take each hole, define a polydata face for it, triangulate() and pick out an arbitrary cell and use its center.

[TysonRayJones]

I really appreciate the potential workarounds, thanks very much! Unfortunately I'm bound to using pyvista, where this seems a genuine bug (or at least, one in VTK beneath).