#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*/


int image_size, head_data1, head_data2 ;
int SetFirstEdgeMapValueByThis ;

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



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

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

  image_size=256;
  head_data1 = 0 ;
  head_data2 = 0 ;
  SetFirstEdgeMapValueByThis = 254 ;
  while((c=getopt(argc-3,argv+3,"s:h:t:f:")) != -1)
    {
      switch(c)
	{
	case 's':
	  image_size=atoi(optarg);
	  break;
	
	case 'h':
	  head_data1=atoi(optarg);
	  break;
	
	case 't':
	  head_data2=atoi(optarg);
	  break;

	case 'f':
	  SetFirstEdgeMapValueByThis=atoi(optarg);
	  break;
	}
    }
  printf(" image_size=%d\n head_data1=%d\n head_data2=%d SetFirstEdgeMapValueByThis=%d\n", image_size, head_data1, head_data2, SetFirstEdgeMapValueByThis) ;

  StackTwoEdgeMaps(argc, argv) ;

  return 0;
}


void show_usage_SHEN()
{
  printf("USAGE: shen3D <inputed_data1> <inputed_data2> <output_data>\n\
\t -s <int>             : image size in x and y dimensions\n\
\t -h <int>             : head size one for image data one\n\
\t -t <int>             : head size two for image data two\n\
\t -f <int>             : the value that the first edge map will be set\n\
");

  printf("The first edge map will keep the same, while the second edge map will be halved and add into the place where the the first edge map has 0 value.\n") ;
  exit(1);
}


void StackTwoEdgeMaps(int argc, char *argv[])
{
  unsigned char ***data1, ***data2 ;
  int   i,j,k,l,m, z_size1, z_size2, z_size ;
  FILE  *fp;
  char  filename[280] ;


  /***** allocate memory and read image data1 *****/
  fp=myopen(argv[1],"rb");
  fseek(fp,0,SEEK_END);
  z_size1=(ftell(fp)-head_data1)/(image_size*image_size);  printf("z_size1=%d\n", z_size1) ;
  rewind(fp);
  
  data1=UCalloc3d(image_size,image_size,z_size1);
  
  fseek(fp,head_data1,SEEK_CUR);
  
  for(k=0;k<z_size1;k++)
    for(i=0;i<image_size;i++)
      fread(data1[k][i],1,image_size,fp);
  fclose(fp);


  /***** allocate memory and read image data1 *****/
  fp=myopen(argv[2],"rb");
  fseek(fp,0,SEEK_END);
  z_size2=(ftell(fp)-head_data2)/(image_size*image_size); printf("z_size2=%d\n", z_size2) ;
  rewind(fp);
  
  data2=UCalloc3d(image_size,image_size,z_size2);
  
  fseek(fp,head_data2,SEEK_CUR);
  
  for(k=0;k<z_size2;k++)
    for(i=0;i<image_size;i++)
      fread(data2[k][i],1,image_size,fp);
  fclose(fp);

  /* the slice numbers in two images may be different */
  z_size = z_size1 ;
  if(z_size2<z_size1) z_size = z_size2 ;

  /* stack */
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  /*  if( data1[k][i][j]==255) data1[k][i][j]=SetFirstEdgeMapValueByThis ;*/
	  if( data1[k][i][j]>=1) data1[k][i][j]=SetFirstEdgeMapValueByThis ;
	  else  /*if( data1[k][i][j]<1 ) April 2003*/
	    data1[k][i][j] = data2[k][i][j] ;  
	}

  /* saving ... */
  printf("attempting to save the results \n");
  sprintf(filename, argv[3]);    printf("%s\n", filename) ; 
  fp=myopen(argv[3],"wb"); 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(data1[k][i],1,image_size,fp);
  fclose(fp);


  /* free */
  UCfree3d(data1, z_size1, image_size) ;
  UCfree3d(data2, z_size2, image_size) ;
}
