Adicionada a biblioteca R deldir em pacotes/rlib/win

Edmar Moretti
1 parent 63df1103
Showing 64 changed files with 4563 additions and 0 deletions Show diff stats
pacotes/rlib/win/deldir/CONTENTS
pacotes/rlib/win/deldir/DESCRIPTION
pacotes/rlib/win/deldir/INDEX
pacotes/rlib/win/deldir/MD5
pacotes/rlib/win/deldir/Meta/Rd.rds
pacotes/rlib/win/deldir/Meta/hsearch.rds
pacotes/rlib/win/deldir/Meta/package.rds
pacotes/rlib/win/deldir/R-ex/deldir.R
pacotes/rlib/win/deldir/R-ex/plot.deldir.R
pacotes/rlib/win/deldir/R-ex/plot.tile.list.R
pacotes/rlib/win/deldir/R-ex/tile.list.R
pacotes/rlib/win/deldir/R/deldir
pacotes/rlib/win/deldir/READ_ME
pacotes/rlib/win/deldir/chtml/deldir.chm
pacotes/rlib/win/deldir/err.list
pacotes/rlib/win/deldir/ex.out
pacotes/rlib/win/deldir/help/AnIndex
pacotes/rlib/win/deldir/help/deldir
pacotes/rlib/win/deldir/help/deldir-internal
pacotes/rlib/win/deldir/help/plot.deldir
@@ -0,0 +1,29 @@
+Entry: deldir-internal
+Aliases: dumpts ind.dup mid.in mnnd get.cnrind acw
+Keywords: internal
+Description: Internal deldir functions
+URL: ../../../library/deldir/html/deldir-internal.html
+
+Entry: deldir
+Aliases: deldir
+Keywords: spatial
+Description: Construct the Delaunay triangulation and the Dirichlet (Voronoi) tessellation of a planar point set.
+URL: ../../../library/deldir/html/deldir.html
+
+Entry: plot.deldir
+Aliases: plot.deldir
+Keywords: 
+Description: Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi) tesselation of a planar point set, as constructed by the function deldir.
+URL: ../../../library/deldir/html/plot.deldir.html
+
+Entry: plot.tile.list
+Aliases: plot.tile.list
+Keywords: hplot
+Description: Plot Dirchlet/Voronoi tiles
+URL: ../../../library/deldir/html/plot.tile.list.html
+
+Entry: tile.list
+Aliases: tile.list
+Keywords: spatial
+Description: Create a list of tiles in a tessellation.
+URL: ../../../library/deldir/html/tile.list.html
@@ -0,0 +1,14 @@
+Package: deldir
+Version: 0.0-7
+Date: 2007-11-03
+Title: Delaunay Triangulation and Dirichlet (Voronoi) Tessellation.
+Author: Rolf Turner <r.turner@auckland.ac.nz>
+Maintainer: Rolf Turner <r.turner@auckland.ac.nz>
+Depends: R (>= 0.99)
+Description: Calculates the Delaunay triangulation and the Dirichlet or
+        Voronoi tessellation (with respect to the entire plane) of a
+        planar point set.
+License: GPL version 2 or newer
+URL: http://www.math.unb.ca/~rolf/
+Packaged: Sat Nov 3 14:49:01 2007; rolf
+Built: R 2.6.0; i386-pc-mingw32; 2007-11-04 12:29:56; windows
@@ -0,0 +1,9 @@
+deldir                  Construct the Delaunay triangulation and the
+                        Dirichlet (Voronoi) tessellation of a planar
+                        point set.
+plot.deldir             Produce a plot of the Delaunay triangulation
+                        and Dirichlet (Voronoi) tesselation of a planar
+                        point set, as constructed by the function
+                        deldir.
+plot.tile.list          Plot Dirchlet/Voronoi tiles
+tile.list               Create a list of tiles in a tessellation.
@@ -0,0 +1,63 @@
+34ef60d93db4e4513ae6b13668340663 *CONTENTS
+ee8fe6e62e3e9bbcb23bc33ec089d31a *DESCRIPTION
+d8d74f53836393d466d4c787015f0bec *INDEX
+e5a045166078a4a7bd4663db86d87395 *Meta/Rd.rds
+cd368d55dede62658f2248198fa90551 *Meta/hsearch.rds
+366a9e3d80f846401b4bbea6cadc152a *Meta/package.rds
+f1e590f87144c3fc1bac5ad7ba9b43c6 *R-ex/deldir.R
+593a215bba63e9347baa5d20f399d843 *R-ex/plot.deldir.R
+21797d2ef6002f35b9539333ad936fe9 *R-ex/plot.tile.list.R
+93422ef22c461433627a006c688284c9 *R-ex/tile.list.R
+3b9c9b8cbab6a88b781daa62777d5f5c *R/deldir
+4e82e5306993afb8f50752201a1801bf *READ_ME
+af071a8a57af7c84d1333c8c402da6a0 *chtml/deldir.chm
+6b9d8dec2bc172cb280aa76cf285b67b *err.list
+b1e64b5aed9484584095ac4e02baea84 *ex.out
+16bd2bedcf7796770e90c3e56d22526e *help/AnIndex
+4d8ead26ee3fc639d613b330b6b4478e *help/deldir
+a84f7d114620c92c0ac0e2ea57265afa *help/deldir-internal
+7881f61ae9a79255ce43f0814a015899 *help/plot.deldir
+4a53379cb1c77bc2cf1c1c425cfbca8c *help/plot.tile.list
+ec37f6d9ceddf4db09dd469dedff211f *help/tile.list
+be17ba586ed30dd13fa15e313a6fd7bf *html/00Index.html
+2ab5a9ca151205217138b518e6d8405d *html/deldir-internal.html
+40f8a12dc94a21820c1ad4ab98886536 *html/deldir.html
+49398bbaf0db6df39fcdcb6f95cf593e *html/plot.deldir.html
+0cc37a506231a768a0bdf3aa3d9c14f5 *html/plot.tile.list.html
+3f9352ec46c00d0fc3f407d3103173a0 *html/tile.list.html
+1f390921c22e235996759df8afeb639a *latex/deldir-internal.tex
+6ad0d68321c5d3be37af79dce2904e68 *latex/deldir.tex
+375f59d2ef90cdfdbc38bb3aec3f45cb *latex/plot.deldir.tex
+e20c7394462a8ef65f11f3f84a2e8819 *latex/plot.tile.list.tex
+6d3ef222fc3e6eff800da2b0adbad24a *latex/tile.list.tex
+797b9359a38e5806fe599e63d173f604 *libs/deldir.dll
+f3c6c1feb137fadf8399c5541034d0c0 *man/deldir.Rd.gz
+32332f6b7cc4f78b3307c77566bb5a89 *ratfor/acchk.r
+5f8b878f079445829c1d9f49e31170bc *ratfor/addpt.r
+93ed44f787438449a98e24be66775d8b *ratfor/adjchk.r
+d53c2495df9c4d1f4d3fcad81e10b801 *ratfor/binsrt.r
+f6e55e5c74ed60b802cb2a827e597603 *ratfor/circen.r
+f54912f73cbfc10232547f97258e6ec6 *ratfor/cross.r
+43bbe616653a6276af389d5af6ff212c *ratfor/delet.r
+1791398b1b1e8db90ef7aca3d901246d *ratfor/delet1.r
+22b1b7d68d94d8d179663d3d7dbee446 *ratfor/delout.r
+6076c5cec0047fe59f8e151eef2041cf *ratfor/delseg.r
+4e775ebb151b3b47656718f9f6aefb8c *ratfor/dirout.r
+ab6664e97b26cd37a8a05ea814eada1d *ratfor/dirseg.r
+77454183e29dd932e494e9f363e7c7a0 *ratfor/dldins.r
+7a0cbe7ba69bc4d3758d7b467378a6a5 *ratfor/inddup.r
+a88b1a39361d1816685cba69561b1733 *ratfor/initad.r
+793f31965fd00d91904dc1145e5a5e95 *ratfor/insrt.r
+8d92aa51712ce4e98089b46daa29614f *ratfor/insrt1.r
+d393dd4b28b863d09b8940bee1e42e97 *ratfor/locn.r
+e6e2b439fea000bc0c4332c1332430d1 *ratfor/master.r
+adaef6c99448a066b58f1a07eb5d883c *ratfor/mnnd.r
+5e7b0ef2ba903423e2b218ca93a9a66a *ratfor/pred.r
+dd315504b6e25408c551b4f2aaf2927f *ratfor/qtest.r
+87629c704f0b27364c80453f5034f86b *ratfor/qtest1.r
+a8636bf7b643eca192877bc38a6f2c87 *ratfor/stoke.r
+6d1ed8f4862aea3c6cd9a28179c51ddf *ratfor/succ.r
+d2d7a614a5ce8ff20006261ac1db36a5 *ratfor/swap.r
+182455cc37ee801235629c44f0277830 *ratfor/testeq.r
+5264af03d65005c360f831afc0b5d160 *ratfor/triar.r
+9afccbeef89e96ce284f18504820e8a0 *ratfor/trifnd.r
@@ -0,0 +1,20 @@
+### Name: deldir
+### Title: Construct the Delaunay triangulation and the Dirichlet (Voronoi)
+###   tessellation of a planar point set.
+### Aliases: deldir
+### Keywords: spatial
+
+### ** Examples
+
+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)
+try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+# Puts dummy points at the corners of the rectangular
+# window, i.e. at (0,0), (10,0), (10,10), and (0,10)
+## Not run: 
+##D try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr')
+## End(Not run)
+# Plots the triangulation which was created (but not the tesselation).
+
+
+
@@ -0,0 +1,24 @@
+### Name: plot.deldir
+### Title: Produce a plot of the Delaunay triangulation and Dirichlet
+###   (Voronoi) tesselation of a planar point set, as constructed by the
+###   function deldir.
+### Aliases: plot.deldir
+
+
+### ** Examples
+
+## Not run: 
+##D try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+##D plot(try)
+##D #
+##D deldir(x,y,list(ndx=4,ndy=4),plot=TRUE,add=TRUE,wl='te',
+##D        col=c(1,1,2,3,4),num=TRUE)
+##D # Plots the tesselation, but does not save the results.
+##D try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr',
+##D               wp='n')
+##D # Plots the triangulation, but not the points, and saves the returned
+##D structure.
+## End(Not run)
+
+
+
@@ -0,0 +1,18 @@
+### Name: plot.tile.list
+### Title: Plot Dirchlet/Voronoi tiles
+### Aliases: plot.tile.list
+### Keywords: hplot
+
+### ** Examples
+
+  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 rainbox(20).
+  plot(w,polycol=ccc,close=TRUE)
+
+
+
@@ -0,0 +1,21 @@
+### Name: tile.list
+### Title: Create a list of tiles in a tessellation.
+### Aliases: tile.list
+### Keywords: spatial
+
+### ** Examples
+
+        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)
+
+
+
+
@@ -0,0 +1,520 @@
+.packageName <- "deldir"
+.First.lib <- function(lib,pkg) {
+	library.dynam("deldir", pkg, lib)
+	ver <- read.dcf(file.path(lib, pkg, "DESCRIPTION"), "Version")
+        cat(paste(pkg, ver, "\n"))
+}
+acw <- function(xxx) {
+xbar <- mean(xxx$x)
+ybar <- mean(xxx$y)
+theta   <- atan2(xxx$y - ybar,xxx$x-xbar)
+theta   <- ifelse(theta > 0, theta, theta + 2 * pi)
+theta.0 <- sort(unique(theta))
+iii     <- match(theta.0, theta)
+xxx$x   <- xxx$x[iii]
+xxx$y   <- xxx$y[iii]
+xxx$bp  <- xxx$bp[iii]
+xxx
+}
+deldir <- function(x,y,dpl=NULL,rw=NULL,eps=1e-9,frac=1e-4,
+                   sort=TRUE,plotit=FALSE,digits=6,...) {
+# Function deldir
+#
+#   Copyright (C) 1996 by T. Rolf Turner
+#
+#   Permission to use, copy, modify, and distribute this software and
+#   its documentation for any purpose and without fee is hereby
+#   granted, provided that the above copyright notice appear in all
+#   copies and that both that copyright notice and this permission
+#   notice appear in supporting documentation.
+#
+# ORIGINALLY PROGRAMMED BY: Rolf Turner in 1987/88, while with the
+# Division of Mathematics and Statistics, CSIRO, Sydney, Australia.
+# Re-programmed by Rolf Turner to adapt the implementation from a
+# stand-alone Fortran program to an S function, while visiting the
+# University of Western Australia, May 1995.  Further revised
+# December 1996.
+# 
+# Function to compute the Delaunay Triangulation (and hence the
+# Dirichlet Tesselation) of a planar point set according to the
+# second (iterative) algorithm of Lee and Schacter, International
+# Journal of Computer and Information Sciences, Vol. 9, No. 3, 1980,
+# pages 219 to 242.
+
+# The triangulation is made to be with respect to the whole plane by
+# `suspending' it from `ideal' points
+# (-R,-R), (R,-R) (R,R), and (-R,R), where R --> infinity.
+
+# It is also enclosed in a finite rectangle (whose boundaries truncate any
+# infinite Dirichlet tiles) with corners (xmin,ymin) etc.  This rectangle
+# is referred to elsewhere as `the' rectangular window.
+
+# If the first argument is a list, extract components x and y:
+if(is.list(x)) {
+	if(all(!is.na(match(c('x','y'),names(x))))) {
+		y <- x$y
+		x <- x$x
+	}
+	else {
+		cat('Error: called with list lacking both x and y elements\n')
+		return()
+	}
+}
+
+# If a data window is specified, get its corner coordinates
+# and truncate the data by this window:
+n <- length(x)
+if(n!=length(y)) stop('data lengths do not match')
+if(!is.null(rw)) {
+	xmin <- rw[1]
+	xmax <- rw[2]
+	ymin <- rw[3]
+	ymax <- rw[4]
+	drop <- (1:n)[x<xmin|x>xmax|y<ymin|y>ymax]
+	if(length(drop)>0) {
+		x <- x[-drop]
+		y <- y[-drop]
+		n <- length(x)
+	}
+}
+
+# If corners of the window are not specified, form them from
+# the minimum and maximum of the data +/- 10%:
+else {
+	xmin <- min(x)
+	xmax <- max(x)
+	ymin <- min(y)
+	ymax <- max(y)
+	xdff <- xmax-xmin
+	ydff <- ymax-ymin
+	xmin <- xmin-0.1*xdff
+	xmax <- xmax+0.1*xdff
+	ymin <- ymin-0.1*ydff
+	ymax <- ymax+0.1*ydff
+	rw   <- c(xmin,xmax,ymin,ymax)
+}
+
+# Add the dummy points:
+if(!is.null(dpl)) {
+	dpts <- dumpts(x,y,dpl,rw)
+	x    <- dpts$x
+	y    <- dpts$y
+}
+
+# Eliminate duplicate points:
+iii <- !ind.dup(x,y,rw,frac)
+ndm <- sum(iii[-(1:n)])
+n   <- sum(iii[1:n])
+x   <- x[iii]
+y   <- y[iii]
+
+# Make space for the total number of points (real and dummy) as
+# well as 4 ideal points and 4 extra corner points which get used
+# (only) by subroutines dirseg and dirout in the ``output'' process
+# (returning a description of the triangulation after it has been
+# calculated):
+npd  <- n + ndm
+ntot <- npd + 4               # ntot includes the 4 ideal points but
+                              # but NOT the 4 extra corners
+x <- c(rep(0,4),x,rep(0,4))
+y <- c(rep(0,4),y,rep(0,4))
+
+# Set up fixed dimensioning constants:
+ntdel  <- 4*npd
+ntdir  <- 3*npd
+
+# Set up dimensioning constants which might need to be increased:
+madj <- max(20,ceiling(3*sqrt(ntot)))
+tadj <- (madj+1)*(ntot+4)
+ndel <- madj*(madj+1)/2
+tdel <- 6*ndel
+ndir <- ndel
+tdir <- 8*ndir
+
+# Call the master subroutine to do the work:
+repeat {
+	tmp <- .Fortran(
+			'master',
+			x=as.double(x),
+			y=as.double(y),
+			sort=as.logical(sort),
+			rw=as.double(rw),
+			npd=as.integer(npd),
+			ntot=as.integer(ntot),
+			nadj=integer(tadj),
+			madj=as.integer(madj),
+			ind=integer(npd),
+			tx=double(npd),
+			ty=double(npd),
+			ilist=integer(npd),
+			eps=as.double(eps),
+			delsgs=double(tdel),
+			ndel=as.integer(ndel),
+			delsum=double(ntdel),
+			dirsgs=double(tdir),
+			ndir=as.integer(ndir),
+			dirsum=double(ntdir),
+			nerror=integer(1),
+			PACKAGE='deldir'
+		)
+	
+# Check for errors:
+	nerror <- tmp$nerror
+	if(nerror < 0) break
+
+	else {
+		if(nerror==4) {
+			cat('nerror =',nerror,'\n')
+			cat('Increasing madj and trying again.\n')
+			madj <- ceiling(1.2*madj)
+			tadj <- (madj+1)*(ntot+4)
+			ndel <- max(ndel,madj*(madj+1)/2)
+			tdel <- 6*ndel
+			ndir <- ndel
+			tdir <- 8*ndir
+			}
+		else if(nerror==14|nerror==15) {
+			cat('nerror =',nerror,'\n')
+			cat('Increasing ndel and ndir and trying again.\n')
+			ndel <- ceiling(1.2*ndel)
+			tdel <- 6*ndel
+	                ndir <- ndel
+	                tdir <- 8*ndir
+		}
+		else {
+			cat('nerror =',nerror,'\n')
+			return(invisible())
+		}
+	}
+}
+
+# Collect up the results for return:
+ndel       <- tmp$ndel
+delsgs     <- round(t(as.matrix(matrix(tmp$delsgs,nrow=6)[,1:ndel])),digits)
+delsum     <- matrix(tmp$delsum,ncol=4)
+del.area   <- sum(delsum[,4])
+delsum     <- round(cbind(delsum,delsum[,4]/del.area),digits)
+del.area   <- round(del.area,digits)
+ndir       <- tmp$ndir
+dirsgs     <- round(t(as.matrix(matrix(tmp$dirsgs,nrow=8)[,1:ndir])),digits)
+dirsgs     <- as.data.frame(dirsgs)
+dirsum     <- matrix(tmp$dirsum,ncol=3)
+dir.area   <- sum(dirsum[,3])
+dirsum     <- round(cbind(dirsum,dirsum[,3]/dir.area),digits)
+dir.area   <- round(dir.area,digits)
+allsum     <- cbind(delsum,dirsum)
+rw         <- round(rw,digits)
+
+# Name the columns of the results:
+dimnames(delsgs) <- list(NULL,c('x1','y1','x2','y2','ind1','ind2'))
+names(dirsgs)    <- c('x1','y1','x2','y2','ind1','ind2','bp1','bp2')
+mode(dirsgs$bp1) <- 'logical'
+mode(dirsgs$bp2) <- 'logical'
+dimnames(allsum) <- list(NULL,c('x','y','n.tri','del.area','del.wts',
+                                'n.tside','nbpt','dir.area','dir.wts'))
+
+# Aw' done!!!
+rslt <- list(delsgs=delsgs,dirsgs=dirsgs,summary=allsum,n.data=n,
+             n.dum=ndm,del.area=del.area,dir.area=dir.area,rw=rw)
+class(rslt) <- 'deldir'
+if(plotit) plot(rslt,...)
+if(plotit) invisible(rslt) else rslt
+}
+dumpts <- function(x,y,dpl,rw) {
+#
+# Function dumpts to append a sequence of dummy points to the
+# data points.
+#
+
+ndm  <- 0
+xd   <- NULL
+yd   <- NULL
+xmin <- rw[1]
+xmax <- rw[2]
+ymin <- rw[3]
+ymax <- rw[4]
+
+# Points on radii of circles emanating from data points:
+if(!is.null(dpl$nrad)) {
+	nrad  <- dpl$nrad # Number of radii from each data point.
+	nper  <- dpl$nper # Number of dummy points per radius.
+	fctr  <- dpl$fctr # Length of each radius = fctr * mean
+                             # interpoint distance.
+	lrad  <- fctr*mnnd(x,y)/nper
+	theta <- 2*pi*(1:nrad)/nrad
+	cs    <- cos(theta)
+	sn    <- sin(theta)
+	xt    <- c(lrad*(1:nper)%o%cs)
+	yt    <- c(lrad*(1:nper)%o%sn)
+	xd    <- c(outer(x,xt,'+'))
+	yd    <- c(outer(y,yt,'+'))
+}
+
+# Ad hoc points passed over as part of dpl:
+if(!is.null(dpl$x)) {
+	xd <- c(xd,dpl$x)
+	yd <- c(yd,dpl$y)
+}
+
+# Delete dummy points outside the rectangular window.
+ndm  <- length(xd)
+if(ndm >0) {
+	drop <- (1:ndm)[xd<xmin|xd>xmax|yd<ymin|yd>ymax]
+	if(length(drop)>0) {
+		xd  <- xd[-drop]
+		yd  <- yd[-drop]
+	}
+}
+
+# Rectangular grid:
+ndx <- dpl$ndx
+okx <- !is.null(ndx) && ndx > 0
+ndy <- dpl$ndy
+oky <- !is.null(ndy) && ndy > 0
+if(okx & oky) {
+	xt  <- if(ndx>1) seq(xmin,xmax,length=ndx) else 0.5*(xmin+xmax)
+	yt  <- if(ndy>1) seq(ymin,ymax,length=ndy) else 0.5*(ymin+ymax)
+	xy <- expand.grid(x=xt,y=yt)
+	xd  <- c(xd,xy$x)
+	yd  <- c(yd,xy$y)
+}
+
+ndm <- length(xd)
+list(x=c(x,xd),y=c(y,yd),ndm=ndm)
+}
+get.cnrind <- function(x,y,rw) {
+x.crnrs <- rw[c(1,2,2,1)]
+y.crnrs <- rw[c(3,3,4,4)]
+M1 <- outer(x,x.crnrs,function(a,b){(a-b)^2})
+M2 <- outer(y,y.crnrs,function(a,b){(a-b)^2})
+MM <- M1 + M2
+apply(MM,2,which.min)
+}
+ind.dup <- function(x,y,rw=NULL,frac=0.0001) {
+#
+# Function ind.dup to calculate the indices of data pairs 
+# which duplicate earlier ones.  (Returns a logical vector;
+# true for such indices, false for the rest.)
+#
+
+if(is.null(rw)) rw <- c(0,1,0,1)
+n <- length(x)
+rslt <- .Fortran(
+		'inddup',
+		x=as.double(x),
+		y=as.double(y),
+		n=as.integer(n),
+		rw=as.double(rw),
+		frac=as.double(frac),
+		dup=logical(n),
+		PACKAGE='deldir'
+	)
+
+rslt$dup
+}
+mid.in <- function(x,y,rx,ry) {
+xm <- 0.5*(x[1]+x[2])
+ym <- 0.5*(y[1]+y[2])
+(rx[1] < xm & xm < rx[2] & ry[1] < ym & ym < ry[2])
+}
+mnnd <- function(x,y) {
+#
+# Function mnnd to calculate the mean nearest neighbour distance
+# between the points whose coordinates are stored in x and y.
+#
+
+n <- length(x)
+if(n!=length(y)) stop('data lengths do not match')
+dmb <- (max(x)-min(x))**2 + (max(y)-min(y))**2
+
+.Fortran(
+	"mnnd",
+	x=as.double(x),
+	y=as.double(y),
+	n=as.integer(n),
+	dmb=as.double(dmb),
+	d=double(1),
+	PACKAGE='deldir'
+	)$d
+}
+plot.deldir <- function(x,add=FALSE,wlines=c('both','triang','tess'),
+                        wpoints=c('both','real','dummy','none'),
+                        number=FALSE,cex=1,nex=1,col=NULL,lty=NULL,
+                        pch=NULL,xlim=NULL,ylim=NULL,xlab='x',ylab='y',...)
+{
+#
+# Function plot.deldir to produce a plot of the Delaunay triangulation
+# and Dirichlet tesselation of a point set, as produced by the
+# function deldir().
+#
+
+wlines  <- match.arg(wlines)
+wpoints <- match.arg(wpoints)
+
+if(is.null(class(x)) || class(x)!='deldir') {
+	cat('Argument is not of class deldir.\n')
+	return(invisible())
+}
+
+col <- if(is.null(col)) c(1,1,1,1,1) else rep(col,length.out=5)
+lty <- if(is.null(lty)) 1:2 else rep(lty,length.out=2)
+pch <- if(is.null(pch)) 1:2 else rep(pch,length.out=2)
+
+plot.del <- switch(wlines,both=TRUE,triang=TRUE,tess=FALSE)
+plot.dir <- switch(wlines,both=TRUE,triang=FALSE,tess=TRUE)
+plot.rl  <- switch(wpoints,both=TRUE,real=TRUE,dummy=FALSE,none=FALSE)
+plot.dum <- switch(wpoints,both=TRUE,real=FALSE,dummy=TRUE,none=FALSE)
+
+delsgs <- x$delsgs
+dirsgs <- x$dirsgs
+n      <- x$n.data
+rw     <- x$rw
+
+if(plot.del) {
+	x1<-delsgs[,1]
+	y1<-delsgs[,2]
+	x2<-delsgs[,3]
+	y2<-delsgs[,4]
+}
+
+if(plot.dir) {
+	u1<-dirsgs[,1]
+	v1<-dirsgs[,2]
+	u2<-dirsgs[,3]
+	v2<-dirsgs[,4]
+}
+
+X<-x$summary[,1]
+Y<-x$summary[,2]
+
+if(!add) {
+	pty.save <- par()$pty
+	on.exit(par(pty=pty.save))
+	par(pty='s')
+	if(is.null(xlim)) xlim <- rw[1:2]
+	if(is.null(ylim)) ylim <- rw[3:4]
+	plot(0,0,type='n',xlim=xlim,ylim=ylim,
+     		xlab=xlab,ylab=ylab,axes=FALSE,...)
+	axis(side=1)
+	axis(side=2)
+}
+
+if(plot.del) segments(x1,y1,x2,y2,col=col[1],lty=lty[1],...)
+if(plot.dir) segments(u1,v1,u2,v2,col=col[2],lty=lty[2],...)
+if(plot.rl) {
+	x.real <- X[1:n]
+	y.real <- Y[1:n]
+	points(x.real,y.real,pch=pch[1],col=col[3],cex=cex,...)
+}
+if(plot.dum) {
+	x.dumm <- X[-(1:n)]
+	y.dumm <- Y[-(1:n)]
+	points(x.dumm,y.dumm,pch=pch[2],col=col[4],cex=cex,...)
+}
+if(number) {
+	xoff <-0.02*diff(range(X))
+	yoff <-0.02*diff(range(Y))
+	text(X+xoff,Y+yoff,1:length(X),cex=nex,col=col[5],...)
+}
+invisible()
+}
+`plot.tile.list` <-
+function (x, verbose = FALSE, close = FALSE, pch = 1, polycol = NA, 
+    showpoints = TRUE, asp = 1, ...) 
+{
+    object <- x
+    if (!inherits(object, "tile.list")) 
+        stop("Argument \"object\" is not of class tile.list.\n")
+    n <- length(object)
+    x.all <- unlist(lapply(object, function(w) {
+        c(w$pt[1], w$x)
+    }))
+    y.all <- unlist(lapply(object, function(w) {
+        c(w$pt[2], w$y)
+    }))
+    x.pts <- unlist(lapply(object, function(w) {
+        w$pt[1]
+    }))
+    y.pts <- unlist(lapply(object, function(w) {
+        w$pt[2]
+    }))
+    rx <- range(x.all)
+    ry <- range(y.all)
+    plot(x.all, y.all, type = "n", asp = asp, xlab = "x", ylab = "y")
+    polycol <- apply(col2rgb(polycol,TRUE),2,
+                     function(x){do.call(rgb,as.list(x/255))})
+    polycol <- rep(polycol, length = length(object))
+    hexbla  <- do.call(rgb,as.list(col2rgb("black",TRUE)/255))
+    hexwhi  <- do.call(rgb,as.list(col2rgb("white",TRUE)/255))
+    ptcol <- ifelse(polycol == hexbla,hexwhi,hexbla)
+    lnwid <- ifelse(polycol == hexbla, 2, 1)
+    for (i in 1:n) {
+        inner <- !any(object[[i]]$bp)
+        if (close | inner) 
+            polygon(object[[i]], col = polycol[i], border = ptcol[i], 
+                lwd = lnwid[i])
+        else {
+            x <- object[[i]]$x
+            y <- object[[i]]$y
+            bp <- object[[i]]$bp
+            ni <- length(x)
+            for (j in 1:ni) {
+                jnext <- if (j < ni) 
+                  j + 1
+                else 1
+                do.it <- mid.in(x[c(j, jnext)], y[c(j, jnext)], 
+                  rx, ry)
+                if (do.it) 
+                  segments(x[j], y[j], x[jnext], y[jnext], col = ptcol[i], 
+                    lwd = lnwid[i])
+            }
+        }
+        if (verbose & showpoints) 
+            points(object[[i]]$pt[1], object[[i]]$pt[2], pch = pch, 
+                col = ptcol[i])
+        if (verbose & i < n) 
+            readline("Go? ")
+    }
+    if (showpoints) 
+        points(x.pts, y.pts, pch = pch, col = ptcol)
+    invisible()
+}
+tile.list <- function (object) 
+{
+    if (!inherits(object, "deldir")) 
+        stop("Argument \"object\" is not of class deldir.\n")
+    rw <- object$rw
+    x.crnrs <- rw[c(1,2,2,1)]
+    y.crnrs <- rw[c(3,3,4,4)]
+    ddd <- object$dirsgs
+    sss <- object$summary
+    npts <- nrow(sss)
+    x <- sss[,"x"]
+    y <- sss[,"y"]
+    i.crnr <- get.cnrind(x,y,rw)
+    rslt <- list()
+    for (i in 1:npts) {
+        m <- as.matrix(rbind(ddd[ddd$ind1 == i, 1:4], ddd[ddd$ind2 == i, 1:4]))
+        bp1 <- c(ddd[ddd$ind1 == i, 7], ddd[ddd$ind2 == i, 7])
+        bp2 <- c(ddd[ddd$ind1 == i, 8], ddd[ddd$ind2 == i, 8])
+        m1 <- cbind(m[, 1:2, drop = FALSE], 0 + bp1)
+        m2 <- cbind(m[, 3:4, drop = FALSE], 0 + bp2)
+        m <- rbind(m1, m2)
+        pt <- sss[i, 1:2]
+        theta <- atan2(m[, 2] - pt[2], m[, 1] - pt[1])
+        theta <- ifelse(theta > 0, theta, theta + 2 * pi)
+        theta.0 <- sort(unique(theta))
+        mm <- m[match(theta.0, theta), ]
+        xx <- mm[, 1]
+        yy <- mm[, 2]
+        bp <- as.logical(mm[, 3])
+# Add corner points if necessary:
+        ii <- i.crnr%in%i
+	xx <- c(xx,x.crnrs[ii])
+	yy <- c(yy,y.crnrs[ii])
+        bp <- c(bp,rep(TRUE,sum(ii)))
+        rslt[[i]] <- acw(list(pt=pt, x = xx, y = yy, bp = bp))
+    }
+    class(rslt) <- "tile.list"
+    rslt
+}
@@ -0,0 +1,186 @@
+
+ Version date: 21 February 2002.
+ This version is simply an adaptation of the Splus version
+ of the package to R.
+
+                       *****************************
+
+ Version date: 14 February 2002.
+ This version supercedes the version dated 24 April 1999.
+
+                       *****************************
+
+ The changes from the version dated 24 April 1999 to the
+ version dated 14 February 2002 were:
+
+	A bug in the procedure for eliminating duplicated points was
+	fixed.  Thanks go to Dr. Berwin Turlach of the Department of
+	Maths and Stats at the University of Western Australia, for
+	spotting this bug.
+
+                       *****************************
+
+ The changes from the version dated 26 October 1998 to the
+ version dated 24 April 1999 were:
+
+	(1) The function mipd(), stored in mipd.sf, and the
+	corresponding Fortran subroutine mipd, stored in mipd.r, have
+	been replaced by mnnd() in mnnd.sf and mnnd in mnnd.r.  The
+	function mipd calculated the mean interpoint distance, to be
+	used in constructing dummy point structures of a certain
+	type.  After some reflection it became apparent that the mean
+	interpoint distance was much too large for the intended
+	purpose, and that a more appropriate value was the ``mean
+	nearest neighbour distance'' which is calculated by the new
+	function.  This new value is now used in constructing dummy
+	point structures.
+
+	Note that the operative result is that the resulting dummy
+	point structures contain many more points than before.  The
+	old value caused large numbers of the dummy points to fall
+	outside the data window and therefore to be clipped.
+
+                       *****************************
+
+ The changes from the version dated 6 December 1996 to the
+ version dated 26 October 1998 were:
+
+	(1) A ratfor/Fortran routine named ``inside'' has been
+	renamed ``dldins'' to avoid conflict with a name built in to
+	some versions of Splus.
+
+	(2) Some minor corrections have been made to dangerous
+	infelicities in a piece of the ratfor/Fortran code.
+
+	(3) The dynamic loading procedure has been changed to use
+	dyn.load.shared so that the package is easily usable on
+	IRIX systems as well as under SunOS/Solaris.
+
+	(4) The package has been adjusted slightly so that it
+	can easily be installed as a section of a library.  In
+	particular, the dynamic loading is now done by the
+	.First.lib() function rather than from within deldir()
+	itself; reference to an environment variable DYN_LOAD_LIB
+	is no longer needed.
+
+                       *****************************
+
+ This package computes and plots the Dirichlet (Voronoi) tesselation
+ and the Delaunay triangulation of a set of of data points and
+ possibly a superimposed ``grid'' of dummy points.
+
+ The tesselation is constructed with respect to the whole plane
+ by suspending it from ideal points at infinity.
+
+ ORIGINALLY PROGRAMMED BY: Rolf Turner in 1987/88, while with the
+ Division of Mathematics and Statistics, CSIRO, Sydney, Australia.
+ Re-programmed by Rolf Turner to adapt the implementation from a
+ stand-alone Fortran program to an S function, while visiting the
+ University of Western Australia, May 1995.  Further revised
+ December 1996, October 1998, April 1999, and February 2002.
+ Adapted to an R package 21 February 2002.
+
+ Current address of the author:
+		Department of Mathematics and Statistics,
+                University of New Brunswick,
+                P.O. Box 4400, Fredericton, New Brunswick,
+                Canada E3B 5A3
+ Email:
+		rolf@math.unb.ca
+
+ 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 the University of
+ Western Australia).  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, is now deprecated in
+ favour of this current package.  Don also collaborated in the
+ preparation of this current package.
+
+ Bill Dunlap of MathSoft Inc. tracked down a bug which was making
+ the deldir() function crash on some systems, and pointed out some
+ other improvements to be made.
+
+ Berwin Turlach of the Department of Maths and Stats at the
+ University of Western Australia pointed out a bug in the procedure
+ for eliminating duplicated points.
+
+                       *****************************
+
+ The man directory, contains, in addition to the R documentation
+ files deldir.Rd and plot.deldir.Rd:
+
+	(a) This READ_ME file.
+	(b) A file err.list, containing a list of meanings of possible
+	    error numbers which could be returned.  NONE of these
+	    errors should ever actually happen except for errors 4, 14,
+	    and 15.  These relate to insufficient dimensioning, and
+	    if they occur, the driver increases the dimensions and
+	    tries again (informing you of this fact).
+	(c) A file ex.out containing a printout of the object returned
+	    by running the example given in the help file for deldir.
+
+ The src directory contains many, many *.f (Fortran) files, which
+ get compiled and dynamically loaded.
+
+ The Fortran code is ponderous --- it was automatically generated
+ from Ratfor code, which was pretty ponderous to start with.   It is
+ quite possibly very kludgy aw well --- i.e. a good programmer could
+ make it ***much*** more efficient I'm sure.  It contains all sorts
+ of checking for anomalous situations which probably can/will never
+ occur.  These checks basically reflect my pessimism and fervent
+ belief in Murphy's Law.
+
+ The program  was also designed with a particular application in mind,
+ in which we wished to superimpose a grid of dummy points onto the
+ actual data points which we were triangulating.  This fact adds
+ slightly to the complication of the code.
+
+                       *****************************
+
+ Here follows a brief description of the package:
+
+ (1) The function deldir computes the Delaunay Triangulation (and
+ hence the Dirichlet Tesselation) of a planar point set according to
+ the second (iterative) algorithm of Lee and Schacter, International
+ Journal of Computer and Information Sciences, Vol. 9, No. 3, 1980,
+ pages 219 to 242.
+
+ The tesselation/triangulation is made to be
+
+             **** with respect to the whole plane ****
+
+ by `suspending' it from `ideal' points (-Inf,-Inf), (Inf,-Inf)
+ (Inf,Inf), and (-Inf,Inf).
+
+ (2) The tesselation/triangulation is also enclosed in a finite
+ rectangle with corners
+
+          (xmin,ymax) * ------------------------ * (xmax,ymax)
+                      |                          |
+                      |                          |
+                      |                          |
+                      |                          |
+                      |                          |
+          (xmin,ymin) * ------------------------ * (xmax,ymin)
+
+ The boundaries of this rectangle truncate some Dirichlet tiles, in
+ particular any infinite ones. This rectangle is referred to
+ elsewhere as `the' rectangular window.
+               ===
+
+ (2) The function plot.deldir is a method for plot.  I.e. it may be
+ invoked simply by typing ``plot(x)'' provided that ``x'' is an
+ object of class ``deldir'' (as produced by the function deldir).
+ The plot (by default) consists of the edges of the Delaunay
+ triangles (solid lines) and the edges of the Dirichlet tiles (dotted
+ lines).  By default the real data points are indicated by circles,
+ and the dummy points are indicated by triangles.
@@ -0,0 +1,49 @@
+
+Error  list:
+===========
+
+nerror =  1: Contradictory adjacency lists.  Error in adjchk.
+
+nerror =  2: Number of points jumbled.  Error in binsrt.
+
+nerror =  3: Vertices of 'triangle' are collinear
+             and vertex 2 is not between 1 and 3.  Error in circen.
+
+nerror =  4: Number of adjacencies too large.  Error in insrt.
+             (Automatically adjusted for in deldir().)
+
+nerror =  5: Adjacency list of i is empty, and so cannot contain j.
+             Error in pred.
+
+nerror =  6: Adjacency list of i does not contain j.  Error in pred.
+
+nerror =  7: Indicator ijk is out of range. (This CAN'T happen!)
+             Error in qtest.
+
+nerror =  8: Fell through all six cases.
+             Something must be totally stuffed up.  Error in stoke.
+
+nerror =  9: Adjacency list of i is empty, and so cannot contain j.
+             Error in succ.
+
+nerror = 10: Adjacency list of i does not contain j.  Error in succ.
+
+nerror = 11: No triangles to find. Error in trifnd.
+
+nerror = 12: Vertices of triangle are collinear.  Error in dirseg.
+
+nerror = 13: Vertices of triangle are collinear.  Error in dirout.
+
+nerror = 14: Number of Delaunay segments exceeds alloted space.
+             Error in delseg.
+             (Automatically adjusted for in deldir().)
+
+nerror = 15: Number of Dirichlet segments exceeds alloted space.
+             Error in dirseg.
+             (Automatically adjusted for in deldir().)
+
+nerror = 16: Line from midpoint to circumcenter does not intersect
+             rectangle boundary; but it HAS to!!! Error in dirseg.
+
+nerror = 17: Line from midpoint to circumcenter does not intersect
+             rectangle boundary; but it HAS to!!! Error in dirout.
@@ -0,0 +1,74 @@
+$delsgs
+      x1 y1   x2   y2 ind1 ind2
+ [1,]  3  3  2.3  2.3    2    1
+ [2,]  7  2  2.3  2.3    3    1
+ [3,]  7  2  3.0  3.0    3    2
+ [4,]  1  5  2.3  2.3    4    1
+ [5,]  1  5  3.0  3.0    4    2
+ [6,]  3  8  3.0  3.0    5    2
+ [7,]  3  8  7.0  2.0    5    3
+ [8,]  3  8  1.0  5.0    5    4
+ [9,]  8  9  7.0  2.0    6    3
+[10,]  8  9  3.0  8.0    6    5
+[11,]  0  0  2.3  2.3    7    1
+[12,]  0  0  7.0  2.0    7    3
+[13,]  0  0  1.0  5.0    7    4
+[14,] 10  0  7.0  2.0    8    3
+[15,] 10  0  8.0  9.0    8    6
+[16,] 10  0  0.0  0.0    8    7
+[17,]  0 10  1.0  5.0    9    4
+[18,]  0 10  3.0  8.0    9    5
+[19,]  0 10  8.0  9.0    9    6
+[20,]  0 10  0.0  0.0    9    7
+[21,] 10 10  8.0  9.0   10    6
+[22,] 10 10 10.0  0.0   10    8
+[23,] 10 10  0.0 10.0   10    9
+
+$dirsgs
+          x1        y1        x2        y2 ind1 ind2   bp1   bp2
+1   1.650000  3.650000  4.560000  0.740000    2    1 FALSE FALSE
+2   4.560000  0.740000  4.512766  0.000000    3    1 FALSE  TRUE
+3   5.750000  5.500000  4.560000  0.740000    3    2 FALSE FALSE
+4   0.000000  2.855556  1.650000  3.650000    4    1  TRUE FALSE
+5   1.650000  3.650000  3.500000  5.500000    4    2 FALSE FALSE
+6   3.500000  5.500000  5.750000  5.500000    5    2 FALSE FALSE
+7   5.750000  5.500000  6.058824  5.705882    5    3 FALSE FALSE
+8   0.500000  7.500000  3.500000  5.500000    5    4 FALSE FALSE
+9   6.058824  5.705882 10.000000  5.142857    6    3 FALSE  TRUE
+10  5.200000 10.000000  6.058824  5.705882    6    5  TRUE FALSE
+11  2.300000  0.000000  0.000000  2.300000    7    1  TRUE  TRUE
+12 10.000000  3.250000  7.833333  0.000000    8    3  TRUE  TRUE
+13  0.000000  7.400000  0.500000  7.500000    9    4  TRUE FALSE
+14  0.500000  7.500000  2.166667 10.000000    9    5 FALSE  TRUE
+15  8.750000 10.000000 10.000000  7.500000   10    6  TRUE  TRUE
+
+$summary
+         x    y n.tri  del.area  del.wts n.tside nbpt  dir.area  dir.wts
+ [1,]  2.3  2.3     4  4.500000 0.045000       4    4  9.092057 0.090921
+ [2,]  3.0  3.0     4  6.050000 0.060500       4    0 10.738500 0.107385
+ [3,]  7.0  2.0     6 18.666667 0.186667       5    4 23.318162 0.233182
+ [4,]  1.0  5.0     5  7.500000 0.075000       4    2  9.394167 0.093942
+ [5,]  3.0  8.0     5 15.000000 0.150000       5    2 18.055637 0.180556
+ [6,]  8.0  9.0     5 16.666667 0.166667       3    4 18.314811 0.183148
+ [7,]  0.0  0.0     4  8.450000 0.084500       1    2  2.645000 0.026450
+ [8,] 10.0  0.0     3 10.500000 0.105000       1    2  3.520833 0.035208
+ [9,]  0.0 10.0     4  7.666667 0.076667       2    2  3.358333 0.033583
+[10,] 10.0 10.0     2  5.000000 0.050000       1    2  1.562500 0.015625
+
+$n.data
+[1] 6
+
+$n.dum
+[1] 4
+
+$del.area
+[1] 100
+
+$dir.area
+[1] 100
+
+$rw
+[1]  0 10  0 10
+
+attr(,"class")
+[1] "deldir"
@@ -0,0 +1,10 @@
+acw	deldir-internal
+deldir	deldir
+dumpts	deldir-internal
+get.cnrind	deldir-internal
+ind.dup	deldir-internal
+mid.in	deldir-internal
+mnnd	deldir-internal
+plot.deldir	plot.deldir
+plot.tile.list	plot.tile.list
+tile.list	tile.list
@@ -0,0 +1,224 @@
+deldir                package:deldir                R Documentation
+
+_C_o_n_s_t_r_u_c_t _t_h_e _D_e_l_a_u_n_a_y _t_r_i_a_n_g_u_l_a_t_i_o_n _a_n_d _t_h_e _D_i_r_i_c_h_l_e_t
+(_V_o_r_o_n_o_i) _t_e_s_s_e_l_l_a_t_i_o_n _o_f _a _p_l_a_n_a_r _p_o_i_n_t _s_e_t.
+
+_D_e_s_c_r_i_p_t_i_o_n:
+
+     This function computes the Delaunay triangulation (and hence the
+     Dirichlet 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.
+
+_U_s_a_g_e:
+
+     deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, frac=0.0001,
+            sort=TRUE, plotit=FALSE, digits=6, ...)
+
+_A_r_g_u_m_e_n_t_s:
+
+     x,y: The coordinates of the point set being triangulated. These
+          can be given by two arguments x and y which are vectors or by
+          a single argument x which is a list with components "x" and
+          "y". 
+
+     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 auxilliary function dumpts().
+           The list may have components:
+
+          ndx: The x-dimension of a rectangular grid; if either ndx or
+          ndy is null, no grid is constructed.
+
+          ndy: The y-dimension of the aforementioned rectangular grid.
+
+          nrad: The number of radii or "spokes", emanating from each
+          data point, along which dummy points are to be added.
+
+          nper: The number of dummy points per spoke.
+
+          fctr: A factor determining the length of each spoke; each
+          spoke is of length equal to fctr times the mean nearest
+          neighbour distance of the data. (This distance is calculated
+          by the auxilliary function mnnd().)
+
+          x: A vector of x-coordinates of "ad hoc" dummy points
+
+          y: A vector of the corresponding y-coordinates of "ad hoc"
+          dummy points
+
+      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. 
+
+    frac: A value specifying the tolerance used in eliminating
+          duplicate points; defaults to 0.0001. Points are considered
+          duplicates if abs(x1-x2) < frac*(xmax-xmin) AND abs(y1-y2) <
+          frac*(ymax-ymin). 
+
+    sort: Logical argument; if 'TRUE' (the default) the data (including
+          dummy points) are sorted into a sequence of "bins" prior to
+          triangulation; this makes the algorithm slightly more
+          efficient.  Normally one would set sort equal to 'FALSE' only
+          if one wished to observe some of the fine detail of the way
+          in which adding a point to a data set affected the
+          triangulation, and therefore wished to make sure that the
+          point in question was added last.  Essentially this argument
+          would get used only in a de-bugging process. 
+
+  plotit: Logical argument; if 'TRUE' a plot of the triangulation and
+          tessellation is produced; the default is 'FALSE'. 
+
+  digits: The number of decimal places to which all numeric values in
+          the returned list should be rounded.  Defaults to 6. 
+
+     ...: Auxilliary arguments add, wlines, wpoints, number, nex, col,
+          lty, pch, xlim, and ylim (and possibly other plotting
+          parameters) may be passed to plot.deldir through "..." if
+          plotit='TRUE'. 
+
+_D_e_t_a_i_l_s:
+
+     This package is a (straightforward) adaptation of the Splus
+     library section ``delaunay'' to R.  That library section is an
+     implementation of the Lee-Schacter algorithm, which was originally
+     written as a stand-alone Fortran program in 1987/88 by Rolf
+     Turner, while with the Division of Mathematics and Statistics,
+     CSIRO, Sydney, Australia.  It was re-written as an Splus function
+     (using dynamically loaded Fortran code), by Rolf Turner while
+     visiting the University of Western Australia, May, 1995.
+
+     Further revisions 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 the University of Western Australia). 
+     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, is now deprecated in
+     favour of the current package ``delaunay'' package.  Don also
+     collaborated in the preparation of that package.
+
+     Further revisions and bug-fixes were made in 1998, 1999, and 2002.
+
+_V_a_l_u_e:
+
+     A list (of class 'deldir'), invisible if plotit='TRUE', with
+     components:
+
+  delsgs: a matrix 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 (x1,y1,x2,y2).  The last two entries
+          are the indices of the two points which are joined. 
+
+  dirsgs: a data frame with 8 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 (x1,y1,x2,y2).  The fifth and
+          sixth entries are the indices of the two points, in the set
+          being triangulated, which are separated by that edge. The
+          seventh and eighth entries are logical values.  The seventh
+          indicates whether the first endpoint of the corresponding
+          edge of a Dirichlet tile is a boundary point (a point on the
+          boundary of the rectangular window).  Likewise for the eighth
+          entry and the second endpoint of the edge. 
+
+ summary: a matrix with 9 columns, and (n.data + n.dum) rows (see
+          below). These rows correspond to the points in the set being
+          triangulated. The columns are named "x" (the x-coordinate of
+          the point), "y" (the y-coordinate), "n.tri" (the number of
+          Delaunay triangles emanating from the point), "del.area" (1/3
+          of the total area of all the Delaunay triangles emanating
+          from the point), "del.wts" (the corresponding entry of the
+          "del.area" column divided by the sum of this column),
+          "n.tside" (the number of sides - within the rectangular
+          window - of the Dirichlet tile surrounding the point), "nbpt"
+          (the number of points in which the Dirichlet tile intersects
+          the boundary of the rectangular window), "dir.area" (the area
+          of the Dirichlet tile surrounding the point), and "dir.wts"
+          (the corresponding entry of the "dir.area" column divided by
+          the sum of this column).  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 "summary" correspond to real
+          points. 
+
+   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
+          "summary" correspond to dummy points. 
+
+del.area: the area of the convex hull of the set of points being
+          triangulated, as formed by summing the "del.area" column of
+          "summary". 
+
+dir.area: the area of the rectangular window enclosing the points being
+          triangulated, as formed by summing the "dir.area" column of
+          "summary". 
+
+      rw: the specification of the corners of the rectangular window
+          enclosing the data, in the order (xmin, xmax, ymin, ymax). 
+
+_S_i_d_e _E_f_f_e_c_t_s:
+
+     If plotit=='TRUE' a plot of the triangulation and/or tessellation
+     is produced or added to an existing plot.
+
+_N_o_t_e:
+
+     If ndx >= 2 and ndy >= 2, then the rectangular window IS the
+     convex hull, and so the values of del.area and dir.area are
+     identical.
+
+_A_u_t_h_o_r(_s):
+
+     Rolf Turner r.turner@auckland.ac.nz <URL:
+     http://www.math.unb.ca/~rolf>
+
+_R_e_f_e_r_e_n_c_e_s:
+
+     Lee, D. T., and Schacter, B. J.  "Two algorithms for constructing
+     a Delaunay triangulation", Int. J. Computer and Information
+     Sciences, Vol. 9, No. 3, 1980, pp. 219 - 242.
+
+     Ahuja, N. and Schacter, B. J. (1983). Pattern Models.  New York:
+     Wiley.
+
+_S_e_e _A_l_s_o:
+
+     plot.deldir
+
+_E_x_a_m_p_l_e_s:
+
+     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)
+     try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+     # Puts dummy points at the corners of the rectangular
+     # window, i.e. at (0,0), (10,0), (10,10), and (0,10)
+     ## Not run: 
+     try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr')
+     ## End(Not run)
+     # Plots the triangulation which was created (but not the tesselation).
+
@@ -0,0 +1,22 @@
+deldir-internal            package:deldir            R Documentation
+
+_I_n_t_e_r_n_a_l _d_e_l_d_i_r _f_u_n_c_t_i_o_n_s
+
+_D_e_s_c_r_i_p_t_i_o_n:
+
+     Internal deldir functions.
+
+_U_s_a_g_e:
+
+     dumpts(x,y,dpl,rw)
+     ind.dup(x,y,rw=NULL,frac=0.0001)
+     mid.in(x,y,rx,ry)
+     mnnd(x,y)
+     get.cnrind(x,y,rw)
+     acw(xxx)
+
+_D_e_t_a_i_l_s:
+
+     These functions are auxilliary and are not intended to be called
+     by the user.
+
@@ -0,0 +1,118 @@
+plot.deldir              package:deldir              R Documentation
+
+_P_r_o_d_u_c_e _a _p_l_o_t _o_f _t_h_e _D_e_l_a_u_n_a_y _t_r_i_a_n_g_u_l_a_t_i_o_n _a_n_d _D_i_r_i_c_h_l_e_t (_V_o_r_o_n_o_i)
+_t_e_s_s_e_l_a_t_i_o_n _o_f _a _p_l_a_n_a_r _p_o_i_n_t _s_e_t, _a_s _c_o_n_s_t_r_u_c_t_e_d _b_y _t_h_e _f_u_n_c_t_i_o_n _d_e_l_d_i_r.
+
+_D_e_s_c_r_i_p_t_i_o_n:
+
+     This is a method for plot.
+
+_U_s_a_g_e:
+
+     ## S3 method for class 'deldir':
+     plot(x,add=FALSE,wlines=c('both','triang','tess'),
+                           wpoints=c('both','real','dummy','none'),
+                           number=FALSE,cex=1,nex=1,col=NULL,lty=NULL,
+                           pch=NULL,xlim=NULL,ylim=NULL,xlab='x',ylab='y',...)
+
+_A_r_g_u_m_e_n_t_s:
+
+       x: An object of class "deldir" as constructed 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) ? 
+
+ wpoints: "which points?".  I.e.  should the real points be plotted
+          (wpoints='real'), should the dummy points be plotted
+          (wpoints='dummy'), should both be plotted (wpoints='both',
+          the default) or should no points be plotted (wpoints='none')? 
+
+  number: Logical argument, defaulting to 'FALSE'; if 'TRUE' then the
+          points plotted will be labelled with their index numbers
+          (corresponding to the row numbers of the matrix "summary" in
+          the output of deldir). 
+
+     cex: The value of the character expansion argument cex to be used
+          with the plotting symbols for plotting the points. 
+
+     nex: The value of the character expansion argument cex to be used
+          by the text function when numbering the points with their
+          indices.  Used only if number='TRUE'. 
+
+     col: the colour numbers for plotting the triangulation, the
+          tesselation, the data points, the dummy points, and the point
+          numbers, in that order; defaults to c(1,1,1,1,1).  If fewer
+          than five numbers are given, they are recycled.  (If more
+          than five numbers are given, the redundant ones are ignored.) 
+
+     lty: the line type numbers for plotting the triangulation and the
+          tesselation, in that order; defaults to 1:2.  If only one
+          value is given it is repeated.  (If more than two numbers are
+          given, the redundant ones are ignored.) 
+
+     pch: the plotting symbols for plotting the data points and the
+          dummy points, in that order; may be either integer or
+          character; defaults to 1:2.  If only one value is given it is
+          repeated.  (If more than two values are given, the redundant
+          ones are ignored.) 
+
+    xlim: the limits on the x-axis.  Defaults to rw[1:2] where rw is
+          the rectangular window specification returned by deldir(). 
+
+    ylim: the limits on the y-axis.  Defaults to rw[3:4] where rw is
+          the rectangular window specification returned by deldir(). 
+
+    xlab: label for the x-axis.  Defaults to 'x'.  Ignored if
+          'add=TRUE'. 
+
+    ylab: label for the y-axis.  Defaults to 'y'.  Ignored if
+          'add=TRUE'. 
+
+     ...: Further plotting parameters to be passed to 'plot()'
+          'segments()' or 'points()'.  Unlikely to be used. 
+
+_D_e_t_a_i_l_s:
+
+     The points in the set being triangulated are plotted with
+     distinguishing symbols.  By default the real points are plotted as
+     circles (pch=1) and the dummy points are plotted as triangles
+     (pch=2).
+
+_S_i_d_e _E_f_f_e_c_t_s:
+
+     A plot of the points being triangulated is produced or added to an
+     existing plot.  As well, the edges of the Delaunay triangles
+     and/or of the Dirichlet tiles are plotted.  By default the
+     triangles are plotted with solid lines (lty=1) and the tiles with
+     dotted lines (lty=2).
+
+_A_u_t_h_o_r(_s):
+
+     Rolf Turner r.turner@auckland.ac.nz <URL:
+     http://www.math.unb.ca/~rolf>
+
+_S_e_e _A_l_s_o:
+
+     'deldir()'
+
+_E_x_a_m_p_l_e_s:
+
+     ## Not run: 
+     try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+     plot(try)
+     #
+     deldir(x,y,list(ndx=4,ndy=4),plot=TRUE,add=TRUE,wl='te',
+            col=c(1,1,2,3,4),num=TRUE)
+     # Plots the tesselation, but does not save the results.
+     try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr',
+                   wp='n')
+     # Plots the triangulation, but not the points, and saves the returned
+     structure.
+     ## End(Not run)
+
@@ -0,0 +1,88 @@
+plot.tile.list            package:deldir            R Documentation
+
+_P_l_o_t _D_i_r_c_h_l_e_t/_V_o_r_o_n_o_i _t_i_l_e_s
+
+_D_e_s_c_r_i_p_t_i_o_n:
+
+     A method for 'plot'.  Plots (sequentially) the tiles associated
+     with each point in the set being tessellated.
+
+_U_s_a_g_e:
+
+     plot.tile.list(x, verbose = FALSE, close=FALSE, pch=1, polycol=NA,
+                    showpoints=TRUE, asp=1, ...)
+     ## S3 method for class 'tile.list':
+     plot(x, verbose = FALSE, close=FALSE, pch=1,
+                              polycol=NA, showpoints=TRUE, asp=1, ...)
+
+_A_r_g_u_m_e_n_t_s:
+
+       x: A list of the tiles in a tessellation, as produced the
+          function 'tile.list()'.
+
+ verbose: Logical scalar; if 'TRUE' the tiles are plotted one at a time
+          (with a ``Go?'' prompt after each) so that the process can be
+          watched.
+
+   close: Logical scalar; if 'TRUE' the outer edges of of the tiles
+          (i.e. the edges of the enclosing rectangle) are drawn. 
+          Otherwise tiles on the periphery of the tessellation are left
+          ``open''.
+
+     pch: The plotting character for plotting the points of the pattern
+          which was tessellated.  Ignored if 'showpoints' is 'FALSE'.
+
+ polycol: Optional vector of integers (or 'NA's); the i-th entry
+          indicates with which colour to fill the i-th tile.  Note that
+          an 'NA' indicates the use of no colour at all.
+
+showpoints: Logical scalar; if 'TRUE' the points of the pattern which
+          was tesselated are plotted.
+
+     asp: The aspect ratio of the plot; integer scalar or 'NA'.  Set
+          this argument equal to 'NA' to allow the data to determine
+          the aspect ratio and hence to make the plot occupy the
+          complete plotting region in both 'x' and 'y' directions. This
+          is inadvisable; see the *Warnings*.
+
+     ...: Optional arguments; not used.  There for consistency with the
+          generic 'plot' function.
+
+_V_a_l_u_e:
+
+     NULL; side effect is a plot.
+
+_W_a_r_n_i_n_g_s:
+
+     The default value for 'verbose' was formerly 'TRUE'; it is now
+     'FALSE'.
+
+     The user is _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' is present ``just because it can be''.
+
+_A_u_t_h_o_r(_s):
+
+     Rolf Turner r.turner@auckland.ac.nz <URL:
+     http://www.math.unb.ca/~rolf>
+
+_S_e_e _A_l_s_o:
+
+     'tile.list()'
+
+_E_x_a_m_p_l_e_s:
+
+       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 rainbox(20).
+       plot(w,polycol=ccc,close=TRUE)
+
@@ -0,0 +1,81 @@
+tile.list               package:deldir               R Documentation
+
+_C_r_e_a_t_e _a _l_i_s_t _o_f _t_i_l_e_s _i_n _a _t_e_s_s_e_l_l_a_t_i_o_n.
+
+_D_e_s_c_r_i_p_t_i_o_n:
+
+     For each point in the set being tessellated produces a list entry
+     describing the Dirichlet/Voronoi tile containing that point.
+
+_U_s_a_g_e:
+
+      tile.list(object) 
+
+_A_r_g_u_m_e_n_t_s:
+
+  object: An object of class 'deldir' as produced by the function
+          'deldir()'.
+
+_V_a_l_u_e:
+
+     A list with one entry for each of the points in the set being
+     tesselated.  Each entry is in turn a list with components
+
+      pt: The coordinates of the point whose tile is being described.
+
+       x: The 'x' coordinates of the vertices of the tile, in
+          anticlockwise order.
+
+       y: The 'y' coordinates of the vertices of the tile, in
+          anticlockwise order.
+
+      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
+
+_A_c_k_n_o_w_l_e_d_g_e_m_e_n_t:
+
+     The author expresses sincere thanks to Majid Yazdani who found and
+     pointed out a serious bug in 'tile.list' in the previous version
+     (0.0-5) of the 'deldir' package.
+
+_W_a_r_n_i_n_g:
+
+     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 the previous version of 'deldir', the corners of
+     the enclosing rectangle are now include 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 when the corresponding polygon is
+     plotted, the result is the intersection of the tile with the
+     enclosing rectangular window.
+
+_A_u_t_h_o_r(_s):
+
+     Rolf Turner r.turner@auckland.ac.nz <URL:
+     http://www.math.unb.ca/~rolf>
+
+_S_e_e _A_l_s_o:
+
+     'deldir()', 'plot.tile.list()'
+
+_E_x_a_m_p_l_e_s:
+
+             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)
+
@@ -0,0 +1,32 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>R: Delaunay Triangulation and Dirichlet (Voronoi) Tessellation.</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" type="text/css" href="../../R.css">
+</head><body>
+<h1>Delaunay Triangulation and Dirichlet (Voronoi) Tessellation. <img class="toplogo" src="../../../doc/html/logo.jpg" alt="[R logo]"></h1>
+
+<hr>
+
+<div align="center">
+<a href="../../../doc/html/packages.html"><img src="../../../doc/html/left.jpg"
+alt="[Package List]" width="30" height="30" border="0"></a>
+<a href="../../../doc/html/index.html"><img src="../../../doc/html/up.jpg"
+alt="[Top]" width="30" height="30" border="0"></a>
+</div>
+
+<h2>Documentation for package &lsquo;deldir&rsquo;</h2>
+
+<h2>Help Pages</h2>
+
+
+<table width="100%">
+<tr><td width="25%"><a href="deldir.html">deldir</a></td>
+<td>Construct the Delaunay triangulation and the Dirichlet (Voronoi) tessellation of a planar point set. </td></tr>
+<tr><td width="25%"><a href="plot.deldir.html">plot.deldir</a></td>
+<td>Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi) tesselation of a planar point set, as constructed by the function deldir. </td></tr>
+<tr><td width="25%"><a href="plot.tile.list.html">plot.tile.list</a></td>
+<td>Plot Dirchlet/Voronoi tiles </td></tr>
+<tr><td width="25%"><a href="tile.list.html">tile.list</a></td>
+<td>Create a list of tiles in a tessellation. </td></tr>
+</table>
+</body></html>
@@ -0,0 +1,41 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>R: Internal deldir functions</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" type="text/css" href="../../R.css">
+</head><body>
+
+<table width="100%" summary="page for deldir-internal {deldir}"><tr><td>deldir-internal {deldir}</td><td align="right">R Documentation</td></tr></table>
+<h2>Internal deldir functions</h2>
+
+
+<h3>Description</h3>
+
+<p>
+Internal deldir functions.
+</p>
+
+
+<h3>Usage</h3>
+
+<pre>
+dumpts(x,y,dpl,rw)
+ind.dup(x,y,rw=NULL,frac=0.0001)
+mid.in(x,y,rx,ry)
+mnnd(x,y)
+get.cnrind(x,y,rw)
+acw(xxx)
+</pre>
+
+
+<h3>Details</h3>
+
+<p>
+These functions are auxilliary and are not intended to be called by
+the user.
+</p>
+
+
+
+<hr><div align="center">[Package <em>deldir</em> version 0.0-7 <a href="00Index.html">Index]</a></div>
+
+</body></html>
@@ -0,0 +1,293 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>R: Construct the Delaunay triangulation and the Dirichlet
+(Voronoi) tessellation of a planar point set.</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" type="text/css" href="../../R.css">
+</head><body>
+
+<table width="100%" summary="page for deldir {deldir}"><tr><td>deldir {deldir}</td><td align="right">R Documentation</td></tr></table>
+<h2>Construct the Delaunay triangulation and the Dirichlet
+(Voronoi) tessellation of a planar point set.</h2>
+
+
+<h3>Description</h3>
+
+<p>
+This function computes the Delaunay triangulation (and hence the
+Dirichlet tesselation) of a planar point set according to the second
+(iterative) algorithm of Lee and Schacter &mdash; see REFERENCES.  The
+triangulation is made to be with respect to the whole plane by
+<code>suspending</code> 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.
+</p>
+
+
+<h3>Usage</h3>
+
+<pre>
+deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, frac=0.0001,
+       sort=TRUE, plotit=FALSE, digits=6, ...)
+</pre>
+
+
+<h3>Arguments</h3>
+
+<table summary="R argblock">
+<tr valign="top"><td><code>x,y</code></td>
+<td>
+The coordinates of the point set being triangulated. These can be
+given by two arguments x and y which are vectors or by a single
+argument x which is a list with components "x" and "y".
+</td></tr>
+<tr valign="top"><td><code>dpl</code></td>
+<td>
+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 auxilliary function dumpts().  The list may have
+components:
+<br>
+ndx: The x-dimension of a rectangular grid; if either ndx or ndy is null,
+no grid is constructed.
+<br>
+ndy: The y-dimension of the aforementioned rectangular grid.
+<br>
+nrad: The number of radii or "spokes", emanating from each data point,
+along which dummy points are to be added.
+<br>
+nper: The number of dummy points per spoke.
+<br>
+fctr: A factor determining the length of each spoke; each spoke is of
+length equal to fctr times the mean nearest neighbour distance of the data.
+(This distance is calculated by the auxilliary function mnnd().)
+<br>
+x: A vector of x-coordinates of "ad hoc" dummy points
+<br>
+y: A vector of the corresponding y-coordinates of "ad hoc" dummy points
+<br>
+</td></tr>
+<tr valign="top"><td><code>rw</code></td>
+<td>
+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.
+</td></tr>
+<tr valign="top"><td><code>eps</code></td>
+<td>
+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.
+</td></tr>
+<tr valign="top"><td><code>frac</code></td>
+<td>
+A value specifying the tolerance used in eliminating duplicate
+points; defaults to 0.0001. Points are considered duplicates if
+abs(x1-x2) &lt; frac*(xmax-xmin) AND abs(y1-y2) &lt; frac*(ymax-ymin).
+</td></tr>
+<tr valign="top"><td><code>sort</code></td>
+<td>
+Logical argument; if <code>TRUE</code> (the default) the data (including dummy
+points) are sorted into a sequence of "bins" prior to triangulation;
+this makes the algorithm slightly more efficient.  Normally one would
+set sort equal to <code>FALSE</code> only if one wished to observe some of the
+fine detail of the way in which adding a point to a data set affected
+the triangulation, and therefore wished to make sure that the point
+in question was added last.  Essentially this argument would get used
+only in a de-bugging process.
+</td></tr>
+<tr valign="top"><td><code>plotit</code></td>
+<td>
+Logical argument; if <code>TRUE</code> a plot of the triangulation and tessellation
+is produced; the default is <code>FALSE</code>.
+</td></tr>
+<tr valign="top"><td><code>digits</code></td>
+<td>
+The number of decimal places to which all numeric values in the
+returned list should be rounded.  Defaults to 6.
+</td></tr>
+<tr valign="top"><td><code>...</code></td>
+<td>
+Auxilliary arguments add, wlines, wpoints, number, nex, col, lty,
+pch, xlim, and ylim (and possibly other plotting parameters) may be
+passed to plot.deldir through "..." if plotit=<code>TRUE</code>.
+</td></tr>
+</table>
+
+<h3>Details</h3>
+
+<p>
+This package is a (straightforward) adaptation of the Splus library
+section ``delaunay'' to R.  That library section is an implementation
+of the Lee-Schacter algorithm, which was originally written as a
+stand-alone Fortran program in 1987/88 by Rolf Turner, while with the
+Division of Mathematics and Statistics, CSIRO, Sydney, Australia.  It
+was re-written as an Splus function (using dynamically loaded Fortran
+code), by Rolf Turner while visiting the University of Western
+Australia, May, 1995.
+</p>
+<p>
+Further revisions 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 the University of Western Australia).  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.
+</p>
+<p>
+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, is now deprecated in
+favour of the current package ``delaunay'' package.  Don also
+collaborated in the preparation of that package.
+</p>
+<p>
+Further revisions and bug-fixes were made in 1998, 1999, and 2002.
+</p>
+
+
+<h3>Value</h3>
+
+<p>
+A list (of class <code>deldir</code>), invisible if plotit=<code>TRUE</code>, with components:
+</p>
+<table summary="R argblock">
+<tr valign="top"><td><code>delsgs</code></td>
+<td>
+a matrix 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 (x1,y1,x2,y2).  The last two entries are the
+indices of the two points which are joined.
+</td></tr>
+<tr valign="top"><td><code>dirsgs</code></td>
+<td>
+a data frame with 8 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 (x1,y1,x2,y2).  The fifth and sixth entries are the indices
+of the two points, in the set being triangulated, which are separated
+by that edge. The seventh and eighth entries are logical values.  The
+seventh indicates whether the first endpoint of the corresponding
+edge of a Dirichlet tile is a boundary point (a point on the boundary
+of the rectangular window).  Likewise for the eighth entry and the
+second endpoint of the edge.
+</td></tr>
+<tr valign="top"><td><code>summary</code></td>
+<td>
+a matrix with 9 columns, and (n.data + n.dum) rows (see below).
+These rows correspond to the points in the set being triangulated.
+The columns are named "x" (the x-coordinate of the point), "y" (the
+y-coordinate), "n.tri" (the number of Delaunay triangles emanating
+from the point), "del.area" (1/3 of the total area of all the
+Delaunay triangles emanating from the point), "del.wts" (the
+corresponding entry of the "del.area" column divided by the sum of
+this column), "n.tside" (the number of sides &mdash; within the
+rectangular window &mdash; of the Dirichlet tile surrounding the point),
+"nbpt" (the number of points in which the Dirichlet tile intersects
+the boundary of the rectangular window), "dir.area" (the area of the
+Dirichlet tile surrounding the point), and "dir.wts" (the
+corresponding entry of the "dir.area" column divided by the sum of
+this column).  Note that the factor of 1/3 associated with the
+del.area column arises because each triangle occurs three times &mdash;
+once for each corner.
+</td></tr>
+<tr valign="top"><td><code>n.data</code></td>
+<td>
+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 "summary" correspond to real points.
+</td></tr>
+<tr valign="top"><td><code>n.dum</code></td>
+<td>
+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 "summary" correspond to dummy
+points.
+</td></tr>
+<tr valign="top"><td><code>del.area</code></td>
+<td>
+the area of the convex hull of the set of points being triangulated,
+as formed by summing the "del.area" column of "summary".
+</td></tr>
+<tr valign="top"><td><code>dir.area</code></td>
+<td>
+the area of the rectangular window enclosing the points being triangulated,
+as formed by summing the "dir.area" column of "summary".
+</td></tr>
+<tr valign="top"><td><code>rw</code></td>
+<td>
+the specification of the corners of the rectangular window enclosing
+the data, in the order (xmin, xmax, ymin, ymax).
+</td></tr>
+</table>
+
+<h3>Side Effects</h3>
+
+<p>
+If plotit==<code>TRUE</code> a plot of the triangulation and/or tessellation is produced
+or added to an existing plot.
+</p>
+
+
+<h3>Note</h3>
+
+<p>
+If ndx &gt;= 2 and ndy &gt;= 2, then the rectangular window IS the convex
+hull, and so the values of del.area and dir.area are identical.
+</p>
+
+
+<h3>Author(s)</h3>
+
+<p>
+Rolf Turner
+<a href="mailto:r.turner@auckland.ac.nz">r.turner@auckland.ac.nz</a>
+<a href="http://www.math.unb.ca/~rolf">http://www.math.unb.ca/~rolf</a>
+</p>
+
+
+<h3>References</h3>
+
+<p>
+Lee, D. T., and Schacter, B. J.  "Two algorithms for constructing a
+Delaunay triangulation", Int. J. Computer and Information
+Sciences, Vol. 9, No. 3, 1980, pp. 219 &ndash; 242.
+</p>
+<p>
+Ahuja, N. and Schacter, B. J. (1983). Pattern Models.  New York: Wiley.
+</p>
+
+
+<h3>See Also</h3>
+
+<p>
+plot.deldir
+</p>
+
+
+<h3>Examples</h3>
+
+<pre>
+x   &lt;- c(2.3,3.0,7.0,1.0,3.0,8.0)
+y   &lt;- c(2.3,3.0,2.0,5.0,8.0,9.0)
+try &lt;- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+# Puts dummy points at the corners of the rectangular
+# window, i.e. at (0,0), (10,0), (10,10), and (0,10)
+## Not run: 
+try &lt;- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr')
+## End(Not run)
+# Plots the triangulation which was created (but not the tesselation).
+</pre>
+
+
+
+<hr><div align="center">[Package <em>deldir</em> version 0.0-7 <a href="00Index.html">Index]</a></div>
+
+</body></html>
@@ -0,0 +1,181 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>R: Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi)
+tesselation of a planar point set, as constructed by the function deldir.</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" type="text/css" href="../../R.css">
+</head><body>
+
+<table width="100%" summary="page for plot.deldir {deldir}"><tr><td>plot.deldir {deldir}</td><td align="right">R Documentation</td></tr></table>
+<h2>Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi)
+tesselation of a planar point set, as constructed by the function deldir.</h2>
+
+
+<h3>Description</h3>
+
+<p>
+This is a method for plot.
+</p>
+
+
+<h3>Usage</h3>
+
+<pre>
+## S3 method for class 'deldir':
+plot(x,add=FALSE,wlines=c('both','triang','tess'),
+                      wpoints=c('both','real','dummy','none'),
+                      number=FALSE,cex=1,nex=1,col=NULL,lty=NULL,
+                      pch=NULL,xlim=NULL,ylim=NULL,xlab='x',ylab='y',...)
+
+</pre>
+
+
+<h3>Arguments</h3>
+
+<table summary="R argblock">
+<tr valign="top"><td><code>x</code></td>
+<td>
+An object of class "deldir" as constructed by the function deldir.
+</td></tr>
+<tr valign="top"><td><code>add</code></td>
+<td>
+logical argument; should the plot be added to an existing plot?
+</td></tr>
+<tr valign="top"><td><code>wlines</code></td>
+<td>
+"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) ?
+</td></tr>
+<tr valign="top"><td><code>wpoints</code></td>
+<td>
+"which points?".  I.e.  should the real points be plotted
+(wpoints='real'), should the dummy points be plotted
+(wpoints='dummy'), should both be plotted (wpoints='both', the
+default) or should no points be plotted (wpoints='none')?
+</td></tr>
+<tr valign="top"><td><code>number</code></td>
+<td>
+Logical argument, defaulting to <code>FALSE</code>; if <code>TRUE</code> then the
+points plotted will be labelled with their index numbers
+(corresponding to the row numbers of the matrix "summary" in the
+output of deldir).
+</td></tr>
+<tr valign="top"><td><code>cex</code></td>
+<td>
+The value of the character expansion argument cex to be used
+with the plotting symbols for plotting the points.
+</td></tr>
+<tr valign="top"><td><code>nex</code></td>
+<td>
+The value of the character expansion argument cex to be used by the
+text function when numbering the points with their indices.  Used only
+if number=<code>TRUE</code>.
+</td></tr>
+<tr valign="top"><td><code>col</code></td>
+<td>
+the colour numbers for plotting the triangulation, the tesselation,
+the data points, the dummy points, and the point numbers, in that
+order; defaults to c(1,1,1,1,1).  If fewer than five numbers are
+given, they are recycled.  (If more than five numbers are given, the
+redundant ones are ignored.)
+</td></tr>
+<tr valign="top"><td><code>lty</code></td>
+<td>
+the line type numbers for plotting the triangulation and the
+tesselation, in that order; defaults to 1:2.  If only one value is
+given it is repeated.  (If more than two numbers are given, the
+redundant ones are ignored.)
+</td></tr>
+<tr valign="top"><td><code>pch</code></td>
+<td>
+the plotting symbols for plotting the data points and the dummy
+points, in that order; may be either integer or character; defaults
+to 1:2.  If only one value is given it is repeated.  (If more than
+two values are given, the redundant ones are ignored.)
+</td></tr>
+<tr valign="top"><td><code>xlim</code></td>
+<td>
+the limits on the x-axis.  Defaults to rw[1:2] where rw is the
+rectangular window specification returned by deldir().
+</td></tr>
+<tr valign="top"><td><code>ylim</code></td>
+<td>
+the limits on the y-axis.  Defaults to rw[3:4] where rw is the
+rectangular window specification returned by deldir().
+</td></tr>
+<tr valign="top"><td><code>xlab</code></td>
+<td>
+label for the x-axis.  Defaults to <code>x</code>.  Ignored if
+<code>add=TRUE</code>.
+</td></tr>
+<tr valign="top"><td><code>ylab</code></td>
+<td>
+label for the y-axis.  Defaults to <code>y</code>.  Ignored if
+<code>add=TRUE</code>.
+</td></tr>
+<tr valign="top"><td><code>...</code></td>
+<td>
+Further plotting parameters to be passed to <code>plot()</code>
+<code>segments()</code> or <code>points()</code>.  Unlikely to be used.
+</td></tr>
+</table>
+
+<h3>Details</h3>
+
+<p>
+The points in the set being triangulated are plotted with distinguishing
+symbols.  By default the real points are plotted as circles (pch=1) and the
+dummy points are plotted as triangles (pch=2).
+</p>
+
+
+<h3>Side Effects</h3>
+
+<p>
+A plot of the points being triangulated is produced or added to
+an existing plot.  As well, the edges of the Delaunay
+triangles and/or of the Dirichlet tiles are plotted.  By default
+the triangles are plotted with solid lines (lty=1) and the tiles
+with dotted lines (lty=2).
+</p>
+
+
+<h3>Author(s)</h3>
+
+<p>
+Rolf Turner
+<a href="mailto:r.turner@auckland.ac.nz">r.turner@auckland.ac.nz</a>
+<a href="http://www.math.unb.ca/~rolf">http://www.math.unb.ca/~rolf</a>
+</p>
+
+
+<h3>See Also</h3>
+
+<p>
+<code><a href="deldir.html">deldir</a>()</code>
+</p>
+
+
+<h3>Examples</h3>
+
+<pre>
+## Not run: 
+try &lt;- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+plot(try)
+#
+deldir(x,y,list(ndx=4,ndy=4),plot=TRUE,add=TRUE,wl='te',
+       col=c(1,1,2,3,4),num=TRUE)
+# Plots the tesselation, but does not save the results.
+try &lt;- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr',
+              wp='n')
+# Plots the triangulation, but not the points, and saves the returned
+structure.
+## End(Not run)
+</pre>
+
+
+
+<hr><div align="center">[Package <em>deldir</em> version 0.0-7 <a href="00Index.html">Index]</a></div>
+
+</body></html>
@@ -0,0 +1,133 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>R: Plot Dirchlet/Voronoi tiles</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" type="text/css" href="../../R.css">
+</head><body>
+
+<table width="100%" summary="page for plot.tile.list {deldir}"><tr><td>plot.tile.list {deldir}</td><td align="right">R Documentation</td></tr></table>
+<h2>Plot Dirchlet/Voronoi tiles</h2>
+
+
+<h3>Description</h3>
+
+<p>
+A method for <code>plot</code>.  Plots (sequentially)
+the tiles associated with each point in the set being tessellated.
+</p>
+
+
+<h3>Usage</h3>
+
+<pre>
+plot.tile.list(x, verbose = FALSE, close=FALSE, pch=1, polycol=NA,
+               showpoints=TRUE, asp=1, ...)
+## S3 method for class 'tile.list':
+plot(x, verbose = FALSE, close=FALSE, pch=1,
+                         polycol=NA, showpoints=TRUE, asp=1, ...)
+</pre>
+
+
+<h3>Arguments</h3>
+
+<table summary="R argblock">
+<tr valign="top"><td><code>x</code></td>
+<td>
+A list of the tiles in a tessellation, as produced
+the function <code><a href="tile.list.html">tile.list</a>()</code>.</td></tr>
+<tr valign="top"><td><code>verbose</code></td>
+<td>
+Logical scalar; if <code>TRUE</code> the tiles are
+plotted one at a time (with a ``Go?'' prompt after each)
+so that the process can be watched.</td></tr>
+<tr valign="top"><td><code>close</code></td>
+<td>
+Logical scalar; if <code>TRUE</code> the outer edges of
+of the tiles (i.e. the edges of the enclosing rectangle)
+are drawn.  Otherwise tiles on the periphery of the
+tessellation are left ``open''.</td></tr>
+<tr valign="top"><td><code>pch</code></td>
+<td>
+The plotting character for plotting the points of the
+pattern which was tessellated.  Ignored if <code>showpoints</code>
+is <code>FALSE</code>.</td></tr>
+<tr valign="top"><td><code>polycol</code></td>
+<td>
+Optional vector of integers (or <code>NA</code>s);
+the <i>i</i>-th entry indicates with which colour to fill
+the <i>i</i>-th tile.  Note that an <code>NA</code> indicates
+the use of no colour at all.</td></tr>
+<tr valign="top"><td><code>showpoints</code></td>
+<td>
+Logical scalar; if <code>TRUE</code> the points of
+the pattern which was tesselated are plotted.</td></tr>
+<tr valign="top"><td><code>asp</code></td>
+<td>
+The aspect ratio of the plot; integer scalar or
+<code>NA</code>.  Set this argument equal to <code>NA</code> to allow the data
+to determine the aspect ratio and hence to make the plot occupy the
+complete plotting region in both <code>x</code> and <code>y</code> directions.
+This is inadvisable; see the <B>Warnings</B>.</td></tr>
+<tr valign="top"><td><code>...</code></td>
+<td>
+Optional arguments; not used.  There for consistency
+with the generic <code>plot</code> function.</td></tr>
+</table>
+
+<h3>Value</h3>
+
+<p>
+NULL; side effect is a plot.</p>
+
+<h3>Warnings</h3>
+
+<p>
+The default value for <code>verbose</code> was formerly <code>TRUE</code>;
+it is now <code>FALSE</code>.
+</p>
+<p>
+The user is <EM>strongly advised</EM> not to set the value of
+<code>asp</code> but rather to leave <code>asp</code> equal to its default
+value of <code>1</code>.  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.)
+</p>
+<p>
+The argument <code>asp</code> is present ``just because it can be''.
+</p>
+
+
+<h3>Author(s)</h3>
+
+<p>
+Rolf Turner
+<a href="mailto:r.turner@auckland.ac.nz">r.turner@auckland.ac.nz</a>
+<a href="http://www.math.unb.ca/~rolf">http://www.math.unb.ca/~rolf</a>
+</p>
+
+
+<h3>See Also</h3>
+
+<p>
+<code><a href="tile.list.html">tile.list</a>()</code>
+</p>
+
+
+<h3>Examples</h3>
+
+<pre>
+  x &lt;- runif(20)
+  y &lt;- runif(20)
+  z &lt;- deldir(x,y,rw=c(0,1,0,1))
+  w &lt;- tile.list(z)
+  plot(w)
+  ccc &lt;- heat.colors(20) # Or topo.colors(20), or terrain.colors(20)
+                         # or cm.colors(20), or rainbox(20).
+  plot(w,polycol=ccc,close=TRUE)
+</pre>
+
+
+
+<hr><div align="center">[Package <em>deldir</em> version 0.0-7 <a href="00Index.html">Index]</a></div>
+
+</body></html>
@@ -0,0 +1,127 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>R: Create a list of tiles in a tessellation.</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" type="text/css" href="../../R.css">
+</head><body>
+
+<table width="100%" summary="page for tile.list {deldir}"><tr><td>tile.list {deldir}</td><td align="right">R Documentation</td></tr></table>
+<h2>Create a list of tiles in a tessellation.</h2>
+
+
+<h3>Description</h3>
+
+<p>
+For each point in the set being tessellated produces a list
+entry describing the Dirichlet/Voronoi tile containing that
+point.
+</p>
+
+
+<h3>Usage</h3>
+
+<pre> tile.list(object) </pre>
+
+
+<h3>Arguments</h3>
+
+<table summary="R argblock">
+<tr valign="top"><td><code>object</code></td>
+<td>
+An object of class <code>deldir</code> as produced
+by the function <code><a href="deldir.html">deldir</a>()</code>.</td></tr>
+</table>
+
+<h3>Value</h3>
+
+<p>
+A list with one entry for each of the points in the set
+being tesselated.  Each entry is in turn a list with
+components
+</p>
+<table summary="R argblock">
+<tr valign="top"><td><code>pt</code></td>
+<td>
+The coordinates of the point whose tile is being described.</td></tr>
+<tr valign="top"><td><code>x</code></td>
+<td>
+The <code>x</code> coordinates of the vertices of the tile, in
+anticlockwise order.</td></tr>
+<tr valign="top"><td><code>y</code></td>
+<td>
+The <code>y</code> coordinates of the vertices of the tile, in
+anticlockwise order.</td></tr>
+<tr valign="top"><td><code>bp</code></td>
+<td>
+Vector of logicals indicating whether the tile vertex is a
+``real'' vertex, or a <EM>boundary point</EM>, i.e. a point where the
+tile edge intersects the boundary of the enclosing rectangle</td></tr>
+</table>
+
+<h3>Acknowledgement</h3>
+
+<p>
+The author expresses sincere thanks to Majid Yazdani who found and
+pointed out a serious bug in <code>tile.list</code> in the previous
+version (0.0-5) of the <code>deldir</code> package.
+</p>
+
+
+<h3>Warning</h3>
+
+<p>
+The set of vertices of each tile may be ``incomplete''.  Only
+vertices which lie within the enclosing rectangle, and ``boundary
+points'' are listed.
+</p>
+<p>
+Note that the enclosing rectangle may be specified by the user
+in the call to <code><a href="deldir.html">deldir</a>()</code>.
+</p>
+<p>
+In contrast to the previous version of <code>deldir</code>, the
+corners of the enclosing rectangle are now include 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 when the corresponding
+polygon is plotted, the result is the intersection of the tile
+with the enclosing rectangular window.
+</p>
+
+
+<h3>Author(s)</h3>
+
+<p>
+Rolf Turner
+<a href="mailto:r.turner@auckland.ac.nz">r.turner@auckland.ac.nz</a>
+<a href="http://www.math.unb.ca/~rolf">http://www.math.unb.ca/~rolf</a>
+</p>
+
+
+<h3>See Also</h3>
+
+<p>
+<code><a href="deldir.html">deldir</a>()</code>, <code><a href="plot.tile.list.html">plot.tile.list</a>()</code>
+</p>
+
+
+<h3>Examples</h3>
+
+<pre>
+        x &lt;- runif(20)
+        y &lt;- runif(20)
+        z &lt;- deldir(x,y)
+        w &lt;- tile.list(z)
+
+        z &lt;- deldir(x,y,rw=c(0,1,0,1))
+        w &lt;- tile.list(z)
+
+        z &lt;- deldir(x,y,rw=c(0,1,0,1),dpl=list(ndx=2,ndy=2))
+        w &lt;- tile.list(z)
+
+</pre>
+
+
+
+<hr><div align="center">[Package <em>deldir</em> version 0.0-7 <a href="00Index.html">Index]</a></div>
+
+</body></html>
@@ -0,0 +1,26 @@
+\HeaderA{deldir-internal}{Internal deldir functions}{deldir.Rdash.internal}
+\aliasA{acw}{deldir-internal}{acw}
+\aliasA{dumpts}{deldir-internal}{dumpts}
+\aliasA{get.cnrind}{deldir-internal}{get.cnrind}
+\aliasA{ind.dup}{deldir-internal}{ind.dup}
+\aliasA{mid.in}{deldir-internal}{mid.in}
+\aliasA{mnnd}{deldir-internal}{mnnd}
+\keyword{internal}{deldir-internal}
+\begin{Description}\relax
+Internal deldir functions.
+\end{Description}
+\begin{Usage}
+\begin{verbatim}
+dumpts(x,y,dpl,rw)
+ind.dup(x,y,rw=NULL,frac=0.0001)
+mid.in(x,y,rx,ry)
+mnnd(x,y)
+get.cnrind(x,y,rw)
+acw(xxx)
+\end{verbatim}
+\end{Usage}
+\begin{Details}\relax
+These functions are auxilliary and are not intended to be called by
+the user.
+\end{Details}
+
@@ -0,0 +1,216 @@
+\HeaderA{deldir}{Construct the Delaunay triangulation and the Dirichlet
+(Voronoi) tessellation of a planar point set.}{deldir}
+\keyword{spatial}{deldir}
+\begin{Description}\relax
+This function computes the Delaunay triangulation (and hence the
+Dirichlet 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
+\code{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.
+\end{Description}
+\begin{Usage}
+\begin{verbatim}
+deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, frac=0.0001,
+       sort=TRUE, plotit=FALSE, digits=6, ...)
+\end{verbatim}
+\end{Usage}
+\begin{Arguments}
+\begin{ldescription}
+\item[\code{x,y}] The coordinates of the point set being triangulated. These can be
+given by two arguments x and y which are vectors or by a single
+argument x which is a list with components "x" and "y".
+
+\item[\code{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 auxilliary function dumpts().  The list may have
+components:
+
+
+ndx: The x-dimension of a rectangular grid; if either ndx or ndy is null,
+no grid is constructed.
+
+
+ndy: The y-dimension of the aforementioned rectangular grid.
+
+
+nrad: The number of radii or "spokes", emanating from each data point,
+along which dummy points are to be added.
+
+
+nper: The number of dummy points per spoke.
+
+
+fctr: A factor determining the length of each spoke; each spoke is of
+length equal to fctr times the mean nearest neighbour distance of the data.
+(This distance is calculated by the auxilliary function mnnd().)
+
+
+x: A vector of x-coordinates of "ad hoc" dummy points
+
+
+y: A vector of the corresponding y-coordinates of "ad hoc" dummy points
+
+
+
+\item[\code{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.
+
+\item[\code{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.
+
+\item[\code{frac}] A value specifying the tolerance used in eliminating duplicate
+points; defaults to 0.0001. Points are considered duplicates if
+abs(x1-x2) < frac*(xmax-xmin) AND abs(y1-y2) < frac*(ymax-ymin).
+
+\item[\code{sort}] Logical argument; if \code{TRUE} (the default) the data (including dummy
+points) are sorted into a sequence of "bins" prior to triangulation;
+this makes the algorithm slightly more efficient.  Normally one would
+set sort equal to \code{FALSE} only if one wished to observe some of the
+fine detail of the way in which adding a point to a data set affected
+the triangulation, and therefore wished to make sure that the point
+in question was added last.  Essentially this argument would get used
+only in a de-bugging process.
+
+\item[\code{plotit}] Logical argument; if \code{TRUE} a plot of the triangulation and tessellation
+is produced; the default is \code{FALSE}.
+
+\item[\code{digits}] The number of decimal places to which all numeric values in the
+returned list should be rounded.  Defaults to 6.
+
+\item[\code{...}] Auxilliary arguments add, wlines, wpoints, number, nex, col, lty,
+pch, xlim, and ylim (and possibly other plotting parameters) may be
+passed to plot.deldir through "\dots" if plotit=\code{TRUE}.
+
+\end{ldescription}
+\end{Arguments}
+\begin{Details}\relax
+This package is a (straightforward) adaptation of the Splus library
+section ``delaunay'' to R.  That library section is an implementation
+of the Lee-Schacter algorithm, which was originally written as a
+stand-alone Fortran program in 1987/88 by Rolf Turner, while with the
+Division of Mathematics and Statistics, CSIRO, Sydney, Australia.  It
+was re-written as an Splus function (using dynamically loaded Fortran
+code), by Rolf Turner while visiting the University of Western
+Australia, May, 1995.
+
+Further revisions 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 the University of Western Australia).  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, is now deprecated in
+favour of the current package ``delaunay'' package.  Don also
+collaborated in the preparation of that package.
+
+Further revisions and bug-fixes were made in 1998, 1999, and 2002.
+\end{Details}
+\begin{Value}
+A list (of class \code{deldir}), invisible if plotit=\code{TRUE}, with components:
+
+\begin{ldescription}
+\item[\code{delsgs}] a matrix 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 (x1,y1,x2,y2).  The last two entries are the
+indices of the two points which are joined.
+
+\item[\code{dirsgs}] a data frame with 8 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 (x1,y1,x2,y2).  The fifth and sixth entries are the indices
+of the two points, in the set being triangulated, which are separated
+by that edge. The seventh and eighth entries are logical values.  The
+seventh indicates whether the first endpoint of the corresponding
+edge of a Dirichlet tile is a boundary point (a point on the boundary
+of the rectangular window).  Likewise for the eighth entry and the
+second endpoint of the edge.
+
+\item[\code{summary}] a matrix with 9 columns, and (n.data + n.dum) rows (see below).
+These rows correspond to the points in the set being triangulated.
+The columns are named "x" (the x-coordinate of the point), "y" (the
+y-coordinate), "n.tri" (the number of Delaunay triangles emanating
+from the point), "del.area" (1/3 of the total area of all the
+Delaunay triangles emanating from the point), "del.wts" (the
+corresponding entry of the "del.area" column divided by the sum of
+this column), "n.tside" (the number of sides --- within the
+rectangular window --- of the Dirichlet tile surrounding the point),
+"nbpt" (the number of points in which the Dirichlet tile intersects
+the boundary of the rectangular window), "dir.area" (the area of the
+Dirichlet tile surrounding the point), and "dir.wts" (the
+corresponding entry of the "dir.area" column divided by the sum of
+this column).  Note that the factor of 1/3 associated with the
+del.area column arises because each triangle occurs three times ---
+once for each corner.
+
+\item[\code{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 "summary" correspond to real points.
+
+\item[\code{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 "summary" correspond to dummy
+points.
+
+\item[\code{del.area}] the area of the convex hull of the set of points being triangulated,
+as formed by summing the "del.area" column of "summary".
+
+\item[\code{dir.area}] the area of the rectangular window enclosing the points being triangulated,
+as formed by summing the "dir.area" column of "summary".
+
+\item[\code{rw}] the specification of the corners of the rectangular window enclosing
+the data, in the order (xmin, xmax, ymin, ymax).
+
+\end{ldescription}
+\end{Value}
+\begin{Section}{Side Effects}
+If plotit==\code{TRUE} a plot of the triangulation and/or tessellation is produced
+or added to an existing plot.
+\end{Section}
+\begin{Note}\relax
+If ndx >= 2 and ndy >= 2, then the rectangular window IS the convex
+hull, and so the values of del.area and dir.area are identical.
+\end{Note}
+\begin{Author}\relax
+Rolf Turner
+\email{r.turner@auckland.ac.nz}
+\url{http://www.math.unb.ca/~rolf}
+\end{Author}
+\begin{References}\relax
+Lee, D. T., and Schacter, B. J.  "Two algorithms for constructing a
+Delaunay triangulation", Int. J. Computer and Information
+Sciences, Vol. 9, No. 3, 1980, pp. 219 -- 242.
+
+Ahuja, N. and Schacter, B. J. (1983). Pattern Models.  New York: Wiley.
+\end{References}
+\begin{SeeAlso}\relax
+plot.deldir
+\end{SeeAlso}
+\begin{Examples}
+\begin{ExampleCode}
+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)
+try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+# Puts dummy points at the corners of the rectangular
+# window, i.e. at (0,0), (10,0), (10,10), and (0,10)
+## Not run: 
+try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr')
+## End(Not run)
+# Plots the triangulation which was created (but not the tesselation).
+\end{ExampleCode}
+\end{Examples}
+
@@ -0,0 +1,113 @@
+\HeaderA{plot.deldir}{Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi)
+tesselation of a planar point set, as constructed by the function deldir.}{plot.deldir}
+\begin{Description}\relax
+This is a method for plot.
+\end{Description}
+\begin{Usage}
+\begin{verbatim}
+## S3 method for class 'deldir':
+plot(x,add=FALSE,wlines=c('both','triang','tess'),
+                      wpoints=c('both','real','dummy','none'),
+                      number=FALSE,cex=1,nex=1,col=NULL,lty=NULL,
+                      pch=NULL,xlim=NULL,ylim=NULL,xlab='x',ylab='y',...)
+
+\end{verbatim}
+\end{Usage}
+\begin{Arguments}
+\begin{ldescription}
+\item[\code{x}] An object of class "deldir" as constructed by the function deldir.
+
+\item[\code{add}] logical argument; should the plot be added to an existing plot?
+
+\item[\code{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) ?
+
+\item[\code{wpoints}] "which points?".  I.e.  should the real points be plotted
+(wpoints='real'), should the dummy points be plotted
+(wpoints='dummy'), should both be plotted (wpoints='both', the
+default) or should no points be plotted (wpoints='none')?
+
+\item[\code{number}] Logical argument, defaulting to \code{FALSE}; if \code{TRUE} then the
+points plotted will be labelled with their index numbers
+(corresponding to the row numbers of the matrix "summary" in the
+output of deldir).
+
+\item[\code{cex}] The value of the character expansion argument cex to be used
+with the plotting symbols for plotting the points.
+
+\item[\code{nex}] The value of the character expansion argument cex to be used by the
+text function when numbering the points with their indices.  Used only
+if number=\code{TRUE}.
+
+\item[\code{col}] the colour numbers for plotting the triangulation, the tesselation,
+the data points, the dummy points, and the point numbers, in that
+order; defaults to c(1,1,1,1,1).  If fewer than five numbers are
+given, they are recycled.  (If more than five numbers are given, the
+redundant ones are ignored.)
+
+\item[\code{lty}] the line type numbers for plotting the triangulation and the
+tesselation, in that order; defaults to 1:2.  If only one value is
+given it is repeated.  (If more than two numbers are given, the
+redundant ones are ignored.)
+
+\item[\code{pch}] the plotting symbols for plotting the data points and the dummy
+points, in that order; may be either integer or character; defaults
+to 1:2.  If only one value is given it is repeated.  (If more than
+two values are given, the redundant ones are ignored.)
+
+\item[\code{xlim}] the limits on the x-axis.  Defaults to rw[1:2] where rw is the
+rectangular window specification returned by deldir().
+
+\item[\code{ylim}] the limits on the y-axis.  Defaults to rw[3:4] where rw is the
+rectangular window specification returned by deldir().
+
+\item[\code{xlab}] label for the x-axis.  Defaults to \code{x}.  Ignored if
+\code{add=TRUE}.
+
+\item[\code{ylab}] label for the y-axis.  Defaults to \code{y}.  Ignored if
+\code{add=TRUE}.
+
+\item[\code{...}] Further plotting parameters to be passed to \code{plot()}
+\code{segments()} or \code{points()}.  Unlikely to be used.
+
+\end{ldescription}
+\end{Arguments}
+\begin{Details}\relax
+The points in the set being triangulated are plotted with distinguishing
+symbols.  By default the real points are plotted as circles (pch=1) and the
+dummy points are plotted as triangles (pch=2).
+\end{Details}
+\begin{Section}{Side Effects}
+A plot of the points being triangulated is produced or added to
+an existing plot.  As well, the edges of the Delaunay
+triangles and/or of the Dirichlet tiles are plotted.  By default
+the triangles are plotted with solid lines (lty=1) and the tiles
+with dotted lines (lty=2).
+\end{Section}
+\begin{Author}\relax
+Rolf Turner
+\email{r.turner@auckland.ac.nz}
+\url{http://www.math.unb.ca/~rolf}
+\end{Author}
+\begin{SeeAlso}\relax
+\code{\LinkA{deldir}{deldir}()}
+\end{SeeAlso}
+\begin{Examples}
+\begin{ExampleCode}
+## Not run: 
+try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
+plot(try)
+#
+deldir(x,y,list(ndx=4,ndy=4),plot=TRUE,add=TRUE,wl='te',
+       col=c(1,1,2,3,4),num=TRUE)
+# Plots the tesselation, but does not save the results.
+try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr',
+              wp='n')
+# Plots the triangulation, but not the points, and saves the returned
+structure.
+## End(Not run)
+\end{ExampleCode}
+\end{Examples}
+
@@ -0,0 +1,81 @@
+\HeaderA{plot.tile.list}{Plot Dirchlet/Voronoi tiles}{plot.tile.list}
+\keyword{hplot}{plot.tile.list}
+\begin{Description}\relax
+A method for \code{plot}.  Plots (sequentially)
+the tiles associated with each point in the set being tessellated.
+\end{Description}
+\begin{Usage}
+\begin{verbatim}
+plot.tile.list(x, verbose = FALSE, close=FALSE, pch=1, polycol=NA,
+               showpoints=TRUE, asp=1, ...)
+## S3 method for class 'tile.list':
+plot(x, verbose = FALSE, close=FALSE, pch=1,
+                         polycol=NA, showpoints=TRUE, asp=1, ...)
+\end{verbatim}
+\end{Usage}
+\begin{Arguments}
+\begin{ldescription}
+\item[\code{x}] A list of the tiles in a tessellation, as produced
+the function \code{\LinkA{tile.list}{tile.list}()}.
+\item[\code{verbose}] Logical scalar; if \code{TRUE} the tiles are
+plotted one at a time (with a ``Go?'' prompt after each)
+so that the process can be watched.
+\item[\code{close}] Logical scalar; if \code{TRUE} the outer edges of
+of the tiles (i.e. the edges of the enclosing rectangle)
+are drawn.  Otherwise tiles on the periphery of the
+tessellation are left ``open''.
+\item[\code{pch}] The plotting character for plotting the points of the
+pattern which was tessellated.  Ignored if \code{showpoints}
+is \code{FALSE}.
+\item[\code{polycol}] Optional vector of integers (or \code{NA}s);
+the \eqn{i}{}-th entry indicates with which colour to fill
+the \eqn{i}{}-th tile.  Note that an \code{NA} indicates
+the use of no colour at all.
+\item[\code{showpoints}] Logical scalar; if \code{TRUE} the points of
+the pattern which was tesselated are plotted.
+\item[\code{asp}] The aspect ratio of the plot; integer scalar or
+\code{NA}.  Set this argument equal to \code{NA} to allow the data
+to determine the aspect ratio and hence to make the plot occupy the
+complete plotting region in both \code{x} and \code{y} directions.
+This is inadvisable; see the \bold{Warnings}.
+\item[\code{...}] Optional arguments; not used.  There for consistency
+with the generic \code{plot} function.
+\end{ldescription}
+\end{Arguments}
+\begin{Value}
+NULL; side effect is a plot.
+\end{Value}
+\begin{Section}{Warnings}
+The default value for \code{verbose} was formerly \code{TRUE};
+it is now \code{FALSE}.
+
+The user is \emph{strongly advised} not to set the value of
+\code{asp} but rather to leave \code{asp} equal to its default
+value of \code{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 \code{asp} is present ``just because it can be''.
+\end{Section}
+\begin{Author}\relax
+Rolf Turner
+\email{r.turner@auckland.ac.nz}
+\url{http://www.math.unb.ca/~rolf}
+\end{Author}
+\begin{SeeAlso}\relax
+\code{\LinkA{tile.list}{tile.list}()}
+\end{SeeAlso}
+\begin{Examples}
+\begin{ExampleCode}
+  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 rainbox(20).
+  plot(w,polycol=ccc,close=TRUE)
+\end{ExampleCode}
+\end{Examples}
+
@@ -0,0 +1,77 @@
+\HeaderA{tile.list}{Create a list of tiles in a tessellation.}{tile.list}
+\keyword{spatial}{tile.list}
+\begin{Description}\relax
+For each point in the set being tessellated produces a list
+entry describing the Dirichlet/Voronoi tile containing that
+point.
+\end{Description}
+\begin{Usage}
+\begin{verbatim} tile.list(object) \end{verbatim}
+\end{Usage}
+\begin{Arguments}
+\begin{ldescription}
+\item[\code{object}] An object of class \code{deldir} as produced
+by the function \code{\LinkA{deldir}{deldir}()}.
+\end{ldescription}
+\end{Arguments}
+\begin{Value}
+A list with one entry for each of the points in the set
+being tesselated.  Each entry is in turn a list with
+components
+
+\begin{ldescription}
+\item[\code{pt}] The coordinates of the point whose tile is being described.
+\item[\code{x}] The \code{x} coordinates of the vertices of the tile, in
+anticlockwise order.
+\item[\code{y}] The \code{y} coordinates of the vertices of the tile, in
+anticlockwise order.
+\item[\code{bp}] Vector of logicals indicating whether the tile vertex is a
+``real'' vertex, or a \emph{boundary point}, i.e. a point where the
+tile edge intersects the boundary of the enclosing rectangle
+\end{ldescription}
+\end{Value}
+\begin{Section}{Acknowledgement}
+The author expresses sincere thanks to Majid Yazdani who found and
+pointed out a serious bug in \code{tile.list} in the previous
+version (0.0-5) of the \code{deldir} package.
+\end{Section}
+\begin{Section}{Warning}
+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 \code{\LinkA{deldir}{deldir}()}.
+
+In contrast to the previous version of \code{deldir}, the
+corners of the enclosing rectangle are now include 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 when the corresponding
+polygon is plotted, the result is the intersection of the tile
+with the enclosing rectangular window.
+\end{Section}
+\begin{Author}\relax
+Rolf Turner
+\email{r.turner@auckland.ac.nz}
+\url{http://www.math.unb.ca/~rolf}
+\end{Author}
+\begin{SeeAlso}\relax
+\code{\LinkA{deldir}{deldir}()}, \code{\LinkA{plot.tile.list}{plot.tile.list}()}
+\end{SeeAlso}
+\begin{Examples}
+\begin{ExampleCode}
+        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)
+
+\end{ExampleCode}
+\end{Examples}
+
@@ -0,0 +1,35 @@
+subroutine acchk(i,j,k,anticl,x,y,ntot,eps)
+# Check whether vertices i, j, k, are in anti-clockwise order.
+# Called by locn, qtest, qtest1.
+
+implicit double precision(a-h,o-z)
+dimension x(-3:ntot), y(-3:ntot), xt(3), yt(3)
+logical anticl
+
+# Create indicator telling which of i, j, and k are ideal points.
+if(i<=0) i1 = 1
+else i1 = 0
+if(j<=0) j1 = 1
+else j1 = 0
+if(k<=0) k1 = 1
+else k1 = 0
+ijk = i1*4+j1*2+k1
+
+# Get the coordinates of vertices i, j, and k. (Pseudo-coordinates for
+# any ideal points.)
+xt(1) = x(i)
+yt(1) = y(i)
+xt(2) = x(j)
+yt(2) = y(j)
+xt(3) = x(k)
+yt(3) = y(k)
+
+# Get the ``normalized'' cross product.
+call cross(xt,yt,ijk,cprd)
+
+# If cprd is positive then (ij-cross-ik) is directed ***upwards*** 
+# and so i, j, k, are in anti-clockwise order; else not.
+if(cprd > eps) anticl = .true.
+else anticl = .false.
+return
+end
@@ -0,0 +1,36 @@
+subroutine addpt(j,nadj,madj,x,y,ntot,eps,nerror)
+# Add point j to the triangulation.
+# Called by master, dirseg.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+logical didswp
+
+# Put the new point in, joined to the vertices of its
+# enclosing triangle.
+call initad(j,nadj,madj,x,y,ntot,eps,nerror)
+if(nerror > 0) return
+
+# Look at each `gap', i.e. pair of adjacent segments
+# emanating from the new point; they form two sides of a
+# quadrilateral; see whether the extant diagonal of this
+# quadrilateral should be swapped with its alternate
+# (according to the LOP: local optimality principle).
+now = nadj(j,1)
+nxt = nadj(j,2)
+ngap = 0
+repeat {
+        call swap(j,now,nxt,didswp,nadj,madj,x,y,ntot,eps,nerror)
+	if(nerror > 0) return
+        n = nadj(j,0)
+        if(!didswp) {         # If no swap of diagonals
+                now  = nxt    # move to the next gap.
+                ngap = ngap+1
+        }
+        call succ(nxt,j,now,nadj,madj,ntot,nerror)
+	if(nerror > 0) return
+}
+until(ngap==n)
+
+return
+end
@@ -0,0 +1,42 @@
+subroutine adjchk(i,j,adj,nadj,madj,ntot,nerror)
+# Check if vertices i and j are adjacent.
+# Called by insrt, delet, trifnd, swap, delseg, dirseg.
+
+dimension nadj(-3:ntot,0:madj)
+logical adj
+
+nerror = -1
+# Check if j is in the adjacency list of i.
+adj = .false.
+ni  = nadj(i,0)
+if(ni>0) {
+        do k = 1,ni {
+                if(j==nadj(i,k)) {
+                        adj = .true.
+                        break
+                }
+        }
+}
+
+# Check if i is in the adjacency list of j.
+nj = nadj(j,0)
+if(nj>0) {
+        do k = 1,nj {
+                if(i==nadj(j,k)) {
+                        if(adj) return # Have j in i's list and i in j's.
+                        else {
+				nerror = 1
+				return
+                        }
+                }
+        }
+}
+
+# If we get to here i is not in j's list.
+if(adj) { # If adj is true, then j IS in i's list.
+	nerror = 1
+	return
+}
+
+return
+end
@@ -0,0 +1,81 @@
+subroutine binsrt(x,y,ntot,rw,npd,ind,tx,ty,ilst,nerror)
+# Sort the data points into bins.
+# Called by master.
+
+implicit double precision(a-h,o-z)
+dimension x(-3:ntot), y(-3:ntot), tx(npd), ty(npd)
+dimension ind(npd), ilst(npd)
+dimension rw(4)
+
+nerror = -1
+kdiv   = 1+dble(npd)**0.25 # Round high.
+xkdiv  = dble(kdiv)
+
+# Dig out the corners of the rectangular window.
+xmin = rw(1)
+xmax = rw(2)
+ymin = rw(3)
+ymax = rw(4)
+
+w = xmax-xmin
+h = ymax-ymin
+
+# Number of bins is to be approx. sqrt(npd); thus number of subdivisions
+# on each side of rectangle is approx. npd**(1/4).
+dw  = w/xkdiv
+dh  = h/xkdiv
+
+# The width of each bin is dw; the height is dh.  We shall move across
+# the rectangle from left to right, then up, then back from right to
+# left, then up, ....  Note that kx counts the divisions from the left,
+# ky counts the divisions from the bottom; kx is incremented by ink, which
+# is +/- 1 and switches sign when ky is incremented; ky is always
+# incremented by 1.
+kx   = 1
+ky   = 1
+ink  = 1
+k    = 0
+do i = 1,npd { ilst(i) = 0 } # Keeps a list of those points already added
+while(ky<=kdiv) {             # to the new list.
+        do i = 1,npd {
+                if(ilst(i)==1) next  # The i-th point has already been added
+                                     # to the new list.
+                # If the i-th point is in the current bin, add it to the list.
+                xt = x(i)
+                yt = y(i)
+                ix = 1+(xt-xmin)/dw
+                if(ix>kdiv) ix = kdiv
+                jy = 1+(yt-ymin)/dh
+                if(jy>kdiv) jy = kdiv
+                if(ix==kx&jy==ky) {
+                        k = k+1
+                        ind(i)  = k  # Index i is the pos'n. of (x,y) in the
+                        tx(k)   = xt # old list; k is its pos'n. in the new one.
+                        ty(k)   = yt
+                        ilst(i) = 1  # Cross the i-th point off the old list.
+                }
+        }
+        # Move to the next bin.
+        kc = kx+ink
+        if((1<=kc)&(kc<=kdiv)) kx = kc
+        else {
+                ky  = ky+1
+                ink = -ink
+        }
+}
+
+# Check that all points from old list have been added to the new,
+# with no spurious additions.
+if(k!=npd) {
+	nerror = 2
+	return
+}
+
+# Copy the new sorted vector back on top of the old ones.
+do i = 1,npd {
+        x(i)   = tx(i)
+        y(i)   = ty(i)
+}
+
+return
+end
@@ -0,0 +1,62 @@
+subroutine circen(i,j,k,x0,y0,x,y,ntot,eps,collin,nerror)
+# Find the circumcentre (x0,y0) of the triangle with
+# vertices (x(i),y(i)), (x(j),y(j)), (x(k),y(k)).
+# Called by qtest1, dirseg, dirout.
+
+implicit double precision(a-h,o-z)
+dimension x(-3:ntot), y(-3:ntot), xt(3), yt(3)
+logical collin
+
+nerror = -1
+
+# Get the coordinates.
+xt(1) = x(i)
+yt(1) = y(i)
+xt(2) = x(j)
+yt(2) = y(j)
+xt(3) = x(k)
+yt(3) = y(k)
+
+# Check for collinearity
+ijk = 0
+call cross(xt,yt,ijk,cprd)
+if(abs(cprd) < eps) collin = .true.
+else collin = .false.
+
+# Form the vector u from i to j, and the vector v from i to k,
+# and normalize them.
+a  = x(j) - x(i)
+b  = y(j) - y(i)
+c  = x(k) - x(i)
+d  = y(k) - y(i)
+c1 = sqrt(a*a+b*b)
+c2 = sqrt(c*c+d*d)
+a  = a/c1
+b  = b/c1
+c  = c/c2
+d  = d/c2
+
+# If the points are collinear, make sure that they're in the right
+# order --- i between j and k.
+if(collin) {
+        alpha = a*c+b*d
+        # If they're not in the right order, bring things to
+        # a shuddering halt.
+        if(alpha>0) {
+		nerror = 3
+		return
+        }
+        # Collinear, but in the right order; think of this as meaning
+        # that the circumcircle in question has infinite radius.
+        return
+}
+
+# Not collinear; go ahead, make my circumcentre.  (First, form
+# the cross product of the ***unit*** vectors, instead of the
+# ``normalized'' cross product produced by ``cross''.)
+crss = a*d - b*c
+x0 = x(i) + 0.5*(c1*d - c2*b)/crss
+y0 = y(i) + 0.5*(c2*a - c1*c)/crss
+
+return
+end
@@ -0,0 +1,81 @@
+subroutine cross(x,y,ijk,cprd)
+implicit double precision(a-h,o-z)
+dimension x(3), y(3)
+# Calculates a ``normalized'' cross product of the vectors joining
+# [x(1),y(1)] to [x(2),y(2)] and [x(3),y(3)] respectively.
+# The normalization consists in dividing by the square of the
+# shortest of the three sides of the triangle.  This normalization is
+# for the purposes of testing for collinearity; if the result is less
+# than epsilon, the smallest of the sines of the angles is less than
+# epsilon.
+
+# Set constants
+zero = 0.d0
+one  = 1.d0
+two  = 2.d0
+four = 4.d0
+
+# Adjust the coordinates depending upon which points are ideal,
+# and calculate the squared length of the shortest side.
+switch(ijk) {
+        case 0: # No ideal points; no adjustment necessary.
+		smin = -one
+		do i = 1,3 {
+			ip = i+1
+			if(ip==4) ip = 1
+			a = x(ip) - x(i)
+			b = y(ip) - y(i)
+			s = a*a+b*b
+			if(smin < zero | s < smin) smin = s
+		}
+        case 1: # Only k ideal.
+		x(2) = x(2) - x(1)
+		y(2) = y(2) - y(1)
+		x(1) = zero
+		y(1) = zero
+		cn    = sqrt(x(2)**2+y(2)**2)
+		x(2) = x(2)/cn
+		y(2) = y(2)/cn
+		smin  = one
+        case 2: # Only j ideal.
+		x(3) = x(3) - x(1)
+		y(3) = y(3) - y(1)
+		x(1) = zero
+		y(1) = zero
+		cn    = sqrt(x(3)**2+y(3)**2)
+		x(3) = x(3)/cn
+		y(3) = y(3)/cn
+		smin  = one
+        case 3: # Both j and k ideal (i not).
+		x(1) = zero
+		y(1) = zero
+		smin  = 2
+        case 4: # Only i ideal.
+		x(3) = x(3) - x(2)
+		y(3) = y(3) - y(2)
+		x(2) = zero
+		y(2) = zero
+		cn    = sqrt(x(3)**2+y(3)**2)
+		x(3) = x(3)/cn
+		y(3) = y(3)/cn
+		smin  = one
+        case 5: # Both i and k ideal (j not).
+		x(2) = zero
+		y(2) = zero
+		smin  = two
+        case 6: # Both i and j ideal (k not).
+		x(3) = zero
+		y(3) = zero
+		smin  = two
+        case 7: # All three points ideal; no adjustment necessary.
+		smin  = four
+}
+
+a = x(2)-x(1)
+b = y(2)-y(1)
+c = x(3)-x(1)
+d = y(3)-y(1)
+
+cprd = (a*d - b*c)/smin
+return
+end
@@ -0,0 +1,20 @@
+subroutine delet(i,j,nadj,madj,ntot,nerror)
+# Delete i and j from each other's adjacency lists.
+# Called by initad, swap.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj)
+logical adj
+
+# First check that they're IN each other's lists.
+call adjchk(i,j,adj,nadj,madj,ntot,nerror)
+if(nerror > 0) return
+
+# Then do the actual deletion if they are.
+if(adj) {
+        call delet1(i,j,nadj,madj,ntot)
+        call delet1(j,i,nadj,madj,ntot)
+}
+
+return
+end
@@ -0,0 +1,19 @@
+subroutine delet1(i,j,nadj,madj,ntot)
+# Delete j from the adjacency list of i.
+# Called by delet.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj)
+
+n    = nadj(i,0)
+do k = 1,n {
+        if(nadj(i,k)==j) { # Find j in the list;
+                           # then move everything back one notch.
+                do kk = k,n-1 { nadj(i,kk) = nadj(i,kk+1) }
+                nadj(i,n) = 0
+                nadj(i,0) = n-1
+                return
+        }
+}
+
+end
@@ -0,0 +1,55 @@
+subroutine delout(delsum,nadj,madj,x,y,ntot,npd,ind,nerror)
+
+# Put a summary of the Delaunay triangles with a vertex at point i,
+# for i = 1, ..., npd, into the array delsum.  Do this in the original
+# order of the points, not the order into which they have been
+# bin-sorted.
+# Called by master.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+dimension delsum(npd,4), ind(npd)
+
+do i1 = 1,npd {
+        area = 0.   # Initialize area of polygon consisting of triangles
+                    # with a vertex at point i.
+
+        # Get the point number, its coordinates and the number of
+        # (real) triangles emanating from it.
+        i = ind(i1)
+        np = nadj(i,0)
+	xi = x(i)
+	yi = y(i)
+        npt = np
+        do k = 1,np {
+                kp = k+1
+                if(kp>np) kp = 1
+                if(nadj(i,k)<=0|nadj(i,kp)<=0) npt = npt-1
+        }
+
+        # For each point in the adjacency list of point i, find its
+        # successor, and the area of the triangle determined by these
+        # three points.
+        do j1 = 1,np {
+                j = nadj(i,j1)
+                if(j<=0) next
+		xj = x(j)
+		yj = y(j)
+                call succ(k,i,j,nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+                if(k<=0) next
+		xk = x(k)
+		yk = y(k)
+                call triar(xi,yi,xj,yj,xk,yk,tmp)
+                # Downweight the area by 1/3, since each
+                # triangle eventually appears 3 times over.
+                area = area+tmp/3.
+        }
+	delsum(i1,1) = xi
+	delsum(i1,2) = yi
+	delsum(i1,3) = npt
+	delsum(i1,4) = area
+}
+
+return
+end
@@ -0,0 +1,40 @@
+subroutine delseg(delsgs,ndel,nadj,madj,x,y,ntot,ind,nerror)
+
+# Output the endpoints of the line segments joining the
+# vertices of the Delaunay triangles.
+# Called by master.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+dimension delsgs(6,1), ind(1)
+logical value
+
+# For each distinct pair of points i and j, if they are adjacent
+# then put their endpoints into the output array.
+npd = ntot-4
+kseg = 0
+do i1 = 2,npd {
+	i = ind(i1)
+        do j1 = 1,i1-1 {
+		j = ind(j1)
+                call adjchk(i,j,value,nadj,madj,ntot,nerror)
+		if(nerror>0) return
+                if(value) {
+			kseg = kseg+1
+			if(kseg > ndel) {
+				nerror = 14
+				return
+			}
+			delsgs(1,kseg) = x(i)
+			delsgs(2,kseg) = y(i)
+			delsgs(3,kseg) = x(j)
+			delsgs(4,kseg) = y(j)
+			delsgs(5,kseg) = i1
+			delsgs(6,kseg) = j1
+                }
+        }
+}
+ndel = kseg
+
+return
+end
@@ -0,0 +1,110 @@
+subroutine dirout(dirsum,nadj,madj,x,y,ntot,npd,rw,ind,eps,nerror)
+
+# Output the description of the Dirichlet tile centred at point
+# i for i = 1, ..., npd.  Do this in the original order of the
+# points, not in the order into which they have been bin-sorted.
+# Called by master.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+dimension dirsum(npd,3), ind(npd), rw(4)
+logical collin, intfnd, bptab, bptcd
+
+# Note that at this point some Delaunay neighbors may be
+# `spurious'; they are the corners of a `large' rectangle in which
+# the rectangular window of interest has been suspended.  This
+# large window was brought in simply to facilitate output concerning
+# the Dirichlet tesselation.  They were added to the triangulation
+# in the routine `dirseg' which ***must*** therefore be called before
+# this routine (`dirout') is called.  (Likewise `dirseg' must be called
+# ***after*** `delseg' and `delout' have been called.)
+
+# Dig out the corners of the rectangular window.
+xmin = rw(1)
+xmax = rw(2)
+ymin = rw(3)
+ymax = rw(4)
+
+do i1 = 1,npd {
+        area = 0. # Initialize the area of the ith tile to zero.
+	nbpt = 0  # Initialize the number of boundary points of
+                  # the ith tile to zero.
+        npt = 0   # Initialize the number of tile boundaries to zero.
+
+        i = ind(i1)
+        np = nadj(i,0)
+	xi = x(i)
+	yi = y(i)
+
+        # Output the point number, its coordinates, and the number of
+        # its Delaunay neighbors == the number of boundary segments in
+        # its Dirichlet tile.
+
+        # For each Delaunay neighbor, find the circumcentres of the
+        # triangles on each side of the segment joining point i to that
+        # neighbor.
+        do j1 = 1,np {
+                j = nadj(i,j1)
+                xj = x(j)
+                yj = y(j)
+                xij = 0.5*(xi+xj)
+                yij = 0.5*(yi+yj)
+                call pred(k,i,j,nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+                call succ(l,i,j,nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+                call circen(i,k,j,a,b,x,y,ntot,eps,collin,nerror)
+		if(nerror>0) return
+                if(collin) {
+			nerror = 13
+			return
+		}
+                call circen(i,j,l,c,d,x,y,ntot,eps,collin,nerror)
+		if(nerror>0) return
+                if(collin) {
+			nerror = 13
+			return
+		}
+
+                # Increment the area of the current Dirichlet
+                # tile (intersected with the rectangular window) by applying
+                # Stokes' Theorem to the segment of tile boundary joining
+                # (a,b) to (c,d).  (Note that the direction is anti-clockwise.)
+                call stoke(a,b,c,d,rw,tmp,sn,eps,nerror)
+		if(nerror > 0) return
+                area = area+sn*tmp
+
+                # If a circumcentre is outside the rectangular window, replace
+                # it with the intersection of the rectangle boundary with the
+                # line joining the circumcentre to the midpoint of
+                # (xi,yi)->(xj,yj).  Then output the number of the current
+                # Delaunay neighbor and the two circumcentres (or the points
+                # with which they have been replaced).
+                call dldins(a,b,xij,yij,ai,bi,rw,intfnd,bptab)
+                if(intfnd) {
+			call dldins(c,d,xij,yij,ci,di,rw,intfnd,bptcd)
+                	if(!intfnd) {
+				nerror = 17
+				return
+			}
+			if(bptab & bptcd) {
+				xm = 0.5*(ai+ci)
+				ym = 0.5*(bi+di)
+				if(xmin<xm&xm<xmax&ymin<ym&ym<ymax) {
+					nbpt = nbpt+2
+					npt  = npt+1
+				}
+			}
+			else {
+				npt = npt + 1
+				if(bptab | bptcd) nbpt = nbpt+1
+			}
+		}
+		dirsum(i1,1) = npt
+		dirsum(i1,2) = nbpt
+		dirsum(i1,3) = area
+	}
+}
+
+return
+end
@@ -0,0 +1,139 @@
+subroutine dirseg(dirsgs,ndir,nadj,madj,x,y,ntot,rw,eps,ind,nerror)
+
+# Output the endpoints of the segments of boundary of Dirichlet
+# tiles.  (Do it economically; each such segment once and only once.)
+# Called by master.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+dimension dirsgs(8,ndir), rw(4), ind(1)
+logical collin, adjace, intfnd, bptab, bptcd, goferit
+
+nerror = -1
+
+# Add in some dummy corner points, outside the actual window.
+# Far enough out so that no resulting tile boundaries intersect the
+# window.
+
+# Note that these dummy corners are needed by the routine `dirout'
+# but will screw things up for `delseg' and `delout'.  Therefore
+# this routine (`dirseg') must be called ***before*** dirout, and
+# ***after*** delseg and delout.
+
+# Dig out the corners of the rectangular window.
+xmin = rw(1)
+xmax = rw(2)
+ymin = rw(3)
+ymax = rw(4)
+
+a = xmax-xmin
+b = ymax-ymin
+c = sqrt(a*a+b*b)
+
+npd  = ntot-4
+nstt = npd+1
+i = nstt
+x(i) = xmin-c
+y(i) = ymin-c
+i = i+1
+x(i) = xmax+c
+y(i) = ymin-c
+i = i+1
+x(i) = xmax+c
+y(i) = ymax+c
+i = i+1
+x(i) = xmin-c
+y(i) = ymax+c
+
+do j = nstt,ntot {
+	call addpt(j,nadj,madj,x,y,ntot,eps,nerror)
+	if(nerror > 0) return
+}
+
+# Put the segments into the array dirsgs.
+
+# For each distinct pair of (genuine) data points, find out if they are
+# adjacent.  If so, find the circumcentres of the triangles lying on each
+# side of the segment joining them.
+kseg = 0
+do i1 = 2,npd {
+	i = ind(i1)
+        do j1 = 1,i1-1 {
+		j = ind(j1)
+                call adjchk(i,j,adjace,nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+                if(adjace) {
+                        xi = x(i)
+                        yi = y(i)
+                        xj = x(j)
+                        yj = y(j)
+                        # Let (xij,yij) be the midpoint of the segment joining
+                        # (xi,yi) to (xj,yj).
+                        xij = 0.5*(xi+xj)
+                        yij = 0.5*(yi+yj)
+                        call pred(k,i,j,nadj,madj,ntot,nerror)
+			if(nerror > 0) return
+                        call circen(i,k,j,a,b,x,y,ntot,eps,collin,nerror)
+			if(nerror > 0) return
+                        if(collin) {
+				nerror = 12
+				return
+			}
+
+                        # If a circumcentre is outside the rectangular window
+                        # of interest, draw a line joining it to the mid-point
+                        # of (xi,yi)->(xj,yj).  Find the intersection of this
+                        # line with the boundary of the window; for (a,b),
+                        # call the point of intersection (ai,bi).  For (c,d),
+                        # call it (ci,di).
+                        call dldins(a,b,xij,yij,ai,bi,rw,intfnd,bptab)
+			if(!intfnd) {
+				nerror = 16
+				return
+			}
+                        call succ(l,i,j,nadj,madj,ntot,nerror)
+			if(nerror > 0) return
+                        call circen(i,j,l,c,d,x,y,ntot,eps,collin,nerror)
+			if(nerror > 0) return
+                        if(collin) {
+				nerror = 12
+				return
+			}
+                        call dldins(c,d,xij,yij,ci,di,rw,intfnd,bptcd)
+			if(!intfnd) {
+				nerror = 16
+				return
+			}
+			goferit = .false.
+			if(bptab & bptcd) {
+				xm = 0.5*(ai+ci)
+				ym = 0.5*(bi+di)
+				if(xmin<xm&xm<xmax&ymin<ym&ym<ymax) {
+					goferit = .true.
+				}
+			}
+			if((!bptab)|(!bptcd)) goferit = .true.
+			if(goferit) {
+				kseg = kseg + 1
+				if(kseg > ndir) {
+					nerror = 15
+					return
+				}
+				dirsgs(1,kseg) = ai
+				dirsgs(2,kseg) = bi
+				dirsgs(3,kseg) = ci
+				dirsgs(4,kseg) = di
+				dirsgs(5,kseg) = i1
+				dirsgs(6,kseg) = j1
+				if(bptab) dirsgs(7,kseg) = 1.d0
+				else dirsgs(7,kseg) = 0.d0
+				if(bptcd) dirsgs(8,kseg) = 1.d0
+				else dirsgs(8,kseg) = 0.d0
+			}
+                }
+        }
+}
+ndir = kseg
+
+return
+end
@@ -0,0 +1,87 @@
+subroutine dldins(a,b,c,d,ai,bi,rw,intfnd,bpt)
+
+# Get a point ***inside*** the rectangular window on the ray from
+# one circumcentre to the next one.  I.e. if the `next one' is
+# inside, then that's it; else find the intersection of this ray with
+# the boundary of the rectangle.
+# Called by dirseg, dirout.
+
+implicit double precision(a-h,o-z)
+dimension rw(4)
+logical intfnd, bpt
+
+# Note that (a,b) is the circumcenter of a Delaunay triangle,
+# and that (c,d) is the midpoint of one of its sides.
+# When `dldins' is called by `dirout' it is possible for (c,d) to
+# lie ***outside*** the rectangular window, and for the ray not to
+# intersect the window at all.  (The point (c,d) might be the midpoint
+# of a Delaunay edge connected to a `fake outer corner', added to facilitate
+# constructing a tiling that completely covers the actual window.)
+# The variable `intfnd' acts as an indicator as to whether such an
+# intersection has been found.
+
+# The variable `bpt' acts as an indicator as to whether the returned
+# point (ai,bi) is a true circumcentre, inside the window (bpt == .false.),
+# or is the intersection of a ray with the boundary of the window
+# (bpt = .true.).
+
+
+intfnd = .true.
+bpt = .true.
+
+# Dig out the corners of the rectangular window.
+xmin = rw(1)
+xmax = rw(2)
+ymin = rw(3)
+ymax = rw(4)
+
+# Check if (a,b) is inside the rectangle.
+if(xmin<=a&a<=xmax&ymin<=b&b<=ymax) {
+        ai = a
+        bi = b
+	bpt = .false.
+        return
+}
+
+# Look for appropriate intersections with the four lines forming
+# the sides of the rectangular window.
+
+# Line 1: x = xmin.
+if(a<xmin) {
+        ai = xmin
+        s = (d-b)/(c-a)
+        t = b-s*a
+        bi = s*ai+t
+        if(ymin<=bi&bi<=ymax) return
+}
+
+# Line 2: y = ymin.
+if(b<ymin) {
+        bi = ymin
+        s = (c-a)/(d-b)
+        t = a-s*b
+        ai = s*bi+t
+        if(xmin<=ai&ai<=xmax) return
+}
+
+# Line 3: x = xmax.
+if(a>xmax) {
+        ai = xmax
+        s = (d-b)/(c-a)
+        t = b-s*a
+        bi = s*ai+t
+        if(ymin<=bi&bi<=ymax) return
+}
+
+# Line 4: y = ymax.
+if(b>ymax) {
+        bi = ymax
+        s = (c-a)/(d-b)
+        t = a-s*b
+        ai = s*bi+t
+        if(xmin<=ai&ai<=xmax) return
+}
+
+intfnd = .false.
+return
+end
@@ -0,0 +1,23 @@
+subroutine inddup(x,y,n,rw,frac,dup)
+implicit double precision(a-h,o-z)
+logical dup(n)
+dimension x(n), y(n), rw(4)
+
+xtol = frac*(rw(2)-rw(1))
+ytol = frac*(rw(4)-rw(3))
+
+dup(1) = .false.
+do i = 2,n {
+	dup(i) = .false.
+	do j = 1,i-1 {
+		dx = abs(x(i)-x(j))
+		dy = abs(y(i)-y(j))
+		if(dx < xtol & dy < ytol) {
+			dup(i) = .true.
+			break
+		}
+	}
+}
+
+return
+end
@@ -0,0 +1,39 @@
+subroutine initad(j,nadj,madj,x,y,ntot,eps,nerror)
+
+# Intial adding-in of a new point j.
+# Called by addpt.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+integer tau(3)
+
+# Find the triangle containing vertex j.
+call trifnd(j,tau,nedge,nadj,madj,x,y,ntot,eps,nerror)
+if(nerror > 0) return
+
+# If the new point is on the edge of a triangle, detach the two
+# vertices of that edge from each other.  Also join j to the vertex
+# of the triangle on the reverse side of that edge from the `found'
+# triangle (defined by tau) -- given that there ***is*** such a triangle.
+if(nedge!=0) {
+        ip = nedge
+        i  = ip-1
+        if(i==0) i = 3 # Arithmetic modulo 3.
+        call pred(k,tau(i),tau(ip),nadj,madj,ntot,nerror)
+	if(nerror > 0) return
+        call succ(kk,tau(ip),tau(i),nadj,madj,ntot,nerror)
+	if(nerror > 0) return
+        call delet(tau(i),tau(ip),nadj,madj,ntot,nerror)
+	if(nerror > 0) return
+        if(k==kk) call insrt(j,k,nadj,madj,x,y,ntot,nerror,eps)
+	if(nerror > 0) return
+}
+
+# Join the new point to each of the three vertices.
+do i = 1,3 {
+        call insrt(j,tau(i),nadj,madj,x,y,ntot,nerror,eps)
+	if(nerror > 0) return
+}
+
+return
+end
@@ -0,0 +1,25 @@
+subroutine insrt(i,j,nadj,madj,x,y,ntot,nerror,eps)
+# Insert i and j into each other's adjacency list.
+# Called by master, initad, swap.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+logical adj
+
+# Check whether i and j are in each other's adjacency lists.
+call adjchk(i,j,adj,nadj,madj,ntot,nerror)
+if(nerror > 0) return
+if(adj) return
+
+# If not, find where in each list they should respectively be.
+call locn(i,j,kj,nadj,madj,x,y,ntot,eps)
+call locn(j,i,ki,nadj,madj,x,y,ntot,eps)
+
+# Put them in each other's lists in the appropriate position.
+call insrt1(i,j,kj,nadj,madj,ntot,nerror)
+if(nerror >0) return
+call insrt1(j,i,ki,nadj,madj,ntot,nerror)
+if(nerror >0) return
+
+return
+end
@@ -0,0 +1,39 @@
+subroutine insrt1(i,j,kj,nadj,madj,ntot,nerror)
+
+# Insert j into the adjacency list of i.
+# Called by insrt.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj)
+
+nerror = -1
+
+# Variable  kj is the index which j ***will***
+# have when it is inserted into the adjacency list of i in
+# the appropriate position.
+
+# If the adjacency list of i had no points just stick j into the list.
+n = nadj(i,0)
+if(n==0) {
+        nadj(i,0) = 1
+        nadj(i,1) = j
+        return
+}
+
+# If the adjacency list had some points, move everything ahead of the
+# kj-th place one place forward, and put j in position kj.
+kk = n+1
+
+if(kk>madj) { # Watch out for over-writing!!!
+	nerror = 4
+	return
+}
+while(kk>kj) {
+        nadj(i,kk) = nadj(i,kk-1)
+        kk = kk-1
+}
+nadj(i,kj) = j
+nadj(i,0)  = n+1
+
+return
+end
@@ -0,0 +1,50 @@
+subroutine locn(i,j,kj,nadj,madj,x,y,ntot,eps)
+
+# Find the appropriate location for j in the adjacency list
+# of i.  This is the index which j ***will*** have when
+# it is inserted into the adjacency list of i in the
+# appropriate place.
+# Called by insrt.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+logical before
+
+n = nadj(i,0)
+
+# If there is nothing already adjacent to i, then j will have place 1.
+if(n==0) {
+        kj = 1
+        return
+}
+
+# Run through i's list, checking if j should come before each element
+# of that list.  (I.e. if i, j, and k are in anti-clockwise order.)
+# If j comes before the kj-th item, but not before the (kj-1)-st, then
+# j should have place kj.
+do ks = 1,n {
+	kj = ks
+        k = nadj(i,kj)
+        call acchk(i,j,k,before,x,y,ntot,eps)
+        if(before) {
+                km = kj-1
+                if(km==0) km = n
+                k = nadj(i,km)
+                call acchk(i,j,k,before,x,y,ntot,eps)
+                if(before) next
+                # If j is before 1 and after n, then it should
+                # have place n+1.
+                if(kj==1) kj = n+1
+                return
+        }
+}
+
+# We've gone right through the list and haven't been before
+# the kj-th item ***and*** after the (kj-1)-st on any occasion.
+# Therefore j is before everything (==> place 1) or after
+# everything (==> place n+1).
+if(before) kj = 1
+else kj = n+1
+
+return
+end
@@ -0,0 +1,84 @@
+subroutine master(x,y,sort,rw,npd,ntot,nadj,madj,ind,tx,ty,ilst,eps,
+                  delsgs,ndel,delsum,dirsgs,ndir,dirsum,nerror)
+
+# Master subroutine, to organize all the others.
+
+implicit double precision(a-h,o-z)
+logical sort
+dimension x(-3:ntot), y(-3:ntot)
+dimension nadj(-3:ntot,0:madj)
+dimension ind(npd), tx(npd), ty(npd), ilst(npd), rw(4)
+dimension delsgs(6,ndel), dirsgs(8,ndir)
+dimension delsum(npd,4), dirsum(npd,3)
+
+# Define one.
+one = 1.d0
+
+# Sort the points into bins, the number of such being approx. sqrt(n).
+if(sort) {
+	call binsrt(x,y,ntot,rw,npd,ind,tx,ty,ilst,nerror)
+	if(nerror > 0) return
+}
+else {
+	do i = 1,npd {
+		ind(i) = i
+	}
+}
+
+# Initialize the adjacency list to 0.
+do i = -3,ntot {
+	do j = 0,madj {
+		nadj(i,j) = 0
+	}
+}
+
+# Put the four ideal points into x and y and the adjacency list.
+# The ideal points are given pseudo-coordinates
+# (-1,-1), (1,-1), (1,1), and (-1,1).  They are numbered as
+#    0       -1      -2         -3
+# i.e. the numbers decrease anticlockwise from the
+# `bottom left corner'.
+x(-3) = -one
+y(-3) =  one
+x(-2) =  one
+y(-2) =  one
+x(-1) =  one
+y(-1) = -one
+x(0)  = -one
+y(0)  = -one
+
+do i = 1,4 {
+        j = i-4
+        k = j+1
+        if(k>0) k = -3
+	call insrt(j,k,nadj,madj,x,y,ntot,nerror,eps)
+        if(nerror>0) return
+}
+
+# Put in the first of the point set into the adjacency list.
+do i = 1,4 {
+        j = i-4
+	call insrt(1,j,nadj,madj,x,y,ntot,nerror,eps)
+        if(nerror>0) return
+}
+
+# Now add the rest of the point set
+do j = 2,npd {
+	call addpt(j,nadj,madj,x,y,ntot,eps,nerror)
+        if(nerror>0) return
+}
+
+# Obtain the description of the triangulation.
+call delseg(delsgs,ndel,nadj,madj,x,y,ntot,ind,nerror)
+if(nerror>0) return
+
+call delout(delsum,nadj,madj,x,y,ntot,npd,ind,nerror)
+if(nerror>0) return
+
+call dirseg(dirsgs,ndir,nadj,madj,x,y,ntot,rw,eps,ind,nerror)
+if(nerror>0) return
+
+call dirout(dirsum,nadj,madj,x,y,ntot,npd,rw,ind,eps,nerror)
+
+return
+end
@@ -0,0 +1,20 @@
+subroutine mnnd(x,y,n,dminbig,dminav)
+implicit double precision(a-h,o-z)
+dimension x(n), y(n)
+
+dminav = 0.d0
+do i = 1,n {
+	dmin = dminbig
+	do j = 1,n {
+		if(i!=j) {
+			d = (x(i)-x(j))**2 + (y(i)-y(j))**2
+			if(d < dmin) dmin = d
+		}
+	}
+	dminav = dminav + sqrt(dmin)
+}
+
+dminav = dminav/n
+
+return
+end
@@ -0,0 +1,34 @@
+subroutine pred(kpr,i,j,nadj,madj,ntot,nerror)
+
+# Find the predecessor of j in the adjacency list of i.
+# Called by initad, trifnd, swap, dirseg, dirout.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj)
+
+nerror = -1
+
+n = nadj(i,0)
+
+# If the adjacency list of i is empty, then clearly j has no predecessor
+# in this adjacency list. Something's wrong; stop.
+if(n==0) {
+	nerror = 5
+	return
+}
+
+# The adjacency list of i is non-empty; search through it until j is found;
+# subtract 1 from the location of j, and find the contents of this new location
+do k = 1,n {
+        if(j==nadj(i,k)) {
+                km = k-1
+                if(km<1) km = n         # Take km modulo n. (The adjacency list
+                kpr = nadj(i,km)        # is circular.)
+                return
+        }
+}
+
+# The adjacency list doesn't contain j.  Something's wrong; stop.
+nerror = 6
+return
+end
@@ -0,0 +1,106 @@
+subroutine qtest(h,i,j,k,shdswp,x,y,ntot,eps,nerror)
+
+# Test whether the LOP is satisified; i.e. whether vertex j is
+# outside the circumcircle of vertices h, i, and k of the
+# quadrilateral.  Vertex h is the vertex being added; i and k are the
+# vertices of the quadrilateral which are currently joined; j is the
+# vertex which is ``opposite'' the vertex being added.  If the LOP is
+# not satisfied, the shdswp ("should-swap") is true, i.e. h and j
+# should be joined, rather than i and k.
+# Called by swap.
+
+implicit double precision(a-h,o-z)
+dimension x(-3:ntot), y(-3:ntot)
+integer h
+logical shdswp
+
+nerror = -1
+
+# Look for ideal points.
+if(i<=0) ii = 1
+else ii = 0
+if(j<=0) jj = 1
+else jj = 0
+if(k<=0) kk = 1
+else kk = 0
+ijk = ii*4+jj*2+kk
+switch(ijk) {
+
+# All three corners other than h (the point currently being
+# added) are ideal --- so h, i, and k are co-linear; so 
+# i and k shouldn't be joined, and h should be joined to j.
+# So swap.  (But this can't happen, anyway!!!)
+        case 7:
+		shdswp = .true.
+                return
+
+# If i and j are ideal, find out which of h and k is closer to the
+# intersection point of the two diagonals, and choose the diagonal
+# emanating from that vertex.  (I.e. if h is closer, swap.)
+# Unless swapping yields a non-convex quadrilateral!!!
+        case 6:
+		ss = 1 - 2*mod(-j,2)
+                xh = x(h)
+                yh = y(h)
+                xk = x(k)
+                yk = y(k)
+                test =(xh*yk+xk*yh-xh*yh-xk*yk)*ss
+                if(test>0) shdswp = .true.
+                else shdswp = .false.
+# Check for convexity:
+		if(shdswp) call acchk(j,k,h,shdswp,x,y,ntot,eps)
+                return
+
+# Vertices i and k are ideal --- can't happen, but if it did, we'd
+# increase the minimum angle ``from 0 to more than 2*0'' by swapping,
+# and the quadrilateral would definitely be convex, so let's say swap.
+        case 5:
+		shdswp = .true.
+                return
+
+# If i is ideal we'd increase the minimum angle ``from 0 to more than
+# 2*0'' by swapping, so just check for convexity:
+        case 4:
+                call acchk(j,k,h,shdswp,x,y,ntot,eps)
+                return
+
+# If j and k are ideal, this is like unto case 6.
+        case 3:
+		ss = 1 - 2*mod(-j,2)
+                xi = x(i)
+                yi = y(i)
+                xh = x(h)
+                yh = y(h)
+                test = (xh*yi+xi*yh-xh*yh-xi*yi)*ss
+                if(test>0.) shdswp = .true.
+                else shdswp = .false.
+# Check for convexity:
+		if(shdswp) call acchk(h,i,j,shdswp,x,y,ntot,eps)
+                return
+
+# If j is ideal we'd decrease the minimum angle ``from more than 2*0
+# to 0'', by swapping; so don't swap.
+        case 2:
+                shdswp = .false.
+                return
+
+# If k is ideal, this is like unto case 4.
+        case 1:
+                call acchk(h,i,j,shdswp,x,y,ntot,eps) # This checks
+                                                      # for convexity.
+                                                      # (Was i,j,h,...)
+                return
+
+# If none of the `other' three corners are ideal, do the Lee-Schacter
+# test for the LOP.
+        case 0:
+                call qtest1(h,i,j,k,x,y,ntot,eps,shdswp,nerror)
+                return
+
+        default:  # This CAN'T happen.
+		nerror = 7
+		return
+
+}
+
+end
@@ -0,0 +1,49 @@
+subroutine qtest1(h,i,j,k,x,y,ntot,eps,shdswp,nerror)
+
+# The Lee-Schacter test for the LOP (all points are
+# real, i.e. non-ideal).  If the LOP is ***not***
+# satisfied then the diagonals should be swapped,
+# i.e. shdswp ("should-swap") is true.
+# Called by qtest.
+
+implicit double precision(a-h,o-z)
+dimension x(-3:ntot), y(-3:ntot)
+integer h
+logical shdswp
+
+# The vertices of the quadrilateral are labelled
+# h, i, j, k in the anticlockwise direction, h
+# being the point of central interest.
+
+# Make sure the quadrilateral is convex, so that
+# it makes sense to swap the diagonal.
+call acchk(i,j,k,shdswp,x,y,ntot,eps)
+if(!shdswp) return
+
+# Get the coordinates of vertices h and j.
+xh = x(h)
+yh = y(h)
+xj = x(j)
+yj = y(j)
+
+# Find the centre of the circumcircle of vertices h, i, k.
+call circen(h,i,k,x0,y0,x,y,ntot,eps,shdswp,nerror)
+if(nerror>0) return
+if(shdswp) return # The points h, i, and k are colinear, so
+                  # the circumcircle has `infinite radius', so
+                  # (xj,yj) is definitely inside.
+
+# Check whether (xj,yj) is inside the circle of centre
+# (x0,y0) and radius r = dist[(x0,y0),(xh,yh)]
+
+a = x0-xh
+b = y0-yh
+r2 = a*a+b*b
+a = x0-xj
+b = y0-yj
+ch = a*a+b*b
+if(ch<r2) shdswp = .true.
+else shdswp = .false.
+
+return
+end
@@ -0,0 +1,165 @@
+subroutine stoke(x1,y1,x2,y2,rw,area,s1,eps,nerror)
+
+# Apply Stokes' theorem to find the area of a polygon;
+# we are looking at the boundary segment from (x1,y1)
+# to (x2,y2), travelling anti-clockwise.  We find the
+# area between this segment and the horizontal base-line
+# y = ymin, and attach a sign s1.  (Positive if the
+# segment is right-to-left, negative if left to right.)
+# The area of the polygon is found by summing the result
+# over all boundary segments.
+
+# Just in case you thought this wasn't complicated enough,
+# what we really want is the area of the intersection of
+# the polygon with the rectangular window that we're using.
+
+# Called by dirout.
+
+implicit double precision(a-h,o-z)
+dimension rw(4)
+logical value
+
+zero   = 0.d0
+nerror = -1
+
+# If the segment is vertical, the area is zero.
+call testeq(x1,x2,eps,value)
+if(value) {
+        area = 0.
+        s1   = 0.
+        return
+}
+
+# Find which is the right-hand end, and which is the left.
+if(x1<x2) {
+        xl = x1
+        yl = y1
+        xr = x2
+        yr = y2
+        s1 = -1.
+}
+else {
+        xl = x2
+        yl = y2
+        xr = x1
+        yr = y1
+        s1 = 1.
+}
+
+# Dig out the corners of the rectangular window.
+xmin = rw(1)
+xmax = rw(2)
+ymin = rw(3)
+ymax = rw(4)
+
+# Now start intersecting with the rectangular window.
+# Truncate the segment in the horizontal direction at
+# the edges of the rectangle.
+slope = (yl-yr)/(xl-xr)
+x  = max(xl,xmin)
+y  = yl+slope*(x-xl)
+xl = x
+yl = y
+
+x  = min(xr,xmax)
+y  = yr+slope*(x-xr)
+xr = x
+yr = y
+
+if(xr<=xmin|xl>=xmax) {
+        area = 0.
+        return
+}
+
+# We're now looking at a trapezoidal region which may or may
+# not protrude above or below the horizontal strip bounded by
+# y = ymax and y = ymin.
+ybot = min(yl,yr)
+ytop = max(yl,yr)
+
+# Case 1; ymax <= ybot:
+# The `roof' of the trapezoid is entirely above the
+# horizontal strip.
+if(ymax<=ybot) {
+        area = (xr-xl)*(ymax-ymin)
+        return
+}
+
+# Case 2; ymin <= ybot <= ymax <= ytop:
+# The `roof' of the trapezoid intersects the top of the
+# horizontal strip (y = ymax) but not the bottom (y = ymin).
+if(ymin<=ybot&ymax<=ytop) {
+        call testeq(slope,zero,eps,value)
+        if(value) {
+                w1 = 0.
+                w2 = xr-xl
+        }
+        else {
+                xit = xl+(ymax-yl)/slope
+                w1 = xit-xl
+                w2 = xr-xit
+                if(slope<0.) {
+                        tmp = w1
+                        w1  = w2
+                        w2  = tmp
+                }
+        }
+        area = 0.5*w1*((ybot-ymin)+(ymax-ymin))+w2*(ymax-ymin)
+        return
+}
+
+# Case 3; ybot <= ymin <= ymax <= ytop:
+# The `roof' intersects both the top (y = ymax) and
+# the bottom (y = ymin) of the horizontal strip.
+if(ybot<=ymin&ymax<=ytop) {
+        xit = xl+(ymax-yl)/slope
+        xib = xl+(ymin-yl)/slope
+        if(slope>0.) {
+                w1 = xit-xib
+                w2 = xr-xit
+        }
+        else {
+                w1 = xib-xit
+                w2 = xit-xl
+        }
+        area = 0.5*w1*(ymax-ymin)+w2*(ymax-ymin)
+        return
+}
+
+# Case 4; ymin <= ybot <= ytop <= ymax:
+# The `roof' is ***between*** the bottom (y = ymin) and
+# the top (y = ymax) of the horizontal strip.
+if(ymin<=ybot&ytop<=ymax) {
+        area = 0.5*(xr-xl)*((ytop-ymin)+(ybot-ymin))
+        return
+}
+
+# Case 5; ybot <= ymin <= ytop <= ymax:
+# The `roof' intersects the bottom (y = ymin) but not
+# the top (y = ymax) of the horizontal strip.
+if(ybot<=ymin&ymin<=ytop) {
+        call testeq(slope,zero,eps,value)
+        if(value) {
+                area = 0.
+                return
+        }
+        xib = xl+(ymin-yl)/slope
+        if(slope>0.) w = xr-xib
+        else w = xib-xl
+        area = 0.5*w*(ytop-ymin)
+        return
+}
+
+# Case 6; ytop <= ymin:
+# The `roof' is entirely below the bottom (y = ymin), so
+# there is no area contribution at all.
+if(ytop<=ymin) {
+        area = 0.
+        return
+}
+
+# Default; all stuffed up:
+nerror = 8
+return
+
+end
@@ -0,0 +1,34 @@
+subroutine  succ(ksc,i,j,nadj,madj,ntot,nerror)
+
+# Find the successor of j in the adjacency list of i.
+# Called by addpt, initad, trifnd, swap, delout, dirseg, dirout.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj)
+
+nerror = -1
+
+n = nadj(i,0)
+
+# If the adjacency list of i is empty, then clearly j has no successor
+# in this adjacency list.  Something's wrong; stop.
+if(n==0) {
+	nerror = 9
+	return
+}
+
+# The adjacency list of i is non-empty; search through it until j is found;
+# add 1 to the location of j, and find the contents of this new location.
+do k = 1,n {
+        if(j==nadj(i,k)) {
+                kp = k+1
+                if(kp>n) kp = 1         # Take kp modulo n. (The adjacency list
+                ksc = nadj(i,kp)        # is circular.)
+                return
+        }
+}
+
+# The adjacency list doesn't contain j.  Something's wrong.
+nerror = 10
+return
+end
@@ -0,0 +1,43 @@
+subroutine swap(j,k1,k2,shdswp,nadj,madj,x,y,ntot,eps,nerror)
+
+# The segment k1->k2 is a diagonal of a quadrilateral
+# with a vertex at j (the point being added to the
+# triangulation).  If the LOP is not satisfied, swap
+# it for the other diagonal.
+# Called by addpt.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
+logical shdswp
+
+
+# If vertices k1 and k2 are not connected there is no diagonal to swap.
+# This could happen if vertices j, k1, and k2 were colinear, but shouldn't.
+call adjchk(k1,k2,shdswp,nadj,madj,ntot,nerror)
+if(nerror > 0) return
+if(!shdswp) return
+
+# Get the other vertex of the quadrilateral.
+call pred(k,k1,k2,nadj,madj,ntot,nerror) # If these aren't the same, then
+if(nerror > 0) return
+call succ(kk,k2,k1,nadj,madj,ntot,nerror) # there is no other vertex.
+if(nerror > 0) return
+if(kk!=k) {
+        shdswp = .false.
+        return
+}
+
+# Check whether the LOP is satisified; i.e. whether
+# vertex k is outside the circumcircle of vertices j, k1, and k2
+call qtest(j,k1,k,k2,shdswp,x,y,ntot,eps,nerror)
+if(nerror > 0) return
+
+# Do the actual swapping.
+if(shdswp) {
+        call delet(k1,k2,nadj,madj,ntot,nerror)
+	if(nerror > 0) return
+        call insrt(j,k,nadj,madj,x,y,ntot,nerror,eps)
+	if(nerror > 0) return
+}
+return
+end
@@ -0,0 +1,32 @@
+subroutine testeq(a,b,eps,value)
+
+# Test for the equality of a and b in a fairly
+# robust way.
+# Called by trifnd, circen, stoke.
+
+implicit double precision(a-h,o-z)
+logical value
+
+# If b is essentially 0, check whether a is essentially zero also.
+# The following is very sloppy!  Must fix it!
+if(abs(b)<=eps) {
+        if(abs(a)<=eps) value = .true.
+        else value = .false.
+        return
+}
+
+# Test if a is a `lot different' from b.  (If it is
+# they're obviously not equal.)  This avoids under/overflow
+# problems in dividing a by b.
+if(abs(a)>10.*abs(b)|abs(a)<0.1*abs(b)) {
+        value = .false.
+        return
+}
+
+# They're non-zero and fairly close; compare their ratio with 1.
+c = a/b
+if(abs(c-1.)<=eps) value = .true.
+else value = .false.
+
+return
+end
@@ -0,0 +1,11 @@
+subroutine triar(x0,y0,x1,y1,x2,y2,area)
+
+# Calculate the area of a triangle with given
+# vertices, in anti clockwise direction.
+# Called by delout.
+
+implicit double precision(a-h,o-z)
+
+area = 0.5*((x1-x0)*(y2-y0)-(x2-x0)*(y1-y0))
+return
+end
@@ -0,0 +1,112 @@
+subroutine trifnd(j,tau,nedge,nadj,madj,x,y,ntot,eps,nerror)
+
+# Find the triangle of the extant triangulation in which
+# lies the point currently being added.
+# Called by initad.
+
+implicit double precision(a-h,o-z)
+dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot), xt(3), yt(3)
+integer tau(3)
+logical adjace
+
+nerror = -1
+
+# The first point must be added to the triangulation before
+# calling trifnd.
+if(j==1) {
+	nerror = 11
+	return
+}
+
+# Get the previous triangle:
+j1     = j-1
+tau(1) = j1
+tau(3) = nadj(j1,1)
+call pred(tau(2),j1,tau(3),nadj,madj,ntot,nerror)
+if(nerror > 0) return
+call adjchk(tau(2),tau(3),adjace,nadj,madj,ntot,nerror)
+if(nerror>0) return
+if(!adjace) {
+        tau(3) = tau(2)
+	call pred(tau(2),j1,tau(3),nadj,madj,ntot,nerror)
+	if(nerror > 0) return
+}
+
+# Move to the adjacent triangle in the direction of the new
+# point, until the new point lies in this triangle.
+1       continue
+ntau  = 0 # This number will identify the triangle to be moved to.
+nedge = 0 # If the point lies on an edge, this number will identify that edge.
+do i = 1,3 {
+        ip = i+1
+        if(ip==4) ip = 1 # Take addition modulo 3.
+
+# Get the coordinates of the vertices of the current side,
+# and of the point j which is being added:
+	xt(1) = x(tau(i))
+	yt(1) = y(tau(i))
+	xt(2) = x(tau(ip))
+	yt(2) = y(tau(ip))
+	xt(3) = x(j)
+	yt(3) = y(j)
+
+# Create indicator telling which of tau(i), tau(ip), and j
+# are ideal points.  (The point being added, j, is ***never*** ideal.)
+	if(tau(i)<=0) i1 = 1
+	else i1 = 0
+	if(tau(ip)<=0) j1 = 1
+	else j1 = 0
+	k1 = 0
+	ijk = i1*4+j1*2+k1
+
+# Calculate the ``normalized'' cross product; if this is positive
+# then the point being added is to the left (as we move along the
+# edge in an anti-clockwise direction).  If the test value is positive
+# for all three edges, then the point is inside the triangle.  Note
+# that if the test value is very close to zero, we might get negative
+# values for it on both sides of an edge, and hence go into an
+# infinite loop.
+	call cross(xt,yt,ijk,cprd)
+	if(cprd >= eps) continue
+	else if(cprd > -eps) nedge = ip
+	else {
+		ntau = ip
+		break
+	}
+}
+
+# We've played ring-around-the-triangle; now figure out the
+# next move:
+switch(ntau) {
+        case 0: # All tests >= 0.; the point is inside; return.
+                return
+
+# The point is not inside; work out the vertices of the triangle to which
+# to move.  Notation: Number the vertices of the current triangle from 1 to 3,
+# anti-clockwise. Then "triangle i+1" is adjacent to the side from vertex i to
+# vertex i+1, where i+1 is taken modulo 3 (i.e. "3+1 = 1").
+        case 1:
+                # move to "triangle 1"
+                #tau(1) = tau(1)
+                tau(2)  = tau(3)
+		call succ(tau(3),tau(1),tau(2),nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+        case 2:
+                # move to "triangle 2"
+                #tau(1) = tau(1)
+                tau(3)  = tau(2)
+		call pred(tau(2),tau(1),tau(3),nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+        case 3:
+                # move to "triangle 3"
+                tau(1)  = tau(3)
+                #tau(2) = tau(2)
+		call succ(tau(3),tau(1),tau(2),nadj,madj,ntot,nerror)
+		if(nerror > 0) return
+}
+
+# We've moved to a new triangle; check if the point being added lies
+# inside this one.
+go to 1
+
+end
...	...	@@ -0,0 +1,29 @@
	1	+Entry: deldir-internal
	2	+Aliases: dumpts ind.dup mid.in mnnd get.cnrind acw
	3	+Keywords: internal
	4	+Description: Internal deldir functions
	5	+URL: ../../../library/deldir/html/deldir-internal.html
	6	+
	7	+Entry: deldir
	8	+Aliases: deldir
	9	+Keywords: spatial
	10	+Description: Construct the Delaunay triangulation and the Dirichlet (Voronoi) tessellation of a planar point set.
	11	+URL: ../../../library/deldir/html/deldir.html
	12	+
	13	+Entry: plot.deldir
	14	+Aliases: plot.deldir
	15	+Keywords:
	16	+Description: Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi) tesselation of a planar point set, as constructed by the function deldir.
	17	+URL: ../../../library/deldir/html/plot.deldir.html
	18	+
	19	+Entry: plot.tile.list
	20	+Aliases: plot.tile.list
	21	+Keywords: hplot
	22	+Description: Plot Dirchlet/Voronoi tiles
	23	+URL: ../../../library/deldir/html/plot.tile.list.html
	24	+
	25	+Entry: tile.list
	26	+Aliases: tile.list
	27	+Keywords: spatial
	28	+Description: Create a list of tiles in a tessellation.
	29	+URL: ../../../library/deldir/html/tile.list.html
...	...
...	...	@@ -0,0 +1,14 @@
	1	+Package: deldir
	2	+Version: 0.0-7
	3	+Date: 2007-11-03
	4	+Title: Delaunay Triangulation and Dirichlet (Voronoi) Tessellation.
	5	+Author: Rolf Turner <r.turner@auckland.ac.nz>
	6	+Maintainer: Rolf Turner <r.turner@auckland.ac.nz>
	7	+Depends: R (>= 0.99)
	8	+Description: Calculates the Delaunay triangulation and the Dirichlet or
	9	+ Voronoi tessellation (with respect to the entire plane) of a
	10	+ planar point set.
	11	+License: GPL version 2 or newer
	12	+URL: http://www.math.unb.ca/~rolf/
	13	+Packaged: Sat Nov 3 14:49:01 2007; rolf
	14	+Built: R 2.6.0; i386-pc-mingw32; 2007-11-04 12:29:56; windows
...	...
...	...	@@ -0,0 +1,9 @@
	1	+deldir Construct the Delaunay triangulation and the
	2	+ Dirichlet (Voronoi) tessellation of a planar
	3	+ point set.
	4	+plot.deldir Produce a plot of the Delaunay triangulation
	5	+ and Dirichlet (Voronoi) tesselation of a planar
	6	+ point set, as constructed by the function
	7	+ deldir.
	8	+plot.tile.list Plot Dirchlet/Voronoi tiles
	9	+tile.list Create a list of tiles in a tessellation.
...	...
...	...	@@ -0,0 +1,63 @@
	1	+34ef60d93db4e4513ae6b13668340663 *CONTENTS
	2	+ee8fe6e62e3e9bbcb23bc33ec089d31a *DESCRIPTION
	3	+d8d74f53836393d466d4c787015f0bec *INDEX
	4	+e5a045166078a4a7bd4663db86d87395 *Meta/Rd.rds
	5	+cd368d55dede62658f2248198fa90551 *Meta/hsearch.rds
	6	+366a9e3d80f846401b4bbea6cadc152a *Meta/package.rds
	7	+f1e590f87144c3fc1bac5ad7ba9b43c6 *R-ex/deldir.R
	8	+593a215bba63e9347baa5d20f399d843 *R-ex/plot.deldir.R
	9	+21797d2ef6002f35b9539333ad936fe9 *R-ex/plot.tile.list.R
	10	+93422ef22c461433627a006c688284c9 *R-ex/tile.list.R
	11	+3b9c9b8cbab6a88b781daa62777d5f5c *R/deldir
	12	+4e82e5306993afb8f50752201a1801bf *READ_ME
	13	+af071a8a57af7c84d1333c8c402da6a0 *chtml/deldir.chm
	14	+6b9d8dec2bc172cb280aa76cf285b67b *err.list
	15	+b1e64b5aed9484584095ac4e02baea84 *ex.out
	16	+16bd2bedcf7796770e90c3e56d22526e *help/AnIndex
	17	+4d8ead26ee3fc639d613b330b6b4478e *help/deldir
	18	+a84f7d114620c92c0ac0e2ea57265afa *help/deldir-internal
	19	+7881f61ae9a79255ce43f0814a015899 *help/plot.deldir
	20	+4a53379cb1c77bc2cf1c1c425cfbca8c *help/plot.tile.list
	21	+ec37f6d9ceddf4db09dd469dedff211f *help/tile.list
	22	+be17ba586ed30dd13fa15e313a6fd7bf *html/00Index.html
	23	+2ab5a9ca151205217138b518e6d8405d *html/deldir-internal.html
	24	+40f8a12dc94a21820c1ad4ab98886536 *html/deldir.html
	25	+49398bbaf0db6df39fcdcb6f95cf593e *html/plot.deldir.html
	26	+0cc37a506231a768a0bdf3aa3d9c14f5 *html/plot.tile.list.html
	27	+3f9352ec46c00d0fc3f407d3103173a0 *html/tile.list.html
	28	+1f390921c22e235996759df8afeb639a *latex/deldir-internal.tex
	29	+6ad0d68321c5d3be37af79dce2904e68 *latex/deldir.tex
	30	+375f59d2ef90cdfdbc38bb3aec3f45cb *latex/plot.deldir.tex
	31	+e20c7394462a8ef65f11f3f84a2e8819 *latex/plot.tile.list.tex
	32	+6d3ef222fc3e6eff800da2b0adbad24a *latex/tile.list.tex
	33	+797b9359a38e5806fe599e63d173f604 *libs/deldir.dll
	34	+f3c6c1feb137fadf8399c5541034d0c0 *man/deldir.Rd.gz
	35	+32332f6b7cc4f78b3307c77566bb5a89 *ratfor/acchk.r
	36	+5f8b878f079445829c1d9f49e31170bc *ratfor/addpt.r
	37	+93ed44f787438449a98e24be66775d8b *ratfor/adjchk.r
	38	+d53c2495df9c4d1f4d3fcad81e10b801 *ratfor/binsrt.r
	39	+f6e55e5c74ed60b802cb2a827e597603 *ratfor/circen.r
	40	+f54912f73cbfc10232547f97258e6ec6 *ratfor/cross.r
	41	+43bbe616653a6276af389d5af6ff212c *ratfor/delet.r
	42	+1791398b1b1e8db90ef7aca3d901246d *ratfor/delet1.r
	43	+22b1b7d68d94d8d179663d3d7dbee446 *ratfor/delout.r
	44	+6076c5cec0047fe59f8e151eef2041cf *ratfor/delseg.r
	45	+4e775ebb151b3b47656718f9f6aefb8c *ratfor/dirout.r
	46	+ab6664e97b26cd37a8a05ea814eada1d *ratfor/dirseg.r
	47	+77454183e29dd932e494e9f363e7c7a0 *ratfor/dldins.r
	48	+7a0cbe7ba69bc4d3758d7b467378a6a5 *ratfor/inddup.r
	49	+a88b1a39361d1816685cba69561b1733 *ratfor/initad.r
	50	+793f31965fd00d91904dc1145e5a5e95 *ratfor/insrt.r
	51	+8d92aa51712ce4e98089b46daa29614f *ratfor/insrt1.r
	52	+d393dd4b28b863d09b8940bee1e42e97 *ratfor/locn.r
	53	+e6e2b439fea000bc0c4332c1332430d1 *ratfor/master.r
	54	+adaef6c99448a066b58f1a07eb5d883c *ratfor/mnnd.r
	55	+5e7b0ef2ba903423e2b218ca93a9a66a *ratfor/pred.r
	56	+dd315504b6e25408c551b4f2aaf2927f *ratfor/qtest.r
	57	+87629c704f0b27364c80453f5034f86b *ratfor/qtest1.r
	58	+a8636bf7b643eca192877bc38a6f2c87 *ratfor/stoke.r
	59	+6d1ed8f4862aea3c6cd9a28179c51ddf *ratfor/succ.r
	60	+d2d7a614a5ce8ff20006261ac1db36a5 *ratfor/swap.r
	61	+182455cc37ee801235629c44f0277830 *ratfor/testeq.r
	62	+5264af03d65005c360f831afc0b5d160 *ratfor/triar.r
	63	+9afccbeef89e96ce284f18504820e8a0 *ratfor/trifnd.r
...	...
...	...	@@ -0,0 +1,20 @@
	1	+### Name: deldir
	2	+### Title: Construct the Delaunay triangulation and the Dirichlet (Voronoi)
	3	+### tessellation of a planar point set.
	4	+### Aliases: deldir
	5	+### Keywords: spatial
	6	+
	7	+### ** Examples
	8	+
	9	+x <- c(2.3,3.0,7.0,1.0,3.0,8.0)
	10	+y <- c(2.3,3.0,2.0,5.0,8.0,9.0)
	11	+try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
	12	+# Puts dummy points at the corners of the rectangular
	13	+# window, i.e. at (0,0), (10,0), (10,10), and (0,10)
	14	+## Not run:
	15	+##D try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr')
	16	+## End(Not run)
	17	+# Plots the triangulation which was created (but not the tesselation).
	18	+
	19	+
	20	+
...	...
...	...	@@ -0,0 +1,24 @@
	1	+### Name: plot.deldir
	2	+### Title: Produce a plot of the Delaunay triangulation and Dirichlet
	3	+### (Voronoi) tesselation of a planar point set, as constructed by the
	4	+### function deldir.
	5	+### Aliases: plot.deldir
	6	+
	7	+
	8	+### ** Examples
	9	+
	10	+## Not run:
	11	+##D try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10))
	12	+##D plot(try)
	13	+##D #
	14	+##D deldir(x,y,list(ndx=4,ndy=4),plot=TRUE,add=TRUE,wl='te',
	15	+##D col=c(1,1,2,3,4),num=TRUE)
	16	+##D # Plots the tesselation, but does not save the results.
	17	+##D try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=TRUE,wl='tr',
	18	+##D wp='n')
	19	+##D # Plots the triangulation, but not the points, and saves the returned
	20	+##D structure.
	21	+## End(Not run)
	22	+
	23	+
	24	+
...	...
...	...	@@ -0,0 +1,18 @@
	1	+### Name: plot.tile.list
	2	+### Title: Plot Dirchlet/Voronoi tiles
	3	+### Aliases: plot.tile.list
	4	+### Keywords: hplot
	5	+
	6	+### ** Examples
	7	+
	8	+ x <- runif(20)
	9	+ y <- runif(20)
	10	+ z <- deldir(x,y,rw=c(0,1,0,1))
	11	+ w <- tile.list(z)
	12	+ plot(w)
	13	+ ccc <- heat.colors(20) # Or topo.colors(20), or terrain.colors(20)
	14	+ # or cm.colors(20), or rainbox(20).
	15	+ plot(w,polycol=ccc,close=TRUE)
	16	+
	17	+
	18	+
...	...
...	...	@@ -0,0 +1,21 @@
	1	+### Name: tile.list
	2	+### Title: Create a list of tiles in a tessellation.
	3	+### Aliases: tile.list
	4	+### Keywords: spatial
	5	+
	6	+### ** Examples
	7	+
	8	+ x <- runif(20)
	9	+ y <- runif(20)
	10	+ z <- deldir(x,y)
	11	+ w <- tile.list(z)
	12	+
	13	+ z <- deldir(x,y,rw=c(0,1,0,1))
	14	+ w <- tile.list(z)
	15	+
	16	+ z <- deldir(x,y,rw=c(0,1,0,1),dpl=list(ndx=2,ndy=2))
	17	+ w <- tile.list(z)
	18	+
	19	+
	20	+
	21	+
...	...
...	...	@@ -0,0 +1,520 @@
	1	+.packageName <- "deldir"
	2	+.First.lib <- function(lib,pkg) {
	3	+ library.dynam("deldir", pkg, lib)
	4	+ ver <- read.dcf(file.path(lib, pkg, "DESCRIPTION"), "Version")
	5	+ cat(paste(pkg, ver, "\n"))
	6	+}
	7	+acw <- function(xxx) {
	8	+xbar <- mean(xxx$x)
	9	+ybar <- mean(xxx$y)
	10	+theta <- atan2(xxx$y - ybar,xxx$x-xbar)
	11	+theta <- ifelse(theta > 0, theta, theta + 2 * pi)
	12	+theta.0 <- sort(unique(theta))
	13	+iii <- match(theta.0, theta)
	14	+xxx$x <- xxx$x[iii]
	15	+xxx$y <- xxx$y[iii]
	16	+xxx$bp <- xxx$bp[iii]
	17	+xxx
	18	+}
	19	+deldir <- function(x,y,dpl=NULL,rw=NULL,eps=1e-9,frac=1e-4,
	20	+ sort=TRUE,plotit=FALSE,digits=6,...) {
	21	+# Function deldir
	22	+#
	23	+# Copyright (C) 1996 by T. Rolf Turner
	24	+#
	25	+# Permission to use, copy, modify, and distribute this software and
	26	+# its documentation for any purpose and without fee is hereby
	27	+# granted, provided that the above copyright notice appear in all
	28	+# copies and that both that copyright notice and this permission
	29	+# notice appear in supporting documentation.
	30	+#
	31	+# ORIGINALLY PROGRAMMED BY: Rolf Turner in 1987/88, while with the
	32	+# Division of Mathematics and Statistics, CSIRO, Sydney, Australia.
	33	+# Re-programmed by Rolf Turner to adapt the implementation from a
	34	+# stand-alone Fortran program to an S function, while visiting the
	35	+# University of Western Australia, May 1995. Further revised
	36	+# December 1996.
	37	+#
	38	+# Function to compute the Delaunay Triangulation (and hence the
	39	+# Dirichlet Tesselation) of a planar point set according to the
	40	+# second (iterative) algorithm of Lee and Schacter, International
	41	+# Journal of Computer and Information Sciences, Vol. 9, No. 3, 1980,
	42	+# pages 219 to 242.
	43	+
	44	+# The triangulation is made to be with respect to the whole plane by
	45	+# `suspending' it from `ideal' points
	46	+# (-R,-R), (R,-R) (R,R), and (-R,R), where R --> infinity.
	47	+
	48	+# It is also enclosed in a finite rectangle (whose boundaries truncate any
	49	+# infinite Dirichlet tiles) with corners (xmin,ymin) etc. This rectangle
	50	+# is referred to elsewhere as `the' rectangular window.
	51	+
	52	+# If the first argument is a list, extract components x and y:
	53	+if(is.list(x)) {
	54	+ if(all(!is.na(match(c('x','y'),names(x))))) {
	55	+ y <- x$y
	56	+ x <- x$x
	57	+ }
	58	+ else {
	59	+ cat('Error: called with list lacking both x and y elements\n')
	60	+ return()
	61	+ }
	62	+}
	63	+
	64	+# If a data window is specified, get its corner coordinates
	65	+# and truncate the data by this window:
	66	+n <- length(x)
	67	+if(n!=length(y)) stop('data lengths do not match')
	68	+if(!is.null(rw)) {
	69	+ xmin <- rw[1]
	70	+ xmax <- rw[2]
	71	+ ymin <- rw[3]
	72	+ ymax <- rw[4]
	73	+ drop <- (1:n)[x<xmin\|x>xmax\|y<ymin\|y>ymax]
	74	+ if(length(drop)>0) {
	75	+ x <- x[-drop]
	76	+ y <- y[-drop]
	77	+ n <- length(x)
	78	+ }
	79	+}
	80	+
	81	+# If corners of the window are not specified, form them from
	82	+# the minimum and maximum of the data +/- 10%:
	83	+else {
	84	+ xmin <- min(x)
	85	+ xmax <- max(x)
	86	+ ymin <- min(y)
	87	+ ymax <- max(y)
	88	+ xdff <- xmax-xmin
	89	+ ydff <- ymax-ymin
	90	+ xmin <- xmin-0.1*xdff
	91	+ xmax <- xmax+0.1*xdff
	92	+ ymin <- ymin-0.1*ydff
	93	+ ymax <- ymax+0.1*ydff
	94	+ rw <- c(xmin,xmax,ymin,ymax)
	95	+}
	96	+
	97	+# Add the dummy points:
	98	+if(!is.null(dpl)) {
	99	+ dpts <- dumpts(x,y,dpl,rw)
	100	+ x <- dpts$x
	101	+ y <- dpts$y
	102	+}
	103	+
	104	+# Eliminate duplicate points:
	105	+iii <- !ind.dup(x,y,rw,frac)
	106	+ndm <- sum(iii[-(1:n)])
	107	+n <- sum(iii[1:n])
	108	+x <- x[iii]
	109	+y <- y[iii]
	110	+
	111	+# Make space for the total number of points (real and dummy) as
	112	+# well as 4 ideal points and 4 extra corner points which get used
	113	+# (only) by subroutines dirseg and dirout in the ``output'' process
	114	+# (returning a description of the triangulation after it has been
	115	+# calculated):
	116	+npd <- n + ndm
	117	+ntot <- npd + 4 # ntot includes the 4 ideal points but
	118	+ # but NOT the 4 extra corners
	119	+x <- c(rep(0,4),x,rep(0,4))
	120	+y <- c(rep(0,4),y,rep(0,4))
	121	+
	122	+# Set up fixed dimensioning constants:
	123	+ntdel <- 4*npd
	124	+ntdir <- 3*npd
	125	+
	126	+# Set up dimensioning constants which might need to be increased:
	127	+madj <- max(20,ceiling(3*sqrt(ntot)))
	128	+tadj <- (madj+1)*(ntot+4)
	129	+ndel <- madj*(madj+1)/2
	130	+tdel <- 6*ndel
	131	+ndir <- ndel
	132	+tdir <- 8*ndir
	133	+
	134	+# Call the master subroutine to do the work:
	135	+repeat {
	136	+ tmp <- .Fortran(
	137	+ 'master',
	138	+ x=as.double(x),
	139	+ y=as.double(y),
	140	+ sort=as.logical(sort),
	141	+ rw=as.double(rw),
	142	+ npd=as.integer(npd),
	143	+ ntot=as.integer(ntot),
	144	+ nadj=integer(tadj),
	145	+ madj=as.integer(madj),
	146	+ ind=integer(npd),
	147	+ tx=double(npd),
	148	+ ty=double(npd),
	149	+ ilist=integer(npd),
	150	+ eps=as.double(eps),
	151	+ delsgs=double(tdel),
	152	+ ndel=as.integer(ndel),
	153	+ delsum=double(ntdel),
	154	+ dirsgs=double(tdir),
	155	+ ndir=as.integer(ndir),
	156	+ dirsum=double(ntdir),
	157	+ nerror=integer(1),
	158	+ PACKAGE='deldir'
	159	+ )
	160	+
	161	+# Check for errors:
	162	+ nerror <- tmp$nerror
	163	+ if(nerror < 0) break
	164	+
	165	+ else {
	166	+ if(nerror==4) {
	167	+ cat('nerror =',nerror,'\n')
	168	+ cat('Increasing madj and trying again.\n')
	169	+ madj <- ceiling(1.2*madj)
	170	+ tadj <- (madj+1)*(ntot+4)
	171	+ ndel <- max(ndel,madj*(madj+1)/2)
	172	+ tdel <- 6*ndel
	173	+ ndir <- ndel
	174	+ tdir <- 8*ndir
	175	+ }
	176	+ else if(nerror==14\|nerror==15) {
	177	+ cat('nerror =',nerror,'\n')
	178	+ cat('Increasing ndel and ndir and trying again.\n')
	179	+ ndel <- ceiling(1.2*ndel)
	180	+ tdel <- 6*ndel
	181	+ ndir <- ndel
	182	+ tdir <- 8*ndir
	183	+ }
	184	+ else {
	185	+ cat('nerror =',nerror,'\n')
	186	+ return(invisible())
	187	+ }
	188	+ }
	189	+}
	190	+
	191	+# Collect up the results for return:
	192	+ndel <- tmp$ndel
	193	+delsgs <- round(t(as.matrix(matrix(tmp$delsgs,nrow=6)[,1:ndel])),digits)
	194	+delsum <- matrix(tmp$delsum,ncol=4)
	195	+del.area <- sum(delsum[,4])
	196	+delsum <- round(cbind(delsum,delsum[,4]/del.area),digits)
	197	+del.area <- round(del.area,digits)
	198	+ndir <- tmp$ndir
	199	+dirsgs <- round(t(as.matrix(matrix(tmp$dirsgs,nrow=8)[,1:ndir])),digits)
	200	+dirsgs <- as.data.frame(dirsgs)
	201	+dirsum <- matrix(tmp$dirsum,ncol=3)
	202	+dir.area <- sum(dirsum[,3])
	203	+dirsum <- round(cbind(dirsum,dirsum[,3]/dir.area),digits)
	204	+dir.area <- round(dir.area,digits)
	205	+allsum <- cbind(delsum,dirsum)
	206	+rw <- round(rw,digits)
	207	+
	208	+# Name the columns of the results:
	209	+dimnames(delsgs) <- list(NULL,c('x1','y1','x2','y2','ind1','ind2'))
	210	+names(dirsgs) <- c('x1','y1','x2','y2','ind1','ind2','bp1','bp2')
	211	+mode(dirsgs$bp1) <- 'logical'
	212	+mode(dirsgs$bp2) <- 'logical'
	213	+dimnames(allsum) <- list(NULL,c('x','y','n.tri','del.area','del.wts',
	214	+ 'n.tside','nbpt','dir.area','dir.wts'))
	215	+
	216	+# Aw' done!!!
	217	+rslt <- list(delsgs=delsgs,dirsgs=dirsgs,summary=allsum,n.data=n,
	218	+ n.dum=ndm,del.area=del.area,dir.area=dir.area,rw=rw)
	219	+class(rslt) <- 'deldir'
	220	+if(plotit) plot(rslt,...)
	221	+if(plotit) invisible(rslt) else rslt
	222	+}
	223	+dumpts <- function(x,y,dpl,rw) {
	224	+#
	225	+# Function dumpts to append a sequence of dummy points to the
	226	+# data points.
	227	+#
	228	+
	229	+ndm <- 0
	230	+xd <- NULL
	231	+yd <- NULL
	232	+xmin <- rw[1]
	233	+xmax <- rw[2]
	234	+ymin <- rw[3]
	235	+ymax <- rw[4]
	236	+
	237	+# Points on radii of circles emanating from data points:
	238	+if(!is.null(dpl$nrad)) {
	239	+ nrad <- dpl$nrad # Number of radii from each data point.
	240	+ nper <- dpl$nper # Number of dummy points per radius.
	241	+ fctr <- dpl$fctr # Length of each radius = fctr * mean
	242	+ # interpoint distance.
	243	+ lrad <- fctr*mnnd(x,y)/nper
	244	+ theta <- 2pi(1:nrad)/nrad
	245	+ cs <- cos(theta)
	246	+ sn <- sin(theta)
	247	+ xt <- c(lrad*(1:nper)%o%cs)
	248	+ yt <- c(lrad*(1:nper)%o%sn)
	249	+ xd <- c(outer(x,xt,'+'))
	250	+ yd <- c(outer(y,yt,'+'))
	251	+}
	252	+
	253	+# Ad hoc points passed over as part of dpl:
	254	+if(!is.null(dpl$x)) {
	255	+ xd <- c(xd,dpl$x)
	256	+ yd <- c(yd,dpl$y)
	257	+}
	258	+
	259	+# Delete dummy points outside the rectangular window.
	260	+ndm <- length(xd)
	261	+if(ndm >0) {
	262	+ drop <- (1:ndm)[xd<xmin\|xd>xmax\|yd<ymin\|yd>ymax]
	263	+ if(length(drop)>0) {
	264	+ xd <- xd[-drop]
	265	+ yd <- yd[-drop]
	266	+ }
	267	+}
	268	+
	269	+# Rectangular grid:
	270	+ndx <- dpl$ndx
	271	+okx <- !is.null(ndx) && ndx > 0
	272	+ndy <- dpl$ndy
	273	+oky <- !is.null(ndy) && ndy > 0
	274	+if(okx & oky) {
	275	+ xt <- if(ndx>1) seq(xmin,xmax,length=ndx) else 0.5*(xmin+xmax)
	276	+ yt <- if(ndy>1) seq(ymin,ymax,length=ndy) else 0.5*(ymin+ymax)
	277	+ xy <- expand.grid(x=xt,y=yt)
	278	+ xd <- c(xd,xy$x)
	279	+ yd <- c(yd,xy$y)
	280	+}
	281	+
	282	+ndm <- length(xd)
	283	+list(x=c(x,xd),y=c(y,yd),ndm=ndm)
	284	+}
	285	+get.cnrind <- function(x,y,rw) {
	286	+x.crnrs <- rw[c(1,2,2,1)]
	287	+y.crnrs <- rw[c(3,3,4,4)]
	288	+M1 <- outer(x,x.crnrs,function(a,b){(a-b)^2})
	289	+M2 <- outer(y,y.crnrs,function(a,b){(a-b)^2})
	290	+MM <- M1 + M2
	291	+apply(MM,2,which.min)
	292	+}
	293	+ind.dup <- function(x,y,rw=NULL,frac=0.0001) {
	294	+#
	295	+# Function ind.dup to calculate the indices of data pairs
	296	+# which duplicate earlier ones. (Returns a logical vector;
	297	+# true for such indices, false for the rest.)
	298	+#
	299	+
	300	+if(is.null(rw)) rw <- c(0,1,0,1)
	301	+n <- length(x)
	302	+rslt <- .Fortran(
	303	+ 'inddup',
	304	+ x=as.double(x),
	305	+ y=as.double(y),
	306	+ n=as.integer(n),
	307	+ rw=as.double(rw),
	308	+ frac=as.double(frac),
	309	+ dup=logical(n),
	310	+ PACKAGE='deldir'
	311	+ )
	312	+
	313	+rslt$dup
	314	+}
	315	+mid.in <- function(x,y,rx,ry) {
	316	+xm <- 0.5*(x[1]+x[2])
	317	+ym <- 0.5*(y[1]+y[2])
	318	+(rx[1] < xm & xm < rx[2] & ry[1] < ym & ym < ry[2])
	319	+}
	320	+mnnd <- function(x,y) {
	321	+#
	322	+# Function mnnd to calculate the mean nearest neighbour distance
	323	+# between the points whose coordinates are stored in x and y.
	324	+#
	325	+
	326	+n <- length(x)
	327	+if(n!=length(y)) stop('data lengths do not match')
	328	+dmb <- (max(x)-min(x))2 + (max(y)-min(y))2
	329	+
	330	+.Fortran(
	331	+ "mnnd",
	332	+ x=as.double(x),
	333	+ y=as.double(y),
	334	+ n=as.integer(n),
	335	+ dmb=as.double(dmb),
	336	+ d=double(1),
	337	+ PACKAGE='deldir'
	338	+ )$d
	339	+}
	340	+plot.deldir <- function(x,add=FALSE,wlines=c('both','triang','tess'),
	341	+ wpoints=c('both','real','dummy','none'),
	342	+ number=FALSE,cex=1,nex=1,col=NULL,lty=NULL,
	343	+ pch=NULL,xlim=NULL,ylim=NULL,xlab='x',ylab='y',...)
	344	+{
	345	+#
	346	+# Function plot.deldir to produce a plot of the Delaunay triangulation
	347	+# and Dirichlet tesselation of a point set, as produced by the
	348	+# function deldir().
	349	+#
	350	+
	351	+wlines <- match.arg(wlines)
	352	+wpoints <- match.arg(wpoints)
	353	+
	354	+if(is.null(class(x)) \|\| class(x)!='deldir') {
	355	+ cat('Argument is not of class deldir.\n')
	356	+ return(invisible())
	357	+}
	358	+
	359	+col <- if(is.null(col)) c(1,1,1,1,1) else rep(col,length.out=5)
	360	+lty <- if(is.null(lty)) 1:2 else rep(lty,length.out=2)
	361	+pch <- if(is.null(pch)) 1:2 else rep(pch,length.out=2)
	362	+
	363	+plot.del <- switch(wlines,both=TRUE,triang=TRUE,tess=FALSE)
	364	+plot.dir <- switch(wlines,both=TRUE,triang=FALSE,tess=TRUE)
	365	+plot.rl <- switch(wpoints,both=TRUE,real=TRUE,dummy=FALSE,none=FALSE)
	366	+plot.dum <- switch(wpoints,both=TRUE,real=FALSE,dummy=TRUE,none=FALSE)
	367	+
	368	+delsgs <- x$delsgs
	369	+dirsgs <- x$dirsgs
	370	+n <- x$n.data
	371	+rw <- x$rw
	372	+
	373	+if(plot.del) {
	374	+ x1<-delsgs[,1]
	375	+ y1<-delsgs[,2]
	376	+ x2<-delsgs[,3]
	377	+ y2<-delsgs[,4]
	378	+}
	379	+
	380	+if(plot.dir) {
	381	+ u1<-dirsgs[,1]
	382	+ v1<-dirsgs[,2]
	383	+ u2<-dirsgs[,3]
	384	+ v2<-dirsgs[,4]
	385	+}
	386	+
	387	+X<-x$summary[,1]
	388	+Y<-x$summary[,2]
	389	+
	390	+if(!add) {
	391	+ pty.save <- par()$pty
	392	+ on.exit(par(pty=pty.save))
	393	+ par(pty='s')
	394	+ if(is.null(xlim)) xlim <- rw[1:2]
	395	+ if(is.null(ylim)) ylim <- rw[3:4]
	396	+ plot(0,0,type='n',xlim=xlim,ylim=ylim,
	397	+ xlab=xlab,ylab=ylab,axes=FALSE,...)
	398	+ axis(side=1)
	399	+ axis(side=2)
	400	+}
	401	+
	402	+if(plot.del) segments(x1,y1,x2,y2,col=col[1],lty=lty[1],...)
	403	+if(plot.dir) segments(u1,v1,u2,v2,col=col[2],lty=lty[2],...)
	404	+if(plot.rl) {
	405	+ x.real <- X[1:n]
	406	+ y.real <- Y[1:n]
	407	+ points(x.real,y.real,pch=pch[1],col=col[3],cex=cex,...)
	408	+}
	409	+if(plot.dum) {
	410	+ x.dumm <- X[-(1:n)]
	411	+ y.dumm <- Y[-(1:n)]
	412	+ points(x.dumm,y.dumm,pch=pch[2],col=col[4],cex=cex,...)
	413	+}
	414	+if(number) {
	415	+ xoff <-0.02*diff(range(X))
	416	+ yoff <-0.02*diff(range(Y))
	417	+ text(X+xoff,Y+yoff,1:length(X),cex=nex,col=col[5],...)
	418	+}
	419	+invisible()
	420	+}
	421	+`plot.tile.list` <-
	422	+function (x, verbose = FALSE, close = FALSE, pch = 1, polycol = NA,
	423	+ showpoints = TRUE, asp = 1, ...)
	424	+{
	425	+ object <- x
	426	+ if (!inherits(object, "tile.list"))
	427	+ stop("Argument \"object\" is not of class tile.list.\n")
	428	+ n <- length(object)
	429	+ x.all <- unlist(lapply(object, function(w) {
	430	+ c(w$pt[1], w$x)
	431	+ }))
	432	+ y.all <- unlist(lapply(object, function(w) {
	433	+ c(w$pt[2], w$y)
	434	+ }))
	435	+ x.pts <- unlist(lapply(object, function(w) {
	436	+ w$pt[1]
	437	+ }))
	438	+ y.pts <- unlist(lapply(object, function(w) {
	439	+ w$pt[2]
	440	+ }))
	441	+ rx <- range(x.all)
	442	+ ry <- range(y.all)
	443	+ plot(x.all, y.all, type = "n", asp = asp, xlab = "x", ylab = "y")
	444	+ polycol <- apply(col2rgb(polycol,TRUE),2,
	445	+ function(x){do.call(rgb,as.list(x/255))})
	446	+ polycol <- rep(polycol, length = length(object))
	447	+ hexbla <- do.call(rgb,as.list(col2rgb("black",TRUE)/255))
	448	+ hexwhi <- do.call(rgb,as.list(col2rgb("white",TRUE)/255))
	449	+ ptcol <- ifelse(polycol == hexbla,hexwhi,hexbla)
	450	+ lnwid <- ifelse(polycol == hexbla, 2, 1)
	451	+ for (i in 1:n) {
	452	+ inner <- !any(object[[i]]$bp)
	453	+ if (close \| inner)
	454	+ polygon(object[[i]], col = polycol[i], border = ptcol[i],
	455	+ lwd = lnwid[i])
	456	+ else {
	457	+ x <- object[[i]]$x
	458	+ y <- object[[i]]$y
	459	+ bp <- object[[i]]$bp
	460	+ ni <- length(x)
	461	+ for (j in 1:ni) {
	462	+ jnext <- if (j < ni)
	463	+ j + 1
	464	+ else 1
	465	+ do.it <- mid.in(x[c(j, jnext)], y[c(j, jnext)],
	466	+ rx, ry)
	467	+ if (do.it)
	468	+ segments(x[j], y[j], x[jnext], y[jnext], col = ptcol[i],
	469	+ lwd = lnwid[i])
	470	+ }
	471	+ }
	472	+ if (verbose & showpoints)
	473	+ points(object[[i]]$pt[1], object[[i]]$pt[2], pch = pch,
	474	+ col = ptcol[i])
	475	+ if (verbose & i < n)
	476	+ readline("Go? ")
	477	+ }
	478	+ if (showpoints)
	479	+ points(x.pts, y.pts, pch = pch, col = ptcol)
	480	+ invisible()
	481	+}
	482	+tile.list <- function (object)
	483	+{
	484	+ if (!inherits(object, "deldir"))
	485	+ stop("Argument \"object\" is not of class deldir.\n")
	486	+ rw <- object$rw
	487	+ x.crnrs <- rw[c(1,2,2,1)]
	488	+ y.crnrs <- rw[c(3,3,4,4)]
	489	+ ddd <- object$dirsgs
	490	+ sss <- object$summary
	491	+ npts <- nrow(sss)
	492	+ x <- sss[,"x"]
	493	+ y <- sss[,"y"]
	494	+ i.crnr <- get.cnrind(x,y,rw)
	495	+ rslt <- list()
	496	+ for (i in 1:npts) {
	497	+ m <- as.matrix(rbind(ddd[ddd$ind1 == i, 1:4], ddd[ddd$ind2 == i, 1:4]))
	498	+ bp1 <- c(ddd[ddd$ind1 == i, 7], ddd[ddd$ind2 == i, 7])
	499	+ bp2 <- c(ddd[ddd$ind1 == i, 8], ddd[ddd$ind2 == i, 8])
	500	+ m1 <- cbind(m[, 1:2, drop = FALSE], 0 + bp1)
	501	+ m2 <- cbind(m[, 3:4, drop = FALSE], 0 + bp2)
	502	+ m <- rbind(m1, m2)
	503	+ pt <- sss[i, 1:2]
	504	+ theta <- atan2(m[, 2] - pt[2], m[, 1] - pt[1])
	505	+ theta <- ifelse(theta > 0, theta, theta + 2 * pi)
	506	+ theta.0 <- sort(unique(theta))
	507	+ mm <- m[match(theta.0, theta), ]
	508	+ xx <- mm[, 1]
	509	+ yy <- mm[, 2]
	510	+ bp <- as.logical(mm[, 3])
	511	+# Add corner points if necessary:
	512	+ ii <- i.crnr%in%i
	513	+ xx <- c(xx,x.crnrs[ii])
	514	+ yy <- c(yy,y.crnrs[ii])
	515	+ bp <- c(bp,rep(TRUE,sum(ii)))
	516	+ rslt[[i]] <- acw(list(pt=pt, x = xx, y = yy, bp = bp))
	517	+ }
	518	+ class(rslt) <- "tile.list"
	519	+ rslt
	520	+}
...	...
...	...	@@ -0,0 +1,186 @@
	1	+
	2	+ Version date: 21 February 2002.
	3	+ This version is simply an adaptation of the Splus version
	4	+ of the package to R.
	5	+
	6	+ *****************************
	7	+
	8	+ Version date: 14 February 2002.
	9	+ This version supercedes the version dated 24 April 1999.
	10	+
	11	+ *****************************
	12	+
	13	+ The changes from the version dated 24 April 1999 to the
	14	+ version dated 14 February 2002 were:
	15	+
	16	+ A bug in the procedure for eliminating duplicated points was
	17	+ fixed. Thanks go to Dr. Berwin Turlach of the Department of
	18	+ Maths and Stats at the University of Western Australia, for
	19	+ spotting this bug.
	20	+
	21	+ *****************************
	22	+
	23	+ The changes from the version dated 26 October 1998 to the
	24	+ version dated 24 April 1999 were:
	25	+
	26	+ (1) The function mipd(), stored in mipd.sf, and the
	27	+ corresponding Fortran subroutine mipd, stored in mipd.r, have
	28	+ been replaced by mnnd() in mnnd.sf and mnnd in mnnd.r. The
	29	+ function mipd calculated the mean interpoint distance, to be
	30	+ used in constructing dummy point structures of a certain
	31	+ type. After some reflection it became apparent that the mean
	32	+ interpoint distance was much too large for the intended
	33	+ purpose, and that a more appropriate value was the ``mean
	34	+ nearest neighbour distance'' which is calculated by the new
	35	+ function. This new value is now used in constructing dummy
	36	+ point structures.
	37	+
	38	+ Note that the operative result is that the resulting dummy
	39	+ point structures contain many more points than before. The
	40	+ old value caused large numbers of the dummy points to fall
	41	+ outside the data window and therefore to be clipped.
	42	+
	43	+ *****************************
	44	+
	45	+ The changes from the version dated 6 December 1996 to the
	46	+ version dated 26 October 1998 were:
	47	+
	48	+ (1) A ratfor/Fortran routine named ``inside'' has been
	49	+ renamed ``dldins'' to avoid conflict with a name built in to
	50	+ some versions of Splus.
	51	+
	52	+ (2) Some minor corrections have been made to dangerous
	53	+ infelicities in a piece of the ratfor/Fortran code.
	54	+
	55	+ (3) The dynamic loading procedure has been changed to use
	56	+ dyn.load.shared so that the package is easily usable on
	57	+ IRIX systems as well as under SunOS/Solaris.
	58	+
	59	+ (4) The package has been adjusted slightly so that it
	60	+ can easily be installed as a section of a library. In
	61	+ particular, the dynamic loading is now done by the
	62	+ .First.lib() function rather than from within deldir()
	63	+ itself; reference to an environment variable DYN_LOAD_LIB
	64	+ is no longer needed.
	65	+
	66	+ *****************************
	67	+
	68	+ This package computes and plots the Dirichlet (Voronoi) tesselation
	69	+ and the Delaunay triangulation of a set of of data points and
	70	+ possibly a superimposed ``grid'' of dummy points.
	71	+
	72	+ The tesselation is constructed with respect to the whole plane
	73	+ by suspending it from ideal points at infinity.
	74	+
	75	+ ORIGINALLY PROGRAMMED BY: Rolf Turner in 1987/88, while with the
	76	+ Division of Mathematics and Statistics, CSIRO, Sydney, Australia.
	77	+ Re-programmed by Rolf Turner to adapt the implementation from a
	78	+ stand-alone Fortran program to an S function, while visiting the
	79	+ University of Western Australia, May 1995. Further revised
	80	+ December 1996, October 1998, April 1999, and February 2002.
	81	+ Adapted to an R package 21 February 2002.
	82	+
	83	+ Current address of the author:
	84	+ Department of Mathematics and Statistics,
	85	+ University of New Brunswick,
	86	+ P.O. Box 4400, Fredericton, New Brunswick,
	87	+ Canada E3B 5A3
	88	+ Email:
	89	+ rolf@math.unb.ca
	90	+
	91	+ The author gratefully acknowledges the contributions, assistance,
	92	+ and guidance of Mark Berman, of D.M.S., CSIRO, in collaboration with
	93	+ whom this project was originally undertaken. The author also
	94	+ acknowledges much useful advice from Adrian Baddeley, formerly of
	95	+ D.M.S. CSIRO (now Professor of Statistics at the University of
	96	+ Western Australia). Daryl Tingley of the Department of Mathematics
	97	+ and Statistics, University of New Brunswick provided some helpful
	98	+ insight. Special thanks are extended to Alan Johnson, of the Alaska
	99	+ Fisheries Science Centre, who supplied two data sets which were
	100	+ extremely valuable in tracking down some errors in the code.
	101	+
	102	+ Don MacQueen, of Lawrence Livermore National Lab, wrote an Splus
	103	+ driver function for the old stand-alone version of this software.
	104	+ That driver, which was available on Statlib, is now deprecated in
	105	+ favour of this current package. Don also collaborated in the
	106	+ preparation of this current package.
	107	+
	108	+ Bill Dunlap of MathSoft Inc. tracked down a bug which was making
	109	+ the deldir() function crash on some systems, and pointed out some
	110	+ other improvements to be made.
	111	+
	112	+ Berwin Turlach of the Department of Maths and Stats at the
	113	+ University of Western Australia pointed out a bug in the procedure
	114	+ for eliminating duplicated points.
	115	+
	116	+ *****************************
	117	+
	118	+ The man directory, contains, in addition to the R documentation
	119	+ files deldir.Rd and plot.deldir.Rd:
	120	+
	121	+ (a) This READ_ME file.
	122	+ (b) A file err.list, containing a list of meanings of possible
	123	+ error numbers which could be returned. NONE of these
	124	+ errors should ever actually happen except for errors 4, 14,
	125	+ and 15. These relate to insufficient dimensioning, and
	126	+ if they occur, the driver increases the dimensions and
	127	+ tries again (informing you of this fact).
	128	+ (c) A file ex.out containing a printout of the object returned
	129	+ by running the example given in the help file for deldir.
	130	+
	131	+ The src directory contains many, many *.f (Fortran) files, which
	132	+ get compiled and dynamically loaded.
	133	+
	134	+ The Fortran code is ponderous --- it was automatically generated
	135	+ from Ratfor code, which was pretty ponderous to start with. It is
	136	+ quite possibly very kludgy aw well --- i.e. a good programmer could
	137	+ make it *much* more efficient I'm sure. It contains all sorts
	138	+ of checking for anomalous situations which probably can/will never
	139	+ occur. These checks basically reflect my pessimism and fervent
	140	+ belief in Murphy's Law.
	141	+
	142	+ The program was also designed with a particular application in mind,
	143	+ in which we wished to superimpose a grid of dummy points onto the
	144	+ actual data points which we were triangulating. This fact adds
	145	+ slightly to the complication of the code.
	146	+
	147	+ *****************************
	148	+
	149	+ Here follows a brief description of the package:
	150	+
	151	+ (1) The function deldir computes the Delaunay Triangulation (and
	152	+ hence the Dirichlet Tesselation) of a planar point set according to
	153	+ the second (iterative) algorithm of Lee and Schacter, International
	154	+ Journal of Computer and Information Sciences, Vol. 9, No. 3, 1980,
	155	+ pages 219 to 242.
	156	+
	157	+ The tesselation/triangulation is made to be
	158	+
	159	+ ** with respect to the whole plane **
	160	+
	161	+ by `suspending' it from `ideal' points (-Inf,-Inf), (Inf,-Inf)
	162	+ (Inf,Inf), and (-Inf,Inf).
	163	+
	164	+ (2) The tesselation/triangulation is also enclosed in a finite
	165	+ rectangle with corners
	166	+
	167	+ (xmin,ymax) * ------------------------ * (xmax,ymax)
	168	+ \| \|
	169	+ \| \|
	170	+ \| \|
	171	+ \| \|
	172	+ \| \|
	173	+ (xmin,ymin) * ------------------------ * (xmax,ymin)
	174	+
	175	+ The boundaries of this rectangle truncate some Dirichlet tiles, in
	176	+ particular any infinite ones. This rectangle is referred to
	177	+ elsewhere as `the' rectangular window.
	178	+ ===
	179	+
	180	+ (2) The function plot.deldir is a method for plot. I.e. it may be
	181	+ invoked simply by typing ``plot(x)'' provided that ``x'' is an
	182	+ object of class ``deldir'' (as produced by the function deldir).
	183	+ The plot (by default) consists of the edges of the Delaunay
	184	+ triangles (solid lines) and the edges of the Dirichlet tiles (dotted
	185	+ lines). By default the real data points are indicated by circles,
	186	+ and the dummy points are indicated by triangles.
...	...