Constrained triangulation wish

[mdsumner]

I would like to have access to constrained triangulation algorithms. This means that index-pairs on a set of coordinates defining segments (i.e. lines) can be input to constrain the result. This is typically used to modify what would be a Delaunay Triangulation (filling the convex hull) to include those lines. Good algorithms can provide provably good triangular meshes and include extra features like 1) insertions of extra Steiner points (segmentization) 2) maximum triangle area 3) maximum internal angles.

Being able to do this means you can triangulate polygons with concavities and/or holes and only get triangles where you had polygon, which is not generally true of most triangulation algorithms.

Extra coolness would come from being able to update an existing mesh by inserting/removing segments, and being able to apply triangulation or update it within local regions of a data set. CGAL has some of this capacity https://doc.cgal.org/latest/Triangulation_2/ The reason that would be useful, beyond raw performance, is that the triangulation result always depends on the coordinate system in use, so while we often use generic crs like "longitude-latitude" or Mercator for general storage and transfer, it will make sense to project individual features to a local space for triangulating, even if those elements then get used in a totally different coordinate system.

My wishlist questions are:

• how feasible is exposing CGAL in R so we can use its constrained triangulation tools?
• can we use the facilities in BOOST for this, perhaps on the back of @Ironholds work behind the https://github.com/ropensci/wicket project?

I can guide the example data sets and preparing them for input from sf, comparing with rgl and RTriangle and so forth, but compiling external libs into R is very hard for me.

Existing libraries

There are many R packages with triangulation, but only two have a constrained version afaik.

• One is RTriangle, that wraps Shewchuk's Triangle lib: https://www.cs.cmu.edu/~quake/triangle.html it's extremely fast and easy to use and offers many bells and whistles, but is precluded from use in commercial contexts.
• The other is rgl with its extrude3d that applies an ear-clipping algorithm, which is a constrained triangulation but has no guarantees of "goodness". Ear clipping can be used to get small triangles by inserting Steiner points and retriangulating, but it'll never be as good as near-Delaunay.

Does anyone know of any others?

• TDA also links the Triangle lib
• INLA has some mesh/triang facilities, but I think "constrained" means something slightly different in this context
• icosa , webglobe and dggridR are related or at least worth considering
• spatialkernel has a copy of Triangle in it too

CGAL, BOOST, JTS, Manifold GIS, and Triangle have tools for the constrained Delaunay algorithms (I have a longer list of libs somewhere, but it's not longer than 10 or so.

Manifold's Radian engine is putting massive effort into this kind of capability, and are quietly developing a super-powered spatially capable database engine with a free viewer, though Windows-only at this stage (mid 2017).

Triangle is available in R package RTriangle, but under a non-commercial license.

CGAL has been built in small form in an R package here: https://github.com/s-u/rcgal
I've compiled and used triangulation with CGAL very basically via Rcpp, but generally the CGAL code is too advanced for me.

geomery provides convhulln which gives the convex hull in nD, and so closed meshes in 3D are available easily.

Manifold is straightforward but it only offers one form of the algorithm, only works on Windows, and is closed source and proprietary. It's good software, but not a serious contender for this problem at the moment.

The best text I know on this topic is:

https://www.crcpress.com/Delaunay-Mesh-Generation/Cheng-Dey-Shewchuk/p/book/9781584887300

[mdsumner]

Examples

Here's what standard Delaunay will do when applied individually to simple features.

p1 <- cbind(x = c(0, 0, 0.75, 1,   0.5, 0.8, 0.69, 0), 
            y = c(0, 1, 1,    0.8, 0.7, 0.6, 0,    0))
p2 <- cbind(x = c(0.2, 0.2, 0.3, 0.5, 0.5, 0.2), 
            y = c(0.2, 0.4, 0.6, 0.4, 0.2, 0.2))
p4 <- cbind(x = c(0.69, 0.8, 1.1, 1.23, 0.69), 
            y = c(0, 0.6, 0.63, 0.3, 0))
library(sf)
g <- data.frame(a = 1:2)
g[["geometry"]] <- st_sfc(st_multipolygon(list(list(p1, p2[nrow(p2):1, ]))), st_multipolygon(list(list(p4))))
g <- st_as_sf(g)

plot(g)
gt <- g
st_geometry(gt)  <- st_triangulate(st_geometry(g))
plot(gt, col = scales::alpha(c("dodgerblue", "firebrick"), 0.3), main = "convex hull each feature")
plot(st_geometry(g), add = TRUE, lwd = 4)

The sf dataset:

The triangles align to the convex hull of the coordinates of each feature, and overlap each other because they are created without knowledge of the other, or of the input edges generally.

In the constrained version the triangles are completely aligned to the original inputs (triangles in holes are trivially filtered out by seeding the algorithm on input, or detecting them afterwards - and it can be done per-feature or per-feature-set, there are pros and cons for different situations).

I don't have concise code for this on hand but this is the result:

Here's a result from rgl for just the first feature, the actual triangles depend on where the paths start and the algorithm being deterministic or not, but it doesn't need segments as input, and comes with easy 3D plotting and so on.

xyna <- rbind(p1[-1, ], NA, p2[-1, ])
rgl_tr <- rgl::triangulate(xyna, random = FALSE)
plot(g, col = "grey")
apply(rgl_tr, 2, function(a) polygon(xyna[a, ], col = "dodgerblue"))

[edzer]

Nice idea! It would be interesting to get boost and cgal exposed to R, not only for all this kind of things but also to be able to combine and/or compare different implementations.

The good news is that you don't need to compile any external libraries into your package since both cgal and boost:geometry are templates libraries that need header files only - rcgal and wicket show how it works: rcgal has them included in the pkg, wicket gets them from BH.

[Ironholds]

Hmn. Michael, are you thinking a standalone library that (say) accepts sf objects, as part of the whole micro-packages thing, or something built into sf/rcgal?

[mdsumner]

@Ironholds I was thinking standalone, since that gives the most options for the developer. In the first instance it's the developers task to get incoming data into shape for the general tool. If it makes sense to build it into another specific package, in whole or in part then I think that's a natural next step, but not the first step. But really that's all part of what I'm asking for advice on as I am pretty clueless about what is sensible in this space.

The comments by @edzer are very helpful and I'll have a look at how it's done in those packages. I have looked at them, it's just that I couldn't take either to the next level for even a simple test. I tried with rcgal, but resorted to a direct build with CGAL here, which should at least show how basic my understanding is:

https://github.com/r-gris/cgalgris

(I kind of forgot about that CGAL experiment, which came before I discovered RTriangle, and I've been using that to understand this topic ever since - now keen to take it further).

[edzer]

@Ironholds in that context, sfcgal seems relevant to explore, too.

If multiple packages would reuse the CGAL headers, it might make sense to package them only once, and provide them in a CGAL headers-only package similar to BH. (This is really a build-time dependency rather than a run-time dependency, but CRAN can't distinguish between the two). The advantage would be that you'd not have multiple CGAL versions distributed over packages (although date of build would matter).

[Ironholds]

Makes sense. I'm happy to take that on (seems simple enough): want it under the r-spatial branding?

[mdsumner]

Yep!

[Ironholds]

Can do; will set it up on my local, transfer over when done.

[mdsumner]

thanks @Ironholds appreciate you doing this!

[mdsumner]

@Ironholds if you need a point of reference I've got this package for applying the RTriangle code to sf:

https://github.com/r-gris/sfdct

The function ct_triangulate returns a GEOMETRYCOLLECTION type, which is what it must return to store neighbouring POLYGON triangles. Features are self-constrained, there's no concept of constraints imposed by neighbours - and it's not clear really what the behaviour would be in that case.

From the example above this works as expected:

plot(ct_triangulate(g))

It doesn't yet work for GEOMETRYCOLLECTION ifself, so you cannot chain triangulations together with successive parameter changes which is something I want to be able to do, but it should be pretty useable. (Available parameters are in RTriangle::triangulate including control over Steiner point insertion, triangle area, internal angle, enforcing Delaunay criterion to the result, and a few others).

In the past I've had similar tools in gris and rangl for sp objects, but without a closed-loop of types as here.

[Ironholds]

Noted! FWIW I was more imagining just packaging the headers up and leaving them be (at least, those headers that are self-contained or depend only on libraries with existing Rcpp header-only libraries). The absence of doing so thus far is more down to current professional time commitments and an ongoing illness than anything else.

[mdsumner]

All good, another heads-up I've found sphere-aware triangulation is in INLA, but not sure yet how it contrasts to other options.

I had also forgotten that we can do 3D-geometry convex hulls with geometry (well nD) and so it's easy to build closed (but not constrained) Delaunay meshes with that. It's a neat example of why this isn't a simple topic to summarize. The INLA option makes it timely to tell the story of what is available across R, there's at least 12 packages to talk about

[Robinlovelace]

Any more recent updates on this @Ironholds ? Just stumbled the http://sfcgal.org/ website that @edzer mentioned above and it looks promising and I see it's already implemented in PostGIS seemingly.

[Ironholds]

Not at the moment :/. I would unfortunately not expect it out of me any time soon - family events, geopolitics and graduate school are sucking up all my free time. If anyone else has the time and will, feel free to go at it!

[mdsumner]

I think I can do it now, I understand the configuration stuff enough and rwinlib works well. I did some very basics at r-gris/cgalgris, so it's not much but enough to start hacking.

[Robinlovelace]

Exciting stuff, please keep us posted!

[mdsumner]

If have you specific things you need that can help experimentation with CGAL and I'd love to start with small-ish real examples when I do try again. The C++ style in CGAL is way more sophisticated than the usual GDAL stuff, it's very head-stretching for me. My key interest will be in performing local triangulations and updating them dynamically, and in having engines for high-fidelity and flexible "overlay" operations (fusing contour layers with polygonal regions and with road networks for example).

RTriangle can already do this of course, and I have tested what I need mostly with rangl. Now mdsumner/sc is a better overall approach with an eye to a future general model (they are just both a mess as they are learning steps for me). I'm exploring ways to represent segments/edges and flip between those and path models and I see this now as a crux foundation that comes up all over the place. The work by @mpadge in osmdata and friends is very much related, and the need for multi-table things in ggplot al a tidygraph and the current limited raster support are clearly related.

(I personally don't have a problem with RTriangle's non-commercial license, it seems perverse to me that it's considered "non open" just because corporations are precluded, it muddies the water for basic research in publically funded institutions. I have no problem with folks choosing not to use non-CC because of their chosen funding arrangement, but there's no reason it should limit research that is otherwise free in how it can be carried out. )

[mdsumner]

This is now progressed well in r-gris/cgalgris, I have all I need to move forward there.

There's a functioning simple features-specific facility in the sfdct CRAN package.

[Robinlovelace]

Great stuff - will aim to check these out when I get a chance - may be a while though.