#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh 'README' <<'END_OF_FILE'
X9/9/94
XREADME file for delaunay.S
XS function and fortran code for delaunay triangulation
Xand dirichlet tesselation. Based on Rolf Turner's
Xdelaunay package in statlib/general. I recommend
Xyou retriev Rolf's package also, for documentation
Xand test cases.
X
X=========================================================
XDon MacQueen
XLawrence Livermore National Lab
Xmacq@llnl.gov
X=========================================================
X#
X# Permission to use, copy, modify, and distribute this software and
X# its documentation for any purpose and without fee is hereby
X# granted, provided that this permission
X# notice appears in all copies and in supporting documentation.
X
Xdistribution consists of four files:
X README
X deldir.r
X deldir.sp
X segments.sp
X
XREADME is this file
X
Xdeldir.r is ratfor (fortran) source code
X contains subroutine deldir, and many others
X
Xdeldir.sp contains an S function deldir that
X calls subroutine deldir
X
Xsegments.sp contains an S function called segments.
X It is an alternate version of Splus version.
X It is optional.
X
XTo install:
X
XMake sure you have an environment variable SHOME defined,
Xand that it points to the top level Splus directory. This is
Xnecessary both for installation and use.
X
XCreate a directory to hold these files, or use an existing
Xdirectory. I use a directory called "delaunay" in the Splus
Xlibrary directory:
X $SHOME/library/delaunay
XIt also needs a .Data subdirectory. In my case:
X $SHOME/library/delaunay/.Data
X
XIf you put them somewhere different than I did, you'll have
Xto edit the function deldir() in the file
Xdeldir.sp so that the argument to dyn.load() points to
Xthe right location.
X
Xcd to wherever you put the files, and execute:
X Splus COMPILE deldir.r
X Splus < deldir.sp
X
XIf you put them in a library as I did, you can use
Xthe function:
X library(delaunay)
Xto make the available. Otherwise, you'll have to attach
Xthe directory you put them in, using attach().
X
XIf you want to use my version of segments(), source it
Xso that it is located earlier in the search path
Xthan the Splus built in version of segments().
X
XI have tested these functions on a Sun Sparcstation LX,
Xrunning SunOS 4.1.3, and Splus 3.2. They obtain the same result as
Xthe sample problems in Rolf Turner's delaunay
Xsubmission to statlib/general.
X
XSorry, I didn't make a help file.
END_OF_FILE
if test 2361 -ne `wc -c <'README'`; then
echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'README~' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'README~'\"
else
echo shar: Extracting \"'README~'\" \(2115 characters\)
sed "s/^X//" >'README~' <<'END_OF_FILE'
X9/9/94
XREADME file for delaunayS
XS function and fortran code for delaunay triangulation
Xand dirichlet tesselation. Based on Rolf Turner's
Xdelaunay package in statlib/general. I recommend
Xyou retriev Rolf's package also, for documentation
Xand test cases.
X
X=========================================================
XDon MacQueen
XLawrence Livermore National Lab
Xmacq@llnl.gov
X=========================================================
X
Xdistribution consists of four files:
X README
X deldir.r
X deldir.sp
X segments.sp
X
XREADME is this file
X
Xdeldir.r is ratfor (fortran) source code
X contains subroutine deldir, and many others
X
Xdeldir.sp contains an S function deldir that
X calls subroutine deldir
X
Xsegments.sp contains an S function called segments.
X It is an alternate version of Splus version.
X It is optional.
X
XTo install:
X
XMake sure you have an environment variable SHOME defined,
Xand that it points to the top level Splus directory. This is
Xnecessary both for installation and use.
X
XCreate a directory to hold these files, or use an existing
Xdirectory. I use a directory called "delaunay" in the Splus
Xlibrary directory:
X $SHOME/library/delaunay
XIt also needs a .Data subdirectory. In my case:
X $SHOME/library/delaunay/.Data
X
XIf you put them somewhere different than I did, you'll have
Xto edit the function deldir() in the file
Xdeldir.sp so that the argument to dyn.load() points to
Xthe right location.
X
Xcd to wherever you put the files, and execute:
X Splus COMPILE deldir.r
X Splus < deldir.sp
X
XIf you put them in a library as I did, you can use
Xthe function:
X library(delaunay)
Xto make the available. Otherwise, you'll have to attach
Xthe directory you put them in, using attach().
X
XIf you want to use my version of segments(), source it
Xso that it is located earlier in the search path
Xthan the Splus built in version of segments().
X
XI have tested these functions on a Sun Sparcstation LX,
Xrunning SunOS 4.1.3, and Splus 3.2. They obtain the same result as
Xthe sample problems in Rolf Turner's delaunay
Xsubmission to statlib/general.
X
XSorry, I didn't make a help file.
END_OF_FILE
if test 2115 -ne `wc -c <'README~'`; then
echo shar: \"'README~'\" unpacked with wrong size!
fi
# end of 'README~'
fi
if test -f 'deldir.r' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'deldir.r'\"
else
echo shar: Extracting \"'deldir.r'\" \(47239 characters\)
sed "s/^X//" >'deldir.r' <<'END_OF_FILE'
X # deldir
X #
X # Copyright (C) 1991 by T. Rolf Turner
X #
X # Permission to use, copy, modify, and distribute this software and
X # its documentation for any purpose and without fee is hereby
X # granted, provided that the above copyright notice appear in all
X # copies and that both that copyright notice and this permission
X # notice appear in supporting documentation.
X #
X # PROGRAMMED BY: Rolf Turner, while with the Division of Mathematics
X # and Statistics, CSIRO, Sydney, Australia.
X #
X # Program to compute the Delaunay Triangulation (and hence the
X # Dirichlet Tesselation) of a planar point set according to the
X # second (iterative) algorithm of Lee and Schacter, International
X # Journal of Computer and Information Sciences, Vol. 9, No. 3, 1980,
X # pages 219 to 242.
X
X # The triangulation is made to be with respect to the whole plane by
X # `suspending' it from `ideal' points
X # (-R,-R), (R,-R) (R,R), and (-R,R), where R --> infinity.
X
X # It is also enclosed in a finite rectangle (whose boundaries truncate any
X # infinite Dirichlet tiles) with corners (xmin,ymin) etc. This rectangle
X # is referred to elsewhere as `the' rectangular window.
X
X # 9/9/94
X # Subroutine version for use with SPlus and dyn.load()
X # Modifications by Donald H. MacQueen
X # Lawrence Livermore National Laboratory, Livermore, CA, USA
X
X # 7/5/94
X # Problems with undefined symbols using dynamic load
X # Eliminate all print statements, replace with cerr='whatever was printed'
X # and set error code ierr=1. Then replace stop statements with
X # return statements. Add check for ierr=1 at beginning of every subroutine
X # and after every subroutine call, to make sure return is propagated
X # back up to the main subroutine deldir, for return to the S calling
X # function
X
X # NOTE: as of 9/9/94, error handling is disabled. I hope to fix this
X # sometime in the future
X
X # Combined all source code into one file, to avoid problems with the order
X # in which object files are loaded. See the warning in section 9.5.3 of S-Plus
X # version 3.2 Programmer's Manual, page 9-18.
X
X # Most modifications have to do with input and output. For example, instead
X # of writing to a file, saving to a matrix. See especially subroutines
X # delseg, delout, dirseg, and dirout.
X
X # 7/8/94
X # add subroutine to initialize all variables in common to zero
X
X subroutine deldir(xin,yin,ndxin,ndyin,xminin,xmaxin,yminin,ymaxin,narg,nptsin,triang,ntseg,tarea,tareatot,polyg,npseg,parea,pareatot,ierro)
X
X include defns.h
X
X include commons.h
X
X dimension triang( 4*(NPTS-1)*NPTS/2) , tarea(NPTS)
X dimension polyg( 4*(NPTS-1)*NPTS/2) , parea(NPTS)
X dimension xin(NPTS),yin(NPTS)
X character*80 errmsg
X
X # initialize variables in common
X call initcom
X # ierr=0
X # cerr=''
X
X # Fill variables in common with variables in subroutine argument
X ndx=ndxin
X ndy=ndyin
X xmin=xminin
X xmax=xmaxin
X ymin=yminin
X ymax=ymaxin
X npts=nptsin
X do i = 1,npts {
X v(i,1)=xin(i)
X v(i,2)=yin(i)
X }
X
X # initialize (getting bizarre behavior when called more than once
X # from within SPlus)
X tareatot=0.
X pareatot=0.
X do i=1,npts {
X tarea(i)=0.
X parea(i)=0.
X }
X do i=1,4*(npts-1)*npts/2 {
X triang(i)=0.
X polyg(i)=0.
X }
X
X # If corners of the window are not specified, form them from
X # the minimum and maximum of the data +/- 10%.
X if(narg<6) {
X call datwin ; if (ierr) { ierro=ierr; errmsg=cerr; return }
X }
X
X # Add in dummy points if required.
X if(narg>=2) call dumpts(ndx,ndy) ; if (ierr) { ierro=ierr; errmsg=cerr; return }
X
X # Sort the points into bins, the number of such being approx. sqrt(npts)
X call binsrt ; if (ierr) { ierro=ierr; errmsg=cerr; return }
X
X # Put the four ideal points into the adjacency list.
X do i = 1,4 {
X j = i-4
X k = j+1
X if(k>0) k = -3
X call insrt(j,k) ; if (ierr) { ierro=ierr; errmsg=cerr; return }
X }
X
X # Put in the first of the point set into the adjacency list.
X do i = 1,4 {
X j = i-4
X call insrt(1,j) ; if (ierr) { ierro=ierr; errmsg=cerr; return }
X }
X
X # Now add the rest of the point set
X do j = 2,npts {
X call addpt(j) ; if (ierr) { ierro=ierr; errmsg=cerr; return }
X }
X
X # Output the triangulation.
X # Order in which called is critical, see dirout comments
X call delseg(triang,ntseg)
X call delout(tarea,tareatot)
X call dirseg(polyg,npseg)
X call dirout(parea,pareatot)
X
X errmsg=cerr
X ierro=ierr
X
X # Aw' done!!!
X return
X end
X subroutine initcom
X # initialize all variables in common to zero
X include commons.h
X do i=-3,NPTS {
X do j=0,NADJ {
X nadj(i,j)=0
X }
X }
X do i=1,NPTS {
X ind(i)=0
X indic(i)=0
X v(i,1)=0.
X v(i,2)=0.
X }
X npts=0
X xmin=0. ; xmax=0. ; ymin=0.; ymax=0.
X ierr=0 ; cerr=''
X return
X end
X # include defns.h
X subroutine acchk(i,j,k,anticl)
X # Check whether vertices i, j, k, are in anti-clockwise order.
X
X logical anticl
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # Create indicator telling which of i, j, and k are ideal points.
X if(i<=0) i1 = 1
X else i1 = 0
X if(j<=0) j1 = 1
X else j1 = 0
X if(k<=0) k1 = 1
X else k1 = 0
X ijk = i1*4+j1*2+k1
X
X # Get the coordinates of vertices i, j, and k. (Pseudo-coordinates for
X # any ideal points.)
X call coords(i,xi,yi) ; # if (ierr) return
X call coords(j,xj,yj) ; # if (ierr) return
X call coords(k,xk,yk) ; # if (ierr) return
X
X # Form the vectors from i to j, and from i to k; how this vector is
X # formed depends upon which points are ideal.
X switch(ijk) {
X case 0: # No ideal points.
X a = xj-xi
X b = yj-yi
X c = xk-xi
X d = yk-yi
X case 1: # Only k ideal.
X a = xj-xi
X b = yj-yi
X c = xk
X d = yk
X case 2: # Only j ideal.
X a = xj
X b = yj
X c = xk-xi
X d = yk-yi
X case 3: # Both j and k ideal (i not).
X a = xj
X b = yj
X c = xk
X d = yk
X case 4: # Only i ideal.
X a = -xi
X b = -yi
X p = xj-xk
X q = yj-yk
X c = a-p
X d = b-q
X case 5: # Both i and k ideal (j not).
X a = -xi
X b = -yi
X c = xk-xi
X d = yk-yi
X case 6: # Both i and j ideal (k not).
X a = xj-xi
X b = yj-yi
X c = -xi
X d = -yi
X case 7: # All three points ideal.
X a = xj-xi
X b = yj-yi
X c = xk-xi
X d = yk-yi
X }
X
X # If (ij x ik) is directed ***upwards*** then i, j, k,
X # are in anti-clockwise order; else not.
X cross = a*d-b*c
X if(cross>0.) anticl = .true.
X else anticl = .false.
X return
X end
X # include defns.h
X subroutine addpt(j)
X # Add point j to the triangulation.
X
X logical didswp
X integer succ
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # Put the new point in, joined to the vertices of its
X # enclosing triangle.
X call initad(j) ; # if (ierr) return
X
X # Look at each `gap', i.e. pair of adjacent segments
X # emanating from the new point; they form two sides of a
X # quadrilateral; see whether the extant diagonal of this
X # quadrilateral should be swapped with its alternate
X # (according to the LOP: local optimality principle).
X now = nadj(j,1)
X nxt = nadj(j,2)
X ngap = 0
X repeat {
X call swap(j,now,nxt,didswp) ; # if (ierr) return
X n = nadj(j,0)
X if(!didswp) { # If no swap of diagonals
X now = nxt # move to the next gap.
X ngap = ngap+1
X }
X nxt = succ(j,now)
X }
X until(ngap==n)
X
X return
X end
X # include defns.h
X subroutine adjchk(i,j,adj)
X # Check if vertices i and j are adjacent.
X
X logical adj
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # Check if j is in the adjacency list of i.
X adj = .false.
X ni = nadj(i,0)
X if(ni>0) {
X do k = 1,ni {
X if(j==nadj(i,k)) {
X adj = .true.
X break
X }
X }
X }
X
X # Check if i is in the adjacency list of j.
X nj = nadj(j,0)
X if(nj>0) {
X do k = 1,nj {
X if(i==nadj(j,k)) {
X if(adj) return # Have j in i's list and i in j's.
X else {
X cerr = 'Contradictory adjacency lists.'
X ierr = 1 ; return
X }
X }
X }
X }
X
X # If we get to here i is not in j's list.
X if(adj) { # If adj is true, then j IS in i's list.
X cerr = 'Contradictory adjacency lists.'
X ierr = 1 ; return
X }
X
X return
X end
X # include defns.h
X subroutine binsrt
X # Sort the data points into bins.
X
X dimension tv(NPTS,2), ilst(NPTS)
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # initialize
X do i=1,npts {tv(i,1)=0. ; tv(i,2)=0.}
X
X w = xmax-xmin
X h = ymax-ymin
X
X # Number of bins is to be approx. sqrt(npts); thus number of subdivisions
X # on each side of rectangle is approx. npts**(1/4).
X kdiv = 1+float(npts)**0.25 # Round high.
X dw = w/float(kdiv)
X dh = h/float(kdiv)
X
X # The width of each bin is dw; the height is dh. We shall move across
X # the rectangle from left to right, then up, then back from right to
X # left, then up, .... Note that kx counts the divisions from the left,
X # ky counts the divisions from the bottom; kx is incremented by ink, which
X # is +/- 1 and switches sign when ky is incremented; ky is always
X # incremented by 1.
X kx = 1
X ky = 1
X ink = 1
X k = 0
X do i = 1,npts { ilst(i) = 0 } # Keeps a list of those points already added
X while(ky<=kdiv) { # to the new list.
X do i = 1,npts {
X if(ilst(i)==1) next # The i-th point has already been added
X # to the new list.
X # If the i-th point is in the current bin, add it to the list.
X x = v(i,1)
X y = v(i,2)
X ix = 1+(x-xmin)/dw
X if(ix>kdiv) ix = kdiv
X jy = 1+(y-ymin)/dh
X if(jy>kdiv) jy = kdiv
X if(ix==kx&jy==ky) {
X k = k+1
X ind(i) = k # Index i is the pos'n. of (x,y) in the
X tv(k,1) = x # old list; k is its pos'n. in the new one.
X tv(k,2) = y
X ilst(i) = 1 # Cross the i-th point off the old list.
X }
X }
X # Move to the next bin.
X kc = kx+ink
X if((1<=kc)&(kc<=kdiv)) kx = kc
X else {
X ky = ky+1
X ink = -ink
X }
X }
X
X # Check that all points from old list have been added to the new,
X # with no spurious additions.
X if(k!=npts) {
X cerr = 'Number of points jumbled in binsrt.'
X ierr = 1 ; return
X }
X
X # Copy the new list back into the old one.
X do i = 1,npts {
X v(i,1) = tv(i,1)
X v(i,2) = tv(i,2)
X }
X return
X end
X # include defns.h
X subroutine circen(i,j,k,x0,y0,collin)
X # Find the circumcentre (x0,y0) of the triangle with
X # vertices (xi,yi), (xj,yj), (xk,yk).
X
X include commons.h
X logical collin
X
X # Check for error code
X # if (ierr) return
X
X # Get the coordinates.
X xi = v(i,1)
X yi = v(i,2)
X xj = v(j,1)
X yj = v(j,2)
X xk = v(k,1)
X yk = v(k,2)
X
X # Check for collinearity of the three vertices.
X fac = (xk-xi)*(yk-yj)-(yk-yi)*(xk-xj)
X call testeq(fac,0.,collin) ; # if (ierr) return
X
X # If they're collinear, make sure that they're in the right
X # order (1, 2, 3).
X if(collin) {
X a = xi-xj
X b = yi-yj
X c = xk-xj
X d = yk-yj
X alpha = (a*c+b*d)/(c*c+d*d)
X # If they're not in the right order, bring things to
X # a shuddering halt.
X if(alpha<=0.|alpha>=1.) {
X cerr = 'Vertices of ''triangle'' are collinear'
X cerr = 'and vertex 2 is not between 1 and 3.'
X cerr = 'Something is wrong'
X ierr = 1 ; return
X }
X # Collinear, but in the right order; think of this as meaning
X # that the circumcircle in question has infinite radius.
X return
X }
X
X # Not collinear; go ahead, make my circumcentre.
X fac = 0.5/fac
X
X a = yk-yj
X b = -yk+yi
X c = -xk+xj
X d = xk-xi
X
X e1 = xk*xk-xi*xi+yk*yk-yi*yi
X e2 = xk*xk-xj*xj+yk*yk-yj*yj
X
X x0 = (a*e1+b*e2)*fac
X y0 = (c*e1+d*e2)*fac
X
X return
X end
X # include defns.h
X subroutine coords(i,a,b)
X # Get the coordinates (a,b) of point i, substituting pseudo-coordinates
X # if i is ideal.
X
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # The ideal points are given pseudo-coordinates (-1,-1), (1,-1)
X # (1,1), and (-1,1). They are numbered anticlockwise from the
X # `bottom left corner' as 0, -1, -2, -3.
X if(i<=0) {
X ia = -i
X a = (-1)**(1+(ia+1)/2)
X b = (-1)**(1+ia/2)
X return
X }
X
X # The point i is real (not ideal); just obtain its coordinates and
X # go home.
X a = v(i,1)
X b = v(i,2)
X
X return
X end
X # include defns.h
X subroutine datwin
X # Rectangular window for the data points has not been specified;
X # form it from the respective maxima and minima of the data +/- 10%.
X
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # Find the maxima and minima in the x and y directions.
X xmin = v(1,1)
X ymin = v(1,2)
X xmax = v(1,1)
X ymax = v(1,2)
X
X do i = 2,npts {
X if(v(i,1) < xmin) xmin = v(i,1)
X if(v(i,2) < ymin) ymin = v(i,2)
X if(v(i,1) > xmax) xmax = v(i,1)
X if(v(i,2) > ymax) ymax = v(i,2)
X }
X
X # Differences to form basis for the 10% expansion.
X xdiff = xmax-xmin
X ydiff = ymax-ymin
X
X # Expand.
X xmin = xmin-0.1*xdiff
X xmax = xmax+0.1*xdiff
X ymin = ymin-0.1*ydiff
X ymax = ymax+0.1*ydiff
X
X return
X end
X # include defns.h
X subroutine delet(i,j)
X # Delete i and j from each other's adjacency lists.
X
X logical adj
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # First check that they're IN each other's lists.
X call adjchk(i,j,adj) ; # if (ierr) return
X
X # Then do the actual deletion if they are.
X if(adj) {
X call delet1(i,j) ; # if (ierr) return
X call delet1(j,i) ; # if (ierr) return
X }
X
X return
X end
X # include defns.h
X subroutine delet1(i,j)
X # Delete j from the adjacency list of i.
X
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X n = nadj(i,0)
X do k = 1,n {
X if(nadj(i,k)==j) { # Find j in the list;
X # then move everything back one notch.
X do kk = k,n-1 { nadj(i,kk) = nadj(i,kk+1) }
X nadj(i,n) = 0
X nadj(i,0) = n-1
X return
X }
X }
X
X end
X # include defns.h
X subroutine delout(tarea,ars)
X # Output the description of the Delaunay triangles with a vertex
X # at point i, for i = 1, ..., npts. Do this in the original order
X # of the points, not the order into which they have been bin-sorted.
X
X # 7/1/94
X # modified to return vector of areas
X
X integer succ
X
X include commons.h
X
X dimension tarea(NPTS)
X
X # Check for error code
X # if (ierr) return
X
X ars = 0. # Initialize area sum to zero.
X do i1 = 1,npts {
X area = 0. # Initialize area of polygon consisting of triangles
X i = ind(i1) # with a vertex at point i.
X np = nadj(i,0)
X call coords(i,xi,yi) ; # if (ierr) return
X
X # Output the point number, its coordinates and the number of
X # (real) triangles emanating from it.
X npt = np
X do k = 1,np {
X kp = k+1
X if(kp>np) kp = 1
X if(nadj(i,k)<=0|nadj(i,kp)<=0) npt = npt-1
X }
X
X # For each point in the adjacency list of point i, find its
X # successor, and the area of the triangle determined by these
X # three points.
X do j1 = 1,np {
X j = nadj(i,j1)
X if(j<=0) next
X call coords(j,xj,yj) ; # if (ierr) return
X k = succ(i,j)
X if(k<=0) next
X call coords(k,xk,yk) ; # if (ierr) return
X call triar(xi,yi,xj,yj,xk,yk,tmp) ; # if (ierr) return
X # Downweight the area by 1/3, since each
X # triangle eventually appears 3 times over.
X area = area+tmp/3.
X }
X tarea(i1)=area
X ars = ars+area
X }
X # Since the area of each triangle has been downweighted by 1/3
X # the sum `ars' is indeed the total area of the convex hull of
X # the point set being triangulated.
X
X return
X end
X # include defns.h
X subroutine delseg(triang,ntseg)
X # Output the endpoints of the line segments joining the
X # vertices of the Delaunay triangles.
X
X # 7/1/94
X # Modified to return matrix
X
X logical value
X
X include commons.h
X
X dimension triang( 4*(NPTS-1)*NPTS/2)
X
X # Check for error code
X # if (ierr) return
X
X n = npts
X
X # For each distinct pair of points i and j, if they are adjacent
X # then output their endpoints.
X ij=0
X do i = 2,n {
X do j = 1,i-1 {
X call adjchk(i,j,value) ; # if (ierr) return
X if(value) {
X ij=ij+4
X triang(ij-3) = v(i,1)
X triang(ij-2) = v(i,2)
X triang(ij-1) = v(j,1)
X triang(ij) = v(j,2)
X }
X }
X }
X ntseg=ij
X return
X end
X# include defns.h
Xsubroutine dirout(parea,ars)
X# Output the description of the Dirichlet tile centred at point
X# i for i = 1, ..., npts. Do this in the original order of the
X# points, not in the order into which they have been bin-sorted.
X
X# 7/1/94 DHM
X# Return vector for use in S function
X
Xinteger pred, succ
Xlogical collin, intfnd, bptab, bptcd
X
Xinclude commons.h
X dimension parea(NPTS)
X
X # Check for error code
X # if (ierr) return
X
X# Note that at this point some Delaunay neighbors may be
X# `spurious'; they are the corners of a `large' rectangle in which
X# the rectangular window of interest has been suspended. This
X# large window was brought in simply to facilitate output concerning
X# the Dirichlet tesselation. They were added to the triangulation
X# in the routine `dirseg' which ***must*** therefore be called before
X# this routine (`dirout') is called. (Likewise `dirseg' must be called
X# ***after*** `delseg' and `delout' have been called.)
X
Xars = 0. # Initialize the sum of the tile areas to zero.
Xdo i1 = 1,npts {
X area = 0. # Initialize the area of the ith tile to zero.
X i = ind(i1)
X np = nadj(i,0)
X npt = np
X
X # Adjust number of tile boundaries for `fake outer corners'.
X do k = 1,np {if(nadj(i,k)>npts) npt = npt-1}
X xi = v(i,1)
X yi = v(i,2)
X
X # Output the point number, its coordinates, and the number of
X # its Delaunay neighbors == the number of boundary segments in
X # its Dirichlet tile.
X
X # For each Delaunay neighbor, find the circumcentres of the
X # triangles on each side of the segment joining point i to that
X # neighbor.
X do j1 = 1,np {
X j = nadj(i,j1)
X xj = v(j,1)
X yj = v(j,2)
X xij = 0.5*(xi+xj)
X yij = 0.5*(yi+yj)
X k = pred(i,j)
X l = succ(i,j)
X call circen(i,k,j,a,b,collin) ; # if (ierr) return
X if(collin) {
X cerr = 'Vertices of ''triangle'' are collinear.'
X ierr = 1 ; return
X }
X call circen(i,j,l,c,d,collin) ; # if (ierr) return
X if(collin) {
X cerr = 'Vertices of ''triangle'' are collinear.'
X ierr = 1 ; return
X }
X
X # Increment the area of the current Dirichlet
X # tile (intersected with the rectangular window) by applying
X # Stoke's Theorem to the segment of tile boundary joining
X # (a,b) to (c,d). (Note that the direction is anti-clockwise.)
X call stoke(a,b,c,d,tmp,sn) ; # if (ierr) return
X area = area+sn*tmp
X
X # If a circumcentre is outside the rectangular window, replace
X # it with the intersection of the rectangle boundary with the
X # line joining the circumcentre to the midpoint of
X # (xi,yi)->(xj,yj). Then output the number of the current
X # Delaunay neighbor and the two circumcentres (or the points
X # with which they have been replaced).
X call inside(a,b,xij,yij,ai,bi,intfnd,bptab) ; # if (ierr) return
X if(intfnd) {
X call inside(c,d,xij,yij,ci,di,intfnd,bptcd) ; # if (ierr) return
X }
X }
X parea(i1) = area
X ars = ars+area # Increment the sum of areas.
X}
X
X50 format(i5,2f10.5,i5)
X51 format(i5,4f10.5)
X511 format(i5,2f10.5,' (bndry pt) ',2f10.5)
X512 format(i5,4f10.5,' (bndry pt) ')
X52 format(f15.5)
X53 format(/f15.5)
X
Xreturn
Xend
X # include defns.h
X subroutine dirseg(polyg,npseg)
X # Output the endpoints of the segments of boundary of Dirichlet
X # tiles. (Do it economically; each such segment once and only once.)
X
X logical collin, adjace, dum1, dum2
X integer pred, succ
X
X include commons.h
X
X dimension polyg( 4*(NPTS-1)*NPTS/2)
X
X # Check for error code
X # if (ierr) return
X
X # Add in some dummy corner points, outside the actual window.
X # Far enough out so that no resulting tile boundaries intersect the
X # window.
X
X # Note that these dummy corners are needed by the routine `dirout'
X # but will screw things up for `delseg' and `delout'. Therefore
X # this routine (`dirseg') must be called ***before*** dirout, and
X # ***after*** delseg and delout.
X
X a = xmax-xmin
X b = ymax-ymin
X c = sqrt(a*a+b*b)
X
X i = npts+1
X v(i,1) = xmin-c
X v(i,2) = ymin-c
X i = i+1
X v(i,1) = xmax+c
X v(i,2) = ymin-c
X i = i+1
X v(i,1) = xmax+c
X v(i,2) = ymax+c
X i = i+1
X v(i,1) = xmin-c
X v(i,2) = ymax+c
X
X nstt = npts+1 # Start and finish of indices for adding
X nfin = npts+4 # the four dummy corner points.
X do j = nstt,nfin {
X call addpt(j) ; # if (ierr) return
X }
X
X # For each distinct pair of (genuine) data points, find out if they are
X # adjacent. If so, find the circumcentres of the triangles lying on each
X # side of the segment joining them.
X ij = 0
X do i = 2,npts {
X do j = 1,i-1 {
X call adjchk(i,j,adjace) ; # if (ierr) return
X if(adjace) {
X xi = v(i,1)
X yi = v(i,2)
X xj = v(j,1)
X yj = v(j,2)
X # Let (xij,yij) be the midpoint of the segment joining
X # (xi,yi) to (xj,yj).
X xij = 0.5*(xi+xj)
X yij = 0.5*(yi+yj)
X k = pred(i,j)
X call circen(i,k,j,a,b,collin) ; # if (ierr) return
X if(collin) {
X cerr = 'Vertices of ''triangle'' are'
X cerr = 'collinear.'
X ierr = 1 ; return
X }
X
X # If a circumcentre is outside the rectangular window
X # of interest, draw a line joining it to the mid-point
X # of (xi,yi)->(xj,yj). Find the intersection of this
X # line with the boundary of the window; for (a,b),
X # call the point of intersection (ai,bi). For (c,d),
X # call it (ci,di).
X call inside(a,b,xij,yij,ai,bi,dum1,dum2) ; # if (ierr) return
X l = succ(i,j)
X call circen(i,j,l,c,d,collin) ; # if (ierr) return
X if(collin) {
X cerr = 'Vertices of ''triangle'' are'
X cerr = 'collinear.'
X ierr = 1 ; return
X }
X call inside(c,d,xij,yij,ci,di,dum1,dum2) ; # if (ierr) return
X ij=ij+4
X polyg(ij-3) = ai
X polyg(ij-2) = bi
X polyg(ij-1) = ci
X polyg(ij) = di
X }
X }
X }
X npseg=ij
X
X return
X end
X # include defns.h
X subroutine dumpts(ndx,ndy)
X # Add in an ndx-by-ndy grid of dummy points.
X
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X if(ndx<=0|ndy<=0) return
X
X # If only one grid-line perp. to the x-axis, made it the midline;
X # else make the lines equi-spaced from side to side.
X if(ndx==1) x0 = 0.5*(xmin+xmax)
X else {
X dx = (xmax-xmin)/(ndx-1)
X x0 = xmin
X }
X
X # Likewise for the grid-lines perp. to the y-axis.
X if(ndy==1) y0 = 0.5*(ymin+ymax)
X else {
X dy = (ymax-ymin)/(ndy-1)
X y0 = ymin
X }
X x = x0
X y = y0
X n = npts
X
X # Run along the grid-lines vertically; there are ndx such
X # vertical lines; add ndy points on each.
X do i = 1,ndx {
X do j = 1,ndy {
X n = n+1
X v(n,1) = x
X v(n,2) = y
X indic(n) = 0
X y = y+dy
X }
X x = x+dx
X y = y0
X }
X
X npts = n
X return
X end
X subroutine initad(j)
X # Intial adding-in of a new point j.
X
X integer tau(3), pred, succ
X
X # Check for error code
X # if (ierr) return
X
X # Find the triangle containing vertex j.
X call trifnd(j,tau,nedge) ; # if (ierr) return
X
X # Check for error code
X # if (ierr) return
X
X # If the new point is on the edge of a triangle, detach the two
X # vertices of that edge from each other. Also join j to the vertex
X # of the triangle on the reverse side of that edge from the `found'
X # triangle (defined by tau) -- given that there ***is*** such a triangle.
X if(nedge!=0) {
X ip = nedge
X i = ip-1
X if(i==0) i = 3 # Arithmetic modulo 3.
X k = pred(tau(i),tau(ip))
X kk = succ(tau(ip),tau(i))
X call delet(tau(i),tau(ip)) ; # if (ierr) return
X if(k==kk) call insrt(j,k) ; # if (ierr) return
X }
X
X # Join the new point to each of the three vertices.
X do i = 1,3 {
X call insrt(j,tau(i)) ; # if (ierr) return
X }
X
X return
X end
X# include defns.h
Xsubroutine inside(a,b,c,d,ai,bi,intfnd,bpt)
X# Get a point ***inside*** the rectangular window on the ray from
X# one circumcentre to the next one. I.e. if the `next one' is
X# inside, then that's it; else find the intersection of this ray with
X# the boundary of the rectangle.
X
Xinclude commons.h
Xlogical intfnd, bpt
X
X # Check for error code
X # if (ierr) return
X
X# Note that (a,b) is the circumcenter of a Delaunay triangle,
X# and that (c,d) is the midpoint of one of its sides.
X# When `inside' is called by `dirout' it is possible for (c,d) to
X# lie ***outside*** the rectangular window, and for the ray not to
X# intersect the window at all. (The point (c,d) might be the midpoint
X# of a Delaunay edge connected to a `fake outer corner', added to facilitate
X# constructing a tiling that completely covers the actual window.)
X# The variable `indfnd' acts as an indicator as to whether such an
X# intersection has been found.
X
X# The variable `bpt' acts as an indicator as to whether the returned
X# point (ai,bi) is a true circumcentre, inside the window (bpt == .false.),
X# or is the intersection of a ray with the boundary of the window
X# (bpt = .true.).
X
Xintfnd = .true.
Xbpt = .true.
X
X# Check if (a,b) is inside the rectangle.
Xif(xmin<=a&a<=xmax&ymin<=b&b<=ymax) {
X ai = a
X bi = b
X bpt = .false.
X return
X}
X
X# Look for appropriate intersections with the four lines forming
X# the sides of the rectangular window.
X
X# Line 1: x = xmin.
Xif(axmax) {
X ai = xmax
X s = (d-b)/(c-a) # See above.
X t = b-s*a
X bi = s*ai+t
X if(yminymax) {
X bi = ymax
X s = (c-a)/(d-b) # See above.
X t = a-s*b
X ai = s*bi+t
X if(xminNADJ) { # Watch out for over-writing!!!
X cerr = 'Number of adjacencies too large.'
X cerr = 'Increase NADJ in defns.h and recompile.'
X ierr = 1 ; return
X }
X while(kk>kj) {
X nadj(i,kk) = nadj(i,kk-1)
X kk = kk-1
X }
X nadj(i,kj) = j
X nadj(i,0) = n+1
X
X return
X end
X # include defns.h
X subroutine locn(i,j,kj)
X # Find the appropriate location for j in the adjacency list
X # of i. This is the index which j ***will*** have when
X # it is inserted into the adjacency list of i in the
X # appropriate place.
X
X logical before
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X n = nadj(i,0)
X
X # If there is nothing already adjacent to i, then j will have place 1.
X if(n==0) {
X kj = 1
X return
X }
X
X # Run through i's list, checking if j should come before each element
X # of that list. (I.e. if i, j, and k are in anti-clockwise order.)
X # If j comes before the kj-th item, but not before the (kj-1)-st, then
X # j should have place kj.
X do ks = 1,n {
X kj = ks
X k = nadj(i,kj)
X call acchk(i,j,k,before) ; # if (ierr) return
X if(before) {
X km = kj-1
X if(km==0) km = n
X k = nadj(i,km)
X call acchk(i,j,k,before) ; # if (ierr) return
X if(before) next
X # If j is before 1 and after n, then it should
X # have place n+1.
X if(kj==1) kj = n+1
X return
X }
X }
X
X # We've gone right through the list and haven't been before
X # the kj-th item ***and*** after the (kj-1)-st on any occasion.
X # Therefore j is before everything (==> place 1) or after
X # everything (==> place n+1).
X if(before) kj = 1
X else kj = n+1
X end
X # include defns.h
X integer function pred(i,j)
X # Find the predecessor of j in the adjacency list of i.
X
X integer pred
X include commons.h
X
X n = nadj(i,0)
X
X # If the adjacency list of i is empty, then clearly j has no predecessor
X # in this adjacency list. Something's wrong; stop.
X if(n==0) {
X # cerr = 'Adjacency of ',i,' is empty,'
X # cerr = 'and so cannot contain ',j,'.'
X # cerr = 'Error in ''pred''.'
X cerr = 'Adjacency problem in function pred()'
X ierr = 1 ; return
X }
X
X # The adjacency list of i is non-empty; search through it until j is found;
X # subtract 1 from the location of j, and find the contents of this new location
X do k = 1,n {
X if(j==nadj(i,k)) {
X km = k-1
X if(km<1) km = n # Take km modulo n. (The adjacency list
X pred = nadj(i,km) # is circular.)
X return
X }
X }
X
X # The adjacency list doesn't contain j. Something's wrong; stop.
X # cerr = 'Adjacency list of ',i,' does not contain ',j,'.'
X # cerr = 'Error in ''pred''.'
X cerr = 'Adjacency problem in function pred()'
X ierr = 1 ; return
X end
X # include defns.h
X subroutine qtest(h,i,j,k,shdswp)
X # Test whether the LOP is satisified; i.e. whether
X # vertex j is outside the circumcircle of vertices h, i, and k
X # of the quadrilateral. (So that vertices i and k should be joined,
X # rather than h and j.) If the LOP is not satisfied, the shdswp
X # ("should-swap") is true.
X
X integer h
X logical shdswp
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # Look for ideal points.
X if(i<=0) ii = 1
X else ii = 0
X if(j<=0) jj = 1
X else jj = 0
X if(k<=0) kk = 1
X else kk = 0
X ijk = ii*4+jj*2+kk
X switch(ijk) {
X
X # If all three corners other than h (the point currently being
X # added) are ideal, then LOP is ***not*** satisfied, i.e. swap.
X case 7:
X shdswp = .true.
X return
X
X # If i and j are ideal, find out which of h and k is closer to the
X # intersection point of the two diagonals, and choose the diagonal
X # emanating from that vertex. (I.e. if k is closer, swap.)
X case 6:
X ja = -j
X xh = v(h,1)
X yh = v(h,2)
X xk = v(k,1)
X yk = v(k,2)
X test =(xh*yk+xk*yh-xh*yh-xk*yk)*(-1)**ja
X if(test>0) shdswp = .true.
X else shdswp = .false.
X return
X
X # If i and k are ideal get minimum angle of pi/4 rather than 0
X # by swapping.
X case 5:
X shdswp = .true.
X return
X
X # If i is ideal get minimum angle > 0 by swapping, but we may
X # not be looking at a convex quadrilateral.
X case 4:
X call acchk(j,k,h,shdswp) ; # if (ierr) return # This checks for convexity.
X return
X
X # If j and k are ideal, this is like unto case 6.
X case 3:
X ja = -j
X xi = v(i,1)
X yi = v(i,2)
X xh = v(h,1)
X yh = v(h,2)
X test = (xh*yi+xi*yh-xh*yh-xi*yi)*(-1)**ja
X if(test>0.) shdswp = .true.
X else shdswp = .false.
X return
X
X # If j is ideal we have a minimum angle > 0, and would get a 0 one
X # by swapping.
X case 2:
X shdswp = .false.
X return
X
X # If k is ideal, this is like unto case 4.
X case 1:
X call acchk(i,j,h,shdswp) ; # if (ierr) return # This checks for convexity.
X return
X
X # If none of the `other' three corners are ideal, do the Lee-Schacter
X # test for the LOP.
X case 0:
X call qtest1(h,i,j,k,shdswp) ; # if (ierr) return
X return
X
X default: # This CAN'T happen.
X cerr = 'Indicator ijk is out of range in qtest.'
X ierr = 1 ; return
X
X }
X
X end
X # include defns.h
X subroutine qtest1(h,i,j,k,shdswp)
X # The Lee-Schacter test for the LOP (all points are
X # real, i.e. non-ideal). If the LOP is ***not***
X # satisfied then the diagonals should be swapped,
X # i.e. shdswp ("should-swap") is true.
X
X integer h
X logical shdswp
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # The vertices of the quadrilateral are labelled
X # h, i, j, k in the anticlockwise direction, h
X # being the point of central interest.
X
X # Make sure the quadrilateral is convex, so that
X # it makes sense to swap the diagonal.
X call acchk(i,j,k,shdswp) ; # if (ierr) return
X if(!shdswp) return
X
X # Get the coordinates of vertices h and j.
X xh = v(h,1)
X yh = v(h,2)
X xj = v(j,1)
X yj = v(j,2)
X
X # Find the centre of the circumcircle of vertices h, i, k.
X call circen(h,i,k,x0,y0,shdswp) ; # if (ierr) return
X if(shdswp) return # The points h, i, and k are colinear, so
X # the circumcircle has `infinite radius', so
X # (xj,yj) is definitely inside.
X
X # Check whether (xj,yj) is inside the circle of centre
X # (x0,y0) and radius r = dist[(x0,y0),(xh,yh)]
X
X a = x0-xh
X b = y0-yh
X r2 = a*a+b*b
X a = x0-xj
X b = y0-yj
X ch = a*a+b*b
X if(ch=xmax) {
X area = 0.
X return
X }
X
X # We're now looking at a trapezoidal region which may or may
X # not protrude above or below the horizontal strip bounded by
X # y = ymax and y = ymin.
X ybot = min(yl,yr)
X ytop = max(yl,yr)
X
X # Case 1; ymax <= ybot:
X # The `roof' of the trapezoid is entirely above the
X # horizontal strip.
X if(ymax<=ybot) {
X area = (xr-xl)*(ymax-ymin)
X return
X }
X
X # Case 2; ymin <= ybot <= ymax <= ytop:
X # The `roof' of the trapezoid intersects the top of the
X # horizontal strip (y = ymax) but not the bottom (y = ymin).
X if(ymin<=ybot&ymax<=ytop) {
X call testeq(slope,0.,value) ; # if (ierr) return
X if(value) {
X w1 = 0.
X w2 = xr-xl
X }
X else {
X xit = xl+(ymax-yl)/slope
X w1 = xit-xl
X w2 = xr-xit
X if(slope<0.) {
X tmp = w1
X w1 = w2
X w2 = tmp
X }
X }
X area = 0.5*w1*((ybot-ymin)+(ymax-ymin))+w2*(ymax-ymin)
X return
X }
X
X # Case 3; ybot <= ymin <= ymax <= ytop:
X # The `roof' intersects both the top (y = ymax) and
X # the bottom (y = ymin) of the horizontal strip.
X if(ybot<=ymin&ymax<=ytop) {
X xit = xl+(ymax-yl)/slope
X xib = xl+(ymin-yl)/slope
X if(slope>0.) {
X w1 = xit-xib
X w2 = xr-xit
X }
X else {
X w1 = xib-xit
X w2 = xit-xl
X }
X area = 0.5*w1*(ymax-ymin)+w2*(ymax-ymin)
X return
X }
X
X # Case 4; ymin <= ybot <= ytop <= ymax:
X # The `roof' is ***between*** the bottom (y = ymin) and
X # the top (y = ymax) of the horizontal strip.
X if(ymin<=ybot&ytop<=ymax) {
X area = 0.5*(xr-xl)*((ytop-ymin)+(ybot-ymin))
X return
X }
X
X # Case 5; ybot <= ymin <= ytop <= ymax:
X # The `roof' intersects the bottom (y = ymin) but not
X # the top (y = ymax) of the horizontal strip.
X if(ybot<=ymin&ymin<=ytop) {
X call testeq(slope,0.,value) ; # if (ierr) return
X if(value) {
X area = 0.
X return
X }
X xib = xl+(ymin-yl)/slope
X if(slope>0.) w = xr-xib
X else w = xib-xl
X area = 0.5*w*(ytop-ymin)
X return
X }
X
X # Case 6; ytop <= ymin:
X # The `roof' is entirely below the bottom (y = ymin), so
X # there is no area contribution at all.
X if(ytop<=ymin) {
X area = 0.
X return
X }
X
X # Default; all stuffed up:
X cerr = 'Fell through all six cases in routine ''stoke''.'
X cerr = 'Something must be totally stuffed up.'
X ierr = 1 ; return
X
X end
X # include defns.h
X integer function succ(i,j)
X # Find the successor of j in the adjacency list of i.
X
X integer succ
X include commons.h
X
X n = nadj(i,0)
X
X # If the adjacency list of i is empty, then clearly j has no successor
X # in this adjacency list. Something's wrong; stop.
X if(n==0) {
X # cerr = 'Adjacency of ',i,' is empty,'
X # cerr = 'and so cannot contain ',j,'.'
X # cerr = 'Error in ''succ''.'
X cerr = 'Adjacency problem in succ()'
X ierr = 1 ; return
X }
X
X # The adjacency list of i is non-empty; search through it until j is found;
X # add 1 to the location of j, and find the contents of this new location.
X do k = 1,n {
X if(j==nadj(i,k)) {
X kp = k+1
X if(kp>n) kp = 1 # Take kp modulo n. (The adjacency list
X succ = nadj(i,kp) # is circular.)
X return
X }
X }
X
X # The adjacency list doesn't contain j. Something's wrong; stop.
X # cerr = 'Adjacency list of ',i,' does not contain ',j,'.'
X # cerr = 'Error in ''succ''.'
X cerr = 'Adjacency problem in succ()'
X ierr = 1 ; return
X end
X # include defns.h
X subroutine swap(j,k1,k2,shdswp)
X # The segment k1->k2 is a diagonal of a quadrilateral
X # with a vertex at j (the point being added to the
X # triangulation). If the LOP is not satisfied, swap
X # it for the other diagonal.
X
X logical shdswp
X integer pred, succ
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # If vertices k1 and k2 are not connected there is no diagonal to swap.
X # This could happen if vertices j, k1, and k2 were colinear, but shouldn't.
X call adjchk(k1,k2,shdswp) ; # if (ierr) return
X if(!shdswp) return
X
X # Get the other vertex of the quadrilateral.
X k = pred(k1,k2) # If these aren't the same, then
X kk = succ(k2,k1) # there is no other vertex.
X if(kk!=k) {
X shdswp = .false.
X return
X }
X
X # Check whether the LOP is satisified; i.e. whether
X # vertex k is outside the circumcircle of vertices j, k1, and k2
X call qtest(j,k1,k,k2,shdswp) ; # if (ierr) return
X
X # Do the actual swapping.
X if(shdswp) {
X call delet(k1,k2) ; # if (ierr) return
X call insrt(j,k) ; # if (ierr) return
X }
X return
X end
X # include defns.h
X subroutine testeq(a,b,value)
X # Test for the equality of a and b in a fairly
X # robust way.
X logical value
X
X # Check for error code
X # if (ierr) return
X
X eps = EPS
X
X # If b is essentially 0, check whether a is essentially zero also.
X if(abs(b)<=eps) {
X if(abs(a)<=eps) value = .true.
X else value = .false.
X return
X }
X
X # Test if a is a `lot different' from b. (If it is
X # they're obviously not equal.) This avoids under/overflow
X # problems in dividing a by b.
X if(abs(a)>10.*abs(b)|abs(a)<0.1*abs(b)) {
X value = .false.
X return
X }
X
X # They're non-zero and fairly close; compare their ratio with 1.
X c = a/b
X if(abs(c-1.)<=eps) value = .true.
X else value = .false.
X
X return
X end
X subroutine triar(x0,y0,x1,y1,x2,y2,area)
X # Calculate the area of a triangle with given
X # vertices, in anti clockwise direction.
X
X # Check for error code
X # if (ierr) return
X
X area = 0.5*((x1-x0)*(y2-y0)-(x2-x0)*(y1-y0))
X return
X end
X # include defns.h
X subroutine trifnd(j,tau,nedge)
X # Find the triangle of the extant triangulation in which
X # lies the point currently being added.
X
X integer tau(3), pred, succ
X logical adjace, onedge
X include commons.h
X
X # Check for error code
X # if (ierr) return
X
X # The first point must be added to the triangulation before
X # calling trifnd.
X if(j==1) {
X cerr = 'No triangles to find.'
X ierr = 1 ; return
X }
X
X # Get the coordinates of point j
X x0 = v(j,1)
X y0 = v(j,2)
X j1 = j-1
X
X # Get the previous triangle
X tau(1) = j1
X tau(3) = nadj(j1,1)
X tau(2) = pred(j1,tau(3))
X call adjchk(tau(2),tau(3),adjace) ; # if (ierr) return
X if(!adjace) {
X tau(3) = tau(2)
X tau(2) = pred(j1,tau(3))
X }
X
X # Move to the adjacent triangle in the direction of the new
X # point, until the new point lies in this triangle.
X 1 continue
X ntau = 0 # This number will identify the triangle to be moved to.
X nedge = 0 # If the point lies on an edge, this number will identify that edge.
X do i = 1,3 {
X ip = i+1
X if(ip==4) ip = 1 # Take addition modulo 3.
X # Get the coordinates of the vertices of the current triangle,
X call coords(tau(i),xi,yi) ; # if (ierr) return
X if(tau(i)<=0) {
X xi = x0+xi
X yi = y0+yi
X }
X call coords(tau(ip),xip,yip) ; # if (ierr) return
X if(tau(ip)<=0) {
X xip = x0+xip
X yip = y0+yip
X }
X
X # If `test' is positive the point being added is to the left
X # (as we move along the edge in an anti-clockwise direction).
X # If this is true for all three edges, then the point is inside
X # the triangle.
X #test = (y0-yi)*(xip-x0)-(x0-xi)*(yip-y0)
X test = (xip-xi)*(y0-yi)-(yip-yi)*(x0-xi)
X call testeq(test,0.,onedge) ; # if (ierr) return
X # If `test' is very close to zero, we might get negative values
X # for it on both sides of an edge, and hence go into an infinite
X # loop.
X if(onedge) nedge = ip
X else if(test<0.) {
X ntau = ip
X break
X }
X }
X switch(ntau) {
X case 0: # All tests >= 0.; the point is inside; return.
X return
X
X # The point is not inside; work out the vertices of the triangle to which
X # to move. Notation: Number the vertices of the current triangle from 1 to 3,
X # anti-clockwise. Then "triangle i+1" is adjacent to the side from vertex i to
X # vertex i+1, where i+1 is taken modulo 3 (i.e. "3+1 = 1").
X case 1:
X # move to "triangle 1"
X #tau(1) = tau(1)
X tau(2) = tau(3)
X tau(3) = succ(tau(1),tau(2))
X case 2:
X # move to "triangle 2"
X #tau(1) = tau(1)
X tau(3) = tau(2)
X tau(2) = pred(tau(1),tau(3))
X case 3:
X # move to "triangle 3"
X tau(1) = tau(3)
X #tau(2) = tau(2)
X tau(3) = succ(tau(1),tau(2))
X }
X
X # We've moved to a new triangle; check if the point being added lies
X # inside this one.
X go to 1
X
X end
END_OF_FILE
if test 47239 -ne `wc -c <'deldir.r'`; then
echo shar: \"'deldir.r'\" unpacked with wrong size!
fi
# end of 'deldir.r'
fi
if test -f 'deldir.sp' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'deldir.sp'\"
else
echo shar: Extracting \"'deldir.sp'\" \(4506 characters\)
sed "s/^X//" >'deldir.sp' <<'END_OF_FILE'
X# 9/9/94
X# Donald H. MacQueen
X# Lawrence Livermore National Laboratory
X# macq@llnl.gov
X#
X# Permission to use, copy, modify, and distribute this software and
X# its documentation for any purpose and without fee is hereby
X# granted, provided that this permission
X# notice appears in all copies and in supporting documentation.
X#
X# S function to perfrom delaunay triangulation and dirichlet tesselation
X# Using fortran code based on the delaunay library by Rolf Turner in statlib
X# See the fortan code for additional information
X
X# Required arguments: x,y two numeric vectors of same length
X# Optional arguments: ndx,ndy
X# Optional arguments: ndx,ndy,xmin,xmax,ymin,ymax
X# See Turner's documentation for description of optional arguments
X
X# Value:
X# The delaunay triangulation as a data frame
X# The dirichlet tesselation as a data frame
X# The numbers of segements in each
X# Vectors of areas of the triangles and polygons
X# The total areas of the triangles and polygons
X# Vectors of weights associated with the triangles and polygons
X# x and y limits useful for plotting the triangles or polygons
X
X# Also included with this distribution is a replacement for the Splus
X# supplied segments() function. This replacement will take a data frame
X# or list argument, and plot segments from its x1, y1, x2, and y2
X# components. Place this function in your search() path ahead of
X# the supplied segments() function if you wish to use it. Then
X# you can plot the triangulation or tesselation like this:
X# tmp <- deldir(x,y)
X# plot(x,y,xlim=tmp$xlim,ylim=tmp$ylim)
X# segments(tmp$triang,col=2)
X# segments(tmp$polyg,col=3)
X# Using the built in segments() the calls are
X# segments(tmp$triang$x1,tmp$triang$y1,tmp$triang$x2,tmp$triang$y2,col=2)
X# and similarly for the polygons.
X
Xdeldir <- function(x,y,ndx=0,ndy=0,xmin=0,xmax=0,ymin=0,ymax=0,narg=0,npts=length(x)) {
X
Xif (!is.loaded(symbol.For("deldir")))
X print(dyn.load( paste(getenv('SHOME'),'/library/delaunay/deldir.o',sep='') ))
X
Xif (any(is.na(c(x,y)))) return('NAs not allowed in x and y')
X
X# check data, set up arguments
Xif (length(y)!=npts) return('Lengths of x and y do not match')
Xif (npts>800) return('Maximum length(x) is 800, increase and recompile')
Xif (all(c(ndx,ndy)>0)) {
X narg <- 2
X if ( xmax > xmin & ymax > ymin ) narg <- 6
X }
X xin <- as.single(x)
X yin <- as.single(y)
X ndxin <- as.integer(ndx)
X ndyin <- as.integer(ndy)
X xminin <- as.single(xmin)
X xmaxin <- as.single(xmax)
X yminin <- as.single(ymin)
X ymaxin <- as.single(ymax)
X narg <- as.integer(narg)
X nptsin <- as.integer(npts)
X triang <- single(4*npts*(npts-1)/2-4)
X ntseg <- integer(1)
X tarea <- single(npts)
X tareatot <- single(1)
X polyg <- single(4*npts*(npts-1)/2-4)
X npseg <- integer(1)
X parea <- single(npts)
X pareatot <- single(1)
X errmsg <- character(1)
X ierro <- integer(1)
X
X retn.full <- .Fortran("deldir",
X xin,yin,
X ndxin,ndyin,
X xminin,xmaxin,yminin,ymaxin,
X narg,nptsin,
X triang=triang,
X ntseg=ntseg,
X tarea=tarea,
X tareatot=tareatot,
X polyg=polyg,
X npseg=npseg,
X parea=parea,
X pareatot=pareatot,
X ierr=ierro
X )
X
Xif (F) {
X retn <- list(triang=matrix(retn.full$triang[1:4*retn.full$ntseg],ncol=4,byrow=T),
X ntseg=retn.full$ntseg/4,
X tarea=retn.full$tarea,
X tareatot=retn.full$tareatot,
X twts=retn.full$tarea/retn.full$tareatot,
X polyg=matrix(retn.full$polyg[1:4*retn.full$npseg],ncol=4,byrow=T),
X npseg=retn.full$npseg/4,
X parea=retn.full$parea,
X pareatot=retn.full$pareatot,
X pwts=retn.full$parea/retn.full$pareatot
X )
X }
X retn <- list(triang=matrix(retn.full$triang[1:retn.full$ntseg],ncol=4,byrow=T),
X ntseg=retn.full$ntseg/4,
X tarea=retn.full$tarea,
X tareatot=retn.full$tareatot,
X twts=retn.full$tarea/retn.full$tareatot,
X polyg=matrix(retn.full$polyg[1:retn.full$npseg],ncol=4,byrow=T),
X npseg=retn.full$npseg/4,
X parea=retn.full$parea,
X pareatot=retn.full$pareatot,
X pwts=retn.full$parea/retn.full$pareatot
X )
X
X retn$xlim <- range(retn$polyg[,1],retn$polyg[,3],xin)
X retn$ylim <- range(retn$polyg[,2],retn$polyg[,4],yin)
X
X dimnames(retn$triang) <- list(NULL,c('x1','y1','x2','y2'))
X dimnames(retn$polyg) <- list(NULL,c('x1','y1','x2','y2'))
X retn$triang <- data.frame(retn$triang)
X retn$polyg <- data.frame(retn$polyg)
X
X retn
X}
X
END_OF_FILE
if test 4506 -ne `wc -c <'deldir.sp'`; then
echo shar: \"'deldir.sp'\" unpacked with wrong size!
fi
# end of 'deldir.sp'
fi
if test -f 'deldir.sp~' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'deldir.sp~'\"
else
echo shar: Extracting \"'deldir.sp~'\" \(4260 characters\)
sed "s/^X//" >'deldir.sp~' <<'END_OF_FILE'
X# 9/9/94
X# Donald H. MacQueen
X# Lawrence Livermore National Laboratory
X# macq@llnl.gov
X#
X# S function to perfrom delaunay triangulation and dirichlet tesselation
X# Using fortran code based on the delaunay library by Rolf Turner in statlib
X# See the fortan code for additional information
X
X# Required arguments: x,y two numeric vectors of same length
X# Optional arguments: ndx,ndy
X# Optional arguments: ndx,ndy,xmin,xmax,ymin,ymax
X# See Turner's documentation for description of optional arguments
X
X# Value:
X# The delaunay triangulation as a data frame
X# The dirichlet tesselation as a data frame
X# The numbers of segements in each
X# Vectors of areas of the triangles and polygons
X# The total areas of the triangles and polygons
X# Vectors of weights associated with the triangles and polygons
X# x and y limits useful for plotting the triangles or polygons
X
X# Also included with this distribution is a replacement for the Splus
X# supplied segments() function. This replacement will take a data frame
X# or list argument, and plot segments from its x1, y1, x2, and y2
X# components. Place this function in your search() path ahead of
X# the supplied segments() function if you wish to use it. Then
X# you can plot the triangulation or tesselation like this:
X# tmp <- deldir(x,y)
X# plot(x,y,xlim=tmp$xlim,ylim=tmp$ylim)
X# segments(tmp$triang,col=2)
X# segments(tmp$polyg,col=3)
X# Using the built in segments() the calls are
X# segments(tmp$triang$x1,tmp$triang$y1,tmp$triang$x2,tmp$triang$y2,col=2)
X# and similarly for the polygons.
X
Xdeldir <- function(x,y,ndx=0,ndy=0,xmin=0,xmax=0,ymin=0,ymax=0,narg=0,npts=length(x)) {
X
Xif (!is.loaded(symbol.For("deldir")))
X print(dyn.load( paste(getenv('SHOME'),'/library/delaunay/deldir.o',sep='') ))
X
Xif (any(is.na(c(x,y)))) return('NAs not allowed in x and y')
X
X# check data, set up arguments
Xif (length(y)!=npts) return('Lengths of x and y do not match')
Xif (npts>800) return('Maximum length(x) is 800, increase and recompile')
Xif (all(c(ndx,ndy)>0)) {
X narg <- 2
X if ( xmax > xmin & ymax > ymin ) narg <- 6
X }
X xin <- as.single(x)
X yin <- as.single(y)
X ndxin <- as.integer(ndx)
X ndyin <- as.integer(ndy)
X xminin <- as.single(xmin)
X xmaxin <- as.single(xmax)
X yminin <- as.single(ymin)
X ymaxin <- as.single(ymax)
X narg <- as.integer(narg)
X nptsin <- as.integer(npts)
X triang <- single(4*npts*(npts-1)/2-4)
X ntseg <- integer(1)
X tarea <- single(npts)
X tareatot <- single(1)
X polyg <- single(4*npts*(npts-1)/2-4)
X npseg <- integer(1)
X parea <- single(npts)
X pareatot <- single(1)
X errmsg <- character(1)
X ierro <- integer(1)
X
X retn.full <- .Fortran("deldir",
X xin,yin,
X ndxin,ndyin,
X xminin,xmaxin,yminin,ymaxin,
X narg,nptsin,
X triang=triang,
X ntseg=ntseg,
X tarea=tarea,
X tareatot=tareatot,
X polyg=polyg,
X npseg=npseg,
X parea=parea,
X pareatot=pareatot,
X ierr=ierro
X )
X
Xif (F) {
X retn <- list(triang=matrix(retn.full$triang[1:4*retn.full$ntseg],ncol=4,byrow=T),
X ntseg=retn.full$ntseg/4,
X tarea=retn.full$tarea,
X tareatot=retn.full$tareatot,
X twts=retn.full$tarea/retn.full$tareatot,
X polyg=matrix(retn.full$polyg[1:4*retn.full$npseg],ncol=4,byrow=T),
X npseg=retn.full$npseg/4,
X parea=retn.full$parea,
X pareatot=retn.full$pareatot,
X pwts=retn.full$parea/retn.full$pareatot
X )
X }
X retn <- list(triang=matrix(retn.full$triang[1:retn.full$ntseg],ncol=4,byrow=T),
X ntseg=retn.full$ntseg/4,
X tarea=retn.full$tarea,
X tareatot=retn.full$tareatot,
X twts=retn.full$tarea/retn.full$tareatot,
X polyg=matrix(retn.full$polyg[1:retn.full$npseg],ncol=4,byrow=T),
X npseg=retn.full$npseg/4,
X parea=retn.full$parea,
X pareatot=retn.full$pareatot,
X pwts=retn.full$parea/retn.full$pareatot
X )
X
X retn$xlim <- range(retn$polyg[,1],retn$polyg[,3],xin)
X retn$ylim <- range(retn$polyg[,2],retn$polyg[,4],yin)
X
X dimnames(retn$triang) <- list(NULL,c('x1','y1','x2','y2'))
X dimnames(retn$polyg) <- list(NULL,c('x1','y1','x2','y2'))
X retn$triang <- data.frame(retn$triang)
X retn$polyg <- data.frame(retn$polyg)
X
X retn
X}
X
END_OF_FILE
if test 4260 -ne `wc -c <'deldir.sp~'`; then
echo shar: \"'deldir.sp~'\" unpacked with wrong size!
fi
# end of 'deldir.sp~'
fi
if test -f 'segments.sp' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'segments.sp'\"
else
echo shar: Extracting \"'segments.sp'\" \(788 characters\)
sed "s/^X//" >'segments.sp' <<'END_OF_FILE'
X# 9/9/94
X# Alternate segments() function, to allow
X# list or data frame argument. Expects
X# the list or data frame to have objects
X# named x1, y1, x2, and y2
X#
X# Permission to use, copy, modify, and distribute this software and
X# its documentation for any purpose and without fee is hereby
X# granted, provided that this permission
X# notice appears in all copies and in supporting documentation.
X
Xsegments <- function(x,...) {
X if ( is.list(x) | is.data.frame(x) ) {
X if ( all(sort(match(c('x1','y1','x2','y2'),names(x)))==1:4) ) {
X .Begin.pic()
X .Cur.pic(.S(segments(x$x1,x$y1,x$x2,x$y2,...), "segments"))
X }
X else return('Data frame or list argument missing x1, y1, x2, or y2')
X }
X else {
X .Begin.pic()
X .Cur.pic(.S(segments(x,...), "segments"))
X }
X }
X
END_OF_FILE
if test 788 -ne `wc -c <'segments.sp'`; then
echo shar: \"'segments.sp'\" unpacked with wrong size!
fi
# end of 'segments.sp'
fi
if test -f 'segments.sp~' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'segments.sp~'\"
else
echo shar: Extracting \"'segments.sp~'\" \(542 characters\)
sed "s/^X//" >'segments.sp~' <<'END_OF_FILE'
X# 9/9/94
X# Alternate segments() function, to allow
X# list or data frame argument. Expects
X# the list or data frame to have objects
X# named x1, y1, x2, and y2
X
Xsegments <- function(x,...) {
X if ( is.list(x) | is.data.frame(x) ) {
X if ( all(sort(match(c('x1','y1','x2','y2'),names(x)))==1:4) ) {
X .Begin.pic()
X .Cur.pic(.S(segments(x$x1,x$y1,x$x2,x$y2,...), "segments"))
X }
X else return('Data frame or list argument missing x1, y1, x2, or y2')
X }
X else {
X .Begin.pic()
X .Cur.pic(.S(segments(x,...), "segments"))
X }
X }
X
END_OF_FILE
if test 542 -ne `wc -c <'segments.sp~'`; then
echo shar: \"'segments.sp~'\" unpacked with wrong size!
fi
# end of 'segments.sp~'
fi
echo shar: End of shell archive.
exit 0