/************************************************************************/
/*                                                                      */
/*         Program: percolate.c                                         */
/*         Purpose: to examine percolation of a phase in a 3-D volume   */
/*                  limited to cubic volumes (X * X * X)                */
/*         Input file: 3-D image, one integer per voxel                 */
/*         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 BURNT 70      /* label for a burnt pixel */
#define SIZE2D 90000 /* size of matrices for holding burning locations */
/* functions defining coordinates for burning in any of three directions */
#define cx(x,y,z,a,b,c) (1-b-c)*x+(1-a-c)*y+(1-a-b)*z
#define cy(x,y,z,a,b,c) (1-a-b)*x+(1-b-c)*y+(1-a-c)*z
#define cz(x,y,z,a,b,c) (1-a-c)*x+(1-a-b)*y+(1-b-c)*z

static short int mic[SYSSIZE][SYSSIZE][SYSSIZE];
long int ntot=0,nperc=0,nburnt=0;
FILE *resfile;

/* routine to assess the connectivity (percolation) of a single phase */
/* Two matrices are used here: one to store the recently burnt locations */
/*                the other to store the newly found burnt locations */
int burn3d(npix,d1,d2,d3)
        int npix;    /* ID of phase to perform burning on */
        int d1,d2,d3;   /* directional flags */
{
        long int ntop,nthrough,ncur,nnew,ntot,nphc;
        int i,inew,j,k,nmatx[SIZE2D],nmaty[SIZE2D],nmatz[SIZE2D];
        int xl,xh,j1,k1,px,py,pz,qx,qy,qz,xcn,ycn,zcn;
        int x1,y1,z1,igood,nnewx[SIZE2D],nnewy[SIZE2D],nnewz[SIZE2D];
        int jnew,icur;
	int bflag;

/* counters for number of pixels of phase accessible from surface #1 */
/* and number which are part of a percolated pathway to surface #2 */
        ntop=0;
	bflag=0;
        nthrough=0;
	nphc=0;

        /* percolation is assessed from top to bottom only */
        /* and burning algorithm is periodic in other two directions */
        /* use of directional flags allow transformation of coordinates */
        /* to burn in direction of choosing (x, y, or z) */
        i=0;

        for(k=0;k<SYSSIZE;k++){
        for(j=0;j<SYSSIZE;j++){

                igood=0;
                ncur=0;
                ntot=0;
                /* Transform coordinates */
                px=cx(i,j,k,d1,d2,d3);
                py=cy(i,j,k,d1,d2,d3);
                pz=cz(i,j,k,d1,d2,d3);
                if(mic [px] [py] [pz]==npix){
                        /* Start a burn front */
                        mic [px] [py] [pz]=BURNT;
                        ntot+=1;
                        ncur+=1;
                        /* burn front is stored in matrices nmat* */
                        /* and nnew* */
                        nmatx[ncur]=i;
                        nmaty[ncur]=j;
                        nmatz[ncur]=k;
                        /* Burn as long as new (fuel) pixels are found */
                        do{
                                nnew=0;
                               	for(inew=1;inew<=ncur;inew++){
                                        xcn=nmatx[inew];
                                       	ycn=nmaty[inew];
                                       	zcn=nmatz[inew];

                                        /* Check all six neighbors */
                                        for(jnew=1;jnew<=6;jnew++){
                                                x1=xcn;
                                                y1=ycn;
                                                z1=zcn;
                                                if(jnew==1){x1-=1;}
                                                if(jnew==2){x1+=1;}
                                                if(jnew==3){y1-=1;}
                                                if(jnew==4){y1+=1;}
                                                if(jnew==5){z1-=1;}
                                                if(jnew==6){z1+=1;}
						/* Periodic in y and */
                                                if(y1>=SYSSIZE){y1-=SYSSIZE;}
                                                else if(y1<0){y1+=SYSSIZE;}	
						/* Periodic in z direction */	
                                                if(z1>=SYSSIZE){z1-=SYSSIZE;}
                                                else if(z1<0){z1+=SYSSIZE;}

			/* Nonperiodic so be sure to remain in the 3-D box */
                                                if((x1>=0)&&(x1<SYSSIZE)){
                                                /* Transform coordinates */
                                                px=cx(x1,y1,z1,d1,d2,d3);
                                                py=cy(x1,y1,z1,d1,d2,d3);
                                                pz=cz(x1,y1,z1,d1,d2,d3);
                                                if(mic [px] [py] [pz]==npix){
                                                   ntot+=1;
                                                   mic [px] [py] [pz]=BURNT;
                                                   nnew+=1;
                                                   if(nnew>=SIZE2D){
                                           printf("error in size of nnew \n");
                                                   }
                                                   nnewx[nnew]=x1;
                                                   nnewy[nnew]=y1;
                                                   nnewz[nnew]=z1;
                                                }
                                        }
                                       	}
                                }
                                if(nnew>0){
                                        ncur=nnew;
                                        /* update the burn front matrices */
                                       	for(icur=1;icur<=ncur;icur++){
                                                nmatx[icur]=nnewx[icur];
                                                nmaty[icur]=nnewy[icur];
                                                nmatz[icur]=nnewz[icur];
                                        }
                                }
                        }while (nnew>0);

                        ntop+=ntot;
                        xl=0;
                        xh=SYSSIZE-1;
                    /* See if current path extends through the microstructure */
                        for(j1=0;j1<SYSSIZE;j1++){
                        for(k1=0;k1<SYSSIZE;k1++){
                                px=cx(xl,j1,k1,d1,d2,d3);
                                py=cy(xl,j1,k1,d1,d2,d3);
                                pz=cz(xl,j1,k1,d1,d2,d3);
                                qx=cx(xh,j1,k1,d1,d2,d3);
                                qy=cy(xh,j1,k1,d1,d2,d3);
                                qz=cz(xh,j1,k1,d1,d2,d3);
                   if((mic [px] [py] [pz]==BURNT)&&(mic [qx] [qy] [qz]==BURNT)){
                                        igood=2;
                                }
                                if(mic [px] [py] [pz]==BURNT){
                                        mic [px] [py] [pz]=BURNT+1;
                                }
                                if(mic [qx] [qy] [qz]==BURNT){
                                        mic [qx] [qy] [qz]=BURNT+1;
                                }
                        }
                        }

                        if(igood==2){
                                nthrough+=ntot;
                        }
               }
        }	
        }
/* return the burnt sites to their original phase values */
        for(k=0;k<SYSSIZE;k++){
        for(j=0;j<SYSSIZE;j++){
        for(i=0;i<SYSSIZE;i++){
                if(mic [i] [j] [k]>=BURNT){
			nphc+=1;
                        mic [i] [j] [k]=npix;
                }
		else if(mic[i][j][k]==npix){
			nphc+=1;
		}
        }
       	}
       	}

        /* Output results to screen and output file */
        printf("Phase ID= %d \n",npix);
        printf("Number accessible from first surface = %ld \n",ntop);
        printf("Number contained in through pathways= %ld \n",nthrough);
	ntot=nphc;
	nperc=nthrough;
	nburnt=ntop;
	if(nthrough>0){
		bflag=1;
	}
        fprintf(resfile,"Phase ID= %d \n",npix);
        fprintf(resfile,"Total number of pixels = %ld \n",ntot);
        fprintf(resfile,"Number accessible from first surface = %ld \n",ntop);
        fprintf(resfile,"Number contained in through pathways= %ld \n",nthrough);
	return(bflag);
}

main()
{
	int ix,iy,iz,valin;
	FILE *infile;
	char filein[80],fileout[80];
	int phasein,garb;

        printf("Enter name of input 3-D image file \n");
	scanf("%s",filein);
        printf("Enter phase ID to examine (0, 1, or 2, etc.) \n");
	scanf("%d",&phasein);
        printf("Enter name of output file to write percolation results to\n");
	scanf("%s",fileout);

	infile=fopen(filein,"r");

	/* Read in the 3-D image, voxel by voxel */
	for(iz=0;iz<SYSSIZE;iz++){
	for(iy=0;iy<SYSSIZE;iy++){
	for(ix=0;ix<SYSSIZE;ix++){
		fscanf(infile,"%d",&valin);
		mic[ix][iy][iz]=valin;
	}
	}
	}

	fclose(infile);
	resfile=fopen(fileout,"w");
	/* Assess percolation in each of the three principal directions */
	garb=burn3d(phasein,1,0,0);
	garb=burn3d(phasein,0,1,0);
	garb=burn3d(phasein,0,0,1);

	fclose(resfile);
}