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