checking triangle orientation

[mdsumner]

With motivation from @brodieG

This example still WIP, but has

• simplest polygon (of two triangles) earcut from a clockwise and anti-clockwise input
• the triangles are anti-clockwise (as is sensible, but different from readme claim)
• plot helpers to provide index of input poly, to compare with earcut(poly)
• plot helper to add arrows of first triangle edge, to aid visual confirm of orientation (and label for area at centre, properly negative)

I'm keen to see any counter examples, does earcut.hpp ever return clockwise-oriented triangles? I don't think so, there are collinear (or visually-collinear at any rate) cases with zero-area, or very small epsilon-positive area - I see that triangulating silicate::inlandwaters, sthng like 44 out of 33000 triangles there.

  ## vectorized index helpers for calculating triangle area
tri_ix <- function(x) {
  offset <- rep(seq(0, nrow(x) - 1, by = 3), each = 3)
  c(3L, 1L, 2L) + offset
}

tri_jx <- function(x) {
  offset <- rep(seq(0, nrow(x) - 1, by = 3), each = 3)
  c(1L, 2L, 3L) + offset
}

## signed triangle area (negative is anti-clockwise orientation)
### https://github.com/hypertidy/silicate/blob/ab7610290e45567c134dd7e48f4ba8199c5bcf59/R/triangles.R
tri_wind <- function(x) {
  ix <- tri_ix(x)
  jx <- tri_jx(x)
  colSums(matrix((x[ix, 1] + x[jx, 1]) * (x[ix, 2] - x[jx, 2]), nrow = 3L))/2
}


library(decido)
## plot polygon as index labels at its vertices
plot_poly <- function(x) {
  xr <- range(x[,1])
  yr <- range(x[,2])
  xr <- xr + c(-1, 1) * diff(xr)/10
  yr <- yr + c(-1, 1) * diff(yr)/10
  
  plot(x, type = "n", xlim = xr, ylim = yr)
  text(x, lab = 1:nrow(x), pos = 4)
}
## plot triangles (input is triplets of xy coordinates)
## with one side an oriented arrow
## signed area of triangle is labelled at its centre
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.35, lty = 2, lwd = 2)
    ## segments the rest
    segments(x[idx[2:3], 1], x[idx[2:3], 2], 
           x[idx[3:4], 1], x[idx[3:4], 2])
    text(mean(x[idx[1:3], 1]), 
         mean(x[idx[1:3], 2]), lab = tri_wind(x[idx[1:3], ]), 
         col = "firebrick", cex = 0.7)
    idx <- idx + 3
    
  }
}


wise <- cbind(c(0, 0, 1, 1), 
              c(0, 1, 1, 0))
awise <- wise[nrow(wise):1, ]

plot_poly(wise)
title("clockwise poly")
plot_tri(wise[earcut(wise), ])

plot_poly(awise)
title("anticlockwise poly")
plot_tri(awise[earcut(awise), ])

## the signed-area should always be negative (sometimes zero or 
## epsilon-postive with a collinear triangle)
tri_wind(awise[earcut(awise), ])
#> [1] -0.5 -0.5

tri_wind(wise[earcut(wise), ])
#> [1] -0.5 -0.5



## 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_poly(cbind(x, y))
plot_tri(cbind(x, y)[ind, ])


sum(abs(tri_wind(cbind(x, y)[ind, ])))
[1] 0.672
sf::st_area(silicate::minimal_mesh[1, ])
[1] 0.672

^{Created on 2020-04-27 by the reprex package (v0.3.0)}

[mdsumner]

Compare rgl, there's a few gotchas for getting into rgl form:

• polygons with holes must be in polygon() NA-split form (this makes indexing a bit confusing as they are relative to the input-NA coords)
• no multi polygons
• winding of inputdoesn't matter
• triangles are anti-clockwise
• turn of the random part of rgl::triangulate or polygon3d

## vectorized index helpers for calculating triangle area
tri_ix <- function(x) {
  offset <- rep(seq(0, nrow(x) - 1, by = 3), each = 3)
  c(3L, 1L, 2L) + offset
}

tri_jx <- function(x) {
  offset <- rep(seq(0, nrow(x) - 1, by = 3), each = 3)
  c(1L, 2L, 3L) + offset
}

## signed triangle area (negative is anti-clockwise orientation)
### https://github.com/hypertidy/silicate/blob/ab7610290e45567c134dd7e48f4ba8199c5bcf59/R/triangles.R
tri_wind <- function(x) {
  ix <- tri_ix(x)
  jx <- tri_jx(x)
  colSums(matrix((x[ix, 1] + x[jx, 1]) * (x[ix, 2] - x[jx, 2]), nrow = 3L))/2
}



## plot polygon as index labels at its vertices
plot_poly <- function(x) {
  xr <- range(x[,1])
  yr <- range(x[,2])
  xr <- xr + c(-1, 1) * diff(xr)/10
  yr <- yr + c(-1, 1) * diff(yr)/10
  
  plot(x, type = "n", xlim = xr, ylim = yr)
  text(x, lab = 1:nrow(x), pos = 4)
}
## plot triangles (input is triplets of xy coordinates)
## with one side an oriented arrow
## signed area of triangle is labelled at its centre
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.35, lty = 2, lwd = 2)
    ## segments the rest
    segments(x[idx[2:3], 1], x[idx[2:3], 2], 
             x[idx[3:4], 1], x[idx[3:4], 2])
    text(mean(x[idx[1:3], 1]), 
         mean(x[idx[1:3], 2]), lab = tri_wind(x[idx[1:3], ]), 
         col = "firebrick", cex = 0.7)
    idx <- idx + 3
    
  }
}


## here we need a object-to-polycoords-na thing
library(silicate)
#> 
#> Attaching package: 'silicate'
#> The following object is masked from 'package:stats':
#> 
#>     filter
x <- minimal_mesh
p <- sc_path(x) %>% dplyr::filter(object_ == object_[1])
xy <- sc_coord(x)[1:sum(p$ncoords_), ]
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
rl <- p %>% dplyr::filter(object_ == object_[1]) %>% 
  pull(ncoords_)
rbind_na <- function(x) {
  head(do.call(rbind, lapply(x, function(a) rbind(head(a, -1), NA))), 
       -1)
}
l <- split(xy, rep(seq_along(rl), rl))
## reverse winding of hole
#l[[2]] <- l[[2]][nrow(l[[2]]):1, ]
pcd <- rbind_na(l)

library(rgl)
## check rgls' take (make sure random = FALSE!)
tri <- rgl::triangulate(pcd$x_, pcd$y_, random = FALSE)
tri_wind(as.matrix(pcd)[tri, ])
#>  [1] -0.0690 -0.1000 -0.0200 -0.0300 -0.0190 -0.0285 -0.0300 -0.0550 -0.0500
#> [10] -0.1125 -0.0625 -0.0955
plot_poly(pcd[!is.na(pcd$x_), ])
plot_tri(as.matrix(pcd)[tri, ])

^{Created on 2020-04-29 by the reprex package (v0.3.0)}

[mdsumner]

You can see the effect of triangle winding in rgl by setting back = "lines"

## rgl polygon with a hole
xy <- structure(c(0, 0, 0.75, 1, 0.5, 0.8, 0.69, NA, 0.2, 0.5, 0.5, 
            0.3, 0.2, 0, 1, 1, 0.8, 0.7, 0.6, 0, NA, 0.2, 0.2, 0.4, 0.6, 
            0.4), .Dim = c(13L, 2L), .Dimnames = list(NULL, c("x_", "y_")))
tri <- 
structure(c(7L, 9L, 1L, 9L, 2L, 1L, 9L, 13L, 2L, 7L, 10L, 9L, 
            7L, 11L, 10L, 7L, 5L, 11L, 5L, 12L, 11L, 5L, 2L, 12L, 2L, 13L, 
            12L, 5L, 3L, 2L, 5L, 4L, 3L, 7L, 6L, 5L), .Dim = c(3L, 12L), nextvert = c(7L, 
                                                                                      1L, 2L, 3L, 4L, 5L, 6L, NA, 13L, 9L, 10L, 11L, 12L))


rgl::clear3d()
## the "back" is triangle-orientation
triangles3d(cbind(xy, 0)[tri, ], back = "lines")

## make some triangles clockwise
tri[, c(1, 3, 5)] <- tri[3:1, c(1, 3, 5)]
rgl::clear3d()
triangles3d(cbind(xy, 0)[tri, ], col = "firebrick", back = "lines")

The triangles are otherwise filled on the other side

[brodieG]

Nice. I'm going to assume going forward that they will always come out a-clockwise out of decido. Will you consider guaranteeing this in the docs (and perhaps with a check for when upgrading earcut.hpp)?

Thanks again!

[mdsumner]

Sure, that should be fine. Are you confident that your clockwise-winded triangle with the vertical side was a false-positive? I'd still like to have that R-polygon, that's a neat example - I couldn't find an easy version of the logo polygon anywhere.

It's been very nice to clean out that wind-order, area, and plot code - that's going to be very handy in lots of places and this has flushed out some confusion for me!

[brodieG]

Yes, pretty sure it was a false positive caused by my assuming that sign() returns only c(1, -1) when it returns 0 for 0, which is me just being an idiot and forgetting and not re-reading docs.

Here is the code to get the R logo in polygons (and the hush hush bit is I'm working on getting it rayrendered for a post of mine, which is how I ended up down this particular rabbit hole).

logo.page <- readLines('https://www.r-project.org/logo/Rlogo.svg')
dat <-regmatches(logo.page, gregexpr("\"[^\"]+\"", logo.page))

get_path <- function(dat) {
  path <- unlist(strsplit(sub('^"(.*)"$', '\\1', dat), " "))
  path <- path[-length(path)]  # drop "Z" command at end
  cmd <- sub("^([A-Z]*).*$", "\\1", path)
  x <- sub("[^0-9]*([0-9.]*),.*$", "\\1", path)
  y <- sub(".*,([0-9.]*).*$", "\\1", path)
  d <- data.frame(i=seq_along(x), x=as.numeric(x), y=as.numeric(y), cmd=cmd)
  d
}
library(gridBezier)

width <- as.numeric(gsub('"', '', dat[[1]][4]))
height <- as.numeric(gsub('"', '', dat[[1]][5]))

dr <- get_path(dat[[13]][[1]])  # R logo
dh <- get_path(dat[[12]][[1]])  # Hula

dr[[2]] <- dr[[2]] / width
dr[[3]] <- dr[[3]] / height
dh[[2]] <- dh[[2]] / width
dh[[3]] <- dh[[3]] / height

# - R Path ---------------------------------------------------------------------

brle <- with(dr, rle(cmd == "C" | cmd == ""))
bends <- with(brle, cumsum(lengths)[values])
bstarts <- with(brle, cumsum(lengths)[!values])

bzs <- Map(
  function(start, end) {
    points <- dr[start:end, c('x', 'y')]
    BezierGrob(points[[1]], points[[2]], stepFn=nSteps(20))
  },
  bstarts,
  bends
)
# BezierPoints returns in inches...

bzsp <- lapply(bzs, BezierPoints)
bzspi <- lapply(bzsp, lapply, unit, 'inches')
bzsp <- lapply(
  bzspi, function(x) list(convertX(x[[1]], 'npc'),convertY(x[[2]], 'npc'))
)
# plot.new(); lapply(bzs, grid.draw)
# reconnect with the line segments

lrle <- with(dr, rle(cmd == "L"))
lends <- with(lrle, cumsum(lengths)[values])
lstarts <- with(lrle, cumsum(lengths)[!values][seq_along(lends)])

lnsp <- Map(
  function(start, end) dr[start:end, c('x', 'y')],
  lstarts,
  lends
)
lns <- lapply(
  lnsp,
  function(x) grid.lines(x[[1]], x[[2]], gp=gpar(col='#000000'))
)
# plot.new(); lapply(lns, grid.draw)

# plot.new()
# for(points in bzsp) {
#   x <- points[[1]]
#   y <- points[[2]]
#   grid.draw(grid.lines(x, y, gp=gpar(col='#000000')))
# }

# reconnect with the line segments
out <- c(bstarts < 33, lstarts < 33)
outside <- c(bzsp, lnsp)[out][order(c(bstarts, lstarts)[out])]
inside <- c(bzsp, lnsp)[!out][order(c(bstarts, lstarts)[!out])]

# plot.new()
# for(points in list(outside)) {
#   x <- unlist(lapply(points, '[[', 1))
#   y <- unlist(lapply(points, '[[', 2))
#   grid.draw(grid.lines(x, y, gp=gpar(col='#000000')))
# }

outx <- unlist(lapply(outside, '[[', 1))
outy <- unlist(lapply(outside, '[[', 2))

inx <- unlist(lapply(inside, '[[', 1))
iny <- unlist(lapply(inside, '[[', 2))

allx <- c(outx, inx)
ally <- -c(outy, iny) + 1

# - Hula -----------------------------------------------------------------------

# This one is all beziers

hbrle <- with(dh, rle(cmd == "C" | cmd == ""))
hbends <- with(hbrle, cumsum(lengths)[values])
hbstarts <- with(hbrle, cumsum(lengths)[!values])

hbzs <- Map(
  function(start, end) {
    points <- dh[start:end, c('x', 'y')]
    BezierGrob(points[[1]], points[[2]], stepFn=nSteps(20))
  },
  hbstarts,
  hbends
)
hbzsp <- lapply(hbzs, BezierPoints)
hbzspi <- lapply(hbzsp, lapply, unit, 'inches')
hbzsp <- lapply(
  hbzspi, function(x) list(convertX(x[[1]], 'npc'),convertY(x[[2]], 'npc'))
)
# plot.new()
# grid.draw(grid.lines(hbzsp[[1]][[1]], hbzsp[[1]][[2]]))
# grid.draw(grid.lines(hbzsp[[2]][[1]], hbzsp[[2]][[2]]))

hulax <- c(hbzsp[[1]][[1]], hbzsp[[2]][[1]])
hulay <- -c(hbzsp[[1]][[2]], hbzsp[[2]][[2]]) + 1

# hind <- decido::earcut(cbind(hulax, hulay), holes=length(hbzsp[[1]][[1]]) + 1)
# decido::plot_ears(cbind(hulax, hulay), hind, col = "#1e64b6")

# - Plotting -------------------------------------------------------------------

# ind <- decido::earcut(cbind(allx, ally), holes=length(outx) + 1)
# decido::plot_ears(cbind(allx, ally), ind, col = "#1e64b6")

# ind <- decido::earcut(cbind(outx, 1-outy))
# decido::plot_ears(cbind(outx, 1-outy), ind, col = "#1e64b6")

[mdsumner]

haha nice manual svg parsing - glad you have done that I baulked!

[brodieG]

Whatever it takes =)

[mdsumner]

I think I got the orientation thing exactly backwards. Signed area should be positive counter-clockwise.

[brodieG]

That's right. Doesn't really matter though, it's pretty obvious once you know about signed areas, so thanks a bunch for showing me that. Kinda like the earcut.hpp docs have the direction backwards. It's important, but really the truly important thing is that it's consistent. You can always spot check and figure out the sign/direction (though, it wouldn't be a bad thing for the docs to be correct).

I don't know how long it would have taken me to run into signed areas on my own given it's not at all my area of expertise.

[mdsumner]

Cool appreciate the feedback, I always need this extra validation to know I'm not off-base, I will check with some local math experts on the general idea and put this in a tiny package.

[mdsumner]

Note to self

• Check why earcut readme claims clockwise
• Update source earcut.hpp and check zero/epsilon-positive cases mentioned above

[mdsumner]

@brodieG I haven't gone any deeper, but it's curious that this example page uses clockwise = positive?

https://rosettacode.org/wiki/Shoelace_formula_for_polygonal_area

It's entirely arbitrary of course but ...

[brodieG]

From my reading, I'm seeing that the absolute value of that equation is positive for clockwise. Due to the use of abs it's hard to know anything about the relationship between CW and CCW and sign. Wolfram is more explicit, and that lines up with what I found when I implemented it for use in rayrender.

[mdsumner]

Oh I see, they're being a bit casual and really calculating "area". Thanks for replying, that link is nice too.

[mdsumner]

Thanks @brodieG I put some of this in the doc, thanks! Needs a much more thorough treatment as part of a bigger story that's developing in R ;)