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

#define Z_THICKINK    1.5 
#define FAIL   0
#define OOKK   1
#define NNO    0
#define YYES   1

#define BG    0
#define CSF   10
#define VN    50
#define GM    150
#define WM    250

#define ALLPLN 1658
#define X_PLN  11
#define Y_PLN  22
#define Z_PLN  33
#define Z_UP   110 /*110: use 3D region growing for z>Z_UP, others you can use 2D region growing */


/* image size */
int  image_size, z_size ;
int  nbrSize ;
Ivector3d    *MyHood ;
int          MyHoodSize, pixelNumIn_MyHood ;
int  GrowingOnPlan ;
int   OutPutTempoaryResult ;


float gaussian(float x, float s)
{
  return(exp((-x*x)/(2*s*s)));
}

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 calculate_hood_byIncreasingRadius(Ivector3d *Hood, int HoodSize, int *pixelNumInHood, float xres,float yres,float zres)
{
  int        i, j, k, l ;
  int        x, y, z, t, r  ;
  int        half ; 


  t = 0 ;  Hood[0].x = 0 ;  Hood[0].y = 0 ;  Hood[0].z = 0 ;
  t++ ;
  half = HoodSize/2;
  for(r=1; r<=half; r++)
    {
      for(z=-r; z<=r; z+=2*r)
	for(x=-r; x<=r; x++)
	  for(y=-r; y<=r; y++)
	    if((x*x+y*y+z*z)<=r*r) /* half*half)*/
	    {
	      Hood[t].x = x/xres ;
	      Hood[t].y = y/yres ;
	      Hood[t].z = z/zres ;
	      
	      t++ ;
	      /*printf("t=%d: (%d,%d,%d)\n", t,x,y,z) ;*/
	    }
      for(x=-r; x<=r; x+=2*r)
	for(z=-r+1; z<=r-1; z++) 
	  for(y=-r; y<=r; y++)
	    if((x*x+y*y+z*z)<=r*r)
	    {
	      Hood[t].x = x/xres ;
	      Hood[t].y = y/yres ;
	      Hood[t].z = z/zres ;
	      
	      t++ ;
	      /*printf("t=%d: (%d,%d,%d)\n", t,x,y,z) ;*/
	    }
      for(y=-r; y<=r; y+=2*r)
	for(z=-r+1; z<=r-1; z++) 
	  for(x=-r+1; x<=r-1; x++) 
	    if((x*x+y*y+z*z)<=r*r)
	    {
	      Hood[t].x = x/xres ;
	      Hood[t].y = y/yres ;
	      Hood[t].z = z/zres ;
	      
	      t++ ;
	      /*printf("t=%d: (%d,%d,%d)\n", t,x,y,z) ;*/
	    }
    }
  (*pixelNumInHood) = t ;
  printf("total %d in the search area!\n", t) ;
}



void OpeningOperationOnBinary3DImage(unsigned char ***outp, int Radius )
{
  unsigned char ***Temp_Img ;
  int  i, j, k ;
  int  l, m, n, index ;
  double TotalWhite ;
  int  NumInSphere, RealNum ;
  double r ;
  int  MyHoodSize ;

  NumInSphere = (2*Radius+1)*(2*Radius+1)*(2*Radius/Z_THICKINK+1) ;
  printf("\n\nopening operation on 3D binary image, using radius =%d\n\n", Radius) ;

  /* spherical mask */
  MyHoodSize = Radius*2+1 ;
  MyHood     = Ivector3dalloc1d(MyHoodSize*MyHoodSize*MyHoodSize) ;
  calculate_hood_byIncreasingRadius(MyHood, MyHoodSize, &pixelNumIn_MyHood, 1., 1., Z_THICKINK) ;


  /* apply for space for temporary image result */
  Temp_Img = UCalloc3d(image_size,image_size,z_size);

  /* find mask ... */
  for(k=0; k<z_size; k++)
    for(i=0; i<image_size; i++)
      for(j=0; j<image_size; j++)
	if( outp[k][i][j]==0 )
	{
	  TotalWhite = 0 ;
	  RealNum = 0 ;
	  for( index=0; index<pixelNumIn_MyHood ; index++)
	    /*for(l=-Radius; l<=Radius; l++)
	    for(m=-Radius; m<=Radius; m++) 
	    for(n=-Radius/Z_THICKINK; n<=Radius/Z_THICKINK; n++)*/
	    {
	      l = MyHood[index].x ;
	      m = MyHood[index].y ;
	      n = MyHood[index].z ;	  
	      if( l*l+m*m+n*n<Radius*Radius &&  i+l>=0 && i+l<image_size && j+m>=0 && j+m<image_size && k+n>=0 && k+n<z_size )
		{
		  TotalWhite += outp[k+n][i+l][j+m] ;
		  RealNum ++ ;
		  
		  /* for speeding up*/
		  if( outp[k+n][i+l][j+m]!=0 ) /* break */
		    index = pixelNumIn_MyHood ;
		}
	    }
	  TotalWhite /= 255.0 ;

	  if( TotalWhite==0 )
	    Temp_Img[k][i][j] = 255 ;
	  else
	    Temp_Img[k][i][j] = 0 ;
	}
	else
	  Temp_Img[k][i][j] = 0 ;

  /* operation ... */
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	outp[k][i][j] = 255-Temp_Img[k][i][j] ;

  
  /* free */
  UCfree3d(Temp_Img, z_size, image_size) ;
  free(MyHood) ;
}


void ClosingOperationOnBinary3DImage(unsigned char ***outp, int Radius )
{
  unsigned char ***Temp_Img ;
  int  i, j, k ;
  int  l, m, n, index ;
  double TotalWhite ;
  int  NumInSphere, RealNum ;
  int  MyHoodSize ;


  NumInSphere = (2*Radius+1)*(2*Radius+1)*(2*Radius/Z_THICKINK+1) ;
  printf("\n\nClosing operation on 3D binary image, using radius =%d\n\n", Radius) ;

  /* spherical mask */
  MyHoodSize = Radius*2+1 ;
  MyHood     = Ivector3dalloc1d(MyHoodSize*MyHoodSize*MyHoodSize) ;
  calculate_hood_byIncreasingRadius(MyHood, MyHoodSize, &pixelNumIn_MyHood, 1., 1., Z_THICKINK) ;


  /* apply for space for temporary image result */
  Temp_Img = UCalloc3d(image_size,image_size,z_size);


  /* find mask ... */
  for(k=0; k<z_size; k++)
    for(i=0; i<image_size; i++)
      for(j=0; j<image_size; j++)
	if( outp[k][i][j]==255 )
	{
	  TotalWhite = 0 ;
	  RealNum = 0 ;
	  for( index=0; index<pixelNumIn_MyHood ; index++)
	    /*for(l=-Radius; l<=Radius; l++)
	      for(m=-Radius; m<=Radius; m++) 
	      for(n=-Radius/Z_THICKINK; n<=Radius/Z_THICKINK; n++)*/
	    {
	      l = MyHood[index].x ;
	      m = MyHood[index].y ;
	      n = MyHood[index].z ;	  
	      if( l*l+m*m+n*n<Radius*Radius &&  i+l>=0 && i+l<image_size && j+m>=0 && j+m<image_size && k+n>=0 && k+n<z_size )
		{
		  TotalWhite += outp[k+n][i+l][j+m] ;
		  RealNum ++ ;
		  
		  
		  /* for speeding up*/
		  if( outp[k+n][i+l][j+m]!=255 ) /* break */
		    index = pixelNumIn_MyHood ;
		}
	    }
	  TotalWhite /= RealNum ;
	  
	  if( TotalWhite==255 )
	    Temp_Img[k][i][j] = 255 ;
	  else
	    Temp_Img[k][i][j] = 0 ;
	}
	else
	  Temp_Img[k][i][j] = 0 ;

  /* operation ... */
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	outp[k][i][j] = Temp_Img[k][i][j] ;

  
  /* free */
  UCfree3d(Temp_Img, z_size, image_size) ;
  free(MyHood) ;
}




void show_usage()
{
  printf("USAGE: extracthead <input_file> <output_file>\n\
\t -n <int>               : the radial size (3)\n\
\t -T <int>               : the size to extend (default: 0)\n\
");
  exit(1);
}


int main(int argc,char *argv[])
{
  int   c, num ;
  FILE  *fp;
  unsigned char ***src ;
  int  i, j, k ;
  char fn[180] ;
  int  OutExtend ;

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


  /* parameters */
  GrowingOnPlan = 0 ;
  nbrSize=3 ;
  OutPutTempoaryResult = YYES ;
  OutExtend=0;
  while((c=getopt(argc-2,argv+2,"n:T:")) != -1)
    {
      switch(c)
	{
	case 'n':
	  nbrSize=atoi(optarg);
	  break;

	case 'T':
	  OutExtend=atoi(optarg);
	  break;

	default:
	  break;
	}
    }

  printf("nbrSize=%d    OutExtend=%d\n", nbrSize, OutExtend) ;

  /* Open image */
   /* image size in xy plane */
   image_size=256;

   fp=myopen(argv[1],"r");
   fseek(fp,0,SEEK_END);
   z_size=ftell(fp)/(image_size*image_size);
   rewind(fp);
   
   /* apply for memory */ 
   src  = UCalloc3d(image_size,image_size,z_size);

   /* read data */ 
   for(k=0;k<z_size;k++)
     for(i=0;i<image_size;i++)
       fread(src[k][i],1,image_size,fp);
   fclose(fp);

   OpeningOperationOnBinary3DImage(src,  nbrSize) ;
   ClosingOperationOnBinary3DImage(src,  nbrSize) ;
   OpeningOperationOnBinary3DImage(src,  1+OutExtend) ;
   sprintf(fn, "%s", argv[2]);
   WriteImg(fn, src, image_size, z_size) ;
   
   return 0;
}

