You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
428 lines
15 KiB
Fortran
428 lines
15 KiB
Fortran
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(256,256),NPX(256),NPY(256),NPCX(256),NPCY(256)
|
|
DIMENSION NPE(256),NPW(256),NPCE(256),NPCW(256),CCFU(256,6)
|
|
C
|
|
CHARACTER*12 CPAR,CINP,COUT,COU2,COU3
|
|
C
|
|
C MATRIX DIMENSION LIMITATION
|
|
C
|
|
NXM=256
|
|
NYM=256
|
|
C
|
|
WRITE (6,400)
|
|
READ (5,100) CPAR
|
|
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) COUT
|
|
READ (1,100) COU2
|
|
READ (1,100) COU3
|
|
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,*) XX,YY,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)
|
|
OPEN (2,FILE=CINP)
|
|
C
|
|
DO 15 I=1,NX
|
|
DO 16 J=1,NY
|
|
READ (2,*) XX,YY,ZVA
|
|
WRITE (1,*) XX,YY,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
|