| Title: | Delaunay Triangulation and Dirichlet (Voronoi) Tessellation |
|---|---|
| Description: | Calculates the Delaunay triangulation and the Dirichlet or Voronoi tessellation (with respect to the entire plane) of a planar point set. Plots triangulations and tessellations in various ways. Clips tessellations to sub-windows. Calculates perimeters of tessellations. Summarises information about the tiles of the tessellation. Calculates the centroidal Voronoi (Dirichlet) tessellation using Lloyd's algorithm. |
| Authors: | Rolf Turner |
| Maintainer: | Rolf Turner <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.0-4 |
| Built: | 2024-10-31 22:18:06 UTC |
| Source: | CRAN |
Calculates the centroidal Voronoi (Dirichlet) tessellation using Lloyd's algorithm.
cvt(object, stopcrit = c("change", "maxit"), tol = NULL, maxit = 100, verbose = FALSE)cvt(object, stopcrit = c("change", "maxit"), tol = NULL, maxit = 100, verbose = FALSE)
object |
An object of class either |
stopcrit |
Text string specifying the stopping criterion for the algorithm. If
this is |
tol |
The tolerance used when the stopping criterion is |
maxit |
The maximum number of iterations to perform when the stopping criterion
is |
verbose |
Logical scalar. If |
The algorithm iteratively tessellates a set of points and then replaces the points with the centroids of the tiles associated with those points. “Eventually” (at convergence) points and the centroids of their associated tiles coincide.
A list with components:
centroids |
A data frame with columns |
tiles |
An object of class |
This function was added to the deldir package at the
suggestion of Dr. Michaël Aupetit, Senior Scientist at the
Qatar Computing Research Institute, Hamad Bin Khalifa University.
Rolf Turner [email protected]
https://en.wikipedia.org/wiki/Lloyd's_algorithm
Lloyd, Stuart P. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory 28 (2), pp. 129–137, doi:10.1109/TIT.1982.1056489.
## Not run: # Takes too long. set.seed(42) x <- runif(20) y <- runif(20) dxy <- deldir(x,y,rw=c(0,1,0,1)) cxy1 <- cvt(dxy,verb=TRUE) plot(cxy1$tiles) with(cxy1$centroids,points(x,y,pch=20,col="red")) cxy2 <- cvt(dxy,stopcrit="m",verb=TRUE) plot(cxy2$tiles) with(cxy2$centroids,points(x,y,pch=20,col="red")) # Visually indistinguishable from the cxy1 plot. # But ... all.equal(cxy1$centroids,cxy2$centroids) # Not quite. cxy3 <- cvt(dxy,stopcrit="m",verb=TRUE,maxit=250) all.equal(cxy1$centroids,cxy3$centroids) # Close, but no cigar. cxy4 <- cvt(dxy,verb=TRUE,tol=1e-14) cxy5 <- cvt(dxy,stopcrit="m",verb=TRUE,maxit=600) all.equal(cxy4$centroids,cxy5$centroids) # TRUE # Takes a lot of iterations or a really small tolerance # to get "good" convergence. But this is almost surely # of no practical importance. txy <- tile.list(dxy) cxy6 <- cvt(txy) all.equal(cxy6$centroids,cxy1$centroids) # TRUE ## End(Not run)## Not run: # Takes too long. set.seed(42) x <- runif(20) y <- runif(20) dxy <- deldir(x,y,rw=c(0,1,0,1)) cxy1 <- cvt(dxy,verb=TRUE) plot(cxy1$tiles) with(cxy1$centroids,points(x,y,pch=20,col="red")) cxy2 <- cvt(dxy,stopcrit="m",verb=TRUE) plot(cxy2$tiles) with(cxy2$centroids,points(x,y,pch=20,col="red")) # Visually indistinguishable from the cxy1 plot. # But ... all.equal(cxy1$centroids,cxy2$centroids) # Not quite. cxy3 <- cvt(dxy,stopcrit="m",verb=TRUE,maxit=250) all.equal(cxy1$centroids,cxy3$centroids) # Close, but no cigar. cxy4 <- cvt(dxy,verb=TRUE,tol=1e-14) cxy5 <- cvt(dxy,stopcrit="m",verb=TRUE,maxit=600) all.equal(cxy4$centroids,cxy5$centroids) # TRUE # Takes a lot of iterations or a really small tolerance # to get "good" convergence. But this is almost surely # of no practical importance. txy <- tile.list(dxy) cxy6 <- cvt(txy) all.equal(cxy6$centroids,cxy1$centroids) # TRUE ## End(Not run)
This function computes the Delaunay triangulation (and hence the Dirichlet or Voronoi tessellation) of a planar point set according to the second (iterative) algorithm of Lee and Schacter — see References. The triangulation is made to be with respect to the whole plane by “suspending” it from so-called ideal points (-Inf,-Inf), (Inf,-Inf) (Inf,Inf), and (-Inf,Inf). The triangulation is also enclosed in a finite rectangular window.
deldir(x, y, z=NULL, rw=NULL, eps=1e-09, sort=TRUE, plot=FALSE, round=TRUE,digits=6, id=NULL, ...)deldir(x, y, z=NULL, rw=NULL, eps=1e-09, sort=TRUE, plot=FALSE, round=TRUE,digits=6, id=NULL, ...)
x, y
|
These arguments specify the coordinates of the point set being
triangulated/tessellated. Argument |
z |
Optional argument specifying “auxiliary” values or
“tags” associated with the respective points. (See
Notes on “tags”.) This argument may be a vector or
factor whose entries constitute these tags, or it may be a text
string naming such a vector or factor. If |
rw |
The coordinates of the corners of the rectangular window enclosing the triangulation, in the order (xmin, xmax, ymin, ymax). Any data points outside this window are discarded. If this argument is omitted, it defaults to values given by the range of the data, plus and minus 10 percent. |
eps |
A value of epsilon used in testing whether a quantity is zero, mainly in the context of whether points are collinear. If anomalous errors arise, it is possible that these may averted by adjusting the value of eps upward or downward. |
sort |
Logical argument; if |
plot |
Logical argument; if |
round |
Logical scalar. Should the data stored in the returned value
be rounded to |
digits |
The number of decimal places to which all numeric values in the
returned list should be rounded. Defaults to 6. Ignored if
|
id |
Optional vector of the same length as The |
... |
Auxiliary arguments |
This package had its origins (way back in the mists of time!) as an Splus library section named “delaunay”. That library section in turn was a re-write of a stand-alone Fortran program written in 1987/88 while the author was with the Division of Mathematics and Statistics, CSIRO, Sydney, Australia. This program was an implementation of the second (iterative) Lee-Schacter algorithm. The stand-alone Fortran program was re-written as an Splus function (which called upon dynamically loaded Fortran code) while the author was visiting the University of Western Australia, May, 1995.
Further revisions were made December 1996. The author gratefully acknowledges the contributions, assistance, and guidance of Mark Berman, of D.M.S., CSIRO, in collaboration with whom this project was originally undertaken. The author also acknowledges much useful advice from Adrian Baddeley, formerly of D.M.S., CSIRO (now Professor of Statistics at Curtin University). Daryl Tingley of the Department of Mathematics and Statistics, University of New Brunswick, provided some helpful insight. Special thanks are extended to Alan Johnson, of the Alaska Fisheries Science Centre, who supplied two data sets which were extremely valuable in tracking down some errors in the code.
Don MacQueen, of Lawrence Livermore National Lab, wrote an Splus driver function for the old stand-alone version of this software. That driver, which was available on Statlib, was deprecated in favour of the Statlib package “delaunay”. Don also collaborated in the preparation of the latter package. It is not clear to me whether the “delaunay” package, or indeed Statlib (or indeed Splus!) still exist.
See the ChangeLog for information about further revisions
and bug-fixes.
A list (of class deldir), invisible if plot=TRUE, with components:
delsgs |
A data frame with 6 columns. The first 4 entries of each row are the
coordinates of the points joined by an edge of a Delaunay triangle,
in the order |
dirsgs |
A data frame with 10 columns. The first 4 entries of each
row are the coordinates of the endpoints of one the edges of a
Dirichlet tile, in the order The ninth and tenth entries, in columns named The entries of columns Note that the entry in column |
summary |
a data frame with 9, 10 or 11 columns and
Note that the factor of 1/3 associated with the del.area column arises because each triangle occurs three times — once for each corner. |
n.data |
the number of points in the set which was triangulated, with any
duplicate points eliminated. It is the same as the number of rows
of |
del.area |
the area of the convex hull of the set of points being triangulated,
as formed by summing the |
dir.area |
the area of the rectangular window enclosing the points being triangulated,
as formed by summing the |
rw |
the specification of the corners of the rectangular window enclosing the data, in the order (xmin, xmax, ymin, ymax). |
ind.orig |
A vector of the indices of the points (x,y) in the
set of coordinates initially supplied to |
If plot=TRUE a plot of the triangulation and/or tessellation
is produced or added to an existing plot.
The protocol for extracting the and coordinates
from the arguments x and y is a bit complicated
and confusing. It is designed to handle a number of different
desiderata and to accommodate various feature-requests that users
have made over the years. Basically the protocol is:
If x is a numeric vector and y is a numeric
vector then x is used as the -coordinates and
y is used as the -coordinates.
If x is a matrix, a data frame, or a generic list),
and y is a numeric vector, then the -coordinates
are sought amongst the components of x and y
is used as the -coordinates.
If x is a matrix, a data frame, or a generic list
and y is not specified or cannot be found, then both the
-coordinates and -coordinates are sought amongst
the components of x.
If x an object of class "ppp" then both the
-coordinates and -coordinates are taken from the
components of x. If y is specified, it is ignored
(with a warning).
If x is a numeric vector and y is not specified
or cannot be found, then an error is thrown.
A few more details:
If x is of class "ppp" then it will definitely
have components named "x" and "y".
If x is a generic list, it must have a
component named "x" (otherwise an error is thrown),
and the -coordinates are set equal to this component.
If y is not specified or cannot be found, then a "y"
component of x is sought. If such a component exists
then the -coordinates are set equal to this component.
Otherwise an error is thrown).
If x is a matrix or a data frame, the protocol gets
a bit more intricate.
If x has a column named "x" then this column
is taken to be the -coordinates.
Otherwise the -coordinates are taken to be the
first column of x that is not named "y" or
znm (where znm is the name of the object providing
the “tags”, if “tags” have been specified).
If there is no such first column (e.g. if there are only
two columns and these have names "y" and znm)
then an error is thrown.
If y is not specified or cannot be found, and
if x has a column named "y" then this column is
taken to be the -coordinates.
Otherwise, in this situation, the -coordinates
are taken to be the first column of x that is
not named "x" or znm and is not equal to the
column previously selected to be the x-coordinates.
If there is no such first column (e.g. if there are only
two columns and these have names "x" and znm),
then an error is thrown.
Got all that? :-) If these instructions seem rather
bewildering (and indeed they are!) just keep things simple and make
calls like deldir(x,y) where x and y are numeric
vectors that have been previously assigned in the global environment.
z
If argument x is a data structure (rather than a numeric
vector) and is not an object of class "ppp" then
z, if specified and not found, is searched for in x.
If x is of class "ppp" then what happens depends
on whether z was specified or left to take its default value
of NULL. In the former case, z takes the specified
value. In the latter case the value of "z" is taken from
the marks of x provided that x is indeed a marked
point pattern and that the marks are atomic (essentially
provided that the marks are not a data frame). Otherwise z
is left NULL, i.e. there are no “tags” associated
with the points.
The “tags” are simply values that are associated in some way
with the data points and hence with the tiles of the tessellation
produced. They DO NOT affect the tessellation. In previous
versions of this package (0.2-10 and earlier) the entries of z
were referred to as “weights”. This terminology has been
changed since it is misleading. The tessellation produced when
a z argument is supplied is the same as is it would be
if there were no z argument (i.e. no “weights”).
The deldir package DOES NOT do weighted tessellation.
It is difficult-to-impossible to determine a priori how much
memory needs to be allocated (in the Fortran code) for storing the
edges of the Delaunay triangles and Dirichlet tiles, and for storing
the “adjacency list” used by the Lee-Schacter algorithm.
In the code, an attempt is made to allocate sufficient storage.
If, during the course of running the algorithm, the amount of
storage turns out to be inadequate, the algorithm is halted, the
storage is incremented, and the algorithm is restarted (with an
informative message). This message may be suppressed by wrapping
the call to deldir() in suppressMessages().
In previous versions of this package, error traps were set in
the underlying Fortran code for 17 different errors. These were
identified by an error number which was passed back up the call stack
and finally printed out by deldir() which then invisibly
returned a NULL value. A glossary of the meanings of the
values of was provided in a file to be found in a file located in the
inst directory (“folder” if you are a Windoze weenie).
This was a pretty shaganappi system. Consequently, as of version
1.2-1, conversion to “proper” error trapping was implemented.
Such error trapping is effected via the rexit() subroutine
which is now available in R. (See “Writing R
Extensions”, section 6.2.1.)
Note that when an error is detected, deldir() now exits with
a genuine error, rather than returning NULL. The glossary
of the meanings of “error numbers” is now irrelevant and
has been removed from the inst directory.
An error trap that merits particular mention was introduced in
version 0.1-16 of deldir. This error trap relates to
“triangle problems”. It was drawn to my attention by Adam
Dadvar (on 18 December, 2018) that in some data sets collinearity
problems may cause the “triangle finding” procedure, used
by the algorithm to successively add new points to a tessellation,
to go into an infinite loop. A symptom of the collinearity is
that the vertices of a putative triangle appear not to be
in anticlockwise order irrespective of whether they are presented
in the order i, j, k or k, j, i. The result of this
anomaly is that the procedure keeps alternating between moving to
“triangle” i, j, k and moving to “triangle”
k, j, i, forever.
The error trap in question is set in trifnd, the triangle
finding subroutine. It detects such occurrences of “clockwise
in either orientation” vertices. The trap causes the deldir()
function to throw an error rather than disappearing into a black
hole.
When an error of the “triangle problems” nature occurs, a
possible remedy is to increase the value of the eps
argument of deldir(). (See the Examples.) There may
conceivably be other problems that lead to infinite loops and so I
put in another error trap to detect whether the procedure has
inspected more triangles than actually exist, and if so to throw
an error.
Note that the strategy of increasing the value of eps
is probably the appropriate response in most (if not all)
of the cases where errors of this nature arise. Similarly this
strategy is probably the appropriate response to the errors
Vertices of “triangle” are collinear and vertex 2 is not between 1 and 3. Error in circen.
Vertices of triangle are collinear. Error in dirseg.
Vertices of triangle are collinear. Error in dirout.
However it is impossible to be sure. The intricacy and numerical delicacy of triangulations is too great for anyone to be able to foresee all the possibilities that could arise.
If there is any doubt as to the appropriateness of the
“increase eps” strategy, users are advised to do
their best to explore the data set, graphically or by other means,
and thereby determine what is actually going on and why problems
are occurring.
The process for determining if points are duplicated
changed between versions 0.1-9 and 0.1-10. Previously there
was an argument frac for this function, which defaulted
to 0.0001. Points were deemed to be duplicates if the difference
in x-coordinates was less than frac times the width
of rw and y-coordinates was less than frac
times the height of rw. This process has been changed to
one which uses duplicated() on the data frame whose
columns are x and y.
As a result it may happen that points which were previously eliminated as duplicates will no longer be eliminated. (And possibly vice-versa.)
The components delsgs and summary of the value
returned by deldir() are now data frames rather than
matrices. The component summary was changed to allow the
“auxiliary” values z to be of arbitrary mode (i.e.
not necessarily numeric). The component delsgs was then
changed for consistency. Note that the other “matrix-like”
component dirsgs has been a data frame since time immemorial.
I would like to express my most warm and sincere thanks to Duncan
Murdoch (Emeritus Professor of Statistics, Western University) for
helping me, with incredible patience and forbearance, to straighten
out my thinking in respect of adjustments that I recently (October
2021) made to the argument processing protocol in the deldir()
function. Duncan provided numerous simple examples to demonstrate
when and how things were going wrong, and patiently explained to
me how I was getting one aspect of the protocol backwards.
Rolf Turner [email protected]
Lee, D. T. and Schacter, B. J. (1980). Two algorithms for constructing a Delaunay triangulation, International Journal of Computer and Information Sciences 9 (3), pp. 219 – 242.
Ahuja, N. and Schacter, B. J. (1983). Pattern Models. New York: Wiley.
plot.deldir(), tile.list(), triang.list()
x <- c(2.3,3.0,7.0,1.0,3.0,8.0) y <- c(2.3,3.0,2.0,5.0,8.0,9.0) # Let deldir() choose the rectangular window. dxy1 <- deldir(x,y) # User chooses the rectangular window. dxy2 <- deldir(x,y,rw=c(0,10,0,10)) # Put "dummy" points at the corners of the rectangular # window, i.e. at (0,0), (10,0), (10,10), and (0,10) xx <- c(x,0,10,10,0) yy <- c(y,0,0,10,10) dxy3 <- deldir(xx,yy,rw=c(0,10,0,10)) # Plot the triangulation created (but not the tessellation). dxy2 <- deldir(x,y,rw=c(0,10,0,10),plot=TRUE,wl="tr") # Example of collinearity error. ## Not run: dniP <- deldir(niProperties) # Throws an error ## End(Not run) dniP <- deldir(niProperties,eps=1e-8) # No error. # Example of using data stored in a data frame. dsw <- deldir(seaweed) # Example where "data" is of class "ppp". dtoy <- deldir(toyPattern) # The "tags", in dtoy$summary$z, are the marks of the toy ppp # object which consists of the letters "a","b","c" and "d". # Artificial example of tags. set.seed(42) trees1 <- sample(c("spruce","birch","poplar","shoe"),20,TRUE) trees2 <- sample(c("fir","maple","larch","palm"),20,TRUE) egDat <- data.frame(x=runif(20),y=runif(20),species=trees1) tagNm <- "species" species <- trees2 dd1 <- deldir(egDat) # No tags. dd2 <- deldir(egDat,z=species) # Uses trees1 as the tags. dd3 <- deldir(egDat,z="species") # Same as dd2. dd4 <- deldir(egDat,z=tagNm) # Same as dd2 and dd3. spec <- species dd5 <- deldir(egDat,z=spec) # Uses trees2 as the tags. # Duncan Murdoch's examples. The deldir() function was not # handling such examples correctly until Duncan kindly set # me on the right path. set.seed(123) dd6 <- deldir(rnorm(32),rnorm(32),rnorm(32)) # x <- cbind(x = 1:10, junk = 11:20) y <- 21:30 z <- 31:40 d7 <- deldir(x=x, y=y, z=z) # # print(d7$summary) reveals that x is 1:10, y is 21:30 # and z is 31:40; x[,"junk"] is ignored as it should be. x <- cbind(x = 1:10, "rnorm(10)" = 11:20) y <- 21:30 z <- 41:50 d8 <- deldir(x=x, y=y, z=rnorm(10)) # # print(d8$summary) reveals that x is 1:10, y is 21:30 and z is a # vector of standard normal values. Even though x has a column with # the name of the z argument i.e. "rnorm(10)" (!!!) the specified # value of z takes precedence over this column (and, as per the usual # R syntax) over the object named "z" in the global environment. # Artificial example of the use of the "id" argument. set.seed(42) x <- runif(50) y <- runif(50) ll <- expand.grid(a=letters[1:10],b=letters[1:10]) aa <- sample(paste0(ll[["a"]],ll[["b"]]),50) dxy.wid <- deldir(x,y,id=aa)x <- c(2.3,3.0,7.0,1.0,3.0,8.0) y <- c(2.3,3.0,2.0,5.0,8.0,9.0) # Let deldir() choose the rectangular window. dxy1 <- deldir(x,y) # User chooses the rectangular window. dxy2 <- deldir(x,y,rw=c(0,10,0,10)) # Put "dummy" points at the corners of the rectangular # window, i.e. at (0,0), (10,0), (10,10), and (0,10) xx <- c(x,0,10,10,0) yy <- c(y,0,0,10,10) dxy3 <- deldir(xx,yy,rw=c(0,10,0,10)) # Plot the triangulation created (but not the tessellation). dxy2 <- deldir(x,y,rw=c(0,10,0,10),plot=TRUE,wl="tr") # Example of collinearity error. ## Not run: dniP <- deldir(niProperties) # Throws an error ## End(Not run) dniP <- deldir(niProperties,eps=1e-8) # No error. # Example of using data stored in a data frame. dsw <- deldir(seaweed) # Example where "data" is of class "ppp". dtoy <- deldir(toyPattern) # The "tags", in dtoy$summary$z, are the marks of the toy ppp # object which consists of the letters "a","b","c" and "d". # Artificial example of tags. set.seed(42) trees1 <- sample(c("spruce","birch","poplar","shoe"),20,TRUE) trees2 <- sample(c("fir","maple","larch","palm"),20,TRUE) egDat <- data.frame(x=runif(20),y=runif(20),species=trees1) tagNm <- "species" species <- trees2 dd1 <- deldir(egDat) # No tags. dd2 <- deldir(egDat,z=species) # Uses trees1 as the tags. dd3 <- deldir(egDat,z="species") # Same as dd2. dd4 <- deldir(egDat,z=tagNm) # Same as dd2 and dd3. spec <- species dd5 <- deldir(egDat,z=spec) # Uses trees2 as the tags. # Duncan Murdoch's examples. The deldir() function was not # handling such examples correctly until Duncan kindly set # me on the right path. set.seed(123) dd6 <- deldir(rnorm(32),rnorm(32),rnorm(32)) # x <- cbind(x = 1:10, junk = 11:20) y <- 21:30 z <- 31:40 d7 <- deldir(x=x, y=y, z=z) # # print(d7$summary) reveals that x is 1:10, y is 21:30 # and z is 31:40; x[,"junk"] is ignored as it should be. x <- cbind(x = 1:10, "rnorm(10)" = 11:20) y <- 21:30 z <- 41:50 d8 <- deldir(x=x, y=y, z=rnorm(10)) # # print(d8$summary) reveals that x is 1:10, y is 21:30 and z is a # vector of standard normal values. Even though x has a column with # the name of the z argument i.e. "rnorm(10)" (!!!) the specified # value of z takes precedence over this column (and, as per the usual # R syntax) over the object named "z" in the global environment. # Artificial example of the use of the "id" argument. set.seed(42) x <- runif(50) y <- runif(50) ll <- expand.grid(a=letters[1:10],b=letters[1:10]) aa <- sample(paste0(ll[["a"]],ll[["b"]]),50) dxy.wid <- deldir(x,y,id=aa)
Create the “dividing chain” of a Dirichlet tessellation. The tessellation must have been created from a set of points having associated “tags”. The dividing chain consists of those edges of Dirichlet tiles which separate points having different values of the given tags.
divchain(x, ...) ## Default S3 method: divchain(x, y, z, ...) ## S3 method for class 'deldir' divchain(x, ...)divchain(x, ...) ## Default S3 method: divchain(x, y, z, ...) ## S3 method for class 'deldir' divchain(x, ...)
x |
Either an object specifying coordinates (in the case of the
|
y |
A numeric vector constituting the |
z |
A vector or factor specifying “auxiliary” values or
“tags”. If this argument is left |
... |
Arguments to be passed to |
An object of class “divchain” consisting of a data frame with columns named “x0”, “y0”, “x1”, “y1”, “v01”, “v02”, “v03”, “v11”, “v12” and “v13”.
The columns named “x0” and “y0” consist of the coordinates of one endpoint of an edge of a Dirichlet tile and the columns named “x1” and “y1” consist of the coordinates of the other endpoint.
The columns named “vij”, i = 0, 1, j = 1, 2, 3, consist
of the indices of the vertices of the Delaunay triangles
whose circumcentres constitute the respective endpoints of the
corresponding edge of a Dirichlet tile. The entries of column
“vi3” may (also) take the values $-1, -2, -3$, and $-4$.
This will be the case if the circumcentre in question lay outside
of the rectangular window rw (see deldir())
enclosing the points being tessellated. In these circumstances the
corresponding endpoint of the tile edge is the intersection of the
line joining the two circumcentres with the boundary of rw,
and the numeric value of the entry of column “vi3” indicates
which side. The numbering follows the convention for numbering
the sides of a plot region in R: 1 for the bottom side,
2 for the left side, 3 for the top side and 4 for the right side.
Note that the triple of vertices uniquely identify the endpoint of the tile edge.
The object has an attribute rw which is equal to
the specification of the rectangular window within which
the triangulation/tessellation in question was constructed.
(See deldir().)
This function was created in response to a question asked
on stackoverflow.com by a user named “Dan”.
Rolf Turner [email protected]
deldir() plot.divchain()
set.seed(42) x <- runif(50) y <- runif(50) z <- factor(kmeans(cbind(x,y),centers=4)$cluster) dc1 <- divchain(x,y,z,rw=c(0,1,0,1)) dxy <- deldir(x,y,z=z,rw=c(0,1,0,1)) dc2 <- divchain(dxy)set.seed(42) x <- runif(50) y <- runif(50) z <- factor(kmeans(cbind(x,y),centers=4)$cluster) dc1 <- divchain(x,y,z,rw=c(0,1,0,1)) dxy <- deldir(x,y,z=z,rw=c(0,1,0,1)) dc2 <- divchain(dxy)
Find which points among a given set are duplicates of others.
duplicatedxy(x, y)duplicatedxy(x, y)
x |
Either a vector of |
y |
A vector of |
Often it is of interest to associate each Dirichlet tile in a
tessellation of a planar point set with the point determining
the tile. This becomes problematic if there are duplicate
points in the set being tessellated/triangulated. Duplicated
points are automatically eliminated “internally” by
deldir(). The association between tiles and the indices
of the original set of points is now preserved by the component
ind.orig of the object returned by deldir().
However confusion could still arise.
If it is of interest to associate Dirichlet tiles with the
points determining them, then it is better to proceed by
eliminating duplicate points to start with. This function
(duplicatedxy()) provides a convenient way of doing so.
A logical vector of length equal to the (original) number
of points being considered, with entries TRUE if the
corresponding point is a duplicate of a point with a smaller
index, and FALSE otherwise.
Which indices will be considered to be indices of duplicated
points (i.e. get TRUE values) will of course depend on
the order in which the points are presented.
The real work is done by the base R function duplicated().
Rolf Turner [email protected]
duplicated(), deldir()
set.seed(42) xy <- data.frame(x=runif(20),y=runif(20)) # Lots of duplicated points. xy <- rbind(xy,xy[sample(1:20,20,TRUE),]) # Scramble. ii <- sample(1:40,40) x <- xy$x[ii] y <- xy$y[ii] # Unduplicate! iii <- !duplicatedxy(x,y) xu <- x[iii] yu <- y[iii] # The i-th tile is determined by (xu[i],yu[i]): dxy <- deldir(xu,yu)set.seed(42) xy <- data.frame(x=runif(20),y=runif(20)) # Lots of duplicated points. xy <- rbind(xy,xy[sample(1:20,20,TRUE),]) # Scramble. ii <- sample(1:40,40) x <- xy$x[ii] y <- xy$y[ii] # Unduplicate! iii <- !duplicatedxy(x,y) xu <- x[iii] yu <- y[iii] # The i-th tile is determined by (xu[i],yu[i]): dxy <- deldir(xu,yu)
Lists the indices (or identifiers if these are provided) of the Delaunay neighbours of each point in the set of points being triangulated/tessellated.
getNbrs(object, interior = NULL)getNbrs(object, interior = NULL)
object |
An object of class |
interior |
Either a rectangle, i.e. a numeric vector of length 4,
If |
If interior is specified then Delaunay neighbours are listed
only for those points which lie in interior. Note that
these neighbours need not themselves lie in interior.
Note also that it is possible for points i and j to
be neighbours even though the “clipped” versions of the tiles
are discontiguous.
A (named) list of length equal to the number points in the set
being triangulated/tessellated, or to the number of such points that
lie in interior if that argument was specified. The names
of the list are the identifiers of the points as specified by
id if id was specified in the call to deldir()
that produced object. If id was not specified, then
the names are the indices of the points, coerced to character mode.
The entries of this list are vectors of the identifiers or indices of the Delaunay neighbours of the point corresponding to that entry.
Be careful about addressing the entries of the list returned
by this function. If id was not specified in the call
to deldir() that produced object, then the names
of this list are the point indices coerced to character mode.
If interior was specified then the name of -th entry
of the list will not in general be i. E.g. given that
xxx is the list returned by this function, xxx[[14]]
will not in general give the Delaunay neighbours of point 14.
Instead, specify xxx[["14"]] or xxx[[id[14]]] where
id is the vector of identifiers supplied in the call to
deldir().
Rolf Turner [email protected]
See deldir() for references.
set.seed(42) x <- runif(60) y <- runif(60) dxy <- deldir(x,y,rw=c(0,1,0,1)) nbrs <- getNbrs(dxy,interior=c(0.2,0.8,0.2,0.8)) nbrs[["14"]] # Correct. nbrs[[14]] # Incorrect. names(nbrs)[14] # The result is 42. # Thus nbrs[[14]] actually gives the Delaunay neighbours of point 42. # Demonstrate that neighbours can have discontiguous clipped tiles. if(require(polyclip)) { x <- c(0.38,0.44,0.04,0.97,0.43,0.96,0.89,0.64,0.97,0.62,0.33,0.35, 0.40,0.78,0.04,0.75,0.68,0.17,0.26,0.51) y <- c(0.68,0.98,0.76,0.57,0.85,0.19,0.27,0.83,0.69,0.24,0.04,0.14, 0.22,0.48,0.20,0.72,0.01,0.38,0.51,0.00) CP <- list(x=c(0.72,0.93,0.76,0.61,-0.03,-0.04,0.41), y=c(0.46,0.76,0.94,1.03,1.01,0.37,0.31)) dxy <- deldir(x,y,rw=c(0,1,0,1)) TL <- tile.list(dxy) plot(TL,labelPts=TRUE) plot(TL[16],clipp=CP,fillcol="orange",labelPts=TRUE,add=TRUE) polygon(CP,border="red") nbrs <- getNbrs(dxy,interior=CP) # Tiles are clipped to CP. # Note that point 14 is a neighbour of point 16, even though their # clipped tiles do not meet. }set.seed(42) x <- runif(60) y <- runif(60) dxy <- deldir(x,y,rw=c(0,1,0,1)) nbrs <- getNbrs(dxy,interior=c(0.2,0.8,0.2,0.8)) nbrs[["14"]] # Correct. nbrs[[14]] # Incorrect. names(nbrs)[14] # The result is 42. # Thus nbrs[[14]] actually gives the Delaunay neighbours of point 42. # Demonstrate that neighbours can have discontiguous clipped tiles. if(require(polyclip)) { x <- c(0.38,0.44,0.04,0.97,0.43,0.96,0.89,0.64,0.97,0.62,0.33,0.35, 0.40,0.78,0.04,0.75,0.68,0.17,0.26,0.51) y <- c(0.68,0.98,0.76,0.57,0.85,0.19,0.27,0.83,0.69,0.24,0.04,0.14, 0.22,0.48,0.20,0.72,0.01,0.38,0.51,0.00) CP <- list(x=c(0.72,0.93,0.76,0.61,-0.03,-0.04,0.41), y=c(0.46,0.76,0.94,1.03,1.01,0.37,0.31)) dxy <- deldir(x,y,rw=c(0,1,0,1)) TL <- tile.list(dxy) plot(TL,labelPts=TRUE) plot(TL[16],clipp=CP,fillcol="orange",labelPts=TRUE,add=TRUE) polygon(CP,border="red") nbrs <- getNbrs(dxy,interior=CP) # Tiles are clipped to CP. # Note that point 14 is a neighbour of point 16, even though their # clipped tiles do not meet. }
A data set taken from an example in the grapherator package. This data set demonstrates handling a data set with duplicated points.
grapherXmplgrapherXmpl
A data frame with 250 observations on the following 2 variables.
xa numeric vector
ya numeric vector
There are 25 duplicated points, so the net number of
observations is 225. These data constitute a structure (named
coordinates) generated internally in the function
addEdgesDelaunay. The call is to be found in the
examples in the help file for the plot.grapherator()
in the grapherator package. The relevant example
initially threw an error, revealing a bug in deldir()
that was triggered when there were duplicated points in the data.
The grapherator package,
https://CRAN.R-project.org/package=grapherator
dgX <- deldir(grapherXmpl) # Now works!!!`dgX <- deldir(grapherXmpl) # Now works!!!`
Produce a summary of a Dirichlet (Voronoi) tessellation in terms of parameters relevant to Lewis's law and Aboav-Weaire's law. Note that “law” in the function name corresponds to “Lewis-Aboav-Weaire”.
lawSummary(object)lawSummary(object)
object |
An object of class |
Tiles are stripped away from the tessellation in “layers”.
Layer 1 consists of “boundary” tiles, i.e. tiles having
at least one vertex on the enclosing rectangle (determined by the
rw argument of deldir()). Layer 2 consists
of tiles which are neighbours of tiles in layer 1 (i.e. tiles
determined by points that are Delaunay neighbours of points
determining the tiles in layer 1). Layer 3 consists of tiles
which are neighbours of tiles in layer 2.
The parameters of interest in respect of the Lewis-Aboav-Weaire summary are then calculated in terms of the tiles that remain after the three layers have been stripped away, which will be referred to as “interior” tiles. These parameters are:
the areas of each of the interior tiles
the number of edges of each of the interior tiles
the number of edges of all neighbouring tiles of each of the interior tiles.
Note that the neighbouring tiles of the interior tiles may include tiles which are not themselves interior tiles (i.e. tiles which are in layer 3).
This function was created at the request of Kai Xu (Fisheries College, Jimei University, Xiamen, Fujian, China 361021).
If no tiles remain after the three layers have been stripped
away, then the returned value is NULL. Otherwise the
returned value is a list with components calculated in terms of
the remaining (“interior”) tiles. These components are:
tile.vertices A list whose entries are data frames
giving the coordinates of the vertices of the interior tiles.
tile.areas A vector of the areas of the interior tiles
in the tessellation in question.
tile.tags A vector or factor whose values are the “tags”
of the interior tiles. The “original” of this object (the
“tags” associated with all of the tiles) is provided
as the z argument to deldir(). The tile.tags
component of the value returned by lawSummary() is present
only if deldir() was called with a (non-NULL) value
of the z argument.
num.edges A vector of the number of edges of each
such tile.
num.nbr.edges A list with a component for each
point, in the set being tessellated, whose corresponding tile
is an interior tile. Each component of this
list is the vector of the number of edges of the interior tiles
determined by points which are Delaunay neighbours of
the point corresponding to the list component in question.
totnum.nbr.edges A vector whose entries consist
of the sums of the vectors in the foregoing list.
The returned list also has attributes as follows:
i1 An integer vector whose entries are in the indices
of the tiles in layer 1.
i2 An integer vector whose entries are in the indices
of the tiles in layer 2.
i3 An integer vector whose entries are in the indices
of the tiles in layer 3.
i.kept An integer vector whose entries are in the indices
of the tiles that are kept, i.e. those that remain after the three layers
have been stripped away.
Rolf Turner [email protected]
# A random pattern: set.seed(42) xy1 <- data.frame(x=runif(400,0,20),y=runif(400,0,20)) dxy1 <- deldir(xy1) ldxy1 <- lawSummary(dxy1) tl1 <- tile.list(dxy1) plot(0,0,type="n",xlim=c(-2,35),ylim=c(0,20),asp=1,xlab="x",ylab="y",bty="l") plot(tl1,showpoints=FALSE,add=TRUE) points(xy1[attr(ldxy1,"i1"),],pch=20,col="yellow") points(xy1[attr(ldxy1,"i2"),],pch=20,col="blue") points(xy1[attr(ldxy1,"i3"),],pch=20,col="green") points(xy1[attr(ldxy1,"i.kept"),],pch=20,col="red") legend("right",pch=20,col=c("yellow","blue","green","red"), legend=c("layer 1","layer 2","layer 3","interior")) # A highly structured pattern (example due to Kai Xu): set.seed(115) x <- c(rep(1:20,10),rep((1:20)+0.5,10)) y <- c(rep(1:10,each=20),rep((1:10)+0.5,each=20))*sqrt(3) a <- runif(400,0,2*pi) b <- runif(400,-1,1) x <- x+0.1*cos(a)*b y <- y+0.1*sin(a)*b xy2 <- data.frame(x,y) dxy2 <- deldir(xy2) ldxy2 <- lawSummary(dxy2) tl2 <- tile.list(dxy2) plot(0,0,type="n",xlim=c(-2,35),ylim=c(0,20),asp=1,xlab="x",ylab="y",bty="l") plot(tl2,showpoints=FALSE,add=TRUE) points(xy2[attr(ldxy2,"i1"),],pch=20,col="yellow") points(xy2[attr(ldxy2,"i2"),],pch=20,col="blue") points(xy2[attr(ldxy2,"i3"),],pch=20,col="green") points(xy2[attr(ldxy2,"i.kept"),],pch=20,col="red") legend("right",pch=20,col=c("yellow","blue","green","red"), legend=c("layer 1","layer 2","layer 3","interior"))# A random pattern: set.seed(42) xy1 <- data.frame(x=runif(400,0,20),y=runif(400,0,20)) dxy1 <- deldir(xy1) ldxy1 <- lawSummary(dxy1) tl1 <- tile.list(dxy1) plot(0,0,type="n",xlim=c(-2,35),ylim=c(0,20),asp=1,xlab="x",ylab="y",bty="l") plot(tl1,showpoints=FALSE,add=TRUE) points(xy1[attr(ldxy1,"i1"),],pch=20,col="yellow") points(xy1[attr(ldxy1,"i2"),],pch=20,col="blue") points(xy1[attr(ldxy1,"i3"),],pch=20,col="green") points(xy1[attr(ldxy1,"i.kept"),],pch=20,col="red") legend("right",pch=20,col=c("yellow","blue","green","red"), legend=c("layer 1","layer 2","layer 3","interior")) # A highly structured pattern (example due to Kai Xu): set.seed(115) x <- c(rep(1:20,10),rep((1:20)+0.5,10)) y <- c(rep(1:10,each=20),rep((1:10)+0.5,each=20))*sqrt(3) a <- runif(400,0,2*pi) b <- runif(400,-1,1) x <- x+0.1*cos(a)*b y <- y+0.1*sin(a)*b xy2 <- data.frame(x,y) dxy2 <- deldir(xy2) ldxy2 <- lawSummary(dxy2) tl2 <- tile.list(dxy2) plot(0,0,type="n",xlim=c(-2,35),ylim=c(0,20),asp=1,xlab="x",ylab="y",bty="l") plot(tl2,showpoints=FALSE,add=TRUE) points(xy2[attr(ldxy2,"i1"),],pch=20,col="yellow") points(xy2[attr(ldxy2,"i2"),],pch=20,col="blue") points(xy2[attr(ldxy2,"i3"),],pch=20,col="green") points(xy2[attr(ldxy2,"i.kept"),],pch=20,col="red") legend("right",pch=20,col=c("yellow","blue","green","red"), legend=c("layer 1","layer 2","layer 3","interior"))
The locations (in longitude and latitude) of a number of properties (land holdings) in Northern Ireland.
data("niProperties")data("niProperties")
A data frame with 240 observations on the following 2 variables.
xA numeric vector of longitudes.
yA numeric vector of latitudes.
These data were kindly provided by Adam Dadvar of the
Cartesian Limited consulting service.
URL: http://www.cartesian.com.
# data(niProperties) # It is unnecessary to use \code{data} since \code{niProperties} is # a "first class object". It is "lazily loaded". plot(niProperties)# data(niProperties) # It is unnecessary to use \code{data} since \code{niProperties} is # a "first class object". It is "lazily loaded". plot(niProperties)
This is a method for plot.
## S3 method for class 'deldir' plot(x,add=FALSE,wlines=c("both","triang","tess"), showpoints=TRUE,labelPts=FALSE,cex=1,lex=1, cmpnt_col=NULL,cmpnt_lty=NULL,pch=1,xlim=NULL, ylim=NULL,axes=FALSE,xlab=if(axes) "x" else "", ylab=if(axes) "y" else"",showrect=FALSE,asp=1,...)## S3 method for class 'deldir' plot(x,add=FALSE,wlines=c("both","triang","tess"), showpoints=TRUE,labelPts=FALSE,cex=1,lex=1, cmpnt_col=NULL,cmpnt_lty=NULL,pch=1,xlim=NULL, ylim=NULL,axes=FALSE,xlab=if(axes) "x" else "", ylab=if(axes) "y" else"",showrect=FALSE,asp=1,...)
x |
An object of class "deldir" as returned by the function deldir. |
add |
logical argument; should the plot be added to an existing plot? |
wlines |
"which lines?". I.e. should the Delaunay triangulation be plotted (wlines="triang"), should the Dirichlet tessellation be plotted (wlines="tess"), or should both be plotted (wlines="both", the default) ? |
showpoints |
Logical scalar; should the points being triangulated/tessellated be plotted? |
labelPts |
Logical argument, defaulting to |
cex |
The value of the character expansion argument cex to be used with the plotting symbols for plotting the points. |
lex |
The value of the character expansion argument cex to be used by the
text function when labelling the points with their identifiers. Used only
if |
cmpnt_col |
A vector or list specifying the colours to be used for plotting
the (up to five) different components of the graphic, explicitly the
triangulation, the tessellation, the data points,
the point labels and the enclosing rectangle ( Colours may be specified as positive integers, corresponding to
the entries of the current If fewer than five colours are given, the missing components are
filled in with |
cmpnt_lty |
A vector or list (of length two) of line types for plotting
the two components of the graphic that consist of lines,
i.e. the triangulation and the tessellation. The types may
consist of integers between 1 and 6, or of approriate text
strings ("solid", "dashed", "dotted", "dotdash", "longdash" or
"twodash"; see the discussion of |
pch |
The plotting symbol for plotting the points. May be either integer or character. |
xlim |
The limits on the x-axis. Defaults to |
ylim |
The limits on the y-axis. Defaults to |
axes |
Logical scalar. Should axes be drawn on the plot? |
xlab |
Label for the x-axis. Defaults to |
ylab |
Label for the y-axis. Defaults to |
showrect |
Logical scalar; show the enclosing rectangle |
asp |
The aspect ratio of the plot; integer scalar or |
... |
Further plotting parameters (e.g. |
A plot of the edges of the Delaunay triangles and/or of the Dirichlet tiles is produced or added to an existing plot. By default the triangles are plotted with solid lines (lty=1) and the tiles with dotted lines (lty=2). In addition the points being triangulated may be plotted.
As of release 1.0-8 the argument number of
plot.deldir() and plot.tile.list()
was changed to labelPts. As a consequence the argument
nex of this function has had its name changed to
lex.
In previous versions of the deldir package, the aspect
ratio was not set. Instead, the shape of the plotting region was
set to be square (pty="s"). This practice was suboptimal.
To reproduce previous behaviour set asp=NA (and possibly
pty="s") in the call to plot.deldir(). Users, unless
they really understand what they are doing and why they are
doing it, are now strongly advised not to set the value of
asp but rather to leave asp equal to its default value
of 1. Any other value may distort the tessellation and destroy
the perpendicular appearance of lines which are indeed perpendicular.
(And conversely may cause lines which are not perpendicular to
appear as if they are.)
Rolf Turner [email protected]
deldir()
x <- c(2.3,3.0,7.0,1.0,3.0,8.0) + 0.5 y <- c(2.3,3.0,2.0,5.0,8.0,9.0) + 0.5 x <- c(x,1,10,10,1) y <- c(y,1,1,10,10) dxy <- deldir(x,y,rw=c(0,11,0,11)) plot(dxy) # Plots the tessellation, but does not save the results: deldir(x,y,rw=c(0,11,0,11),plot=TRUE, wl="tess",cmpnt_col=c(1,1,2,3,4),labelPts=TRUE) # Plots the triangulation, but not the tessellation or the points, # and saves the returned structure: dxy <- deldir(x,y,rw=c(0,11,0,11),plot=TRUE,wl="triang",wp="none") # Plot everything: plot(dxy,cmpnt_col=c("orange","green","red","blue","black"),cmpnt_lty=1, labelPts=TRUE,lex=1.5,pch=c(19,9),showrect=TRUE,axes=TRUE) # Complicated example from He Huang: # "Linguistic distances". vldm <- c(308.298557,592.555483,284.256926,141.421356,449.719913, 733.976839,591.141269,282.842712,1.414214,732.562625) ldm <- matrix(nrow=5,ncol=5) ldm[row(ldm) > col(ldm)] <- vldm ldm[row(ldm) <= col(ldm)] <- 0 ldm <- (ldm + t(ldm))/2 rownames(ldm) <- LETTERS[1:5] colnames(ldm) <- LETTERS[1:5] # Data to be triangulated. id <- c("A","B","C","D","E") x <- c(0.5,1,1,1.5,2) y <- c(5.5,3,7,6.5,5) dat_Huang <- data.frame(id = id, x = x, y = y) # Form the triangulation/tessellation. dH <- deldir(dat_Huang) # Plot the triangulation with line widths proportional # to "linguistic distances". all_col <- rainbow(100) i <- pmax(1,round(100*vldm/max(vldm))) use_col <- all_col[i] ind <- as.matrix(dH$delsgs[,c("ind1","ind2")]) lwv <- ldm[ind] lwv <- 10*lwv/max(lwv) plot(dH,wlines="triang",col=use_col,showpoints=FALSE, lw=lwv,asp=NA) with(dH,text(x,y,id,cex=1.5))x <- c(2.3,3.0,7.0,1.0,3.0,8.0) + 0.5 y <- c(2.3,3.0,2.0,5.0,8.0,9.0) + 0.5 x <- c(x,1,10,10,1) y <- c(y,1,1,10,10) dxy <- deldir(x,y,rw=c(0,11,0,11)) plot(dxy) # Plots the tessellation, but does not save the results: deldir(x,y,rw=c(0,11,0,11),plot=TRUE, wl="tess",cmpnt_col=c(1,1,2,3,4),labelPts=TRUE) # Plots the triangulation, but not the tessellation or the points, # and saves the returned structure: dxy <- deldir(x,y,rw=c(0,11,0,11),plot=TRUE,wl="triang",wp="none") # Plot everything: plot(dxy,cmpnt_col=c("orange","green","red","blue","black"),cmpnt_lty=1, labelPts=TRUE,lex=1.5,pch=c(19,9),showrect=TRUE,axes=TRUE) # Complicated example from He Huang: # "Linguistic distances". vldm <- c(308.298557,592.555483,284.256926,141.421356,449.719913, 733.976839,591.141269,282.842712,1.414214,732.562625) ldm <- matrix(nrow=5,ncol=5) ldm[row(ldm) > col(ldm)] <- vldm ldm[row(ldm) <= col(ldm)] <- 0 ldm <- (ldm + t(ldm))/2 rownames(ldm) <- LETTERS[1:5] colnames(ldm) <- LETTERS[1:5] # Data to be triangulated. id <- c("A","B","C","D","E") x <- c(0.5,1,1,1.5,2) y <- c(5.5,3,7,6.5,5) dat_Huang <- data.frame(id = id, x = x, y = y) # Form the triangulation/tessellation. dH <- deldir(dat_Huang) # Plot the triangulation with line widths proportional # to "linguistic distances". all_col <- rainbow(100) i <- pmax(1,round(100*vldm/max(vldm))) use_col <- all_col[i] ind <- as.matrix(dH$delsgs[,c("ind1","ind2")]) lwv <- ldm[ind] lwv <- 10*lwv/max(lwv) plot(dH,wlines="triang",col=use_col,showpoints=FALSE, lw=lwv,asp=NA) with(dH,text(x,y,id,cex=1.5))
Plot the dividing chain of a Dirichlet tessellation. The tessellation must have been created from a set of points having associated categorical “tags”. The dividing chain consists of those edges of Dirichlet tiles which separate points having different values of the given tags.
## S3 method for class 'divchain' plot(x, add = FALSE, ...)## S3 method for class 'divchain' plot(x, add = FALSE, ...)
x |
An object of class “divchain”. See |
add |
Logical scalar. It |
... |
Graphical parameters such as |
None.
This function was created in response to a question asked
on stackoverflow.com by a user named “Dan”.
Rolf Turner [email protected]
divchain() divchain.default()
divchain.deldir() deldir()
set.seed(42) x <- runif(50) y <- runif(50) z <- factor(kmeans(cbind(x,y),centers=4)$cluster) dc <- divchain(x,y,z,rw=c(0,1,0,1)) plot(dc,lwd=2,col="blue",bty="o")set.seed(42) x <- runif(50) y <- runif(50) z <- factor(kmeans(cbind(x,y),centers=4)$cluster) dc <- divchain(x,y,z,rw=c(0,1,0,1)) plot(dc,lwd=2,col="blue",bty="o")
A method for plot. Plots (sequentially)
the tiles associated with each point in the set being tessellated.
## S3 method for class 'tile.list' plot(x, verbose = FALSE, close = FALSE, pch = 1, fillcol = getCol(x,warn=warn), col.pts=NULL, col.lbls=NULL,border=NULL, showpoints = !labelPts, add = FALSE, asp = 1, clipp=NULL, xlab = "x", ylab = "y", main = "", axes=TRUE,warn=TRUE, labelPts=FALSE,adj=NULL,...)## S3 method for class 'tile.list' plot(x, verbose = FALSE, close = FALSE, pch = 1, fillcol = getCol(x,warn=warn), col.pts=NULL, col.lbls=NULL,border=NULL, showpoints = !labelPts, add = FALSE, asp = 1, clipp=NULL, xlab = "x", ylab = "y", main = "", axes=TRUE,warn=TRUE, labelPts=FALSE,adj=NULL,...)
x |
A list of the tiles in a tessellation, as produced
the function |
verbose |
Logical scalar; if |
close |
Logical scalar; if |
pch |
The plotting character (or vector of plotting
characters) with which to plot the points of the pattern which
was tessellated. Ignored if |
fillcol |
Optional vector (possibly of length 1, i.e. a scalar) whose
entries can be interpreted as colours by |
col.pts |
Optional vector like unto |
col.lbls |
Optional vector like unto |
border |
A scalar that can be interpreted as a colour by |
showpoints |
Logical scalar; if |
add |
Logical scalar; should the plot of the tiles be added to an existing plot? |
asp |
The aspect ratio of the plot; integer scalar or
|
clipp |
An object specifying a polygon to which the tessellation being plotted should be clipped. It should consist either of:
If this argument is provided then the plot of the tessellation
is “clipped” to the polygon specified by |
xlab |
Label for the |
ylab |
Label for the |
main |
A title for the plot (used only if |
axes |
Logical scalar. Should axes (with a marked numeric scale)
be plotted? This argument is used only if |
warn |
Logical scalar passed to the internal function |
labelPts |
Logical scalar; if |
adj |
The “adjustment” argument to |
... |
Optional arguments; may be passed to |
The list of tiles being plotted. This will be the input
list of tiles specified by argument x, or this list
clipped to the polygon clipp if the latter was specified.
As of release 1.0-8 the argument number of
plot.deldir() and plot.tile.list()
was changed to labelPts. As a consequence the argument
col.num of this function has had its name changed to
col.lbls.
The behaviour of this function with respect to
“clipping” has changed substantially since the previous
release of deldir, i.e. 1.1-0. The argument clipwin
has been re-named clipp (“p” for “polygon”).
Clipping is now effected via the new package polyclip.
The spatstat package is no longer used. The argument
use.gpclib has been eliminated, since gpclib (which
used to be called upon by spatstat has been superseded by
polyclip which has an unrestrictive license.
As of release 0.1-1 of the deldir package, the
argument fillcol to this function replaces the old
argument polycol, but behaves somewhat differently.
The argument showrect which was present in versions
of this function prior to release 0.1-1 has been eliminated.
It was redundant.
As of release 0.1-1 the col.pts argument might
behave somewhat differently from how it behaved in the past.
The arguments border, clipp, and warn
are new as of release 0.1-1.
Users, unless they really understand what they are
doing and why they are doing it, are strongly advised
not to set the value of asp but rather to leave asp
equal to its default value of 1. Any other value distorts
the tesselation and destroys the perpendicular appearance of
lines which are indeed perpendicular. (And conversely can cause
lines which are not perpendicular to appear as if they are.)
The argument asp was added at the request of Zubin Dowlaty
(who presumably knows what he's doing!).
If clipp is not NULL and showpoints
is TRUE then it is possible that some of the points
“shown” will not fall inside any of the plotted tiles.
(This will happen if the parts of the tiles in which they fall
have been “clipped” out.) If a tile is clipped out
completely then the point which determines that tile is
not plotted irrespective of the value of showpoints.
If the z components of the entries of x
exist but cannot all be interpreted as colours then the internal
function getCol() returns NA. If warn is
TRUE then a warning is issued. The function getCol()
will also return NA (no warning is issued in this case)
if there aren't any z components. This circumstance
will arise if no z argument was supplied in the call to
deldir(). An NA value of fillcol results
(as is indicated by the argument list entry for fillcol)
in (all of) the tiles being left unfilled.
The change from argument polycol to argument
fillcol, and the resulting change in the way in which
plotted tiles are filled with colours, was made as a result of
a request from Chris Triggs. Likewise the argument clipp
was added due to a request from Chris Triggs.
Rolf Turner [email protected]
deldir(), tile.list(),
triang.list(), plot.triang.list()
set.seed(42) x <- runif(20) y <- runif(20) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) plot(w) ccc <- heat.colors(20) # Or topo.colors(20), or terrain.colors(20) # or cm.colors(20), or rainbow(20). plot(w,fillcol=ccc,close=TRUE) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) cul <- rainbow(10)[c(1,7,3:6,2,8:10)] # Rearranging colours to improve # the contrast between contiguous tiles. plot(w,clipp=CP,showpoints=FALSE,fillcol=cul) } plot(w,labelPts=TRUE,col.lbls="red") plot(w,labelPts=TRUE,col.lbls="red",cex=0.5) plot(w,showpoints=TRUE,labelPts=TRUE,col.pts="green",col.lbls="red") plot(w,axes=FALSE,xlab="",ylab="",close=TRUE)set.seed(42) x <- runif(20) y <- runif(20) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) plot(w) ccc <- heat.colors(20) # Or topo.colors(20), or terrain.colors(20) # or cm.colors(20), or rainbow(20). plot(w,fillcol=ccc,close=TRUE) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) cul <- rainbow(10)[c(1,7,3:6,2,8:10)] # Rearranging colours to improve # the contrast between contiguous tiles. plot(w,clipp=CP,showpoints=FALSE,fillcol=cul) } plot(w,labelPts=TRUE,col.lbls="red") plot(w,labelPts=TRUE,col.lbls="red",cex=0.5) plot(w,showpoints=TRUE,labelPts=TRUE,col.pts="green",col.lbls="red") plot(w,axes=FALSE,xlab="",ylab="",close=TRUE)
A method for plot. Plots the triangles of
a Delaunay triangulation of a set of points in the plane.
## S3 method for class 'triang.list' plot(x, showrect = FALSE, add = FALSE, xlab = "x", ylab = "y", main = "", asp = 1, rectcol="black", ...)## S3 method for class 'triang.list' plot(x, showrect = FALSE, add = FALSE, xlab = "x", ylab = "y", main = "", asp = 1, rectcol="black", ...)
x |
An object of class “triang.list” as produced by
|
showrect |
Logical scalar; show the enclosing rectangle |
add |
Logical scalar; should the plot of the triangles be added to an existing plot? |
xlab |
Label for the |
ylab |
Label for the |
main |
A title for the plot (used only if |
asp |
The aspect ratio of the plot; integer scalar or
|
rectcol |
Text string or integer specifying the colour in which the enclosing
rectangle should be plotted. Ignored unless |
... |
Arguments passed to |
None. This function has the side effect of producing (or adding to) a plot.
Users are strongly advised not to set the value of
asp (unless they really know what they are doing) but rather
to leave asp equal to its default value of 1.
Any other value distorts the tesselation and destroys the
perpendicular appearance of lines which are indeed perpendicular.
(And conversely can cause lines which are not perpendicular to
appear as if they are.)
The argument asp was added at the request of Zubin
Dowlaty (who presumably knows what he is doing!).
Rolf Turner [email protected]
deldir(), plot.triang.list(),
tile.list(), plot.tile.list()
set.seed(42) x <- runif(20) y <- runif(20) d <- deldir(x,y) ttt <- triang.list(d) plot(ttt,border="red",showrect=TRUE,rectcol="green") sss <- tile.list(d) plot(sss) plot(ttt,add=TRUE,border="blue",showrect=TRUE,rectcol="red")set.seed(42) x <- runif(20) y <- runif(20) d <- deldir(x,y) ttt <- triang.list(d) plot(ttt,border="red",showrect=TRUE,rectcol="green") sss <- tile.list(d) plot(sss) plot(ttt,add=TRUE,border="blue",showrect=TRUE,rectcol="red")
Prints a very brief description of an object of class "deldir"
as returned by deldir().
## S3 method for class 'deldir' print(x,digits=NULL,...)## S3 method for class 'deldir' print(x,digits=NULL,...)
x |
A Delaunay triangulation and Dirichlet (Voronoi) tessellation
of a set of points (object of class |
digits |
Integer scalar. The number of digits to which to round the
numeric information before printing. Note this may be give
negative values. (See |
... |
Not used. |
This is a method for the generic print() function.
Rolf Turner [email protected]
print()
set.seed(42) x <- rnorm(200,0,4) y <- rnorm(200,0,4) dxy1 <- deldir(x,y) dxy2 <- deldir(x,y,rw=c(-12,12,-11,11)) dxy1 dxy2 print(dxy1,digits=4)set.seed(42) x <- rnorm(200,0,4) y <- rnorm(200,0,4) dxy1 <- deldir(x,y) dxy2 <- deldir(x,y,rw=c(-12,12,-11,11)) dxy1 dxy2 print(dxy1,digits=4)
Print a reasonably readable summary of an object of class
tileInfo as produced by the tileInfo() function.
## S3 method for class 'tileInfo' print(x, digits = 4, npl = 6, ...)## S3 method for class 'tileInfo' print(x, digits = 4, npl = 6, ...)
x |
An object of class |
digits |
Integer scalar. The (maximum) number of decimal digits to which the output is to be printed. |
npl |
Integer scalar. “Number per line”. It specifies the
(maximum) number of values per line. Used (only) when printing
the edge lengths component of |
... |
Not used. Present for compatibility with the generic
|
The list produced by tileInfo() is a bit messy and
hard to comprehend, especially if there is a large number
of tiles. This print method produces a screen display which
is somewhat more perspicuous.
There are four components to the display:
A list, each entry of which is the vector of edge lengths
of the tile. Each edge length is formatted to have a number
of digits specified by the digits argument. Each list
entry may be displayed over a number of lines. The first of these
lines is prefixed by an “informative” string indicating
the point that determines the tile whose edge lengths are being
printed. The string is formed from the identifier of the point.
See deldir(), plot.deldir() and
getNbrs(). The identifier may consist essentially
of the index of the point in the sequence of points that is
being tessellated.
Succeeding lines, corresponding to the same list entry, are prefixed with a number of blanks so as to produce an aesthetically pleasing alignment.
A table of the edge counts of the tiles.
A simple print out of the areas of the tiles (rounded
to a maximum of digits decimal digits).
A simple print out of the perimeters of the tiles (rounded
to a maximum of digits decimal digits).
This screen display is for “looking at” only. In order
to do further calculations on the output of tileInfo it
is necessary to delve into the bowels of x and extract
the relevant bits.
None.
Rolf Turner [email protected]
tileInfo()
set.seed(179) x <- runif(100) y <- runif(100) dxy <- deldir(x,y,rw=c(0,1,0,1)) ixy1 <- tileInfo(dxy) print(ixy1) ixy2 <- tileInfo(dxy,bndry=TRUE) print(ixy2) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) ixy3 <- tileInfo(dxy,clipp=CP) options(width=120) # And enlarge the console window. print(ixy3) # 33 tiles are retained. print(ixy3$perimeters$perComps) # The tiles for points 9 and 94 have # been split into two components. }set.seed(179) x <- runif(100) y <- runif(100) dxy <- deldir(x,y,rw=c(0,1,0,1)) ixy1 <- tileInfo(dxy) print(ixy1) ixy2 <- tileInfo(dxy,bndry=TRUE) print(ixy2) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) ixy3 <- tileInfo(dxy,clipp=CP) options(width=120) # And enlarge the console window. print(ixy3) # 33 tiles are retained. print(ixy3$perimeters$perComps) # The tiles for points 9 and 94 have # been split into two components. }
A data frame whose columns are the coordinates of the centroids of the cells in a seaweed frond. The points are estimates-by-eye of where the centroids of the cells occur.
data("seaweed")data("seaweed")
A data frame with 266 observations on the following 2 variables.
xThe -coordinates of the cell centroids.
yThe -coordinates of the cell centroids.
These data were kindly supplied by Dr. John Bothwell of the Department of Biosciences, Durham University. The data were collected by Kevin Yun and Georgia Campbell, members of Dr. Bothwell's research group.
# data(seaweed) # It is unnecessary to use \code{data} since \code{seaweed} is # a "first class object". It is "lazily loaded". dsw <- deldir(seaweed) isw <- tileInfo(dsw) # Expand the width of the terminal window. options(width=120) isw tsw <- tile.list(dsw) plot(tsw,labelPts=TRUE,col.lbls="red",cex=0.5,adj=0.5)# data(seaweed) # It is unnecessary to use \code{data} since \code{seaweed} is # a "first class object". It is "lazily loaded". dsw <- deldir(seaweed) isw <- tileInfo(dsw) # Expand the width of the terminal window. options(width=120) isw tsw <- tile.list(dsw) plot(tsw,labelPts=TRUE,col.lbls="red",cex=0.5,adj=0.5)
Given a list of Dirichlet tiles, as produced by tile.list(),
produces a data frame consisting of the centroids of those tiles.
tile.centroids(tl)tile.centroids(tl)
tl |
A list of the tiles (produced by |
A data frame with two columns named x and y.
Each row of this data frame constitutes the centroid of one
of the Dirichlet tiles.
Rolf Turner [email protected]
URL http://en.wikipedia.org/wiki/Centroid
set.seed(42) x <- runif(20) y <- runif(20) d <- deldir(x,y) l <- tile.list(d) g <- tile.centroids(l) plot(l,close=TRUE) points(g,pch=20,col="red")set.seed(42) x <- runif(20) y <- runif(20) d <- deldir(x,y) l <- tile.list(d) g <- tile.centroids(l) plot(l,close=TRUE) points(g,pch=20,col="red")
For each point in the set being tessellated produces a list entry describing the Dirichlet/Voronoi tile containing that point.
tile.list(object,minEdgeLength=NULL,clipp=NULL)tile.list(object,minEdgeLength=NULL,clipp=NULL)
object |
An object of class |
minEdgeLength |
Positive numeric scalar specifying the minimum length that
an edge of a tile may have. It is used to eliminate edges
that are effectively of zero length, which can cause tiles
to be “invalid”. This argument defaults to
|
clipp |
An object specifying a polygon to which the tessellation, whose tiles are being determined, should be clipped. It should consist either of:
If this argument is provided then the tiles in the list that
is produced are “clipped” to the polygon specified by
|
A list with one entry for each of the points in the set being
tessellated, or for each of the tiles that are retained after clipping
if clipp is not NULL. Each entry is in turn a list
with a number of components. These components always include:
ptNum |
The index of the point in the original sequence of points
that is being tessellated. Note that if a point is one of a set
of duplicated points then |
pt |
The coordinates of the point whose tile is being described. |
area |
The area of the tile. |
If the tile in question has not been subdivided by the clipping process then the list components also include:
x |
The |
y |
The |
bp |
Vector of logicals indicating whether the tile vertex is a “real” vertex, or a boundary point, i.e. a point where the tile edge intersects the boundary of the enclosing rectangle. |
If the tile in question has been subdivided then the list
does not have the foregoing three components but rather has a
component tileParts which is in turn a list of length equal
to the number of parts into which the tile was subdivided. Each
component of tileParts is yet another list with four
components x, y, bp and area as described above
and as are appropriate for the part in question.
z |
The “auxiliary value” or “tag” associated
with the |
The author expresses sincere thanks to Majid Yazdani who found and
pointed out a serious bug in tile.list in a previous version
(0.0-5) of the deldir package.
The set of vertices of each tile may be “incomplete”. Only vertices which lie within the enclosing rectangle, and “boundary points” are listed.
Note that the enclosing rectangle may be specified by the user
in the call to deldir().
In contrast to some earlier versions of deldir, the corners
of the enclosing rectangle are now included as vertices of tiles.
I.e. a tile which in fact extends beyond the rectangular window
and contains a corner of that window will have that corner added
to its list of vertices. Thus the corresponding polygon is the
intersection of the tile with the enclosing rectangle.
Rolf Turner [email protected]
deldir(), plot.tile.list()
triang.list() plot.triang.list()
set.seed(42) x <- runif(20) y <- runif(20) z <- deldir(x,y) w <- tile.list(z) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) z <- deldir(x,y,rw=c(0,1,0,1),dpl=list(ndx=2,ndy=2)) w <- tile.list(z) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) wc <- tile.list(z,clipp=CP) # 10 tiles are retained; the third tile, # corresponding to point 6, is # subdivided into two parts. # Determine the tiles on the border of a clipping region. # Example due to Huimin Wang. set.seed(112) x <- runif(100) y <- runif(100) dxy <- deldir(x,y) txy <- tile.list(dxy) chind <- chull(x,y) bdry <- list(x=x[chind],y=y[chind]) ctxy <- tile.list(dxy,clipp=bdry) inside <- lapply(ctxy,function(tile,bdry){insidePoly(tile$x,tile$y,bdry)}, bdry=bdry) border <- sapply(inside,function(x){any(!x) | any(attr(x,"on.boundary"))}) plot(ctxy[border],main="Border tiles") }set.seed(42) x <- runif(20) y <- runif(20) z <- deldir(x,y) w <- tile.list(z) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) z <- deldir(x,y,rw=c(0,1,0,1),dpl=list(ndx=2,ndy=2)) w <- tile.list(z) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) wc <- tile.list(z,clipp=CP) # 10 tiles are retained; the third tile, # corresponding to point 6, is # subdivided into two parts. # Determine the tiles on the border of a clipping region. # Example due to Huimin Wang. set.seed(112) x <- runif(100) y <- runif(100) dxy <- deldir(x,y) txy <- tile.list(dxy) chind <- chull(x,y) bdry <- list(x=x[chind],y=y[chind]) ctxy <- tile.list(dxy,clipp=bdry) inside <- lapply(ctxy,function(tile,bdry){insidePoly(tile$x,tile$y,bdry)}, bdry=bdry) border <- sapply(inside,function(x){any(!x) | any(attr(x,"on.boundary"))}) plot(ctxy[border],main="Border tiles") }
Calculates the area of a Dirichlet tile, applying a discrete version of Stoke's theorem.
tileArea(x, y, rw)tileArea(x, y, rw)
x |
The |
y |
The |
rw |
A vector of length 4 specifying the rectangular window in
which the relevant tessellation was constructed. See
|
The heavy lifting is done by the Fortran subroutine stoke()
which is called by the .Fortran() function.
A positive scalar.
Rolf Turner [email protected]
set.seed(42) x <- runif(20) y <- runif(20) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) with(w[[1]],tileArea(x,y,rw=z$rw)) sapply(w,function(x,rw){tileArea(x$x,x$y,attr(w,"rw"))}) x <- c(0.613102,0.429294,0.386023,0.271880,0.387249,0.455900,0.486101) y <- c(0.531978,0.609665,0.597780,0.421738,0.270596,0.262953,0.271532) # The vertices of the Dirichlet tile for point 6. tileArea(x,y,rw=c(0,1,0,1)) tileArea(x,y,rw=c(-1,2,-3,4)) # Same as above.set.seed(42) x <- runif(20) y <- runif(20) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) with(w[[1]],tileArea(x,y,rw=z$rw)) sapply(w,function(x,rw){tileArea(x$x,x$y,attr(w,"rw"))}) x <- c(0.613102,0.429294,0.386023,0.271880,0.387249,0.455900,0.486101) y <- c(0.531978,0.609665,0.597780,0.421738,0.270596,0.262953,0.271532) # The vertices of the Dirichlet tile for point 6. tileArea(x,y,rw=c(0,1,0,1)) tileArea(x,y,rw=c(-1,2,-3,4)) # Same as above.
Produces a summary of information about the tiles in an
object of class deldir as produced by the function
deldir().
tileInfo(object, bndry = FALSE, clipp=NULL)tileInfo(object, bndry = FALSE, clipp=NULL)
object |
An object of class |
bndry |
Logical scalar. If |
clipp |
An object specifying a polygon to which the tiles of
the tessellation should be clipped. See |
An object of class "tileInfo" which consists of a
list with components:
indivTiles |
This is itself a list. If The entries of
|
allEdgeCounts |
An integer vector of the number of edges for each of the tiles. |
tabEdgeCounts |
A table of |
allEdgeLengths |
A vector of all of the tile edge lengths;
a catenation of the |
Areas |
A vector of the areas of the tiles. |
uniqueEdgeLengths |
A vector of the lengths of the tiles edges
with the duplicates (which occur in |
perimeters |
A list, as returned by |
There is a print() method for class "tileInfo" which
produces a convenient display of the information returned by this
function.
Thanks to Krisztina Konya of Ruhr-University Bochum, who provided an example illustrating the need for an error trap in the setting in which all tiles are boundary tiles.
Rolf Turner [email protected]
deldir() tile.list() print.tileInfo()
tilePerim()
set.seed(42) x <- runif(20) y <- runif(20) dxy <- deldir(x,y,rw=c(0,1,0,1)) ixy1 <- tileInfo(dxy) ixy2 <- tileInfo(dxy,bndry=TRUE) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) ixy3 <- tileInfo(dxy,clipp=CP) # 10 tiles are retained; the third tile, # corresponding to point 6, is # subdivided into two parts. }set.seed(42) x <- runif(20) y <- runif(20) dxy <- deldir(x,y,rw=c(0,1,0,1)) ixy1 <- tileInfo(dxy) ixy2 <- tileInfo(dxy,bndry=TRUE) if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) ixy3 <- tileInfo(dxy,clipp=CP) # 10 tiles are retained; the third tile, # corresponding to point 6, is # subdivided into two parts. }
Calculates the perimeters of all of the Dirichlet (Voronoi) tiles in a tessellation of a set of planar points. Also calculates the sum and the mean of these perimeters.
tilePerim(object,inclbdry=TRUE)tilePerim(object,inclbdry=TRUE)
object |
An object of class |
inclbdry |
Logical scalar. Should boundary segments (edges of tiles
at least one of whose endpoints lies on the enclosing
rectangle |
A list with components
perimeters |
A vector consisting of the values of the perimeters of the Dirichlet tiles in the tessellation. |
totalPerim |
The sum of |
meanPerim |
The mean of |
perComps |
A list whose entries are vectors consisting of the “components” of the perimeters of each tile. If/when the tiles are clipped, some tiles may be subdivided by the clipping into discontiguous parts. The components referred to above are the perimeters of this parts. If no subdivision has occurred then the vector in question has a single entry equal to the perimeter of the corresponding tile. If subdivision has occurred then the perimeter of the tile is the sum of the perimeters of the components. |
Function added at the request of Haozhe Zhang.
Rolf Turner [email protected]
tile.list(), plot.tile.list()
x <- runif(20) y <- runif(20) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) p1 <- tilePerim(w) p0 <- tilePerim(w,inclbdry=FALSE) p1$totalPerim - p0$totalPerim # Get 4 = the perimeter of rw. ss <- apply(as.matrix(z$dirsgs[,1:4]),1, function(x){(x[1]-x[3])^2 + (x[2]-x[4])^2}) 2*sum(sqrt(ss)) - p0$totalPerim # Get 0; in tilePerim() each interior # edge is counted twice. if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) wc <- tile.list(z,clipp=CP) p2 <- tilePerim(wc) # Doesn't matter here if inclbdry is TRUE or FALSE. p2$perComps[["pt.6"]] # The tile for point 6 has got subdivided into # two parts, a tetrahedron and a triangle. cul <- rainbow(10)[c(1,7,3:6,2,8:10)] # Rearranging colours to improve # the contrast between contiguous tiles. plot(wc,labelPts=TRUE,fillcol=cul) }x <- runif(20) y <- runif(20) z <- deldir(x,y,rw=c(0,1,0,1)) w <- tile.list(z) p1 <- tilePerim(w) p0 <- tilePerim(w,inclbdry=FALSE) p1$totalPerim - p0$totalPerim # Get 4 = the perimeter of rw. ss <- apply(as.matrix(z$dirsgs[,1:4]),1, function(x){(x[1]-x[3])^2 + (x[2]-x[4])^2}) 2*sum(sqrt(ss)) - p0$totalPerim # Get 0; in tilePerim() each interior # edge is counted twice. if(require(polyclip)) { CP <- list(x=c(0.49,0.35,0.15,0.20,0.35,0.42, 0.43,0.62,0.46,0.63,0.82,0.79), y=c(0.78,0.86,0.79,0.54,0.58,0.70, 0.51,0.46,0.31,0.20,0.37,0.54)) wc <- tile.list(z,clipp=CP) p2 <- tilePerim(wc) # Doesn't matter here if inclbdry is TRUE or FALSE. p2$perComps[["pt.6"]] # The tile for point 6 has got subdivided into # two parts, a tetrahedron and a triangle. cul <- rainbow(10)[c(1,7,3:6,2,8:10)] # Rearranging colours to improve # the contrast between contiguous tiles. plot(wc,labelPts=TRUE,fillcol=cul) }
A simulated object of class "ppp" provided for use
in an example illustrating the application of deldir()
to point pattern objects.
toyPatterntoyPattern
An object of class "ppp" consisting of a simulated marked
point pattern. Entries include
x |
Cartesian |
y |
Cartesian |
marks |
factor with levels "a","b","c","d"
|
Simulated.
dtoy <- deldir(toyPattern) # "Tags" are the marks of the pattern. set.seed(42) dtoy.nt <- deldir(toyPattern,z=round(runif(59),2)) # Tags are numeric.dtoy <- deldir(toyPattern) # "Tags" are the marks of the pattern. set.seed(42) dtoy.nt <- deldir(toyPattern,z=round(runif(59),2)) # Tags are numeric.
From an object of class “deldir” produces a list of the Delaunay triangles in the triangulation of a set of points in the plane.
triang.list(object)triang.list(object)
object |
An object of class “deldir” as produced by |
A list each of whose components is a or
data frame corresponding to one of the
Delaunay triangles specified by “object”. The rows of each
such data frame correspond to the vertices of the corresponding
Delaunay triangle. The columns are:
ptNum (the index of the point in the original sequence
of points that is being triangulated. Note that if a point is
one of a set of duplicated points then ptNum is
the first of the indices of the points in this set.)
x (the -coordinate of the vertex)
y (the -coordinate of the vertex)
z (the “auxiliary value” or “tag”
z associated with the vertex; present only if such values
were supplied in the call to deldir())
The returned value has an attribute “rw” consisting of the enclosing rectangle of the triangulation.
There may not actually be any triangles determined by
object, in which case this function returns an empty
list with an "rw" attribute. See Examples.
The code of this function was taken more-or-less directly from code written by Adrian Baddeley for the “delaunay()” function in the “spatstat” package.
Rolf Turner [email protected]
deldir(), plot.triang.list(),
tile.list(), plot.tile.list()
set.seed(42) x <- runif(20) y <- runif(20) z <- sample(1:100,20) d1 <- deldir(x,y,z=z) t1 <- triang.list(d1) # A "triangulation" with no triangles! d2 <- deldir(x=1:10,y=11:20) plot(d2) t2 <- triang.list(d2) plot(t2,showrect=TRUE,rectcol="blue") # Pretty boring!set.seed(42) x <- runif(20) y <- runif(20) z <- sample(1:100,20) d1 <- deldir(x,y,z=z) t1 <- triang.list(d1) # A "triangulation" with no triangles! d2 <- deldir(x=1:10,y=11:20) plot(d2) t2 <- triang.list(d2) plot(t2,showrect=TRUE,rectcol="blue") # Pretty boring!
Lists the indices of the vertices of each Delaunay triangle in
the triangulation of a planar point set. The indices are listed
(in increasing numeric order) as the rows of an
matrix where is the number of Delaunay triangles in the
triangulation.
triMat(object)triMat(object)
object |
An object of class |
This function was suggested by Robin Hankin of the School of Mathematical and Computing Sciences at Auckland University of Technology.
An matrix where is the number
of Delaunay triangles in the triangulation specified by object.
The row consists of the indices (in the original
list of points being triangulated) of vertices of the
Delaunay triangle. The indices are listed in increasing numeric
order in each row.
Earlier versions of this function (prior to release 0.1-14
of deldir) could sometimes give incorrect results.
This happened if the union of three contiguous Delaunay triangles
happened to constitute another triangle. This latter triangle
would appear in the list of triangles produced by triMat()
but is not itself a Delaunay triangle. The updated version
of triMat() now checks for this possibility and gives
(I think!) correct results.
Many thanks to Jay Call, who pointed out this bug to me.
Rolf Turner [email protected]
deldir()
triang.list()
plot.triang.list()
# These are the data used by Jay Call to illustrate the bug # that appeared in a previous incarnation of triMat. xy <- data.frame( x = c(0.048,0.412,0.174,0.472,0.607,0.565,0.005,0.237,0.810,0.023), y = c(0.512,0.928,0.955,0.739,0.946,0.134,0.468,0.965,0.631,0.782) ) dxy <- deldir(xy) M <- triMat(dxy) plot(dxy,wlines="triang",num=TRUE,axes=FALSE,cmpnt_col=c(1,1,1,1,2,1)) # The triangle with vertices {4,5,8} was listed in the output of # the previous (buggy) version of triMat(). It is NOT a Delaunay # triangle and hence should NOT be listed.# These are the data used by Jay Call to illustrate the bug # that appeared in a previous incarnation of triMat. xy <- data.frame( x = c(0.048,0.412,0.174,0.472,0.607,0.565,0.005,0.237,0.810,0.023), y = c(0.512,0.928,0.955,0.739,0.946,0.134,0.468,0.965,0.631,0.782) ) dxy <- deldir(xy) M <- triMat(dxy) plot(dxy,wlines="triang",num=TRUE,axes=FALSE,cmpnt_col=c(1,1,1,1,2,1)) # The triangle with vertices {4,5,8} was listed in the output of # the previous (buggy) version of triMat(). It is NOT a Delaunay # triangle and hence should NOT be listed.
Example solute plume concentration data set
associated with the GWSDAT (“GroundWater
Spatiotemporal Data Analysis Tool”) project
https://protect-au.mimecast.com/s/demRC91WzLH6qo3TorzN7?domain=gwsdat.net.
The deldir package is used in this project as part of a
numerical routine to estimate plume quantities (mass, average
concentration and centre of mass).
volTriPointsvolTriPoints
A data frame with 232 observations on the following 3 variables.
xThe x-coordinates of the centres of mass
of the plumes.
yThe y-coordinates of the centres of mass
of the plumes.
zThe plume concentrations.
This data set played a critical role in revealing a bug in the
Fortran code underlying the deldir() function.
These data were kindly provided by Wayne W. Jones of the GWSDAT
project. The data set was originally named Vol.Tri.Points;
the name was changed so as to be more consistent with my usual
naming conventions.
Jones, W. R., Spence, M. J., Bowman, A. W., Evers, L. and Molinari, D. A. (2014). A software tool for the spatiotemporal analysis and reporting of groundwater monitoring data. Environmental Modelling & Software 55, pp. 242–249.
dvtp <- deldir(volTriPoints) plot(dvtp)dvtp <- deldir(volTriPoints) plot(dvtp)
Finds the Dirichlet/Voronoi tile, of a tessellation produced
by deldir(), that contains a given point.
which.tile(x, y, tl)which.tile(x, y, tl)
x |
The |
y |
The |
tl |
A tile list, as produced by the function |
Just minimises the distance from the point in question to the points of the pattern determining the tiles.
An integer equal to the index of the tile in which the given point lies.
Rolf Turner [email protected]
set.seed(42) x <- runif(20,0,100) y <- runif(20,0,100) dxy <- deldir(x,y) txy <- tile.list(dxy) i <- which.tile(30,50,txy) # The value of i here is 14. plot(txy,showpoints=FALSE) text(x,y,labels=1:length(txy),col="red") points(30,50,pch=20,col="blue")set.seed(42) x <- runif(20,0,100) y <- runif(20,0,100) dxy <- deldir(x,y) txy <- tile.list(dxy) i <- which.tile(30,50,txy) # The value of i here is 14. plot(txy,showpoints=FALSE) text(x,y,labels=1:length(txy),col="red") points(30,50,pch=20,col="blue")