Nothing
#' @title 2D Delaunay triangulation
#' @description Performs a (constrained) Delaunay triangulation of a set of
#' 2d points.
#'
#' @param points a numeric matrix with two or three columns (three colums for
#' an elevated Delaunay triangulation)
#' @param edges the edges for the constrained Delaunay triangulation,
#' an integer matrix with two columns; \code{NULL} for no constraint
#' @param elevation Boolean, whether to perform an elevated Delaunay
#' triangulation (also known as 2.5D Delaunay triangulation)
#'
#' @return A list. There are three possibilities.
#' #' \itemize{
#' \item \strong{If the dimension is 2} and \code{edges=NULL},
#' the returned value is a list with three fields:
#' \code{vertices}, \code{mesh} and \code{edges}.
#' The \code{vertices} field contains the given vertices.
#' The \code{mesh} field is an object of
#' class \code{\link[rgl]{mesh3d}}, ready for plotting with the
#' \strong{rgl} package.
#' The \code{edges} field provides the indices of the vertices of the
#' edges, given by the first two columns of a three-columns integer
#' matrix.
#' The third column, named \code{border}, only contains some
#' zeros and some ones; a border (exterior) edge is labelled by a
#' \code{1}.
#' \item \strong{If the dimension is 2} and \code{edges} is not
#' \code{NULL}, the returned value is a list with four fields:
#' \code{vertices}, \code{mesh}, \code{edges}, and \code{constraints}.
#' The \code{vertices} field contains the vertices of the triangulation.
#' They coincide with the given vertices if the constraint edges do not
#' intersect; otherwise there are the intersections in addition to the
#' given vertices.
#' The \code{mesh} and \code{edges} fields are similar to the previous
#' case, the unconstrained Delaunay triangulation.
#' The \code{constraints} field is an integer matrix with
#' two columns, it represents the constraint edges. They are not the
#' same as the ones provided by the user if these ones intersect.
#' If they do not intersect, then in general these are the same,
#' but not always, in some rare corner cases.
#' \item \strong{If} \code{elevation=TRUE}, the returned value is a list with
#' five fields: \code{vertices}, \code{mesh}, \code{edges},
#' \code{volume}, and \code{surface}.
#' The \code{vertices} field contains the given vertices.
#' The \code{mesh} field is an object of class
#' \code{\link[rgl]{mesh3d}}, ready for plotting with the
#' \strong{rgl} package. The \code{edges} field is similar to the
#' previous cases.
#' The \code{volume} field provides a number, the sum of the volumes
#' under the Delaunay triangles, that is to say the total volume
#' under the triangulated surface. Finally, the \code{surface} field
#' provides the sum of the areas of all triangles, thereby
#' approximating the area of the triangulated surface.
#' }
#'
#' @note The triangulation can depend on the order of the points; this is
#' shown in the examples.
#'
#' @importFrom rgl tmesh3d
#' @importFrom Rvcg vcgGetEdge
#'
#' @export
#'
#' @examples library(RCDT)
#' # random points in a square ####
#' set.seed(314)
#' library(uniformly)
#' square <- rbind(
#' c(-1, 1), c(1, 1), c(1, -1), c(-1, -1)
#' )
#' pts_in_square <- runif_in_cube(10L, d = 2L)
#' pts <- rbind(square, pts_in_square)
#' del <- delaunay(pts)
#' opar <- par(mar = c(0, 0, 0, 0))
#' plotDelaunay(
#' del, type = "n", xlab = NA, ylab = NA, asp = 1,
#' fillcolor = "random", luminosity = "light", lty_edges = "dashed"
#' )
#' par(opar)
#'
#' # the order of the points matters ####
#' # the Delaunay triangulation is not unique in general;
#' # it can depend on the order of the points
#' points <- cbind(
#' c(1, 2, 1, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 4, 3, 4),
#' c(1, 1, 2, 1, 2, 3, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4)
#' )
#' del <- delaunay(points)
#' opar <- par(mar = c(0, 0, 0, 0))
#' plotDelaunay(
#' del, type = "p", pch = 19, xlab = NA, ylab = NA, axes = FALSE,
#' asp = 1, lwd_edges = 2, lwd_borders = 3
#' )
#' par(opar)
#' # now we randomize the order of the points
#' set.seed(666L)
#' points2 <- points[sample.int(nrow(points)), ]
#' del2 <- delaunay(points2)
#' opar <- par(mar = c(0, 0, 0, 0))
#' plotDelaunay(
#' del2, type = "p", pch = 19, xlab = NA, ylab = NA, axes = FALSE,
#' asp = 1, lwd_edges = 2, lwd_borders = 3
#' )
#' par(opar)
#'
#' # a constrained Delaunay triangulation: outer and inner dodecagons ####
#' # points
#' nsides <- 12L
#' angles <- seq(0, 2*pi, length.out = nsides+1L)[-1L]
#' points <- cbind(cos(angles), sin(angles))
#' points <- rbind(points, points/1.5)
#' # constraint edges
#' indices <- 1L:nsides
#' edges_outer <- cbind(
#' indices, c(indices[-1L], indices[1L])
#' )
#' edges_inner <- edges_outer + nsides
#' edges <- rbind(edges_outer, edges_inner)
#' # constrained Delaunay triangulation
#' del <- delaunay(points, edges)
#' # plot
#' opar <- par(mar = c(0, 0, 0, 0))
#' plotDelaunay(
#' del, type = "n", fillcolor = "yellow", lwd_borders = 2, asp = 1,
#' axes = FALSE, xlab = NA, ylab = NA
#' )
#' par(opar)
#'
#' # another constrained Delaunay triangulation: a face ####
#' V <- read.table(
#' system.file("extdata", "face_vertices.txt", package = "RCDT")
#' )
#' E <- read.table(
#' system.file("extdata", "face_edges.txt", package = "RCDT")
#' )
#' del <- delaunay(
#' points = as.matrix(V)[, c(2L, 3L)], edges = as.matrix(E)[, c(2L, 3L)]
#' )
#' opar <- par(mar = c(0, 0, 0, 0))
#' plotDelaunay(
#' del, type="n", col_edges = NULL, fillcolor = "salmon",
#' col_borders = "black", col_constraints = "purple",
#' lwd_borders = 3, lwd_constraints = 3,
#' asp = 1, axes = FALSE, xlab = NA, ylab = NA
#' )
#' par(opar)
delaunay <- function(points, edges = NULL, elevation = FALSE){
stopifnot(isBoolean(elevation))
if(!is.matrix(points) || !is.numeric(points)){
stop(
"The `points` argument must be a numeric matrix with two or three columns.",
call. = TRUE
)
}
dimension <- ncol(points)
if(!is.element(dimension, c(2L, 3L))){
stop(
"The `points` argument must be a numeric matrix with two or three columns.",
call. = TRUE
)
}
if(any(is.na(points))){
stop("Missing values are not allowed.", call. = TRUE)
}
if(elevation){
if(dimension != 3L){
stop(
"To get an elevated Delaunay tessellation (`elevation=TRUE`), ",
"you have to provide three-dimensional points.",
call. = TRUE
)
}
x <- points[, 1L]
y <- points[, 2L]
x <- (x - min(x)) / diff(range(x))
y <- (y - min(y)) / diff(range(y))
o <- order(round(x+y, 6L), y-x)
xy <- cbind(x, y)[o, ]
if(anyDuplicated(xy)){
stop("There are some duplicated points.", call. = TRUE)
}
Triangles <- Rcpp_delaunay(t(xy))
vertices <- points[o, ]
mesh <- tmesh3d(
vertices = t(vertices),
indices = Triangles,
homogeneous = FALSE
)
volumes_and_areas <- apply(Triangles, 2L, function(trgl){
trgl <- vertices[trgl, ]
c(
volume_under_triangle(trgl[, 1L], trgl[, 2L], trgl[, 3L]),
triangleArea(trgl[1L, ], trgl[2L, ], trgl[3L, ])
)
})
out <- list(
"vertices" = vertices,
"mesh" = mesh,
"edges" = `colnames<-`(
as.matrix(vcgGetEdge(mesh))[, -3L], c("v1", "v2", "border")
),
"volume" = sum(volumes_and_areas[1L, ]),
"surface" = sum(volumes_and_areas[2L, ])
)
attr(out, "elevation") <- TRUE
class(out) <- "delaunay"
return(out)
}
if(anyDuplicated(points)){
stop("There are some duplicated points.", call. = TRUE)
}
storage.mode(points) <- "double"
tpoints <- t(points)
if(is.null(edges)){
triangles <- Rcpp_delaunay(tpoints)
mesh <- tmesh3d(
vertices = rbind(tpoints, 0),
indices = triangles
)
Edges <- `colnames<-`(
as.matrix(vcgGetEdge(mesh))[, c(1L, 2L, 4L)], c("v1", "v2", "border")
)
out <- list(
"vertices" = points,
"mesh" = mesh,
"edges" = Edges
)
}else{
if(!is.matrix(edges) || !is.numeric(edges) || ncol(edges) != 2L){
stop(
"The `edges` argument must be an integer matrix with two columns.",
call. = TRUE
)
}
if(any(is.na(edges))){
stop("Missing values in `edges` are not allowed.", call. = TRUE)
}
storage.mode(edges) <- "integer"
stopifnot(all(edges >= 1L))
stopifnot(all(edges <= nrow(points)))
edges <- t(apply(edges, 1L, sort))
if(anyDuplicated(edges)){
stop("There are some duplicated constraint edges.", call. = TRUE)
}
if(any(edges[, 1L] == edges[, 2L])){
stop("There are some invalid constraint edges.", call. = TRUE)
}
cpp <- Rcpp_constrained_delaunay(tpoints, t(edges))
triangles <- cpp[["triangles"]]
storage.mode(triangles) <- "integer"
Vertices <- cpp[["vertices"]]
if(ncol(triangles) == 0L){
mesh <- NULL
Edges <- `colnames<-`(
matrix(integer(0L), nrow = 0L, ncol = 3L), c("v1", "v2", "border")
)
}else{
mesh <- tmesh3d(
vertices = rbind(Vertices, 0),
indices = triangles
)
Edges <- `colnames<-`(
as.matrix(vcgGetEdge(mesh))[, c(1L, 2L, 4L)], c("v1", "v2", "border")
)
}
borderEdges <- cpp[["borderEdges"]]
storage.mode(borderEdges) <- "integer"
out <- list(
"mesh" = mesh,
"vertices" = t(Vertices),
"edges" = Edges,
"constraints" = t(borderEdges)
)
attr(out, "constrained") <- TRUE
}
class(out) <- "delaunay"
out
}
#' @title Area of Delaunay triangulation
#' @description Computes the area of a region subject to Delaunay triangulation.
#'
#' @param del an output of \code{\link{delaunay}} executed with
#' \code{elevation=FALSE}
#'
#' @return A number, the area of the region triangulated by the Delaunay
#' triangulation.
#' @export
#'
#' @examples library(RCDT)
#' # random points in a square ####
#' set.seed(666L)
#' library(uniformly)
#' square <- rbind(
#' c(-1, 1), c(1, 1), c(1, -1), c(-1, -1)
#' )
#' pts <- rbind(square, runif_in_cube(8L, d = 2L))
#' del <- delaunay(pts)
#' delaunayArea(del)
#'
#' # a constrained Delaunay triangulation: outer and inner squares ####
#' innerSquare <- rbind( # the hole
#' c(-1, 1), c(1, 1), c(1, -1), c(-1, -1)
#' ) # area: 4
#' outerSquare <- 2*innerSquare # area: 16
#' points <- rbind(innerSquare, outerSquare)
#' edges_inner <- rbind(c(1L, 2L), c(2L, 3L), c(3L, 4L), c(4L, 1L))
#' edges_outer <- edges_inner + 4L
#' edges <- rbind(edges_inner, edges_outer)
#' del <- delaunay(points, edges = edges)
#' delaunayArea(del) # 16-4
delaunayArea <- function(del){
stopifnot(inherits(del, "delaunay"))
if(isTRUE(attr(del, "elevation"))){
stop(
"This function is not conceived for elevated Delaunay triangulations.",
call. = TRUE
)
}
mesh <- del[["mesh"]]
if(is.null(mesh)){
stop(
"This Delaunay triangulation is empty.", call. = TRUE
)
}
triangles <- mesh[["it"]]
vertices <- del[["vertices"]]
ntriangles <- ncol(triangles)
areas <- numeric(ntriangles)
for(i in 1L:ntriangles){
points <- vertices[triangles[, i], ]
areas[i] <- triangleArea(points[1L, ], points[2L, ], points[3L, ])
}
sum(areas)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.