C
C Program TRGL to divide polygons on a curved surface into triangles, C right-handed with respect to the surface normals C C Version: 5.40 C Date: 2000, January 21 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: klimes@seis.karlov.mff.cuni.cz C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Data specifying input files: C VRTX='string'... Name of the file with vertices of the polygons. C Description of file VRTX C Default: VRTX='vrtx.out' C PLGN='string'... Name of the file describing the polygons. C Description of file PLGN C Default: PLGN='plgn.out' C Data specifying output file: C TRGL='string'... Name of the file describing the triangles. C Description of file TRGL C Default: TRGL='trgl.out' C C C Input file VRTX with the vertices: C (1) None to several strings terminated by / (a slash) C (2) For each vertex data (2.1): C (2.1) 'NAME',X1,X2,X3,Z1,Z2,Z3,/ C 'NAME'... Name of the vertex. Not considered. May be blank. C X1,X2,X3... Coordinates of the vertex. C Z1,Z2,Z3... Normal to the surface at the vertex. C /... None to several values terminated by a slash. C (3) / or end of file. C C C Input file PLGN with the polygons: C (1) For each polygon data (1.1): C (1.1) I1,I2,...,IN,/ C I1,I2,...,IN... Indices of N vertices of the polygon. C The vertices in file VRTX are indexed by positive integers C according to their order. C /... List of vertices must be terminated by a slash. C (2) / or end of file. C C C Output file TRGL with the triangles: C (1) For each polygon data (1.1): C (1.1) I1,I2,I3,/ C I1,I2,I3... Indices of 3 vertices of the triangle, right-handed C with respect to the given surface normals. C The vertices in file VRTX are indexed by positive integers C according to their order. C /... List of vertices is terminated by a slash. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C....................................................................... C EXTERNAL LENGTH,ERROR,WARN,RSEP1,RSEP3T INTEGER LENGTH C C....................................................................... C C Filenames and parameters: CHARACTER*80 FSEP,FVRTX,FPLGN,FTRGL INTEGER LU,IUNDEF,MVRTX PARAMETER (LU=1,IUNDEF=-999999,MVRTX=1000) C Input data: CHARACTER*9 FORMAT CHARACTER*80 TEXT C Other variables: INTEGER NVRTX,NPLGN,NTRGL,NBIGON,N,I,I1,I2,I3,J1,J2,J3,K1,K2,N1,N2 REAL X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,A1,A2,A3,D,DMAX,DMIN,HAND,S,C,AUX C C MVRTX...Maximum number of vertices of one polygon. C NVRTX...Number of storage locations for the vertices, i.e. 6 times C the number of vertices. C NPLGN...Last storage location for polygons, i.e. NVRTX + 4 times C the number of future triangles. Each polygon with N C vertices occupies 4*(N-2) locations because it will be C split into N-2 triangles. 4*(N-2)-N locations following C N vertex indices are filled wth zeros. C NTRGL...Last storage location with triangles, i.e. NVRTX + 4 times C the number of created triangles. C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+TRGL: Enter input filename: ' FSEP=' ' READ (*,*) FSEP WRITE(*,'(A)') '+TRGL: Working... ' C C Reading all data from the SEP file into the memory: IF (FSEP.NE.' ') THEN CALL RSEP1(LU,FSEP) ELSE C TRGL-11 CALL ERROR('TRGL-11: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C C Reading input and output filenames: CALL RSEP3T('VRTX',FVRTX,'vrtx.out') CALL RSEP3T('PLGN',FPLGN,'plgn.out') CALL RSEP3T('TRGL',FTRGL,'trgl.out') C C Reading vertices: OPEN(LU,FILE=FVRTX) READ(LU,*) (TEXT,I=1,20) NVRTX=0 10 CONTINUE IF(NVRTX+6.GT.MRAM) THEN C TRGL-01 CALL ERROR('TRGL-01: Too small array RAM') END IF TEXT='$' READ(LU,*,END=18) TEXT,(RAM(I),I=NVRTX+1,NVRTX+6) IF(TEXT.EQ.'$') THEN GO TO 18 END IF C Normalizing the normal AUX=SQRT(RAM(NVRTX+4)**2+RAM(NVRTX+5)**2+RAM(NVRTX+6)**2) RAM(NVRTX+4)=RAM(NVRTX+4)/AUX RAM(NVRTX+5)=RAM(NVRTX+5)/AUX RAM(NVRTX+6)=RAM(NVRTX+6)/AUX NVRTX=NVRTX+6 GO TO 10 18 CONTINUE CLOSE(LU) C C Reading polygons: DO 19 I=NVRTX+1,MRAM IRAM(I)=0 19 CONTINUE OPEN(LU,FILE=FPLGN) NPLGN=NVRTX NBIGON=0 20 CONTINUE IRAM(NPLGN+1)=IUNDEF READ(LU,*,END=29) (IRAM(I),I=NPLGN+1,MIN0(NPLGN+MVRTX+1,MRAM)) IF(IRAM(NPLGN+1).EQ.IUNDEF) THEN GO TO 29 END IF DO 21 I=NPLGN+1,MIN0(NPLGN+MVRTX+1,MRAM) IF(IRAM(I).LE.0) THEN C Number of polygon vertices N=I-1-NPLGN GO TO 22 ELSE IF(IRAM(I).GT.NVRTX/6) THEN C TRGL-02 WRITE(TEXT,'(A,I6)') 'TRGL-02: Wrong vertex index:',IRAM(I) CALL ERROR(TEXT(1:LENGTH(TEXT))) END IF 21 CONTINUE C TRGL-03 CALL ERROR('TRGL-03: Too many vertices in polygons') 22 CONTINUE IF(N.LT.3) THEN C TRGL-04 CALL ERROR('TRGL-04: Polygon of less than 3 vertices') END IF IF(NPLGN+4*(N-2).GT.MRAM) THEN C TRGL-05 CALL ERROR('TRGL-05: Too small array RAM') END IF C Checking vertex indices: DO 24 I2=NPLGN+1,NPLGN+N DO 23 I1=I2+1,NPLGN+N IF(IRAM(I2).EQ.IRAM(I1)) THEN C TRGL-06 WRITE(TEXT,'(A,I6)') * 'TRGL-06: The same vertex twice in a polygon:',IRAM(I2) CALL ERROR(TEXT(1:LENGTH(TEXT))) C All vertices of a polygon must have different indices. END IF 23 CONTINUE 24 CONTINUE C Check for identical indices: 25 CONTINUE DO 27 I2=NPLGN+1,NPLGN+N J1=6*(IRAM(I2)-1) IF(I2.LT.NPLGN+N) THEN J2=6*(IRAM(I2+1)-1) ELSE J2=6*(IRAM(NPLGN+1)-1) END IF IF(RAM(J1+1).EQ.RAM(J2+1).AND. * RAM(J1+2).EQ.RAM(J2+2).AND. * RAM(J1+3).EQ.RAM(J2+3)) THEN IF(RAM(J1+1).NE.RAM(J2+1).OR. * RAM(J1+2).NE.RAM(J2+2).OR. * RAM(J1+3).NE.RAM(J2+3)) THEN C TRGL-10 WRITE(TEXT,'(A,3I6)') * 'TRGL-10: Different normals at identical points:', * IRAM(I2),IRAM(I2+1) CALL ERROR(TEXT(1:LENGTH(TEXT))) C Two subsequent points of a polygon have the same C coordinates but different normals. END IF DO 26 I1=I2+1,NPLGN+N-1 IRAM(I1)=IRAM(I1+1) 26 CONTINUE IRAM(NPLGN+N)=0 N=N-1 GO TO 25 END IF 27 CONTINUE IF(N.LT.3) THEN NBIGON=NBIGON+1 GO TO 20 END IF C Leaving the space for N-2 triangles, each terminated by zero: DO 28 I=NPLGN+N+1,NPLGN+4*(N-2) IRAM(I)=0 28 CONTINUE NPLGN=NPLGN+4*(N-2) GO TO 20 29 CONTINUE IF(NBIGON.GT.0) THEN C TRGL-53 WRITE(TEXT,'(A,3I6)') * 'TRGL-53: Number of polygons with less than 3 vertices:', * NBIGON CALL WARN(TEXT(1:LENGTH(TEXT))) C Polygons with less than 3 vertices are created from polygons C with subsequent vertices of identical coordinates. END IF CLOSE(LU) C C Dividing polygons into triangles: IF(NVRTX.EQ.NPLGN) THEN C There are no polygons GO TO 49 END IF NTRGL=NVRTX C Loop over polygons 30 CONTINUE C Determining number N of polygon vertices DO 31 I=NTRGL+1,NPLGN IF(IRAM(I).LE.0) THEN N=I-1-NTRGL GO TO 32 END IF 31 CONTINUE 32 CONTINUE IF(N.EQ.3) THEN C Current polygon is a triangle, proceeding to the next polygon NTRGL=NTRGL+4 IF(NTRGL.LT.NPLGN) THEN GO TO 30 ELSE C All polygons are divided into triangles GO TO 49 END IF END IF C C Dividing the polygon into 2 polygons: K1=0 DMIN=1.0E20 C Loop over polygon vertices C (a) Determination of the handedness of the polygon HAND=0. DO 35 I2=1,N I3=MOD(I2,N)+1 IF(I2.GT.1) THEN I1=I2-1 ELSE I1=N END IF J1=6*(IRAM(NTRGL+I1)-1) J2=6*(IRAM(NTRGL+I2)-1) J3=6*(IRAM(NTRGL+I3)-1) X1=RAM(J2+1)-RAM(J1+1) X2=RAM(J2+2)-RAM(J1+2) X3=RAM(J2+3)-RAM(J1+3) Y1=RAM(J3+1)-RAM(J2+1) Y2=RAM(J3+2)-RAM(J2+2) Y3=RAM(J3+3)-RAM(J2+3) Z1=X2*Y3-X3*Y2 Z2=X3*Y1-X1*Y3 Z3=X1*Y2-X2*Y1 S=Z1*RAM(J2+4)+Z2*RAM(J2+5)+Z3*RAM(J2+6) C=X1*Y1+X2*Y2+X3*Y3 HAND=HAND+ATAN2(S,C) 35 CONTINUE C (b) Selection of the triangle to be separated DO 37 I2=1,N I3=MOD(I2,N)+1 IF(I2.GT.1) THEN I1=I2-1 ELSE I1=N END IF J1=6*(IRAM(NTRGL+I1)-1) J2=6*(IRAM(NTRGL+I2)-1) J3=6*(IRAM(NTRGL+I3)-1) X1=RAM(J2+1)-RAM(J1+1) X2=RAM(J2+2)-RAM(J1+2) X3=RAM(J2+3)-RAM(J1+3) Y1=RAM(J3+1)-RAM(J2+1) Y2=RAM(J3+2)-RAM(J2+2) Y3=RAM(J3+3)-RAM(J2+3) Z1=X2*Y3-X3*Y2 Z2=X3*Y1-X1*Y3 Z3=X1*Y2-X2*Y1 S=Z1*RAM(J2+4)+Z2*RAM(J2+5)+Z3*RAM(J2+6) IF(S*HAND.GT.0.) THEN C Normal Z to the separation plane X1=RAM(J3+1)-RAM(J1+1) X2=RAM(J3+2)-RAM(J1+2) X3=RAM(J3+3)-RAM(J1+3) Y1=RAM(J3+4)+RAM(J1+4) Y2=RAM(J3+5)+RAM(J1+5) Y3=RAM(J3+6)+RAM(J1+6) Z1=X2*Y3-X3*Y2 Z2=X3*Y1-X1*Y3 Z3=X1*Y2-X2*Y1 C Point on the diagonal Y1=RAM(J1+1) Y2=RAM(J1+2) Y3=RAM(J1+3) C Vertex to be separated D=Z1*(RAM(J2+1)-Y1)+Z2*(RAM(J2+2)-Y2)+Z3*(RAM(J2+3)-Y3) DMAX=0. IF(D*HAND.GT.0) THEN C Loop over opposite vertices DO 36 I=1,N IF(I.NE.I1.AND.I.NE.I2.AND.I.NE.I3) THEN J2=6*(IRAM(NTRGL+I)-1) AUX=Z1*(RAM(J2+1)-Y1)+Z2*(RAM(J2+2)-Y2) * +Z3*(RAM(J2+3)-Y3) DMAX=AMAX1(-AUX*SIGN(1.,HAND),DMAX) END IF 36 CONTINUE AUX=AMIN1(ABS(D),DMAX) IF(AUX.GT.0.) THEN AUX=ABS((RAM(J3+4)-RAM(J1+4))*X1 * +(RAM(J3+5)-RAM(J1+5))*X2 * +(RAM(J3+6)-RAM(J1+6))*X3) * *SQRT(X1*X1+X2*X2+X3*X3)/AUX IF(K1.EQ.0.OR.AUX.LT.DMIN) THEN K1=MIN0(I1,I3) K2=MAX0(I1,I3) DMIN=AUX END IF END IF END IF END IF 37 CONTINUE C ^^^^^^^^ C Loop over polygon diagonals to find the best one for division *521 DO 42 I1=1,N-2 * DO 42 I1=1,N * J1=6*(IRAM(NTRGL+I1)-1) C521 Loop over all polygon diagonals *521 DO 41 I2=I1+2,N+MIN0(I1-2,0) C-NEW *** C Limiting the loop to a single diagonal per vertex * I2=MOD(I1+1,N)+1 * J2=6*(IRAM(NTRGL+I2)-1) C Normal Z to the separation plane * X1=RAM(J2+1)-RAM(J1+1) * X2=RAM(J2+2)-RAM(J1+2) * X3=RAM(J2+3)-RAM(J1+3) * Y1=RAM(J2+4)+RAM(J1+4) * Y2=RAM(J2+5)+RAM(J1+5) * Y3=RAM(J2+6)+RAM(J1+6) * Z1=X2*Y3-X3*Y2 * Z2=X3*Y1-X1*Y3 * Z3=X1*Y2-X2*Y1 C Point on the diagonal * Y1=RAM(J1+1) * Y2=RAM(J1+2) * Y3=RAM(J1+3) C Vertex to be separated * I3=MOD(I1,N)+1 * J3=6*(IRAM(NTRGL+I3)-1) * D=Z1*(RAM(J3+1)-Y1)+Z2*(RAM(J3+2)-Y2)+Z3*(RAM(J3+3)-Y3) * DMAX=0. * IF(D.NE.0) THEN C Loop over opposite vertices * DO 35 I=1,N * IF(I.NE.I1.AND.I.NE.I2.AND.I.NE.I3) THEN * J3=6*(IRAM(NTRGL+I)-1) * AUX=Z1*(RAM(J3+1)-Y1)+Z2*(RAM(J3+2)-Y2) * * +Z3*(RAM(J3+3)-Y3) * IF(AUX*D.GE.0.) THEN C Vertex I3 is not separated from opposite vertices * GO TO 41 * END IF * DMAX=AMAX1(ABS(AUX),DMAX) * END IF * 35 CONTINUE * AUX=SQRT((RAM(J2+4)-RAM(J1+4))**2 * * +(RAM(J2+5)-RAM(J1+5))**2 * * +(RAM(J2+6)-RAM(J1+6))**2) * * *(X1*X1+X2*X2+X3*X3)/AMIN1(ABS(D),DMAX) * IF(K1.EQ.0.OR.AUX.LT.DMIN) THEN * K1=MIN0(I1,I2) * K2=MAX0(I1,I2) * DMIN=AUX * END IF * END IF * 41 CONTINUE * 42 CONTINUE C ^^^^^^^^ IF(K1.EQ.0) THEN C TRGL-09 WRITE(TEXT,'(A,7I6)') * 'TRGL-09: Polygon cannot be divided:', * (IRAM(I),I=NTRGL+1,NTRGL+MIN0(7,N)) IF(N.GT.7) THEN TEXT(78:80)='...' END IF CALL ERROR(TEXT(1:LENGTH(TEXT))) C No polygon vertex can be separated from all opposite polygon C vertices by the plane determined by the diagonal connecting C the two adjacent vertices and the sum of the normals in the C adjacent vertices. END IF C v---521---v C Dividing polygon along diagonal between K1-th and K2-th vertices N1=K1+N-K2+1 N2=K2-K1+1 C Doubling K2-th and K1-th vertices DO 43 I=NTRGL+N,NTRGL+K2,-1 IRAM(I+1)=IRAM(I) 43 CONTINUE DO 44 I=NTRGL+N+1,NTRGL+K1,-1 IRAM(I+1)=IRAM(I) 44 CONTINUE C Moving vertices K2,K2+1,...,N C between vertices 1,2,...,K1 and K1,K1+1,...,K2 DO 46 I2=K2,N I=IRAM(NTRGL+N+2) DO 45 I1=NTRGL+N+1,NTRGL+K1+1,-1 IRAM(I1+1)=IRAM(I1) 45 CONTINUE IRAM(NTRGL+K1+1)=I 46 CONTINUE C Moving the second polygon to the proper location DO 47 I=NTRGL+N2,NTRGL+1,-1 IRAM(I+4*(N1-2))=IRAM(I+N1) 47 CONTINUE C Storing zeros between the polygons DO 48 I=NTRGL+4*(N1-2),NTRGL+N1+1,-1 IRAM(I)=0 48 CONTINUE GO TO 30 49 CONTINUE C C Making triangles right-handed with respect to the normals: DO 51 I=NVRTX+1,NPLGN,4 IF(IRAM(I+3).NE.0) THEN C TRGL-07 CALL ERROR('TRGL-07: Triangle not terminated by 0') C This errror should not appear. Contact the author. END IF J1=6*(IRAM(I )-1) J2=6*(IRAM(I+1)-1) J3=6*(IRAM(I+2)-1) X1=RAM(J2+1)-RAM(J1+1) X2=RAM(J2+2)-RAM(J1+2) X3=RAM(J2+3)-RAM(J1+3) Y1=RAM(J3+1)-RAM(J2+1) Y2=RAM(J3+2)-RAM(J2+2) Y3=RAM(J3+3)-RAM(J2+3) Z1=X2*Y3-X3*Y2 Z2=X3*Y1-X1*Y3 Z3=X1*Y2-X2*Y1 IF(Z1.EQ.0..AND.Z2.EQ.0..AND.Z3.EQ.0.) THEN C TRGL-51 WRITE(TEXT,'(A,3I6)') * 'TRGL-51: Triangle has shrunk into a line:', * IRAM(I),IRAM(I+1),IRAM(I+2) CALL WARN(TEXT(1:LENGTH(TEXT))) C The sides of a triangle are parallel. C Marking the triangle linear in order not to write it into the C output file: IRAM(I+3)=-1 ELSE A1=Z1*RAM(J1+4)+Z2*RAM(J1+5)+Z3*RAM(J1+6) A2=Z1*RAM(J2+4)+Z2*RAM(J2+5)+Z3*RAM(J2+6) A3=Z1*RAM(J3+4)+Z2*RAM(J3+5)+Z3*RAM(J3+6) IF(A1.LT.0..AND.A2.LT.0..AND.A3.LT.0.) THEN C Changing left-handed triangle into right-handed one J2=IRAM(I+1) J3=IRAM(I+2) IRAM(I+1)=J3 IRAM(I+2)=J2 ELSE IF(A1.LE.0..OR.A2.LE.0..OR.A3.LE.0.) THEN AUX=X1*(X1+Y1)+Y1*Y1+X2*(X2+Y2)+Y2*Y2+X3*(X3+Y3)+Y3*Y3 AUX=SQRT(Z1*Z1+Z2*Z2+Z3*Z3)/AUX IF(AUX.LE.0.000010) THEN C TRGL-52 WRITE(TEXT,'(A,3I6)') * 'TRGL-52: Triangle is too narrow:', * IRAM(I),IRAM(I+1),IRAM(I+2) CALL WARN(TEXT(1:LENGTH(TEXT))) C Triangle is too narrow, i.e., the area of the triangle C is too small compared with the sum of the squares of its C sides. C Marking the triangle linear in order not to write it into C the output file: IRAM(I+3)=-1 ELSE C TRGL-08 WRITE(TEXT,'(A,3I6)') * 'TRGL-08: Wrong normals at the triangle vertices:', * IRAM(I),IRAM(I+1),IRAM(I+2) CALL ERROR(TEXT(1:LENGTH(TEXT))) C Normals at the triangle vertices do not point at the same C side of the triangle. END IF END IF END IF 51 CONTINUE C C Writing triangles: OPEN(LU,FILE=FTRGL) FORMAT='(3(I0,A))' I=INT(ALOG10(FLOAT(NVRTX)+0.5))+1 FORMAT(5:5)=CHAR(ICHAR('0')+I) DO 61 I=NVRTX+1,NPLGN,4 IF(IRAM(I+3).EQ.0) THEN WRITE(LU,FORMAT) IRAM(I),' ',IRAM(I+1),' ',IRAM(I+2),' /' END IF 61 CONTINUE CLOSE(LU) C IF(NVRTX.EQ.NPLGN) THEN C There are no polygons WRITE(*,'(A)') '+TRGL: No triangles. ' ELSE WRITE(*,'(A)') '+TRGL: Done. ' END IF STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.for C C======================================================================= C