/************************************************************************/
/* */
/* Program: rand3d.c */
/* Purpose: To read in a 3-D image and output */
/* a filtered image in which one phase from the original */
/* image has been separated into two distinct phases */
/* Programmer: Dale P. Bentz */
/* NIST */
/* 100 Bureau Drive Mail Stop 8621 */
/* Gaithersburg, MD 20899-0001 */
/* Phone: (301) 975-5865 */
/* E-mail: dale.bentz@.nist.gov */
/* */
/* Created: From Fortran version 2/00 */
/************************************************************************/
#include <stdio.h>
#include <math.h>
/* Random number generator */
#define IA 16807
#define IM 2147483647
#define IQ 127773
#define IR 2836
#define NTAB 32
#define EPS (1.2E-07)
#define MAX(a,b) (a>b)?a:b
#define MIN(a,b) (a<b)?a:b
double ran1(idum)
int *idum;
{
int j,k;
static int iv[NTAB],iy=0;
void nrerror();
static double NDIV = 1.0/(1.0+(IM-1.0)/NTAB);
static double RNMX = (1.0-EPS);
static double AM = (1.0/IM);
if ((*idum <= 0) || (iy == 0)) {
*idum = MAX(-*idum,*idum);
for(j=NTAB+7;j>=0;j--) {
k = *idum/IQ;
*idum = IA*(*idum-k*IQ)-IR*k;
if(*idum < 0) *idum += IM;
if(j < NTAB) iv[j] = *idum;
}
iy = iv[0];
}
k = *idum/IQ;
*idum = IA*(*idum-k*IQ)-IR*k;
if(*idum<0) *idum += IM;
j = iy*NDIV;
iy = iv[j];
iv[j] = *idum;
return MIN(AM*iy,RNMX);
}
#undef IA
#undef IM
#undef IQ
#undef IR
#undef NTAB
#undef EPS
#undef MAX
#undef MIN
#define SYSIZE 100
main(){
int *seed,iseed,syssize;
float snew,s2,ss,sdiff,phasein,phaseout,xtmp,ytmp;
static float normm[SYSIZE+1][SYSIZE+1][SYSIZE+1];
static float res[SYSIZE+1][SYSIZE+1][SYSIZE+1];
static int mask[SYSIZE+1][SYSIZE+1][SYSIZE+1];
static float filter [32][32][32];
int done,r[61];
static float s[61],xr[61],sum[502];
float val2;
float t1,t2,x1,x2,u1,u2,xrad,resmax,resmin;
float pi,xtot,filval,radius,xpt,sect,sumtot,vcrit;
int valin,r1,r2,i1,i2,i3,i,j,k,j1,k1;
int ido,iii,jjj,ix,iy,iz,index;
char filen[80],filem[80];
FILE *infile,*outfile;
/* Get user input */
pi=3.1415926;
printf("Enter random number seed (negative integer) \n");
scanf("%d",&iseed);
printf("%d\n",iseed);
seed=(&iseed);
printf("Enter existing phase assignment for matching\n");
scanf("%f",&phasein);
printf("%f\n",phasein);
printf("Enter phase assignment to be created by program\n");
scanf("%f",&phaseout);
printf("%f\n",phaseout);
printf("Enter name of cement microstructure image file\n");
scanf("%s",filen);
printf("%s\n",filen);
/* Read in the original microstructure image file */
infile=fopen(filen,"r");
for(k=1;k<=SYSIZE;k++){
for(j=1;j<=SYSIZE;j++){
for(i=1;i<=SYSIZE;i++){
fscanf(infile,"%d",&valin);
mask[i][j][k]=valin;
}
}
}
fclose(infile);
/* Create the Gaussian noise image 2 elements at a time */
i1=i2=i3=1;
for(i=1;i<=500000;i++){
u1=ran1(seed);
u2=ran1(seed);
t1=2.*pi*u2;
t2=sqrt(-2.*log(u1));
x1=cos(t1)*t2;
x2=sin(t1)*t2;
normm[i1][i2][i3]=x1;
i1+=1;
if(i1>SYSIZE){
i1=1;
i2+=1;
if(i2>SYSIZE){
i2=1;
i3+=1;
}
}
normm[i1][i2][i3]=x2;
i1+=1;
if(i1>SYSIZE){
i1=1;
i2+=1;
if(i2>SYSIZE){
i2=1;
i3+=1;
}
}
}
/* Now perform the convolution of */
/* the Gaussian noise image with the correlation filter */
printf("Enter filename to read in autocorrelation from\n");
scanf("%s",filen);
printf("%s\n",filen);
infile=fopen(filen,"r");
fscanf(infile,"%d",&ido);
printf("Number of points in correlation file is %d \n",ido);
for(i=1;i<=ido;i++){
fscanf(infile,"%d %f",&valin,&val2);
r[i]=valin;
s[i]=val2;
xr[i]=(float)r[i];
}
fclose(infile);
ss=s[1];
s2=ss*ss;
/* Load up the convolution (filter) matrix with the correlation function*/
sdiff=ss-s2;
for(i=0;i<31;i++){
iii=i*i;
for(j=0;j<31;j++){
jjj=j*j;
for(k=0;k<31;k++){
xtmp=(float)(iii+jjj+k*k);
radius=sqrt(xtmp);
r1=(int)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;
}
}
}
/* Now filter the image maintaining periodic boundaries */
resmax=0.0;
resmin=1.0;
for(i=1;i<=SYSIZE;i++){
for(j=1;j<=SYSIZE;j++){
for(k=1;k<=SYSIZE;k++){
res[i][j][k]=0.0;
if((float)mask[i][j][k]==phasein){
for(ix=1;ix<=31;ix++){
i1=i+ix-1;
if(i1<1){i1+=SYSIZE;}
else if(i1>SYSIZE){i1-=SYSIZE;}
for(iy=1;iy<=31;iy++){
j1=j+iy-1;
if(j1<1){j1+=SYSIZE;}
else if(j1>SYSIZE){j1-=SYSIZE;}
for(iz=1;iz<=31;iz++){
k1=k+iz-1;
if(k1<1){k1+=SYSIZE;}
else if(k1>SYSIZE){k1-=SYSIZE;}
res[i][j][k]+=normm[i1][j1][k1]*filter[ix][iy][iz];
}
}
}
if(res[i][j][k]>resmax){resmax=res[i][j][k];}
if(res[i][j][k]<resmin){resmin=res[i][j][k];}
}
}
}
}
/* Now threshold the image */
printf("Input desired threshold phase fraction\n");
scanf("%f",&xpt);
printf("%f\n",xpt);
/* Divide the filtered image values into 500 separate bins */
sect=(resmax-resmin)/500.;
printf("SECT is %f\n",sect);
printf("RESMAX and RESMIN are %f and %f\n",resmax,resmin);
for(i=1;i<=500;i++){
sum[i]=0.0;
}
xtot=0.0;
for(i=1;i<=SYSIZE;i++){
for(j=1;j<=SYSIZE;j++){
for(k=1;k<=SYSIZE;k++){
if((float)mask[i][j][k]==phasein){
xtot+=1;
index=1+(int)((res[i][j][k]-resmin)/sect);
if(index>500){index=500;}
sum[index]+=1.0;
}
}
}
}
/* Determine which bin to choose for correct thresholding */
sumtot=vcrit=0.0;
done=0;
for(i=1;((i<=500)&&(done==0));i++){
sumtot+=sum[i]/xtot;
if(sumtot>xpt){
ytmp=(float)i;
vcrit=resmin+(resmax-resmin)*(ytmp-0.5)/500.;
done=1;
}
}
printf("Critical volume fraction is %f\n",vcrit);
/* Write out the final resultant image */
printf("Enter name of new cement microstructure image file\n");
scanf("%s",filem);
printf("%s\n",filem);
outfile=fopen(filem,"w");
for(k=1;k<=SYSIZE;k++){
for(j=1;j<=SYSIZE;j++){
for(i=1;i<=SYSIZE;i++){
if((float)mask[i][j][k]==phasein){
if(res[i][j][k]>vcrit){
res[i][j][k]=phaseout;
}
else{
res[i][j][k]=phasein;
}
}
else{
res[i][j][k]=(float)mask[i][j][k];
}
fprintf(outfile,"%2d\n",(int)res[i][j][k]);
}
}
}
fclose(outfile);
}