chrisp@bilby.cs.uwa.oz.au (Chris Pudney) (03/26/91)
G'day, I require an algorithm for plotting the contour lines on a surface. The surfaces are basically 2-D potential fields (f:R^2 -> R) which I would like to depict graphically using contour lines (or perhaps I should use greyscale/colour?) Could you provide me with pointers to a contour plotting algorithm. Public domain (ftp'able) source would be most desirable. If I get a large enough response I will post a compilation of the information I receive. AtDhVaAnNkCsE. -------------------------------------------------------------------------------- Chris Pudney Department of Computer Science, PHONE: (local: 09) (int'l: +61 9) 380 3455 University of Western Australia, Nedlands, Western AUSTRALIA, 6009. FAX: (local: 09) (int'l: +61 9) 382 1688 E-MAIL: chrisp@cs.uwa.oz.au -- -------------------------------------------------------------------------------- Chris Pudney Department of Computer Science, PHONE: (local: 09) (int'l: +61 9) 380 3455
jdm5548@tamsun.tamu.edu (James Darrell McCauley) (03/27/91)
This was posted by marwk@levels.sait.edu.au on Tue Dec 11 22:01:20 1990. Darrell -- James Darrell McCauley (jdm5548@diamond.tamu.edu, jdm5548@tamagen.bitnet) Spatial Analysis Lab, Department of Agricultural Engineering, Texas A&M University, College Station, Texas 77843-2117, USA --------------------------- These are the references I already had before posting to NEWS. CARTOGRAPHIC JOURNAL: 1977, Vol 14, Pt. 2: "Computer ...", by Yoeli, P. 1979, 16, 1: "Contouring ...", Elfick, M.H. 1980, 17, 2: "Triangulation ...", by McCullagh et al. "Construction of a digital terrain model using bicubic polynomials", A.C. Jones, Mathematical Project G, 1978. (This has FORTRAN code but I do not know if I can make it public yet - I need to request permission for this.) --------------------------- P.G.Hamer@stl.stc.co.uk The ACM algorithm collection should have something to start you off. The mail server netlib provides quite a few of them, start by looking in the index toms. Algorithm 626 Tricp might be close to what you are looking for. It is in fortran; if this hurts you netlib also provides f2c a fortran to C converter. --------------------------- AVINOAM%HUJIVMS.BITNET@CUNYVM.CUNY.EDU After using Akima's method (or any other method) for genereting a regular grid, you can use the following algorithm: contour plotting Collected Algorithms, from the ACM Transactions on Mathematical Software. algorithm 531 ACM Transactions on Mathematical Software volume 4 issue number 3 p. 290 The program follows. I did not check it my self (maybe you should change some of the machine dependent constants, like word length. see the paper for details.) You can get the program for netlib yourself, if you want... Good luck. Avinoam@hujivms, avinoam@shum From: HUJICC::BITNET%"netlibd@research.att.com" 17-FEB-1990 09:46:49.96 To: avinoam@HUJIVMS <= netlibd@research.att.com (M NETLIBD.MAIL ASCII) CC: Subj: Re:Subject: send 531 from toms Disclaimer. This on-line distribution of algorithms is not an official service of the Association for Computing Machinery (ACM). This service is being provided on a trial basis and may be discontinued in the future. ACM is not responsible for any transmission errors in algorithms received by on-line distribution. Official ACM versions of algorithms may be obtained from: IMSL, Inc. 2500 ParkWest Tower One 2500 CityWest Blvd. Houston, TX 77042-3020 (713) 782-6060 Columns 73-80 have been deleted, trailing blanks removed, and program text mapped to lower case. Assembly language programs have sometimes been deleted and standard machine constant and error handling routines removed from the individual algorithms and collected into the libraries "core" and "slatec". In a few cases, this trimming was overzealous, and will be repaired when we can get back the original files. We have tried to incorporate published Remarks; if we missed something, please contact Eric Grosse, 201-582-5828, research!ehg or ehg@research.att.com. The material in this library is copyrighted by the ACM, which grants general permission to distribute provided the copies are not made for direct commercial advantage. For details of the copyright and dissemination agreement, consult the current issue of TOMS. Careful! Anything free comes with no guarantee. *** from netlib, Fri Feb 16 07:08:08 EST 1990 *** subroutine fill0(bitmap, n) c c fill the first n bits of bitmap with zeroes. c integer bitmap(1), n c data nbpw /35/ c nbpw is the minimum number of significant bits per word used c by integer arithmetic. this is usually one less than the c actual number of bits per word, but an important exception is c the cdc-6000 series of machines, where nbpw should be 48. c loop = n/nbpw nblw = mod(n,nbpw) if (loop.eq.0) go to 20 do 10 i=1,loop bitmap(i) = 0 10 continue 20 if (nblw.ne.0) bitmap(loop+1) = mod(bitmap(loop+1),2**(nbpw-nblw) * ) return end subroutine mark1(bitmap, n) c c put a one in the nth bit of bitmap. c integer bitmap(1), n c data nbpw /35/ c nbpw is the minimum number of significant bits per word used c by integer arithmetic. this is usually one less than the c actual number of bits per word, but an important exception is c the cdc-6000 series of machines, where nbpw should be 48. c nword = (n-1)/nbpw nbit = mod(n-1,nbpw) i = 2**(nbpw-nbit-1) bitmap(nword+1) = bitmap(nword+1) + i*(1-mod(bitmap(nword+1)/i,2) * ) return end function iget(bitmap, n) c c iget=0 if the nth bit of bitmap is zero, else iget is one. c integer bitmap(1), n c data nbpw /35/ c nbpw is the minimum number of significant bits per word used c by integer arithmetic. this is usually one less than the c actual number of bits per word, but an important exception is c the cdc-6000 series of machines, where nbpw should be 48. c nword = (n-1)/nbpw nbit = mod(n-1,nbpw) iget = mod(bitmap(nword+1)/2**(nbpw-nbit-1),2) return end subroutine gcontr(z, nrz, nx, ny, cv, ncv, zmax, bitmap, draw) c c this subroutine draws a contour through equal values of an array. c c ***** formal arguments *********************************** c c z is the array for which contours are to be drawn. the elements c of z are assumed to lie upon the nodes of a topologically c rectangular coordinate system - e.g. cartesian, polar (except c the origin), etc. c c nrz is the number of rows declared for z in the calling program. c c nx is the limit for the first subscript of z. c c ny is the limit for the second subscript of z. c c cv are the values of the contours to be drawn. c c ncv is the number of contour values in cv. c c zmax is the maximum value of z for consideration. a value of c z(i,j) greater than zmax is a signal that that point and the c grid line segments radiating from that point to it's neighbors c are to be excluded from contouring. c c bitmap is a work area large enough to hold 2*nx*ny*ncv bits. it c is accessed by low-level routines, which are described below. c let j be the number of useful bits in each word of bitmap, c as determined by the user machine and implementation of c the bitmap manipulation subprograms described below. then c the number of words required for the bitmap is the floor of c (2*nx*ny*ncv+j-1)/j. c c draw is a user-provided subroutine used to draw contours. c the calling sequence for draw is: c c call draw (x,y,iflag) c let nx = integer part of x, fx = fractional part of x. c then x should be interpreted such that increases in nx c correspond to increases in the first subscript of z, and c fx is the fractional distance from the abscissa corresponding c to nx to the abscissa corresponding to nx+1, c and y should be interpreted similarly for the second c subscript of z. c the low-order digit of iflag will have one of the values: c 1 - continue a contour, c 2 - start a contour at a boundary, c 3 - start a contour not at a boundary, c 4 - finish a contour at a boundary, c 5 - finish a closed contour (not at a boundary). c note that requests 1, 4 and 5 are for pen-down c moves, and that requests 2 and 3 are for pen-up c moves. c 6 - set x and y to the approximate 'pen' position, using c the notation discussed above. this call may be c ignored, the result being that the 'pen' position c is taken to correspond to z(1,1). c iflag/10 is the contour number. c c ***** external subprograms ******************************* c c draw is the user-supplied line drawing subprogram described above. c draw may be sensitive to the host computer and to the plot device. c fill0 is used to fill a bitmap with zeroes. call fill0 (bitmap,n) c fills the first n bits of bitmap with zeroes. c mark1 is used to place a 1 in a specific bit of the bitmap. c call mark1 (bitmap,n) puts a 1 in the nth bit of the bitmap. c iget is used to determine the setting of a particular bit in the c bitmap. i=iget(bitmap,n) sets i to zero if the nth bit of the c bitmap is zero, and sets i to one if the nth bit is one. c fill0, mark1 and iget are machine sensitive. c c ****************************************************************** c real z(nrz,1), cv(1) integer bitmap(1) integer l1(4), l2(4), ij(2) c c l1 and l2 contain limits used during the spiral search for the c beginning of a contour. c ij stores subcripts used during the spiral search. c integer i1(2), i2(2), i3(6) c c i1, i2 and i3 are used for subscript computations during the c examination of lines from z(i,j) to it's neighbors. c real xint(4) c c xint is used to mark intersections of the contour under c consideration with the edges of the cell being examined. c real xy(2) c c xy is used to compute coordinates for the draw subroutine. c equivalence (l2(1),imax), (l2(2),jmax), (l2(3),imin), * (l2(4),jmin) equivalence (ij(1),i), (ij(2),j) equivalence (xy(1),x), (xy(2),y) c data l1(3) /-1/, l1(4) /-1/ data i1 /1,0/, i2 /1,-1/, i3 /1,0,0,1,1,0/ c l1(1) = nx l1(2) = ny dmax = zmax c c set the current pen position. the default position corresponds c to z(1,1). c x = 1.0 y = 1.0 call draw(x, y, 6) icur = max0(1,min0(int(x),nx)) jcur = max0(1,min0(int(y),ny)) c c clear the bitmap c call fill0(bitmap, 2*nx*ny*ncv) c c search along a rectangular spiral path for a line segment having c the following properties: c 1. the end points are not excluded, c 2. no mark has been recorded for the segment, c 3. the values of z at the ends of the segment are such that c one z is less than the current contour value, and the c other is greater than or equal to the current contour c value. c c search all boundaries first, then search interior line segments. c note that the interior line segments near excluded points may be c boundaries. c ibkey = 0 10 i = icur j = jcur 20 imax = i imin = -i jmax = j jmin = -j idir = 0 c direction zero is +i, 1 is +j, 2 is -i, 3 is -j. 30 nxidir = idir + 1 k = nxidir if (nxidir.gt.3) nxidir = 0 40 i = iabs(i) j = iabs(j) if (z(i,j).gt.dmax) go to 140 l = 1 c l=1 means horizontal line, l=2 means vertical line. 50 if (ij(l).ge.l1(l)) go to 130 ii = i + i1(l) jj = j + i1(3-l) if (z(ii,jj).gt.dmax) go to 130 assign 100 to jump c the next 15 statements (or so) detect boundaries. 60 ix = 1 if (ij(3-l).eq.1) go to 80 ii = i - i1(3-l) jj = j - i1(l) if (z(ii,jj).gt.dmax) go to 70 ii = i + i2(l) jj = j + i2(3-l) if (z(ii,jj).lt.dmax) ix = 0 70 if (ij(3-l).ge.l1(3-l)) go to 90 80 ii = i + i1(3-l) jj = j + i1(l) if (z(ii,jj).gt.dmax) go to 90 if (z(i+1,j+1).lt.dmax) go to jump, (100, 280) 90 ix = ix + 2 go to jump, (100, 280) 100 if (ix.eq.3) go to 130 if (ix+ibkey.eq.0) go to 130 c now determine whether the line segment is crossed by the contour. ii = i + i1(l) jj = j + i1(3-l) z1 = z(i,j) z2 = z(ii,jj) do 120 icv=1,ncv if (iget(bitmap,2*(nx*(ny*(icv-1)+j-1)+i-1)+l).ne.0) go to 120 if (cv(icv).le.amin1(z1,z2)) go to 110 if (cv(icv).le.amax1(z1,z2)) go to 190 110 call mark1(bitmap, 2*(nx*(ny*(icv-1)+j-1)+i-1)+l) 120 continue 130 l = l + 1 if (l.le.2) go to 50 140 l = mod(idir,2) + 1 ij(l) = isign(ij(l),l1(k)) c c lines from z(i,j) to z(i+1,j) and z(i,j+1) are not satisfactory. c continue the spiral. c 150 if (ij(l).ge.l1(k)) go to 170 ij(l) = ij(l) + 1 if (ij(l).gt.l2(k)) go to 160 go to 40 160 l2(k) = ij(l) idir = nxidir go to 30 170 if (idir.eq.nxidir) go to 180 nxidir = nxidir + 1 ij(l) = l1(k) k = nxidir l = 3 - l ij(l) = l2(k) if (nxidir.gt.3) nxidir = 0 go to 150 180 if (ibkey.ne.0) return ibkey = 1 go to 10 c c an acceptable line segment has been found. c follow the contour until it either hits a boundary or closes. c 190 iedge = l cval = cv(icv) if (ix.ne.1) iedge = iedge + 2 iflag = 2 + ibkey xint(iedge) = (cval-z1)/(z2-z1) 200 xy(l) = float(ij(l)) + xint(iedge) xy(3-l) = float(ij(3-l)) call mark1(bitmap, 2*(nx*(ny*(icv-1)+j-1)+i-1)+l) call draw(x, y, iflag+10*icv) if (iflag.lt.4) go to 210 icur = i jcur = j go to 20 c c continue a contour. the edges are numbered clockwise with c the bottom edge being edge number one. c 210 ni = 1 if (iedge.lt.3) go to 220 i = i - i3(iedge) j = j - i3(iedge+2) 220 do 250 k=1,4 if (k.eq.iedge) go to 250 ii = i + i3(k) jj = j + i3(k+1) z1 = z(ii,jj) ii = i + i3(k+1) jj = j + i3(k+2) z2 = z(ii,jj) if (cval.le.amin1(z1,z2)) go to 250 if (cval.gt.amax1(z1,z2)) go to 250 if (k.eq.1) go to 230 if (k.ne.4) go to 240 230 zz = z1 z1 = z2 z2 = zz 240 xint(k) = (cval-z1)/(z2-z1) ni = ni + 1 ks = k 250 continue if (ni.eq.2) go to 260 c c the contour crosses all four edges of the cell being examined. c choose the lines top-to-left and bottom-to-right if the c interpolation point on the top edge is less than the interpolation c point on the bottom edge. otherwise, choose the other pair. this c method produces the same results if the axes are reversed. the c contour may close at any edge, but must not cross itself inside c any cell. c ks = 5 - iedge if (xint(3).lt.xint(1)) go to 260 ks = 3 - iedge if (ks.le.0) ks = ks + 4 c c determine whether the contour will close or run into a boundary c at edge ks of the current cell. c 260 l = ks iflag = 1 assign 280 to jump if (ks.lt.3) go to 270 i = i + i3(ks) j = j + i3(ks+2) l = ks - 2 270 if (iget(bitmap,2*(nx*(ny*(icv-1)+j-1)+i-1)+l).eq.0) go to 60 iflag = 5 go to 290 280 if (ix.ne.0) iflag = 4 290 iedge = ks + 2 if (iedge.gt.4) iedge = iedge - 4 xint(iedge) = xint(ks) go to 200 c end dimension z(51,51), c(10), work(1680) c dimension of work is large enough to contain c 2*(dimension of c)*(total dimension of z) useful bits. see the c bitmap routines accessed by gcontr. real mu external draw common /cur/ xcur, ycur data c(1), c(2), c(3), c(4), c(5) /3.05,3.2,3.5,3.50135,3.6/ data c(6), c(7), c(8), c(9), c(10) /3.766413,4.0,4.130149,5.0, * 10.0/ data nx /51/, ny /51/, nf /10/ data xmin /-2.0/, xmax /2.0/, ymin /-2.0/, ymax /2.0/, mu /0.3/ dx = (xmax-xmin)/float(nx-1) dy = (ymax-ymin)/float(ny-1) xcur = 1.0 ycur = 1.0 if (mod(nx,2).ne.0) ycur = float(ny) if (mod(ny,2).ne.0) xcur = float(nx) x = xmin - dx do 20 i=1,nx y = ymin - dy x = x + dx do 10 j=1,ny y = y + dy z(i,j) = (1.0-mu)*(2.0/sqrt((x-mu)**2+y**2)+(x-mu)**2+y**2) * + mu*(2.0/sqrt((x+1.0-mu)**2+y**2)+(x+1.0-mu)**2+y**2) 10 continue 20 continue call gcontr(z, 51, nx, ny, c, nf, 1.e6, work, draw) stop end real z(51,51), c(10), cval(10), mu integer work(1680), l(10), clab(10) c dimension of work is large enough to contain c 2*(dimension of c)*(total dimension of z) useful bits. see the c bitmap routines accessed by gcontr. external draw common /gctcom/ xcur, ycur, xl, yl, cval, clab, nch data c(1), c(2), c(3), c(4), c(5) /3.05,3.2,3.5,3.50135,3.6/ data c(6), c(7), c(8), c(9), c(10) /3.766413,4.0,4.130149,5.0, * 10.0/ data l(1), l(2), l(3), l(4), l(5) /1ha,1hb,1hc,1hd,1he/ data l(6), l(7), l(8), l(9), l(10) /1hf,1hg,1hh,1hi,1hj/ data nx /51/, ny /51/, nf /10/, nxg /5/, nyg /5/ data xmin /-2.0/, xmax /2.0/, ymin /-2.0/, ymax /2.0/, mu /0.3/ data xlen /8.0/, ylen /8.0/ c initialize plotting subroutines. call plots dx = (xmax-xmin)/float(nx-1) dy = (ymax-ymin)/float(ny-1) xl = xlen/float(nx) yl = ylen/float(ny) xcur = 1.0 ycur = 1.0 if (mod(nx,2).ne.0) ycur = float(ny) if (mod(ny,2).ne.0) xcur = float(nx) x = xmin - dx do 20 i=1,nx y = ymin - dy x = x + dx do 10 j=1,ny y = y + dy c evaluate function to be plotted. z(i,j) = (1.0-mu)*(2.0/sqrt((x-mu)**2+y**2)+(x-mu)**2+y**2) * + mu*(2.0/sqrt((x+1.0-mu)**2+y**2)+(x+1.0-mu)**2+y**2) 10 continue 20 continue do 30 i=1,nf cval(i) = c(i) clab(i) = l(i) 30 continue nch = 1 c pen up move to below lower left corner of page. c this call works differently on different machines. you may c need to change it. call plot(0.0, -11.0, -3) c pen up move to 1 inch above lower left corner of page. call plot(0.0, 1.0, -3) sx = 8.0/float(nxg) sy = 8.0/float(nxg) c draw a grid. call cgrid(1, nxg, sx, 0.0, 0.0, nyg, sy, 0.0, 0.0) c draw the contour plots. call gcontr(z, 51, nx, ny, cval, nf, 1.0e6, work, draw) xx = 9.0 yy = 8.0 c write a table of contour labels and values. call symbol(xx, yy+0.14, 0.07, 10hcontour id, 0.0, 10) do 40 i=1,nf call symbol(xx, yy, 0.07, l(i), 0.0, 2) call number(xx+0.12, yy, 0.07, c(i), 0.0, 5) yy = yy - 0.14 40 continue c pen up move to below lower right corner of page. c this call works differently on different machines. you may need c to change it, or you may not need it. call plot(10.0, -11.0, -3) c reduce picture size, plot end of file information. c the end of file information may not be available at all sites. c if not available, change the next two statements to comments. call factor(0.3) call plot(0.0, 0.0, 999) stop c end subroutine draw(x, y, iflag) c this subroutine uses calcomp plot routines to draw lines for the c contour plotting routine gcontr. real cval(10) integer clab(10) common /gctcom/ xcur, ycur, xl, yl, cval, clab, nch data iblank /1h / ih = iflag/10 il = iflag - 10*ih if (il.eq.6) go to 40 ipen = 2 if (il.eq.2) ipen = 3 if (il.eq.3) ipen = 3 xcur = x ycur = y xx = (x-1.0)*xl yy = (y-1.0)*yl call plot(xx, yy, ipen) if (il.lt.2) go to 30 if (il.gt.4) go to 30 if (nch.lt.1) go to 30 if (clab(ih).eq.iblank) go to 30 if (clab(ih).ne.0) go to 10 call number(xx, yy-0.03, 0.07, cval(ih), 0.0, -1) go to 20 10 call symbol(xx, yy-0.03, 0.07, clab(ih), 0.0, nch) 20 call plot(xx, yy, 3) 30 return 40 x = xcur y = ycur return c end subroutine draw(x, y, iflag) c c do output for gcontr. c integer print common /cur/ xcur, ycur data print /6/ c print is the system printer fortran i/o unit number. icont = iflag/10 jump = mod(iflag,10) go to (10, 20, 30, 40, 50, 60), jump 10 write (print,99999) icont, x, y go to 70 20 write (print,99998) icont, x, y go to 70 30 write (print,99997) icont, x, y go to 70 40 write (print,99996) icont, x, y go to 70 50 write (print,99995) icont, x, y go to 70 60 write (print,99994) x = xcur y = ycur 70 return 99999 format (17h continue contour, i3, 3h to, 1p2e14.7) 99998 format (14h start contour, i3, 19h on the boundary at, 1p2e14.7) 99997 format (14h start contour, i3, 19h in the interior at, 1p2e14.7) 99996 format (15h finish contour, i3, 19h on the boundary at, 1p2e14.7) 99995 format (15h finish contour, i3, 19h in the interior at, 1p2e14.7) 99994 format (33h request for current pen position) end subroutine cgrid(nopt, nx, sx, xs, xf, ny, sy, ys, yf) c c subroutine which draws a frame around the plot and draws c either tick marks or grid lines. c c parameters: nopt -- =0, draw ticks only c =1, draw grid lines c =2, draw grid lines to edge of frame. c nx -- number of intervals in x direction c sx -- spacing in inches between tick marks or grid lines c along the x axis c xs -- location of first tick or grid line on x axis c xf -- location of right edge of frame c ny -- number of intervals in y direction c sy -- spacing in inches between tick marks or grid lines c along the y axis c ys -- location of first tick or grid line on y axis c yf -- location of top edge of frame c assumptions: nx, sx, ny, sy all positive. c the lower left-hand corner of the frame is drawn at (0,0) c if xs<0, use 0; if ys<0, use 0 c if xf<=0, use nx*sx; if yf<=0, use ny*sy. c xinc = sx yinc = sy xlgth = float(nx)*sx ylgth = float(ny)*sy xmin = amax1(xs,0.0) ymin = amax1(ys,0.0) xmax = amax1(xf,xlgth+xmin) ymax = amax1(yf,ylgth+ymin) c c draw frame. c call plot(0.0, 0.0, 3) call plot(xmax, 0.0, 2) call plot(xmax, ymax, 2) call plot(0.0, ymax, 2) call plot(0.0, 0.0, 2) if (nopt.ne.0) go to 130 c c draw tick marks. c do 120 j=1,4 go to (10, 50, 20, 40), j 10 x2 = 0.0 if (xmin.ne.0.0) x2 = xmin - sx y2 = 0.0 go to 30 20 xinc = -sx x2 = xmin + xlgth + sx if (xmax.eq.xmin+xlgth) x2 = xmax y2 = ymax 30 y1 = y2 y2 = y2 + sign(0.125,xinc) n = nx if (abs(xmax-xmin-xlgth)+abs(xmin)) 70, 80, 70 40 yinc = -sy y2 = ymin + ylgth + sy if (ymax.eq.ymin+ylgth) y2 = ymax x2 = 0.0 go to 60 50 y2 = 0.0 if (ymin.ne.0.0) y2 = ymin - sy x2 = xmax 60 x1 = x2 n = ny x2 = x2 - sign(0.125,yinc) if (abs(ymax-ymin-ylgth)+abs(ymin)) 70, 80, 70 70 n = n + 1 80 do 110 i=1,n if (mod(j,2).eq.0) go to 90 x2 = x2 + xinc x1 = x2 go to 100 90 y2 = y2 + yinc y1 = y2 100 call plot(x1, y1, 3) call plot(x2, y2, 2) 110 continue 120 continue go to 240 c c draw grid lines c 130 x1 = xmin x2 = xmin + xlgth if (nopt.ne.2) go to 140 x1 = 0.0 x2 = xmax 140 y1 = ymin - sy n = ny + 1 if (ymax.eq.ymin+ylgth) n = n - 1 if (ymin.ne.0.0) go to 150 y1 = 0.0 n = n - 1 150 if (n.le.0) go to 170 j = 1 do 160 i=1,n j = -j y1 = y1 + sy call plot(x1, y1, 3) call plot(x2, y1, 2) xx = x1 x1 = x2 x2 = xx 160 continue 170 y1 = ymin + ylgth y2 = ymin if (nopt.ne.2) go to 180 y1 = ymax y2 = 0.0 180 n = nx + 1 if (j.lt.0) go to 200 x1 = xmin - sx if (xmax.eq.xmin+xlgth) n = n - 1 if (xmin.ne.0.0) go to 190 x1 = 0.0 n = n - 1 190 if (n.le.0) go to 240 xinc = sx go to 220 200 x1 = xmin + xlgth + sx if (xmin.eq.0.0) n = n - 1 if (xmax.ne.xlgth+xmin) go to 210 n = n - 1 x1 = xmax 210 xinc = -sx 220 do 230 i=1,n x1 = x1 + xinc call plot(x1, y1, 3) call plot(x1, y2, 2) xx = y1 y1 = y2 y2 = xx 230 continue 240 return c end --------------------------- dynamics@janus.Berkeley.EDU I have no idea what method is best, but the June 1987 issue of BYTE had a clever algorithm that I used. The basic idea is pretty simple, and I'll try to explain it. They gave an implementation in basic as I recall. Consider your typical grid square (it doesn't really need to be a square): 1 2 Each of these has a height attached to it. 3 4 Calculate the height at the middle of the square (at the 5 below) as just the average of the other four heights: 1 2 5 3 4 Now you have a set of four triangles in each square. For example, one is bounded by the lines connecting the heights at points 3, 4 and 5. The contours are then the intersections of the contour planes (i.e. z=z_0) with these triangles. So you just find all these line segments and plot them. Hope this helps. Actually, if you find out about a better algorithm, I'd be interested in hearing about it. chris goedde dynamics@janus.berkeley.edu --------------------------- Rich@rice.edu Hi Ray, I located the code but time won't permit me to look at it. I just grabbed what was there without even trying to compile it. I put the code on qed.rice.edu (128.42.4.38) in pub/contour.tar.Z. Most of it is in fortran. If you can help clean it up a little please do. I've got several other peices I'd like to release now, so I can't spend the time required to clean it up. I have supported plot3d which is derived from the same sources. If contour doesn't work right away, take a look at plot3d.tar.Z in the same directory. I'm sure it will be easier to compile. I've worked on the C front end to the fortran code for plot3d, and the C front end for the fortran contour plot code is similar. If you start working on it I can help. I'm looking for testers for a simple 3d equivalent of unix graph. Yell if your'e interested. You need unix and g++ to compile it. Rich --------------------------- mccalpin@perelandra.cms.udel.edu "John D. McCalpin" I suspect that the easiest approach is that used in the NCAR (National Center for Atmospheric Research) CONREC package. For randomly spaced data, one first interpolates to a fixed grid. Then you simply walk through the grid and compare the current point with its neighbors. If the contour level falls in between the values, simply linearly interpolate to find the location in space of the contour and follow that contour around until it either leaves the domain or closes on itself. --------------------------- --------------------------- -------------------------------------------- froncio@caip.rutgers.edu (Andy Froncioni) Well, it's not the most sophisticated way, but you could try o triangulating the points (*) and defining a linear (planar) interpolation within each triangular cell. The interpolation must be continuous across cell boundaries, but that's not difficult to do. I'll provide you with a simple one (**). o query each cell about its surface intersection with the given contour value. The result is either <nil> or a segment. Each cell of interest therefore coughs up a segment. (***) o plot each segment (in any order) and you'll have a contour line, since the segments are guaranteed to match up at cell boundaries. ----------- (*) Note 1: Make sure that the trianglulation orientation is counterclockwise. That is, you are to contruct ordered triplets of points such that the triplets define triangles with the vertices in counterclockwise order. If you have a convex shape, you can use voronoi triangulation. You can get a voronoi triangulation program from netlib (send a one-line e-mail message to netlib@research.att.com with the query 'send sweep2 from voronoi') which will perform this step automatically for you. ----------- (**) Note 2: An interpolation you can use within a cell whose vertex coordinates are ( (x1,y1) ; (x2,y2) ; (x3,y3) ) and with corresponding surface values (f1,f2,f3) is: f(x,y) = f1*z1(x,y) + f2*z2(x,y) + f3*z3(x,y) where: z1(x,y) = ( (x2*y3-x3*y2) + x*(y2-y3) + y*(x3-x2) ) / (2*A) z2(x,y) = ( (y1*x3-y3*x1) + x*(y3-y1) + y*(x1-x3) ) / (2*A) z3(x,y) = ( (x1*y2-x2*y1) + x*(y1-y2) + y*(x2-x1) ) / (2*A) and (2*A) is the twice the area of the triangle: (2*A) = (x2*y3-x3*y2) - x1*(y3-y2) + y1*(x2-x2) ------------- (***) Note 3: Actually, for the linear interpolation, it is easiest to just query each triangle edge for an intersection point. A cell which gets "cut" by the contour value will generate 2 such cut-points. This determines your inter- section segment in the cell. In such a case, the interpolation need only be defined linearly along each unique segment. ------------- >References would help, but I would like to know which method is best and/or >easiest to implement. In 3d, there is an algorithm called 'Marching Cubes' which effectively does the same thing as I've described: Lorensen, W. and Cline, H. "Marching Cubes: A High-Resolution 3D Surface Construction Algorithm." Computer Graphics, 21(4): 163-170, August, 1987. For details of the 2D stuff, above, I'd be glad to clarify anything that I've managed to confuse you about... Andy Andy Froncioni froncio@caip.rutgers.edu -------------------------------------------- ahiggins@pequod.cso.uiuc.edu (Andrew Higgins) > >Well, it's not the most sophisticated way, but you could try > >Lorensen, W. and Cline, H. "Marching Cubes: A High-Resolution 3D Surface >Construction Algorithm." Computer Graphics, 21(4): 163-170, August, 1987. >Andy > >Andy Froncioni >froncio@caip.rutgers.edu I'd also recommend: Lawson,C.L. "Generation of a Triangular Grid with Application to Contour Plotting." Tech. Memo. 299, Sect. 914, Jet Propulsion Lab., Calif. Instit. Tech., Pasadena, Calif., Feb. 1972. Akima, Hiroshi "A Method of Bivariate Interpolation and Smooth Surface Fitting for Irregularly Distributed Data Points." ACM Transactions on Mathematical Software, Vol. 4, No. 2, June 1978, Pages 148-159. This combined Lawson (for triangulation) and Akima (for interpolation) method forms the basis of both NCAR and PVI DI-3000's contouring package for irregular or randomly spaced data. Be forwarned, if you have lots (a very precise term) of data, this method is extremely S-L-O-W. -- Andrew J. Higgins -------------------------------------------- wahba@stat.wisc.edu (Grace Wahba) Getting Better Contour Plots with S and GCVPACK by Douglas M. Bates, Fred Reames and Grace Wahba University of Wisconsin-Madison Statistics Dept TR865 April, 1990 We show how to obtain aesthetically pleasing contour plots with thin plate splines by using GCVPACK (available thru netlib) and New S or any similar code that plots contours from points on a grid. Available while supplies last by sending your mailing address to gao@stat.wisc.edu with S-GCVPACK in the subject. hanche@imf.unit.no (Harald Hanche-Olsen) You may want to take the probabilistic, or stocahstic process, approach and look up a technique called Kriging. It has been used in (inter/extra)polating values from measurements at discrete points in geophysical applications. Sorry, but I have no references available. The keyword Kriging might get you somewhere at your library, though...
pdbourke@ccu1.aukuni.ac.nz (Bourke) (03/29/91)
chrisp@bilby.cs.uwa.oz.au (Chris Pudney) writes: >G'day, >I require an algorithm for plotting the contour lines on a surface. The >surfaces are basically 2-D potential fields (f:R^2 -> R) which I would like >to depict graphically using contour lines (or perhaps I should use >greyscale/colour?) >Could you provide me with pointers to a contour plotting algorithm. Public >domain (ftp'able) source would be most desirable. If I get a large enough >response I will post a compilation of the information I receive. A few years ago I wrote an article for BYTE magazine on just such an algorithm, in fact from memory the example I used was of a potential field. I don't have the article here (it's at home somewhere along with the source code) but it appeared I think in the June/July issue in either 1987 or 1988. It was titled CONREC - A Contouring Algorithm. Briefly - it contoured any rectangular grid of points. To contour a surface described by a function one would first evaluate the function on a rectangular mesh (not necessarily a regular spaced mesh) and pass that data onto the contouring program. One of the main features of the program was that it was very efficient, ie: even version in BASIC were quite fast. Have a look in past issues of BYTE, if you can't find it I can send the code in BASIC or FORTRAN (maybe C).