/************************************************************************/
/*                                                                      */
/*         Program: clusterid.c                                         */
/*         Purpose: to identify clusters ("particles") in a 3-D         */
/*                  segmented volume and remove those smaller than      */
/*                  a user-specified volume                             */
/*                  Only applicable to a cubic (X * X * X) volume       */ 
/*         Input file: segemented 3-D image, one integer per line per   */
/*                      voxel						*/
/*         Output file: binary 3-D image, with values of 1 for particle */
/*                      (cluster) voxels and 0 otherwise		*/
/*         Programmer: Dale P. Bentz                                    */
/*                     National Institute of Standards and Technology   */
/*                     100 Bureau Drive Stop 8621                       */
/*                     Gaithersburg, MD 20899-8621                      */
/*                     Phone: (301) 975-5865                            */
/*                     Fax: (301) 990-6891                              */
/*                     E-mail: dale.bentz@nist.gov                      */
/*                                                                      */
/*         Date: 2001                                                   */
/*                                                                      */
/************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define SYSSIZE 300
#define SYSSIZEP1 301    /* SYSSIZE plus one */
#define NUMPARTS 50000  /* maximum number of possible clusters */
#define POROSITY 0
#define ANHYDROUS 1
#define PRODUCTS 2
#define BURNT 8000     /* indicates burnt pixel during cluster identification */

/* array phase stores the microstructure representation */
static unsigned short int phase [SYSSIZEP1] [SYSSIZEP1] [SYSSIZEP1];
static int particle [SYSSIZEP1] [SYSSIZEP1] [SYSSIZEP1];
long int npixels[NUMPARTS];
short int pkeep[NUMPARTS];

/* routine to count phase fractions (porosity and solids) */
void phcount()
{
	long int npore,nanhydrous,nproducts;
	int ix,iy,iz;

	npore=0;
	nanhydrous=nproducts=0;
	/* check all voxels in the 3-D system */
	for(iz=0;iz<SYSSIZE;iz++){
	for(iy=0;iy<SYSSIZE;iy++){
	for(ix=0;ix<SYSSIZE;ix++){

		if(phase [ix] [iy] [iz]==POROSITY){
			npore+=1;
		}
		else if(phase[ix][iy][iz]==ANHYDROUS){
			nanhydrous+=1;
		}
		else if(phase[ix][iy][iz]==PRODUCTS){
			nproducts+=1;
		}
	}
	}
	}

	printf("Pores are: %ld \n",npore);
	printf("Anhydrous are: %ld \n",nanhydrous);
	printf("Products are: %ld \n",nproducts);
}


/*	routine burncer to find "separate" particles in 3-D image */
/*      particles have a phase ID input by user */ 
void burncer()
{
	int ix1,iy1,iz1,phkey;
	int ix,iy,iz,iw;
	int ic,ixx,iyy,izz,ncount,n1,n2;
	int nmat1;
	long int xsum,ysum,zsum;
	float xcent,ycent,zcent,save,ssum,s2sum;
	int xl1 [8000],yl1 [8000],zl1 [8000],xl2 [8000],yl2 [8000],zl2 [8000];
	long int bcount,ccount,tcount,iclust;

	printf("Enter phase ID (1- anhydrous, 2-products, etc.) for particles \n");
        scanf("%d",&phkey);
        printf("%d\n",phkey);

/*	bcount tallies pixels which are accessible from source */
	bcount=0;
	tcount=0;
	iclust=0; /* counter for number of partciles found */

/*	check all pixels for particles */
	
	for(iz=0;iz<SYSSIZE;iz++){
	for(iy=0;iy<SYSSIZE;iy++){
	for(ix=0;ix<SYSSIZE;ix++){

	tcount+=1;
	ccount=0;
	ncount=0;
	nmat1=0;
        xsum=ysum=zsum=0;

/*	if current pixel is particle, then flood at this site */
	if(phase [ix] [iy] [iz]==phkey){
		iw=1;
		iclust+=1;
                if(iclust>=NUMPARTS){
                  printf("Too many particles, increase value of NUMPARTS! \n");
                  exit(1);
                }
                particle [ix] [iy] [iz]=iclust;

/*	store current pixel location  */
		nmat1+=1;
		xsum+=ix;
		ysum+=iy;
		zsum+=iz;
		xl1[nmat1]=ix;
		yl1[nmat1]=iy;
		zl1[nmat1]=iz;

/*       convert pixel ID to burnt */
		phase [ix] [iy] [iz]=BURNT;
		ncount+=1;
		ccount+=1;
		bcount+=1;

/* while new sites reamin for propagation */
	while(ncount>0){

		n2=0;
		nmat1=0;
/*	check all new site pixels for possible flood propagation */
	for(n1=1;n1<=ncount;n1++){

/*	read in location of next pixel site to be checked */
	if(iw==1){
		ixx=xl1[n1];
		iyy=yl1[n1];
		izz=zl1[n1];
		}
	if(iw==2){
		ixx=xl2[n1];
		iyy=yl2[n1];
		izz=zl2[n1];
		}

/*	check 6 neighbors for flood propagation  */
	for(ic=1;ic<=6;ic++){
		ix1=ixx;
		iy1=iyy;
		iz1=izz;

	switch (ic) {
		case 1:
			ix1=ixx-1;
			break;
		case 2:
			ix1=ixx+1;
			break;
		case 3:
			iy1=iyy-1;
			break;
		case 4:
			iy1=iyy+1;
			break;
		case 5:
			iz1=izz-1;
			break;
		case 6:
			iz1=izz+1;
			break;
		default:
			break;
			}

/*	if new site is in original box, then */
/*	check for propagation	*/
/*      note- no periodic boundaries are assumed */
	if((ix1>=0)&&(ix1<SYSSIZE)&&(iy1>=0)&&(iy1<SYSSIZE)&&(iz1>=0)&&(iz1<SYSSIZE)){

/*	if flow is into particle space, store location in new */
/*	source array					*/
	if(phase [ix1] [iy1] [iz1]==phkey){
		n2+=1;
		if(n2>=8000){
			printf("Error in burning too many new sources \n");
			exit(1);
		}
		ccount+=1;
		bcount+=1;
		xsum+=ix1;
		ysum+=iy1;
		zsum+=iz1;
                particle[ix1][iy1][iz1]=iclust;
		if(iw==1){
			xl2[n2]=ix1;
			yl2[n2]=iy1;
			zl2[n2]=iz1;
			}
		if(iw==2){
			xl1[n2]=ix1;
			yl1[n2]=iy1;
			zl1[n2]=iz1;
			}
		phase [ix1] [iy1] [iz1]=BURNT;
		}
		}

	}
	}
/*	get ready for next cycle of the propagation	*/
	ncount=n2;
	iw+=1;
	if(iw>2){iw=1;}
	
	}

/*      Compute particle centroid */
	xcent=(float)xsum/(float)ccount;
	ycent=(float)ysum/(float)ccount;
	zcent=(float)zsum/(float)ccount;
	printf("Cluster %ld  Size %ld \n",iclust,ccount);
	printf("Source= (%d, %d, %d) \n",ix,iy,iz);
	printf("Centroid= (%f, %f, %f) \n",xcent,ycent,zcent);
        npixels[iclust]=ccount;
	/* Tabulate statistics for particle greater than 2 voxels in volume */
	if(ccount>=2){
		ssum+=(float)ccount;
		s2sum+=((float)ccount*(float)ccount);
		}

		}

	}
	}
	}

	printf("Total pixels checked- %ld  Clusters found- %ld  \n",tcount,bcount);
	save=0.0;
	if(ssum>0.0){
		save=s2sum/ssum;
	}
	printf("ssum count is %f \n",ssum);
	printf("Site weighted average cluster size- %f \n",save);

/*	Reconvert burnt areas to original phase */
	for(iz=0;iz<SYSSIZE;iz++){
	for(iy=0;iy<SYSSIZE;iy++){
	for(ix=0;ix<SYSSIZE;ix++){
		if(phase [ix] [iy] [iz]==BURNT){
			phase [ix] [iy] [iz]=phkey;
		}

	}
	}
	}

}

void readmic()
{
	FILE *infile;
	char filen[80];
	int iout,i1,i2,i3;
	long int porcnt,totcnt,gelcnt,solidcnt;

	porcnt=gelcnt=totcnt=solidcnt=0;
	printf("Enter name of file to read in \n");
	scanf("%s",filen);
	infile=fopen(filen,"r");
	
	for(i3=0;i3<SYSSIZE;i3++){
	for(i2=0;i2<SYSSIZE;i2++){
	for(i1=0;i1<SYSSIZE;i1++){
		fscanf(infile,"%d",&iout);
		totcnt+=1;
		phase [i1] [i2] [i3]=iout;
		if(phase [i1] [i2] [i3]==POROSITY){
			porcnt+=1;
		}
		if(phase [i1] [i2] [i3]==ANHYDROUS){
			solidcnt+=1;
		}
		if(phase [i1] [i2] [i3]==PRODUCTS){
			gelcnt+=1;
		}
                particle[i1][i2][i3]=0;
	}
	}
	}
	printf("Count for porosity is %ld out of %ld \n",porcnt,totcnt);
	printf("Count for products is %ld and for anhydrous %ld \n",gelcnt,solidcnt);
	fclose(infile);
}

void savemic()
{
	FILE *outfile;
        char fileout[80];
	int iout,i1,i2,i3,tv;
	long int porcnt,totcnt,sizecrit,nclusters=0;

        printf("Enter name of file to save microstructure to \n");
        scanf("%s",fileout);
        printf("%s\n",fileout);

	printf("Enter minimum size of clusters to keep \n");
        scanf("%ld",&sizecrit);
        printf("%ld\n",sizecrit);

	porcnt=totcnt=0;
	outfile=fopen(fileout,"w");
	
	/* Examine each voxel in the 3-D microstructure */
	for(i3=0;i3<SYSSIZE;i3++){
	for(i2=0;i2<SYSSIZE;i2++){
	for(i3=0;i3<SYSSIZE;i3++){
		totcnt+=1;
		iout=0;
                tv=particle [i1][i2][i3];
                /* Determine if voxel is part of a large enough cluster */
		if((tv!=0)&&(npixels[tv]>=sizecrit)){
			iout=1;
			porcnt+=1;
                        if(pkeep[tv]==0){
                           pkeep[tv]=1;
                           nclusters+=1;
                        }
		}
		fprintf(outfile,"%d\n",iout);
	}
	}
	}
	printf("Count for new solids is %ld out of %ld \n",porcnt,totcnt);
	printf("Count for clusters kept is %ld  \n",nclusters);
	fclose(outfile);
}

main()
{
	int i,nseed,menuch;

	/* Initialize arrays */	
        for(i=0;i<NUMPARTS;i++){
          npixels[i]=0;
          pkeep[i]=0;
        }

	printf("Enter random number seed: \n");
	scanf("%d",&nseed);
	printf("%d \n",nseed);
	srandom(nseed);

	menuch=7;
	while(menuch!=0){
		printf("Enter choice: \n");
		printf("0) Exit program \n");
		printf("1) Read in microstructure from file \n");
		printf("2) Measure phase fractions \n");
		printf("3) Check particle sizes \n");
		printf("4) Output particle microstructure to file \n");
		scanf("%d",&menuch);
		printf("%d \n",menuch);

		switch(menuch){
			case 1:
				readmic();
				break;
			case 2:
				phcount();
				break;
                        case 3:
				burncer();
				break;
			case 4:
				savemic();
				break;
			default: 
				break;
		}
	}

}