#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include <strings.h>
#include <mvcd.h>
#include <cres.h>
#include <matrixSHEN.h>  /*by SHEN*/


#define WM      250
#define GM      150
#define CSF      10
#define NNO       0
#define YYES      1


typedef struct 
{ 
  unsigned char v1 ;
  unsigned char v2 ; 
  unsigned char v3 ; 
} mask ; 


int  image_size;
mask MaskValue ;
char SubjSegImgFile[80] ;


void show_usage_SHEN() ;
void MaskoutImg(int argc, char *argv[]) ;
void Maskout2ByteImg(int argc, char *argv[]) ;


int main(int argc,char *argv[])
{
  int    c,num;
  extern char *optarg;
  int v1,v2,v3;
  int ByteImg ;
 

  num=4;
 
  if(argc<num)
    show_usage_SHEN() ;

  image_size=256;
  MaskValue.v1 = 0 ;
  MaskValue.v2 = 0 ;
  MaskValue.v3 = 0 ;
  ByteImg = NNO ;

  while((c=getopt(argc-3,argv+3,"s:v:b")) != -1)
    {
      switch(c)
	{
	case 's':
	  image_size=atoi(optarg);
	  break;
	
	case 'v':
	  sscanf(optarg, "%d,%d,%d", &v1, &v2, &v3) ;
	  MaskValue.v1 = v1 ; MaskValue.v2 = v2 ; MaskValue.v3 = v3 ;
	  break ;

	case 'b':
	  ByteImg = YYES ;
	  break;
	}
    }
  printf(" image_size=%d\n", image_size) ;
  printf(" MaskValue.v1=%d MaskValue.v2=%d\n MaskValue.v3=%d\n", MaskValue.v1, MaskValue.v2, MaskValue.v3) ;


  if( ByteImg==YYES )
    Maskout2ByteImg(argc, argv) ;
  else
    MaskoutImg(argc, argv) ;
  return 0;
}


void show_usage_SHEN()
{
  printf("USAGE: Maskout3Dimg <inputed_data> <mask img> <output image> \n\
\t -s <int>             : image size in x and y dimensions\n\
\t -v <int,int,int>     : voxels with these three values in the mask image will be used remove the content in the input image in the same voxels. CSF=10 \n\
\t -b                   : the input images are at two bytes\n\
");
  exit(1);
}

void WriteImg(char filename[80], unsigned char ***data, int image_size, int z_size)
{
  FILE  *fp;
  int   i, k ;

  /* write the smoothed image */
  fp=myopen(filename,"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(data[k][i],1,image_size,fp);
  fclose(fp);
}

void Write2ByteImg(char filename[80], short ***data, int image_size, int z_size)
{
  FILE  *fp;
  int   i, k ;

  /* write the smoothed image */
  fp=myopen(filename,"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(data[k][i],2,image_size,fp);
  fclose(fp);
}

void MaskoutImg(int argc, char *argv[])
{
  unsigned char ***InputImg, ***MaskImg ;
  int   i,j,k,l,m, z_size ;
  FILE  *fp;
  char  filename[80] ;


  /***** allocate memory and read input image *****/
  fp=myopen(argv[1],"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);  printf("z_size=%d\n", z_size) ;
  rewind(fp);
  
  InputImg = UCalloc3d(image_size,image_size,z_size);
  MaskImg  = UCalloc3d(image_size,image_size,z_size);
    
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fread(InputImg[k][i],1,image_size,fp);
  fclose(fp);


  /***** allocate memory and read mask image *****/
  fp=myopen(argv[2],"r");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fread(MaskImg[k][i],1,image_size,fp);
  fclose(fp);

  /* check difference */
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	if( MaskImg[k][i][j]==MaskValue.v1 || MaskImg[k][i][j]==MaskValue.v2 || MaskImg[k][i][j]==MaskValue.v3)
	  InputImg[k][i][j]=0;

  /* write result */
  WriteImg(argv[3], InputImg, image_size, z_size) ;


  /* free */
  UCfree3d(InputImg, z_size, image_size) ;
  UCfree3d(MaskImg,  z_size, image_size) ;
}



void Maskout2ByteImg(int argc, char *argv[])
{
  short ***InputImg, ***MaskImg ;
  int   i,j,k,l,m, z_size ;
  FILE  *fp;
  char  filename[80] ;


  /***** allocate memory and read input image *****/
  fp=myopen(argv[1],"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size*2);  printf("z_size=%d\n", z_size) ;
  rewind(fp);
  
  InputImg = Salloc3d(image_size,image_size,z_size);
  MaskImg  = Salloc3d(image_size,image_size,z_size);
    
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fread(InputImg[k][i],2,image_size,fp);
  fclose(fp);


  /***** allocate memory and read mask image *****/
  fp=myopen(argv[2],"r");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fread(MaskImg[k][i],2,image_size,fp);
  fclose(fp);

  /* check difference */
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	if( MaskImg[k][i][j]==MaskValue.v1 || MaskImg[k][i][j]==MaskValue.v2 || MaskImg[k][i][j]==MaskValue.v3)
	  InputImg[k][i][j]=0;

  /* write result */
  Write2ByteImg(argv[3], InputImg, image_size, z_size) ;


  /* free */
  Sfree3d(InputImg, z_size, image_size) ;
  Sfree3d(MaskImg,  z_size, image_size) ;
}
