/************************************************************************/
/*                                                                      */
/*         Program: median.c                                            */
/*         Purpose: to perform a 3x3x3 or 5x5x5 median filter on a      */
/*                  3-D segmented image volume                          */
/*         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 XSIZE 300
#define YSIZE 300
#define ZSIZE 300

/* Global image array */
static short int image [XSIZE] [YSIZE] [ZSIZE];

short int findmed(ix,iy,iz,imedsize)
int ix,iy,iz,imedsize;
{
	short int temp,valback,values[126];
	int numtot,itrue,jtrue,id,jd,kd,idx,idy,idz;
	int numhalf,i,j;

	/* First store all of the values in the cubic box */
	numtot=0;
	for(id=-imedsize;id<=imedsize;id++){
		idx=ix+id;
		if((idx>=0)&&(idx<XSIZE)){
			itrue=1;
		}
		else{
			itrue=0;
		}
	for(jd=-imedsize;((jd<=imedsize)&&(itrue==1));jd++){
		idy=iy+jd;
		if((idy>=0)&&(idy<YSIZE)){
			jtrue=1;
		}
		else{
			jtrue=0;
		}
	for(kd=-imedsize;((kd<=imedsize)&&(itrue==1)&&(jtrue==1));kd++){
		idz=iz+kd;
		if((itrue)&&(jtrue)&&(idz>=0)&&(idz<ZSIZE)){
			values[numtot]=image[idx][idy][idz];
			numtot+=1;
		}
	}
	}
	}

	numhalf=numtot/2;
	/* Twin loop sort of values to find median */
	for(i=0;i<(numtot-1);i++){
	for(j=(i+1);j<numtot;j++){
		if(values[j]<values[i]){
			temp=values[j];
			values[j]=values[i];
			values[i]=temp;
		}
	}
	}
        /* Store median value and return it to the calling routine */
	valback=values[numhalf];
	return(valback);
}

main(){
	FILE *infile,*outfile;
	char filein[80],fileout[80];
	int i,j,k,medsize;
	int xsize,ysize,zsize;
	short int valin,valout;


	/* Get user input */
	printf("Enter name of file to read from \n");
	scanf("%s",filein);
	printf("%s\n",filein);
	infile=fopen(filein,"r");

	printf("Enter name of file to create \n");
	scanf("%s",fileout);
	printf("%s\n",fileout);
	outfile=fopen(fileout,"w");
        
        do{
	        printf("Enter size of median filter (3 or 5) \n");
       		scanf("%d",&medsize);
       		printf("%d\n",medsize);
        }while ((medsize!=3)&&(medsize!=5));

	if(medsize==3){
		medsize=1;
	}
	else{
		medsize=2;
	}

	/* Read in the initial microstructure */
	for(k=0;k<zsize;k++){
	for(j=0;j<ysize;j++){
	for(i=0;i<xsize;i++){
		fscanf(infile,"%d",&valin);
		image[i][j][k]=valin;
	}
	}
	}
	fclose(infile);

	/* Now in a second pass apply the median filter */
        /* and output the resultant image to user-specified file */
	for(k=0;k<zsize;k++){
	for(j=0;j<ysize;j++){
	for(i=0;i<xsize;i++){

		valout=findmed(i,j,k,medsize);
		fprintf(outfile,"%d\n",valout);
                fflush(outfile);
	}
	}
	}
	fclose(outfile);

}