Voronoi function

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(print,delaunay)
+export(Voronoi)
 export(delaunay)
 export(mesh2d)
 export(plotDelaunay2D)
@@ -20,5 +21,6 @@ importFrom(rgl,material3d)
 importFrom(rgl,points3d)
 importFrom(rgl,tmesh3d)
 importFrom(rgl,triangles3d)
+importFrom(sets,pair)
 importFrom(utils,combn)
 useDynLib(delaunay, .registration=TRUE)

R/voronoi.R

@@ -0,0 +1,132 @@
+voronoiEdgeFromDelaunayEdge <- function(edgeIndex, Edges, Circumcenters) {
+  faces <- unlist(Edges[edgeIndex, c("f1", "f2")])
+  face1 <- faces[1L]
+  face2 <- faces[2L]
+  c1 <- Circumcenters[face1, ]
+  if(is.na(face2)){
+    return(cbind(c1, NA_real_))
+  }
+  c2 <- Circumcenters[face2, ]
+  if(isTRUE(all.equal(c1, c2))){
+    return(NULL)
+  }
+  cbind(c1, c2)
+}
+
+vertexNeighborEdges <- function(vertexId, Edges) {
+  i1 <- Edges[["i1"]]
+  i2 <- Edges[["i2"]]
+  which((i1 == vertexId) | (i2 == vertexId))
+}
+
+VoronoiCell <- function(mesh, facetsQuotienter, edgeTransformer){
+  Edges <- mesh[["edges"]]
+  Faces <- mesh[["faces"]]
+  Circumcenters <- Faces[, c("ccx", "ccy")]
+  function(vertexId){
+    delaunayEdges <- facetsQuotienter(
+      unname(vertexNeighborEdges(vertexId, Edges))
+    )
+    voronoiEdges <- Filter(
+      Negate(is.null), lapply(delaunayEdges, function(dedge){
+        voronoiEdgeFromDelaunayEdge(dedge, Edges, Circumcenters)
+      })
+    )
+    nedges <- length(voronoiEdges)
+    bounded <- nedges > 0L
+    while(nedges > 0L){
+      bounded <- bounded && !anyNA(voronoiEdges[[nedges]])
+      nedges <- nedges - 1L
+    }
+    cell <- edgeTransformer(voronoiEdges)
+    attr(cell, "bounded") <- bounded
+    cell
+  }
+}
+
+#' @importFrom sets pair
+#' @noRd
+zip <- function(matrix, list) {
+  lapply(seq_len(nrow(matrix)), function(i) {
+    pair(site = matrix[i, ], cell = list[[i]])
+  })
+}
+
+Voronoi0 <- function(cellGetter, mesh) {
+  vertices <- mesh[["vertices"]]
+  nvertices <- nrow(vertices)
+  bounded <- logical(nvertices)
+  L <- lapply(1L:nvertices, function(i) {
+    cell <- cellGetter(i)
+    bounded[i] <<- attr(cell, "bounded")
+    cell
+  })
+  out <- zip(vertices, L)
+  attr(out, "nbounded") <- sum(bounded)
+  out
+}
+
+#' @title Voronoï diagram
+#' @description Returns the Voronoï tessellation corresponding to a 2D 
+#'   Delaunay triangulation. If the Delaunay triangulation is constrained, 
+#'   the output is not a true constrained Voronoï tessellation.
+#' 
+#' @param triangulation a 2D Delaunay triangulation obtained with the 
+#'   \code{\link{delaunay}} function
+#'
+#' @return The Voronoï diagram given as list of pairs; each pair is made of 
+#'   one site (a vertex of the Delaunay triangulation) and one Voronoï cell, 
+#'   a polygon given as numeric matrix with two columns (that can be plotted 
+#'   with \code{\link[graphics:polygon]{polygon}}). 
+#' @export
+#' @seealso \code{\link{plotVoronoi}}
+#'
+#' @examples
+#' library(delaunay)
+#' # make the vertices
+#' nsides <- 6L
+#' angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
+#' outer_points <- cbind(cos(angles), sin(angles))
+#' inner_points <- outer_points / 4
+#' nsides <- 12L
+#' angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
+#' middle_points <- cbind(cos(angles), sin(angles)) / 2
+#' points <- rbind(outer_points, inner_points, middle_points)
+#' angles <- angles + pi/24
+#' middle_points <- cbind(cos(angles), sin(angles)) / 3
+#' points <- rbind(points, middle_points)
+#' middle_points <- cbind(cos(angles), sin(angles)) / 1.5
+#' points <- rbind(points, middle_points)
+#' # constraint edges
+#' indices <- 1L:6L
+#' edges <- cbind(
+#'   indices, c(indices[-1L], indices[1L])
+#' )
+#' edges <- rbind(edges, edges + 6L)
+#' ## | constrained Delaunay triangulation 
+#' del <- delaunay(points, constraints = edges)
+#' opar <- par(mar = c(0,0,0,0))
+#' plotDelaunay2D(
+#'   del, type = "n", xlab = NA, ylab = NA, axes = FALSE, asp = 1,
+#'   luminosity = "dark", col_borders = "black", lwd_borders = 3
+#' )
+#' par(opar)
+#' ## | corresponding Voronoï diagram
+#' vor <- Voronoi(del)
+Voronoi <- function(triangulation) {
+  if(!inherits(triangulation, "delaunay")){
+    stop(
+      "The argument `triangulation` must be an output of the `delaunay` function.",
+      call. = TRUE
+    )
+  }
+  dimension <- attr(triangulation, "dimension")
+  if(dimension != 2){
+    stop(
+      sprintf("Invalid dimension (%s instead of 2).", dimension),
+      call. = TRUE
+    )
+  }
+  mesh <- triangulation[["mesh"]]
+  Voronoi0(VoronoiCell(mesh, identity, identity), mesh)  
+}

inst/essais/voronoi2.R

@@ -1,31 +1,28 @@
 library(delaunay)
 library(cgalMeshes)
 
+# make the vertices
 nsides <- 6L
 angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
 outer_points <- cbind(cos(angles), sin(angles))
 inner_points <- outer_points / 4
-
 nsides <- 12L
 angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
 middle_points <- cbind(cos(angles), sin(angles)) / 2
 points <- rbind(outer_points, inner_points, middle_points)
-
 angles <- angles + pi/24
 middle_points <- cbind(cos(angles), sin(angles)) / 3
 points <- rbind(points, middle_points)
 middle_points <- cbind(cos(angles), sin(angles)) / 1.5
 points <- rbind(points, middle_points)
-
-
 # constraint edges
 indices <- 1L:6L
 edges <- cbind(
   indices, c(indices[-1L], indices[1L])
 )
 edges <- rbind(edges, edges + 6L)
 # constrained Delaunay triangulation
-d <- delaunay(points, constraints = edges)
+del <- delaunay(points, constraints = edges)
 
 #####
 

inst/essais/voronoi4.R

@@ -79,9 +79,10 @@ VoronoiCell <- function(mesh, facetsQuotienter, edgeTransformer){
   }
 }
 
+#' @importFrom sets pair
 zip <- function(matrix, list) {
   lapply(seq_len(nrow(matrix)), function(i) {
-    sets::pair(site = matrix[i, ], cell = list[[i]])
+    pair(site = matrix[i, ], cell = list[[i]])
   })
 }
 
@@ -99,8 +100,68 @@ Voronoi0 <- function(cellGetter, mesh) {
   out
 }
 
-Voronoi <- function(tesselation) {
-  mesh <- tesselation[["mesh"]]
+#' @title Voronoï diagram
+#' @description Returns the Voronoï tessellation corresponding to a 2D 
+#'   Delaunay triangulation. If the Delaunay triangulation is constrained, 
+#'   the output is not a true constrained Voronoï tessellation.
+#' 
+#' @param triangulation a 2D Delaunay triangulation obtained with the 
+#'   \code{\link{delaunay}} function
+#'
+#' @return The Voronoï diagram given as list of pairs; each pair is made of 
+#'   one site (a vertex of the Delaunay triangulation) and one Voronoï cell, 
+#'   a polygon given as numeric matrix with two columns (that can be plotted 
+#'   with \code{\link[graphics:polygon]{polygon}}). 
+#' @export
+#' @seealso \code{\link{plotVoronoi}}
+#'
+#' @examples
+#' library(delaunay)
+#' # make the vertices
+#' nsides <- 6L
+#' angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
+#' outer_points <- cbind(cos(angles), sin(angles))
+#' inner_points <- outer_points / 4
+#' nsides <- 12L
+#' angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
+#' middle_points <- cbind(cos(angles), sin(angles)) / 2
+#' points <- rbind(outer_points, inner_points, middle_points)
+#' angles <- angles + pi/24
+#' middle_points <- cbind(cos(angles), sin(angles)) / 3
+#' points <- rbind(points, middle_points)
+#' middle_points <- cbind(cos(angles), sin(angles)) / 1.5
+#' points <- rbind(points, middle_points)
+#' # constraint edges
+#' indices <- 1L:6L
+#' edges <- cbind(
+#'   indices, c(indices[-1L], indices[1L])
+#' )
+#' edges <- rbind(edges, edges + 6L)
+#' ## | constrained Delaunay triangulation 
+#' del <- delaunay(points, constraints = edges)
+#' opar <- par(mar = c(0,0,0,0))
+#' plotDelaunay2D(
+#'   del, type = "n", xlab = NA, ylab = NA, axes = FALSE, asp = 1,
+#'   luminosity = "dark", col_borders = "black", lwd_borders = 3
+#' )
+#' par(opar)
+#' ## | corresponding Voronoï diagram
+#' vor <- Voronoi(del)
+Voronoi <- function(triangulation) {
+  if(!inherits(triangulation, "delaunay")){
+    stop(
+      "The argument `triangulation` must be an output of the `delaunay` function.",
+      call. = TRUE
+    )
+  }
+  dimension <- attr(triangulation, "dimension")
+  if(dimension != 2){
+    stop(
+      sprintf("Invalid dimension (%s instead of 2).", dimension),
+      call. = TRUE
+    )
+  }
+  mesh <- triangulation[["mesh"]]
   Voronoi0(VoronoiCell(mesh, identity, identity), mesh)  
 }