home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Hack-Phreak Scene Programs
/
cleanhpvac.zip
/
cleanhpvac
/
FAQSYS18.ZIP
/
FAQS.DAT
/
CONREC.TXT
< prev
next >
Wrap
Text File
|
1991-05-27
|
24KB
|
634 lines
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 DMIN = D(I,J) ELSE DMIN = D(I,J+1)
1580 IF (D(I+1,J)<DMIN) THEN DMIN = D(I+1,J)
1590 IF (D(I+1,J+1)<DMIN) THEN DMIN = D(I+1,J+1)
1600 REM Find the highest vertex in the rectangle
1610 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 (DMAX<Z(0) OR DMIN>Z(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)<DMIN) OR (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 ===================================================================