#! /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