home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / pc / source / contour.f < prev    next >
Encoding:
Text File  |  1992-04-19  |  9.4 KB  |  251 lines
  1.  
  2. The program ACORD is taken, with a slight revision, from the article 
  3. "ACORD: AUTOMATIC COUNTOURING OF RAW DATA, Computers & Geosciences, 
  4. vol. 8, no. 1, p. 97-101, 1982", by D.F. Watson.  It was written in 
  5. 1979 in an ancient version of FORTRAN to run on the Control Data 
  6. Cyber main-frame.  It accepts a set of bivariate data {x, y, f(x,y)} 
  7. and finds the Delaunay triangulation.  Then the trace of each contour 
  8. level on each triangle is output, that is, piecewise linear triangle-
  9. based interpolation.  This program is notable only for the concise 
  10. triangulation algorithm, which is better than O(N**2).  The arrays 
  11. will take a maximum of 1200 data plus allowing for three control 
  12. vertices.  The data set that follows, included for reference, is 
  13. from "Davis, J.C., 1986, Statistics and data analysis in geology, 
  14. 2ed., Wiley, p.362" and is prefaced by the number of data, the number 
  15. of contours, the input format, and the contour levels, in that order.  
  16. Most data sets will require perturbation to avoid degenerate 
  17. configuations so a small random value is randomly added or subtracted 
  18. to or from each x and y coordinate.
  19. Discussion or (gentle) criticsm is invited.
  20.  
  21. Dave Watson                   Internet: watson@maths.uwa.oz.au
  22. Department of Mathematics            
  23. The University of Western Australia       Tel: (61 9) 380 3359
  24. Nedlands, WA 6009  Australia.             FAX: (61 9) 380 1028
  25. --------------------------- snipity-snip -------------------------------
  26.       PROGRAM ACORD
  27. C.....ACORD: Automatic Contouring Of Raw Data - by D.F.Watson
  28. C.....READ DATA POINTS AND FORM ALL 3-TUPLES S.T. NO OTHER POINT
  29. C     LIES WITHIN THAT 3-TUPLE'S CIRCUMCIRCLE
  30. C.....PNT HOLDS THE DATA POINTS AND TETR CARRIES THE CIRCUMCIRCLE
  31. C     CENTRE AND RADIUS SQUARED FOR EACH 3-TUPLE
  32. C.....ITETR HOLDS THE DATA POINT INDICES IN INPUT ORDER OF EACH 3-TUPLE
  33. C.....ISTACK IS A LAST-IN-FIRST-OUT STACK OF INDICES OF VACANT 3-TUPLES
  34. C.....KTETR IS A TEMPORARY LIST OF EDGES OF DELETED 3-TUPLES
  35. C.....ID IS POINTER TO ISTACK, AND JT IS POINTER TO TETR AND ITETR
  36.       REAL PNT(1203,3), TETR(2401,3), XPNT(3,3), CONT(50), DET(2,3)
  37.       INTEGER ITETR(2401,3), ISTACK(2401), KTETR(50,2), ITEMP(3,2)
  38.       CHARACTER*32 FORM
  39.       DATA ITEMP, XPNT/1,1,2,2,3,3,-1.,5.,-1.,-1.,-1.,5.,2.,2.,18./
  40.       DATA NTETR,ID,XMIN,YMIN,XMAX,YMAX/1,2,1.E37,1.E37,-1.E37,-1.E37/
  41. C......................................................................|
  42. C.....OPEN FILE FOR DATA INPUT
  43.       OPEN(UNIT=7, FILE="jdavis.dat", ERR=99)
  44. C.....READ NO. OF DATA POINTS, NO. OF CONTOURS, DATA FORMAT
  45.       READ(7,200) NDATA, LC, FORM
  46.   200 FORMAT(I8,I5,A13)
  47. C.....READ THE CONTOUR VALUES
  48.       READ(7,FORM) (CONT(I),I=1,LC)
  49.       RNUM = RAND(1)
  50. C.....READ THE DATA
  51.       DO 5 I = 1,NDATA
  52.         READ(7,FORM) PNT(I,1), PNT(I,2), PNT(I,3)
  53.         RNUM = RAND(RNUM)
  54.         PNT(I,1) = PNT(I,1) + 0.00001 * (RNUM - 0.5)
  55.         RNUM = RAND(RNUM)
  56.         PNT(I,2) = PNT(I,2) + 0.00001 * (RNUM - 0.5)
  57.         IF(PNT(I,1) .GT. XMAX) XMAX = PNT(I,1)
  58.         IF(PNT(I,1) .LT. XMIN) XMIN = PNT(I,1)
  59.         IF(PNT(I,2) .GT. YMAX) YMAX = PNT(I,2)
  60.     5   IF(PNT(I,2) .LT. YMIN) YMIN = PNT(I,2)
  61.       CLOSE(7)
  62.       DATAX = XMAX - XMIN
  63.       Y = YMAX - YMIN
  64.       IF(Y .GT. DATAX) DATAX = Y
  65. C.....NORMALIZE THE DATA
  66.       DO 7 I = 1,NDATA
  67.         PNT(I,1) = (PNT(I,1) - XMIN)/DATAX
  68.     7   PNT(I,2) = (PNT(I,2) - YMIN)/DATAX
  69.       DO 8 I = 1,3
  70.         ITETR(1,I) = I + NDATA
  71.         TETR(1,I) = XPNT(I,3)
  72.         DO 8 J = 1,2
  73.     8     PNT(I + NDATA,J) = XPNT(I,J)
  74.       J = NDATA * 2 + 1
  75.       DO 9 I = 2,J
  76.     9   ISTACK(I) = I
  77. C......................................................................|
  78.       DO 50 NUC = 1,NDATA
  79.         KM = 0
  80. C.......LOOP THRU THE ESTABLISHED 3-TUPLES
  81.         DO 30 JT = 1,NTETR
  82. C.........TEST IF NEW DATA POINT IS WITHIN THE JT CIRCUMCIRCLE
  83.           DSQ = TETR(JT,3) - (PNT(NUC,1) - TETR(JT,1))**2
  84.           IF(DSQ .LT. 0.0) GO TO 30
  85.           DSQ = DSQ - (PNT(NUC,2) - TETR(JT,2))**2
  86.           IF(DSQ .LT. 0.0) GO TO 30
  87. C.........DELETE THIS 3-TUPLE BUT SAVE ITS EDGES
  88.           ID = ID - 1
  89.           ISTACK(ID) = JT
  90. C.........ADD EDGES TO KTETR BUT DELETE IF ALREADY LISTED
  91.           DO 28 I = 1,3
  92.             L1 = ITEMP(I,1)
  93.             L2 = ITEMP(I,2)
  94.             IF(KM .LE. 0) GO TO 26
  95.             KMT = KM
  96.             DO 24 J = 1,KMT
  97.               IF(ITETR(JT,L1) .NE. KTETR(J,1)) GO TO 24
  98.               IF(ITETR(JT,L2) .NE. KTETR(J,2)) GO TO 24
  99.               KM = KM - 1
  100.               IF(J .GT. KM) GO TO 28
  101.               DO 20 K = J,KM
  102.                 K1 = K + 1
  103.                 DO 20 L = 1,2
  104.    20             KTETR(K,L) = KTETR(K1,L)
  105.                 GO TO 28
  106.    24         CONTINUE
  107.    26       KM = KM + 1
  108.             KTETR(KM,1) = ITETR(JT,L1)
  109.             KTETR(KM,2) = ITETR(JT,L2)
  110.    28       CONTINUE
  111.    30     CONTINUE
  112. C
  113. C.......FORM NEW 3-TUPLE3
  114.         DO 48 I = 1,KM
  115.           KT = ISTACK(ID)
  116.           ID = ID + 1
  117. C.........CALCULATE THE CIRCUMCIRCLE CENTER AND RADIUS
  118. C         SQUARED OF POINTS KTETR(I,*) AND PLACE IN TETR(KT,*)
  119.           DO 44 JZ = 1,2
  120.             I2 = KTETR(I,JZ)
  121.             DET(JZ,1) = PNT(I2,1) - PNT(NUC,1)
  122.             DET(JZ,2) = PNT(I2,2) - PNT(NUC,2)
  123.    44       DET(JZ,3) = DET(JZ,1) * (PNT(I2,1) + PNT(NUC,1))/2.0
  124.      +                + DET(JZ,2) * (PNT(I2,2) + PNT(NUC,2))/2.0
  125.           DD = DET(1,1) * DET(2,2) - DET(1,2) * DET(2,1)
  126.           TETR(KT,1) = (DET(1,3) * DET(2,2) - DET(2,3) * DET(1,2))/DD
  127.           TETR(KT,2) = (DET(1,1) * DET(2,3) - DET(2,1) * DET(1,3))/DD
  128.           TETR(KT,3) = (PNT(NUC,1) - TETR(KT,1))**2
  129.      +               + (PNT(NUC,2) - TETR(KT,2))**2
  130.           ITETR(KT,1) = KTETR(I,1)
  131.           ITETR(KT,2) = KTETR(I,2)
  132.    48     ITETR(KT,3) = NUC
  133.    50   NTETR = NTETR + 2
  134. C......................................................................|
  135. C.....SET OUTPUT SIZE FACTOR
  136.       SIZE = 1.0
  137.       OPEN(UNIT=8, FILE="contsegm.dat", ERR=99)
  138. C
  139.       DO 90 JT = 1,NTETR
  140.         IF(ITETR(JT,1) .GT. NDATA .OR. TETR(JT,3) .GT. 1) GO TO 90
  141. C.......FIND CONTOUR INTERSECTIONS
  142.         ABIT = 0.0
  143.         IF(PNT(ITETR(JT,1),3) .EQ. PNT(ITETR(JT,2),3) .OR.
  144.      +     PNT(ITETR(JT,1),3) .EQ. PNT(ITETR(JT,3),3) .OR.
  145.      +     PNT(ITETR(JT,2),3) .EQ. PNT(ITETR(JT,3),3)) ABIT = 1.E-10
  146.         TOP = AMAX1(PNT(ITETR(JT,1),3),PNT(ITETR(JT,2),3),
  147.      +              PNT(ITETR(JT,3),3))
  148.         BOT = AMIN1(PNT(ITETR(JT,1),3),PNT(ITETR(JT,2),3),
  149.      +              PNT(ITETR(JT,3),3))
  150.         DO 87 JC = 1,LC
  151.           IF(CONT(JC) .GT. TOP .OR. CONT(JC) .LT. BOT) GO TO 87
  152.           CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,2),3)
  153.      +                   - PNT(ITETR(JT,1),3) + ABIT)
  154.           IF(CZ .LE. 0.0 .OR. CZ .GE. 1.0) GO TO 83
  155.           X1 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,2),1)
  156.      +       -  PNT(ITETR(JT,1),1)) * CZ) * SIZE
  157.           Y1 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,2),2)
  158.      +       -  PNT(ITETR(JT,1),2)) * CZ) * SIZE
  159.    82     CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,3),3)
  160.      +                   - PNT(ITETR(JT,1),3) + ABIT)
  161.           IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 84
  162.           X2 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,3),1)
  163.      +       -  PNT(ITETR(JT,1),1)) * CZ) * SIZE
  164.           Y2 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,3),2)
  165.      +       -  PNT(ITETR(JT,1),2)) * CZ) * SIZE
  166.           GO TO 85
  167.    83     CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,3),3)
  168.      +                   - PNT(ITETR(JT,1),3) + ABIT)
  169.           IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 87
  170.           X1 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,3),1)
  171.      +       -  PNT(ITETR(JT,1),1)) * CZ) * SIZE
  172.           Y1 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,3),2)
  173.      +       -  PNT(ITETR(JT,1),2)) * CZ) * SIZE
  174.    84     CZ = (CONT(JC) - PNT(ITETR(JT,2),3))/(PNT(ITETR(JT,3),3)
  175.      +                   - PNT(ITETR(JT,2),3) + ABIT)
  176.           IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 87
  177.           X2 = (PNT(ITETR(JT,2),1) + (PNT(ITETR(JT,3),1)
  178.      +       -  PNT(ITETR(JT,2),1)) * CZ) * SIZE
  179.           Y2 = (PNT(ITETR(JT,2),2) + (PNT(ITETR(JT,3),2)
  180.      +       -  PNT(ITETR(JT,2),2)) * CZ) * SIZE
  181.    85     CONTINUE
  182. C         WRITE A CONTOUR SEGMENT TO FILE
  183. C         "3" means 'pen-up' and "2" means 'pen-down"
  184.           WRITE(8,400) X1, Y1, "3", X2, Y2, "2"
  185.   400     FORMAT(2F10.4,A5,2F10.4,A5)
  186.    87     CONTINUE
  187.    90   CONTINUE
  188.       CLOSE(8)
  189.    99 STOP
  190.       END
  191. ------------------------------- snipity-snip  ------------------------------
  192.       52   11 (3F5.0)
  193.      700  725  750
  194.      775  800  825
  195.      850  875  900
  196.      925  950
  197.      0.3  6.1  870
  198.      1.4  6.2  793
  199.      2.4  6.1  755
  200.      3.6  6.2  690
  201.      5.7  6.2  800
  202.      1.6  5.2  800
  203.      2.9  5.1  730
  204.      3.4  5.3  728
  205.      3.4  5.7  710
  206.      4.8  5.6  780
  207.      5.3  5.0  804
  208.      6.2  5.2  855
  209.      0.2  4.3  830
  210.      0.9  4.2  813
  211.      2.3  4.8  762
  212.      2.5  4.5  765
  213.      3.0  4.5  740
  214.      3.5  4.5  765
  215.      4.1  4.6  760
  216.      4.9  4.2  790
  217.      6.3  4.3  820
  218.      0.9  3.2  855
  219.      1.7  3.8  812
  220.      2.4  3.8  773
  221.      3.7  3.5  812
  222.      4.5  3.2  827
  223.      5.2  3.2  805
  224.      6.3  3.4  840
  225.      0.3  2.4  890
  226.      2.0  2.7  820
  227.      3.8  2.3  873
  228.      6.3  2.2  875
  229.      0.6  1.7  873
  230.      1.5  1.8  865
  231.      2.1  1.8  841
  232.      2.1  1.1  862
  233.      3.1  1.1  908
  234.      4.5  1.8  855
  235.      5.5  1.7  850
  236.      5.7  1.0  882
  237.      6.2  1.0  910
  238.      0.4  0.5  940
  239.      1.4  0.6  915
  240.      1.4  0.1  890
  241.      2.1  0.7  880
  242.      2.3  0.3  870
  243.      3.1  0.0  880
  244.      4.1  0.8  960
  245.      5.4  0.4  890
  246.      6.0  0.1  860
  247.      5.7  3.0  830
  248.      3.6  6.0  705
  249. -------------------------------------------------
  250.  
  251.