#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include <strings.h>
#include <mvcd.h>
#include <cres.h>

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

int image_size, z_size ;
int nbrSize, Z_nbrSize ;


/* main part of region growing*/
/*void NeighboringSimilarPointsSearch( long int *PointNum, int ***Detected_Region, int order, unsigned char ***src, unsigned char ***status, 
				     int i, int j, int k, int Greythreshold ) 
{
  int kk,ii,jj; 
  int k_plus_nbrSize, k_minus_nbrSize, i_plus_nbrSize, i_minus_nbrSize, j_plus_nbrSize, j_minus_nbrSize, x, y;


  k_plus_nbrSize =k+nbrSize;
  k_minus_nbrSize=k-nbrSize;
  i_plus_nbrSize =i+nbrSize;
  i_minus_nbrSize=i-nbrSize;
  j_plus_nbrSize =j+nbrSize;
  j_minus_nbrSize=j-nbrSize;

  if (k_plus_nbrSize>=z_size)     k_plus_nbrSize=z_size-1;
  if (k_minus_nbrSize<0)          k_minus_nbrSize=0;
  if (j_plus_nbrSize>=image_size) j_plus_nbrSize=image_size-1;
  if (j_minus_nbrSize<0)          j_minus_nbrSize=0;
  if (i_plus_nbrSize>=image_size) i_plus_nbrSize=image_size-1;
  if (i_minus_nbrSize<0)          i_minus_nbrSize=0;

  if(status[k][i][j]==FAIL) 
    {
      status[k][i][j] = OOKK ;

      (*PointNum) ++ ;      
      Detected_Region[k][i][j] = order ;

      for(kk=k_minus_nbrSize; kk<=k_plus_nbrSize; kk++)
	for(ii=i_minus_nbrSize; ii<=i_plus_nbrSize; ii++)
	  for(jj=j_minus_nbrSize; jj<=j_plus_nbrSize; jj++)
	    {
	      if( !(kk==k && jj==j && ii==i) && src[kk][ii][jj]>Greythreshold && status[kk][ii][jj]==FAIL )
		{
		  NeighboringSimilarPointsSearch( PointNum, Detected_Region, order, src, status, ii, jj, kk, Greythreshold ) ;
		}
	    }
    }
}*/
Ivector3d  *Stack ; 
unsigned char ***VoxelInStack ;
void NeighboringSimilarPointsSearch( long int *PointNum, int ***Detected_Region, int order, unsigned char ***src, unsigned char ***status, 
				     int i, int j, int k, int Greythreshold ) 
{
  int kk,ii,jj; 
  int k_plus_nbrSize, k_minus_nbrSize, i_plus_nbrSize, i_minus_nbrSize, j_plus_nbrSize, j_minus_nbrSize, x, y;

  /* stack */
  long int   Pointer ;


  /* add the first searched point */
  Pointer = 0 ; 
  Stack[Pointer].x = i; 
  Stack[Pointer].y = j; 
  Stack[Pointer].z = k; 


  do{
    i = Stack[Pointer].x ;
    j = Stack[Pointer].y ;
    k = Stack[Pointer].z ;
    Pointer -- ;

    status[k][i][j] = OOKK ;
    (*PointNum) ++ ;   
    Detected_Region[k][i][j] = order ;

    k_plus_nbrSize =k+Z_nbrSize; /* nbrSize; 2005 Jun */
    k_minus_nbrSize=k-Z_nbrSize; /* nbrSize; 2005 Jun */
    i_plus_nbrSize =i+nbrSize; 
    i_minus_nbrSize=i-nbrSize;
    j_plus_nbrSize =j+nbrSize;
    j_minus_nbrSize=j-nbrSize;

    if (k_plus_nbrSize>=z_size)     k_plus_nbrSize=z_size-1;
    if (k_minus_nbrSize<0)          k_minus_nbrSize=0;
    if (j_plus_nbrSize>=image_size) j_plus_nbrSize=image_size-1;
    if (j_minus_nbrSize<0)          j_minus_nbrSize=0;
    if (i_plus_nbrSize>=image_size) i_plus_nbrSize=image_size-1;
    if (i_minus_nbrSize<0)          i_minus_nbrSize=0;

    /* now we can check immediately adjacent points to see if they too could be added to the track */
    for(kk=k_minus_nbrSize; kk<=k_plus_nbrSize; kk++)
      for(ii=i_minus_nbrSize; ii<=i_plus_nbrSize; ii++)
	for(jj=j_minus_nbrSize; jj<=j_plus_nbrSize; jj++)
	  {
	    if( !(kk==k && jj==j && ii==i) && src[kk][ii][jj]>Greythreshold && status[kk][ii][jj]==FAIL && VoxelInStack[kk][ii][jj]==NNO )
	      {
		Pointer ++ ;
		Stack[Pointer].x = ii ;
		Stack[Pointer].y = jj; 
		Stack[Pointer].z = kk; 

		VoxelInStack[kk][ii][jj]= YYES ;
	      }
	  }
  }while(Pointer>=0) ;
}


void region_growing(int argc, char *argv[], int nbrSize, int Greythreshold, int RegionPointNum)
{
  unsigned char ***src, ***status; 
  int  ***Detected_Region, order ;
  int  i,j,k, ii,jj,kk ;
  long int PointNum ;
  FILE *fp;
  int  *PointNumInCertainOrder ;


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


  /* 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);
  status = UCalloc3d(image_size,image_size,z_size);
  Detected_Region   = Ialloc3d(image_size,image_size,z_size);

  PointNumInCertainOrder = Ialloc1d(image_size*image_size*1) ;

  /* 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);

  /* reset ... */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	status[k][i][j] = FAIL ;

  /* Global variables: for growing */	
   VoxelInStack = UCalloc3d(image_size,image_size,z_size);
   Stack = Ivector3dalloc1d(image_size*image_size*z_size) ;
   /* status for representing whether current voxel already in stack? */
   for(kk=0; kk<z_size; kk++)
     for(ii=0; ii<image_size; ii++)
       for(jj=0; jj<image_size; jj++)
	 VoxelInStack[kk][ii][jj] = NNO ;


  /* region growing ... */
  order = 1 ;
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  if(src[k][i][j]>Greythreshold && status[k][i][j]==FAIL) 
	    {
	      PointNum = 0 ;

	      printf("(%d,%d,%d) = %d\n", i,j,k, src[k][i][j]) ;

	      NeighboringSimilarPointsSearch( &PointNum, Detected_Region, order, src, status, i, j, k, Greythreshold ) ;
	      
	      if(PointNum>RegionPointNum )
		printf("PointNum=%d\n", PointNum) ;
	      PointNumInCertainOrder[order] = PointNum ;

	      order ++ ;
	    }
	}

  /* remove small connected regions */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  order = Detected_Region[k][i][j] ;
	  if( order>=1 && PointNumInCertainOrder[order]<RegionPointNum )
	    src[k][i][j] = 0 ;
	}

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


  /* free */
  UCfree3d(status, z_size, image_size) ;
  UCfree3d(src, z_size, image_size) ;
  Ifree3d(Detected_Region, z_size, image_size) ;
  free(PointNumInCertainOrder) ;

  /* Global variables: free */
  free(Stack);
  UCfree3d(VoxelInStack, z_size, image_size) ;
}



void show_usage()
{
  printf("USAGE: dsphere <src_file> <out_file>\n\
\t -n <int>               : the size defined for connection (1)\n\
\t -z <int>               : different search neighbor size for z direction \n\
\t -g <int>               : threshold for image grey-level value (default: 128)\n\
\t -t <int>               : keep the region with the point number large than (300)\n\
");
  exit(1);
}


int main(int argc,char *argv[])
{
  int   c, num ;
  int   Greythreshold, RegionPointNum ;
  int   separate_Z_nbrSize ;

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

  nbrSize=1 ;
  Greythreshold = 128 ; 
  RegionPointNum = 300 ;
  separate_Z_nbrSize = nbrSize ;

  while((c=getopt(argc-2,argv+2,"n:s:g:t:z:")) != -1)
    {
      switch(c)
	{
	case 'n':
	  nbrSize=atoi(optarg);
	  break;

	case 'g':
	  Greythreshold=atoi(optarg);
	  break;

	case 't':
	  RegionPointNum=atoi(optarg);
	  break;

	case 'z':
	  separate_Z_nbrSize=atoi(optarg);
	  break;

	default:
	  break;
	}
    }
  printf(" nbrSize=%d\n Greythreshold=%d\n RegionPointNum=%d\n \n", nbrSize, Greythreshold, RegionPointNum) ;

  Z_nbrSize = separate_Z_nbrSize ;
  printf("Z_nbrSize=%d\n", Z_nbrSize) ;

  region_growing(argc, argv, nbrSize, Greythreshold, RegionPointNum) ;

  return 0;
}
