gpclib::tristrip() alternatives

[bastistician]

polyCub.exact.Gauss() is based on polygon triangulation, currently via gpclib::tristrip():

polyCub/R/polyCub.exact.Gauss.R

Lines 91 to 93 in 39ef504

	## triangulation: tristrip() returns a list where each element is a
	## coordinate matrix of vertices of triangles
	triangleSets <- gpclib::tristrip(polyregion)

However, the R package gpclib is unmaintained since 2013 and requires special treatment because of its non-standard license. Sooner or later we might have to switch to an alternative implementation (with support for non-convex polygons).

Simple non-convex example polygon:

xypoly <- list(x = c(3.9, 3.8, 3.7, 3.5, 3.4, 3.5, 3.7, 3.8, 3.8, 3.7,
                     3.7, 3.5, 3.3, 2, 2, 2.7, 2.7, 2.9, 3, 3.3, 3.9),
               y = c(0.7, 1.1, 1.3, 1.7, 1.8, 1.9, 2.1, 2.3, 2.5, 2.8, 3,
                     3.2, 3.3, 3.3, 0.7, 0.7, 1.7, 1.7, 1.5, 0.7, 0.6))
plot(xypoly, type="l")

• gpclib::tristrip() (current approach):

library("gpclib")
xygpc <- new("gpc.poly", pts = list(c(xypoly, list(hole=FALSE))))
xystrip <- tristrip(xygpc)
plot(xygpc)
lapply(xystrip, lines, lty = 2)
## this gives quite a lot of triangles (over each of which we need to integrate)
sum(sapply(xystrip, function(x) nrow(x) - 2))  # 29 triangles

• interp::tri.mesh() (relatively fast Delaunay triangulation):

xymesh <- interp::tri.mesh(xypoly)
plot(xymesh); lines(xypoly, col=2, lwd=2)
nrow(interp::triangles(xymesh))  # 32 triangles
## BUT: still need to find interior triangles

• decido::earcut() (wraps earcut.hpp, which implements a modified ear slicing algorithm):

xyear <- decido::earcut(xypoly)
decido::plot_ears(xypoly, xyear); lines(xypoly, col=2, lwd=2)
dim(xyear) <- c(3, length(xyear)/3)
ncol(xyear)  # 19 triangles

• RTriangle::triangulate() (non-commercial license):

xymat <- as.matrix(as.data.frame(xypoly))
## create line graph object, setting segments explicitly (default would use convex hull)
lg <- RTriangle::pslg(xymat, S = cbind(seq_len(nrow(xymat)), c(seq_len(nrow(xymat))[-1L],1L)))
tri <- RTriangle::triangulate(lg)
plot(tri)
nrow(tri$T)  # also 19 triangles

Runtimes:

timings <- microbenchmark::microbenchmark(
  gpclib = tristrip(xygpc),
  interp = interp::tri.mesh(xypoly),  # BUT: still need to find interior triangles
  decido = matrix(decido::earcut(xypoly), nrow = 3),
  RTriangle = RTriangle::triangulate(RTriangle::pslg(xymat, S = cbind(seq_len(nrow(xymat)), c(seq_len(nrow(xymat))[-1L],1L)))),
  times = 300)
timings

Unit: microseconds
      expr min  lq mean median  uq max neval  cld
    gpclib  64  81   89     90  95 176   300 a   
    interp  95 109  117    117 123 187   300  b  
    decido  17  26   30     31  33  61   300    d
 RTriangle 149 165  173    172 178 282   300   c 

Switch to decido?

• Pro: decido::earcut() ...
  • is almost twice as fast as gpclib::tristrip()
  • produces fewer triangles!
• Con: according to its README, earcut.hpp ...

  doesn't guarantee correctness of triangulation, but attempts to always produce acceptable results for practical data like geographical shapes

Alternative: RTriangle