The program ACORD is taken, with a slight revision, from the article "ACORD: AUTOMATIC COUNTOURING OF RAW DATA, Computers & Geosciences, vol. 8, no. 1, p. 97-101, 1982", by D.F. Watson. It was written in 1979 in an ancient version of FORTRAN to run on the Control Data Cyber main-frame. It accepts a set of bivariate data {x, y, f(x,y)} and finds the Delaunay triangulation. Then the trace of each contour level on each triangle is output, that is, piecewise linear triangle- based interpolation. This program is notable only for the concise triangulation algorithm, which is better than O(N**2). The arrays will take a maximum of 1200 data plus allowing for three control vertices. The data set that follows, included for reference, is from "Davis, J.C., 1986, Statistics and data analysis in geology, 2ed., Wiley, p.362" and is prefaced by the number of data, the number of contours, the input format, and the contour levels, in that order. Most data sets will require perturbation to avoid degenerate configuations so a small random value is randomly added or subtracted to or from each x and y coordinate. Discussion or (gentle) criticsm is invited. Dave Watson Internet: watson@maths.uwa.oz.au Department of Mathematics The University of Western Australia Tel: (61 9) 380 3359 Nedlands, WA 6009 Australia. FAX: (61 9) 380 1028 --------------------------- snipity-snip ------------------------------- PROGRAM ACORD C.....ACORD: Automatic Contouring Of Raw Data - by D.F.Watson C.....READ DATA POINTS AND FORM ALL 3-TUPLES S.T. NO OTHER POINT C LIES WITHIN THAT 3-TUPLE'S CIRCUMCIRCLE C.....PNT HOLDS THE DATA POINTS AND TETR CARRIES THE CIRCUMCIRCLE C CENTRE AND RADIUS SQUARED FOR EACH 3-TUPLE C.....ITETR HOLDS THE DATA POINT INDICES IN INPUT ORDER OF EACH 3-TUPLE C.....ISTACK IS A LAST-IN-FIRST-OUT STACK OF INDICES OF VACANT 3-TUPLES C.....KTETR IS A TEMPORARY LIST OF EDGES OF DELETED 3-TUPLES C.....ID IS POINTER TO ISTACK, AND JT IS POINTER TO TETR AND ITETR REAL PNT(1203,3), TETR(2401,3), XPNT(3,3), CONT(50), DET(2,3) INTEGER ITETR(2401,3), ISTACK(2401), KTETR(50,2), ITEMP(3,2) CHARACTER*32 FORM DATA ITEMP, XPNT/1,1,2,2,3,3,-1.,5.,-1.,-1.,-1.,5.,2.,2.,18./ DATA NTETR,ID,XMIN,YMIN,XMAX,YMAX/1,2,1.E37,1.E37,-1.E37,-1.E37/ C......................................................................| C.....OPEN FILE FOR DATA INPUT OPEN(UNIT=7, FILE="jdavis.dat", ERR=99) C.....READ NO. OF DATA POINTS, NO. OF CONTOURS, DATA FORMAT READ(7,200) NDATA, LC, FORM 200 FORMAT(I8,I5,A13) C.....READ THE CONTOUR VALUES READ(7,FORM) (CONT(I),I=1,LC) RNUM = RAND(1) C.....READ THE DATA DO 5 I = 1,NDATA READ(7,FORM) PNT(I,1), PNT(I,2), PNT(I,3) RNUM = RAND(RNUM) PNT(I,1) = PNT(I,1) + 0.00001 * (RNUM - 0.5) RNUM = RAND(RNUM) PNT(I,2) = PNT(I,2) + 0.00001 * (RNUM - 0.5) IF(PNT(I,1) .GT. XMAX) XMAX = PNT(I,1) IF(PNT(I,1) .LT. XMIN) XMIN = PNT(I,1) IF(PNT(I,2) .GT. YMAX) YMAX = PNT(I,2) 5 IF(PNT(I,2) .LT. YMIN) YMIN = PNT(I,2) CLOSE(7) DATAX = XMAX - XMIN Y = YMAX - YMIN IF(Y .GT. DATAX) DATAX = Y C.....NORMALIZE THE DATA DO 7 I = 1,NDATA PNT(I,1) = (PNT(I,1) - XMIN)/DATAX 7 PNT(I,2) = (PNT(I,2) - YMIN)/DATAX DO 8 I = 1,3 ITETR(1,I) = I + NDATA TETR(1,I) = XPNT(I,3) DO 8 J = 1,2 8 PNT(I + NDATA,J) = XPNT(I,J) J = NDATA * 2 + 1 DO 9 I = 2,J 9 ISTACK(I) = I C......................................................................| DO 50 NUC = 1,NDATA KM = 0 C.......LOOP THRU THE ESTABLISHED 3-TUPLES DO 30 JT = 1,NTETR C.........TEST IF NEW DATA POINT IS WITHIN THE JT CIRCUMCIRCLE DSQ = TETR(JT,3) - (PNT(NUC,1) - TETR(JT,1))**2 IF(DSQ .LT. 0.0) GO TO 30 DSQ = DSQ - (PNT(NUC,2) - TETR(JT,2))**2 IF(DSQ .LT. 0.0) GO TO 30 C.........DELETE THIS 3-TUPLE BUT SAVE ITS EDGES ID = ID - 1 ISTACK(ID) = JT C.........ADD EDGES TO KTETR BUT DELETE IF ALREADY LISTED DO 28 I = 1,3 L1 = ITEMP(I,1) L2 = ITEMP(I,2) IF(KM .LE. 0) GO TO 26 KMT = KM DO 24 J = 1,KMT IF(ITETR(JT,L1) .NE. KTETR(J,1)) GO TO 24 IF(ITETR(JT,L2) .NE. KTETR(J,2)) GO TO 24 KM = KM - 1 IF(J .GT. KM) GO TO 28 DO 20 K = J,KM K1 = K + 1 DO 20 L = 1,2 20 KTETR(K,L) = KTETR(K1,L) GO TO 28 24 CONTINUE 26 KM = KM + 1 KTETR(KM,1) = ITETR(JT,L1) KTETR(KM,2) = ITETR(JT,L2) 28 CONTINUE 30 CONTINUE C C.......FORM NEW 3-TUPLE3 DO 48 I = 1,KM KT = ISTACK(ID) ID = ID + 1 C.........CALCULATE THE CIRCUMCIRCLE CENTER AND RADIUS C SQUARED OF POINTS KTETR(I,*) AND PLACE IN TETR(KT,*) DO 44 JZ = 1,2 I2 = KTETR(I,JZ) DET(JZ,1) = PNT(I2,1) - PNT(NUC,1) DET(JZ,2) = PNT(I2,2) - PNT(NUC,2) 44 DET(JZ,3) = DET(JZ,1) * (PNT(I2,1) + PNT(NUC,1))/2.0 + + DET(JZ,2) * (PNT(I2,2) + PNT(NUC,2))/2.0 DD = DET(1,1) * DET(2,2) - DET(1,2) * DET(2,1) TETR(KT,1) = (DET(1,3) * DET(2,2) - DET(2,3) * DET(1,2))/DD TETR(KT,2) = (DET(1,1) * DET(2,3) - DET(2,1) * DET(1,3))/DD TETR(KT,3) = (PNT(NUC,1) - TETR(KT,1))**2 + + (PNT(NUC,2) - TETR(KT,2))**2 ITETR(KT,1) = KTETR(I,1) ITETR(KT,2) = KTETR(I,2) 48 ITETR(KT,3) = NUC 50 NTETR = NTETR + 2 C......................................................................| C.....SET OUTPUT SIZE FACTOR SIZE = 1.0 OPEN(UNIT=8, FILE="contsegm.dat", ERR=99) C DO 90 JT = 1,NTETR IF(ITETR(JT,1) .GT. NDATA .OR. TETR(JT,3) .GT. 1) GO TO 90 C.......FIND CONTOUR INTERSECTIONS ABIT = 0.0 IF(PNT(ITETR(JT,1),3) .EQ. PNT(ITETR(JT,2),3) .OR. + PNT(ITETR(JT,1),3) .EQ. PNT(ITETR(JT,3),3) .OR. + PNT(ITETR(JT,2),3) .EQ. PNT(ITETR(JT,3),3)) ABIT = 1.E-10 TOP = AMAX1(PNT(ITETR(JT,1),3),PNT(ITETR(JT,2),3), + PNT(ITETR(JT,3),3)) BOT = AMIN1(PNT(ITETR(JT,1),3),PNT(ITETR(JT,2),3), + PNT(ITETR(JT,3),3)) DO 87 JC = 1,LC IF(CONT(JC) .GT. TOP .OR. CONT(JC) .LT. BOT) GO TO 87 CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,2),3) + - PNT(ITETR(JT,1),3) + ABIT) IF(CZ .LE. 0.0 .OR. CZ .GE. 1.0) GO TO 83 X1 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,2),1) + - PNT(ITETR(JT,1),1)) * CZ) * SIZE Y1 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,2),2) + - PNT(ITETR(JT,1),2)) * CZ) * SIZE 82 CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,3),3) + - PNT(ITETR(JT,1),3) + ABIT) IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 84 X2 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,3),1) + - PNT(ITETR(JT,1),1)) * CZ) * SIZE Y2 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,3),2) + - PNT(ITETR(JT,1),2)) * CZ) * SIZE GO TO 85 83 CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,3),3) + - PNT(ITETR(JT,1),3) + ABIT) IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 87 X1 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,3),1) + - PNT(ITETR(JT,1),1)) * CZ) * SIZE Y1 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,3),2) + - PNT(ITETR(JT,1),2)) * CZ) * SIZE 84 CZ = (CONT(JC) - PNT(ITETR(JT,2),3))/(PNT(ITETR(JT,3),3) + - PNT(ITETR(JT,2),3) + ABIT) IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 87 X2 = (PNT(ITETR(JT,2),1) + (PNT(ITETR(JT,3),1) + - PNT(ITETR(JT,2),1)) * CZ) * SIZE Y2 = (PNT(ITETR(JT,2),2) + (PNT(ITETR(JT,3),2) + - PNT(ITETR(JT,2),2)) * CZ) * SIZE 85 CONTINUE C WRITE A CONTOUR SEGMENT TO FILE C "3" means 'pen-up' and "2" means 'pen-down" WRITE(8,400) X1, Y1, "3", X2, Y2, "2" 400 FORMAT(2F10.4,A5,2F10.4,A5) 87 CONTINUE 90 CONTINUE CLOSE(8) 99 STOP END ------------------------------- snipity-snip ------------------------------ 52 11 (3F5.0) 700 725 750 775 800 825 850 875 900 925 950 0.3 6.1 870 1.4 6.2 793 2.4 6.1 755 3.6 6.2 690 5.7 6.2 800 1.6 5.2 800 2.9 5.1 730 3.4 5.3 728 3.4 5.7 710 4.8 5.6 780 5.3 5.0 804 6.2 5.2 855 0.2 4.3 830 0.9 4.2 813 2.3 4.8 762 2.5 4.5 765 3.0 4.5 740 3.5 4.5 765 4.1 4.6 760 4.9 4.2 790 6.3 4.3 820 0.9 3.2 855 1.7 3.8 812 2.4 3.8 773 3.7 3.5 812 4.5 3.2 827 5.2 3.2 805 6.3 3.4 840 0.3 2.4 890 2.0 2.7 820 3.8 2.3 873 6.3 2.2 875 0.6 1.7 873 1.5 1.8 865 2.1 1.8 841 2.1 1.1 862 3.1 1.1 908 4.5 1.8 855 5.5 1.7 850 5.7 1.0 882 6.2 1.0 910 0.4 0.5 940 1.4 0.6 915 1.4 0.1 890 2.1 0.7 880 2.3 0.3 870 3.1 0.0 880 4.1 0.8 960 5.4 0.4 890 6.0 0.1 860 5.7 3.0 830 3.6 6.0 705 -------------------------------------------------