Logo Search packages:      
Sourcecode: felt version File versions  Download package

vispoly.f

C
C     The following code was excerpted from: vispol.f
C
      SUBROUTINE VISPOL(XEYE,YEYE,NVRT,XC,YC,NVIS,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER NVIS,NVRT
      INTEGER IVIS(0:NVRT)
      DOUBLE PRECISION XEYE,YEYE
      DOUBLE PRECISION XC(0:NVRT),YC(0:NVRT)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: Compute the visibility polygon VP from an eyepoint in
C        the interior or blocked exterior of a simple polygon P or
C        on the boundary of a simply connected polygonal region P.
C      In the latter case, the interior angles at all vertices must
C        be strictly between 0 and 2*PI.
C
C     Input parameters:
C        XEYE,YEYE - coordinates of eyepoint; must be a simple vertex
C              if it lies on the boundary (i.e. occurs only once)
C      NVRT - upper subscript of XC, YC (approx number of vertices)
C      XC(0:NVRT),YC(0:NVRT) - If eyepoint is interior or blocked
C              exterior then arrays contain coordinates in CCW or CW
C              order, respectively, with (XC(0),YC(0)) = (XC(NVRT),
C              YC(NVRT)); (XC(0),YC(0)) is a vertex visible from
C              (XEYE,YEYE), e.g. as computed by routine ROTIPG.
C              If eyepoint is a vertex of P then arrays contain
C            coordinates in CCW order; (XC(0),YC(0)) is successor
C              vertex of (XEYE,YEYE); (XC(NVRT),YC(NVRT)) is
C              predecessor vertex of (XEYE,YEYE).
C
C     Updated parameters:
C        XC(0:NVIS),YC(0:NVIS) - vertices of VP in CCW order;
C              if eyepoint is interior or blocked exterior then
C              (XC(0),YC(0)) = (XC(NVIS),YC(NVIS)), else (XC(0),YC(0))
C              and (XC(NVIS),YC(NVIS)) are the successor and
C              predecessor vertices of (XEYE,YEYE) in VP
C
C     Output parameters:
C      NVIS - upper subscript of XC, YC on output (approx number
C              of vertices of VP); NVIS <= NVRT
C        IVIS(0:NVIS) - contains information about the vertices of VP
C              w.r.t. the vertices of P; IVIS(I) = K if (XC(I),YC(I))
C              is the vertex of index K in the input polygon; IVIS(I)
C              = -K if (XC(I),YC(I)) is on the interior of the edge
C              joining vertices of index K-1 and K in input polygon
C
C     Note about algorithm:
C        On input, XC and YC contain vertex coordinates of P. During
C        the algorithm, part of XC, YC is used as a stack, which, on
C        output, contains the vertex coordinates of VP. The stack
C        vertices overwrite the input vertices as the input vertices
C        are scanned. Elements of IVIS are set when vertices are added
C        to the stack; these values may have +NV or -NV added to them
C        to indicate that stack point has same angle as previous one.
C
C     Reference: 
C        B. Joe and R. B. Simpson, BIT 27 (1987), pp. 458-473.
C
C     Abnormal return:
C        IERR is set to 206, 207, 208, 209, or 210
C
C     Routines called:
C        LRLINE, VPLEFT, VPRGHT, VPSCNA, VPSCNB, VPSCNC, VPSCND
C
      INTEGER CUR,IERR,NV,OPER,TOP
      DOUBLE PRECISION XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GERROR/ IERR
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GERROR/,/GVPVAR/
C
C     Variables in common block GVPVAR:
C        NV - NVRT
C        OPER - operation code 1 to 7 for LEFT, RIGHT, SCANA, SCANB,
C              SCANC, SCAND, FINISH
C        CUR - index of current vertex of P in XC, YC arrays
C        TOP - index of top vertex of stack in XC, YC arrays
C              (TOP <= CUR is always satisfied)
C        XE,YE - XEYE,YEYE
C        XW,YW - coordinates of point on last or second-last edge
C              processed (needed for routines VPSCNB, VPSCNC, VPSCND)
C        BEYE - .TRUE. iff eyepoint is on boundary
C
      INTEGER I,LR,LRLINE
C
      BEYE = XC(0) .NE. XC(NVRT) .OR. YC(0) .NE. YC(NVRT)
      NV = NVRT
      XE = XEYE
      YE = YEYE
      IVIS(0) = 0
      CUR = 1
      IF (.NOT. BEYE) GO TO 20
   10 CONTINUE
       LR = LRLINE(XC(NV-1),YC(NV-1),XE,YE,XC(NV),YC(NV),0.0D0)
       IF (LR .EQ. 0) THEN
          NV = NV - 1
          GO TO 10
       ENDIF
C
   20 CONTINUE
         LR = LRLINE(XC(CUR),YC(CUR),XE,YE,XC(0),YC(0),0.0D0)
       IF (LR .EQ. 0) THEN
          CUR = CUR + 1
          GO TO 20
       ENDIF
      IF (LR .EQ. -1) THEN
       OPER = 1
       IF (CUR .EQ. 1) THEN
            TOP = 1
            IVIS(1) = CUR
       ELSE IF (BEYE) THEN
          TOP = 1
          XC(0) = XC(CUR-1)
          YC(0) = YC(CUR-1)
          IVIS(0) = CUR - 1
          XC(1) = XC(CUR)
          YC(1) = YC(CUR)
          IVIS(1) = CUR
       ELSE
          TOP = 2
          XC(1) = XC(CUR-1)
          YC(1) = YC(CUR-1)
          IVIS(1) = CUR - 1 + NV
          XC(2) = XC(CUR)
          YC(2) = YC(CUR)
          IVIS(2) = CUR
       ENDIF
      ELSE
       OPER = 3
       TOP = 0
       IF (BEYE .AND. CUR. GT. 1) THEN
          XC(0) = XC(CUR-1)
          YC(0) = YC(CUR-1)
          IVIS(0) = CUR - 1
       ENDIF
      ENDIF
C
C     Angular displacement of stack points are in nondecreasing order,
C     with at most two consecutive points having the same displacement.
C
   30 CONTINUE
       IF (OPER .EQ. 1) THEN
          CALL VPLEFT(XC,YC,IVIS)
       ELSE IF (OPER .EQ. 2) THEN
          CALL VPRGHT(XC,YC,IVIS)
       ELSE IF (OPER .EQ. 3) THEN
          CALL VPSCNA(XC,YC,IVIS)
       ELSE IF (OPER .EQ. 4) THEN
          CALL VPSCNB(XC,YC,IVIS)
       ELSE IF (OPER .EQ. 5) THEN
          CALL VPSCNC(XC,YC,IVIS)
       ELSE
          CALL VPSCND(XC,YC,IVIS)
       ENDIF
       IF (IERR .NE. 0) THEN
          NVIS = TOP
          RETURN
       ENDIF
      IF (OPER .LE. 6) GO TO 30
C
C     Add or subtract NV from those IVIS values which are used to
C     indicate that stack point has same angle as previous one.
C
      DO 40 I = 1,TOP
       IF (IVIS(I) .GT. NV) THEN
          IVIS(I) = IVIS(I) - NV
       ELSE IF (IVIS(I) .LT. -NV) THEN
          IVIS(I) = IVIS(I) + NV
       ENDIF
   40 CONTINUE
      NVIS = TOP
      END
C
C     The following code was excerpted from: visvrt.f
C
      SUBROUTINE VISVRT(ANGSPC,XEYE,YEYE,NVIS,XC,YC,IVIS,MAXN,NVSVRT,
     $   THETA)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXN,NVIS,NVSVRT
      INTEGER IVIS(0:MAXN)
      DOUBLE PRECISION ANGSPC,XEYE,YEYE
      DOUBLE PRECISION THETA(0:MAXN),XC(0:MAXN),YC(0:MAXN)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: Determine a list of visible vertices, ordered by
C        increasing "polar angle", on the boundary of the visibilty
C        polygon from boundary eyepoint (XEYE,YEYE). This list
C        includes the vertices of visibility polygon such that a
C        line segment from (XEYE,YEYE) to vertex lies in interior
C        of polygon, as well as extra points on edges which subtend
C        an angle >= 2*ANGSPC at (XEYE,YEYE). These extra points are
C        at an equal angular spacing of >= ANGSPC and < 2*ANGSPC. The
C        successor and predecessor of eyepoint are included in list.
C
C     Input parameters:
C        ANGSPC - angle spacing parameter in radians which controls
C              how many extra points become visible vertices
C        XEYE,YEYE - coordinates of boundary eyepoint
C        NVIS - (number of vertices of visibility polygon) - 2
C        XC(0:NVIS),YC(0:NVIS) - the coordinates of the vertices of
C              visibility polygon in CCW order; (XC(0),YC(0)) and
C              (XC(NVIS),YC(NVIS)) are the successor and predecessor
C              vertices of eyepoint in visibility polygon; at most 2
C              consecutive vertices have same polar angle wrt eyepoint
C        IVIS(0:NVIS) - contains information about the vertices of
C              XC, YC arrays with respect to the original polygon from
C              which visibility polygon is computed; if IVIS(I) >= 0
C              then (XC(I),YC(I)) has index I in original polygon;
C              if IVIS(I) < 0 then (XC(I),YC(I)) is on the edge
C              ending at vertex of index -IVIS(I) in original polygon;
C              indexing starts at 0 from successor of eyepoint
C        MAXN - upper bound on NVSVRT; should be at least
C              NVIS + INT(PHI/ANGSPC) where PHI is the interior
C              angle at (XEYE,YEYE)
C
C     Updated parameters:
C        XC(0:NVSVRT),YC(0:NVSVRT) - coordinates of visible vertices
C              which overwrite the input coordinates
C        IVIS(0:NVSVRT) - contains information about the output
C              vertices of XC, YC arrays as described above for input
C
C     Output parameters:
C        NVSVRT - (number of visible vertices) - 1
C        THETA(0:NVSVRT) - polar angles of visible vertices wrt (XEYE,
C              YEYE) at origin and (XC(0),YC(0)) on positive x-axis
C
C     Routines called:
C        ANGLE, LRLINE
C
      DOUBLE PRECISION PI,TOL
      COMMON /GCONST/ PI,TOL
      SAVE /GCONST/
C
      INTEGER CUR,I,IND,K,LR,LRLINE,N,TOP
      DOUBLE PRECISION ALPHA,ANG,ANG1,ANG2,ANGDIF,ANGLE,ANGSP2
      DOUBLE PRECISION COSANG,DIFF,DX,DY,NUMER,R,SINANG
C
C     Shift input vertices right, and possibly remove first and last
C     vertices due to collinearity with eyepoint.
C
      ANGSP2 = 2.0D0*ANGSPC
      CUR = MAXN + 1
      N = MAXN
      DO 10 I = NVIS,0,-1
       CUR = CUR - 1
       XC(CUR) = XC(I)
       YC(CUR) = YC(I)
       IVIS(CUR) = IVIS(I)
   10 CONTINUE
      LR = LRLINE(XC(CUR+1),YC(CUR+1),XEYE,YEYE,XC(CUR),YC(CUR),0.0D0)
      IF (LR .GE. 0) THEN
       CUR = CUR + 1
       XC(0) = XC(CUR)
       YC(0) = YC(CUR)
       IVIS(0) = IVIS(CUR)
      ENDIF
      LR = LRLINE(XC(N-1),YC(N-1),XEYE,YEYE,XC(N),YC(N),0.0D0)
      IF (LR .LE. 0) N = N - 1
      ALPHA = ATAN2(YC(0)-YEYE,XC(0)-XEYE)
      ANG2 = 0.0D0
      THETA(0) = 0.0D0
      TOP = 0
      CUR = CUR + 1
C
C     Process edge from vertices of indices CUR-1, CUR.
C
   20 CONTINUE
       ANG1 = ANG2
       ANG2 = ANGLE(XC(CUR),YC(CUR),XEYE,YEYE,XC(0),YC(0))
         ANGDIF = ANG2 - ANG1
         IF (ANGDIF .LE. TOL) THEN
          DIFF = ((XC(CUR) - XEYE)**2 + (YC(CUR) - YEYE)**2) -
     $         ((XC(CUR-1) - XEYE)**2 + (YC(CUR-1) - YEYE)**2)
          IF (DIFF .LT. 0.0D0) THEN
             XC(TOP) = XC(CUR)
             YC(TOP) = YC(CUR)
             IVIS(TOP) = IVIS(CUR)
             THETA(TOP) = ANG2
          ENDIF
         ELSE
          IF (ANGDIF .GE. ANGSP2) THEN
             K = INT(ANGDIF/ANGSPC)
             IND = -ABS(IVIS(CUR))
             ANGDIF = ANGDIF/DBLE(K)
             DX = XC(CUR) - XC(CUR-1)
             DY = YC(CUR) - YC(CUR-1)
             NUMER = (XC(CUR) - XEYE)*DY - (YC(CUR) - YEYE)*DX
             DO 30 I = 1,K-1
                TOP = TOP + 1
              THETA(TOP) = ANG1 + DBLE(I)*ANGDIF
              ANG = THETA(TOP) + ALPHA
                COSANG = COS(ANG)
                SINANG = SIN(ANG)
                R = NUMER/(DY*COSANG - DX*SINANG)
                XC(TOP) = R*COSANG + XEYE
                YC(TOP) = R*SINANG + YEYE
                IVIS(TOP) = IND
   30          CONTINUE
          ENDIF
          TOP = TOP + 1
          XC(TOP) = XC(CUR)
          YC(TOP) = YC(CUR)
          IVIS(TOP) = IVIS(CUR)
          THETA(TOP) = ANG2
         ENDIF
         CUR = CUR + 1
      IF (CUR .LE. N) GO TO 20
      NVSVRT = TOP
      END
C
C     The following code was excerpted from: vornbr.f
C
      SUBROUTINE VORNBR(XEYE,YEYE,NVRT,XC,YC,NVOR,IVOR,XVOR,YVOR)
      IMPLICIT LOGICAL (A-Z)
      INTEGER NVOR,NVRT
      INTEGER IVOR(0:NVRT)
      DOUBLE PRECISION XEYE,YEYE
      DOUBLE PRECISION XC(0:NVRT),XVOR(0:NVRT),YC(0:NVRT),YVOR(0:NVRT)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: Determine the Voronoi neighbours of (XEYE,YEYE) from a
C        list of vertices which are in increasing "polar angle" order.
C        The Voronoi neighbours are a sublist of this list. The
C        Voronoi polygon is restricted to the sector formed from the
C        the edges joining (XEYE,YEYE) to the first and last vertices
C        of this list. Each Voronoi neighbour corresponds to an edge
C        of the Voronoi polygon.
C
C     Input parameters:
C        XEYE,YEYE - coordinates of eyepoint
C        NVRT - (number of vertices in list) - 1
C        XC(0:NVRT),YC(0:NVRT) - vertex coordinates from which
C              Voronoi neighbours are determined; (XC(0),YC(0)),...,
C              (XC(NVRT),YC(NVRT)) are in increasing angular
C              displacement order w.r.t. (XEYE,YEYE)
C
C     Output parameters:
C        NVOR - (number of Voronoi neighbours) - 1 [<= NVRT]
C        IVOR(0:NVOR) - indices of Voronoi neighbours in XC, YC
C               arrays; 0 <= IVOR(0) < ... < IVOR(NVOR) <= NVRT
C
C     Working parameters:
C        XVOR(0:NVRT),YVOR(0:NVRT) - arrays for storing the vertex
C              coordinates of the Voronoi polygon
C
C     Abnormal return:
C        IERR is set to 212
C
C     Routines called:
C        LRLINE
C
      INTEGER IERR
      DOUBLE PRECISION PI,TOL
      COMMON /GERROR/ IERR
      COMMON /GCONST/ PI,TOL
      SAVE /GERROR/,/GCONST/
C
      INTEGER IM,K,LR,LRLINE,M
      DOUBLE PRECISION A11,A12,A21,A22,B1,B2,DET,TOLABS,XI,YI
C
      K = 1
      M = 0
      IVOR(0) = 0
      XVOR(0) = (XEYE + XC(0))*0.5D0
      YVOR(0) = (YEYE + YC(0))*0.5D0
C
C     Beginning of main loop
C
   10 CONTINUE
      IF (K .GT. NVRT) GO TO 20
C
C        Determine the intersection of the perpendicular bisectors
C        of edges from (XEYE,YEYE) to (XC(K),YC(K)) and from
C        (XEYE,YEYE) to (XC(IM),YC(IM)).
C
       IM = IVOR(M)
       A11 = XC(K) - XEYE
       A12 = YC(K) - YEYE
       A21 = XC(IM) - XEYE
       A22 = YC(IM) - YEYE
       TOLABS = TOL*MAX(ABS(A11),ABS(A12),ABS(A21),ABS(A22))
       DET = A11*A22 - A21*A12
       IF (ABS(DET) .LE. TOLABS) THEN
          IERR = 212
          RETURN
       ENDIF
       B1 = (A11**2 + A12**2)*0.5D0
       B2 = (A21**2 + A22**2)*0.5D0
       XI = (B1*A22 - B2*A12)/DET
       YI = (B2*A11 - B1*A21)/DET
C
C        Determine whether (XVOR(M+1),YVOR(M+1)) is to the left of or
C        on the directed line from (XEYE,YEYE) to (XVOR(M),YVOR(M)).
C
       XVOR(M+1) = XI + XEYE
       YVOR(M+1) = YI + YEYE
       LR = LRLINE(XVOR(M+1),YVOR(M+1),XEYE,YEYE,XVOR(M),YVOR(M),
     1      0.0D0)
       IF (LR .LE. 0) THEN
          M = M + 1
          IVOR(M) = K
          K = K + 1
       ELSE IF (M .GT. 0) THEN
          M = M - 1
       ELSE
C
C           Determine the intersection of edge from (XEYE,YEYE) to
C           (XC(0),YC(0)) and the perpendicular bisector of the edge
C           from (XEYE,YEYE) to (XC(K),YC(K)).
C
          A11 = XC(K) - XEYE
          A12 = YC(K) - YEYE
          A21 = YC(0) - YEYE
          A22 = XEYE - XC(0)
          TOLABS = TOL*MAX(ABS(A11),ABS(A12),ABS(A21),ABS(A22))
          DET = A11*A22 - A21*A12
          IF (ABS(DET) .LE. TOLABS) THEN
             IERR = 212
             RETURN
          ENDIF
          B1 = (A11**2 + A12**2)*0.5D0
          B2 = 0.0D0
          XI = (B1*A22 - B2*A12)/DET
          YI = (B2*A11 - B1*A21)/DET
          XVOR(M) = XI + XEYE
          YVOR(M) = YI + YEYE
          IVOR(M) = K
          K = K + 1
       ENDIF
      GO TO 10
C
C     The following short loop determines which vertices at the end
C     of list are not Voronoi neighbours.
C
   20 CONTINUE
       LR = LRLINE(XVOR(M),YVOR(M),XEYE,YEYE,XC(NVRT),YC(NVRT),0.0D0)
       IF (LR .GE. 0) GO TO 30
       M = M - 1
      IF (M .GE. 0) GO TO 20
   30 CONTINUE
      NVOR = M
      END
C
C     The following code was excerpted from: vpleft.f
C
      SUBROUTINE VPLEFT(XC,YC,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IVIS(0:*)
      DOUBLE PRECISION XC(0:*),YC(0:*)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: This routine is called by routine VISPOL for the LEFT
C        operation (OPER = 1).
C
C     Input and updated parameters:
C        XC,YC,IVIS - see comments in routine VISPOL
C
C     Routines called:
C        LRLINE, XEDGE
C
      INTEGER CUR,NV,OPER,TOP
      DOUBLE PRECISION XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GVPVAR/
C
      INTEGER J,LR,LR1,LR2,LRLINE
      DOUBLE PRECISION XU,YU
      LOGICAL INTSCT
C
C     EYE-V(CUR-1)-V(CUR) is a left turn, S(TOP) = V(CUR), TOP <= CUR,
C     S(TOP-1) = V(CUR-1) or on interior of edge V(CUR-1)-V(CUR).
C
   10 CONTINUE
      IF (CUR .EQ. NV) THEN
       OPER = 7
       RETURN
      ENDIF
      IF (.NOT. BEYE .AND. TOP .LE. 2) GO TO 20
C
C     Check if angular displacement of stack chain >= 2*PI or
C     interior angle at boundary viewpoint.
C
      CALL XEDGE(1,XE,YE,XC(NV),YC(NV),XC(TOP-1),YC(TOP-1),XC(TOP),
     $   YC(TOP),XU,YU,INTSCT)
      IF (INTSCT) THEN
       OPER = 4
       XW = XC(CUR)
       YW = YC(CUR)
         LR = LRLINE(XC(TOP),YC(TOP),XE,YE,XC(NV),YC(NV),0.0D0)
       IF (LR .EQ. -1) THEN
          XC(TOP) = XU
          YC(TOP) = YU
          IVIS(TOP) = -CUR
       ENDIF
       RETURN
      ENDIF
C
C     Process next edge.
C
   20 CONTINUE
      LR = LRLINE(XC(CUR+1),YC(CUR+1),XE,YE,XC(CUR),YC(CUR),0.0D0)
      IF (LR .EQ. -1) THEN
       CUR = CUR + 1
       TOP = TOP + 1
       XC(TOP) = XC(CUR)
       YC(TOP) = YC(CUR)
       IVIS(TOP) = CUR
      ELSE
       J = CUR + 1
         LR1 = LRLINE(XC(J),YC(J),XC(TOP-1),YC(TOP-1),XC(CUR),YC(CUR),
     $      0.0D0)
       IF (LR1 .EQ. 1) THEN
          OPER = 3
          CUR = J
       ELSE
          IF (LR .EQ. 1) THEN
             LR2 = 1
             GO TO 40
          ENDIF
   30       CONTINUE
             J = J + 1
             LR2 = LRLINE(XC(J),YC(J),XE,YE,XC(CUR),YC(CUR),0.0D0)
          IF (LR2 .EQ. 0) GO TO 30
   40       CONTINUE
          IF (LR2 .EQ. -1) THEN
             TOP = TOP + 1
             XC(TOP) = XC(J-1)
             YC(TOP) = YC(J-1)
             IVIS(TOP) = J - 1 + NV
             TOP = TOP + 1
             XC(TOP) = XC(J)
             YC(TOP) = YC(J)
             IVIS(TOP) = J
          ELSE
             OPER = 2
          ENDIF
          CUR = J
       ENDIF
      ENDIF
C
C     This test avoids extra subroutine calls.
C
      IF (OPER .EQ. 1) GO TO 10
      END
C
C     The following code was excerpted from: vprght.f
C
      SUBROUTINE VPRGHT(XC,YC,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IVIS(0:*)
      DOUBLE PRECISION XC(0:*),YC(0:*)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: This routine is called by routine VISPOL for the RIGHT
C        operation (OPER = 2).
C
C     Input and updated parameters:
C        XC,YC,IVIS - see comments in routine VISPOL
C
C     Abnormal return:
C        IERR is set to 206
C
C     Routines called:
C        LRLINE, XEDGE
C
      INTEGER CUR,IERR,NV,OPER,TOP
      DOUBLE PRECISION XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GERROR/ IERR
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GERROR/,/GVPVAR/
C
      INTEGER CASE,J,LR,LR1,LR2,LRLINE
      DOUBLE PRECISION XU,YU
      LOGICAL INTSCT
C
C     EYE-V(CUR-1)-V(CUR) is a right turn, EYE-S(TOP)-V(CUR) is a right
C     turn, EYE-S(TOP-1)-S(TOP) is a left turn, TOP < CUR, S(TOP) =
C     V(CUR-1) and S(TOP-1)-S(TOP)-V(CUR) is a left turn or S(TOP) is
C     not on edge V(CUR-1)-V(CUR) and V(CUR-1)-V(CUR) intersects
C     EYE-S(TOP).
C     Pop points from stack. If BEYE, it is not possible for
C     (XC(CUR),YC(CUR)) to be identical to any stack points.
C
   10 CONTINUE
      CASE = 0
      J = TOP
   20 CONTINUE
       IF (ABS(IVIS(J)) .LE. NV) THEN
          LR = LRLINE(XC(CUR),YC(CUR),XE,YE,XC(J-1),YC(J-1),0.0D0)
          IF (LR .EQ. -1) THEN
             CASE = 1
          ELSE IF (LR .EQ. 0) THEN
             IF (ABS(IVIS(J-1)) .LE. NV) THEN
              J = J - 1
              CASE = 2
             ELSE IF ((XC(J-2) - XE)**2 + (YC(J-2) - YE)**2 .GE.
     $            (XC(J-1) - XE)**2 + (YC(J-1) - YE)**2) THEN
              J = J - 2
              CASE = 2
             ELSE
              CASE = -1
             ENDIF
          ENDIF
       ELSE IF (CASE .EQ. -1) THEN
          IF ((XC(J-1) - XE)**2 + (YC(J-1) - YE)**2 .GE.
     $         (XC(CUR) - XE)**2 + (YC(CUR) - YE)**2) THEN
             J = J - 1
             CASE = 2
          ELSE
             XW = XC(CUR)
             YW = YC(CUR)
             CASE = 3
          ENDIF
       ELSE
          CALL XEDGE(0,XC(CUR-1),YC(CUR-1),XC(CUR),YC(CUR),
     $         XC(J-1),YC(J-1),XC(J),YC(J),XW,YW,INTSCT)
          IF (INTSCT) CASE = 3
       ENDIF
       IF (CASE .GT. 0) GO TO 30
       J = J - 1
      IF (J .GE. 1) GO TO 20
C
C     Error from no more edges in stack.
C
      IERR = 206
      RETURN
C
C     Process next edge.
C
   30 CONTINUE
      IF (CASE .EQ. 3) THEN
       OPER = 6
       TOP = J - 1
      ELSE
       TOP = J
       XW = XC(CUR-1)
       YW = YC(CUR-1)
       IF (CASE .EQ. 1) THEN
          CALL XEDGE(1,XE,YE,XC(CUR),YC(CUR),XC(TOP-1),YC(TOP-1),
     $         XC(TOP),YC(TOP),XU,YU,INTSCT)
          XC(TOP) = XU
          YC(TOP) = YU
          IVIS(TOP) = -ABS(IVIS(TOP))
       ENDIF
       LR = LRLINE(XC(CUR+1),YC(CUR+1),XE,YE,XC(CUR),YC(CUR),0.0D0)
       IF (LR .EQ. 1) THEN
          CUR = CUR + 1
       ELSE
          J = CUR + 1
          LR1 = LRLINE(XC(J),YC(J),XW,YW,XC(CUR),YC(CUR),0.0D0)
          IF (LR1 .EQ. -1) THEN
             OPER = 5
             CUR = J
          ELSE
             IF (LR .EQ. -1) THEN
              LR2 = -1
              GO TO 50
             ENDIF
   40          CONTINUE
              J = J + 1
              LR2 = LRLINE(XC(J),YC(J),XE,YE,XC(CUR),YC(CUR),0.0D0)
             IF (LR2 .EQ. 0) GO TO 40
   50          CONTINUE
             IF (LR2 .EQ. -1) THEN
                OPER = 1
                TOP = TOP + 1
                XC(TOP) = XC(J-1)
                YC(TOP) = YC(J-1)
                IVIS(TOP) = J - 1 + NV
                TOP = TOP + 1
                XC(TOP) = XC(J)
                YC(TOP) = YC(J)
                IVIS(TOP) = J
             ENDIF
             CUR = J
          ENDIF
       ENDIF
      ENDIF
C
C     This test avoids extra subroutine calls.
C
      IF (OPER .EQ. 2) GO TO 10
      END
C
C     The following code was excerpted from: vpscna.f
C
      SUBROUTINE VPSCNA(XC,YC,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IVIS(0:*)
      DOUBLE PRECISION XC(0:*),YC(0:*)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: This routine is called by routine VISPOL for the SCANA
C        operation (OPER = 3).
C
C     Input and updated parameters:
C        XC,YC,IVIS - see comments in routine VISPOL
C
C     Abnormal return:
C        IERR is set to 207
C
C     Routines called:
C        LRLINE, XEDGE
C
      INTEGER CUR,IERR,NV,OPER,TOP
      DOUBLE PRECISION XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GERROR/ IERR
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GERROR/,/GVPVAR/
C
      INTEGER CASE,J,K,LR,LR1,LR2,LR3,LRLINE
      LOGICAL INTSCT
C
C     EYE-V(CUR-1)-V(CUR) is a right turn or forward move, S(TOP) =
C     V(CUR-1) or EYE-S(TOP)-V(CUR-1) is a forward move and TOP = 0,
C     TOP < CUR; S(TOP-1)-S(TOP)-V(CUR) is a right turn if TOP >= 1
C     or EYE-S(TOP)-V(CUR) is a right turn if TOP = 0.
C     If BEYE, it is possible that (XC(TOP),YC(TOP)) is a non-simple
C     vertex but any edge incident on this vertex encountered during
C     scan must be invisible from (XE,YE).
C
      K = CUR
   10 CONTINUE
       IF (XC(K+1) .EQ. XC(TOP) .AND. YC(K+1) .EQ. YC(TOP)) THEN
          K = K + 2
       ELSE
          CALL XEDGE(1,XE,YE,XC(TOP),YC(TOP),XC(K),YC(K),XC(K+1),
     $         YC(K+1),XW,YW,INTSCT)
          IF (INTSCT) THEN
             LR = LRLINE(XC(K+1),YC(K+1),XE,YE,XC(K),YC(K),0.0D0)
             IF (LR .EQ. 1) THEN
              IF ((XC(TOP) - XE)**2 + (YC(TOP) - YE)**2 .GE.
     $               (XW - XE)**2 + (YW - YE)**2) THEN
                 IF (TOP .GT. 0) THEN
                    CASE = 1
                  GO TO 20
                 ENDIF
              ELSE
                   LR1 = LRLINE(XC(K),YC(K),XE,YE,XC(TOP),YC(TOP),
     $                  0.0D0)
                 IF (LR1 .EQ. -1) THEN
                    CASE = 2
                    GO TO 20
                 ENDIF
              ENDIF
             ELSE
                LR1 = LRLINE(XC(K+1),YC(K+1),XE,YE,XC(TOP),YC(TOP),
     $               0.0D0)
              IF (LR1 .EQ. -1) THEN
                 CASE = 3
                 GO TO 20
              ENDIF
             ENDIF
          ENDIF
          K = K + 1
       ENDIF
      IF (K .LT. NV) GO TO 10
C
C     Error from unsuccessful scan.
C
      IERR = 207
      RETURN
C
C     Process current edge.
C
   20 CONTINUE
      IF (CASE .EQ. 3) THEN
       OPER = 1
       CUR = K + 1
       LR = LRLINE(XC(K),YC(K),XE,YE,XC(TOP),YC(TOP),0.0D0)
       TOP = TOP + 1
       IF (LR .EQ. 0) THEN
          XC(TOP) = XC(K)
          YC(TOP) = YC(K)
          IVIS(TOP) = K + NV
       ELSE
          XC(TOP) = XW
          YC(TOP) = YW
          IVIS(TOP) = -(K + 1 + NV)
       ENDIF
       TOP = TOP + 1
       XC(TOP) = XC(CUR)
       YC(TOP) = YC(CUR)
       IVIS(TOP) = CUR
      ELSE IF (CASE .EQ. 1) THEN
       CUR = K + 1
       LR = LRLINE(XC(CUR),YC(CUR),XE,YE,XC(TOP),YC(TOP),0.0D0)
       IF (LR .EQ. 1) THEN
          OPER = 2
       ELSE
          J = CUR + 1
          LR1 = LRLINE(XC(J),YC(J),XE,YE,XC(CUR),YC(CUR),0.0D0)
          LR2 = LRLINE(XC(J),YC(J),XC(K),YC(K),XC(CUR),YC(CUR),0.0D0)
          IF (LR1 .LE. 0 .AND. LR2 .EQ. -1) THEN
             OPER = 5
             XW = XC(K)
             YW = YC(K)
             CUR = J
          ELSE
             IF (LR1 .NE. 0) THEN
              LR3 = LR1
              GO TO 40
             ENDIF
   30          CONTINUE
              J = J + 1
              LR3 = LRLINE(XC(J),YC(J),XE,YE,XC(CUR),YC(CUR),0.0D0)
             IF (LR3 .EQ. 0) GO TO 30
   40          CONTINUE
             IF (LR3 .EQ. 1) THEN
              OPER = 2
             ELSE
              OPER = 1
              TOP = TOP + 1
              XC(TOP) = XC(J-1)
              YC(TOP) = YC(J-1)
              IVIS(TOP) = J - 1 + NV
              TOP = TOP + 1
              XC(TOP) = XC(J)
              YC(TOP) = YC(J)
              IVIS(TOP) = J
             ENDIF
             CUR = J
          ENDIF
       ENDIF
      ELSE
       OPER = 6
       CUR = K + 1
       LR = LRLINE(XC(CUR),YC(CUR),XE,YE,XC(TOP),YC(TOP),0.0D0)
       IF (LR .EQ. 0) THEN
          XW = XC(CUR)
          YW = YC(CUR)
       ENDIF
      ENDIF
      END
C
C     The following code was excerpted from: vpscnb.f
C
      SUBROUTINE VPSCNB(XC,YC,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IVIS(0:*)
      DOUBLE PRECISION XC(0:*),YC(0:*)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: This routine is called by routine VISPOL for the SCANB
C        operation (OPER = 4).
C
C     Input and updated parameters:
C        XC,YC,IVIS - see comments in routine VISPOL
C
C     Abnormal return:
C        IERR is set to 208
C
C     Routines called:
C        LRLINE, XEDGE
C
      INTEGER CUR,IERR,NV,OPER,TOP
      DOUBLE PRECISION PI,TOL,XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GERROR/ IERR
      COMMON /GCONST/ PI,TOL
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GERROR/,/GCONST/,/GVPVAR/
C
      INTEGER K,LR,LR1,LRLINE
      DOUBLE PRECISION TOLABS,XU,YU
      LOGICAL INTSCT
C
C     EYE-V(CUR-1)-V(CUR) is a left turn, S(TOP) = V(CUR) or S(TOP) is
C     on interior of edge V(CUR-1)-V(CUR), TOP <= CUR, S(TOP) has
C     angular displacement of 2*PI or interior angle at boundary eye.
C     (XW,YW) is the input version of (XC(CUR),YC(CUR)).
C     If BEYE, it is possible that (XC(TOP),YC(TOP)) is a non-simple
C     point but any edge containing this point encountered during scan
C     must be invisible from (XE,YE), except for 1 case where K = CUR.
C
      TOLABS = TOL*((XC(NV) - XC(TOP))**2 + (YC(NV) - YC(TOP))**2)
      K = CUR
      IF (IVIS(TOP) .LT. 0 .OR. K + 1 .EQ. NV) GO TO 10
      LR = LRLINE(XC(K+1),YC(K+1),XE,YE,XC(TOP),YC(TOP),0.0D0)
      LR1 = LRLINE(XC(K+1),YC(K+1),XC(TOP-1),YC(TOP-1),XC(TOP),YC(TOP),
     $   0.0D0)
      IF (LR .EQ. 1 .AND. LR1 .EQ. -1) THEN
       OPER = 2
       CUR = K + 1
       RETURN
      ELSE
       K = K + 1
      ENDIF
C
   10 CONTINUE
       IF (K + 1 .EQ. NV) THEN
          OPER = 7
          CUR = NV
          TOP = TOP + 1
          XC(TOP) = XC(NV)
          YC(TOP) = YC(NV)
          IVIS(TOP) = NV
          RETURN
       ELSE
          IF (K .EQ. CUR) THEN
             CALL XEDGE(0,XC(NV),YC(NV),XC(TOP),YC(TOP),XW,YW,
     $            XC(K+1),YC(K+1),XU,YU,INTSCT)
          ELSE
             CALL XEDGE(0,XC(NV),YC(NV),XC(TOP),YC(TOP),XC(K),YC(K),
     $            XC(K+1),YC(K+1),XU,YU,INTSCT)
          ENDIF
          IF (INTSCT) THEN
             IF ((XC(TOP) - XU)**2 + (YC(TOP) - YU)**2 .LE. TOLABS)
     $            GO TO 20
             LR = LRLINE(XC(K+1),YC(K+1),XE,YE,XC(NV),YC(NV),0.0D0)
             IF (LR .EQ. 1) THEN
              OPER = 2
              CUR = K + 1
              RETURN
             ENDIF
          ENDIF
   20       CONTINUE
          K = K + 1
       ENDIF
      IF (K .LT. NV) GO TO 10
C
C     Error from unsuccessful scan.
C
      IERR = 208
      END
C
C     The following code was excerpted from: vpscnc.f
C
      SUBROUTINE VPSCNC(XC,YC,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IVIS(0:*)
      DOUBLE PRECISION XC(0:*),YC(0:*)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: This routine is called by routine VISPOL for the SCANC
C        operation (OPER = 5).
C
C     Input and updated parameters:
C        XC,YC,IVIS - see comments in routine VISPOL
C
C     Abnormal return:
C        IERR is set to 209
C
C     Routines called:
C        LRLINE, XEDGE
C
      INTEGER CUR,IERR,NV,OPER,TOP
      DOUBLE PRECISION XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GERROR/ IERR
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GERROR/,/GVPVAR/
C
      INTEGER J,K,LR,LR1,LR2,LRLINE
      DOUBLE PRECISION XP,XU,YP,YU
      LOGICAL INTSCT
C
C     EYE-V(CUR-1)-V(CUR) is a left turn or forward move, EYE-V(CUR-2)-
C     V(CUR-1) is a right turn, V(CUR-2)-V(CUR-1)-V(CUR) is a left turn,
C     TOP < CUR-1, W = V(CUR-2), S(TOP) is not on V(CUR-1)-V(CUR), EYE-
C     S(TOP)-V(CUR-1) is a backward move, EYE-S(TOP-1)-S(TOP) is a left
C     turn. If BEYE, it is possible that V(CUR-1) is a non-simple point,
C     but intersection at (XC(TOP),YC(TOP)) cannot occur.
C
      XP = XC(CUR-1)
      YP = YC(CUR-1)
      K = CUR
   10 CONTINUE
       IF (XC(K+1) .EQ. XP .AND. YC(K+1) .EQ. YP) THEN
          GO TO 40
       ELSE IF (XC(K) .EQ. XP .AND. YC(K) .EQ. YP) THEN
          J = K + 1
          LR = LRLINE(XC(J),YC(J),XE,YE,XP,YP,0.0D0)
          LR1 = LRLINE(XC(J),YC(J),XW,YW,XP,YP,0.0D0)
          IF (LR .LE. 0 .AND. LR1 .EQ. -1) GO TO 40
          IF (LR .NE. 0) THEN
             LR2 = LR
             GO TO 30
          ENDIF
   20       CONTINUE
             J = J + 1
             LR2 = LRLINE(XC(J),YC(J),XE,YE,XP,YP,0.0D0)
          IF (LR2 .EQ. 0) GO TO 20
   30       CONTINUE
          IF (LR2 .EQ. 1) THEN
             OPER = 2
          ELSE
             OPER = 1
             TOP = TOP + 1
             XC(TOP) = XC(J-1)
             YC(TOP) = YC(J-1)
             IVIS(TOP) = J - 1 + NV
             TOP = TOP + 1
             XC(TOP) = XC(J)
             YC(TOP) = YC(J)
             IVIS(TOP) = J
          ENDIF
          CUR = J
          RETURN
       ELSE
          CALL XEDGE(0,XP,YP,XC(TOP),YC(TOP),XC(K),YC(K),XC(K+1),
     $         YC(K+1),XU,YU,INTSCT)
          IF (INTSCT) THEN
             LR = LRLINE(XC(K+1),YC(K+1),XE,YE,XP,YP,0.0D0)
             IF (LR .EQ. 1) THEN
              OPER = 2
              CUR = K + 1
              RETURN
             ENDIF
          ENDIF
       ENDIF
   40    CONTINUE
       K = K + 1
      IF (K .LT. NV) GO TO 10
C
C     Error from unsuccessful scan.
C
      IERR = 209
      END
C
C     The following code was excerpted from: vpscnd.f
C
      SUBROUTINE VPSCND(XC,YC,IVIS)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IVIS(0:*)
      DOUBLE PRECISION XC(0:*),YC(0:*)
C
C     Written and copyright by:
C        Barry Joe, Dept. of Computing Science, Univ. of Alberta
C        Edmonton, Alberta, Canada  T6G 2H1
C        Phone: (403) 492-5757      Email: barry@cs.ualberta.ca
C
C     Purpose: This routine is called by routine VISPOL for the SCAND
C        operation (OPER = 6).
C
C     Input and updated parameters:
C        XC,YC,IVIS - see comments in routine VISPOL
C
C     Abnormal return:
C        IERR is set to 210
C
C     Routines called:
C        LRLINE, XEDGE
C
      INTEGER CUR,IERR,NV,OPER,TOP
      DOUBLE PRECISION XE,XW,YE,YW
      LOGICAL BEYE
      COMMON /GERROR/ IERR
      COMMON /GVPVAR/ NV,OPER,CUR,TOP,XE,YE,XW,YW,BEYE
      SAVE /GERROR/,/GVPVAR/
C
      INTEGER K,LR,LR1,LR2,LRLINE
      DOUBLE PRECISION XP,XU,YP,YU
      LOGICAL INTSCT
C
C     EYE-V(CUR-1)-V(CUR) is a right turn, S(TOP) is a V vertex not on
C     V(CUR-1)-V(CUR), TOP < CUR, W is intersection of V(CUR-1)-V(CUR)
C     and ray EYE-S(TOP), EYE-S(TOP)-W is a forward move, and
C     EYE-S(TOP-1)-S(TOP) is a left turn if TOP >= 1.
C     If BEYE, it is possible that (XW,YW) is a non-simple point,
C     but intersection at (XC(TOP),YC(TOP)) cannot occur.
C
      XP = XC(CUR-1)
      YP = YC(CUR-1)
      K = CUR
   10 CONTINUE
       CALL XEDGE(0,XW,YW,XC(TOP),YC(TOP),XC(K),YC(K),XC(K+1),
     $      YC(K+1),XU,YU,INTSCT)
       IF (INTSCT) THEN
          LR = LRLINE(XC(K+1),YC(K+1),XE,YE,XC(K),YC(K),0.0D0)
          LR1 = LRLINE(XC(K+1),YC(K+1),XE,YE,XC(TOP),YC(TOP),0.0D0)
          IF (LR .EQ. -1 .AND. LR1 .EQ. -1) THEN
             IF (XC(K) .NE. XW .OR. YC(K) .NE. YW) GO TO 20
             LR2 = LRLINE(XC(K+1),YC(K+1),XP,YP,XW,YW,0.0D0)
             IF (LR2 .EQ. -1) GO TO 30
   20          CONTINUE
             OPER = 1
             CUR = K + 1
             LR2 = LRLINE(XC(K),YC(K),XE,YE,XC(TOP),YC(TOP),0.0D0)
             TOP = TOP + 1
             IF (LR2 .EQ. 0) THEN
              XC(TOP) = XC(K)
              YC(TOP) = YC(K)
              IVIS(TOP) = K + NV
             ELSE
              XC(TOP) = XU
              YC(TOP) = YU
              IVIS(TOP) = -(K + 1 + NV)
             ENDIF
             TOP = TOP + 1
             XC(TOP) = XC(CUR)
             YC(TOP) = YC(CUR)
             IVIS(TOP) = CUR
             RETURN
          ENDIF
       ENDIF
   30    CONTINUE
       K = K + 1
      IF (K .LT. NV) GO TO 10
C
C     Error from unsuccessful scan.
C
      IERR = 210
      END

Generated by  Doxygen 1.6.0   Back to index