---
title: "decido - polygon triangulation by ear clipping"
author: "Michael D. Sumner"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 7
vignette: >
  %\VignetteIndexEntry{decido}
  %\VignetteEngine{knitr::knitr}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Constrained triangulation of polygons

Decido aims to provide very fast constrained triangulations for polygons using a robust C++ library created by Mapbox. 

Constrained triangulation by ear clipping is an important low-cost method for decomposing a polygon into triangles, used for 3D visualization and many analytical techniques. Here *constrained* means **shape-preserving**, in the sense that every edge of the polygon will be included ^[1] as the edge of a triangle in the result. 

Ear clipping (or cutting) must be applied to a *path-based* polygon that consists of only one island, with zero or more holes. Any multiple island polygons must be triangulated separately by this method. 

The better-known *Delaunay triangulation* is generally not shape-preserving as it works only on points, without knowledge of input line segments or polygon edges. Shape-preserving Delaunay or more commonly *near-Delaunay* triangulation is performed on a set of edges rather than paths, and can be
performed on multiple polygons at once. Holes require careful trimming of triangles after the decomposition is performed in this approach.  (See the [anglr package](https://CRAN.r-project.org/package=anglr) for the `DEL0()` function for high-quality triangles, with control over size, angle, and the Delaunay criterion). 


## Example

This is a basic example of triangulating a single-ring polygon. The
output is a vector of triplet indices defining each triangle. 


```{r example}
library(decido)
x <- c(0, 0, 0.75, 1, 0.5, 0.8, 0.69)
y <- c(0, 1, 1, 0.8, 0.7, 0.6, 0)
(ind <- earcut(cbind(x, y)))


plot_ears(cbind(x, y), ind)
```

Support for holes is provided by the argument `holes`. The values are the starting index of each hole, here in R's 1-based convention. 

```{r}
## polygon with a hole
x <- c(0, 0, 0.75, 1, 0.5, 0.8, 0.69,
     0.2, 0.5, 0.5, 0.3, 0.2)
y <- c(0, 1, 1, 0.8, 0.7, 0.6, 0,
     0.2, 0.2, 0.4, 0.6, 0.4)
ind <- earcut(cbind(x, y), holes = 8)
plot_ears(cbind(x, y), ind)
```


The hole-specification is a little subtle, since usually R's functions
(polygon and polypath, and others) expect NA values to separate paths. 


Notice how the hole begins at index 8, hence `holes = 8` above, and `holes = c(8, 13)` below. 

```{r}
plot_ears(cbind(x, y), ind, col = "grey", border = NA)
text(x, y, labels = seq_along(x), pos = 2)

```

The method is subtle. 

This example adds a third polygon, a second hole in the island. 

```{r}
## add another hole
x <- c(0, 0, 0.75, 1, 0.5, 0.8, 0.69,
     0.2, 0.5, 0.5, 0.3, 0.2,
      0.15, 0.23, 0.2)
y <- c(0, 1, 1, 0.8, 0.7, 0.6, 0,
     0.2, 0.2, 0.4, 0.6, 0.4,
      0.65, 0.65, 0.81)
ind <- earcut(cbind(x, y), holes = c(8, 13))
plot_ears(cbind(x, y), ind, col = "grey")

```


## Performance for developers

There is a headers API for decido. 

```{r headers, eval = FALSE, include = TRUE}
## this code is not run in the vignette
## but can be run as is with decido installed
library(Rcpp)

cppFunction(
  depends = "decido"
  , includes = '#include "decido/decido.hpp"'
  , code = '
    Rcpp::IntegerVector earcut0( SEXP polygon ) {
      return decido::api::earcut( polygon );
    }
  '
)

poly <- list(matrix(c(0,0,0,1,1,1,1,0,0,0), ncol = 2, byrow = T))
earcut0( poly )
# [1] 1 4 3 3 2 1
```


## Orientation

Triangles are returned in *counter-clockwise* order, as is proper and good in a world where area is positive (I'm kidding, it's just a convention that folks use ... see `fortunes::fortune("illogical")`). 

First a function so I can show this with a plot. 

```{r triplot}
## plot triangles (input is triplets of xy coordinates)
## with one side an oriented arrow
plot_tri <- function(x, add = TRUE) {
  if (!add) plot(x, asp = 1, type = "n")
  idx <- c(1:3, 1)
  for (i in seq_len(nrow(x)/3)) {
    #print(idx)
    ## plot only the first arrow
    arrows(x[idx[1], 1], x[idx[1], 2], 
           x[idx[2], 1], x[idx[2], 2], length = 0.25, lwd = 2, col = "firebrick")
    ## segments the rest
    segments(x[idx[2:3], 1], x[idx[2:3], 2], 
           x[idx[3:4], 1], x[idx[3:4], 2])
  
    idx <- idx + 3
    
  }
}

```


See how each triangle has its first edge (the first two indexes) shown as an arrow. They're all counter-clockwise. I've checked this on very large sets of real-world inputs, you can get 0-area polygons and *very small negative-epsilon area* polygons when the 3 points are collinear or so close to it that numeric precision fails us. 

```{r orientation}
plot_tri(cbind(x, y)[ind, ], add = FALSE)
```

(See [this discussion](https://github.com/hypertidy/decido/issues/17) for more on the topic where I was exactly backwards on the counter-clockwise thing, but airing the topic was very helpful). 

Orientation is not really important in planar contexts, though you'll find very little consensus on this. (Would it not be handy if all holes in polygons had *negative area* ...).  It really matters in 3D because the orientation has a more definite meaning when the triangles are representing surfaces. 

## Comparing paths and triangles

This example defines a much simpler shape, the minimal shape able to
be decomposed to triangles (and not stay a triangle). 

A quadrilateral, with two holes that are open to each other allows the use of the same data, and we can tweak whether we wanted one hole or two. This is an important example used for validating early versions of this package. 

```{r}
x <- c(0, 0, 1, 1,
       0.4, 0.2, 0.2, 0.4,
       0.6, 0.8, 0.8, 0.6
)
y <- c(0, 1, 1, 0,
       0.2, 0.2, 0.4, 0.4,
       0.6, 0.6, 0.4, 0.4
)
ind <- decido::earcut(cbind(x, y), holes = c(5, 9))
plot_ears(cbind(x, y), ind, col = "grey")
title("triangle plot, two holes")
plot_holes(cbind(x, y), holes = c(5, 9), col = "grey")
title("path plot, two holes")

ind <- decido::earcut(cbind(x, y), holes = 5)
plot_ears(cbind(x, y), ind, col = "grey")
title("triangle plot, two holes as one")
plot_holes(cbind(x, y), holes = 5, col = "grey")
title("path plot, two holes as one")
```

## A geographic example

For good measure we include a geographic example, a triangulation
of the mainland part of Tasmania from the `oz` package. 


```{r}
library(oz)
oz_ring <- oz::ozRegion(states = FALSE)
ring <- oz_ring$lines[[6]]
indices <- earcut(ring[c("x", "y")])
plot_ears(cbind(ring$x, ring$y), indices)
```


## Variants

The actual triangulation obtained depends on where the polygons start. This is complicated, because of the sheer number of possible variants, combinations of starting points among the island and its holes. 

Important "edge" cases are degeneracies, holes touching the island or each other, duplicated edges, intersecting edges, zero-length edges, holes actually external to islands, already existing triangles, and existing quadrilaterals (amongst others). We aren't going to explore those here, and we are reasonably confident that Mapbox has been presented with a rich enough pool of variant polygons to make its library pretty robust. There's no uniquely "correct" here either, different systems and standards will apply different rules and allow or choose differently. 


First a function to "rotate" coordinates to different start/end. 

```{r}
vecrot <- function(x, k) {
  if (k < 0 || k > length(x)) warning("k out of bounds of 'x' index")
  k <- k %% length(x)
  ## from DescTools::VecRot
  rep(x, times = 2)[(length(x) - k + 1):(2 * length(x) - k)]
}
```

Now plot each possible variant of defining the polygon ring by traversing the boundary anti-clockwise from 
a different starting vertex, shown as a filled point symbol. 

```{r}
x <- c(0, 0, 0.75, 1, 0.5, 0.8, 0.69)
y <- c(0, 1, 1, 0.8, 0.7, 0.6, 0)
(ind <- earcut(cbind(x, y)))

plot_ears(cbind(x, y), ind)
points(x[1], y[1], pch = 19, col = "firebrick")
title("original")
op <- par(mfrow = c(2, 3))
for (rot in head(seq_along(x), -1)) {
xx <- vecrot(x, rot); yy <- vecrot(y, rot)
ind <- earcut(cbind(xx, yy))
plot_ears(cbind(xx, yy), ind)
title(sprintf("rot %i", rot))
points(xx[1], yy[1], pch = 19, col = "firebrick")
}
par(op)


```


Compare this to constrained Delaunay triangulation from the `RTriangle` package (not illustrated here 
but the result is that the *constrained* triangulation is different from all variants above, and the *conforming* triangulation inserts two points not included in the original data - these are called Steiner points). 

```R
library(sfdct)
library(sf)
xsf <- st_sfc(st_polygon(list(cbind(x, y)[c(seq_along(x), 1), ])))
plot(ct_triangulate(xsf, Y = TRUE)[1], col = NA)

## now a conforming Delaunay triangulation
plot(ct_triangulate(xsf, D = TRUE)[1], col = NA)

```

## Performance

Compare timing of C++ versus JS implementations, the original C++ decido was born from the Javascript version 'rearcut'. With version 0.3.0 decido gained a C++ headers API thanks to David Cooley. 


```R

ring <- 
remotes::install_github("hypertidy/rearcut")

library(Rcpp)
## some people love to see the world burn by
## ignoring indexes and working purely in coordinates ...
cppFunction(
  depends = "decido"
  , includes = '#include "decido/decido.hpp"'
  , code = '
    Rcpp::IntegerVector earcut0( SEXP polygon ) {
      return decido::api::earcut( polygon );
    }
  '
)

dim(ring <- decido:::ant_cont )

rbenchmark::benchmark(js = rearcut::earcut(ring), 
                      decido = decido::earcut(ring),
                      headers_decido = earcut0(list(ring)))
                      
                      
#            test replications elapsed relative user.self sys.self user.child sys.child
# 2         decido          100   0.304    1.924     0.305    0.000          0         0
# 3 headers_decido          100   0.158    1.000     0.158    0.000          0         0
# 1             js          100   6.876   43.519     6.818    0.148          0         0
```


[1]: See [https://hypertidy.github.io/sfdct/reference/ct_triangulate.html#examples] of `st_triangulate` vs `ct_triangulate` on a set of polygons, the unconstrained version does not include all the polygon edges, the boundaries are the convex of *each polygon* and so each triangle collection is overlapping. 


```{r, include=FALSE}
# , and it's also not the only way to do it, the grid package uses a grouping vector rather than a sparse index like this. The spatial packages sp and sf explicitly use structural hierarchies rather than more abstract specifications. The fortify-approach in ggplot2 is more like the grid one. A sparse representation is closer to what is needed for topological operations and visualization, consider that when we have triangles there are no need for "holes", we can identify which triangles will be plotted and how (or not), and an index into the vertices available becomes a key efficiency feature. (See [silicate](https://github.com/hypertidy/silicate) for a lot more on this topic). 
```