C CONNEC2D.FOR VER.21-FEB-2002 C C PROGRAM FOR CONNECTIVITY ANALYSIS OF A 2D INDICATOR MAP C C GIVEN AN INDICATOR FIELDS (VALUES 0 AND 1 ONLY), THE RANDOM SET C WITH VALUES 1 IS ANALYSED FOR CONNECTIVITY. C C INPUT PARAMETER FILE WITH C C IPHA : 0 OR 1 FOR CONNECTIVITY ANALYSIS OF PHASE 0 OR 1. C ICON : 4 FOR 4-CONNECTIVITY AND 8 FOR 8-CONNECTIVITY C CINP : INPUT FILE WITH INDICATOR VARIABLE (VALUS 0/1 ONLY) C NX NY : NUMBER OF POINTS IN X AND Y C DX DY : GRID DIMENSIONS IN X AND Y C COUT : OUTPUT FILE WITH STATISTICS C COU2 : OUTPUT FILE WITH CONNECTED COMPONENTS C COU3 : OUTPUT FILE WITH CONNECTIVITY FUNCITON C C DX,DY ARE ONLY USED FOR CONVERTING GRID UNITS TO REAL UNITS. C NX : NUMBER OF POINTS ALONG THE X DIRECTION C NY : NUMBER OF POINTS ALONG THE Y DIRECTION C C IND(I,J) : 0/1 INDICATOR VALUE C AT THE TERMINATION OF THE PROGRAM THE MATRIX IND(I,J) CONTAINS C THE CONNECTED COMPONENTS WITH VALUES 1,2,3,4,5,... FOR FIRST, C SECOND, THIRD, ... ETC CONNECTED COMMPONENTS. C THE VALUE 0 REMAINS 0. C COUT : OUTPUT FILE WITH STATISTICS AND CONNECTIVITY FUNCTION. C COU2 : OUTPUT FILE WITH THE CONNECTED COMPONENTS. C COU3 : OUTPUT FILE WITH CONNECTIVITY FUNCTION C C FOR THE CONNECTIVITY FUNCTION: C NPX(K): NUMBER OF PAIRS OF VALUES SEPARATED A DISTANCE K ALONG C THE X DIRECTION AND THAT BELONG TO FACIES 1. C NPCX(K): NUMBER OF THE PREVIOUS VALUES THAT ARE CONNECTED. C NPY(K): NUMBER OF PAIRS OF VALUES SEPARATED A DISTANCE K ALONG C THE Y DIRECTION AND THAT BELONG TO FACIES 1. C NPCY(K): NUMBER OF THE PREVIOUS VALUES THAT ARE CONNECTED. C C C PPHA : PROPORTION OF FACIES 1. C NCC : NUMBER OF CONNECTED COMPONENTS. C RISM : MEAN CONNECTED COMPONENT SIZE IN PIXELS. C RSME : MRAN SIZE IN REAL UNITS. C RTOT : MEAN SIZE RELATIVE TO SIZE OF FACIES 1. C RXME : MEAN LENGTH ALONG THE X DIRECTION. C RYME : MEAN LENGTH ALONG THE Y DIRECTION. C ICOM : NUMBER OF THE LARGEST COMPONENT. C ISMA : SIZE IN PIXELS (OR LARGEST COMPONENT). C RSIZ : SIZE RELATIVE TO SIZE OF FACIES 1. C IXMA : MAXIMUM LENGTH ALONG X. C IYMA : MAXIMUM LENGTH ALONG Y. C ISMI : SIZE OF SMALLEST COMPONENT. C IXMI : MINIMUM LENGTH ALONG X. C IYMI : MINIMUM LENGTH ALONG Y. C IPX : NUMBER OF PERCOLATING COMPONENTS ALONG X. C IPY : NUMBER OF PERCOLATING COMPONENTS ALONG Y. C CCFU(.,.) : CONNECTIVITY FUNCTION. C C CCFU(.,1): CONNECTIVITY FUNCTION ALONG X (E-W). C CCFU(.,2): CONNECTIVITY FUNCTION ALONG Y (N-S). C CCFU(.,3): CONNECTIVITY FUNCTION ALONG NE-SW DIRECTION. C CCFU(.,4): CONNECTIVITY FUNCTION ALONG NW-SE DIRECTION. C CCFU(.,5): MEAN CONNECTIVITY FUNCTION ALONG X AND Y. C CCFU(.,6): MEAN CONNECTIVITY FUNCTION ALONG THE DIAGONALS. C C C ---------------------------------------------------------------- C DIMENSION IND(4096,4096) C,NPX(256),NPY(256),NPCX(256),NPCY(256) C DIMENSION NPE(256),NPW(256),NPCE(256),NPCW(256),CCFU(256,6) C CHARACTER*12 CPAR,CINP,COU2 C C MATRIX DIMENSION LIMITATION C NXM=4096 NYM=4096 C CPAR="coninput.txt" OPEN (1,FILE=CPAR) READ (1,*) IPHA IF (IPHA.NE.0.AND.IPHA.NE.1) THEN WRITE (6,*)'THE FIRST LINE IN THE PARAMETER FILE' WRITE (6,*)'MUST BE A 0 OR A 1' WRITE (6,*)'0 FOR CONNECTIVITY ANALYSIS OF PHASE 0' WRITE (6,*)'1 FOR CONNECTIVITY ANALYSIS OF PHASE 1' WRITE (6,*)'THE ACTUAL VALUE IN THE PARAMETER FILE IS: ',IPHA STOP END IF READ (1,*) ICON IF (ICON.NE.4.AND.ICON.NE.8) THEN WRITE (6,*)'THE SECOND LINE IN THE PARAMETER FILE' WRITE (6,*)'MUST BE A 4 OR A 8' WRITE (6,*)'4 MEANS 4-CONNECTIVITY' WRITE (6,*)'8 MEANS 8-CONNECTIVITY' WRITE (6,*)'THE ACTUAL VALUE IN THE PARAMETER FILE IS: ',ICON STOP END IF READ (1,100) CINP READ (1,*) NX,NY READ (1,*) DX,DY READ (1,100) COU2 CLOSE (1) IF (NX.GT.NXM) THEN WRITE (6,*)'MAXIMUM NX IS ',NXM WRITE (6,*)'ACTUAL VALUE IS ',NX STOP END IF IF (NY.GT.NYM) THEN WRITE (6,*)'MAXIMUM NY IS ',NYM WRITE (6,*)'ACTUAL VALUE IS ',NY STOP END IF C C READING EXPERIMENTAL INDICATOR DATA 2D FIELD C OPEN (1,FILE=CINP) JPHA=0 DO 1 I=1,NX DO 2 J=1,NY READ (1,*) VAL IND(I,J)=INT(VAL) IF (IPHA.EQ.0) THEN IF (INT(VAL).EQ.0) THEN IND(I,J)=1 ELSE IND(I,J)=0 END IF END IF IF (IND(I,J).NE.1.AND.IND(I,J).NE.0) THEN WRITE (6,*)'EXPERIMENTAL DATA MUST BE INDICATOR DATA 0/1' WRITE (6,*)'ACTUAL VALUE : ',IND(I,J) STOP END IF IF (IPHA.EQ.0) THEN IF (INT(VAL).EQ.0) THEN IND(I,J)=1 ELSE IND(I,J)=0 END IF END IF IF (IND(I,J).EQ.1) JPHA=JPHA+1 2 CONTINUE 1 CONTINUE CLOSE (1) IARE=NX*NY PPHA=FLOAT(JPHA)/IARE WRITE (6,*) IF (ICON.EQ.4) WRITE (6,*)'4-CONNECTIVITY ANALYSIS' IF (ICON.EQ.8) WRITE (6,*)'8-CONNECTIVITY ANALYSIS' WRITE (6,*) WRITE (6,*)'PROPORTION OF PHASE 1 IS: ',PPHA C C C LOOKING FOR CONNECTED COMPONENTS C NCC=0 C 5 DO 6 I=1,NX DO 7 J=1,NY IF (IND(I,J).EQ.1) THEN IAI=I IAJ=J GOTO 8 END IF 7 CONTINUE 6 CONTINUE C NO MORE CONNECTED COMPONENTS SENDING CONTROL TO LABEL 20 GOTO 20 C C A NEW COMPONEND HAS BEEN FOUND, INCREASE NUMBER OF C CONNECTED COMPONENTS BY 1 C C 8 NCC=NCC+1 C C LOOKING FOR COMPONENT NCC WITH LABEL NNC+1 C IND(IAI,IAJ)=NCC+1 DO 9 I=IAI,NX DO 10 J=IAJ,NY IF ((I-1).GT.0) THEN IF (IND(I-1,J).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I-1,J)=NCC+1 END IF END IF IF (I+1.LE.NX) THEN IF (IND(I+1,J).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I+1,J)=NCC+1 END IF END IF IF (J-1.GT.0) THEN IF (IND(I,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I,J-1)=NCC+1 END IF END IF IF (J+1.LE.NY) THEN IF (IND(I,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I,J+1)=NCC+1 END IF END IF IF (ICON.EQ.8) THEN IF ((I-1).GT.0.AND.(J-1).GT.0) THEN IF (IND(I-1,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I-1,J-1)=NCC+1 END IF END IF IF ((I+1).LE.NX.AND.(J-1).GT.0) THEN IF (IND(I+1,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I+1,J-1)=NCC+1 END IF END IF IF ((I-1).GT.0.AND.(J+1).LE.NY) THEN IF (IND(I-1,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I-1,J+1)=NCC+1 END IF END IF IF ((I+1).LE.NX.AND.(J+1).LE.NY) THEN IF (IND(I+1,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I+1,J+1)=NCC+1 END IF END IF END IF 10 CONTINUE 9 CONTINUE DO 40 I=IAI,1,-1 DO 41 J=IAJ,1,-1 IF ((I-1).GT.0) THEN IF (IND(I-1,J).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I-1,J)=NCC+1 END IF END IF IF (I+1.LE.NX) THEN IF (IND(I+1,J).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I+1,J)=NCC+1 END IF END IF IF (J-1.GT.0) THEN IF (IND(I,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I,J-1)=NCC+1 END IF END IF IF (J+1.LE.NY) THEN IF (IND(I,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I,J+1)=NCC+1 END IF END IF IF (ICON.EQ.8) THEN IF ((I-1).GT.0.AND.(J-1).GT.0) THEN IF (IND(I-1,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I-1,J-1)=NCC+1 END IF END IF IF ((I+1).LE.NX.AND.(J-1).GT.0) THEN IF (IND(I+1,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I+1,J-1)=NCC+1 END IF END IF IF ((I-1).GT.0.AND.(J+1).LE.NY) THEN IF (IND(I-1,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I-1,J+1)=NCC+1 END IF END IF IF ((I+1).LE.NX.AND.(J+1).LE.NY) THEN IF (IND(I+1,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) IND(I+1,J+1)=NCC+1 END IF END IF END IF 41 CONTINUE 40 CONTINUE C 62 NNN=0 DO 60 I=1,NX DO 61 J=1,NY IF ((I-1).GT.0) THEN IF (IND(I-1,J).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I-1,J)=NCC+1 NNN=NNN+1 END IF END IF END IF IF (I+1.LE.NX) THEN IF (IND(I+1,J).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I+1,J)=NCC+1 NNN=NNN+1 END IF END IF END IF IF (J-1.GT.0) THEN IF (IND(I,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I,J-1)=NCC+1 NNN=NNN+1 END IF END IF END IF IF (J+1.LE.NY) THEN IF (IND(I,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I,J+1)=NCC+1 NNN=NNN+1 END IF END IF END IF IF (ICON.EQ.8) THEN IF ((I-1).GT.0.AND.(J-1).GT.0) THEN IF (IND(I-1,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I-1,J-1)=NCC+1 NNN=NNN+1 END IF END IF END IF IF ((I+1).LE.NX.AND.(J-1).GT.0) THEN IF (IND(I+1,J-1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I+1,J-1)=NCC+1 NNN=NNN+1 END IF END IF END IF IF ((I-1).GT.0.AND.(J+1).LE.NY) THEN IF (IND(I-1,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I-1,J+1)=NCC+1 NNN=NNN+1 END IF END IF END IF IF ((I+1).LE.NX.AND.(J+1).LE.NY) THEN IF (IND(I+1,J+1).EQ.1) THEN IF (IND(I,J).EQ.NCC+1) THEN IND(I+1,J+1)=NCC+1 NNN=NNN+1 END IF END IF END IF END IF 61 CONTINUE 60 CONTINUE IF (NNN.GT.0) GOTO 62 GOTO 5 C C THE CONNECTED COMPONENT IS THE COMPONENT WITH VALUES C IA(I,J)=NCC, WE SUBSTRACT 1 FROM THE LEBEL C 20 DO 50 I=1,NX DO 51 J=1,NY IF (IND(I,J).NE.0) IND(I,J)=IND(I,J)-1 51 CONTINUE 50 CONTINUE C C C OUTPUT FILE WITH CONNECTED COMPONENTS C OPEN (1,FILE=COU2) C OPEN (2,FILE=CINP) C DO 15 I=1,NX DO 16 J=1,NY C READ (2,*) ZVA WRITE (1,*) IND(I,J) 16 CONTINUE 15 CONTINUE C CLOSE (1) C 100 FORMAT (A) 110 FORMAT (1X,'INPUT INDICATOR FILE : ',A12) 120 FORMAT (1X,'NX : ',I6) 130 FORMAT (1X,'NY : ',I6) 140 FORMAT (1X,'DX : ',F10.3) 150 FORMAT (1X,'DY : ',F10.3) 160 FORMAT (1X,'PROPORTION OF FACIES ',I1,' : ',F10.4) 170 FORMAT (/1X,'CONNECTED COMPONENTS STATISTICS'// * /1X,'NUMBER OF CONNECTED COMPONENTS : ',I6) 175 FORMAT (/1X,'AVERAGES'/1X,8('=')) 176 FORMAT (/1X,'MINIMA'/1X,6('=')) 177 FORMAT (/1X,'MAXIMA'/1X,6('=')) 178 FORMAT (/1X,'PERCOLATION'/1X,11('=')) 180 FORMAT (1X,'MEAN SIZE IN PIXELS : ',F10.4) 190 FORMAT (1X,'MEAN SIZE REAL UNITS : ',F12.4) 200 FORMAT (1X,'MEAN SIZE RELATIVE TO TOTAL AREA OF FACIES ',I1, * ' : ',F10.4) 210 FORMAT (1X,'MEAN LENGTH ALONG X (IN PIXELS) : ',F10.4) 220 FORMAT (1X,'MINIMUM LENGTH ALONG X (IN PIXELS) : ',I6) 230 FORMAT (1X,'MAXIMUM LENGTH ALONG X (IN PIXELS) : ',I6) 235 FORMAT (1X,'THE LARGEST COMPONENT IS NUMBER : ',I6) 236 FORMAT (1X,'WITH MAXIMUM SIZE IN PIXELS : ',I6) 237 FORMAT (1X,'AND RELATIVE TO TOTAL AREA OF FACIES ' * ,I1,' : ',F10.4) 238 FORMAT (1X,'SIZE IN PIXELS OF SMALLEST COMPONENT : ',I6) 240 FORMAT (1X,'MEAN LENGTH ALONG Y (IN PIXELS) : ',F10.4) 250 FORMAT (1X,'MIMIMUN LENGTH ALONG Y (IN PIXELS) : ',I6) 260 FORMAT (1X,'MAXIMUM LENGTH ALONG Y (IN PIXELS) : ',I6) 270 FORMAT (1X,'NUMBER OF PERCOLATING COMPONENTS IN X : ',I6) 280 FORMAT (1X,'NUMBER OF PERCOLATING COMPONENTS IN Y : ',I6) 290 FORMAT (/1X,'OUTPUT FILE WITH CONNECTED COMPONENTS : ',A12) 295 FORMAT (/1X,'OUTPUT FILE WITH CONNECTIVIY FUNCTION : ',A12) 300 FORMAT (/1X,'CONNECTIVITY FUNCTION'/1X,21('=') * //1X,'ALONG THE X DIRECTION') 310 FORMAT (F10.3,F12.6,2I12) 320 FORMAT (/1X,'ALONG THE Y DIRECTION') 321 FORMAT (/1X,'ALONG THE NE-SW DIRECTION') 322 FORMAT (/1X,'ALONG THE NW-SE DIRECTION') 323 FORMAT (/1X,'AVERAGE ALONG THE DIAGONALS') 330 FORMAT (/1X,'AVERAGE ALONG X AND Y') 400 FORMAT (25(/),1X,'CONNEC2D PROGRAM VER. 1.0'/1X,27('=')/// * 1X,'INPUT PARAMETER FILE ---> ',$) 410 FORMAT (/1X,'OUTPUT FILE WITH STATISTICS AND CONNECTIVITY ' * 'FUNCTION: ',A12) 420 FORMAT (/1X,'OUTPUT FILE WITH CONNECTED COMPONENTS: ',A12) 425 FORMAT (/1X,'OUTPUT FILE WITH CONNECTIVITY FUNCTION: ',A12) 430 FORMAT (1X,'OUTPUT RESULTS OF CONNEC2D'//1X'4-CONNECTIVITY' * ' ANALYSIS'/1X,'PARAMETER FILE WAS : ',A12//) 440 FORMAT (1X,'OUTPUT RESULTS OF CONNEC2D'//1X'4-CONNECTIVITY' * ' ANALYSIS'/1X,'PARAMETER FILE WAS : ',A12//) END