#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<mvcd.h>
typedef unsigned char datatype;
void showUsage();
main(int argc, char **argv)
{
int x = 0, y = 0, z =0, k = 0, i = 0, j =0;
int xs = 256, ys = 256, zs = 124; /*input dimension*/
int oxs = 256, oys = 256, ozs = 124; /*output dimension*/
float xs_res = 0.9375, ys_res = 0.9375, zs_res = 1.5; /*input spacing */
float oxs_res = 0.9375, oys_res = 0.9375, ozs_res = 1.5; /*output spacing */
int xs_origin = 0, ys_origin = 0, zs_origin = 0; /* input origin */
int oxs_origin = 0, oys_origin = 0, ozs_origin = 0; /* output origin */
int xs_interpType = 0, ys_interpType = 0, zs_interpType = 0; /* interplation type */ 
char filename[256];
float distance[9]; /* distance to the eight neighbor */
int fileSize = 0;
FILE *fp;

int c; /* getopt to fetch options */
extern char *optarg;

int verboseFlag = 0;

if(argc < 3){
	showUsage();
	exit(1);
}

while((c=getopt(argc-2,argv+2,"d:D:r:R:o:O:t:v"))!=-1)
{
	switch(c)
	{
		case 'v':
			verboseFlag = 1;
			break;
		case 'd':
			sscanf(optarg,"%d,%d,%d",&xs,&ys,&zs);
			if(verboseFlag)
				printf("input image volume dimension (x,y,z)=(%d,%d,%d)\n",xs,ys,zs);
			break;
		case 'D':
			sscanf(optarg,"%d,%d,%d",&oxs,&oys,&ozs);
			if(verboseFlag)
				printf("output image volume dimension (x,y,z)=(%d,%d,%d)\n",oxs,oys,ozs);
			break;
		case 'r':
			sscanf(optarg,"%f,%f,%f",&xs_res,&ys_res,&zs_res);
			if(verboseFlag)
				printf("input image volume resolution (x,y,z)=(%f,%f,%f)\n",xs_res,ys_res,zs_res);
			break;
		case 'R':
			sscanf(optarg,"%f,%f,%f",&oxs_res,&oys_res,&ozs_res);
			if(verboseFlag)
				printf("output image volume resolution (x,y,z)=(%f,%f,%f)\n",oxs_res,oys_res,ozs_res);
			break;
		case 'o':
			sscanf(optarg,"%d,%d,%d",&xs_origin,&ys_origin,&zs_origin);
			if(verboseFlag)
				printf("input image volume origin (x,y,z)=(%d,%d,%d)\n",xs_origin,ys_origin,zs_origin);
			break;
		case 'O':
			sscanf(optarg,"%d,%d,%d",&oxs_origin,&oys_origin,&ozs_origin);
			if(verboseFlag)
				printf("output image volume origin (x,y,z)=(%d,%d,%d)\n",oxs_origin,oys_origin,ozs_origin);
			break;
		case 't':
			sscanf(optarg,"%d,%d,%d",&xs_interpType,&ys_interpType,&zs_interpType);
			if(verboseFlag)
				printf("interplation method is (x,y,z) = (%d,%d,%d)\n",xs_interpType,ys_interpType,zs_interpType);
			break;
		case '?':
			showUsage();
			exit(1);
			break;
	}
}

datatype ***orgImg;
datatype ***outImg;


sprintf(filename,argv[1]);
fp=myopen(filename,"r");
fseek(fp,0,SEEK_END);
fileSize = ftell(fp);
zs = fileSize / (xs * ys * sizeof(datatype));
rewind(fp);
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);
}

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

if(2 == sizeof(datatype)){
	outImg = (datatype ***)USalloc3d(oxs,oys,ozs);
}else if(1 == sizeof(datatype)){
	outImg = (datatype ***)UCalloc3d(oxs,oys,ozs);
}

for(z=0; z<ozs;z++) {
	for(x=0; x<oxs; x++ ) {
		for(y=0;y<oys;y++) {
			outImg[z][x][y] = 0;
		}
	}
}

if(verboseFlag)
	printf("memory allocated successfully\n");

// the coordinate that lies on the original physical space
float f_x_coord = 0.0;
float f_y_coord = 0.0;
float f_z_coord = 0.0;
// the lineared value at x,y,z direction
float yy0 = 0.0;
float yy1 = 0.0;
float zz0 = 0.0;
float zz1 = 0.0;

/* successfully mapped using the center for resample
*/
if(xs_interpType == ys_interpType == zs_interpType == 0){
	if(verboseFlag)printf("using nearest neighbor\n");
	for(z = -ozs/2; z < ozs/2; z++){
		for(x = -oxs/2; x < oxs/2; x++){
			for(y = -oys/2; y < oys/2; y++){
				f_x_coord = xs/2 + (x * oxs_res + oxs_origin - xs_origin)/xs_res;
				f_y_coord = ys/2 + (y * oys_res + oys_origin - ys_origin)/ys_res;
				f_z_coord = zs/2 + (z * ozs_res + ozs_origin - zs_origin)/zs_res;
				if(f_x_coord < 0 || f_y_coord < 0 || f_z_coord < 0 || f_x_coord >= xs || f_y_coord >= ys || f_z_coord >= zs){
					outImg[z + ozs/2][x + oxs/2][y + oys/2] = 0;
				}else{
					outImg[z + ozs/2][x + oxs/2][y + oys/2] = orgImg[(datatype)f_z_coord][(datatype)f_x_coord][(datatype)f_y_coord];
				}
			}
		}
	}
}else if(xs_interpType == ys_interpType == zs_interpType == 1){
	if(verboseFlag)printf("using tri-linear interpolation\n ");
	for(z = -ozs/2; z < ozs/2; z++){
		printf(".");
		for(x = -oxs/2; x < oxs/2; x++){
			for(y = -oys/2; y < oys/2; y++){
				f_x_coord = 0.999999*xs/2 + (x * oxs_res + oxs_origin - xs_origin)/xs_res;
				f_y_coord = 0.999999*ys/2 + (y * oys_res + oys_origin - ys_origin)/ys_res;
				f_z_coord = 0.999999*zs/2 + (z * ozs_res + ozs_origin - zs_origin)/zs_res;
				if((int)floorf(f_x_coord) < 0 || (int)floorf(f_y_coord) < 0 || (int)floorf(f_z_coord) < 0 || (int)ceilf(f_x_coord) >= xs || (int)ceilf(f_y_coord) >= ys || (int)ceilf(f_z_coord) >= zs){
					outImg[z + ozs/2][x + oxs/2][y + oys/2] = 0;
				}else{
					distance[1] = f_y_coord - floorf(f_y_coord);
					distance[2] = ceilf(f_y_coord) - f_y_coord;
					distance[3] = f_x_coord - floorf(f_x_coord);
					distance[4] = ceilf(f_x_coord) - f_x_coord;
					distance[5] = f_z_coord - floorf(f_z_coord);
					distance[6] = ceilf(f_z_coord) - f_z_coord;
					
					yy0 = (distance[2] * orgImg[(int)floorf(f_z_coord)][(int)floorf(f_x_coord)][(int)floorf(f_y_coord)]) / (distance[2] + distance[1]) + (distance[1] * orgImg[(int)floorf(f_z_coord)][(int)floorf(f_x_coord)][(int)ceilf(f_y_coord)]) / (distance[2] + distance[1]);
					yy1 = (distance[2] * orgImg[(int)floorf(f_z_coord)][(int)ceilf(f_x_coord)][(int)floorf(f_y_coord)]) / (distance[2] + distance[1]) + (distance[1] * orgImg[(int)floorf(f_z_coord)][(int)ceilf(f_x_coord)][(int)ceilf(f_y_coord)]) / (distance[2] + distance[1]);
					zz0 = distance[4] * yy0 / (distance[4] + distance[3]) + distance[3] * yy1 / (distance[4] + distance[3]);
					yy0 = (distance[2] * orgImg[(int)ceilf(f_z_coord)][(int)floorf(f_x_coord)][(int)floorf(f_y_coord)]) / (distance[2] + distance[1]) + (distance[1] * orgImg[(int)ceilf(f_z_coord)][(int)floorf(f_x_coord)][(int)ceilf(f_y_coord)]) / (distance[2] + distance[1]);
					yy1 = (distance[2] * orgImg[(int)ceilf(f_z_coord)][(int)ceilf(f_x_coord)][(int)floorf(f_y_coord)]) / (distance[2] + distance[1]) + (distance[1] * orgImg[(int)ceilf(f_z_coord)][(int)ceilf(f_x_coord)][(int)ceilf(f_y_coord)]) / (distance[2] + distance[1]);
					zz1 = distance[4] * yy0 / (distance[4] + distance[3]) + distance[3] * yy1 / (distance[4] + distance[3]);

					
					outImg[z + ozs/2][x + oxs/2][y + oys/2] = distance[6] * zz0 / (distance[6] + distance[5]) + distance[5] * zz1 / (distance[6] + distance[5]);
					
				}
			}
		}
	}
}else{
	printf("please use the same interpolation method on 3 dimension or this interpolation method has not been applicated");
	exit(EXIT_FAILURE);
}
if(verboseFlag)
	printf("output image writing!\n");

sprintf(filename,argv[2]);
fp=myopen(filename,"w");
for(k=0; k < ozs; k++){
	for(i = 0; i < oxs; i++){
		fwrite(outImg[k][i],sizeof(datatype),oys,fp);
	}
}
fclose(fp);

//the end of program
}

void showUsage(){
	printf("This program was used to resample an input image volume to other resolution or dimension or both, please be noted that all the operation was done in physical space in stead of pixel/grid space\n");
	printf("USAGE:  resample inputImg outputImg\n\
\t   -d <int>,<int>,<int>       : intput volume dimension (row,col,depth)(256,256,124) \n\
\t   -D <int>,<int>,<int>       : output volume dimension (row,col,depth)(256,256,124) \n\
\t   -r <float>,<float>,<float> : intput volume resolution (row,col,depth)(0.9375,0.9375,1.5)\n\
\t   -R <float>,<float>,<float> : output volume resolution (row,col,depth)(0.9375,0.9375,1.5) \n\
\t   -o <int>,<int>,<int>       : intput volume origin(row,col,depth)(0,0,0)\n\
\t   -O <int>,<int>,<int>       : output volume origin(row,col,depth)(0,0,0) \n\
\t   -t <int>,<int>,<int>       : interplation method of row col depth(0),0(nearest neiber,DEFAULT) 1(tri-linear) 2(8-pt sinc) 3(B-spline)(currently please use the same interpolation method on these 3 directions\n\
");
exit(-1);
}