#include<stdio.h>
#include<stdlib.h>
#include<math.h>

typedef unsigned char datatype;

void showUsage(char *);

int main(int argc, char **argv)
{
  FILE *fp;
  float sum,scale=1.0;
  int volNumber=0;
  int N,k,i,j,s,zs=124,xs=256,ys=256;
  int c;
  int verbose;
  datatype *tempImg;
  double *outImg;
  extern char *optarg;
  
  if(argc<3){
	  showUsage(argv[0]);
	  exit(EXIT_FAILURE);
  }

  while((c=getopt(argc-2,argv+2,"s:n:x:y:z:v"))!=-1)
    {
      switch(c)
	{
	case 's':
	  scale=atof(optarg);
	  break;
	case 'x':
	  xs=atoi(optarg);
	  break;
	case 'y':
	  ys=atoi(optarg);
	  break;
	case 'z':
	  zs=atoi(optarg);
	  break;
	case 'n':
	  volNumber = atoi(optarg);
	  break;
	case 'v':
	  verbose = 1;
	  break;
	case '?':
	  showUsage(argv[0]);
	  exit(EXIT_FAILURE) ;
	  break;
	}
    }
  if(verbose)
	  printf("Scale=%f,xs=%d, ys=%d, zs=%d, volNumber=%d\n",scale,xs,ys,zs,volNumber);
  outImg = (double *)malloc( xs * ys * zs * sizeof(double));
  tempImg = (datatype *)malloc( xs * ys * zs * sizeof(datatype));
  for(j = 0; j< xs*ys*zs;j++){
	  outImg[j] = 0;
	  tempImg[j] = 0;
  }
  
  fp=fopen(argv[1],"r");
  for(i = 0; i < volNumber; i++){
	  fread(tempImg,sizeof(datatype), xs*ys*zs,fp);
	  if(verbose)
		  printf("finished volume number%d\n",i);
	  for(j = 0; j< xs*ys*zs;j++)
		  outImg[j] += tempImg[j];
  } 
fclose(fp);
for(j = 0; j< xs*ys*zs;j++){
	outImg[j] /= volNumber;
	outImg[j] *= scale;
	tempImg[j] = outImg[j];
}
if(verbose)
	printf("Writing output\n");

fp=fopen(argv[2],"w");
fwrite(tempImg,sizeof(datatype),xs*ys*zs,fp);
close(fp);

return 0;
}

void showUsage(char *progname){
printf("USAGE: %s inputVolStack outputVol [options]\n\
\t  -x <int>               : height size,default 256 \n\
\t  -y <int>               : width size,default 256 \n\
\t  -z <int>               : z size,default 124 \n\
\t  -n <int>               : volumes number to average( part volumes also good,do not biger than the actual volume numbers )\n\
\t  -s <float>             : Scaling factor\n\
",progname);
exit(EXIT_FAILURE);
}
