St_delaunaytriangles polygon with holes

Question

I'm using PostGIS to decompose a polygon with a few holes into triangles via st_delaunaytriangles. Unfortunately, it does not return the right result with holes, because the triangles being outputted cross the section with holes (I want the holes to be empty, i.e. no triangle should enter any of these holes).

I attempted to split the polygons with a line segment to create two polygons without a hole, but this did not work as expected.

Is there a way to decompose a polygon with holes into triangles, while forcing the triangles to not enter the holes?

Answer 1

The ST_DelaunayTriangles function is based on points in the input and doesn't look at the already available linework in your geometry. Hence, it will not only ignore the holes, the resulting triangles may also cross your original polygon. What you are looking for is a constrained delaunay triangulation. Luckily, there is a hidden functionality in SFCGAL that does precisely that, and faster than the geos version. Make sure you get postgis including sfcgal and add the functionality with:

CREATE OR REPLACE FUNCTION public.st_triangulate2dz(geometry)
  RETURNS geometry AS
'$libdir/postgis-2.3', 'sfcgal_triangulate'
  LANGUAGE c IMMUTABLE STRICT
  COST 100;

(make sure you set the correct postgis version).

Please note: the constrained delaunay triangulation will still fill-up the interior and exterior of the polygon, which means you need to remove the non-overlapping polygons after triangulation.

Answer 2

I can't fit all of this in a comment, so I'm adding it as an addendum to @tilt's answer.

A working example using st_triangulate2dz to reproduce the constrained Delaunay triangulation (CDT) example on the Triangle website:

WITH a_collection_to_triangulate AS
(SELECT st_geometryfromtext('MultiPolygon (((0.8 -0.7924, 0.5792 -0.7924, 0.5792 -0.7764, 0.6216 -0.7716, 0.6336 -0.7628, 0.6392 -0.7444, 0.6208 -0.6844, 0.5872 -0.6044, 0.3608 -0.6044, 0.3192 -0.7068, 0.312 -0.7396, 0.3184 -0.7612, 0.3344 -0.7716, 0.3712 -0.7764, 0.3712 -0.7924, 0.2 -0.7924, 0.2 -0.7764, 0.22 -0.7732, 0.2456 -0.7564, 0.2776 -0.702, 0.4888 -0.2076, 0.5048 -0.2076, 0.7408 -0.7396, 0.756 -0.7612, 0.7744 -0.7724, 0.8 -0.7764, 0.8 -0.7924),(0.4736 -0.3308, 0.3744 -0.57, 0.5744 -0.57, 0.4736 -0.3308)))') AS geom
),
b_triangulation_dumped_rebuillt_as_simple AS --Triangulate then dump and rebuild as "simple features" to avoid errors with other ST functions
(SELECT (st_dump(st_makepolygon(st_exteriorring((st_dump(st_triangulate2dz(ST_makevalid(geom)))).geom)))).geom AS geom
FROM a_collection_to_triangulate)
--Constrained, exterior removed triangulation
SELECT b_triangulation_dumped_rebuillt_as_simple.geom
FROM b_triangulation_dumped_rebuillt_as_simple
INNER JOIN a_collection_to_triangulate
ON st_intersects(b_triangulation_dumped_rebuillt_as_simple.geom, a_collection_to_triangulate.geom)
AND NOT st_touches(b_triangulation_dumped_rebuillt_as_simple.geom, a_collection_to_triangulate.geom)