A naive implementation of Visvalingam's line simplification in R

2022-05-07

Visvalingham’s Line Simplification Algorithm

  • Given a number of points defining a line
  • Find the areas of the triangles defined by each set of 3 continuous points
  • Remove the point which is the centre of the triangle with the minimum area
  • Recalculate areas (as two will have changed when removing the point)

Naive R implementation

  • Full set of all triangle areas is calculated prior to the remove of each point
    • A priority queue would help here

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#' Visvalingham Line Simplification
#'
#' @param x,y coordinates
#' @param n number of points to keep
#'
#' @importFrom utils head tail
#' @export
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
visvalingam <- function(x, y, n) {

  # Sanity Check
  stopifnot(length(x) == length(y))
  stopifnot (n <= length(x) && n >= 2)
  if (length(x) == 2) {
    return(list(x=x, y=y))
  }

  # Remove points
  for (i in seq_len(length(x) - n)) {
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Find areas
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    x1 <- head(x, -2)
    x2 <- head(x[-1], -1)
    x3 <- tail(x, -2)

    y1 <- head(y, -2)
    y2 <- head(y[-1], -1)
    y3 <- tail(y, -2)

    a0 <- x1 - x2
    a1 <- x3 - x2

    b0 <- y1 - y2
    b1 <- y3 - y2

    tri_areas <- abs(a0 * b1 - a1 * b0)  # / 2

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Find minimim area triangle and remove its centre vertex from all pts
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    min_tri <- which.min(tri_areas)
    
    x <- x[-(min_tri + 1)]
    y <- y[-(min_tri + 1)]
  }

  list(x = x, y = y)
}

Demo

This demo shows how visvalingam can be used to tidy up a noisy drawing of a square

library(grid)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Create a messy square
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
set.seed(2)
N <- 20
x <- c(seq(0, 1, length.out = N), rep(1, N), seq(1, 0, length.out = N), rep(0, N)) / 2 + 0.25
y <- c(rep(0, N), seq(0, 1, length.out = N), rep(1, N), seq(1, 0, length.out = N)) / 2 + 0.25

x <- x + runif(length(x), -0.02, 0.02)
y <- y + runif(length(x), -0.02, 0.02)

x[1] <- y[1] <- x[N*4] <- y[N*4] <- 0.25

grid::grid.newpage()
grid::grid.lines(x, y)



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Capture the simplification step by step.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
xs <- list()
ys <- list()

for (n in (4*N):5) {
  res <- visvalingam(x, y, n)
  xs <- c(xs, list(res$x))
  ys <- c(ys, list(res$y))
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Save frames 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for (i in seq_along(xs)) {
  filename <- sprintf("jnk/%03i.png", i)
  png(filename)
    grid.lines(xs[[i]], ys[[i]], gp = gpar(lwd= 4))
  dev.off()
}
# ffmpeg -y -framerate 15 -pattern_type glob -i '*.png' -c:v libx264 -pix_fmt yuv420p -s 800x800 anim.mp4