Delaunay triangulation and Dirichlet tessellation
This function computes the Delaunay triangulation (and hence the
Dirichlet or Voronoi tesselation) 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. A set of dummy
points may be added, in various ways, to the set of data points
being triangulated.
deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, sort=TRUE, plot=FALSE,
round=TRUE,digits=6, z=NULL, zdum=NULL, suppressMsge=FALSE, ...)x,y |
These arguments specify the coordinates of the point set being
triangulated or tessellated. These can be given by two separate
arguments x and y which are vectors or by a single argument x
which is either a data frame or a generic list, possibly one of
class If If there is a column named “z” and if the argument If |
dpl |
A list describing the structure of the dummy points to be added to the data being triangulated. The addition of these dummy points is effected by the auxiliary function dumpts(). The list may have components:
|
rw |
The coordinates of the corners of the rectangular window enclosing the triangulation, in the order (xmin, xmax, ymin, ymax). Any data points (including dummy 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
|
z |
An optional vector of “auxiliary” values or “weights”
associated with the respective points. (NOTE: These
“weights” are values associated with the points and hence
with the tiles of the tessellation produced. They DO NOT
affect the tessellation, i.e. the tessellation produced is the
same as is it would be if there were no weights. The If |
zdum |
Values of |
suppressMsge |
Logical scalar indicating whether a message (alerting the user to
changes from previous versions of |
... |
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 that 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 nineth 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 real (as opposed to dummy) points in the set which was
triangulated, with any duplicate points eliminated. The first n.data
rows of |
n.dum |
the number of dummy points which were added to the set being triangulated,
with any duplicate points (including any which duplicate real points)
eliminated. The last n.dum 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 (as data points or as dummy
points) to |
If ndx >= 2 and ndy >= 2, then the rectangular window IS the convex
hull, and so the values of del.area and dir.area (if the latter is
not NULL) are identical.
If plot=TRUE a plot of the triangulation and/or tessellation
is produced or added to an existing plot.
It is difficult-to-impossible to determine a priori
how much memory needs to be allocated 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 inadquate, 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 and 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
conceiveably 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 the appropriateness of the “increase
eps” strategy, the user is advised to do his or her 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.
A message alerting the user to the foregoing two items is printed
out the first time that deldir() is called with
suppressMsge=FALSE in a given session. In succeeding
calls to deldir() in the same session, no message is printed.
(I.e. the “alerting” message is printed at most once
in any given session.)
The “alerting” message is not produced via the
warning() function, so suppressWarnings() will
not suppress its appearance. To effect such suppression
(necessary only on the first call to deldir() in a
given session) one must set the suppressMsge argument of
deldir equal to TRUE.
If any dummy points are created, and if a vector z, of
“auxiliary” values or “weights” associated with the
points being triangulated, is supplied, then it is up to the user to
supply the corresponding auxiliary values or weights associated with
the dummy points. These values should be supplied as zdum.
If zdum is not supplied then the auxiliary values or weights
associated with the dummy points are all taken to be missing values
(i.e. NA).
Rolf Turner r.turner@auckland.ac.nz
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)
dxy3 <- deldir(x,y,dpl=list(ndx=2,ndy=2),rw=c(0,10,0,10))
# Plot the triangulation created (but not the tesselation).
## Not run:
dxy2 <- deldir(x,y,rw=c(0,10,0,10),plot=TRUE,wl='tr')
## End(Not run)
# Auxiliary values associated with points; 4 dummy points to be
# added so 4 dummy "z-values" provided.
z <- c(1.63,0.79,2.84,1.56,0.22,1.07)
zdum <- rep(42,4)
dxy4 <- deldir(x,y,dpl=list(ndx=2,ndy=2),rw=c(0,10,0,10),z=z,zdum=zdum)
# Example of collinearity error.
## Not run:
dniP <- deldir(niProperties) # Throws an error
## End(Not run)
dniP <- deldir(niProperties,eps=1e-8) # No error.Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.