4

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?

  • It's not possible in simple features, generally and not in this GEOS function. See here for some discussion: rpubs.com/cyclemumner/sc-rationale PostGIS probably has a *constrained (near-)Delaunay triangulation" method, you'll find it in CGAL, and the Triangle library and a few others. Manifold 8.0 calls it Constrained Triangulation (or ToTrianglesAdv). There's an R package sfdct to drive Triangle for the simple features case if you can use that. There are other constrained non-near Delaunay triangulation methods such as ear-clipping. – mdsumner May 23 '17 at 9:52
  • Near duplicate: gis.stackexchange.com/questions/232052/… – mdsumner May 23 '17 at 9:53
  • @mdsumner thank you ! I did post a link below of notes by Ben Discoe. It includes lots of libraries (so far none have handled all my test cases), except for sfdct which is GPL and also Triangle (also GPL). Manifold looks interesting though ! – DaynaJuliana May 24 '17 at 9:02
10

-UPDATE- Since Postgis3 this function is available by default when SFCGAL is enabled as: ST_ConstrainedDelauneyTriangles


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.

  • Cool, it's good to have this answered with the SFCGAL link and the connection to the other q. I think it occurs relatively regularly so linking to here will be good. – mdsumner May 23 '17 at 10:39
  • Wouldn't this just be using CGAL under the hood? Why not just use CGAL directly? I would love to mark this as answered but because of the restrictive license it doesn't fit my use case :/ Using postgis or shapely is what I am looking for (though I appreciate the thoughtful and correct reply!) – DaynaJuliana May 24 '17 at 0:08
  • @DaynaJuliana: CGAL and SFCGAL are LGPL. Is that too restrictive for you? Postgis itself is GPL. – tilt May 24 '17 at 7:34
  • @tilt are you sure? I thought Postgis is open source? CGAL's trangulation is also GPL (CGAL itself is split between LGPL and GPL). I was looking, and what I've seen online is most people are looking for a python or C implementation. I've tried poly2tri but it crahes on complex polygons – DaynaJuliana May 24 '17 at 8:51
  • I did find this as well, not sure of current state vterrain.org/Implementation/Libs/triangulate.html – DaynaJuliana May 24 '17 at 8:55
0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.