marwk@levels.sait.edu.au (12/12/90)
Repost of summary of contouring methods - I tried replying to keith@rses...
but it is not recognised at this site.
A summary of the relevant responses appears below; I will evaluate the
methods as time permits.
If anyone else evaluates these responses could they please let us all have the
benefit of their evaluations?
I thank all for their responses.
Ray
---------------------------
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...