Constrained Delaunay Triangulation with Refinement

Binding to the 'CDT' header-only C++ library: constrained and conforming Delaunay triangulation, Ruppert refinement by area or minimum angle (Ruppert, 1995, ), resolution of crossing constraints, and topological region depth per triangle. Aimed at replacing 'RTriangle' in 'anglr'.

cdtr

Constrained Delaunay triangulation with refinement for R, binding the header-only CDT library (MPL-2.0). It takes the same input as a planar straight line graph (vertices, segments as index pairs, optional per-vertex attributes) and returns vertices, triangles, and a per-triangle depth that says where each triangle sits relative to the constraints, with no point-in-polygon step.

remotes::install_github("hypertidy/cdtr")

The pipeline, on a toy

Vertices and segments in, vertices and triangles out. Here a square with a square hole: two closed rings of segments.

library(cdtr)
x <- c(0, 1, 1, 0, 0.25, 0.75, 0.75, 0.25)
y <- c(0, 0, 1, 1, 0.25, 0.25, 0.75, 0.75)
s0 <- 1:8; s1 <- c(2, 3, 4, 1, 6, 7, 8, 5)
r <- cdt_triangulate(x, y, s0, s1)
str(r[c("P", "T", "S", "depth")])
#> List of 4
#>  $ P    : num [1:8, 1:2] 0 1 1 0 0.25 0.75 0.75 0.25 0 0 ...
#>  $ T    : int [1:10, 1:3] 8 6 8 7 6 5 8 6 2 2 ...
#>  $ S    : int [1:8, 1:2] 1 3 1 5 6 2 7 5 2 4 ...
#>  $ depth: int [1:10] 2 1 1 1 2 1 1 1 1 1

depth is the number of constraint boundaries crossed on the way in from outside: 0 is outside everything, and for nested rings odd is inside, even is a hole. By default triangles outside the outermost boundary are already dropped (erase = "outer", which is what Triangle’s -p does); erase = "holes" also drops the even depths.

plot_mesh <- function(m, main = "") {
  e <- cdt_edges(m)
  plot(m$P, asp = 1, axes = FALSE, xlab = "", ylab = "", pch = ".", main = main)
  segments(m$P[e[, 1], 1], m$P[e[, 1], 2], m$P[e[, 2], 1], m$P[e[, 2], 2], col = "grey50")
  segments(m$P[m$S[, 1], 1], m$P[m$S[, 1], 2], m$P[m$S[, 2], 1], m$P[m$S[, 2], 2], lwd = 2)
}
op <- par(mfrow = c(1, 3), mar = c(0.5, 0.5, 2, 0.5))
plot_mesh(r, "constrained")
plot_mesh(cdt_triangulate(x, y, s0, s1, max_area = 0.01), "max_area = 0.01")
plot_mesh(cdt_triangulate(x, y, s0, s1, max_area = 0.01, min_angle = 25, erase = "holes"),
          "+ min_angle = 25, holes erased")

par(op)

Refinement reports what it could not do

Ruppert refinement cannot always meet the bound (two constraints meeting at a sharp angle, say). Rather than silently stopping, the result carries the counts.

r2 <- cdt_triangulate(x, y, s0, s1, max_area = 0.01, min_angle = 25)
r2$unrefined
#>   criterion shortEdgeTriangles circumcenterOutside circumcenterOnVertex
#> 1      area                  0                   0                    0
#> 2     angle                  0                   0                    0
#>   sharpFixedCorner shortEdges splitVertexInvalid
#> 1                0          0                  0
#> 2                0          0                  0
r2$min_edge_length
#> [1] 0.025

min_edge_length is the floor below which the refiner gives up on an edge; by default it is derived from the input (a fraction of the median segment length, capped by the area target) so refinement does not cascade into sharp input corners. Pass 0 for never-give-up.

Real polygons

Any polygon source reduces to the same two tables. With sf’s North Carolina counties, rings become vertex runs and consecutive index pairs; shared boundaries dedupe to one segment.

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
rings <- unlist(lapply(sf::st_geometry(nc), function(g) lapply(unclass(g), `[[`, 1)), recursive = FALSE)
xy <- do.call(rbind, lapply(rings, function(m) m[-nrow(m), ]))
n <- vapply(rings, function(m) nrow(m) - 1L, 1L)
start <- rep(cumsum(c(0L, n[-length(n)])), n)
s0 <- seq_along(xy[, 1])
s1 <- start + (s0 - start) %% rep(n, n) + 1L
m <- cdt_triangulate(xy[, 1], xy[, 2], s0, s1, max_area = 0.01)
table(m$depth)
#> 
#>    1    3    5 
#> 2100 1119  521
op <- par(mar = rep(0, 4))
plot_mesh(m)

par(op)

The counties tile the state, so there are no holes, but a county two boundaries in from the coast has depth 3, and so on: depth is the peel layer, a count of constraint crossings from outside, not an inside/outside flag in general. It is computed topologically from the mesh alone, before anything is erased, which makes it a cheap first cut for classifying triangles by region: nested rings get odd/even, a coverage gets a boundary-distance ordering, and only ties within a layer need a geometric lookup.

Attributes ride along

Per-vertex values (an elevation, a time) are carried onto the new Steiner vertices by linear interpolation from the input mesh, the way Triangle’s PA works.

z <- 2 * xy[, 1] + 3 * xy[, 2]
ma <- cdt_triangulate_attr(xy[, 1], xy[, 2], s0, s1, PA = cbind(z = z), max_area = 0.01)
range(ma$PA[, "z"] - (2 * ma$P[, 1] + 3 * ma$P[, 2]))
#> [1] -1.421085e-14  1.421085e-14

What this is for

cdtr is the triangulation primitive for the hypertidy mesh packages, where it replaces RTriangle (Shewchuk’s Triangle, non-commercial licence). The inst/benchmarks/ scripts compare the two on anglr’s Tasmanian contour and cadastre data; the short version is that constrained-only output is identical, refinement honours max_area where Triangle may give up, and Triangle’s refiner is faster per Steiner point.