/* This progarm was revised  by Deepthi.Koka@uphs.upenn.edu in Aug,2008 to fix a bug in the -l option of this code that was broken on olympus.  */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<mvcd.h>
#include<termio.h>

typedef unsigned char datatype;
#define JAKOBLABELTOTAL 104

//#define DEBUG

#define NOTHING 0
#define VOLSUM 1
#define HISTINFO 2
#define INTENSITYSUM 3
#define LABELVOLSUM 4

void showUsage();

main(int argc, char **argv)
{
int i,j,k,x,y,z,xs=256,ys=256,zs=124;
float xs_res = 0.9375,ys_res = 0.9375,zs_res = 1.5;
datatype ***orgImg;
unsigned char ***outImg;
char filename[256];
int op = 0, sum = 0,nonZero = 0;
int c;
int verbose = 0, batchProcessing = 0;
int labelIndex = 0;
double labelVolSum = 0.0;
// input image format
int byteNumber = 1;
int fileSize=0,hdrSize = 0;
int hdrFlag = 0,outputFlag = 0;
char outFilename[FILENAME_MAX];

/* for termal information */
int fd,termCols = 80;
struct winsize s_winsize;

float histUpperThresh = 100.0;
int lowerThresh = 0.0;
extern char *optarg;
FILE *fp;

/*int labelsToPick[JAKOBLABELTOTAL];
float labelVolToPick[JAKOBLABELTOTAL]; */
int labelsToPick;
float labelVolToPick = 0.0;
int labelsToPickCnt = 0;
char *pch;

if(argc<2){
	showUsage();
	exit(1);
}
/* Initialization */
/*for(i = 0; i < JAKOBLABELTOTAL; i++) {
	labelVolToPick[i] = 0.0;
}*/
while((c=getopt(argc-1,argv+1,"x:y:z:sSh:l:o:vbV:"))!=-1)
    {
      switch(c)
	{
	case 'x':
	  xs=atoi(optarg);
	  break;
	case 'y':
	  ys=atoi(optarg);
	  break;
	case 'z':
	  zs=atoi(optarg);
	  hdrFlag = 1;
	  break;
	case 'S':
	  op=VOLSUM;
	  break;
	case 's':
	  op=INTENSITYSUM;
	  break;
	case 'h':
	  op=HISTINFO;
	  sscanf(optarg,"%f,%d\n",&histUpperThresh,&lowerThresh);
	  break;
	case 'o':
	  outputFlag = 1;
	  sprintf(outFilename,optarg);
	  printf("%c",*optarg);
	  printf("output file is in byte mode with name %s\n",outFilename);
	  break;
	case 'l':
	  op=LABELVOLSUM;
	  labelIndex=atoi(optarg);
	  /*pch = strtok(optarg,",");
	  labelsToPick[0] = atoi(pch);
	  while(pch != NULL){
		  pch = strtok(NULL,",");
		  labelsToPickCnt++;
		  labelsToPick[labelsToPickCnt] = atoi(pch);
	        
	  }
	  if(verbose){
		  for(i = 0; i < labelsToPickCnt; i++)
			  printf("labels you want to pick is %d\n",labelsToPick[i]);
	  }*/
	  labelsToPick=labelIndex;
	  break;
	case 'V':
	  sscanf(optarg,"%f,%f,%f\n",&xs_res,&ys_res,&zs_res);
	  if(verbose)printf("resolution (%f,%f,%f)\n",zs_res,ys_res,xs_res);
	  break;
	case 'v':
	  verbose=1;
	  break;
#ifdef DEBUG
	  printf("calculate the sum of label volume\n");
#endif
	case 'b':
	  batchProcessing = 1;
	  break;
	case '?':
	  showUsage();
	  exit(1);
	  break;
	}
    }
byteNumber = sizeof(datatype);

sprintf(filename,argv[1]);
fp=myopen(filename,"r");
fseek(fp,0,SEEK_END);
fileSize = ftell(fp);
if(hdrFlag){
	hdrSize = fileSize - xs*ys*zs*sizeof(datatype);
}else{
	zs = fileSize / (xs * ys * sizeof(datatype));
	hdrSize = 0;
}

fseek(fp,hdrSize,SEEK_SET);
#ifdef DEBUG
printf("Input image para: xs=%d,ys=%d,zs=%d\n",xs,ys,zs);
#endif
/*
if(2 == sizeof(datatype)){
	orgImg = (datatype ***)USalloc3d(xs,ys,zs); 
}else if(1 == sizeof(datatype)){
	orgImg = (datatype ***)UCalloc3d(xs,ys,zs);
}else{
	printf("error\n");
	exit(-1);
}
*/

orgImg = (datatype ***)MyAlloc3d(sizeof(datatype), xs, ys, zs);
outImg = (unsigned char ***)MyAlloc3d(sizeof(unsigned char), xs, ys, zs);

for(z=0; z<zs;z++) {
    for(x=0; x<xs; x++ ) {
      for(y=0;y<ys;y++) {
        orgImg[z][x][y] = 0;
      }
    }
}
#ifdef DEBUG
	printf("memory allocated successfully\n");
#endif

for(k=0; k<zs; k++)
	for(i = 0; i<xs; i++){
		fread(orgImg[k][i],sizeof(datatype),ys,fp);
	}
fclose(fp);

if( op == VOLSUM ){
	nonZero = 0;
	for(z=0; z<zs;z++) {
		for(x=0; x<xs; x++ ) {
			for(y=0;y<ys;y++) {
				if(0 != orgImg[z][x][y])nonZero++;
			}
		}
	}
	//printf("sum of %d nonZero vols is %d\n",nonZero,sum);
	printf("%d\n",nonZero);
}
if( op == INTENSITYSUM ){
	nonZero = 0;
	for(z=0; z<zs;z++) {
		for(x=0; x<xs; x++ ) {
			for(y=0;y<ys;y++) {
				if(0 != orgImg[z][x][y])
					nonZero +=orgImg[z][x][y];
			}
		}
	}
	//printf("sum of %d nonZero vols is %d\n",nonZero,sum);
	if(verbose){
		printf("sum of nonZero Vols intensity is %d\n",nonZero);
	}else{
		printf("%d\n",nonZero);
	}
}
if( op == LABELVOLSUM){
	for(z=0; z<zs;z++) {
		for(x=0; x<xs; x++ ) {
			for(y=0;y<ys;y++) {
				if(orgImg[z][x][y]){
					/*for(i = 0; i < labelsToPickCnt; i++){*/
						if(orgImg[z][x][y] == labelsToPick) {
							labelVolToPick += xs_res * ys_res * zs_res;
						}
					/*}*/
				}
			}
		}
	}
#ifdef DEBUG
	printf("sum of label %d is %1.1f\n",labelIndex,labelVolSum);
#endif
	/*for(i = 0; i < labelsToPickCnt; i++){
		printf("%1.1f ",labelVolToPick[i]);
		printf("\n");
	}*/
	if(verbose){
		printf("sum of label - %d  is : %1.1f\n",labelsToPick,labelVolToPick);
	}else{
		printf("%1.1f\n",labelVolToPick);
	}
	
}

if(op == HISTINFO){
	int histLen = (1<<8*sizeof(datatype));
	int hist[histLen];
	int histAcc = 0;
	int histUpperThreshold_val = 0;
	int j = 0;
	int signedShortFlag = 0;
	
	for(i = 0; i < histLen; i++){
		hist[i] = 0;
	}
	if(verbose)printf("byteNumber=%d,histLen=%d\n",byteNumber,histLen);
	for(z=0; z<zs;z++) {
		for(x=0; x<xs; x++ ) {
			for(y=0;y<ys;y++) {
				if(orgImg[z][x][y] < 0){
					signedShortFlag += 1;
				}
			}
		}
	}
	if(verbose)printf("signed short input points below zero is %d\n",signedShortFlag);
	/*
	if(signedShortFlag){
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0;y<ys;y++) {
					hist[orgImg[z][x][y]+32767]++;
					nonZero++;
				}
			}
		}
	}else{*/
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0;y<ys;y++) {
					if( orgImg[z][x][y] > lowerThresh){
						hist[orgImg[z][x][y]]++;
						nonZero++;
					}
				}
			}
		}
/*	}*/
	for(i = 1; i < histLen; i++){
//		printf("HERE hist[%d]=%d,histAcc=%d,nonZero=%d \n",i,hist[i],histAcc,nonZero);
		histAcc += hist[i];
		if(histAcc >= (nonZero * 1.0 * histUpperThresh/100.0)){
			histUpperThreshold_val = i;
		       	i += histLen;
		}
	}
	if(verbose)printf("the value at %f%% of the hist is %d\n",histUpperThresh,histUpperThreshold_val); 
	/* output the byte image */
	if(outputFlag){
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0;y<ys;y++) {
					if( orgImg[z][x][y] >= lowerThresh && orgImg[z][x][y] <= histUpperThreshold_val*1.2){ /* give some space betwennthe scale */
						outImg[z][x][y] = (orgImg[z][x][y] * 255.0)/(histUpperThreshold_val*1.2);
					}else if (orgImg[z][x][y]<lowerThresh){
						outImg[z][x][y] = 0;
					}else{
						outImg[z][x][y] = 255;
					}
				}
			}
		}
		fp = myopen(outFilename,"w");
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				fwrite(outImg[z][x],sizeof(unsigned char),ys,fp);
			}
		}
		fclose(fp);
	}
	/* printout the histgram */
	ioctl(fd,TIOCGWINSZ, &s_winsize);
	termCols = s_winsize.ws_col; /* get terminal colrs information */
	int histPeakNum = 0;
	for(i = 1; i < histLen; i++){
		if(histPeakNum < hist[i])
			histPeakNum = hist[i];
	}
	int histDisplayCols = termCols - 15;
	if(0 == batchProcessing)printf("for the below histogram, each point represents %d points\n",histPeakNum/histDisplayCols);	
	for(i = 1; i < histUpperThreshold_val; i++){
		if(hist[i] > histPeakNum/(histDisplayCols)){
			if(0 == batchProcessing)printf("%5d\033[33m",i);
			if(0 == batchProcessing)
				for(j = 0; j < hist[i]; j += histPeakNum/histDisplayCols){
					printf(".");
				}
			if(0 == batchProcessing)printf("\033[0m%d\n",hist[i]);
		}
	}
	/* using gnuplot for ploting */
	FILE *fp1 = fopen("hist.dat","w+");
	if(NULL == fp1)printf("ERROR open file hist.dat\n");
	fprintf(fp1,"%d %d\n",0,0);
	for(i = 1; i < histUpperThreshold_val; i++){
		if(hist[i] > histPeakNum/100)
			fprintf(fp1,"%d %d\n",i,hist[i]);
	}
	fclose(fp1);
	if(0 == batchProcessing){
		FILE *pPlot = popen("/usr/freeware/bin/gnuplot>dump1","w");
		fprintf(pPlot, "%s\n","set terminal x11;set title 'histgrams of the input image volume'\n");
		fflush(pPlot);
		fprintf(pPlot, "%s\n","set xlabel 'Intensity';set ylabel 'Points'\n");
		fflush(pPlot);
		fprintf(pPlot, "%s\n","set autoscale y\n");
		fflush(pPlot);
		//fprintf(pPlot,"%s","plot 'hist.dat' with lines\\\n");
		fprintf(pPlot,"%s\n","plot 'hist.dat' with impulses,'hist.dat' with lines");
		fflush(pPlot);
		printf("press any key to quit.....\n");
		getchar();
		system("rm hist.dat");
		pclose(pPlot);
		/* end gnuplot */
	}
}
//end of program
}

void showUsage(){
   printf("USAGE:\033[40;31m imageinfo inputImg [options]\n\
\t  \033[31;32m -x  <int>               : x size \n\
\t  \033[31;34m -y  <int>               : y size \n\
\t  \033[31;35m -z  <int>               : z size \n\
\t  \033[31;36m -V <float,float,float>  : Voxel spacing in mm(x_res,y_res,z_res)\n\
\t  \033[31;33m -o <filename> 	        : must be used with -h to rescale two byte to one byte account perc of the hist  \n\
\t  \033[31;36m -S  	                : prints the sum of vol interms of # of voxels \n\
\t  \033[31;36m -s  	                : prints the sum of vol intensity \n\
\t  \033[31;37m -l  <int>               : prints the sum of certain label(value=<int>) in mm^3 \n\
\t  \033[31;32m -h  <float>,<int>       : prints the perc(99.00 default) of the hist \n\
\t  \033[31;33m -b                     : batch processing, no visual histograms\n\
\033[0m");
exit(1);
}
