/* 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 <unistd.h>
#include<math.h>
#include<mvcd.h>
#include<string.h>

typedef unsigned char datatype;
#define PI 3.14159
#define MAXGRAY 256
#define JAKOBLABELTOTAL 104
#define SUBTRACT -1
#define NOTHING 0
#define ADD 1
#define REPLACE 2
#define MASK 3
#define PICKLABEL 4
#define AND 5
#define EMPTY 6
#define FILLHOLE 7
#define CHGSLICENUM 8
#define INVERSE 9
#define MASKVAL 10
#define THRESH 11
#define SHIFTING 12
#define PADDING 13
#define SCALING 14
#define CUTTING 15
#define PICKING 16
#define MASKHOLE 17
#define SHIFTVAL 18
#define MERGE 19
#define ORIENTATION 20
#define FLIP 21
#define DELETELABEL 22
#define MASKOUT 23
#define EQUALIZATION 24
#define BIASCORRECTION 25
#define DIVIDE 26
#define MASKAVERAGE 27
#define SMOOTHING 28
#define MAXTOZERO 29

void endianSwap(char dataFormat, unsigned char * rawStream);
void endianSwap(char dataFormat, unsigned char * rawStream)
{
	unsigned char t;
	switch(dataFormat){
		case 'b':
			break;
		case 'u':
			t = *rawStream;
			*rawStream = rawStream[1];
			rawStream[1] = t;
			break;
		case 's':
			t = *rawStream;
			*rawStream = rawStream[1];
			rawStream[1] = t;
			break;
		case 'f':
			t = *rawStream;
			*rawStream = rawStream[3];
			rawStream[3] = t;
			t = rawStream[1];
			rawStream[1] = rawStream[2];
			rawStream[2] = t;
			break;
		default:
			break;
	}
}

void fillFilterArray(float *kernel1DArray,float res, float sigma,int filterWidth)
{
	int i,j,k;
	float sum = 0.0;
	for(i = 0; i < filterWidth; i++){
		kernel1DArray[i] = exp(-(i-(filterWidth-1)/2)*(i-(filterWidth-1)/2)/(2*sigma*sigma))/(sigma*sqrt(2*PI));
	}
	/* normalise the kernel array */
	for(i = 0; i < filterWidth; i++){
		sum += kernel1DArray[i];
	}
	for(i = 0; i < filterWidth; i++){
		kernel1DArray[i] /= sum;
	}
	printf("\n");

}
void fillHole(datatype ***, int , datatype ***, int, int, int);
void showUsage(char **argv);
void showDebugUsage(char **argv);

int main(int argc, char **argv)
{
datatype c_maxValue = 0;	
int i,j,k,x,y,z,xs=256,ys=256,zs=124,oxs=256,oys=256,ozs=124;
float x_res = 0.9375, y_res = 0.9375, z_res = 1.5;
int op = 0,maskVal,fillHoleSize = 1, sliceChgNum = 0;
datatype ***orgImg,***subImg,***outImg;
float *floatOutput;
char filename[256];
unsigned char shiftDirection;
int shiftAmount = 0,paddingAmountX = 0, paddingAmountY = 0;
int cuttingAmountX = 0, cuttingAmountY = 0,cuttingOffsetX = 0, cuttingOffsetY = 0;
int c;
int pick_low = 0,pick_high = zs, pick_jump = 1;
int shiftVal = 0;
int lowThreshold = 0, highThreshold = c_maxValue, threshedVal[3] = {0,-1,0};
char orientation[6];
int orientParaNum = 0;
int flipDirection = -1;
int labelVal=-1;
int verboseFlag = 0;
int hdrFlag = 0;
int fileSize = 0;
int maskout = 0; /* if maskout = -1 ,pick the non mask part */
int endianType = 1;
char dataFormat = 'b';
int valsToMask=-1;
float subtractScale = 1.0;
float divideScale = 1.0;
float scales = 1.0;

int labelsToPick;
int labelsToPickCnt = 0;
char *pch;
float fTemp1,fTemp2;
float kernelSize = 0.0;
long ucnt=0,lcnt=0;

extern char *optarg;
FILE *fp, *fp1;

if(argc<4){
	showUsage(argv);
	exit(1);
}

while((c=getopt(argc-3,argv+3,"x:y:z:Z:vV:m:M:f:F:c:C:T:s:S:p:P:w:aAbrR:l:ne:E:ihg:O:d:D:Q10"))!=-1)
    {
      switch(c)
	{
	case 'x':
	  xs=atoi(optarg);
	  break;
	case 'y':
	  ys=atoi(optarg);
	  break;
	case 'z':
	  zs=atoi(optarg);
	  hdrFlag = 1;
	  break;
	case 'Z':
	  op = SCALING;
	  scales = atof(optarg);
	  if(scales < 0.0){
		  printf("please input zoom factor: ");
		  scanf("%f",&scales);
	  }
	  scales = 75.0/scales;
	  break;
	case 'V':
	  sscanf(optarg,"%f,%f,%f",&x_res,&y_res,&z_res);
	  printf("resolutino of x,y,z is (%f,%f,%f)\n",x_res,y_res,z_res);
	  break;
	case 'v':
	  verboseFlag = 1;
	  break;
	case 'a':
	  op = ADD;
	  break;
	case 'A':
	  op = MERGE;
	  break;
	case '0':
	  op = MAXTOZERO;
	  break;
	case 's':
	  op = SUBTRACT;
	  subtractScale = atof(optarg);
	  break;
	case 'd':
	  op = DIVIDE;
	  divideScale = atof(optarg);
	  if(divideScale < 0.0){
		  printf("please input the delta age: ");
		  scanf("%f",&divideScale);
	  }
	  if(verboseFlag)printf("divideScale=%f\n",divideScale);
	  break;
	case 'Q':
	  op = EQUALIZATION;
	  break;
	case 'S':
	  sscanf(optarg,"%c,%d",&shiftDirection,&shiftAmount);
	  printf("%c direction shifting amount is %d\n",shiftDirection,shiftAmount);
	  op = SHIFTING;
	  break;
	case 'm':
	  op = MASK;
	  maskout = atoi(optarg);
	  break;
	case 'w':
	  op = SHIFTVAL;
	  shiftVal = atoi(optarg);
	  hdrFlag=1;
	  break;
	case 'T':
	  op = THRESH;
	  sscanf(optarg,"%d,%d,%d,%d,%d",&lowThreshold,&highThreshold,&threshedVal[0],&threshedVal[1],&threshedVal[2]);
	  break;
	case 'M':
	  op = MASKVAL;
	  sscanf(optarg,"%d,%d",&valsToMask,&maskVal);
	  break;
	case 'f':
	  op = FILLHOLE;
	  fillHoleSize = atoi(optarg);
	  break;
	case 'C':
	  op = CHGSLICENUM;
	  sliceChgNum = atoi(optarg);
	  printf("change slice number to %d\n",sliceChgNum);
	  break;
	case 'P':
	  op = PICKING;
	  sscanf(optarg,"%d,%d,%d",&pick_low,&pick_high,&pick_jump);
	  sliceChgNum = (int)ceilf(((float)(pick_high - pick_low +1)/(float)pick_jump));
	  break;
	case 'r':
	  op = REPLACE;
	  break;
	case 'R':
	  op = MASKAVERAGE;
	  maskVal = atoi(optarg);
	  break;
	case 'O':
	  op = ORIENTATION;
	  orientParaNum=sscanf(optarg,"%s",orientation);
	  if(verboseFlag)printf("orientParaNum is %d\n",orientParaNum);
	  break;
	case 'g':
	  op = SMOOTHING;
	  kernelSize = atof(optarg);
	  if(verboseFlag)printf("gaussian kernel size if %f mm",kernelSize);
	  break;
	case 'i':
	  op = INVERSE;
	  break;
	case 'l':
	  op = PICKLABEL;
	  labelVal = atoi(optarg);
	  printf("labelVal is : %d ",labelVal);
	  /*pch = strtok(optarg,",");
	  labelsToPick[0] = atoi(pch);
	  //pch = strtok(NULL,",");
	  //printf("H %s H\n",pch);fflush(stdout);
	  while(pch != NULL){
		  pch = strtok(NULL,",");
		  labelsToPickCnt++;
		  labelsToPick[labelsToPickCnt] = atoi(pch);
	  printf("JUNK %d,",labelsToPickCnt);fflush(stdout);
	  pch = strtok(NULL,",");
	  }
	  printf("JUNK");fflush(stdout);
	  if(verboseFlag){
		  for(i = 0; i < labelsToPickCnt; i++)
			  printf("labels you want to pick is %d\n",labelsToPick[i]);fflush(stdout);
	  } */
	  labelsToPick=labelVal;
	  if(verboseFlag) printf("label you want to pick is %d\n",labelsToPick);
	  break;
	case 'e':
	  op = EMPTY;
	  sliceChgNum = atoi(optarg);
	  break;
	case 'E':
	  endianType = atoi(optarg);
	  break;
	case 'n':
	  op = AND;
	  break;
	case 'h':
	  op = MASKHOLE;
	  break;
	case 'p':
	  op = PADDING;
	  sscanf(optarg,"%d,%d",&paddingAmountX,&paddingAmountY);
	  break;
	case 'c':
	  op = CUTTING;
	  sscanf(optarg,"%d,%d,%d,%d",&cuttingAmountX,&cuttingAmountY,&cuttingOffsetX,&cuttingOffsetY);
	  break;
	case 'F':
	  op = FLIP;
	  flipDirection = atoi(optarg);
	  break;
	case 'D':
	  op = DELETELABEL;
	  labelVal = atoi(optarg);
	  break;
	case 'b':
	  op = BIASCORRECTION;
	  printf("Thsi is bias correction using fsl generated bias field, the second input will be used as the bias field in float mode\n");
	  break;
	case '?':
	  showUsage(argv);
	  exit(1);
	  break;
	}
    }
if(1 == sizeof(datatype)){
	c_maxValue = 255;
	dataFormat = 'b';
}else if(4 == sizeof(datatype)){
	c_maxValue = (datatype)0xFFFFFFFF;
	dataFormat = 'f';
}else if( (32768 == (1<< (8*sizeof(datatype) - 1))) && (datatype)((1<< (8 * sizeof(datatype) - 1)) | 0x40000000 ) > 0){
	c_maxValue = (datatype)65535;
	dataFormat = 'u';
}else{
	c_maxValue = (datatype)32767;
	dataFormat = 's';
}
if(verboseFlag){
	printf("c_maxValue = %d,dataFormat is %c\n",c_maxValue,dataFormat);
	fflush(stdout);
}
sprintf(filename,argv[1]);
fp=myopen(filename,"r");
fseek(fp,(long)0,SEEK_END);
fileSize = ftell(fp);
if(!hdrFlag){
	zs = fileSize/(xs * ys * sizeof(datatype));
	rewind(fp);
}
if(verboseFlag)printf("filesize is %d, zs=%d,calced zs = %d",fileSize,zs,fileSize/(xs * ys * sizeof(datatype)));

if( 0 != fileSize%(xs * ys * zs *sizeof(datatype))){
	printf("Input Image has header in it with size of %d\n",fileSize-(zs * xs * ys * sizeof(datatype)));
	fseek(fp,fileSize-(zs * xs * ys * sizeof(datatype)),0);
}

if(verboseFlag){
	printf("Input image para: xs=%d,ys=%d,zs=%d\n",xs,ys,zs);
	printf("operation is %d\n",op);
}
/*
if(2 == sizeof(datatype)){
	orgImg = (datatype ***)USalloc3d(xs,ys,zs);
	subImg = (datatype ***)USalloc3d(xs,ys,zs);
}else{
	orgImg = (datatype ***)UCalloc3d(xs,ys,zs); 
	subImg = (datatype ***)UCalloc3d(xs,ys,zs); 
}
*/
orgImg = (datatype ***)MyAlloc3d(sizeof(datatype), xs, ys, zs);
subImg = (datatype ***)MyAlloc3d(sizeof(datatype), xs, ys, zs);
if(DIVIDE == op || SUBTRACT == op){
	floatOutput = malloc(xs*ys*zs*sizeof(float)); 
	for(i = 0; i < xs*ys*zs; i++){
		floatOutput[i] = 0.0;
	}
}
if(NULL == orgImg){
	printf("memory allocation casted 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;
        subImg[z][x][y] = 0;
      }
    }
}
if(verboseFlag){
	printf("memory allocated successfully\n");
}

unsigned char *dataFormatEndianTypeConversionArray;
datatype *dataFormatEndianTypeConversionArray_casted;
sprintf(filename,argv[2]);
fp1=myopen(filename,"r");
if( !endianType ){
	if(verboseFlag)
		printf("small endian!\n");
	dataFormatEndianTypeConversionArray = calloc(xs * ys *zs * sizeof(datatype), sizeof(unsigned char));
	if( NULL == dataFormatEndianTypeConversionArray){
		printf(" Error in dataFormatEndianTypeConversionArray memory allocation\n");
		exit(-1);
	}
	if( xs*ys*zs*sizeof(datatype) != fread(dataFormatEndianTypeConversionArray, sizeof(unsigned char),xs*ys*zs*sizeof(datatype),fp)){
		printf(" Error in dataFormatEndianTypeConversionArray data reading\n");
		exit(-1);
	}
	for(j = 0; j < xs * ys * zs * sizeof(datatype); j += sizeof(datatype)){
		endianSwap(dataFormat,&dataFormatEndianTypeConversionArray[j]);
	}
	dataFormatEndianTypeConversionArray_casted = (datatype *)dataFormatEndianTypeConversionArray;
	for(z=0; z<zs;z++) {
		for(x=0; x<xs; x++ ) {
			for(y=0;y<ys;y++) {
				orgImg[z][x][y] = dataFormatEndianTypeConversionArray_casted[z*xs*ys + x*ys + y];
			}
		}
	}
	/* read in the second input image volume */
	if( xs*ys*zs*sizeof(datatype) != fread(dataFormatEndianTypeConversionArray, sizeof(unsigned char),xs*ys*zs*sizeof(datatype),fp1)){
		printf(" Error in dataFormatEndianTypeConversionArray data reading\n");
		exit(-1);
	}
	for(j = 0; j < xs * ys * zs * sizeof(datatype); j += sizeof(datatype)){
		endianSwap(dataFormat,&dataFormatEndianTypeConversionArray[j]);
	}
	dataFormatEndianTypeConversionArray_casted = (datatype *)dataFormatEndianTypeConversionArray;
	for(z=0; z<zs;z++) {
		for(x=0; x<xs; x++ ) {
			for(y=0;y<ys;y++) {
				subImg[z][x][y] = dataFormatEndianTypeConversionArray_casted[z*xs*ys + x*ys + y];
			}
		}
	}
	free(dataFormatEndianTypeConversionArray);
}else{
	for(k=0; k<zs; k++)
		for(i = 0; i<xs; i++){
			fread(orgImg[k][i],sizeof(datatype),ys,fp);
		}
	for(k=0; k<zs; k++)
		for(i = 0; i<xs; i++){
			fread(subImg[k][i],sizeof(datatype),ys,fp1);
		}
}


fclose(fp);
fclose(fp1);

oxs = xs;
oys = ys;
ozs = (sliceChgNum>0?sliceChgNum:zs);

/* when padding required */
if(PADDING == op){
	printf(" padding image volume x,y for %d,%d more\n",paddingAmountX,paddingAmountY);
	oys = ys + paddingAmountY;
	oxs = xs + paddingAmountX;
}
if(CUTTING == op){
	printf(" cutting image volume x,y for (%d,%d) less,with offset of (%d,%d)\n",cuttingAmountX,cuttingAmountY,cuttingOffsetX,cuttingOffsetY);
	oys = ys - cuttingAmountY;
	oxs = xs - cuttingAmountX;
}

/* output memory allocate */
if(verboseFlag){
printf("output image para oxs=%d,oys=%d, ozs=%d\n",oxs,oys,ozs);
}

if( ORIENTATION == op){
	int temp[3];
	temp[0] = zs; 
	temp[1] = xs; 
	temp[2] = ys;
	i = 0;
	int temp_x,temp_y,temp_z;
	char tempxyz[3];
	int xyzDirection[3] ;
	xyzDirection[0] = 0;
	xyzDirection[1] = 0;
	xyzDirection[2] = 0;
	for( k = 0; k < 3; k++){
		if( '-' == orientation[i] ){
			switch(orientation[i+1]){
				case 'x':
					xs = temp[k];
					tempxyz[k] = 'x';
					xyzDirection[k] = 1;
					break;
				case 'y':
					ys = temp[k];
					tempxyz[k] = 'y';
					xyzDirection[k] = 1;
					break;
				case 'z':
					zs = temp[k];
					tempxyz[k] = 'z';
					xyzDirection[k] = 1;
					break;
				default:
					printf("error passing orientation parameters\n");
					exit(EXIT_FAILURE);
					break;
			}
			i++;
		}else{
			switch(orientation[i]){
				case 'x':
					xs = temp[k];
					tempxyz[k] = 'x';
					break;
				case 'y':
					ys = temp[k];
					tempxyz[k] = 'y';
					break;
				case 'z':
					zs = temp[k];
					tempxyz[k] = 'z';
					break;
				default:
					printf("error passing orientation parameters\n");
					exit(EXIT_FAILURE);
					break;
			}
		}
		
	i++;
	}
	oxs = xs;
	oys = ys;
	ozs = zs;
	printf("tempxyz is from (z,x,y) to (%c,%c,%c)\n",tempxyz[0],tempxyz[1],tempxyz[2]);
	printf("rotated result of z,x,y are(%d,%d,%d)\n",zs,xs,ys);
	outImg = (datatype ***)MyAlloc3d(sizeof(datatype), xs, ys, zs);
	if(NULL == outImg){
		printf("memory allocation casted error\n");
		exit(-1);
	}
	if(verboseFlag)printf("outImge mem allocated\n");
	for(z = 0; z < zs; z++) {
		for(x=0; x < xs; x++ ) {
			for(y=0; y < ys; y++) {
				if('x' == tempxyz[0]){
					temp_z = xyzDirection[0] > 0?(xs -x-1):x;
				}else if('y' == tempxyz[0]){
					temp_z = xyzDirection[0] > 0?(ys -y-1):y;
				}else if('z' == tempxyz[0]){
					temp_z = xyzDirection[0] > 0?(zs -z-1):z;
				}else {
					printf("temp_x error\n");
					exit(-1);
				}

				if('x' == tempxyz[1]){
					temp_x = xyzDirection[1] > 0?(xs -x-1):x;
				}else if('y' == tempxyz[1]){
					temp_x = xyzDirection[1] > 0?(ys -y-1):y;
				}else if('z' == tempxyz[1]){
					temp_x = xyzDirection[1] > 0?(zs -z-1):z;
				}else {
					printf("temp_y error\n");
					exit(-1);
				}
				if('x' == tempxyz[2]){
					temp_y = xyzDirection[2] > 0?(xs -x-1):x;
				}else if('y' == tempxyz[2]){
					temp_y = xyzDirection[2] > 0?(ys -y-1):y;
				}else if('z' == tempxyz[2]){
					temp_y = xyzDirection[2] > 0?(zs -z-1):z;
				}else {
					printf("temp_z error\n");
					exit(-1);
				}
				 outImg[z][x][y] = (datatype)orgImg[temp_z][temp_x][temp_y];
			}
		}
	}
}else {
	outImg = (datatype ***)MyAlloc3d(sizeof(datatype), oxs, oys, (sliceChgNum > 0?sliceChgNum:ozs));
	if(NULL == outImg){
		printf("memory allocation casted error\n");
		exit(-1);
	}
	/*
	if(2 == sizeof(datatype)){
		outImg = (datatype ***)USalloc3d(oxs,oys,(sliceChgNum > 0?sliceChgNum:zs));
	}else{
		outImg = (datatype ***)UCalloc3d(oxs,oys,(sliceChgNum > 0?sliceChgNum:zs)); 
	}
	*/
	for(z=0; z<(sliceChgNum > 0?sliceChgNum:ozs);z++) {
		for(x=0; x<oxs; x++ ) {
			for(y=0;y<oys;y++) {
				outImg[z][x][y] = 0;
			}
		}
	}
}
/*
for(z=0; z<(sliceChgNum>zs?zs:sliceChgNum);z++) {
  for(x=0; x<xs; x++ ) {
     for(y=0; y<ys; y++) {
          outImg[z][x][y] = orgImg[z][x][y];
     }
  }
}
*/

/* do nothing */
for(z=0;z< (sliceChgNum >= zs?zs:sliceChgNum);z++){
	for(x=0;x<oxs;x++){
		for(y=0;y<oys;y++){
			outImg[z][x][y] = orgImg[z][x][y];
		}
	}
}

if( FILLHOLE == op){
	for(z=0;z<zs;z++)
		for(x=0;x<xs;x++)
			for(y=0;y<ys;y++)
				outImg[z][x][y] = orgImg[z][x][y];
       fillHole(orgImg, fillHoleSize, outImg, xs, ys, zs);
}
/* test area */

/* test area */

if( SMOOTHING == op ){
	float *kernel1DArray;
	int kernelWidth;
	float sigma;
	float sum;
	sigma = kernelSize / sqrt(log(4));
	if(verboseFlag)printf("sigma isze is %f\n",sigma);
	kernelWidth = 2*kernelSize/y_res+1;
	kernel1DArray = calloc(kernelWidth,sizeof(float));
	fillFilterArray(kernel1DArray,y_res,sigma,2*kernelSize/y_res+1);
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				if(subImg[z][x][y] < 0){
					orgImg[z][x][y] = 0;
				}
			}
		}
	}
	if(verboseFlag){
		printf("the filter array is:\n ");
		for(i = 0; i < kernelWidth; i++){
			printf("%f ",kernel1DArray[i]);
		}
	}
	printf("\nsmoothing y direction with kwidth=%d ",kernelWidth);
	for(z=0;z<zs;z++){if(z%10 == 0){printf(".");fflush(stdout);}
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				outImg[z][x][y] = 0;
				for(i = -(kernelWidth-1)/2; i <= (kernelWidth-1)/2; i++){
					if( y+i < 0 || y+i >= ys )
						;
					else
						outImg[z][x][y] += orgImg[z][x][y+i] * kernel1DArray[i+(kernelWidth-1)/2];
				}
			}
		}
	}
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				orgImg[z][x][y] = outImg[z][x][y];
			}
		}
	}
	free(kernel1DArray);
	kernelWidth = 2*kernelSize/x_res+1;
	kernel1DArray = calloc(kernelWidth,sizeof(float));
	fillFilterArray(kernel1DArray,x_res,sigma,2*kernelSize/x_res+1);
	printf("\nsmoothing x direction with kwidth=%d",kernelWidth);
	if(verboseFlag){
		printf("the filter array is:\n ");
		for(i = 0; i < kernelWidth; i++){
			printf("%f ",kernel1DArray[i]);
		}
	}
	for(z=0;z<zs;z++){if(z%10 == 0){printf(".");fflush(stdout);}
		for(y=0;y<ys;y++){
			for(x=0;x<xs;x++){
				outImg[z][x][y] = 0;
				for(i = -(kernelWidth-1)/2; i <= (kernelWidth-1)/2; i++){
					if( x+i < 0 || x+i >= xs )break;
					else
						outImg[z][x][y] += orgImg[z][x+i][y] * kernel1DArray[i+(kernelWidth-1)/2];
				}
			}
		}
	}
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				orgImg[z][x][y] = outImg[z][x][y];
			}
		}
	}
	free(kernel1DArray);
	kernelWidth = 2*kernelSize/z_res+1;
	kernel1DArray = calloc(kernelWidth,sizeof(float));
	fillFilterArray(kernel1DArray,z_res,sigma,2*kernelSize/z_res+1);
	printf("\nsmoothing z direction with kwidth=%d",kernelWidth);
	if(verboseFlag){
		printf("the filter array is:\n ");
		for(i = 0; i < kernelWidth; i++){
			printf("%f ",kernel1DArray[i]);
		}
	}
	for(x=0;x<xs;x++){if(x%10 == 0){printf(".");fflush(stdout);}
		for(y=0;y<ys;y++){
			for(z=0;z<zs;z++){
				outImg[z][x][y] = 0;
				for(i = -(kernelWidth-1)/2; i <= (kernelWidth-1)/2; i++){
					if( z+i < 0 || z+i >= zs )break;
					else
						outImg[z][x][y] += orgImg[z+i][x][y] * kernel1DArray[i+(kernelWidth-1)/2];
				}
			}
		}
	}
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				orgImg[z][x][y] = outImg[z][x][y];
			}
		}
	}
}
if( CHGSLICENUM == op){
	for(z=0;z<(sliceChgNum<zs?sliceChgNum:zs);z++)
		for(x=0;x<xs;x++)
			for(y=0;y<ys;y++)
				outImg[z][x][y] = orgImg[z][x][y];
}

if( DELETELABEL == op){
	for(z=0;z<zs;z++)
		for(x=0;x<xs;x++)
			for(y=0;y<ys;y++)
				if(orgImg[z][x][y] == labelVal)
					outImg[z][x][y] = 0;
				else
					outImg[z][x][y] = orgImg[z][x][y];
}
if( SCALING == op){
	for(z=0;z<zs;z++)
		for(x=0;x<xs;x++)
			for(y=0;y<ys;y++)
				outImg[z][x][y] = orgImg[z][x][y]*scales;
}

if( PICKLABEL == op){
	for(z=0;z<zs;z++)
		for(x=0;x<xs;x++)
			for(y=0;y<ys;y++)
				if(orgImg[z][x][y] != 0){
					/*for(i = 0; i < labelsToPickCnt; i++) { */
						if(orgImg[z][x][y] == labelsToPick)
							outImg[z][x][y] = orgImg[z][x][y];
							/*break; */
							 else 
							outImg[z][x][y] = 0;
							/*}*/
				}
}
if( EQUALIZATION == op){
	int i, R;
	int left[MAXGRAY], width[MAXGRAY];
	long total, Hsum, Havg, histo[MAXGRAY];
	
	printf("MAXGRAY = %d\n",MAXGRAY);
	for(i=0; i < MAXGRAY; i++) 
		histo[i] = 0;
	for(z = 0; z<zs;z++)
		for(x=0;x<xs;x++)
			for(y=0;y<ys;y++){
				histo[(int)orgImg[z][x][y]]++;
				if(orgImg[z][x][y])total++;
			}
	histo[0] = 0;
	R    = 0;
	Hsum = 0;
	Havg = total / MAXGRAY;

	for(i=0; i < MAXGRAY; i++) {
		left[i] = R;
		Hsum += histo[i];
		while(Hsum>Havg && R<MAXGRAY){
			Hsum -= Havg;
			R++;
		}
		width[i] = R - left[i] + 1;
	}
	for(z = 0; z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				if(width[(int)orgImg[z][x][y]] == 1) 
					outImg[z][x][y] = left[(int)orgImg[z][x][y]];
				else {
					R = ((rand()&0x7fff)*width[(int)orgImg[z][x][y]]) >> 15;
					outImg[z][x][y] = left[(int)orgImg[z][x][y]] + R;
				}
			}
		}
	}/*
	for(z = 0; z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				outImg[z][x][y] = outImg[z][x][y] * 255.0 / (1.0 * left[MAXGRAY-1]);
			}
		}
	}
	*/
}
if( PICKING == op){
	  if(pick_high>zs){
		  pick_high = zs;
		  sliceChgNum = (int)ceilf(((float)(pick_high - pick_low )/(float)pick_jump));
		  ozs = sliceChgNum;
	  }
	if(verboseFlag)printf("picking slices %d %d %d %d\n",pick_low,pick_high,pick_jump,sliceChgNum);
	for(z = 0; z < sliceChgNum; z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				/*outImg[z][x][y] = orgImg[zs-1-(pick_low+pick_jump*z)][x][y];
				*/outImg[z][x][y] = orgImg[pick_low+pick_jump*z][x][y];
			}
		}
	}
}
if ( MASKAVERAGE == op){
	int voxelCnt = 0;
	float sum = 0.0;
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				if(subImg[z][x][y] == maskVal && orgImg != 0){
					sum += orgImg[z][x][y];
					voxelCnt++;
				}
			}
		}
	}
	/* this is for GM normalizing 
	sum /= voxelCnt;
	printf("%f\n",75.0/sum);
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				outImg[z][x][y] = (float)orgImg[z][x][y] * 75.0/ sum;
			}
		}
	}
	 GM normalizing done */
	/* this is for masked nbrPt */
	sum /= voxelCnt;
	if(verboseFlag)printf("%d ",voxelCnt);
	printf("%f\n",sum);
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				outImg[z][x][y] = (datatype)orgImg[z][x][y] * 1.0;
			}
		}
	}
	/* masked nbrPt averaging done */ 
}
			
if( SHIFTING == op){
	if('H' == shiftDirection){printf("Horizonal shifting %d",shiftAmount);
		for(z=0;z<zs;z++)
			for(x=0;x<xs;x++)
				for(y=0;y<ys;y++)
					if(y-shiftAmount >= 0 && y-shiftAmount < ys )
						outImg[z][x][y] = orgImg[z][x][y-shiftAmount];
					else if(y-shiftAmount >= ys)
						outImg[z][x][y] = orgImg[z][x][y-shiftAmount-ys];
					else
						outImg[z][x][y] = orgImg[z][x][y-shiftAmount+ys];
	}else if ('V' == shiftDirection){printf("vertical shifting %d",shiftAmount);
		for(z=0;z<zs;z++)
	            for(x=0;x<xs;x++)
	                for(y=0;y<ys;y++)
			if(x-shiftAmount >= 0 && x-shiftAmount < xs )		
	                        outImg[z][x][y] = orgImg[z][x-shiftAmount][y];
			else if (x-shiftAmount >= xs)
				outImg[z][x][y] = orgImg[z][x-shiftAmount-xs][y];
			else
				outImg[z][x][y] = orgImg[z][x-shiftAmount+xs][y];
	}else if('D' == shiftDirection){printf("the program will decide the horizontal shifting\n");
		shiftAmount = getShiftAmount(orgImg, xs, ys, zs);
		if(shiftAmount){
			printf("Calculated horizontal shifting is %d\n",shiftAmount);
		}else{
			printf("I did not think you need shift the image volume, so no shifting was did to the original image volume\n");
		}
		for(z=0;z<zs;z++)
			for(x=0;x<xs;x++)
				for(y=0;y<ys;y++)
					if(y-shiftAmount >= 0 && y-shiftAmount < ys )
						outImg[z][x][y] = orgImg[z][x][y-shiftAmount];
					else if(y-shiftAmount >= ys)
						outImg[z][x][y] = orgImg[z][x][y-shiftAmount-ys];
					else
						outImg[z][x][y] = orgImg[z][x][y-shiftAmount+ys];
	}
}
if( PADDING == op){
	if(paddingAmountY/2*paddingAmountX/2 < 0){
		printf(" padding value must greater than zero!!\n");
		exit(1);
	}
	if(verboseFlag)
		printf(" padding !!\n");
	for(z=0;z<zs;z++)
		for(x = paddingAmountX/2; x < xs + paddingAmountX/2; x++)
			for(y = paddingAmountY/2 ; y < ys + paddingAmountY/2 ; y++){
				outImg[z][x][y] = orgImg[z][x-paddingAmountX/2][y-paddingAmountY/2];
			}
}

if( CUTTING == op){
	if(cuttingAmountX/2 * cuttingAmountY/2 <0){
		printf(" cutting value must greater than zero!!\n");
		exit(EXIT_FAILURE);
	}
	if(cuttingAmountX/2-cuttingOffsetX < 0 || y+cuttingAmountY/2-cuttingOffsetY < 0){
		printf(" cutting and offset value not matching\n");
		exit(EXIT_FAILURE);
	}
	for(z=0;z<zs;z++)
		for(x = 0; x < oxs; x++)
			for(y = 0 ; y < oys ; y++)
				outImg[z][x][y] = orgImg[z][x+cuttingAmountX/2-cuttingOffsetX][y+cuttingAmountY/2-cuttingOffsetY];
}

if( MERGE == op){
	for(z=0;z<zs;z++)
		for(x=0;x < oxs;x++)
			 for(y=0;y < oys;y++)
				 if (orgImg[z][x][y]*subImg[z][x][y] == 0)
					 outImg[z][x][y] = orgImg[z][x][y]+subImg[z][x][y];
				 else
					 outImg[z][x][y] = (orgImg[z][x][y]+subImg[z][x][y])/2;
}
if( BIASCORRECTION == op){
	float *biasfield;
	biasfield=malloc(xs*ys*zs*sizeof(float));
	for(k=0;k<xs*ys*zs;k++){
		biasfield[k] = 0.0;
	}
	sprintf(filename,argv[2]);
	fp=myopen(filename,"r");
	for(k=0;k<xs*ys*zs;k++){
		fread(biasfield,sizeof(float),zs*xs*ys,fp);
	}
	fclose(fp);
	for(z=0;z<zs;z++)
		for(x=0;x < oxs;x++)
			for(y=0;y < oys;y++)
				outImg[z][x][y] = (int)(orgImg[z][x][y] * biasfield[z*xs*ys + x * ys + ys]);

}

if( THRESH == op){
	 for(z=0;z<zs;z++)
		 for(x=0;x < oxs;x++)
			 for(y=0;y < oys;y++)
				 if (orgImg[z][x][y] <= lowThreshold)
					 outImg[z][x][y] = (threshedVal[0]<0?orgImg[z][x][y]:threshedVal[0]);
				 else if (orgImg[z][x][y] <= highThreshold)
					 outImg[z][x][y] = (threshedVal[1]<0?orgImg[z][x][y]:threshedVal[1]);
				 else if (orgImg[z][x][y] > highThreshold)
					 outImg[z][x][y] = (threshedVal[2]<0?orgImg[z][x][y]:threshedVal[2]);
				 else{
					 printf("Error in parsing threshold,exist\n");
					 exit(-1);
				 }
}

if( FLIP == op){
	if(0 == flipDirection){
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0; y<ys; y++) {
					outImg[z][x][y] = orgImg[z][xs-x-1][y];
				}
			}
		}
	}else if( 1 == flipDirection){
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0; y<ys; y++) {
					outImg[z][x][y] = orgImg[z][x][ys-y-1];
				}
			}
		}
	}else if( 2 == flipDirection){
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0; y<ys; y++) {
					outImg[z][x][y] = orgImg[zs-z-1][x][y];
				}
			}
		}
	}else if( 5 == flipDirection){
		ucnt = 0;lcnt=0;
		int histLen = (1<<8*sizeof(datatype));
		int hist[histLen];
		int nonZero = 0;
		int histAcc = 0;
		int lowerThresh = 0;
		for(i = 0; i < histLen; i++){
			hist[i] = 0;
		}
		for(z=0; z<zs;z++) {
			for(x=0; x<xs; x++ ) {
				for(y=0;y<ys;y++) {
					if( orgImg[z][x][y] > 5){
						hist[(int)orgImg[z][x][y]]++;
						nonZero++;
					}
				}
			}
		}
		for(i = 1; i < histLen; i++){
			histAcc += hist[i];
			if(histAcc >= (nonZero * 0.98)){
				lowerThresh = i;
				i += histLen;
			}
		}
				
		for(z = 0; z<30; z++){
			for(x=(int)(xs-xs*2.0/3.0); x<(int)(xs-xs/3.0); x++ ) {
				for(y=0; y<ys; y++) {
					if(orgImg[z][x][y]>lowerThresh)ucnt++;
				}
			}
		}
		for(z = zs-30; z<zs; z++){
			for(x=(int)(xs-xs*2.0/3.0); x<(int)(xs-xs/3.0); x++ ) {
				for(y=0; y<ys; y++) {
					if(orgImg[z][x][y]>lowerThresh)lcnt++;
				}
			}
		}
		if(verboseFlag)printf("%d,%d with thres %d\n",lcnt,ucnt,lowerThresh);
		if(ucnt>lcnt){
			printf("1\n");
			for(z=0; z<zs;z++) {
				for(x=0; x<xs; x++ ) {
					for(y=0; y<ys; y++) {
						outImg[z][x][y] = orgImg[zs-z-1][x][y];
					}
				}
			}
		}else{
			for(z=0; z<zs;z++) {
				for(x=0; x<xs; x++ ) {
					for(y=0; y<ys; y++) {
						outImg[z][x][y] = orgImg[z][x][y];
					}
				}
			}
		}
	}else{
		printf("ERROR in flipping\n");
		exit(1);
	}
}

if( MASKVAL == op){
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				if(valsToMask == orgImg[z][x][y] && subImg[z][x][y] > 0){
					outImg[z][x][y] = maskVal;
				}else{
					outImg[z][x][y] = orgImg[z][x][y];
				}
			}
		}
	}
}
if( ADD == op){
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
//				if(subImg[z][x][y] > 0)
					outImg[z][x][y] = orgImg[z][x][y]>subImg[z][x][y]?orgImg[z][x][y]:subImg[z][x][y];
//				else
//					outImg[z][x][y] = orgImg[z][x][y];
//				outImg[z][x][y] = ( (orgImg[z][x][y] + subImg[z][x][y]) < c_maxValue?(orgImg[z][x][y] + subImg[z][x][y]):c_maxValue);
			}
		}
	}
}
if( EMPTY == op){
	 for(z=0;z<sliceChgNum;z++){
		 for(x=0;x<xs;x++){
			 for(y=0;y<ys;y++){
				 outImg[z][x][y] = 0;
			 }
		 }
	 }
}
if( MAXTOZERO == op){
	float fmax = 0.0;
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
//				printf(" %f",orgImg[z][x][y]);
				if(orgImg[z][x][y] > fmax){
					fmax = orgImg[z][x][y];
				}
			}
		}
	}
	/*
	printf("fmax = %f\n",fmax);
	for(z=0;z<zs;z++){
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				if(fmax - orgImg[z][x][y] <= 100.000001){
					outImg[z][x][y] = 0;
				}else{
					outImg[z][x][y] = orgImg[z][x][y];
				}
			}
		}
	}*/
}
for(z=0;z<zs;z++){
	for(x=0;x<xs;x++){
		for(y=0;y<ys;y++){
			if( REPLACE == op) 
				outImg[z][x][y] = (subImg[z][x][y] > 0?subImg[z][x][y]:orgImg[z][x][y]);
			else if( MASK == op)
				if(-1 == maskout){
					outImg[z][x][y] = (((orgImg[z][x][y] > 0) && (subImg[z][x][y] > 0))?0:orgImg[z][x][y]);
				} else {
					outImg[z][x][y] = (((orgImg[z][x][y] > 0) && (subImg[z][x][y] > 0))?orgImg[z][x][y]:0);
				}
			else if( MASKHOLE == op){
				if((orgImg[z][x][y] == 0) && (subImg[z][x][y] > 0)){
					outImg[z][x][y] = subImg[z][x][y];
				}else{
					outImg[z][x][y] = orgImg[z][x][y];
				}
			}
			else if( SHIFTVAL == op)
				outImg[z][x][y] = ((shiftVal <= orgImg[z][x][y])?orgImg[z][x][y]-shiftVal:orgImg[z][x][y]);
			else if( INVERSE == op)
				outImg[z][x][y] = c_maxValue - orgImg[z][x][y];
			else if( AND == op){
				outImg[z][x][y] = (subImg[z][x][y] > 0?orgImg[z][x][y]:0);
			}else if( SUBTRACT == op){
			//	outImg[z][x][y] = (int)(abs((orgImg[z][x][y] - subImg[z][x][y]) * subtractScale)); 
				if(/*0 != subImg[z][x][y] && */orgImg[z][x][y] != 0){
					floatOutput[z*xs*ys+x*ys+y] = (orgImg[z][x][y] - subImg[z][x][y]) * subtractScale;
				}else{
					floatOutput[z*xs*ys+x*ys+y] = 0.0;
				}
/*				outImg[z][x][y] = abs(orgImg[z][x][y] - subImg[z][x][y]);
 */
/*				outImg[z][x][y] = (orgImg[z][x][y] - subImg[z][x][y]) > 0? abs(orgImg[z][x][y] - subImg[z][x][y]):0; */
			}else if( DIVIDE == op) {
				if(0 != subImg[z][x][y] && orgImg[z][x][y] != 0){
					//outImg[z][x][y] = (datatype)(100.0 * ((orgImg[z][x][y]-subImg[z][x][y])*1.0/subImg[z][x][y]*1.0)/divideScale);
					floatOutput[z*xs*ys+x*ys+y] = (float)(100.0 * ((orgImg[z][x][y]-subImg[z][x][y])*1.0/subImg[z][x][y]*1.0)/divideScale);
					/*
					if(floatOutput[z*xs*ys+x*ys+y] < 0){
						printf("org[%d][%d][%d] = %d,sub[%d][%d][%d] = %d,floatOutput=%f\n",z,x,y,orgImg[z][x][y],z,x,y,subImg[z][x][y],floatOutput[z*xs*ys+x*ys+y]);
					}*/
				//	if(outImg[z][x][y] >0)printf("%f ",outImg[z][x][y]);
					/*
					fTemp1=orgImg[z][x][y] * 1.0;
					fTemp2=subImg[z][x][y] * 1.0;
					outImg[z][x][y] = (datatype)(100.0 * ((fTemp1/fTemp2)-1)/divideScale);
					*///printf("%f ",((fTemp1/fTemp2-1.0))/divideScale);
				}
				else
					floatOutput[z*xs*ys+x*ys+y] = 0.0;
			}
		}
	}
}
if(DIVIDE == op){
	printf("org[z][138][90] = %d,sub[2][138][90] = %d,floatOutput=%f\n",orgImg[87][148][120],subImg[87][148][120],floatOutput[87*148*120]);
}
/****** save output image *****/
if(verboseFlag){
	printf("save output image\n");
}
if( ORIENTATION == op){
	sprintf(filename,argv[3]);
	fp=myopen(filename,"w");
	for(k=0;k<zs;k++){
		for(i=0;i<xs;i++){
			fwrite(outImg[k][i],sizeof(datatype),ys,fp);
		}
	}
	fclose(fp);
}else if(DIVIDE == op || SUBTRACT == op){
	sprintf(filename,argv[3]);
	fp=myopen(filename,"w");
	fwrite(floatOutput,sizeof(float),zs*xs*ys,fp);
	fclose(fp);
}else {
	sprintf(filename,argv[3]);

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

return 0;

}

int getShiftAmount(datatype ***orgImg, int xs, int ys, int zs)
{
	int shiftAmountY = ys;
	const c_maxValue = 1<< (8*sizeof(datatype) - 1);
	int min= 2 *c_maxValue*xs;
	int x=0,y=0,z=0,sum=0;
	for(z=zs/2; z<zs/2+1;z++) {
		for(y=ys/2; y<ys; y++ ) {
			sum = 0;
			for(x=0;x<xs;x++) {
				sum +=orgImg[z][x][y];
			}
			if(sum < min){
				min = sum;
				shiftAmountY = y;
			}
		}
	}
	if(ys-shiftAmountY < 2)
		return 0;
	else
		return ys-shiftAmountY;
}

void fillHole(datatype ***orgImg, int fillHoleSize, datatype ***outImg, int xs, int ys, int zs)
{
	int x,y,z,curLabel = 1;
	int c_maxValue = 0;
	datatype ***flagImg;

	if(1 == sizeof(datatype)){
		c_maxValue = 255;
	}else if(4 == sizeof(datatype)){
		c_maxValue = (datatype)0xFFFFFFFF;
	}else if( 32768 == (1<< (8 * sizeof(datatype)-1))){
		c_maxValue = (datatype)65535;
	}else{
		c_maxValue = (datatype)32767;
	}

	flagImg = (datatype ***)MyAlloc3d(sizeof(datatype), xs, ys, zs);
	if(NULL == flagImg){
		printf("memory allocation casted error\n");
		exit(-1);
	}
	/*
	if(2 == sizeof(datatype)){
		flagImg = (datatype ***)USalloc3d(xs,ys,zs);
	}else{
		flagImg = (datatype ***)UCalloc3d(xs,ys,zs);
	}
	*/
	for(z = 0;z<zs;z++){
		for(x = 0; x<xs;x++){
			for(y =0; y<ys;y++){
				flagImg[z][x][y] = 0; /*init*/
			}
		}
	}
	

	/*mark all the holes*/
	
	for(z = 0;z<zs;z++){
		for(y =0; y<ys;y++){
			if(( 0 == orgImg[z][0][y] ) && (0 == flagImg[z][0][y-1] )){
			}
		}
	}
			
	for(z = 0;z<zs;z++){
		for(x = 0; x<xs;x++){
			for(y =0; y<ys;y++){
				if((0 == flagImg[z][x][y]) && (0 == orgImg[z][x][y]) && (x - 1 > 0) && (y -1 > 0) && (x < c_maxValue) && (y < c_maxValue)&& ( 0 ==flagImg[z][x][y-1]) && (0 ==flagImg[z][x-1][y-1]) && (0 ==flagImg[z][x-1][y]) && (0 ==flagImg[z][x-1][y+1])) {
					flagImg[z][x][y] = ++curLabel;
				}else if((0 == flagImg[z][x][y]) && (0 == orgImg[z][x][y]) && (x - 1 > 0) && (y -1 > 0) && (x < c_maxValue) && (y < c_maxValue)&& ( flagImg[z][x][y-1] != 0)){flagImg[z][x][y] = flagImg[z][x][y-1];
				}else if((0 == flagImg[z][x][y]) && (0 == orgImg[z][x][y]) && (x - 1 > 0) && (y -1 > 0) && (x < c_maxValue) && (y < c_maxValue)&& ( flagImg[z][x-1][y-1] != 0)){ flagImg[z][x][y] = flagImg[z][x-1][y-1]; 
				}else if((0 == flagImg[z][x][y]) && (0 == orgImg[z][x][y]) && (x - 1 > 0) && (y -1 > 0) && (x < c_maxValue) && (y < c_maxValue)&& ( flagImg[z][x-1][y] != 0)){flagImg[z][x][y] = flagImg[z][x-1][y];
				}else if((0 == flagImg[z][x][y]) && (0 == orgImg[z][x][y]) && (x - 1 > 0) && (y -1 > 0) && (x < c_maxValue) && (y < c_maxValue)&& ( flagImg[z][x-1][y+1] != 0)){flagImg[z][x][y] = flagImg[z][x-1][y+1];
				}
			}
		}
	}
}

void showUsage(char **argv) {
   printf("USAGE: %s orgImg subImg outpImg (if only one input, duplicate it as second input) [options] \n\
\t  -x <int>               : height size:256\n\
\t  -y <int>               : width size:256 \n\
\t  -z <int>               : slice number, if given, generally there are headers in the volume \n\
\t  -a                     : add the first two images\n\
\t  -A                     : merge the first two images(if overlap, divided by 2)\n\
\t  -s <scale>             : subtract the frist two images and scale the result\n\
\t  -S H[VD],<int>         : shift the input image H(orizontally),V(ertically) Left(negtive),Right(positive),Up(Negtive),Down(Positive) using circle shifting, or D(ecide) by the program)\n\
\t  -m <int>               : mask the first input, if -1, keep the non-mask part, else the org part  \n\
\t  -M <valTomask>,<maskedval> : mask the value(s) to assigned value \n\
\t  -h                     : mask the image hole to mask value \n\
\t  -T <int lowThresh>,<int highThresh>,<int lowerVal>,<int midVal>,<int highVal>\n\
\t                         : threshold the image to by <lowThresh:0> <highThresh:255>\n\
\t                           and assign val( -1 to keep original value),default is 0,255,0,-1,0,\n\
\t                           x<=lowThresh x<=highThresh, x>highThresh \n\
\t  -r                     : replace the first image using the 2nd one \n\
\t  -O <zxy>      : change the orientation of the input image as zxy,xyz,xzy,...[-]means image reading direction\n\
\t  -l <int>               : pick one label in the input image \n\
\t  -n                     : pixel AND-like operation \n\
\t  -e <int>               : create an empty image volume with slice number of <int>\n\
\t  -E                     : input image volume endian type(1-big endian, 0-little endian)\n\
\t  -f <int>               : fill holes in binarized image \n\
\t  -F <int>               : flip the input volume in direction row(0),column(1),depth(2)\n\
\t  -C <int>               : change the slice number to # from the end\n\
\t  -P <int>,<int>,<int>   : pick the slices from pick_low(0) to pick_high(include border slices)using jump\n\
\t  -i                     : reverse the intensity(255-i) \n\
\t  -p <int>,<int>	       : padding the input image in height,width of <int>, <int> more \n\
\t  -c <int>,<int>	       : cutting the input image in height,width of <int>, <int> less \n\
\t  -d <int>,<int>	       : downsampling the input image in height,width with rate of <int>,<int> \n\
\t  -D <int>               : delete label <int>\n",argv[0]);
exit(1);

}
