Text and code for a contouring algorithm, wrote it a long time ago but it is still useful. If you want the image see the June or July 1987 Byte article. CONREC - A Contouring Subroutine Written by Paul D. Bourke (MSc.) Introduction This article introduces a straightforward method of contouring some surface. The algorithm is implemented in a subroutine written in BASIC. Contouring aids in visualizing three dimensional surfaces on a two dimensional medium (on paper or in this case a computer graphics screen). Two most common applications are displaying topological features of an area on a map or the air pressure on a weather map. In all cases some parameter is plotted as a function of two variables, the longitude and latitude or x and y axis. One problem with computer contouring is the process is usually CPU intensive and the algorithms often using advanced mathematical techniques. CONREC To do contouring in software you need to describe the data surface and the contour levels you want to have drawn. The software given this information must call the algorithm that calculate sthe line segments that make up a contour curve and then plot these line segments on whatever graphics device is available. CONREC satisfies the above description, it is relatively simple to implement, very reliable, and does not require sophisticated programming techniques or a high level of mathematics to understand how it works. The input parameters to the CONREC subroutine are as follows : - The number of horizontal and vertical data points designated iub and jub. - The number of contouring levels, nc. - A one dimensional array z(0:nc-1) that seves as a list of the contour levels in increasing order. (The order of course can be relaxed if the program will sort the levels) - A two dimensional array d(0:iub,0:jub) that contains the description of the data array to be contoured. Each element of the array is a sample of the surface being studied at a point (x,y) - Two, one dimensional arrays x(0:iub) and y(0:jub) which contain the horizontal and vertical coordinates of each sample point. This allows for a rectangular grid of samples. Diagram 0 illustates some of the above input parameters. The contouring subroutine CONREC does not assume anything about the device that will be used to plot the contours. It instead expects a user written subroutine called VECOUT. CONREC calls VECOUT with the horizontal and vertical coordinates of the start and end coordinates of a line segment along with the contour level for that line segment. In the simplest case this is very similar the the usual LINE (x1,y1)-(x2,y2) command in BASIC. See examples. Algorithm As already mentioned the samples of the three dimensional surface are stored in a two dimensional real array. This rectangular grid is considered four points at a time, namely the rectangle d(i,j), d(i+1,j), d(i,j+1), and d(i+1,j+1). The centre of each rectangle is assigned a value corresponding to the average values of each of the four vertices. Each rectangle is in turn divided into four triangular regions by cutting along the diagonals. Each of these triangular planes may be bisected by horizontal contour plane. The intersection of these two planes is a straight line segment, part of the the contour curve at that contour height. Depending on the value of a contour level with respect to the height at the vertices of a triangle, certain types of contour lines are drawn. The 10 possible cases which may occur are summarised below (also see table 1) a) All the vertices lie below the contour level. b) Two vertices lie below and one on the contour level. c) Two vertices lie below and one above the contour level. d) One vertex lies below and two on the contour level. e) One vertex lies below, one on and one above the contour level. f) One vertex lies below and two above the contour level. g) Three vertices lie on the contour level. h) Two vertices lie on and one above the contour level. i) One vertex lies on and two above the contour level. j) All the vertices lie above the contour level. In cases a, b, i and j the two planes do not intersect, ie: no line need be drawn. For cases d and h the two planes intersect along an edge of the triangle and therefore line is drawn between the two vertices that lie on the contour level. Case e requires that a line be drawn from the vertex on the contour level to a point on the opposite edge. This point is determined by the intersection of the contour level with the straight line between the other two vertices. Cases c and f are the most common situations where the line is drawn from one edge to another edge of the triangle. The last possibility or case g above has no really satisfactory solution and fortunately will occur rarely with real arithmetic. Diagram 2 summarises the possible line orientations. Example As a simple example consider one triangle with vertices labelled m1,m2 and m3 with heights 0,2 and 3 respectively (See diagram 3) To calculate where a contour line at a height of 1 should be drawn, it can be seen that this is case f described above. Level 1 intersects line segment m1-m2 half the way along and it intersects line segment m1-m3 one third of the way along. A line segment is drawn between these two points. Subroutine In summary, CONREC takes each rectangle of adjacent data points and splits it into 4 triangles after chooseing the height at the centre of the rectangle. For each of the triangles the line segment resulting from the intersection with each contour plane. A routine is then called with the starting and stopping coordinates of the line segment. Listing 1 shows the CONREC subroutine written in MicroSoft Basic (version 2 binary) for the Macintosh. An attempt is made at optimization by checking first to see if there are any contour levels within the present rectangle and second that there are some contour levels within the present triangle. The indices i and j are used to step through each rectangle in turn, k refers to each contour level and m to the four triangles in each rectangle. Diagram 1 gives some of the notation used for identifying the rectangles and triangles in the subroutine. Note that for large arrays the whole data array need be stored in memory . Since the algorithm is a local one only requiring 4 points at a time, the data for each rectangle could be read from disk when required. Improvements You can add features such as labeling the contour and axes, drawing the contour lines with different colours or in the case of a black and white display, drawing the contour lines with different mark-space ratio's. The contours may be labeled with the contour values either by software or by adding them later by hand. Examples of CONREC Example 1 is a program that uses CONREC to draw the contour of a mathematical function. The function is given by 1 f ( x , y ) = sin ( ( x2 + y2 )1/2 ) + ---------------- ( ( x - c )2 + y2 )1/2 Example 2 is taken from physics. The potential a distance r from the origin from two charges q1 and q2, located at plus and minus c is given by q1 q2 V( x , y ) = const * ( --- - --- ) r1 r2 where r1 and r2 are the distances from the point ( x , y ) to the charge q1 and q2 respectively. Contouring this function gives rise to equipotential surfaces, Similar functions also arise in the case of gravitational potentials. Example 3 is a contour diagram of the function 1 f( x , y ) = ------------------------------------------------- ((x2+(y-0.842)(y+0.842))2 + (x(y-0.842) + x(y-0.842))2 Example 4 is a contour diagram of cowboy hat shaped function f( x , y ) = sin ( x2 + y2 ) The examples given here are written in Basic and are therefore rather slow. Example 1 (31x31 data array and requesting 10 contour levals) took about 5 minutes to draw. ****** Demo FORTRAN-77 program ******* program demo implicit none c c The follwoing is a simplistic application of the CONREC routine. c A mathematical function is evaluated over a regular grid of points c on a computer raster graphics screen. c c Paul D. Bourke c integer*4 pxmin,pymin,pxmax,pymax parameter (pxmin = 10,pymin = 10,pxmax = 460,pymax = 300) integer*4 nx,ny parameter (nx = 100,ny = 50) integer*4 nc parameter (nc = 10) c real*4 d(0:nx,0:ny),x(0:nx),y(0:ny),z(1:nc) real*4 x1,y1,x2,y2 real*4 zmax,zmin integer*4 i,j c c Create an artificial data surface and calculate the c surface bounds for choosing the contour levels. c zmin = 1.0e30 zmax = -1.0e30 do 200 i=0,nx do 100 j=0,ny d(i,j) = i * j zmin = min(zmin,d(i,j)) zmax = max(zmax,d(i,j)) 100 continue 200 continue c c Set coordinates in Y array suitable for c automatic plotting on the graphics screen c do 300 j=0,ny y(j) = pymax - j * (pymax - pymin) / float(ny) 300 continue c c Set coordinates in X array suitable for c automatic plotting on the graphics screen c do 400 i=0,nx x(i) = i * (pxmax - pxmin) / float(nx) + pxmin 400 continue c c Set a full contingent of contour levels c do 500 i=1,nc z(i) = i * (zmax - zmin) / (nc + 1) 500 continue c c Draw a border around the contour plot c x1 = pxmin y1 = pymin x2 = pxmax y2 = pymax call vecout(x1,y1,x1,y2,0.0) call vecout(x1,y2,x2,y2,0.0) call vecout(x2,y2,x2,y1,0.0) call vecout(x2,y1,x1,y1,0.0) c c Call the contouring routine c call conrec(d,0,nx,0,ny,x,y,nc,z) pause c end ****** FORTRAN-77 CODE ******* c c====================================================================== c c CONREC is a contouring subroutine for rectangularily spaced data. c c It emits calls to a line drawing subroutine supplied by the user c which draws a contour map corresponding to real*4data on a randomly c spaced rectangular grid. The coordinates emitted are in the same c units given in the x() and y() arrays. c c Any number of contour levels may be specified but they must be c in order of increasing value. c c subroutine conrec(d,ilb,iub,jlb,jub,x,y,nc,z) c real*4 d(ilb:iub,jlb:jub) ! matrix of data to contour c integer ilb,iub,jlb,jub ! index bounds of data matrix c real*4 x(ilb:iub) ! data matrix column coordinates c real*4 y(jlb,jub) ! data matrix row coordinates c integer nc ! number of contour levels c real*4 z(1:nc) ! contour levels in increasing order c subroutine conrec(d,ilb,iub,jlb,jub,x,y,nc,z) real*4 d(ilb:iub,jlb:jub) integer ilb,iub,jlb,jub real*4 x(ilb:iub) real*4 y(jlb:jub) integer nc real*4 z(1:nc) c c Local declarations c real*4 h(0:4) integer sh(0:4) real*4 xh(0:4),yh(0:4) integer im(1:4),jm(1:4) integer case integer castab(-1:1,-1:1,-1:1) integer p1,p2 c c Data c data im/0,1,1,0/ data jm/0,0,1,1/ data castab/0,0,9, 0,1,5, 7,4,8, 1 0,3,6, 2,3,2, 6,3,0, 2 8,4,7, 5,1,0, 9,0,0/ c c Use statement functions for the line intersections c xsect(p1,p2) = (h(p2)*xh(p1)-h(p1)*xh(p2))/(h(p2)-h(p1)) ysect(p1,p2) = (h(p2)*yh(p1)-h(p1)*yh(p2))/(h(p2)-h(p1)) c c Scan the arrays, top down, left to right within rows c 20 do 100 j=jub-1,jlb,-1 do 90 i=ilb,iub-1 dmin = min(d(i,j),d(i,j+1),d(i+1,j),d(i+1,j+1)) dmax = max(d(i,j),d(i,j+1),d(i+1,j),d(i+1,j+1)) if (dmax.ge.z(1) .and. dmin.le.z(nc)) then do 80 k=1,nc if (z(k).ge.dmin .and. z(k).le.dmax) then do 22 m=4,0,-1 if (m.gt.0) then h(m)=d(i+im(m),j+jm(m))-z(k) xh(m)=x(i+im(m)) yh(m)=y(j+jm(m)) else h(0)=0.25*(h(1)+h(2)+h(3)+h(4)) xh(0)=0.5*(x(i)+x(i+1)) yh(0)=0.5*(y(j)+y(j+1)) endif if (h(m).gt.0.0) then sh(m)=+1 else if (h(m).lt.0.0) then sh(m)=-1 else sh(m)=0 endif 22 continue c c Note: at this stage the relative heights of the corners and the c centre are in the h array, and the corresponding coordinates are c in the xh and yh arrays. The centre of the box is indexed by 0 c and the 4 corners by 1 to 4 as shown below. c Each triangle is then indexed by the parameter m, and the 3 c vertices of each triangle are indexed by parameters m1,m2,and m3. c It is assumed that the centre of the box is always vertex 2 though c this isimportant only when all 3 vertices lie exactly on the same c contour level, in which case only the side of the box is drawn. c c c vertex 4 +-------------------+ vertex 3 c | \ / | c | \ m-3 / | c | \ / | c | \ / | c | m=2 X m=2 | the centre is vertex 0 c | / \ | c | / \ | c | / m=1 \ | c | / \ | c vertex 1 +-------------------+ vertex 2 c c c c Scan each triangle in the box c do 60 m=1,4 m1=m m2=0 if (m.ne.4) then m3=m+1 else m3=1 endif case = castab(sh(m1),sh(m2),sh(m3)) if (case.ne.0) then goto (31,32,33,34,35,36,37,38,39),case c c Case 1 - Line between vertices 1 and 2 c 31 x1=xh(m1) y1=yh(m1) x2=xh(m2) y2=yh(m2) goto 40 c c Case 2 - Line between vertices 2 and 3 c 32 x1=xh(m2) y1=yh(m2) x2=xh(m3) y2=yh(m3) goto 40 c c Case 3 - Line between vertices 3 and 1 c 33 x1=xh(m3) y1=yh(m3) x2=xh(m1) y2=yh(m1) goto 40 c c Case 4 - Line between vertex 1 and side 2-3 c 34 x1=xh(m1) y1=yh(m1) x2=xsect(m2,m3) y2=ysect(m2,m3) goto 40 c c Case 5 - Line between vertex 2 and side 3-1 c 35 x1=xh(m2) y1=yh(m2) x2=xsect(m3,m1) y2=ysect(m3,m1) goto 40 c c Case 6 - Line between vertex 3 and side 1-2 c 36 x1=xh(m3) y1=yh(m3) x2=xsect(m1,m2) y2=ysect(m1,m2) goto 40 c c Case 7 - Line between sides 1-2 and 2-3 c 37 x1=xsect(m1,m2) y1=ysect(m1,m2) x2=xsect(m2,m3) y2=ysect(m2,m3) goto 40 c c Case 8 - Line between sides 2-3 and 3-1 c 38 x1=xsect(m2,m3) y1=ysect(m2,m3) x2=xsect(m3,m1) y2=ysect(m3,m1) goto 40 c c Case 9 - Line between sides 3-1 and 1-2 c 39 x1=xsect(m3,m1) y1=ysect(m3,m1) x2=xsect(m1,m2) y2=ysect(m1,m2) goto 40 40 call vecout(x1,y1,x2,y2,z(k)) endif 60 continue endif 80 continue endif 90 continue 100 continue return end c c====================================================================== c c This is a sample vector output routine. For a local environment c either replace the VECOUT call in the main line, or better c replace the contents of this subroutine between the *'s shown. c c There is often the requirement to distinguish each contour c line with a different colour or a different line style. This c can be done in many ways using the contour values z for a c particular line segment. c subroutine vecout(x1,y1,x2,y2,z) implicit none real*4 x1,y1,x2,y2,z c c***** Replace from here c c The following should be ignored since it is specific to c the version of FORTRAN running on the Macintosh microcomputer c on which this particular example was written. c INTEGER LINETO PARAMETER (LINETO=Z'89109000') INTEGER MOVETO PARAMETER (MOVETO=Z'89309000') call toolbx(MOVETO,nint(x1),nint(y1)) call toolbx(LINETO,nint(x2),nint(y2)) c c***** To here c return end ****** IBM BASIC CODE ******* 1000 REM ==================================================================== 1010 REM Documentation 1020 REM ============= 1030 REM 1040 REM Originally written in FORTRAN-77 1050 REM Translated to BASICA on the IBM-PC microcomputer 1060 REM By Paul D. Bourke, August 1987 1070 REM 1080 REM Input variables to CONREC 1090 REM d(0:iub,0:jub) matrix for the data surface height values 1100 REM iub, jub index bounds of the data array 1110 REM x(0:iub) data array for column coordinates 1120 REM y(0:jub) data array for row coordinates 1130 REM nc number of contour levels 1140 REM z(0:nc-1) contour levels in increasing order 1150 REM It is assumed that the logical variables TRUE and FALSE are 1160 REM either available or have been defined in the main line. 1170 REM 1180 REM NOTE: All arrays have a minimum index of 0 1190 REM 1200 REM Local declarations for CONREC 1210 REM ============================= 1220 DIM H(4) ' relative heights of the box above contour 1230 DIM ISH(4) ' sign of h() 1240 DIM XH(4) ' x coordinates of box 1250 DIM YH(4) ' y coordinates of box 1260 DIM IM(3) ' mapping from vertex numbers to x offsets 1270 IM(0)=0 : IM(1)=1 : IM(2)=1 : IM(3)=0 1280 DIM JM(3) ' mapping from vertex numbers to y offsets 1290 JM(0)=0 : JM(1)=0 : JM(2)=1 : JM(3)=1 1300 DIM CASTAB(2,2,2) 1310 DATA 0,0,8,0,2,5,7,6,9, 0,3,4,1,3,1,4,3,0, 9,6,7,5,2,0,8,0,0 1320 FOR K=0 TO 2 1330 FOR J=0 TO 2 1340 FOR I=0 TO 2 1350 READ CASTAB(K,J,I) 1360 NEXT I 1370 NEXT J 1380 NEXT K 1390 REM 1400 REM Check the input parameters for validity 1410 REM PRMERR should be tested by the calling routine and 1420 REM MSG$ printed if PRMERR is TRUE 1430 REM 1440 PRMERR = FALSE 1450 IF (IUB<=0 OR JUB<=0) THEN PRMERR = TRUE 1460 IF (NC<=0) THEN PRMERR = TRUE 1470 FOR K=1 TO NC-1 1480 IF (Z(K)<=Z(K-1)) THEN PRMERR = TRUE 1490 NEXT K 1500 IF (PRMERR) THEN MSG$ = "Error in input parameters" : RETURN 1510 REM 1520 REM Scan the array, top down, left to right 1530 REM ======================================= 1540 FOR J=JUB-1 TO 0 STEP -1 1550 FOR I=0 TO IUB-1 1560 REM Find the lowest vertex in the rectangle 1570 IF (D(I,J)D(I,J+1)) THEN DMAX = D(I,J) ELSE DMAX = D(I,J+1) 1620 IF (D(I+1,J)>DMAX) THEN DMAX = D(I+1,J) 1630 IF (D(I+1,J+1)>DMAX) THEN DMAX = D(I+1,J+1) 1640 IF (DMAXZ(NC-1)) THEN GOTO 2700 1650 REM 1660 REM Draw each contour within this box 1670 FOR K=0 TO NC-1 1680 IF ((Z(K)DMAX)) THEN GOTO 2690 1690 FOR M=4 TO 0 STEP -1 1700 IF M>0 THEN GOTO 1720 1710 IF M=0 THEN GOTO 1770 1720 REM m > 0 1730 H(M) = D(I+IM(M-1),J+JM(M-1))-Z(K) 1740 XH(M) = X(I+IM(M-1)) 1750 YH(M) = Y(J+JM(M-1)) 1760 GOTO 1820 1770 REM m = 0 1780 H(0) = (H(1)+H(2)+H(3)+H(4))/4 1790 XH(0) = (X(I)+X(I+1))/2 1800 YH(0) = (Y(J)+Y(J+1))/2 1810 GOTO 1820 1820 IF H(M) > 0 THEN ISH(M) = 2 1830 IF H(M) = 0 THEN ISH(M) = 1 1840 IF H(M) < 0 THEN ISH(M) = 0 1850 NEXT M 1860 REM 1870 REM Scan each triangle in the box 1880 FOR M=1 TO 4 1890 M1 = M 1900 M2 = 0 1910 M3 = M+1 1920 IF (M3 = 5) THEN M3 = 1 1930 REM 1940 REM Select the type of line/plane intersection 1950 CASE = CINT(CASTAB(ISH(M1),ISH(M2),ISH(M3))) 1960 IF CASE = 0 THEN GOTO 2680 1970 IF CASE = 1 THEN GOTO 2060 1980 IF CASE = 2 THEN GOTO 2120 1990 IF CASE = 3 THEN GOTO 2180 2000 IF CASE = 4 THEN GOTO 2240 2010 IF CASE = 5 THEN GOTO 2300 2020 IF CASE = 6 THEN GOTO 2360 2030 IF CASE = 7 THEN GOTO 2420 2040 IF CASE = 8 THEN GOTO 2480 2050 IF CASE = 9 THEN GOTO 2540 2060 REM Case 1: Line between vertices m1 and m2 2070 X1 = XH(M1) 2080 Y1 = YH(M1) 2090 X2 = XH(M2) 2100 Y2 = YH(M2) 2110 GOTO 2650 2120 REM Case 2: Line between vertices m2 and m3 2130 X1 = XH(M2) 2140 Y1 = YH(M2) 2150 X2 = XH(M3) 2160 Y2 = YH(M3) 2170 GOTO 2650 2180 REM Case 3: Line between vertices m3 and m1 2190 X1 = XH(M3) 2200 Y1 = YH(M3) 2210 X2 = XH(M1) 2220 Y2 = YH(M1) 2230 GOTO 2650 2240 REM Case 4: Line between vertex m1 and side m2-m3 2250 X1 = XH(M1) 2260 Y1 = YH(M1) 2270 X2 = (H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2)) 2280 Y2 = (H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2)) 2290 GOTO 2650 2300 REM Case 5: Line between vertex m2 and side m3-m1 2310 X1 = XH(M2) 2320 Y1 = YH(M2) 2330 X2 = (H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3)) 2340 Y2 = (H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3)) 2350 GOTO 2650 2360 REM Case 6: Line between vertex m3 and side m1-m2 2370 X1 = XH(M3) 2380 Y1 = YH(M3) 2390 X2 = (H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1)) 2400 Y2 = (H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1)) 2410 GOTO 2650 2420 REM CASE 7: Line between sides m1-m2 and m2-m3 2430 X1 = (H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1)) 2440 Y1 = (H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1)) 2450 X2 = (H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2)) 2460 Y2 = (H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2)) 2470 GOTO 2650 2480 REM Case 8: Line between sides m2-m3 and m3-m1 2490 X1 = (H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2)) 2500 Y1 = (H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2)) 2510 X2 = (H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3)) 2520 Y2 = (H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3)) 2530 GOTO 2650 2540 REM Case 9: Line between sides m3-m1 and m1-m2 2550 X1 = (H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3)) 2560 Y1 = (H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3)) 2570 X2 = (H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1)) 2580 Y2 = (H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1)) 2590 REM This is where a contour line segment of value z(i) 2600 REM is drawn between points (x1,y1) and (x2,y2) 2610 REM The exact details can vary depending on the particular 2620 REM graphics capabilities available. 2630 REM If colour is available then each contour level could 2640 REM be drawn in a new colour 2650 line (x1,y1)-(x2,y2) 2660 NEXT M 2670 NEXT K 2680 NEXT I 2690 NEXT J 2700 RETURN 2710 REM ===================================================================