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

decomp2d.f

C
C     The following code was excerpted from: cvdec2.f
C
      SUBROUTINE CVDEC2(ANGSPC,ANGTOL,NVC,NPOLG,NVERT,MAXVC,MAXHV,MAXPV,
     $   MAXIW,MAXWK,VCL,REGNUM,HVL,PVL,IANG,IWK,WK)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXHV,MAXIW,MAXPV,MAXVC,MAXWK,NPOLG,NVC,NVERT
      INTEGER HVL(MAXHV),IWK(MAXIW),PVL(4,MAXPV),REGNUM(MAXHV)
      DOUBLE PRECISION ANGSPC,ANGTOL,IANG(MAXPV),VCL(2,MAXVC),WK(MAXWK)
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: Decompose general polygonal region (which is decomposed
C        into simple polygons on input) into convex polygons using
C        vertex coordinate list, head vertex list, and polygon vertex
C        list data structures.
C
C     Input parameters:
C        ANGSPC - angle spacing parameter in radians used in controlling
C              vertices to be considered as an endpoint of a separator
C        ANGTOL - angle tolerance parameter in radians used in accepting
C              separator(s)
C        NVC - number of vertex coordinates or positions used in VCL
C              array
C        NPOLG - number of polygonal subregions or positions used in
C              HVL array
C        NVERT - number of polygon vertices or positions used in PVL
C              array
C        MAXVC - maximum size available for VCL array, should be >=
C              number of vertex coordinates required for decomposition 
C        MAXHV - maximum size available for HVL, REGNUM arrays, should
C              be >= number of polygons required for decomposition
C        MAXPV - maximum size available for PVL, IANG arrays; should be
C              >= number of polygon vertices required for decomposition
C        MAXIW - maximum size available for IWK array; should be about
C              3 times maximum number of vertices in any polygon
C        MAXWK - maximum size available for WK array; should be about
C              5 times maximum number of vertices in any polygon
C        VCL(1:2,1:NVC) - vertex coordinate list
C        REGNUM(1:NPOLG) - region numbers
C        HVL(1:NPOLG) - head vertex list
C        PVL(1:4,1:NVERT),IANG(1:NVERT) - polygon vertex list and
C              interior angles; see routine DSPGDC for more details
C        [Note: The data structures should be as output from routine
C              SPDEC2.]
C
C     Updated parameters:
C        NVC,NPOLG,NVERT,VCL,REGNUM,HVL,PVL,IANG
C
C     Working parameters:
C        IWK(1:MAXIW) - integer work array
C        WK(1:MAXWK) - double precision work array
C
C     Abnormal return:
C        IERR is set to 3, 4, 5, 6, 7, 206, 207, 208, 209, 210, or 212
C
C     Routines called:
C        INSED2, RESVRT
C
      INTEGER IERR
      DOUBLE PRECISION PI,TOL
      COMMON /GERROR/ IERR
      COMMON /GCONST/ PI,TOL
      SAVE /GERROR/,/GCONST/
C
      INTEGER V,W1,W2
      DOUBLE PRECISION PIPTOL
C
C     For each reflex vertex, resolve it with one or two separators
C     and update VCL, HVL, PVL, IANG.
C
      PIPTOL = PI + TOL
      V = 1
   10 CONTINUE
      IF (V .GT. NVERT) RETURN
       IF (IANG(V) .GT. PIPTOL) THEN
            CALL RESVRT(V,ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,MAXIW,
     $         MAXWK,VCL,PVL,IANG,W1,W2,IWK,WK)
          IF (IERR .NE. 0) RETURN
          CALL INSED2(V,W1,NPOLG,NVERT,MAXHV,MAXPV,VCL,REGNUM,HVL,
     $         PVL,IANG)
          IF (IERR .NE. 0) RETURN
            IF (W2 .GT. 0) CALL INSED2(V,W2,NPOLG,NVERT,MAXHV,MAXPV,
     $         VCL,REGNUM,HVL,PVL,IANG)
          IF (IERR .NE. 0) RETURN
         ENDIF
      V = V + 1
      GO TO 10
      END
C
C     The following code was excerpted from: dsmcpr.f
C
      SUBROUTINE DSMCPR(NHOLE,NVBC,VCL,MAXHV,MAXPV,MAXHO,NVC,NPOLG,
     $   NVERT,NHOLA,REGNUM,HVL,PVL,IANG,HOLV)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXHO,MAXHV,MAXPV,NHOLA,NHOLE,NPOLG,NVC,NVERT
      INTEGER HVL(NHOLE+1),HOLV(MAXHO),NVBC(NHOLE+1),PVL(4,MAXPV)
      INTEGER REGNUM(1)
      DOUBLE PRECISION IANG(MAXPV),VCL(2,*)
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: Initialize the polygonal decomposition data structure
C        given a multiply-connected polygonal region with 1 outer
C        boundary curve and 0 or more inner boundary curves of holes.
C
C     Input parameters:
C        NHOLE - number of holes in region
C        NVBC(1:NHOLE+1) - number of vertices per boundary curve; first
C              boundary curve is the outer boundary of the region
C        VCL(1:2,1:NVC) - vertex coordinates of boundary curves in CCW
C              order; NVC = NVBC(1) + ... + NVBC(NHOLE+1); positions 1
C              to NVBC(1) of VCL contain the vertex coordinates of the
C              outer boundary in CCW order; positions NVBC(1)+1 to
C              NVBC(1)+NVBC(2) contain the vertex coordinates of the
C              first hole boundary in CCW order, etc.
C        MAXHV - maximum size available for HVL, REGNUM arrays, should
C              be >= NHOLE + 1
C        MAXPV - maximum size available for PVL, IANG arrays; should be
C               >= NVC
C        MAXHO - maximum size available for HOLV array; should be
C               >= NHOLE*2
C
C     Output parameters:
C        NVC - number of vertex coordinates, set to sum of NVBC(I)
C        NPOLG - number of polygonal subregions, set to 1
C        NVERT - number of vertices in PVL, set to NVC
C        NHOLA - number of attached holes, set to 0
C        REGNUM(1:1) - region number of only subregion, set to 1
C        [Note: Above 4 parameters are for consistency with DSPGDC.]
C        HVL(1:NHOLE+1) - head vertex list; first entry is the head
C              vertex (index in PVL) of outer boundary curve; next
C              NHOLE entries contain the head vertex of a hole
C        PVL(1:4,1:NVC),IANG(1:NVC) - polygon vertex list and interior
C              angles; vertices of outer boundary curve are in CCW order
C              followed by vertices of each hole in CW hole; vertices
C              of each polygon are in a circular linked list; see
C              routine DSPGDC for more details of this data structure
C        HOLV(1:NHOLE*2) - indices in PVL of top and bottom vertices of
C              holes; first (last) NHOLE entries are for top (bottom)
C              vertices; top (bottom) vertices are sorted in decreasing
C              (increasing) lexicographic (y,x) order of coordinates
C
C     Abnormal return:
C        IERR is set to 2, 4, or 5
C
C     Routines called:
C        ANGLE, HOLVRT
C
      INTEGER IERR
      COMMON /GERROR/ IERR
      SAVE /GERROR/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER I,IV,IVS,J,LV,LVP,LVS,NV,NVS
      DOUBLE PRECISION ANGLE
C
      NVC = 0
      DO 10 I = 1,NHOLE+1
       NVC = NVC + NVBC(I)
   10 CONTINUE
      NPOLG = 1
      NVERT = NVC
      NHOLA = 0
      REGNUM(1) = 1
      IF (NHOLE + 1 .GT. MAXHV) THEN
       IERR = 4
       RETURN
      ELSE IF (NVC .GT. MAXPV) THEN
       IERR = 5
       RETURN
      ELSE IF (NHOLE + NHOLE .GT. MAXHO) THEN
       IERR = 2
       RETURN
      ENDIF
C
C     Initialize HVL, PVL arrays.
C
   20 CONTINUE
      HVL(1) = 1
      NV = NVBC(1)
      DO 30 I = 1,NV
       PVL(LOC,I) = I
       PVL(POLG,I) = 1
       PVL(SUCC,I) = I + 1
       PVL(EDGV,I) = 0
   30 CONTINUE
      PVL(SUCC,NV) = 1
      DO 50 J = 1,NHOLE
       HVL(J+1) = NV + 1
       NVS = NV + NVBC(J+1)
         DO 40 I = NV+1,NVS
          PVL(LOC,I) = I
          PVL(POLG,I) = 1
          PVL(SUCC,I) = I - 1
          PVL(EDGV,I) = 0
   40    CONTINUE
         PVL(SUCC,NV+1) = NVS
       NV = NVS
   50 CONTINUE
C
C     Initialize IANG array.
C
      DO 70 I = 1,NHOLE+1
       J = HVL(I)
       LVP = PVL(LOC,J)
       IV = PVL(SUCC,J)
       LV = PVL(LOC,IV)
   60    CONTINUE
          IVS = PVL(SUCC,IV)
          LVS = PVL(LOC,IVS)
          IANG(IV) = ANGLE(VCL(1,LVP),VCL(2,LVP),VCL(1,LV),VCL(2,LV),
     $         VCL(1,LVS),VCL(2,LVS))
          IF (IV .EQ. J) GO TO 70
          LVP = LV
          IV = IVS
          LV = LVS
       GO TO 60
   70 CONTINUE
C
C     Initialize HOLV array.
C
      IF (NHOLE .GT. 0) CALL HOLVRT(NHOLE,VCL,HVL(2),PVL,HOLV)
      END
C
C     The following code was excerpted from: edght.f
C
      SUBROUTINE EDGHT(A,B,V,N,HTSIZ,MAXEDG,HDFREE,LAST,HT,EDGE,W)
      IMPLICIT LOGICAL (A-Z)
      INTEGER A,B,HDFREE,HTSIZ,LAST,MAXEDG,N,V,W
      INTEGER EDGE(4,MAXEDG),HT(0:HTSIZ-1)
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: Search in hash table HT for record in EDGE containing
C        key (A,B).
C
C     Input parameters:
C        A,B - vertex indices, > 0, of edge (also key of hash table)
C        V - value associated with edge
C        N - upper bound on A, B
C        HTSIZ - size of hash table HT
C        MAXEDG - maximum size available for EDGE array
C        HDFREE - head pointer to linked list of free entries of EDGE
C              array due to deletions
C        LAST - index of last entry used in EDGE array
C        HT(0:HTSIZ-1) - hash table of head pointers (direct chaining
C              with ordered lists is used)
C        EDGE(1:4,1:MAXEDG) - entries of hash table records;
C              EDGE(1,I) = MIN(A,B); EDGE(2,I) = MAX(A,B);
C              EDGE(3,I) = V; EDGE(4,I) = link
C        [Note: Before first call to this routine, HDFREE, LAST, and
C              entries of HT should be set to 0.]
C
C     Updated parameters:
C        HDFREE,LAST - at least one of these is updated
C        HT,EDGE - if key with A,B is found then this record is deleted
C              from hash table, else record is inserted in hash table
C
C     Output parameters:
C        W - EDGE(3,INDEX), where INDEX is index of record, if found;
C              else 0
C
C     Abnormal return:
C        IERR is set to 1
C
      INTEGER IERR
      COMMON /GERROR/ IERR
      SAVE /GERROR/
C
      INTEGER AA,BB,BPTR,K,NEWP,PTR
C
      IF (A .LT. B) THEN
         AA = A
         BB = B
      ELSE
         AA = B
         BB = A
      ENDIF
      K = MOD(AA*N + BB, HTSIZ)
      BPTR = -1
      PTR = HT(K)
   10 CONTINUE
      IF (PTR .NE. 0) THEN
       IF (EDGE(1,PTR) .GT. AA) THEN
          GO TO 20
       ELSE IF (EDGE(1,PTR) .EQ. AA) THEN
          IF (EDGE(2,PTR) .GT. BB) THEN
             GO TO 20
          ELSE IF (EDGE(2,PTR) .EQ. BB) THEN
             IF (BPTR .EQ. -1) THEN
              HT(K) = EDGE(4,PTR)
             ELSE
              EDGE(4,BPTR) = EDGE(4,PTR)
             ENDIF
             EDGE(4,PTR) = HDFREE
             HDFREE = PTR
             W = EDGE(3,PTR)
             RETURN
          ENDIF
       ENDIF
       BPTR = PTR
       PTR = EDGE(4,PTR)
       GO TO 10
      ENDIF
   20 CONTINUE
      IF (HDFREE .GT. 0) THEN
       NEWP = HDFREE
       HDFREE = EDGE(4,HDFREE)
      ELSE
       LAST = LAST + 1
       NEWP = LAST
       IF (LAST .GT. MAXEDG) THEN
          IERR = 1
          RETURN
       ENDIF
      ENDIF
      IF (BPTR .EQ. -1) THEN
       HT(K) = NEWP
      ELSE
       EDGE(4,BPTR) = NEWP
      ENDIF
      EDGE(1,NEWP) = AA
      EDGE(2,NEWP) = BB
      EDGE(3,NEWP) = V
      EDGE(4,NEWP) = PTR
      W = 0
      END
C
C     The following code was excerpted from: fndsep.f
C
      SUBROUTINE FNDSEP(ANGAC1,XR,YR,NVRT,XC,YC,IVIS,THETA,NV,IV,
     $   VCL,PVL,IANG,ANGSEP,I1,I2,WKANG)
      IMPLICIT LOGICAL (A-Z)
      INTEGER I1,I2,NV,NVRT
      INTEGER IV(0:NV),IVIS(0:NVRT),PVL(4,*)
      DOUBLE PRECISION ANGAC1,ANGSEP,XR,YR,IANG(*),THETA(0:NVRT)
      DOUBLE PRECISION VCL(2,*),WKANG(0:NV),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: Find 1 or 2 separators which can resolve reflex vertex
C        (XR,YR) using a max-min angle criterion from list of vertices
C        in increasing polar angle w.r.t. reflex vertex. Preference
C        is given to 1 separator.
C
C     Input parameters:
C        ANGAC1 - angle tolerance parameter used for preference
C              in accepting one separator
C        XR,YR - coordinates of reflex vertex
C        NVRT - (number of vertices) - 1
C        XC(0:NVRT), YC(0:NVRT) - vertex coordinates of possible
C              endpoints of a separator
C        IVIS(0:NVRT) - contains information about the vertices of
C              XC, YC arrays w.r.t. the polygon vertex list; if
C              IVIS(I) > 0 then vertex (XC(I),YC(I)) has index IVIS(I)
C              in PVL; if IVIS(I) < 0 then vertex (XC(I),YC(I)) is on
C              the edge joining vertices with indices -IVIS(I) and
C              SUCC(-IVIS(I)) in PVL
C        THETA(0:NVRT) - polar angles of vertices in increasing order;
C              THETA(NVRT) is the interior angle of reflex vertex;
C              THETA(I), I >= 0, is the polar angle of (XC(I),YC(I))
C              w.r.t. reflex vertex
C        NV - (number of vertices to be considered as endpoint of a
C              separator) - 1
C        IV(0:NV) - indices of vertices in XC, YC arrays to be
C              considered as endpoint of a separator; angle between
C              consecutive vertices is assumed to be < 180 degrees
C        VCL(1:2,1:*) - vertex coordinate list
C        PVL(1:4,1:*),IANG(1:*) - polygon vertex list, interior angles
C
C     Output parameters:
C        ANGSEP - minimum of the 4 or 7 angles at the boundary
C              resulting from 1 or 2 separators, respectively
C        I1,I2 - indices of endpoints of separators in XC, YC arrays;
C              I2 = -1 if there is only one separator, else I1 < I2
C
C     Working parameters:
C        WKANG(0:NV) - working array for angles
C
C     Routines called:
C        MINANG
C
      DOUBLE PRECISION PI,TOL
      COMMON /GCONST/ PI,TOL
      SAVE /GCONST/
C
      INTEGER I,II,K,L,M,NL,NR,P,Q,R
      DOUBLE PRECISION ANG,ANGSP2,MINANG,PHI
C
C     Determine the vertices in the inner cone - indices P to Q.
C
      I = 0
      P = -1
      PHI = THETA(NVRT) - PI + TOL
   10 CONTINUE
      IF (P .GE. 0) GO TO 20
       IF (THETA(IV(I)) .GE. PHI) THEN
          P = I
       ELSE
          I = I + 1
       ENDIF
      GO TO 10
   20 CONTINUE
      I = NV
      Q = -1
      PHI = PI - TOL
   30 CONTINUE
      IF (Q .GE. 0) GO TO 40
       IF (THETA(IV(I)) .LE. PHI) THEN
          Q = I
       ELSE
          I = I - 1
       ENDIF
      GO TO 30
   40 CONTINUE
C
C     Use the max-min angle criterion to find the best separator
C     in inner cone.
C
      ANGSEP = 0.0D0
      DO 50 I = P,Q
       K = IV(I)
       ANG = MINANG(XR,YR,XC(K),YC(K),IVIS(K),THETA(K),THETA(NVRT),
     $      VCL,PVL,IANG)
       IF (ANG .GT. ANGSEP) THEN
          ANGSEP = ANG
          II = IV(I)
       ENDIF
   50 CONTINUE
      ANGSP2 = ANGSEP
      IF (ANGSEP .GE. ANGAC1) GO TO 110
C
C     If the best separator in inner cone is not 'good' enough,
C     use max-min angle criterion to try to find a better pair
C     of separators from the right and left cones.
C
      NR = 0
      NL = 0
      DO 60 R = 0,P-1
         WKANG(R) = 0.0D0
         IF (THETA(IV(R)) .GT. ANGSEP) THEN
          K = IV(R)
          ANG = MINANG(XR,YR,XC(K),YC(K),IVIS(K),THETA(K),THETA(NVRT),
     $         VCL,PVL,IANG)
          IF (ANG .GT. ANGSEP) THEN
             NR = NR + 1
             WKANG(R) = ANG
          ENDIF
         ENDIF
   60 CONTINUE
      IF (NR .EQ. 0) GO TO 110
      PHI = THETA(NVRT) - ANGSEP
      DO 70 L = Q+1,NV
         WKANG(L) = 0.0D0
         IF (THETA(IV(L)) .LT. PHI) THEN
          K = IV(L)
          ANG = MINANG(XR,YR,XC(K),YC(K),IVIS(K),THETA(K),THETA(NVRT),
     $         VCL,PVL,IANG)
          IF (ANG .GT. ANGSEP) THEN
             NL = NL + 1
             WKANG(L) = ANG
          ENDIF
         ENDIF
   70 CONTINUE
      IF (NL .EQ. 0) GO TO 110
C
C     Check all possible pairs for the best pair of separators
C     in the right and left cones.
C
      M = NV
      DO 100 R = P-1,0,-1
         IF (M .GT. Q .AND. WKANG(R) .GT. ANGSP2) THEN
          PHI = THETA(IV(R))
   80       CONTINUE
          IF (M .GT. Q .AND. (WKANG(M) .LE. ANGSP2 .OR.
     $         THETA(IV(M)) - PHI .GT. PI - TOL)) THEN
             M = M - 1
             GO TO 80
          ENDIF
          DO 90 L = Q+1,M
             IF (WKANG(L) .GT. ANGSP2) THEN
                ANG = MIN(THETA(IV(L)) - PHI, WKANG(R), WKANG(L))
                IF (ANG .GT. ANGSP2) THEN
                 ANGSP2 = ANG
                 I1 = IV(R)
                 I2 = IV(L)
                ENDIF
             ENDIF
   90       CONTINUE
         ENDIF
  100 CONTINUE
C
C     Choose 1 or 2 separators based on max-min angle criterion or
C     ANGAC1 parameter.
C
  110 CONTINUE
      IF (ANGSP2 .LE. ANGSEP) THEN
       I1 = II
       I2 = -1
      ELSE
       ANGSEP = ANGSP2
      ENDIF
      END
C
C     The following code was excerpted from: holvrt.f
C
      SUBROUTINE HOLVRT(NHOLE,VCL,HVL,PVL,HOLV)
      IMPLICIT LOGICAL (A-Z)
      INTEGER NHOLE
      INTEGER HOLV(NHOLE*2),HVL(NHOLE),PVL(4,*)
      DOUBLE PRECISION VCL(2,*)
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 top and bottom vertices of holes in polygonal
C        region(s), and sort top vertices in decreasing (y,x) order
C        and bottom vertices in increasing (y,x) order.
C
C     Input parameters:
C        NHOLE - number of holes in region(s)
C        VCL(1:2,1:*) - vertex coordinate list
C        HVL(1:NHOLE) - head vertex list; HVL(I) is index in PVL of
C              head vertex of Ith hole
C        PVL(1:4,1:*) - polygon vertex list; see routine DSPGDC
C
C     Output parameters:
C        HOLV(1:NHOLE*2) - indices in PVL of top and bottom vertices of
C              holes; first (last) NHOLE entries are for top (bottom)
C              vertices; top (bottom) vertices are sorted in decreasing
C              (increasing) lexicographic (y,x) order of coordinates
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER HV,I,IMAX,IMIN,IV,J,LV,NHP1
      DOUBLE PRECISION X,XMAX,XMIN,Y,YMAX,YMIN
C
C     Determine top and bottom vertices of holes.
C
      IMIN = HVL(1)
      IMAX = HVL(1)
      XMIN = VCL(1,PVL(LOC,HVL(1)))
      YMIN = VCL(2,PVL(LOC,HVL(1)))
      XMAX = XMIN
      YMAX = YMIN
      DO 20 I = 1,NHOLE
       HV = HVL(I)
       IV = HV
   10    CONTINUE
          LV = PVL(LOC,IV)
          IF (IV .EQ. HV) THEN
             IMIN = IV
             IMAX = IV
             XMIN = VCL(1,LV)
             YMIN = VCL(2,LV)
             XMAX = XMIN
             YMAX = YMIN
          ELSE
             X = VCL(1,LV)
             Y = VCL(2,LV)
             IF (Y .LT. YMIN .OR. Y .EQ. YMIN .AND. X .LT. XMIN) THEN
              IMIN = IV
              XMIN = X
              YMIN = Y
             ELSE IF (Y .GT. YMAX .OR. Y .EQ. YMAX .AND. X .GT. XMAX)
     $         THEN
              IMAX = IV
              XMAX = X
              YMAX = Y
             ENDIF
          ENDIF
          IV = PVL(SUCC,IV)
       IF (IV .NE. HV) GO TO 10
       HOLV(I) = IMAX
       HOLV(I+NHOLE) = IMIN
   20 CONTINUE
C
C     Use linear insertion sort to sort the top vertices of holes
C     in decreasing (y,x) order, then bottom vertices in increasing
C     (y,x) order. It is assumed NHOLE is small.
C
      DO 40 I = 2,NHOLE
       HV = HOLV(I)
       LV = PVL(LOC,HV)
       X = VCL(1,LV)
       Y = VCL(2,LV)
       J = I
   30    CONTINUE
          IV = HOLV(J-1)
          LV = PVL(LOC,IV)
          IF (Y .GT. VCL(2,LV) .OR. Y .EQ. VCL(2,LV) .AND.
     $         X .GT. VCL(1,LV)) THEN
             HOLV(J) = IV
             J = J - 1
             IF (J .GT. 1) GO TO 30
          ENDIF
       HOLV(J) = HV
   40 CONTINUE
C
      NHP1 = NHOLE + 1
      DO 60 I = NHP1+1,NHOLE+NHOLE
       HV = HOLV(I)
       LV = PVL(LOC,HV)
       X = VCL(1,LV)
       Y = VCL(2,LV)
       J = I
   50    CONTINUE
          IV = HOLV(J-1)
          LV = PVL(LOC,IV)
          IF (Y .LT. VCL(2,LV) .OR. Y .EQ. VCL(2,LV) .AND.
     $         X .LT. VCL(1,LV)) THEN
             HOLV(J) = IV
             J = J - 1
             IF (J .GT. NHP1) GO TO 50
          ENDIF
       HOLV(J) = HV
   60 CONTINUE
      END
C
C     The following code was excerpted from: insed2.f
C
      SUBROUTINE INSED2(V,W,NPOLG,NVERT,MAXHV,MAXPV,VCL,REGNUM,HVL,
     $   PVL,IANG)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXHV,MAXPV,NPOLG,NVERT,V,W
      INTEGER HVL(MAXHV),PVL(4,MAXPV)
      INTEGER REGNUM(MAXHV)
      DOUBLE PRECISION IANG(MAXPV),VCL(2,*)
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: Insert edge joining vertices V, W into head vertex
C        list and polygon vertex list data structures.
C
C     Input parameters:
C        V,W - indices in PVL of vertices which are the endpoints
C              of an edge to be added to decomposition
C        NPOLG - number of positions used in HVL array
C        NVERT - number of positions used in PVL array
C        MAXHV - maximum size available for HVL array
C        MAXPV - maximum size available for PVL array
C        VCL(1:2,1:*) - vertex coordinate list
C        REGNUM(1:NPOLG) - region numbers
C        HVL(1:NPOLG) - head vertex list
C        PVL(1:4,1:NVERT),IANG(1:NVERT) - polygon vertex list and
C              interior angles
C
C     Updated parameters:
C        REGNUM,HVL,PVL,IANG
C
C     Abnormal return:
C        IERR is set to 4 or 5
C
C     Routines called:
C        ANGLE
C
      INTEGER IERR,IPRT,MSGLVL
      COMMON /GERROR/ IERR
      COMMON /GPRINT/ IPRT,MSGLVL
      SAVE /GERROR/,/GPRINT/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER I,L,LV,LW,VV,WW
      DOUBLE PRECISION ANGLE
C
      IF (NPOLG .GE. MAXHV) THEN
       IERR = 4
       RETURN
      ELSE IF (NVERT+2 .GT. MAXPV) THEN
       IERR = 5
       RETURN
      ENDIF
C
C     Split linked list of vertices of the polygon containing vertices
C     V and W into two linked list of vertices of polygons with common
C     edge joining V and W.
C
      NVERT = NVERT + 2
      VV = NVERT - 1
      WW = NVERT
      LV = PVL(LOC,V)
      LW = PVL(LOC,W)
      PVL(LOC,VV) = LV
      PVL(LOC,WW) = LW
      PVL(POLG,WW) = PVL(POLG,V)
      PVL(SUCC,VV) = PVL(SUCC,V)
      PVL(SUCC,WW) = PVL(SUCC,W)
      PVL(SUCC,V) = WW
      PVL(SUCC,W) = VV
      PVL(EDGV,VV) = PVL(EDGV,V)
      PVL(EDGV,WW) = PVL(EDGV,W)
      PVL(EDGV,V) = W
      PVL(EDGV,W) = V
      IF (PVL(EDGV,VV) .GT. 0) PVL(EDGV,PVL(EDGV,VV)) = VV
      IF (PVL(EDGV,WW) .GT. 0) PVL(EDGV,PVL(EDGV,WW)) = WW
      L = PVL(LOC,PVL(SUCC,VV))
      IANG(VV) = ANGLE(VCL(1,LW),VCL(2,LW),VCL(1,LV),VCL(2,LV),
     $   VCL(1,L),VCL(2,L))
      IANG(V) = IANG(V) - IANG(VV)
      L = PVL(LOC,PVL(SUCC,WW))
      IANG(WW) = ANGLE(VCL(1,LV),VCL(2,LV),VCL(1,LW),VCL(2,LW),
     $   VCL(1,L),VCL(2,L))
      IANG(W) = IANG(W) - IANG(WW)
      NPOLG = NPOLG + 1
      I = VV
   10 CONTINUE
       PVL(POLG,I) = NPOLG
       I = PVL(SUCC,I)
      IF (I .NE. VV) GO TO 10
      HVL(PVL(POLG,V)) = V
      HVL(NPOLG) = VV
      REGNUM(NPOLG) = REGNUM(PVL(POLG,V))
C
      END
C
C     The following code was excerpted from: insvr2.f
C
      SUBROUTINE INSVR2(XI,YI,WP,NVC,NVERT,MAXVC,MAXPV,VCL,PVL,IANG,W)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXPV,MAXVC,NVC,NVERT,PVL(4,MAXPV),W,WP
      DOUBLE PRECISION IANG(MAXPV),VCL(2,MAXVC),XI,YI
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: Insert point (XI,YI) into vertex coordinate list and
C        polygon vertex list data structures.
C
C     Input parameters:
C        XI,YI - coordinates of point to be inserted
C        WP - index of vertex in PVL which is to be the predecessor
C              vertex of the inserted vertex
C        NVC - number of positions used in VCL array
C        NVERT - number of positions used in PVL array
C        MAXVC - maximum size available for VCL array
C        MAXPV - maximum size available for PVL array
C        VCL(1:2,1:NVC) - vertex coordinate list
C        PVL(1:4,1:NVERT),IANG(1:NVERT) - polygon vertex list and
C              interior angles
C
C     Updated parameters:
C        NVC,NVERT,VCL,PVL,IANG
C
C     Output parameter:
C        W - index of inserted vertex in PVL
C
C     Abnormal return:
C        IERR is set to 3 or 5
C
      INTEGER IERR
      DOUBLE PRECISION PI,TOL
      COMMON /GERROR/ IERR
      COMMON /GCONST/ PI,TOL
      SAVE /GERROR/,/GCONST/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER WS,WW,WWP,WWS
C
      IF (NVC .GE. MAXVC) THEN
       IERR = 3
       RETURN
      ELSE IF (NVERT+2 .GT. MAXPV) THEN
       IERR = 5
       RETURN
      ENDIF
C
C     Update linked list of vertices of polygon containing vertex WP.
C
      NVC = NVC + 1
      VCL(1,NVC) = XI
      VCL(2,NVC) = YI
      WS = PVL(SUCC,WP)
      NVERT = NVERT + 1
      W = NVERT
      PVL(LOC,W) = NVC
      PVL(POLG,W) = PVL(POLG,WP)
      PVL(SUCC,WP) = W
      PVL(SUCC,W) = WS
      IANG(W) = PI
      PVL(EDGV,W) = PVL(EDGV,WP)
C
C     If edge containing (XI,YI) is shared by another polygon,
C     then also update linked list of vertices of that polygon.
C
      IF (PVL(EDGV,WP) .GT. 0) THEN
       WWS = PVL(EDGV,WP)
         WWP = PVL(SUCC,WWS)
         NVERT = NVERT + 1
         WW = NVERT
         PVL(LOC,WW) = NVC
         PVL(POLG,WW) = PVL(POLG,WWS)
         PVL(SUCC,WWS) = WW
         PVL(SUCC,WW) = WWP
         IANG(WW) = PI
         PVL(EDGV,WP) = WW
       PVL(EDGV,WW) = WP
       PVL(EDGV,WWS) = W
      ENDIF
      END
C
C     The following code was excerpted from: jnhole.f
C
      SUBROUTINE JNHOLE(ITOPHV,ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,
     $   MAXIW,MAXWK,VCL,HVL,PVL,IANG,IWK,WK)
      IMPLICIT LOGICAL (A-Z)
      INTEGER ITOPHV,MAXIW,MAXPV,MAXVC,MAXWK,NVC,NVERT
      INTEGER HVL(*),IWK(MAXIW),PVL(4,MAXPV)
      DOUBLE PRECISION ANGSPC,ANGTOL,IANG(MAXPV),VCL(2,MAXVC),WK(MAXWK)
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: Join hole boundary to boundary of polygon containing hole
C        by finding a cut edge originating from the top vertex of hole
C        to a point on outer polygon boundary above it.
C
C     Input parameters:
C        ITOPHV - index in PVL of top vertex of hole
C        ANGSPC - angle spacing parameter used in controlling the
C              vertices to be considered as an endpoint of a separator
C        ANGTOL - angle tolerance parameter used in accepting
C              separator(s)
C        NVC - number of positions used in VCL array
C        NVERT - number of positions used in PVL array
C        MAXVC - maximum size available for VCL array
C        MAXPV - maximum size available for PVL array
C        MAXIW - maximum size available for IWK array; should be about
C              3 times number of vertices in outer polygon
C        MAXWK - maximum size available for WK array; should be about
C              5 times number of vertices in outer polygon
C        VCL(1:2,1:NVC) - vertex coordinate list
C        HVL(1:*) - head vertex list
C        PVL(1:4,1:NVERT),IANG(1:NVERT) - polygon vertex list and
C              interior angles
C
C     Updated parameters:
C        NVC,NVERT,VCL,PVL,IANG
C
C     Working parameters:
C        IWK(1:MAXIW) - integer work array
C        WK(1:MAXWK) - double precision work array
C
C     Abnormal return:
C        IERR is set to 3, 5, 6, 7, 206 to 210, 212, 218, or 219
C
C     Routines called:
C        ANGLE, RESVRT
C
      INTEGER IERR,IPRT,MSGLVL
      DOUBLE PRECISION PI,TOL
      COMMON /GERROR/ IERR
      COMMON /GCONST/ PI,TOL
      COMMON /GPRINT/ IPRT,MSGLVL
      SAVE /GERROR/,/GCONST/,/GPRINT/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER HV,ILFT,IPOLY,IRGT,IV,IVS,L,LV,LW,SUCCIL,SUCCIR
      INTEGER VP,VR,VS,VV,W,WW
      DOUBLE PRECISION ANGLE,DY,S,SLFT,SRGT,XINT,XLFT,XRGT,XT,XV,XVS
      DOUBLE PRECISION YLFT,YRGT,YT,YTMTOL,YTPTOL,YV,YVS
C
      IF (NVC+3 .GT. MAXVC) THEN
       IERR = 3
       RETURN
      ELSE IF (NVERT+5 .GT. MAXPV) THEN
       IERR = 5
       RETURN
      ENDIF
C
C     Determine 'closest' vertices on outer boundary which are to the
C     left and right of top vertex of hole and on the horizontal line
C     through top vertex. The two closest vertices must be on edges
C     which intersect the horizontal line and are partially above the
C     line. Ties are broken (in the case of a vertex on a cut edge)
C     by choosing the vertex on the edge of maximum or minimum dx/dy
C     slope depending on whether the vertex is to the left or right
C     of top vertex, respectively.
C
      IPOLY = PVL(POLG,ITOPHV)
      LV = PVL(LOC,ITOPHV)
      XT = VCL(1,LV)
      YT = VCL(2,LV)
      DY = 0.0D0
      HV = HVL(IPOLY)
      IV = HV
      YV = VCL(2,PVL(LOC,IV))
      SLFT = 0.0
      SRGT = 0.0
   10 CONTINUE
       IV = PVL(SUCC,IV)
       YVS = VCL(2,PVL(LOC,IV))
       DY = MAX(DY,ABS(YVS-YV))
       YV = YVS
      IF (IV .NE. HV) GO TO 10
      YTMTOL = YT - TOL*DY
      YTPTOL = YT + TOL*DY
      ILFT = 0
      IRGT = 0
      XLFT = 0.0D0
      XRGT = 0.0D0
      HV = HVL(IPOLY)
      IV = HV
      LV = PVL(LOC,IV)
      XV = VCL(1,LV)
      YV = VCL(2,LV)
   20 CONTINUE
       IVS = PVL(SUCC,IV)
       LV = PVL(LOC,IVS)
       XVS = VCL(1,LV)
       YVS = VCL(2,LV)
       IF (YV .LE. YTPTOL .AND. YVS .GT. YTPTOL) THEN
          IF (YV .GE. YTMTOL) THEN
             IF (XV .GT. XT) THEN
              IF (XV .LT. XRGT .OR. IRGT .EQ. 0) THEN
                 IRGT = IV
                 XRGT = XV
                 YRGT = YV
                 SRGT = (XVS - XV)/(YVS - YV)
              ELSE IF (XV .EQ. XRGT) THEN
                 S = (XVS - XV)/(YVS - YV)
                 IF (S .LT. SRGT) THEN
                  IRGT = IV
                    YRGT = YV
                  SRGT = S
                 ENDIF
              ENDIF
             ENDIF
          ELSE
             XINT = (YT - YV)*(XVS - XV)/(YVS - YV) + XV
             IF (XINT .GT. XT) THEN
              IF (XINT .LT. XRGT .OR. IRGT .EQ. 0) THEN
                 IRGT = IV
                 XRGT = XINT
                 YRGT = YT
              ENDIF
             ENDIF
          ENDIF
       ELSE IF (YV .GT. YTPTOL .AND. YVS .LE. YTPTOL) THEN
          IF (YVS .GE. YTMTOL) THEN
             IF (XVS .LT. XT) THEN
              IF (XVS .GT. XLFT .OR. ILFT .EQ. 0) THEN
                 ILFT = IV
                 XLFT = XVS
                 YLFT = YVS
                 SLFT = (XVS - XV)/(YVS - YV)
              ELSE IF (XVS .EQ. XLFT) THEN
                 S = (XVS - XV)/(YVS - YV)
                 IF (S .GT. SLFT) THEN
                  ILFT = IV
                    YLFT = YVS
                  SLFT = S
                 ENDIF
              ENDIF
             ENDIF
          ELSE
             XINT = (YT - YV)*(XVS - XV)/(YVS - YV) + XV
             IF (XINT .LT. XT) THEN
              IF (XINT .GT. XLFT .OR. ILFT .EQ. 0) THEN
                 ILFT = IV
                 XLFT = XINT
                 YLFT = YT
              ENDIF
             ENDIF
          ENDIF
       ENDIF
       IV = IVS
       XV = XVS
       YV = YVS
      IF (IV .NE. HV) GO TO 20
C
      IF (ILFT .EQ. 0 .OR. IRGT .EQ. 0) THEN
       IERR = 218
       RETURN
      ENDIF
C
C     Temporarily modify PVL to pass the subregion 'above' top vertex
C     of hole to routine RESVRT. The top vertex is the reflex vertex
C     passed to RESVRT (in the temporary subregion, it has interior
C     angle PI). This causes one separator to be chosen by RESVRT
C     and its other endpoint is above the top vertex.
C
      SUCCIL = PVL(SUCC,ILFT)
      SUCCIR = PVL(SUCC,IRGT)
      VCL(1,NVC+2) = XLFT
      VCL(2,NVC+2) = YLFT
      VCL(1,NVC+3) = XRGT
      VCL(2,NVC+3) = YRGT
      VP = NVERT + 3
      VR = NVERT + 4
      VS = NVERT + 5
      IANG(VR) = ANGLE(XLFT,YLFT,XT,YT,XRGT,YRGT)
      IF (IANG(VR) .LT. PI - TOL .OR. IANG(VR) .GT. PI + TOL) THEN
       IERR = 219
       RETURN
      ENDIF
      PVL(LOC,VP) = NVC + 2
      PVL(POLG,VP) = IPOLY
      PVL(SUCC,VP) = VR
      PVL(EDGV,VP) = 0
      PVL(LOC,VR) = PVL(LOC,ITOPHV)
      PVL(POLG,VR) = IPOLY
      PVL(SUCC,VR) = VS
      PVL(EDGV,VR) = 0
      PVL(LOC,VS) = NVC + 3
      PVL(POLG,VS) = IPOLY
      PVL(SUCC,VS) = SUCCIR
      PVL(EDGV,VS) = PVL(EDGV,IRGT)
      PVL(SUCC,ILFT) = VP
      LV = PVL(LOC,ILFT)
      IANG(VP) = ANGLE(VCL(1,LV),VCL(2,LV),XLFT,YLFT,XT,YT)
      LV = PVL(LOC,SUCCIR)
      IANG(VS) = ANGLE(XT,YT,XRGT,YRGT,VCL(1,LV),VCL(2,LV))
      W = 0
      CALL RESVRT(VR,ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,MAXIW,MAXWK,
     $   VCL,PVL,IANG,W,WW,IWK,WK)
C
C     Remove temporary modification to PVL. There are three cases
C     depending on where the endpoint of separator is located:
C     successor of closest vertex to the right of top vertex,
C     predecessor of closest vertex to the left of top vertex,
C     or neither of these.
C
      IF (PVL(SUCC,VS) .EQ. W) THEN
       PVL(SUCC,ILFT) = SUCCIL
       PVL(SUCC,IRGT) = W
       PVL(EDGV,IRGT) = PVL(EDGV,VS)
         IF (PVL(EDGV,IRGT) .GT. 0) PVL(EDGV,PVL(EDGV,IRGT)) = IRGT
      ELSE IF (PVL(SUCC,ILFT) .EQ. W) THEN
       PVL(SUCC,W) = SUCCIL
      ELSE
       PVL(SUCC,ILFT) = SUCCIL
      ENDIF
      IF (IERR .NE. 0) RETURN
C
C     Update PVL with cut edge, i.e. join linked lists of vertices
C     of the hole polygon and the outer boundary polygon into one
C     linked list of vertices by adding the cut edge from the top
C     vertex of hole to the vertex on the outer boundary.
C
      NVERT = NVERT + 2
      VV = NVERT - 1
      WW = NVERT
      LV = PVL(LOC,ITOPHV)
      LW = PVL(LOC,W)
      PVL(LOC,VV) = LV
      PVL(LOC,WW) = LW
      PVL(POLG,VV) = IPOLY
      PVL(POLG,WW) = IPOLY
      PVL(SUCC,VV) = PVL(SUCC,ITOPHV)
      PVL(SUCC,WW) = PVL(SUCC,W)
      PVL(SUCC,ITOPHV) = WW
      PVL(SUCC,W) = VV
      PVL(EDGV,VV) = PVL(EDGV,ITOPHV)
      PVL(EDGV,WW) = PVL(EDGV,W)
      PVL(EDGV,ITOPHV) = W
      PVL(EDGV,W) = ITOPHV
      IF (PVL(EDGV,VV) .GT. 0) PVL(EDGV,PVL(EDGV,VV)) = VV
      IF (PVL(EDGV,WW) .GT. 0) PVL(EDGV,PVL(EDGV,WW)) = WW
      L = PVL(LOC,PVL(SUCC,VV))
      IANG(VV) = ANGLE(VCL(1,LW),VCL(2,LW),VCL(1,LV),VCL(2,LV),
     $   VCL(1,L),VCL(2,L))
      IANG(ITOPHV) = IANG(ITOPHV) - IANG(VV)
      L = PVL(LOC,PVL(SUCC,WW))
      IANG(WW) = ANGLE(VCL(1,LV),VCL(2,LV),VCL(1,LW),VCL(2,LW),
     $   VCL(1,L),VCL(2,L))
      IANG(W) = IANG(W) - IANG(WW)
C
      END
C
C     The following code was excerpted from: minang.f
C
      DOUBLE PRECISION FUNCTION MINANG(XR,YR,XS,YS,IND,ALPHA,THETA,
     $   VCL,PVL,IANG)
      IMPLICIT LOGICAL (A-Z)
      INTEGER IND,PVL(4,*)
      DOUBLE PRECISION ALPHA,IANG(*),THETA,VCL(2,*),XR,XS,YR,YS
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 minimum of the 4 angles at the boundary
C        resulting from using edge joining vertices (XR,YR) and
C        (XS,YS) as a separator.
C
C     Input parameters:
C        XR,YR - coordinates of reflex vertex
C        XS,YS - coordinates of other endpoint of possible separator
C        IND - if positive then (XS,YS) has index IND in PVL; else
C              (XS,YS) is on edge joining vertices with indices -IND
C              and SUCC(-IND) in PVL
C        ALPHA - polar angle of (XS,YS) w.r.t. (XR,YR)
C        THETA - interior angle at reflex vertex
C        VCL(1:2,1:*) - vertex coordinate list
C        PVL(1:4,1:*),IANG(1:*) - polygon vertex list, interior angles
C
C     Returned function value:
C        MINANG - minimum of the 4 angles in radians
C
C     Routines called:
C        ANGLE
C
      DOUBLE PRECISION PI,TOL
      COMMON /GCONST/ PI,TOL
      SAVE /GCONST/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER J,L
      DOUBLE PRECISION ANG,ANGLE,BETA1
C
      IF (IND .GT. 0) THEN
       J = PVL(SUCC,IND)
       ANG = IANG(IND)
      ELSE
       J = PVL(SUCC,-IND)
       ANG = PI
      ENDIF
      L = PVL(LOC,J)
      BETA1 = ANGLE(XR,YR,XS,YS,VCL(1,L),VCL(2,L))
      MINANG = MIN(ALPHA, THETA - ALPHA, ANG - BETA1, BETA1)
      END
C
C     The following code was excerpted from: resvrt.f
C
      SUBROUTINE RESVRT(VR,ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,MAXIW,
     $   MAXWK,VCL,PVL,IANG,W1,W2,IWK,WK)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXIW,MAXPV,MAXVC,MAXWK,NVC,NVERT,VR,W1,W2
      INTEGER IWK(MAXIW),PVL(4,MAXPV)
      DOUBLE PRECISION ANGSPC,ANGTOL,IANG(MAXPV),VCL(2,MAXVC),WK(MAXWK)
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: Resolve reflex vertex of simply connected polygon with
C        one or two separators. The reflex vertex must be a 'simple'
C        vertex of the polygon.
C
C     Input parameters:
C        VR - index in PVL of reflex vertex
C        ANGSPC - angle spacing parameter used in controlling the
C              vertices to be considered as an endpoint of a separator
C        ANGTOL - angle tolerance parameter used in accepting
C              separator(s)
C        NVC - number of positions used in VCL array
C        NVERT - number of positions used in PVL array
C        MAXVC - maximum size available for VCL array
C        MAXPV - maximum size available for PVL array
C        MAXIW - maximum size available for IWK array; should be about
C              3 times number of vertices in polygon
C        MAXWK - maximum size available for WK array; should be about
C              5 times number of vertices in polygon
C        VCL(1:2,1:NVC) - vertex coordinate list
C        PVL(1:4,1:NVERT),IANG(1:NVERT) - polygon vertex list and
C              interior angles
C
C     Updated parameters:
C        NVC,NVERT,VCL,PVL,IANG - are updated if INSVR2 is called
C
C     Output parameters:
C        W1 - index in PVL of vertex which is the endpoint of separator
C              in inner cone or right cone w.r.t. reflex vertex
C        W2 - 0 if there is only one separator; else index in PVL of
C              vertex which is endpoint of 2nd separator in left cone
C
C     Working parameters:
C        IWK(1:MAXIW) - integer work array
C        WK(1:MAXWK) - double precision work array
C
C     Abnormal return:
C        IERR is set to 3, 5, 6, 7, 206, 207, 208, 209, 210, or 212
C
C     Routines called:
C        FNDSEP, INSVR2, VISPOL, VISVRT, VORNBR
C
      INTEGER IERR
      COMMON /GERROR/ IERR
      SAVE /GERROR/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER I,I1,I2,IVIS,IVOR,IVRT,L,MAXN,NVIS,NVOR,NVRT,NVSVRT
      INTEGER THETA,V,WKANG,XC,XVOR,YC,YVOR
      DOUBLE PRECISION ANGSEP,XR,YR
C
C     Determine number of vertices in polygon containing reflex vertex.
C
      NVRT = 0
      V = VR
   10 CONTINUE
         V = PVL(SUCC,V)
         IF (V .NE. VR) THEN
          NVRT = NVRT + 1
          GO TO 10
         ENDIF
      MAXN = NVRT + INT(IANG(VR)/ANGSPC)
      L = PVL(LOC,VR)
      XR = VCL(1,L)
      YR = VCL(2,L)
C
C     Set up work arrays for routine VISPOL, and determine whether there
C     is enough workspace. XC, YC are d.p. arrays of length NVRT in WK,
C     used for the coordinates of the polygon containing the reflex
C     vertex. MAXN positions are reserved for XC, YC since this is the
C     maximum space required by routine VISVRT. IVIS is an integer array
C     of length MAXN in IWK. IVRT is an integer array of length NVRT in
C     IWK used temporarily for storing indices of vertices in PVL.
C
      IF (MAXN + NVRT .GT. MAXIW) THEN
       IERR = 6
       RETURN
      ELSE IF (MAXN + MAXN .GT. MAXWK) THEN
       IERR = 7
       RETURN
      ENDIF
      IVIS = 1
      IVRT = IVIS + MAXN
      XC = 1
      YC = XC + MAXN
      V = PVL(SUCC,VR)
      DO 20 I = 0,NVRT-1
       L = PVL(LOC,V)
       WK(XC+I) = VCL(1,L)
       WK(YC+I) = VCL(2,L)
       IWK(IVRT+I) = V
       V = PVL(SUCC,V)
   20 CONTINUE
      CALL VISPOL(XR,YR,NVRT-1,WK(XC),WK(YC),NVIS,IWK(IVIS))
      IF (IERR .NE. 0) RETURN
C
C     XC, YC now contain visibility polygon coordinates. Update MAXN
C     and set up d.p. array THETA of length MAXN in WK for routine
C     VISVRT. Elements of IVIS are changed to indices of PVL after call.
C
      MAXN = MAXN - NVRT + NVIS + 1
      THETA = YC + MAXN
      IF (THETA + MAXN - 1 .GT. MAXWK) THEN
       IERR = 7
       RETURN
      ENDIF
      CALL VISVRT(ANGSPC,XR,YR,NVIS,WK(XC),WK(YC),IWK(IVIS),MAXN-1,
     $   NVSVRT,WK(THETA))
      WK(THETA+NVSVRT) = IANG(VR)
      DO 30 I = IVIS,IVIS+NVSVRT
       L = IWK(I)
       IF (L .GE. 0) THEN
          IWK(I) = IWK(IVRT+L)
       ELSE
          IWK(I) = -IWK(IVRT-L-1)
       ENDIF
   30 CONTINUE
C
C     XC, YC now contain coord. of visible vertices to be considered
C     as an endpoint of a separator. Set up work arrays for routine
C     VORNBR. Integer array IVOR and d.p. arrays XVOR, YVOR, each of
C     length NVSVRT+1, are added at the end of IWK and WK arrays.
C
      IVOR = IVIS + NVSVRT + 1
      XVOR = THETA + NVSVRT + 1
      YVOR = XVOR + NVSVRT + 1
      IF (IVOR + NVSVRT .GT. MAXIW) THEN
       IERR = 6
       RETURN
      ELSE IF (YVOR + NVSVRT .GT. MAXWK) THEN
       IERR = 7
       RETURN
      ENDIF
      CALL VORNBR(XR,YR,NVSVRT,WK(XC),WK(YC),NVOR,IWK(IVOR),WK(XVOR),
     $   WK(YVOR))
      IF (IERR .NE. 0) RETURN
C
C     Set up d.p. array WKANG of length NVOR+1 <= NVSVRT+1 in WK for
C     routine FNDSEP. Only Voronoi neighbours are considered as an
C     endpoint of a separator in the first call to FNDSEP. If the
C     minimum angle created at the boundary by the separator(s) is too
C     small, then a second call is made to FNDSEP in which all visible
C     vertices are considered as an endpoint of a separator.
C
      WKANG = XVOR
      IF (IWK(IVOR+NVOR) .EQ. NVSVRT) NVOR = NVOR - 1
      IF (IWK(IVOR) .EQ. 0) THEN
       IVOR = IVOR + 1
       NVOR = NVOR - 1
      ENDIF
      CALL FNDSEP(ANGTOL+ANGTOL,XR,YR,NVSVRT,WK(XC),WK(YC),IWK(IVIS),
     $   WK(THETA),NVOR,IWK(IVOR),VCL,PVL,IANG,ANGSEP,I1,I2,WK(WKANG))
      IF (ANGSEP .LT. ANGTOL) THEN
       IVRT = IVIS + NVSVRT + 1
       DO 40 I = 1,NVSVRT-1
          IWK(IVRT+I-1) = I
   40    CONTINUE
         CALL FNDSEP(ANGTOL+ANGTOL,XR,YR,NVSVRT,WK(XC),WK(YC),IWK(IVIS),
     $      WK(THETA),NVSVRT-2,IWK(IVRT),VCL,PVL,IANG,ANGSEP,I1,I2,
     $      WK(WKANG))
      ENDIF
C
C     Insert endpoint(s) of separator(s) in vertex coordinate list and
C     polygon vertex list data structures, if they are not yet there.
C
      IF (I2 .EQ. -1) THEN
       W2 = 0
      ELSE IF (IWK(IVIS+I2) .LT. 0) THEN
       CALL INSVR2(WK(XC+I2),WK(YC+I2),-IWK(IVIS+I2),NVC,NVERT,MAXVC,
     $      MAXPV,VCL,PVL,IANG,W2)
       IF (IERR .NE. 0) RETURN
      ELSE
       W2 = IWK(IVIS+I2)
      ENDIF
      IF (IWK(IVIS+I1) .LT. 0) THEN
       CALL INSVR2(WK(XC+I1),WK(YC+I1),-IWK(IVIS+I1),NVC,NVERT,MAXVC,
     $      MAXPV,VCL,PVL,IANG,W1)
       IF (IERR .NE. 0) RETURN
      ELSE
       W1 = IWK(IVIS+I1)
      ENDIF
      END
C
C     The following code was excerpted from: spdec2.f
C
      SUBROUTINE SPDEC2(ANGSPC,ANGTOL,NVC,NPOLG,NVERT,NHOLE,NHOLA,MAXVC,
     $   MAXHV,MAXPV,MAXIW,MAXWK,HOLV,VCL,REGNUM,HVL,PVL,IANG,IWK,WK)
      IMPLICIT LOGICAL (A-Z)
      INTEGER MAXHV,MAXIW,MAXPV,MAXVC,MAXWK,NHOLA,NHOLE,NPOLG,NVC,NVERT
      INTEGER HOLV(*),HVL(MAXHV),IWK(MAXIW),PVL(4,MAXPV),REGNUM(MAXHV)
      DOUBLE PRECISION ANGSPC,ANGTOL,IANG(MAXPV),VCL(2,MAXVC),WK(MAXWK)
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: Decompose general polygonal region with interfaces and
C        holes into simple polygons using vertex coordinate list,
C        head vertex list, and polygon vertex list data structures.
C
C     Input parameters:
C        ANGSPC - angle spacing parameter in radians used in controlling
C              vertices to be considered as an endpoint of a separator
C        ANGTOL - angle tolerance parameter in radians used in accepting
C              separator(s)
C        NVC - number of vertex coordinates or positions used in VCL
C              array
C        NPOLG - number of polygonal subregions or positions used in
C              HVL array
C        NVERT - number of polygon vertices or positions used in PVL
C              array
C        NHOLE - number of holes and hole interfaces
C        NHOLA - number of 'attached' holes; these holes are attached
C              to the outer boundary of a subregion through vertices
C              or cut interfaces and have their edges in consecutive
C              order on the boundary
C        MAXVC - maximum size available for VCL array, should be >=
C              number of vertex coordinates required for decomposition
C        MAXHV - maximum size available for HVL, REGNUM arrays, should
C              be >= number of polygons required for decomposition
C        MAXPV - maximum size available for PVL, IANG arrays; should be
C              >= number of polygon vertices required for decomposition
C        MAXIW - maximum size available for IWK array; should be about
C              3 times maximum number of vertices in any polygon
C        MAXWK - maximum size available for WK array; should be about
C              5 times maximum number of vertices in any polygon
C        HOLV(1:NHOLE*2+NHOLA) - indices in PVL of bottom or top vertex
C              of holes; first (next) NHOLE entries are for top (bottom)
C              vertices of holes and hole interfaces, with top (bottom)
C              vertices sorted in decreasing (increasing) lexicographic
C              (y,x) order of coord; last NHOLA entries are for attached
C              holes; if bottom vertex of attached hole is a simple
C              vertex of boundary curve containing the hole then entry
C              contains index of bottom vertex otherwise entry contains
C              index of top vertex (which is simple)
C        VCL(1:2,1:NVC) - vertex coordinate list
C        REGNUM(1:NPOLG) - region numbers
C        HVL(1:NPOLG) - head vertex list
C        PVL(1:4,1:NVERT),IANG(1:NVERT) - polygon vertex list and
C              interior angles; see routine DSPGDC for more details
C        [Note: The data structures should be as output from routines
C              DSMCPR or DSPGDC.]
C
C     Updated parameters:
C        NVC,NPOLG,NVERT,VCL,REGNUM,HVL,PVL,IANG
C
C     Working parameters:
C        IWK(1:MAXIW) - integer work array
C        WK(1:MAXWK) - double precision work array
C
C     Abnormal return:
C        IERR is set to 3, 4, 5, 6, 7, 206 to 210, 212, 218, or 219
C
C     Routines called:
C        INSED2, JNHOLE, RESVRT
C
      INTEGER IERR
      DOUBLE PRECISION PI,TOL
      COMMON /GERROR/ IERR
      COMMON /GCONST/ PI,TOL
      SAVE /GERROR/,/GCONST/
C
      INTEGER EDGV,LOC,POLG,SUCC
      PARAMETER (LOC = 1, POLG = 2, SUCC = 3, EDGV = 4)
C
      INTEGER I,J,P,VR,W1,W2
      DOUBLE PRECISION PIPTOL
      LOGICAL CI,CJ
C
C     For each simple hole, find cut edge from top vertex of hole to
C     a point on the outer boundary above top vertex, and update
C     VCL, HVL, PVL, IANG.
C
      PIPTOL = PI + TOL
      DO 10 I = 1,NHOLE
       CALL JNHOLE(HOLV(I),ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,MAXIW,
     $      MAXWK,VCL,HVL,PVL,IANG,IWK,WK)
       IF (IERR .NE. 0) RETURN
   10 CONTINUE
C
C     Resolve remaining vertices in HOLV array if they are reflex
C     vertices. These vertices may no longer be reflex if they are the
C     endpoint of a cut edge from the top vertex of another hole or
C     of a previous separator.
C
      DO 20 I = NHOLE+1,NHOLE+NHOLE+NHOLA
       VR = HOLV(I)
       IF (IANG(VR) .GT. PIPTOL) THEN
            CALL RESVRT(VR,ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,MAXIW,
     $         MAXWK,VCL,PVL,IANG,W1,W2,IWK,WK)
          IF (IERR .NE. 0) RETURN
          CALL INSED2(VR,W1,NPOLG,NVERT,MAXHV,MAXPV,VCL,REGNUM,HVL,
     $         PVL,IANG)
          IF (IERR .NE. 0) RETURN
            IF (W2 .GT. 0) CALL INSED2(VR,W2,NPOLG,NVERT,MAXHV,MAXPV,
     $         VCL,REGNUM,HVL,PVL,IANG)
          IF (IERR .NE. 0) RETURN
         ENDIF
   20 CONTINUE
      IF (NHOLA .EQ. 0) RETURN
C
C     Check that polygons are simple. If polygon is simply-connected and
C     not simple then find a simple reflex vertex in polygon to resolve.
C
      P = 1
   30 CONTINUE
      IF (P .GT. NPOLG) RETURN
       I = HVL(P)
   40    CONTINUE
          IF (PVL(POLG,PVL(EDGV,I)) .EQ. P) GO TO 50
          I = PVL(SUCC,I)
       IF (I .NE. HVL(P)) GO TO 40
       P = P + 1
       GO TO 30
   50    CONTINUE
       CI = .TRUE.
   60    CONTINUE
          J = PVL(SUCC,I)
          CJ = (PVL(POLG,PVL(EDGV,J)) .EQ. P)
       IF (CI .OR. CJ .OR. IANG(J) .LE. PIPTOL) THEN
          I = J
          CI = CJ
          GO TO 60
       ENDIF
       VR = J
         CALL RESVRT(VR,ANGSPC,ANGTOL,NVC,NVERT,MAXVC,MAXPV,MAXIW,
     $      MAXWK,VCL,PVL,IANG,W1,W2,IWK,WK)
       IF (IERR .NE. 0) RETURN
       CALL INSED2(VR,W1,NPOLG,NVERT,MAXHV,MAXPV,VCL,REGNUM,HVL,
     $      PVL,IANG)
       IF (IERR .NE. 0) RETURN
         IF (W2 .GT. 0) CALL INSED2(VR,W2,NPOLG,NVERT,MAXHV,MAXPV,
     $      VCL,REGNUM,HVL,PVL,IANG)
       IF (IERR .NE. 0) RETURN
      GO TO 30
      END

Generated by  Doxygen 1.6.0   Back to index