#!/bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #!/bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# READ_ME
# README
# First.lib.sf
# deldir.sf
# dumpts.sf
# mnnd.sf
# ind.dup.sf
# plot.deldir.sf
# print.deldir.sf
# deldir.d
# plot.deldir.d
# print.deldir.d
# Makefile
# install.help
# acchk.r
# addpt.r
# adjchk.r
# binsrt.r
# circen.r
# cross.r
# delet.r
# delet1.r
# delout.r
# delseg.r
# dirout.r
# dirseg.r
# inddup.r
# initad.r
# inside.r
# insrt.r
# insrt1.r
# locn.r
# master.r
# mnnd.r
# pred.r
# qtest.r
# qtest1.r
# stoke.r
# succ.r
# swap.r
# testeq.r
# triar.r
# trifnd.r
# err.list
# ex.out
# This archive created: Thu Feb 14 10:25:13 2002
export PATH; PATH=/bin:$PATH
if test -f 'READ_ME'
then
echo shar: over-writing existing file "'READ_ME'"
fi
cat << \SHAR_EOF > 'READ_ME'
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.
*****************************
LIBRARY SECTION OF FUNCTIONS TO COMPUTE AND PLOT THE DIRICHLET
TESSELLATION AND DELAUNAY TRIANGULATION OF A SET OF DATA POINTS AND
POSSIBLY A SUPERIMPOSED GRID OF DUMMY POINTS.
THE TRIANGULATION IS CONSTRUCTED WITH REPECT TO THE WHOLE PLANE
BY SUSPENDING IT FROM IDEAL POINTS AT INFINITY.
Copyright (C) 2002 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, October 1998, April 1999, and 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.
*****************************
This archive contains:
(a) This READ_ME file.
(b) A brief file named README (no ``_'') as required to
form a library section.
(c) The S code of the functions .First.lib(), deldir(),
dumpts(), mnnd(), ind.dup(), plot.deldir(), and
print.deldir(), in files First.lib.sf, deldir.sf,
dumpts.sf, mnnd.sf(), ind.dup.sf, plot.deldir.sf, and
print.deldir.sf
(d) Documentation files deldir.d, plot.deldir.d, and print.deldir.d
(e) A Makefile.
(f) A shell script install.help, used by the Makefile.
(g) Many, many *.r files of ratfor subroutines.
(h) 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).
(i) A file ex.out containing a printout of the object returned
by running the example given in the help file for deldir.
The ratfor code is ponderous and quite possibly kludgy -- 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 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).
It 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). The data points are distinguished (by default) by the real
data points being indicated by circles, and the dummy points being
indicated by triangles.
(3) The function print.deldir (a "method" for print) is included to
avoid a minor inconvenience. The output of the function deldir is a
list one of whose components is a data frame. The print.list
function has the default "quote=T" which puts quote-marks into the
printed version of the data frame, which looks messy. The function
print.deldir averts this annoying feature.
(4) To get started:
(a) Place all the files from this archive (unpack this archive)
in a directory which is to form a ``section'' of a library.
(Either your own personal library, or the system library if
you are installing this package for use by multiple users on the
system.)
(b) Modify the `Makefile' for compatibility with your local
system. In particular you may need to change the definition
of the SPLUS macro, depending on the name and location of the
Splus executable on your system.
(c) Type ``make'' or ``make load'' which creates a ``shared
library'' object file called ``deldirld.so'' to be dynamically
loaded using dyn.load.shared.
(d) Make sure that the shell script ``install.help'' is
executable (e.g. execute ``chmod +x install.help'').
(e) Type ``make install'' to install the S functions .First.lib,
deldir, dumpts, mnnd, ind.dup, and plot.deldir in .Data, and
the help files for deldir, plot.deldir, and print.deldir in
.Data/.Help
(f) You can remove unnecessary *.o files by typing ``make clean''.
(g) Start an S session in whatever directory you wish to work
in. Issue the approprate library() command, of the form
library(lib.loc=library_name,section=section_name)
For example if the Delaunay package is installed as a
section in a directory
/home/users/melvin/mylib/delaunay
where ``mylib'' is your personal library, then the command
would be
library(lib.loc="/home/users/melvin/mylib",section="delaunay")
To avoid having to repeatedly type the longish value of the
lib.loc string you can assign a value to lib.loc in frame 0:
assign("lib.loc","/home/users/melvin/mylib",where=0)
and this in turn can be done in your .First function so that
it happens automatically whenever you start Splus in the
given directory. Then you can just type
library("delaunay")
to access this library.
(h) Type
help(deldir)
and/or
help(plot.deldir)
to obtain instructions on how to invoke these functions.
SHAR_EOF
if test -f 'README'
then
echo shar: over-writing existing file "'README'"
fi
cat << \SHAR_EOF > 'README'
Functions to implement Delaunay triangulation and Dirichlet (Voronoi)
tessellation of planar point sets
=====================================================================
functions: brief descriptions:
deldir calculate the triangulation and tessellation
dumpts create a set of dummy points if desired
mnnd auxilliary --- mean nearest neighbour distance
ind.dup auxilliary --- find indices of duplicate points
plot.deldir plot the triangulation and tessellation
print.deldir print the output of deldir
SHAR_EOF
if test -f 'First.lib.sf'
then
echo shar: over-writing existing file "'First.lib.sf'"
fi
cat << \SHAR_EOF > 'First.lib.sf'
.First.lib <- function(lib.loc,section) {
path <- paste(lib.loc,section,"deldirld.so",sep="/")
if(!is.loaded(symbol.For('acchk')))
dyn.load.shared(path)
invisible()
}
SHAR_EOF
if test -f 'deldir.sf'
then
echo shar: over-writing existing file "'deldir.sf'"
fi
cat << \SHAR_EOF > 'deldir.sf'
deldir <- function(x,y,dpl=NULL,rw=NULL,eps=1e-9,frac=1e-4,
sort=T,plotit=F,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)[xxmax|yymax]
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)
)
# 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
}
SHAR_EOF
if test -f 'dumpts.sf'
then
echo shar: over-writing existing file "'dumpts.sf'"
fi
cat << \SHAR_EOF > 'dumpts.sf'
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)[xdxmax|ydymax]
if(length(drop)>0) {
xd <- xd[-drop]
yd <- yd[-drop]
}
}
# Rectangular grid:
if(!is.null(dpl$ndx)) {
ndx <- dpl$ndx
ndy <- dpl$ndy
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)
xd <- c(xd,rep(xt,ndy))
yd <- c(yd,rep(yt,rep(ndx,ndy)))
}
ndm <- length(xd)
list(x=c(x,xd),y=c(y,yd),ndm=ndm)
}
SHAR_EOF
if test -f 'mnnd.sf'
then
echo shar: over-writing existing file "'mnnd.sf'"
fi
cat << \SHAR_EOF > 'mnnd.sf'
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)
)$d
}
SHAR_EOF
if test -f 'ind.dup.sf'
then
echo shar: over-writing existing file "'ind.dup.sf'"
fi
cat << \SHAR_EOF > 'ind.dup.sf'
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)
)
rslt$dup
}
SHAR_EOF
if test -f 'plot.deldir.sf'
then
echo shar: over-writing existing file "'plot.deldir.sf'"
fi
cat << \SHAR_EOF > 'plot.deldir.sf'
plot.deldir <- function(object,add=F,wlines=c('both','triang','tess'),
wpoints=c('both','real','dummy','none'),
number=F,cex=0.5,nex=0.75,
col=NULL,lty=NULL,pch=NULL,xlim=NULL,ylim=NULL)
{
#
# 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(object)) || class(object)!='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=T,triang=T,tess=F)
plot.dir <- switch(wlines,both=T,triang=F,tess=T)
plot.rl <- switch(wpoints,both=T,real=T,dummy=F,none=F)
plot.dum <- switch(wpoints,both=T,real=F,dummy=T,none=F)
delsgs <- object$delsgs
dirsgs <- object$dirsgs
n <- object$n.data
rw <- object$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<-object$summary[,1]
y<-object$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='',ylab='',axes=F)
axis(side=1)
axis(side=2)
axes(xlab='x',ylab='y',axes=F)
}
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()
}
SHAR_EOF
if test -f 'print.deldir.sf'
then
echo shar: over-writing existing file "'print.deldir.sf'"
fi
cat << \SHAR_EOF > 'print.deldir.sf'
print.deldir <- function(x) {
#
# Function print.deldir --- a ``method'' for print.
# Needed only to keep quotes from being put into the output
# in printing the ``dirseg'' component of the list returned
# by deldir().
#
print.list(x, quote = F)
}
SHAR_EOF
if test -f 'deldir.d'
then
echo shar: over-writing existing file "'deldir.d'"
fi
cat << \SHAR_EOF > 'deldir.d'
.BG
.FN deldir
.TL
Construct the Delaunay triangulation and the Dirichlet
(or Voronoi) tessellation of a planar point set.
.DN
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.
.CS
deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, frac=0.0001,
sort=T, plot=F, digits=6, ...)
.RA
.AG 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".
.OA
.AG 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:
.PP
ndx: The x-dimension of a rectangular grid; if either ndx or ndy is null,
no grid is constructed.
.PP
ndy: The y-dimension of the aforementioned rectangular grid.
.PP
nrad: The number of radii or "spokes", emanating from each data point,
along which dummy points are to be added.
.PP
nper: The number of dummy points per spoke.
.PP
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().)
.PP
x: A vector of x-coordinates of "ad hoc" dummy points
.PP
y: A vector of the corresponding y-coordinates of "ad hoc" dummy points
.PP
.AG 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.
.AG 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.
.AG 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).
.AG 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.
.AG plot
Logical argument; if TRUE a plot of the triangulation and tessellation
is produced; the default is FALSE.
.AG digits
The number of decimal places to which all numeric values in the
returned list should be rounded. Defaults to 6.
.AG ...
Auxilliary arguments add, wlines, wpoints, number, nex, col, lty,
pch, xlim, and ylim may be passed to plot.deldir through "..." if
plot=T.
.RT
A list (of class `deldir'), invisible if plot=T, with components:
.RC 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.
.RC 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.
.RC 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.
.RC 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.
.RC 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.
.RC 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".
.RC dir.area
the area of the rectangular window enclosing the points being triangulated,
as formed by summing the "dir.area" column of "summary".
.RC rw
the specification of the corners of the rectangular window enclosing
the data, in the order (xmin, xmax, ymin, ymax).
.SH NOTE
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.
.SE
If plot==T a plot of the triangulation and/or tessellation is produced
or added to an existing plot.
.DT
The code to effect this implementation of the Lee-Schacter algorithm
was originally written in 1987/88 by Rolf Turner, while with the
Division of Mathematics and Statistics, CSIRO, Sydney, Australia. It
was re-written to adapt the implementation from a stand-alone Fortran
program to an Splus function, 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 this current package. Don also collaborated in the
preparation of this current package.
.SH REFERENCES
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.
.SA
plot.deldir
.EX
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)
try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=T,wl='tr')
# Plots the triangulation which was created (but not the tesselation).
.KW
.WR
SHAR_EOF
if test -f 'plot.deldir.d'
then
echo shar: over-writing existing file "'plot.deldir.d'"
fi
cat << \SHAR_EOF > 'plot.deldir.d'
.BG
.FN plot.deldir
.TL
Produce a plot of the Delaunay triangulation and Dirichlet (Voronoi)
tesselation of a planar point set, as constructed by the function deldir.
.DN
This is a method for plot.
.CS
plot.deldir(object, add=F, wlines='both', wpoints='both',
number=F, cex=0.5, nex=0.75, col=NULL, lty=NULL,
pch=NULL, xlim=NULL, ylim=NULL)
plot(object, add=F, wlines='both', wpoints='both',
number=F, cex=0.5, nex=0.75, col=NULL, lty=NULL,
pch=NULL, xlim=NULL, ylim=NULL)
.RA
.AG object
An object of class "deldir" as constructed by the function deldir.
.OA
.AG add
logical argument; should the plot be added to an existing plot?
.AG 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) ?
.AG 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')?
.AG 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).
.AG cex
The value of the character expansion argument cex to be used
with the plotting symbols for plotting the points.
.AG 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=T.
.AG 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.)
.AG 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.)
.AG 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.)
.AG xlim
the limits on the x-axis. Defaults to rw[1:2] where rw is the
rectangular window specification returned by deldir().
.AG ylim
the limits on the y-axis. Defaults to rw[3:4] where rw is the
rectangular window specification returned by deldir().
.SE
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).
.DT
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).
.SA
deldir
.EX
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=T,add=T,wl='te',
col=c(1,1,2,3,4),num=T)
# Plots the tesselation, but does not save the results.
try <- deldir(x,y,list(ndx=2,ndy=2),c(0,10,0,10),plot=T,wl='tr',
wp='n')
# Plots the triangulation, but not the points, and saves the returned
structure.
.KW
.WR
SHAR_EOF
if test -f 'print.deldir.d'
then
echo shar: over-writing existing file "'print.deldir.d'"
fi
cat << \SHAR_EOF > 'print.deldir.d'
.BG
.FN print.deldir
.TL
print an object of class "deldir" as produced by the function
deldir
.DN
a method for print. Needed only to keep quotes from being put into
the output in printing the ``dirseg'' component of objects of class
"deldir"
.CS
print.deldir(x)
print(x)
.RA
.AG x
an object of class "deldir" as returned by the funtion deldir.
.RT
x, with the invisible flag set to prevent reprinting.
.SE
x is printed (component by component).
.SA
deldir, plot.deldir
.EX
try <- deldir(x,y,dp.str=list(ndx=2,ndy=2))
print(try)
.KW
.WR
SHAR_EOF
if test -f 'Makefile'
then
echo shar: over-writing existing file "'Makefile'"
fi
cat << \SHAR_EOF > 'Makefile'
SPLUS = /usr/local/bin/splus
OBJ = acchk.o addpt.o adjchk.o binsrt.o circen.o cross.o delet.o \
delet1.o delout.o delseg.o dirout.o dirseg.o inddup.o initad.o \
inside.o insrt.o insrt1.o locn.o master.o mnnd.o pred.o qtest.o \
qtest1.o stoke.o succ.o swap.o testeq.o triar.o trifnd.o
SFUNS = First.lib.sf deldir.sf dumpts.sf ind.dup.sf mnnd.sf \
plot.deldir.sf print.deldir.sf
HELPS = deldir.d plot.deldir.d print.deldir.d
RM = -rm
load: deldirld.so
deldirld.so: $(OBJ)
$(SPLUS) SHLIB -o deldirld.so $(OBJ)
install: sfuns help
sfuns: $(SFUNS) .Data
cat $(SFUNS) | $(SPLUS)
help: $(HELPS) .Data/.Help
install.help $(HELPS)
.Data:
-mkdir $@
.Data/.Help: .Data
-mkdir $@
clean:
$(RM) -f $(OBJ) core
SHAR_EOF
if test -f 'install.help'
then
echo shar: over-writing existing file "'install.help'"
fi
cat << \SHAR_EOF > 'install.help'
#! /bin/csh
foreach help ($argv)
cp $help .Data/.Help/$help:r
end
SHAR_EOF
chmod +x 'install.help'
if test -f 'acchk.r'
then
echo shar: over-writing existing file "'acchk.r'"
fi
cat << \SHAR_EOF > 'acchk.r'
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
SHAR_EOF
if test -f 'addpt.r'
then
echo shar: over-writing existing file "'addpt.r'"
fi
cat << \SHAR_EOF > 'addpt.r'
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
SHAR_EOF
if test -f 'adjchk.r'
then
echo shar: over-writing existing file "'adjchk.r'"
fi
cat << \SHAR_EOF > 'adjchk.r'
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
SHAR_EOF
if test -f 'binsrt.r'
then
echo shar: over-writing existing file "'binsrt.r'"
fi
cat << \SHAR_EOF > 'binsrt.r'
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
SHAR_EOF
if test -f 'circen.r'
then
echo shar: over-writing existing file "'circen.r'"
fi
cat << \SHAR_EOF > 'circen.r'
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
SHAR_EOF
if test -f 'cross.r'
then
echo shar: over-writing existing file "'cross.r'"
fi
cat << \SHAR_EOF > 'cross.r'
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
SHAR_EOF
if test -f 'delet.r'
then
echo shar: over-writing existing file "'delet.r'"
fi
cat << \SHAR_EOF > 'delet.r'
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
SHAR_EOF
if test -f 'delet1.r'
then
echo shar: over-writing existing file "'delet1.r'"
fi
cat << \SHAR_EOF > 'delet1.r'
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
SHAR_EOF
if test -f 'delout.r'
then
echo shar: over-writing existing file "'delout.r'"
fi
cat << \SHAR_EOF > 'delout.r'
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
SHAR_EOF
if test -f 'delseg.r'
then
echo shar: over-writing existing file "'delseg.r'"
fi
cat << \SHAR_EOF > 'delseg.r'
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
SHAR_EOF
if test -f 'dirout.r'
then
echo shar: over-writing existing file "'dirout.r'"
fi
cat << \SHAR_EOF > 'dirout.r'
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 'dirseg.r'
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 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
SHAR_EOF
if test -f 'inddup.r'
then
echo shar: over-writing existing file "'inddup.r'"
fi
cat << \SHAR_EOF > 'inddup.r'
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
SHAR_EOF
if test -f 'initad.r'
then
echo shar: over-writing existing file "'initad.r'"
fi
cat << \SHAR_EOF > 'initad.r'
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
SHAR_EOF
if test -f 'inside.r'
then
echo shar: over-writing existing file "'inside.r'"
fi
cat << \SHAR_EOF > 'inside.r'
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(axmax) {
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
SHAR_EOF
if test -f 'insrt.r'
then
echo shar: over-writing existing file "'insrt.r'"
fi
cat << \SHAR_EOF > 'insrt.r'
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
SHAR_EOF
if test -f 'insrt1.r'
then
echo shar: over-writing existing file "'insrt1.r'"
fi
cat << \SHAR_EOF > 'insrt1.r'
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
SHAR_EOF
if test -f 'locn.r'
then
echo shar: over-writing existing file "'locn.r'"
fi
cat << \SHAR_EOF > 'locn.r'
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
SHAR_EOF
if test -f 'master.r'
then
echo shar: over-writing existing file "'master.r'"
fi
cat << \SHAR_EOF > 'master.r'
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
SHAR_EOF
if test -f 'mnnd.r'
then
echo shar: over-writing existing file "'mnnd.r'"
fi
cat << \SHAR_EOF > 'mnnd.r'
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
SHAR_EOF
if test -f 'pred.r'
then
echo shar: over-writing existing file "'pred.r'"
fi
cat << \SHAR_EOF > 'pred.r'
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
SHAR_EOF
if test -f 'qtest.r'
then
echo shar: over-writing existing file "'qtest.r'"
fi
cat << \SHAR_EOF > 'qtest.r'
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
SHAR_EOF
if test -f 'qtest1.r'
then
echo shar: over-writing existing file "'qtest1.r'"
fi
cat << \SHAR_EOF > 'qtest1.r'
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 'stoke.r'
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=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
SHAR_EOF
if test -f 'succ.r'
then
echo shar: over-writing existing file "'succ.r'"
fi
cat << \SHAR_EOF > 'succ.r'
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
SHAR_EOF
if test -f 'swap.r'
then
echo shar: over-writing existing file "'swap.r'"
fi
cat << \SHAR_EOF > 'swap.r'
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
SHAR_EOF
if test -f 'testeq.r'
then
echo shar: over-writing existing file "'testeq.r'"
fi
cat << \SHAR_EOF > 'testeq.r'
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
SHAR_EOF
if test -f 'triar.r'
then
echo shar: over-writing existing file "'triar.r'"
fi
cat << \SHAR_EOF > 'triar.r'
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
SHAR_EOF
if test -f 'trifnd.r'
then
echo shar: over-writing existing file "'trifnd.r'"
fi
cat << \SHAR_EOF > 'trifnd.r'
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
SHAR_EOF
if test -f 'err.list'
then
echo shar: over-writing existing file "'err.list'"
fi
cat << \SHAR_EOF > 'err.list'
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.
SHAR_EOF
if test -f 'ex.out'
then
echo shar: over-writing existing file "'ex.out'"
fi
cat << \SHAR_EOF > 'ex.out'
$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 F F
2 4.560000 0.740000 4.512766 0.000000 3 1 F T
3 5.750000 5.500000 4.560000 0.740000 3 2 F F
4 0.000000 2.855556 1.650000 3.650000 4 1 T F
5 1.650000 3.650000 3.500000 5.500000 4 2 F F
6 3.500000 5.500000 5.750000 5.500000 5 2 F F
7 5.750000 5.500000 6.058824 5.705882 5 3 F F
8 0.500000 7.500000 3.500000 5.500000 5 4 F F
9 6.058824 5.705882 10.000000 5.142857 6 3 F T
10 5.200000 10.000000 6.058824 5.705882 6 5 T F
11 2.300000 0.000000 0.000000 2.300000 7 1 T T
12 10.000000 3.250000 7.833333 0.000000 8 3 T T
13 0.000000 7.400000 0.500000 7.500000 9 4 T F
14 0.500000 7.500000 2.166667 10.000000 9 5 F T
15 8.750000 10.000000 10.000000 7.500000 10 6 T T
$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
SHAR_EOF
# End of shell archive
exit 0