COMMON NCOUNT DIMENSION XX(2,25),IN(25),IH(25) DIMENSION X(25),Y(25) INTEGER IWORK(50) INTEGER IL(50) DATA IWORK/50*0/ NCOUNT=0 C NCOUNT IS TOTAL NUMBER OF POINTS PASSED TO SPLIT READ(5,1)N 1 FORMAT(I5) WRITE(6,1)N N1=N+1 DO 2 I=1,N J=N1-I 2 IN(J)=I C ARRAY IN CONTAINS INDICES 1-N IN REVERSE ORDER DO 3 I=1,N 3 READ(5,4)XX(1,I),XX(2,I) 4 FORMAT(2F10.5) DO 5 I=1,N J=IN(I) 5 WRITE(6,4)XX(1,J),XX(2,J) DO 10 M=4,N CALL CONVEX(N,XX,M,IN,IWORK,IWORK(N+1),IH,NHULL,IL) IK=IL(1) DO 6 I=1,NHULL J=IH(IK) X(I)=XX(1,J) Y(I)=XX(2,J) 6 IK=IL(IK) WRITE(6,7)M,NHULL,NCOUNT 7 FORMAT(12H0SAMPLE SIZE ,I5,9H VERTICES ,I5,6H SPLIT ,I5) DO 8 I=1,NHULL 8 WRITE(6,9)X(I),Y(I) 9 FORMAT(1X,2F10.5) 10 CONTINUE STOP END SUBROUTINE SPLIT(N,X,M,IN,II,JJ,S,IABV,NA,MAXA,IBEL, 1 NB,MAXB) C THIS SUBROUTINE TAKES THE M POINTS OF ARRAY X WHOSE C SUBSCRIPTS ARE IN ARRAY IN AND PARTITIONS THEM BY THE C LINE JOINING THE TWO POINTS IN ARRAY X WHOSE SUBSCRIPTS C ARE II AND JJ. THE SUBSCRIPTS OF THE POINTS ABOVE THE C LINE ARE PUT INTO ARRAY IABV, AND THE SUBSCRIPTS OF THE C POINTS BELOW ARE PUT INTO ARRAY IBEL. NA AND NB ARE, C RESPECTIVELY, THE NUMBER OF POINTS ABOVE THE LINE AND THE C NUMBER BELOW. MAXA AND MAXB ARE THE SUBSCRIPTS FOR ARRAY C X OF THE POINT FURTHEST ABOVE THE LINE AND THE POINT C FURTHEST BELOW, RESPECTIVELY. IF EITHER SUBSET IS NULL C THE CORRESPONDING SUBSCRIPT (MAXA OR MAXB) IS SET TO ZERO C FORMAL PARAMETERS C INPUT C N INTEGER TOTAL NUMBER OF DATA POINTS C X REAL ARRAY (2,N) (X,Y) CO-ORDINATES OF THE DATA C M INTEGER NUMBER OF POINTS IN INPUT SUBSET C IN INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE C POINTS IN THE INPUT SUBSET C II INTEGER SUBSCRIPT FOR ARRAY X OF ONE POINT C ON THE PARTITIONING LINE C JJ INTEGER SUBSCRIPT FOR ARRAY X OF ANOTHER C POINT ON THE PARTITIONING LINE C S INTEGER SWITCH TO DETERMINE OUTPUT. REFER C TO COMMENTS BELOW C OUTPUT C IABV INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE C POINTS ABOVE THE PARTITIONING LINE C NA INTEGER NUMBER OF ELEMENTS IN IABV C MAXA INTEGER SUBSCRIPT FOR ARRAY X OF POINT C FURTHEST ABOVE THE LINE. SET TO C ZERO IF NA IS ZERO C IBEL INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE C POINTS BELOW THE PARTITIONING LINE C NB INTEGER NUMBER OF ELEMENTS IN IBEL C MAXB INTEGER SUBSCRIPT FOR ARRAY X OF POINT C FURTHEST BELOW THE LINE. SET TO C ZERO IF NB IS ZERO DIMENSION X(2,N) DIMENSION IN(M),IABV(M),IBEL(M) INTEGER S C IF S = 2 DONT SAVE IBEL,NB,MAXB. C IF S =-2 DONT SAVE IABV,NA,MAXA. C OTHERWISE SAVE EVERYTHING C IF S IS POSITIVE THE ARRAY BEING PARTITIONED IS ABOVE C THE INITIAL PARTITIONING LINE. IF IT IS NEGATIVE, THEN C THE SET OF POINTS IS BELOW. LOGICAL T T=.FALSE. C CHECK TO SEE IF THE LINE IS VERTICAL IF(X(1,JJ).NE.X(1,II))GOTO 1 XT=X(1,II) DIR=SIGN(1.,X(2,JJ)-X(2,II))*SIGN(1.,FLOAT(S)) T=.TRUE. GOTO 2 1 A=(X(2,JJ)-X(2,II))/(X(1,JJ)-X(1,II)) B=X(2,II)-A*X(1,II) 2 UP=0. NA=0 MAXA=0 DOWN=0. NB=0 MAXB=0 DO 6 I=1,M IS=IN(I) IF(T)GOTO 3 Z=X(2,IS)-A*X(1,IS)-B GOTO 4 3 Z=DIR*(X(1,IS)-XT) 4 IF(Z.LE.0.)GOTO 5 C THE POINT IS ABOVE THE LINE IF(S.EQ.-2)GOTO 6 NA=NA+1 IABV(NA)=IS IF(Z.LT.UP)GOTO 6 UP=Z MAXA=NA GOTO 6 5 IF(S.EQ.2)GOTO 6 IF(Z.GE.0.)GOTO 6 C THE POINT IS BELOW THE LINE NB=NB+1 IBEL(NB)=IS IF(Z.GT.DOWN)GOTO 6 DOWN=Z MAXB=NB 6 CONTINUE RETURN END SUBROUTINE CONVEX(N,X,M,IN,IA,IB,IH,NH,IL) C THIS SUBROUTINE DETERMINES WHICH OF THE M POINTS OF ARRAY C X WHOSE SUBSCRIPTS ARE IN ARRAY IN ARE VERTICES OF THE C MINIMUM AREA CONVEX POLYGON CONTAINING THE M POINTS. THE C SUBSCRIPTS OF THE VERTICES ARE PLACED IN ARRAY IH IN THE C ORDER THEY ARE FOUND. NH IS THE NUMBER OF ELEMENTS IN C ARRAY IH AND ARRAY IL. ARRAY IL IS A LINKED LIST GIVING C THE ORDER OF THE ELEMENTS OF ARRAY IH IN A COUNTER C CLOCKWISE DIRECTION. THIS ALGORITHM CORRESPONDS TO A C PREORDER TRAVERSAL OF A CERTAIN BINARY TREE. EACH VERTEX C OF THE BINARY TREE REPRESENTS A SUBSET OF THE M POINTS. C AT EACH STEP THE SUBSET OF POINTS CORRESPONDING TO THE C CURRENT VERTEX OF THE TREE IS PARTITIONED BY A LINE C JOINING TWO VERTICES OF THE CONVEX POLYGON. THE LEFT SON C VERTEX IN THE BINARY TREE REPRESENTS THE SUBSET OF POINTS C ABOVE THE PARTITIONING LINE AND THE RIGHT SON VERTEX, THE C SUBSET BELOW THE LINE. THE LEAVES OF THE TREE REPRESENT C EITHER NULL SUBSETS OR SUBSETS INSIDE A TRIANGLE WHOSE C VERTICES COINCIDE WITH VERTICES OF THE CONVEX POLYGON. C FORMAL PARAMETERS C INPUT C N INTEGER TOTAL NUMBER OF DATA POINTS C X REAL ARRAY (2,N) (X,Y) CO-ORDINATES OF THE DATA C M INTEGER NUMBER OF POINTS IN THE INPUT SUBSET C IN INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE POINTS C IN THE INPUT SUBSET C WORK AREA C IA INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF LEFT SON C SUBSETS. SEE COMMENTS AFTER DIMENSION C STATEMENTS C IB INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF RIGHT SON C SUBSETS C OUTPUT C IH INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE C VERTICES OF THE CONVEX HULL C NH INTEGER NUMBER OF ELEMENTS IN ARRAY IH AND C ARRAY IL. SAME AS NUMBER OF VERTICES C OF THE CONVEX POLYGON C IL INTEGER ARRAY (M) A LINKED LIST GIVING IN ORDER IN A C COUNTER-CLOCKWISE DIRECTION THE C ELEMENTS OF ARRAY IH DIMENSION X(2,N) DIMENSION IN(M),IA(M),IB(M),IH(M),IL(M) C THE UPPER END OF ARRAY IA IS USED TO STORE TEMPORARILY C THE SIZES OF THE SUBSETS WHICH CORRESPOND TO RIGHT SON C VERTICES, WHILE TRAVERSING DOWN THE LEFT SONS WHEN ON THE C LEFT HALF OF THE TREE, AND TO STORE THE SIZES OF THE LEFT C SONS WHILE TRAVERSING THE RIGHT SONS(DOWN THE RIGHT HALF) LOGICAL MAXE,MINE IF(M.EQ.1)GOTO 22 IL(1)=2 IL(2)=1 KN=IN(1) KX=IN(2) IF(M.EQ.2)GOTO 21 MP1=M+1 MIN=1 MX=1 KX=IN(1) MAXE=.FALSE. MINE=.FALSE. C FIND TWO VERTICES OF THE CONVEX HULL FOR THE INITIAL C PARTITION DO 6 I=2,M J=IN(I) IF(X(1,J)-X(1,KX))3,1,2 1 MAXE=.TRUE. GOTO 3 2 MAXE=.FALSE. MX=I KX=J 3 IF(X(1,J)-X(1,KN))5,4,6 4 MINE=.TRUE. GOTO 6 5 MINE=.FALSE. MIN=I KN=J 6 CONTINUE C IF THE MAX AND MIN ARE EQUAL, ALL M POINTS LIE ON A C VERTICAL LINE IF(KX.EQ.KN)GOTO 18 C IF MAXE (OR MINE) HAS THE VALUE TRUE THERE ARE SEVERAL C MAXIMA (OR MINIMA) WITH EQUAL FIRST COORDINATES IF(MAXE.OR.MINE)GOTO 23 7 IH(1)=KX IH(2)=KN NH=3 INH=1 NIB=1 MA=M IN(MX)=IN(M) IN(M)=KX MM=M-2 IF(MIN.EQ.M)MIN=MX IN(MIN)=IN(M-1) IN(M-1)=KN C BEGIN BY PARTITIONING THE ROOT OF THE TREE CALL SPLIT(N,X,MM,IN,IH(1),IH(2),0,IA,MB,MXA,IB,IA(MA), 1 MXBB) C FIRST TRAVERSE THE LEFT HALF OF THE TREE C START WITH THE LEFT SON 8 NIB=NIB+IA(MA) MA=MA-1 9 IF(MXA.EQ.0)GOTO 11 IL(NH)=IL(INH) IL(INH)=NH IH(NH)=IA(MXA) IA(MXA)=IA(MB) MB=MB-1 NH=NH+1 IF(MB.EQ.0)GOTO 10 ILINH=IL(INH) CALL SPLIT(N,X,MB,IA,IH(INH),IH(ILINH),1,IA,MBB,MXA, 1 IB(NIB),IA(MA),MXB) MB=MBB GOTO 8 C THEN THE RIGHT SON 10 INH=IL(INH) 11 INH=IL(INH) MA=MA+1 NIB=NIB-IA(MA) IF(MA.GE.M)GOTO 12 IF(IA(MA).EQ.0)GOTO 11 ILINH=IL(INH) C ON THE LEFT SIDE OF THE TREE, THE RIGHT SON OF A RIGHT SON C MUST REPRESENT A SUBSET OF POINTS WHICH IS INSIDE A C TRIANGLE WITH VERTICES WHICH ARE ALSO VERTICES OF THE C CONVEX POLYGON AND HENCE THE SUBSET MAY BE NEGLECTED. CALL SPLIT(N,X,IA(MA),IB(NIB),IH(INH),IH(ILINH),2,IA, 1 MB,MXA,IB(NIB),MBB,MXB) IA(MA)=MBB GOTO 9 C NOW TRAVERSE THE RIGHT HALF OF THE TREE 12 MXB=MXBB MA=M MB=IA(MA) NIA=1 IA(MA)=0 C START WITH THE RIGHT SON 13 NIA=NIA+IA(MA) MA=MA-1 14 IF(MXB.EQ.0)GOTO 16 IL(NH)=IL(INH) IL(INH)=NH IH(NH)=IB(MXB) IB(MXB)=IB(MB) MB=MB-1 NH=NH+1 IF(MB.EQ.0)GOTO 15 ILINH=IL(INH) CALL SPLIT(N,X,MB,IB(NIB),IH(INH),IH(ILINH),-1,IA(NIA), 1 IA(MA),MXA,IB(NIB),MBB,MXB) MB=MBB GOTO 13 C THEN THE LEFT SON 15 INH=IL(INH) 16 INH=IL(INH) MA=MA+1 NIA=NIA-IA(MA) IF(MA.EQ.MP1)GOTO 17 IF(IA(MA).EQ.0)GOTO 16 ILINH=IL(INH) C ON THE RIGHT SIDE OF THE TREE, THE LEFT SON OF A LEFT SON C MUST REPRESENT A SUBSET OF POINTS WHICH IS INSIDE A C TRIANGLE WITH VERTICES WHICH ARE ALSO VERTICES OF THE C CONVEX POLYGON AND HENCE THE SUBSET MAY BE NEGLECTED. CALL SPLIT(N,X,IA(MA),IA(NIA),IH(INH),IH(ILINH),-2, 1 IA(NIA),MBB,MXA,IB(NIB),MB,MXB) GOTO 14 17 NH=NH-1 RETURN C ALL THE SPECIAL CASES ARE HANDLED DOWN HERE C IF ALL THE POINTS LIE ON A VERTICAL LINE 18 KX=IN(1) KN=IN(1) DO 20 I=1,M J=IN(I) IF(X(2,J).LE.X(2,KX))GOTO 19 MX=I KX=J 19 IF(X(2,J).GE.X(2,KN))GOTO 20 MIN=I KN=J 20 CONTINUE IF(KX.EQ.KN)GOTO 22 C IF THERE ARE ONLY TWO POINTS 21 IH(1)=KX IH(2)=KN NH=3 IF((X(1,KN).EQ.X(1,KX)).AND.(X(2,KN).EQ.X(2,KX)))NH=2 GOTO 17 C IF THERE IS ONLY ONE POINT 22 NH=2 IH(1)=IN(1) IL(1)=1 GOTO 17 C MULTIPLE EXTREMES ARE HANDLED HERE C IF THERE ARE SEVERAL POINTS WITH THE (SAME) LARGEST C FIRST COORDINATE 23 IF(.NOT.MAXE)GOTO 25 DO 24 I=1,M J=IN(I) IF(X(1,J).NE.X(1,KX))GOTO 24 IF(X(2,J).LE.X(2,KX))GOTO 24 MX=I KX=J 24 CONTINUE C IF THERE ARE SEVERAL POINTS WITH THE (SAME) SMALLEST C FIRST COORDINATE 25 IF(.NOT.MINE)GOTO 7 DO 26 I=1,M J=IN(I) IF(X(1,J).NE.X(1,KN))GOTO 26 IF(X(2,J).GE.X(2,KN))GOTO 26 MIN=I KN=J 26 CONTINUE GOTO 7 END 20 2.0 0.0 1.73 -1.0 1.0 1.73 0.0 2.0 0.1 0.1 -1.0 -1.73 0.2 -0.2 -1.73 1.0 -0.3 0.3 0.0 -2.0 -0.4 -0.4 -2.0 0.0 0.5 0.5 1.73 1.0 0.6 -0.6 -1.0 1.73 -0.7 0.7 -1.73 -1.0 -0.8 -0.8 1.0 -1.73