voronoifortune (version 1.0)

voronoi: Voronoi Tessellation and Delaunay Triangulation

Description

Voronoi tessellation and Delaunay triangulation are performed simultaneously with the Fortune (1987) algorithm.

Usage

voronoi(X, sorted = FALSE, debug = FALSE)
# S3 method for voronoi
print(x, ...)

Value

voronoi() returns an object of class "voronoi" with four elements:

Triplets

a three-column matrix of integers giving the triangles of the Delaunay triangulation (indices from the original data);

Vertices

a two-column matrix of reals giving the coordinates of the vertices of the Voronoi tessellation;

Edges

a two-column matrix of integers giving the edges of the tessellation (row indices of the vertices in the previous matrix);

Lines

a description of the previous edges (some of them are semi-infinite indicated by a 0 in the matrix Edges; this is used to draw the dashed lines by plot.voronoi).

Arguments

X

a two-column matrix with the coordinates.

sorted

a logical: are the coordinates already sorted in increasing order first by y, then by x? (See order.) If TRUE, this shortens the running time.

debug

a logical: if TRUE, some details of the computation are printed for debugging.

x

an object of class "voronoi".

...

(unused).

Author

Emmanuel Paradis, Steven Fortune

References

Fortune, S. (1987) A sweepline algorithm for Voronoi diagrams. Algorithmica, 2, 153--174. tools:::Rd_expr_doi("10.1007/BF01840357").

See Also

plot.voronoi

Examples

n <- 100L
xx <- runif(n)
yy <- runif(n)
X <- cbind(xx, yy)
(res <- voronoi(X))
str(res)

### show the circumcircle of each Delaunay triangle ###

n <- 10L
X <- cbind(runif(n), runif(n))
res <- voronoi(X)

## if 12 windows not enough par(ask) is set to TRUE
op <- par(mfrow = c(3, 4), ask = interactive())
## to show the circle as disk:
col <- rgb(.5, .5, 1, .3)

## for each triangle:
for (i in 1:nrow(res$Triplets)) {
    plot(res, X = X, voronoi = FALSE, main = i)
    ## get the 3 vertices of the triangle:
    P <- res$Triplets[i, ]
    ## center the coordinates on the 1st vertex:
    ## (B and C are the new coordinates of the 2nd
    ## and 3rd vertices)
    A <- X[P[1], ]
    B <- X[P[2], ] - A
    C <- X[P[3], ] - A
    ## the coordinates of the center of the circumcircle are:
    ## (https://en.wikipedia.org/wiki/Circumcircle)
    cr <- c(C[2] * sum(B^2) - B[2] * sum(C^2),
            B[1] * sum(C^2) - C[1] * sum(B^2))
    cr <- cr / (2 * (B[1] * C[2] - B[2] * C[1]))
    ## since the circle intersects with the new origin,
    ## the radius is simply:
    r <- sqrt(sum(cr^2))
    ## translate back to the original coordinate system:
    cr <- cr + A
    ## draw the circumcircle:
    symbols(cr[1], cr[2], circles = r, inches = FALSE, add = TRUE,
            fg = col, bg = col)
    ## test numerically that no points are inside the circumcircle:
    ## 1/ build a matrix with the coordinates of the center
    ## and all the vertices:
    M <- rbind(cr, X)
    ## 2/ compute the Euclidean distances, we keep only the n first
    ## values which are the distances from the center to the n vertices:
    Dis <- dist(M)[1:n]
    ## 3/ all other points than the 3 in P should be farther
    ## to the center:
    stopifnot(all(Dis[-P] > r))
}

par(op)