Next: Listing for stat3d.c
Up: Computer programs for
Previous: Listing for genpart3d.c
C
C Program RAND3D.F: to create a 100*100*100
C microstructure based on filtering Gaussian
C noise with autocorrelation filter of a 2-D
C image
C Programmer: Dale Bentz (dale.bentz@nist.gov)
C Date: 1993
C 1995: Modified to filter random particle images
C for cement hydration model
C
INTEGER SEED,SYSSIZE,IRES
REAL SNEW,S2,SS,SDIFF,PHASEIN,PHASEOUT,XTMP,YTMP
REAL NORMM(100,100,100),RES(100,100,100)
REAL MASK(100,100,100)
REAL FILTER(31,31,31)
INTEGER R(60)
REAL S(60),XR(60),SUM(501)
REAL T1,T2,X1,X2,U1,U2,XRAD,RESMAX,RESMIN
INTEGER R1,R2,I1,I2,I3,I,J,K,J1,K1
INTEGER IDO,III,JJJ,IX,IY,IZ,INDEX
REAL PI,XTOT,FILVAL,RADIUS,XPT,SECT,SUMTOT,VCRIT
CHARACTER*80 FILEN,FILEM
PI=3.1415926
SYSSIZE=100
WRITE(6,*)'Input random seed'
READ(5,*)SEED
WRITE(6,*)SEED
C
C In this version, only a portion of the original 3-D image is filtered
C so that for instance, silicates may be separated into C3S and C2S
C without modification of the aluminate and gypsum phases
C
WRITE(6,*)'Input existing phase assignment for matching'
READ(5,*)PHASEIN
WRITE(6,*)PHASEIN
WRITE(6,*)'Input phase assignment to be created'
READ(5,*)PHASEOUT
WRITE(6,*)PHASEOUT
C
C Read in mask image (particles) from file
C
WRITE(6,*)'Enter name of cement microstructure image file'
READ(5,*)FILEM
OPEN(13,FILE=FILEM,STATUS='OLD')
DO 19 I=1,SYSSIZE
DO 19 J=1,SYSSIZE
DO 19 K=1,SYSSIZE
19 READ(13,*)MASK(I,J,K)
CLOSE(13)
C
C Create the gaussian noise image
C by appropriately modifying uniform noise
C image
C Source: Law, A.M, and Kelton, W.D.,
C Simulation Modeling and Analysis, McGraw-Hill, 1982.
C
I1=1
I2=1
I3=1
DO 70 I=1,500000
U1=RAN1(SEED)
U2=RAN1(SEED)
T1=2.*PI*U2
T2=(-2.*LOG(U1))**0.5
X1=COS(T1)*(T2)
X2=SIN(T1)*(T2)
NORMM(I1,I2,I3)=X1
I1=I1+1
IF(I1.GT.SYSSIZE) THEN
I1=1
I2=I2+1
IF(I2.GT.SYSSIZE) THEN
I2=1
I3=I3+1
ENDIF
ENDIF
NORMM(I1,I2,I3)=X2
I1=I1+1
IF(I1.GT.SYSSIZE) THEN
I1=1
I2=I2+1
IF(I2.GT.SYSSIZE) THEN
I2=1
I3=I3+1
ENDIF
ENDIF
70 CONTINUE
CLOSE(14)
C
C Now convolve with the autocorrelation function of the goal 2-D image
C
WRITE(6,*)'Enter filename to read in autocorrelation from'
READ(5,*)FILEN
OPEN(8,FILE=FILEN,STATUS='OLD')
READ(8,*)IDO
WRITE(6,*)IDO
DO 95 I=1,IDO
READ(8,*)R(I),S(I)
XR(I)=R(I)
95 CONTINUE
CLOSE(8)
SS=S(1)
S2=SS*SS
C
C Load up the convolution matrix (31*31*31 in size)
C
OPEN(10,FILE='FILTER3D.IMG')
SDIFF=SS-S2
DO 14 I=0,30
III=I*I
DO 13 J=0,30
JJJ=J*J
DO 12 K=0,30
XTMP=III+JJJ+K*K
C
C Use Linear interpolation to estimate S(x,y,z) from available S(r)
C
RADIUS=SQRT(XTMP)
R1=RADIUS+1
R2=R1+1
XRAD=RADIUS+1-R1
FILVAL=S(R1)+(S(R2)-S(R1))*XRAD
FILTER(I+1,J+1,K+1)=(FILVAL-S2)/(SDIFF)
12 WRITE(10,*)FILTER(I+1,J+1,K+1)
13 CONTINUE
14 CONTINUE
CLOSE(10)
C
C Now filter the image maintaining periodic boundaries
C
C Store minimum (RESMIN) and maximum (RESMAX) values for later binning
C
RESMAX=0.0
RESMIN=1.0
DO 97 I=1,SYSSIZE
DO 97 J=1,SYSSIZE
DO 97 K=1,SYSSIZE
RES(I,J,K)=0.0
DO 98 IX=1,31
I1=I+IX-1
C
C Periodic boundaries
C
IF(I1.LT.1) THEN
I1=I1+SYSSIZE
ELSE IF(I1.GT.SYSSIZE) THEN
I1=I1-SYSSIZE
ENDIF
DO 98 IY=1,31
J1=J+IY-1
IF(J1.LT.1) THEN
J1=J1+SYSSIZE
ELSE IF(J1.GT.SYSSIZE) THEN
J1=J1-SYSSIZE
ENDIF
DO 98 IZ=1,31
K1=K+IZ-1
IF(K1.LT.1) THEN
K1=K1+SYSSIZE
ELSE IF(K1.GT.SYSSIZE) THEN
K1=K1-SYSSIZE
ENDIF
RES(I,J,K)=RES(I,J,K)+NORMM(I1,J1,K1)*FILTER(IX,IY,IZ)
98 CONTINUE
C
C Check for min/max if pixel is part of proper phase
C
IF(MASK(I,J,K).EQ.PHASEIN) THEN
IF(RES(I,J,K).GT.RESMAX) THEN
RESMAX=RES(I,J,K)
ENDIF
IF(RES(I,J,K).LT.RESMIN) THEN
RESMIN=RES(I,J,K)
ENDIF
ENDIF
97 CONTINUE
C
C Now threshold the image
C
WRITE(6,*)'Input desired threshold phase fraction'
READ(5,*)XPT
WRITE(6,*)XPT
C
C Bin the resultant (filtered) values into 500 bins
C so that adequate binarization (thresholding) can be performed
C
SECT=(RESMAX-RESMIN)/500.0
WRITE(6,*)'SECT is',SECT
WRITE(6,*)'RESMAX and RESMIN are ',RESMAX,RESMIN
DO 34 I=1,500
34 SUM(I)=0.0
C
C Fill in the 500-bin histogram
C
XTOT=0.0
DO 33 I=1,SYSSIZE
DO 33 J=1,SYSSIZE
DO 33 K=1,SYSSIZE
IF(MASK(I,J,K).EQ.PHASEIN) THEN
XTOT=XTOT+1.0
INDEX=1+(RES(I,J,K)-RESMIN)/SECT
IF(INDEX.GT.500) THEN
INDEX=500
ENDIF
SUM(INDEX)=SUM(INDEX)+1.0
ENDIF
33 CONTINUE
C
C Determine which bin to choose for correct thresholding
C
SUMTOT=0.0
DO 35 I=1,500
SUMTOT=SUMTOT+SUM(I)/(XTOT)
IF(SUMTOT.GT.XPT) THEN
YTMP=I
VCRIT=RESMIN+(RESMAX-RESMIN)*(YTMP-0.5)/500.
GOTO 36
ENDIF
35 CONTINUE
36 WRITE(6,*)'VCRIT is',VCRIT
IRES=0
C
C Perform the segmentation and output the resultant image
C
OPEN(9,FILE='IMAGE3D.OUT')
DO 32 I=1,SYSSIZE
DO 32 J=1,SYSSIZE
DO 32 K=1,SYSSIZE
C
C Only set values that were originally the phase of choice
C
IF(MASK(I,J,K).EQ.PHASEIN) THEN
IF(RES(I,J,K).GT.VCRIT) THEN
RES(I,J,K)=PHASEOUT
ELSE
RES(I,J,K)=PHASEIN
IRES=IRES+1
ENDIF
ELSE
RES(I,J,K)=MASK(I,J,K)
ENDIF
WRITE(9,*)INT(RES(I,J,K))
32 CONTINUE
CLOSE(9)
WRITE(6,*)'Solids are ',FLOAT(IRES)/1000000.
STOP
END
FUNCTION RAN1(IDUM)
C
C Portable random number generator, RAN1
C To generate uniform random deviates between 0 and 1.0
C From: Computers in Physics, Vol. 6, (5), 1992, 522-524
C Press and Teukolsky
C
C CAll with idum set to a negative number to initialize
C RNMAX should approximate the largest floating value that
C is less than 1.0
INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV
REAL RAN1,AM,EPS,RNMAX
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
+ NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2E-7,RNMAX=1.-EPS)
INTEGER J,K,IV(NTAB),IY
SAVE IV,IY
DATA IV /NTAB*0/, IY /0/
IF(IDUM.LE.0.OR.IY.EQ.0) THEN
IF(IDUM.LE.0) THEN
IDUM=MAX(-IDUM,1)
ENDIF
DO 11 J=NTAB+8,1,-1
K=IDUM/IQ
IDUM=IA*(IDUM-K*IQ)-IR*K
IF(IDUM.LT.0) IDUM=IDUM+IM
IF(J.LE.NTAB) IV(J)=IDUM
11 CONTINUE
IY=IV(1)
ENDIF
K=IDUM/IQ
IDUM=IA*(IDUM-K*IQ)-IR*K
IF(IDUM.LT.0) IDUM=IDUM+IM
J=1+IY/NDIV
IY=IV(J)
IV(J)=IDUM
RAN1=MIN(AM*IY,RNMAX)
RETURN
END