Delaunay triangulation

Delaunay triangulation (or tessellation) of a set of points.

Usage

delaunay(
  points,
  atinfinity = FALSE,
  degenerate = FALSE,
  exteriorEdges = FALSE,
  elevation = FALSE
)

Arguments

points: the points given as a matrix, one point per row
atinfinity: Boolean, whether to include a point at infinity
degenerate: Boolean, whether to include degenerate tiles
exteriorEdges: Boolean, for dimension 3 only, whether to return the exterior edges (see below)
elevation: Boolean, only for three-dimensional points; if TRUE, the function performs an elevated Delaunay triangulation (also called 2.5D Delaunay triangulation), using the third coordinate of a point as its elevation; see the example

Value

If the function performs an elevated Delaunay tessellation, then the returned value is a list with four fields: mesh, edges, volume, and surface. The mesh field is an object of class mesh3d, ready for plotting with the rgl package. The edges field is an integer matrix which provides the indices of the vertices of the edges, and an indicator of whether an edge is a border edge; this matrix is obtained with vcgGetEdge. The volume field provides the sum of the volumes under the Delaunay triangles, that is to say the total volume under the triangulated surface. Finally, the surface field provides the sum of the areas of the Delaunay triangles, thus this is an approximate value of the area of the surface that is triangulated. The elevated Delaunay tessellation is built with the help of the interp package. Otherwise, the function returns the Delaunay tessellation with many details, in a list. This list contains five fields:

vertices: the vertices (or sites) of the tessellation; these are the points passed to the function
tiles: the tiles of the tessellation (triangles in dimension 2, tetrahedra in dimension 3)
tilefacets: the facets of the tiles of the tessellation
mesh: a 'rgl' mesh (mesh3d object)
edges: a two-columns integer matrix representing the edges, each row represents an edge; the two integers of a row are the indices of the two points which form the edge.

In dimension 3, the list contains an additional field exteriorEdges if you set exteriorEdges = TRUE. This is the list of the exterior edges, represented as Edge3 objects. This field is involved in the function plotDelaunay3D. The vertices field is a list with the following fields:

id: the id of the vertex; this is nothing but the index of the corresponding point passed to the function
neighvertices: the ids of the vertices of the tessellation connected to this vertex by an edge
neightilefacets: the ids of the tile facets this vertex belongs to
neightiles: the ids of the tiles this vertex belongs to

The tiles field is a list with the following fields:

id: the id of the tile
simplex: a list describing the simplex (that is, the tile); this list contains four fields: vertices, a hash giving the simplex vertices and their id, circumcenter, the circumcenter of the simplex, circumradius, the circumradius of the simplex, and volume, the volume of the simplex
facets: the ids of the facets of this tile
neighbors: the ids of the tiles adjacent to this tile
family: two tiles have the same family if they share the same circumcenter; in this case the family is an integer, and the family is NA for tiles which do not share their circumcenter with any other tile
orientation: 1 or -1, an indicator of the orientation of the tile

The tilefacets field is a list with the following fields:

id: the id of this tile facet
subsimplex: a list describing the subsimplex (that is, the tile facet); this list is similar to the simplex list of tiles
facetOf: one or two ids, the id(s) of the tile this facet belongs to
normal: a vector, the normal of the tile facet
offset: a number, the offset of the tile facet

Note

The package provides the functions plotDelaunay2D to plot a 2D Delaunay tessellation and plotDelaunay3D to plot a 3D Delaunay tessellation. But there is no function to plot an elevated Delaunay tessellation; the examples show how to plot such a Delaunay tessellation.

Examples

library(tessellation)
points <- rbind(
 c(0.5,0.5,0.5),
 c(0,0,0),
 c(0,0,1),
 c(0,1,0),
 c(0,1,1),
 c(1,0,0),
 c(1,0,1),
 c(1,1,0),
 c(1,1,1)
)
del <- delaunay(points)
del$vertices[[1]]
#> $id
#> [1] 1
#> 
#> $point
#> [1] 0.5 0.5 0.5
#> 
#> $neighvertices
#> [1] 2 3 4 5 6 7 8 9
#> 
#> $neightilefacets
#>  [1]  1  3  4  5  7  9 10 11 13 15 16 17 20 21 23 25 26 29
#> 
#> $neightiles
#>  [1]  1  2  3  4  5  6  7  8  9 10 11 12
#> 
del$tiles[[1]]
#> $vertices
#> [1] 8 1 6 2
#> 
#> $id
#> [1] 1
#> 
#> $simplex
#> $simplex$vertices
#> <hash> containing 4 key-value pair(s).
#>   1 : 0.5 0.5 0.5
#>   2 : 0 0 0
#>   6 : 1 0 0
#>   8 : 1 1 0
#> 
#> $simplex$circumcenter
#> [1]  0.50  0.50 -0.25
#> 
#> $simplex$circumradius
#> [1] 0.75
#> 
#> $simplex$volume
#> [1] 0.08333333
#> 
#> 
#> $facets
#> [1] 1 2 3 4
#> 
#> $neighbors
#> [1] 3 2 5
#> 
#> $family
#> [1] 0
#> 
#> $orientation
#> [1] -1
#> 
del$tilefacets[[1]]
#> $id
#> [1] 1
#> 
#> $subsimplex
#> $subsimplex$vertices
#> <hash> containing 3 key-value pair(s).
#>   1 : 0.5 0.5 0.5
#>   2 : 0 0 0
#>   6 : 1 0 0
#> 
#> $subsimplex$circumcenter
#> [1] 0.500 0.125 0.125
#> 
#> $subsimplex$circumradius
#> [1] 0.5303301
#> 
#> $subsimplex$volume
#> [1] 0.3535534
#> 
#> 
#> $facetOf
#> [1] 1 3
#> 
#> $normal
#> [1]  0.0000000  0.7071068 -0.7071068
#> 
#> $offset
#> [1] 0
#> 

# an elevated Delaunay tessellation ####
f <- function(x, y){
  dnorm(x) * dnorm(y)
}
x <- y <- seq(-5, 5, length.out = 50)
grd <- expand.grid(x = x, y = y) # grid on the xy-plane
points <- as.matrix(transform( # data (x_i, y_i, z_i)
  grd, z = f(x, y)
))
del <- delaunay(points, elevation = TRUE)
del[["volume"]] # close to 1, as expected
#> [1] 0.9999988
# plotting
library(rgl)
mesh <- del[["mesh"]]
open3d(windowRect = c(100, 100, 612, 356), zoom = 0.6)
aspect3d(1, 1, 20)
shade3d(mesh, color = "limegreen")
wire3d(mesh)


# another elevated Delaunay triangulation, to check the correctness of
#   the calculated surface ####
library(Rvcg)
library(rgl)
cap <- vcgSphericalCap(angleRad = pi/2, subdivision = 3)
open3d(windowRect = c(100, 100, 612, 356), zoom = 0.6)
shade3d(cap, color = "lawngreen")
wire3d(cap)

# exact value of the surface of the spherical cap:
R <- 1
h <- R * (1 - sin(pi/2/2))
2 * pi * R * h
#> [1] 1.840302
# our approximation:
points <- t(cap$vb[-4, ]) # the points on the spherical cap
del <- delaunay(points, elevation = TRUE)
del[["surface"]]
#> [1] 1.832746
# try to increase `subdivision` in `vcgSphericalCap` to get a
#   better approximation of the true value
# note that 'Rvcg' returns the same result as ours:
vcgArea(cap)
#> [1] 1.832746
# let's check the volume as well:
pi * h^2 * (R - h/3) # true value
#> [1] 0.2431939
del[["volume"]]
#> [1] 0.2405314
# there's a warning with 'Rvcg':
tryCatch(vcgVolume(cap), warning = function(w) message(w))
#> Warning: Mesh is not watertight! USE RESULT WITH CARE!
#> Error in invokeRestart("muffleWarning"): no 'restart' 'muffleWarning' found
suppressWarnings({vcgVolume(cap)})
#> [1] 0.2405313

Usage

Arguments

Value

Note

See also

Examples