/**
 * @file  extractventricle.cxx
 * @brief Originally part of HAMMER.
 *
 * Copyright (c) 2001, 2012, 2013 University of Pennsylvania.
 *
 * This file is part of DTI-DROID.
 *
 * DTI-DROID is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * DTI-DROID is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with DTI-DROID.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Contact: SBIA Group <sbia-software at uphs.upenn.edu>
 */

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

#include "mvcd.h"
#include "matrixSHEN.h"
#include "cres.h"

#include <dtidroid/basis.h>


// acceptable in .cxx file
using namespace dtidroid;


// ===========================================================================
// global variables and constants
// ===========================================================================

#define Z_THICKINK    2.0  /* 1.6 1.4 1.5 2.0*/
#define ORIENT_SCALE  24
#define FAIL   0
#define OOKK   1

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

int imageSz ;

// ===========================================================================
// help
// ===========================================================================

// ---------------------------------------------------------------------------
void show_usage()
{
  print_version("extractventricle", "", "");
  printf("\n");
  printf("USAGE: extractVentricle <input_file> <output_file>\n\
\t -a                     : half the grey-level valus of image\n\
\t -i                     : double the grey-level valus of image\n\
\t -h                     : reverse the balck point in image to be white\n\
\t -b <float>             : blur 3D image by Guassian filter (1.0)\n\
\t -p <float, int, int>   : blur 3D image by shape-based Guassian filter (1.0) and smoothing time (2), with strongness of intensity to focus (6=20/3). Here, the boundary voxels will not blurred by bacground\n\
\t -s <int>               : extract 3D image points with the grey-level value equal to (50)\n\
\t -l <int>,<int>         : extract 3D image points with the grey-level value around (50) with variation of 10\n\
\t -t <int>               : threshold the input image (100)\n\
\t -r <int>               : reflect 3D image by y=image_size/2 (128)\n\
\t -n <float>             : Guassian noise with variance (1.0) added to the input file\n\
\t -w <int>               : Remove nosie in the 3D edge map (edge: 1x1x1)\n\
\t -e <float>             : extract and scale 3D part at whole image center, by pixel-wise sampling (1.8)\n\
\t -d <float>             : extract and scale 3D part at (x,y,z) or image center, by linear interpolation (1.8)\n\
\t -g <float>             : Lapalacian edge detection, without smoothing (10)\n\
\t -f <float>             : Canny-liked edge detection, with Guassina smoothing (sigma: 1.5)\n\
\t -c <float>             : 3D detection by 2D Canny method (2.0)\n\
\t -k <float>             : 3D detection by 2D Canny method (2.0) at three directions\n\
\t -m                     : mapping intensities\n\
\t -C <int>,<int>,<int>   : crop one part of image around the image center (lengthX, lengthY, lengthZ)\n\
\t -U <int>               : smooth size \n\
\t -S <float, float, int> : smaple down image at rate sx,sz,image_size(i.e. 256)\n\
\t -A <int>               : number of image samples, used to calculate average image\n\
\t -R <int>               : remove content in slice ? \n\
\t -N                     : mapping intensities to 50, 150 or 250 based on nearest distance\n\
\t -Z <char>,<int><float> : crop and zoom certain slice, such as x,128,3.0\n\
\t -L <int, float>        : level to increase on each pixel (i.e. 10), and gamma value (1.0)\n\
\t -P <int>,<int>         : to pick up intensity range <min, max> and set them to 255, others 0\n\
\t -Q <int>,<int>         : to pick up intensity range <min, max> \n\
\t -B                     : save coordinates only in WM,GM,CSF\n\
\t -I <int>               : increase slice number to <int>\n\
\t -M                     : find max and min intensity from 3D image\n\
\t -H                     : total number of WM and GM in brain image\n\
\t -D                     : histogram\n\
\t -V                     : total intensities in the image, used to compute total changes in certain part\n\
\t -W <int>               : Get Volume Of Voxels With A Certain Label\n\
\t -F <int>               : Increase a certain value for each non-zero voxel; default: 0\n\
");
}

// ===========================================================================
// auxiliary functions
// ===========================================================================

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

void Smoothing3DimageWithShapeBasedGuassian(int argc, char *argv[], float Smooth_Sigma, int SmoothTimes, int FocusIntensity)
{
  unsigned char ***src ;
  int  i,j,k,l, ii,jj,kk, image_size, z_size ;
  FILE *fp;
  float ***Density, ***outp ;\
  float src_nonZeroNum, src_TotalIntensities, Current_nonZeroNum, Current_TotalIntensities ;

  /* variables for Guassian filetering */
  int    filter_width;               /* length of 1-D gaussian mask */
  float a,b,c,d;        /* mask generation intermediate vars*/
  float ti,tj,tk,total,weights;
  float  r ;
  int   efficient ;
  float ***G, TotalWeight ;
  int   num ;
  float filterCoefficientThreshold ;



  filterCoefficientThreshold = 0.1 ; /* 0.005 */
  printf("\n\nStarting ...\nSmooth_Sigma=%f\nfilterCoefficientThreshold=%f\n", Smooth_Sigma, filterCoefficientThreshold) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp   = Falloc3d(image_size,image_size,z_size);
  Density= Falloc3d(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);

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

	 
  /* calc coeffs for 1-dimensional G, dG/dn and Delta-squared G filters */
  for(i=0; i<100; i++) 	
    {
      a=gaussian(i, Smooth_Sigma);
      if(a<filterCoefficientThreshold) {  filter_width=i; break; }
    }
  filter_width -- ;
  printf("filter_width=%d\n", filter_width) ;


  G = Falloc3d(filter_width*2+1, filter_width*2+1, filter_width*2+1); 
  for(ii=-filter_width; ii<=filter_width; ii++)
    for(jj=-filter_width; jj<=filter_width; jj++)
      for(kk=-filter_width; kk<=filter_width; kk++)
	{
	  r = sqrt(ii*ii+jj*jj+kk*kk) ;
	  
	  a=gaussian(r, Smooth_Sigma);
	  if(a>filterCoefficientThreshold || r<2.)
	    G[ii+filter_width][jj+filter_width][kk+filter_width]=a/(6.283185*Smooth_Sigma*Smooth_Sigma);
	}

  /* normalize ... */
  weights = 0 ;
  for(ii=-filter_width; ii<=filter_width; ii++)
    for(jj=-filter_width; jj<=filter_width; jj++)
      for(kk=-filter_width; kk<=filter_width; kk++)
	weights += G[ii+filter_width][jj+filter_width][kk+filter_width] ; 

  for(ii=-filter_width; ii<=filter_width; ii++)
    for(jj=-filter_width; jj<=filter_width; jj++)
      for(kk=-filter_width; kk<=filter_width; kk++)
	G[ii+filter_width][jj+filter_width][kk+filter_width] /= weights ; 

  /* print */
  for(ii=-filter_width; ii<=filter_width; ii++)
    {
      for(jj=-filter_width; jj<=filter_width; jj++)
	{
	  for(kk=-filter_width; kk<=filter_width; kk++)
	    printf(" %f ", G[ii+filter_width][jj+filter_width][kk+filter_width]) ;
	  printf("\n") ;
	}
      printf("\n") ;
    }



  /* beginning ... */  
  /* produce x-, y- and z- convolutions with gaussian filter */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	Density[k][i][j] = src[k][i][j] ;
  for( num=0; num<SmoothTimes; num++ )
    {
      printf("The %dth smoothing time ...\n", num) ;
      for( i=0; i<image_size; i++ )
	{
	  if(i%1==0) printf("i=%d\n", i) ;
	  
	  for( j=0; j<image_size; j++ )
	    for( k=0; k<z_size; k++ )
	      {
		if( src[k][i][j]>FocusIntensity )
		  {
		    total = 0 ;
		    TotalWeight = 0 ;
		    for(ii=-filter_width; ii<=filter_width; ii++)
		      for(jj=-filter_width; jj<=filter_width; jj++)
			for(kk=-filter_width; kk<=filter_width; kk++)
			  {
			    if(k+kk>=0 && k+kk<z_size && i+ii>=0 && i+ii<image_size && j+jj>=0 && j+jj<image_size )
			      {
				if( src[k+kk][i+ii][j+jj]>FocusIntensity ) 
				  {
				    TotalWeight += G[ii+filter_width][jj+filter_width][kk+filter_width] ;
				    total +=  G[ii+filter_width][jj+filter_width][kk+filter_width]*Density[k+kk][i+ii][j+jj] ; 
				  }
			      }
			  }
		    total /= TotalWeight ;
		    
		    outp[k][i][j] = total ;
		  }
		else
		  outp[k][i][j] = Density[k][i][j] ;
	      }
	}
      
      /* switching ... */
      for( i=0; i<image_size; i++ )
	for( j=0; j<image_size; j++ )
	  for( k=0; k<z_size; k++ )
	    Density[k][i][j] = outp[k][i][j] ;
    }
  

  /* total intensities normalization */
  src_nonZeroNum = 0 ;
  src_TotalIntensities = 0 ;
  Current_nonZeroNum = 0 ;
  Current_TotalIntensities = 0 ;
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  if( outp[k][i][j]>0 )
	    {
	      Current_nonZeroNum ++ ;
	      Current_TotalIntensities += outp[k][i][j] ;
	    }
	  if( src[k][i][j]>0 )
	    {
	      src_nonZeroNum ++ ;
	      src_TotalIntensities += src[k][i][j] ;
	    }
	}
   printf("original: average=%f, num=%f\n",     src_TotalIntensities/src_nonZeroNum,     src_nonZeroNum) ;
   printf("current:  average=%f, num=%f\n", Current_TotalIntensities/Current_nonZeroNum, Current_nonZeroNum) ;
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	outp[k][i][j] *= (src_TotalIntensities/Current_TotalIntensities) ;

  
  /* write the smoothed image */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  if( outp[k][i][j]>255. ) 
	    src[k][i][j] = 255 ;
	  else
	    src[k][i][j] = (unsigned char)(outp[k][i][j]+0.5) ;
	}
  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(src, z_size, image_size) ;
  Ffree3d(outp, z_size, image_size) ;
  Ffree3d(Density, z_size, image_size) ;
}


void Smoothing3DimageWithGuassian(int argc, char *argv[], float Smooth_Sigma)
{
  unsigned char ***src, ***outp ;
  int  i,j,k,l, ii,jj,kk, image_size, z_size ;
  FILE *fp;

  /* variables for Guassian filetering */
  int    filter_width;               /* length of 1-D gaussian mask */
  float a,b,c,d;        /* mask generation intermediate vars*/
  float ti,tj,tk,total,weights;
  float  r ;
  int efficient ;
  float ***G ;



  printf("Smooth_Sigma=%f\n", Smooth_Sigma) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp=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);

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

	 
  /* calc coeffs for 1-dimensional G, dG/dn and Delta-squared G filters */
  for(i=0; i<100; i++) 	
    {
      a=gaussian(i, Smooth_Sigma);
      if(a<0.005) {  filter_width=i; break; }
    }
  filter_width -- ;
  printf("filter_width=%d\n", filter_width) ;


  G = Falloc3d(filter_width*2+1, filter_width*2+1, filter_width*2+1); 
  for(ii=-filter_width; ii<=filter_width; ii++)
    for(jj=-filter_width; jj<=filter_width; jj++)
      for(kk=-filter_width; kk<=filter_width; kk++)
	{
	  r = sqrt(ii*ii+jj*jj+kk*kk) ;
	  
	  a=gaussian(r, Smooth_Sigma);
	  if(a>0.005 || r<2.)
	    G[ii+filter_width][jj+filter_width][kk+filter_width]=a/(6.283185*Smooth_Sigma*Smooth_Sigma);
	}

  /* normalize ... */
  weights = 0 ;
  for(ii=-filter_width; ii<=filter_width; ii++)
    for(jj=-filter_width; jj<=filter_width; jj++)
      for(kk=-filter_width; kk<=filter_width; kk++)
	weights += G[ii+filter_width][jj+filter_width][kk+filter_width] ; 

  for(ii=-filter_width; ii<=filter_width; ii++)
    for(jj=-filter_width; jj<=filter_width; jj++)
      for(kk=-filter_width; kk<=filter_width; kk++)
	G[ii+filter_width][jj+filter_width][kk+filter_width] /= weights ; 

  /* print */
  for(ii=-filter_width; ii<=filter_width; ii++)
    {
      for(jj=-filter_width; jj<=filter_width; jj++)
	{
	  for(kk=-filter_width; kk<=filter_width; kk++)
	    printf(" %f ", G[ii+filter_width][jj+filter_width][kk+filter_width]) ;
	  printf("\n") ;
	}
      printf("\n") ;
    }
  
  /* produce x-, y- and z- convolutions with gaussian filter */
  for( i=0; i<image_size; i++ )
    {
      if(i%1==0) printf("i=%d\n", i) ;

      for( j=0; j<image_size; j++ )
	for( k=0; k<z_size; k++ )
	  {
	    total = 0 ;
	    for(ii=-filter_width; ii<=filter_width; ii++)
	      for(jj=-filter_width; jj<=filter_width; jj++)
		for(kk=-filter_width; kk<=filter_width; kk++)
		  {
		    if(k+kk>=0 && k+kk<z_size && i+ii>=0 && i+ii<image_size && j+jj>=0 && j+jj<image_size )
		      total +=  G[ii+filter_width][jj+filter_width][kk+filter_width]*src[k+kk][i+ii][j+jj] ; 
		  }

	    if( total>255. ) total=255. ;
	    
	    outp[k][i][j] = (unsigned char)(total+0.5) ;	  
	  }
    }


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


  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(outp, z_size, image_size) ;
 }



void ThresholdingInput3DImage(int argc, char *argv[], int BinaryThreshold)
{
  unsigned char ***src, ***outp ;
  int  i,j,k, image_size, z_size ;
  FILE *fp;


  printf("\n\nBinarize 3D image. thresholding=%d\n\n", BinaryThreshold) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);

  /* thresholding */
  for(k=0;k<z_size;k++)
    {
      for(i=0;i<image_size;i++)
	for(j=0;j<image_size;j++)
	  if(src[k][i][j]>BinaryThreshold) outp[k][i][j]=255 ;
	  else                             outp[k][i][j]=0 ;
    }


  /* write the edge map */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);

  /* free */
  UCfree3d(outp, z_size, image_size) ;
  UCfree3d(src,  z_size, image_size) ;
}


double hypotenuse(double x,double y)
{
    if (x==0.0 && y==0.0) return(0.0);
    else return(hypot(x,y));
}


void canny_core(double s,int cols,int rows,unsigned char *data,unsigned char *derivative_mag,unsigned char *magnitude,unsigned char *orientation)
{
  int      filter_width;            /* length of 1-D gaussian mask */
  float    *gsmooth_x,*gsmooth_y;
  float    *derivative_x,*derivative_y;
  int      i,k,n;                   /* counters */
  int      t;                       /* temp. grad magnitude variable */
  double   a,b,c,d,g0;              /* mask generation intermediate vars*/
  double   ux,uy;
  double   t1,t2;
  double   G[20],dG[20],D2G[20];    /*gaussian & derivative filter masks*/
  double   gc,gn,gs,gw,ge,gnw,gne,gsw,gse;
  int      picsize,jstart,jlimit;
  int      ilimit;
  register int jfactor;
  int      kfactor1,kfactor2;
  int      kfactor;
  register int cindex,nindex,sindex,windex,eindex,nwindex,neindex,swindex,seindex;
  int      low=1,high=255;          /* tracker hysteresis parameters */
  int      mag_overflow_count=0;    /* used to measure how oft mag array overflows */
  int      mag_underflow_count=0;   /* used to measure how oft mag array underflows */

  double weightG, weightDG ;  /*normalizing ...*/



  picsize=cols*rows;    /* picture area */

  /* calc coeffs for 1-dimensional G, dG/dn and Delta-squared G filters */
  for(n=0; n<20; ++n)
    {
      a=gaussian((double)n,s);
      if (a>0.005 || n<2)
        {
	  b=gaussian((double)n-0.5,s);
	  c=gaussian((double)n+0.5,s);
	  d=gaussian((double)n,s*0.5);
	  /*fprintf(stderr,"a,b,c: %lf,%lf,%lf\n",a,b,c);*/
	  G[n]=(a+b+c)/3/(6.283185*s*s);
	  dG[n]=c-b;
	  D2G[n]=1.6*d-a; /* DOG */
	  /* fprintf(stderr,"G[%d]: %lf\n",n,G[n]);
	  fprintf(stderr,"dG[%d]: %lf\n",n,dG[n]); 
	  fprintf(stderr,"D2G[%d]: %lf\n",n,D2G[n]);*/	  
        }
      else break;
    }
  filter_width=n; /* printf("n=%d\n", n); */


  /*normalizing ...*/
  weightG  = 0.5*G[0] ; 
  weightDG = 0 ; 
  for(k=1;k<filter_width;++k)
    {
      weightG  +=  G[k] ; 
      weightDG += fabs(dG[k]) ;
    }
  for(k=0;k<filter_width;++k)
    {
      G[k]  /= weightG ; 
      dG[k] /= weightDG ; 
      /*fprintf(stderr,"G[%d]: %lf\n",k,G[k]);
	fprintf(stderr,"dG[%d]: %lf\n",k,dG[k]);*/
    }


  /*fprintf(stderr,"canny_core: smooth pic\n");*/
  /* allocate space for gaussian smoothing arrays */
  if ((gsmooth_x=(float *)calloc(picsize,sizeof(float)))==(float *)NULL)
    {
      fprintf(stderr,"can't alloc gsmooth_x\n");
      exit(0);
    }
  if ((gsmooth_y=(float *)calloc(picsize,sizeof(float)))==(float *)NULL)
    {
      fprintf(stderr,"can't alloc gsmooth_y\n");
      exit(0);
    }


  /* produce x- and y- convolutions with gaussian */
  ilimit=cols-(filter_width-1);
  jstart=cols*(filter_width-1);
  jlimit=cols*(rows-(filter_width-1));
  for (i=filter_width-1;i<ilimit;++i)
    {
      for(jfactor=jstart;
	  jfactor<jlimit;
	  jfactor+=cols)
        {
	  cindex=i+jfactor;
	  t1=data[cindex]*G[0];
	  t2=t1;
	  for(k=1,kfactor1=cindex-cols,
                kfactor2=cindex+cols;
	      k<filter_width;
	      k++,
                kfactor1-=cols,
                kfactor2+=cols)
            {
	      t1+=G[k]*(data[kfactor1]+
			data[kfactor2]);
	      t2+=G[k]*(data[cindex-k]+
			data[cindex+k]);
            }
	  gsmooth_x[cindex]=(float)t1;
	  gsmooth_y[cindex]=(float)t2;
        }
    }
    
  /* allocate space for gradient arrays */
  /*fprintf(stderr,"canny_core: find grad\n");*/


  if ((derivative_x=(float *)calloc(picsize,sizeof(float)))==(float *)NULL)
    {
      fprintf(stderr,"can't alloc x\n");
      exit(0);
    }
  /* produce x and y convolutions with derivative of gaussian */

  for (i=filter_width-1;i<ilimit;++i)
    {
      for(jfactor=jstart;
	  jfactor<jlimit;
	  jfactor+=cols)
        {
	  t1=0;
	  cindex=i+jfactor;
	  for (k=1;k<filter_width;++k)
	    t1+=dG[k]*(gsmooth_x[cindex-k]-
                       gsmooth_x[cindex+k]);
	  derivative_x[cindex] = (float)t1;
        }
    }
  free(gsmooth_x);
  if ((derivative_y=(float *)calloc(picsize,sizeof(float)))==(float *)NULL)
    {
      fprintf(stderr,"can't alloc y\n");
      exit(0);
    }
  
  for (i=n;i<cols-n;++i)
    {
      for(jfactor=jstart;jfactor<jlimit;jfactor+=cols)
        {
	  t2=0;
	  cindex=i+jfactor;
	  for (k=1,kfactor=cols;
	       k<filter_width;
	       k++,kfactor+=cols)
	    t2+=dG[k]*(gsmooth_y[cindex-kfactor]-gsmooth_y[cindex+kfactor]);
	  derivative_y[cindex] = (float)t2;
        }
    }
  free(gsmooth_y);
    

  /* non-maximum suppression (4 cases for orientation of line of max slope) */
  /* fprintf(stderr,"canny_core: non-maximum suppression\n"); */
  ilimit=cols-filter_width;
  jstart=cols*filter_width;
  jlimit=cols*(rows-filter_width);

  for (i=filter_width;i<ilimit;++i)
    {
      for (jfactor=jstart;
	   jfactor<jlimit;
	   jfactor+=cols)
        {
				/* calculate current indeces */
	  cindex=i+jfactor;
	  nindex=cindex-cols;
	  sindex=cindex+cols;
	  windex=cindex-1;
	  eindex=cindex+1;
	  nwindex=nindex-1;
	  neindex=nindex+1;
	  swindex=sindex-1;
	  seindex=sindex+1;
	  ux=derivative_x[cindex];
	  uy=derivative_y[cindex];
	  gc=hypotenuse(ux,uy);
	  /* scale gc to fit into an unsigned char array */
	  t=(int)(gc*1.5); /* 3.0 July 26, 1999 */
	  /*fprintf(stderr,"canny_core: i,j=(%d,%d), t=%lf\n",i,jfactor/cols,t);*/
	  derivative_mag[cindex]=(t<256 ? t : 255);
	  gn=hypotenuse(derivative_x[nindex],derivative_y[nindex]);
	  gs=hypotenuse(derivative_x[sindex],derivative_y[sindex]);
	  gw=hypotenuse(derivative_x[windex],derivative_y[windex]);
	  ge=hypotenuse(derivative_x[eindex],derivative_y[eindex]);
	  gne=hypotenuse(derivative_x[neindex],derivative_y[neindex]);
	  gse=hypotenuse(derivative_x[seindex],derivative_y[seindex]);
	  gsw=hypotenuse(derivative_x[swindex],derivative_y[swindex]);
	  gnw=hypotenuse(derivative_x[nwindex],derivative_y[nwindex]);
	  if (ux*uy>0)
            {
	      if(fabs(ux)<fabs(uy))
                {
		  if((g0=fabs(uy*gc))
                     < fabs(ux*gse+(uy-ux)*gs) ||
                     g0<=fabs(ux*gnw+(uy-ux)*gn))
                    continue;
                }
	      else
                {
		  if((g0=fabs(ux*gc))
                     < fabs(uy*gse+(ux-uy)*ge) ||
                     g0<=fabs(uy*gnw+(ux-uy)*gw))
                    continue;
                }
            }
	  else
            {
	      if(fabs(ux)<fabs(uy))
                {
		  if((g0=fabs(uy*gc))
                     < fabs(ux*gne-(uy+ux)*gn) ||
                     g0<=fabs(ux*gsw-(uy+ux)*gs))
		    continue;
                }
	      else
                {
		  if((g0=fabs(ux*gc))
                     < fabs(uy*gne-(ux+uy)*ge) ||
                     g0<=fabs(uy*gsw-(ux+uy)*gw))
		    continue;
                }
            }
	  /* seems to be a good scale factor */
	  magnitude[cindex]=derivative_mag[cindex];
	  /* pi*40 ~= 128 - direction is (thought of as) a signed byte */
	  orientation[cindex]=(unsigned char)((PI+atan2(ux, uy))*ORIENT_SCALE);
/*		if (orientation[cindex]>(3.1416*ORIENT_SCALE)) 
			orientation[cindex]=orientation[cindex];
		if (orientation[cindex]<0)
			orientation[cindex]=orientation[cindex];
	  float t1,t2,t3,t4;
	  unsigned char ut1,ut2,ut3,ut4;

	  t1=atan2(1,2);
	  t2=atan2(-1,2);
	  t3=atan2(1,-2);
	  t4=atan2(-1,-2);
	  	  ut1=(unsigned char)(t1*ORIENT_SCALE);	 
	  ut2=(unsigned char)(t2*ORIENT_SCALE);	
	  ut3=(unsigned char)(t3*ORIENT_SCALE);	
	  ut4=(unsigned char)(t4*ORIENT_SCALE);	 
*/ 
	   }
    } 

  free(derivative_x);
  free(derivative_y);
}


void Cubic_edgeDetectionBy2DcannyMethod(int argc, char *argv[], double CannyGuassianSigma)
{
  /* Canny */
  unsigned char *XY_img, *XY_magnitude, *XY_orient, *XY_derivative_mag; /* mag of del G before non-maximum suppression */
  int  low=10, high=16 ; /* for track */

  /* 3D */
  unsigned char ***src, ***outp ;
  int  i,j,k, image_size, z_size ;
  FILE *fp;
  int ii,jj,kk ;
  int GreyValue ;


  printf("\n\n2D Canny edge detector for 3D application, with Guassian sigma=%f\n\n", CannyGuassianSigma) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);


  /* Canny part: memory applciation*/
  XY_img           = (unsigned char *)malloc(image_size*image_size*sizeof(unsigned char));
  XY_magnitude     = (unsigned char *)calloc(image_size*image_size,1) ;
  XY_orient        = (unsigned char *)calloc(image_size*image_size,1) ;
  XY_derivative_mag= (unsigned char *)calloc(image_size*image_size,1) ;


  /* edge detection*/
   for(k=0;k<z_size;k++)
   {
     if(k%10==0)
       printf("%dth slice\n", k) ;

     for(i=0;i<image_size;i++)
       for(j=0;j<image_size;j++)
	 XY_img[i*image_size+j] = src[k][i][j] ;

     /* Canny */
     canny_core(CannyGuassianSigma, image_size, image_size, XY_img, XY_derivative_mag, XY_magnitude, XY_orient);

     for(i=0;i<image_size;i++)
       for(j=0;j<image_size;j++)
	 {
	   outp[k][i][j] = XY_magnitude[i*image_size+j] ;
	   XY_magnitude[i*image_size+j] = 0 ; /* reset */
	 }
   }


  /* write the edge map */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);

 
  /* 2D: free */
   free(XY_img) ;
   free(XY_orient) ;
   free(XY_derivative_mag) ;
   free(XY_magnitude) ;
  /* 3D: free */
   UCfree3d(outp, z_size, image_size) ;
   UCfree3d(src,  z_size, image_size) ;
}



void Cubic_edgeDetectionBythree2DcannyMethod(int argc, char *argv[], double CannyGuassianSigma)
{
  /* Canny */
  unsigned char *XY_img, *XY_magnitude, *XY_orient, *XY_derivative_mag; /* mag of del G before non-maximum suppression */
  int  low=10, high=16 ; /* for track */

  /* 3D */
  unsigned char ***src, ***outp ;
  int  i,j,k, image_size, z_size ;
  FILE *fp;
  int ii,jj,kk ;
  int GreyValue ;


  printf("\n\nthree 2D Canny edge detector for 3D application, with Guassian sigma=%f\n\n", CannyGuassianSigma) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);


  /* Z dimension: */
   printf("\n\nZ dimension\n") ;
   /* Canny part--memory applciation*/
    XY_img           = (unsigned char*)malloc(image_size*image_size*sizeof(unsigned char));
    XY_magnitude     = (unsigned char *)calloc(image_size*image_size,1) ;
    XY_orient        = (unsigned char *)calloc(image_size*image_size,1) ;
    XY_derivative_mag= (unsigned char *)calloc(image_size*image_size,1) ;

    /* edge detection*/
    for(k=0;k<z_size;k++)
      {
	if(k%10==0)
	  printf("%dth slice\n", k) ;
	
	for(i=0;i<image_size;i++)
	  for(j=0;j<image_size;j++)
	    XY_img[i*image_size+j] = src[k][i][j] ;
	
	/* Canny */
	canny_core(CannyGuassianSigma, image_size, image_size, XY_img, XY_derivative_mag, XY_magnitude, XY_orient);
	
	for(i=0;i<image_size;i++)
	  for(j=0;j<image_size;j++)
	    {
	      outp[k][i][j] = XY_magnitude[i*image_size+j] ;
	      XY_magnitude[i*image_size+j] = 0 ; /* reset */
	    }
      }
    /* 2D: free */
     free(XY_img) ;
     free(XY_orient) ;
     free(XY_derivative_mag) ;
     free(XY_magnitude) ;

  /* X dimension: */
   printf("\n\nX dimension\n") ;
   /* Canny part--memory applciation*/
    XY_img           = (unsigned char *)malloc(z_size*image_size*sizeof(unsigned char));
    XY_magnitude     = (unsigned char *)calloc(z_size*image_size,1) ;
    XY_orient        = (unsigned char *)calloc(z_size*image_size,1) ;
    XY_derivative_mag= (unsigned char *)calloc(z_size*image_size,1) ;

    /* edge detection*/
    for(i=0;i<image_size;i++)
      {
	if(i%10==0)
	  printf("%dth slice\n", i) ;
	
	for(k=0;k<z_size;k++)
	  for(j=0;j<image_size;j++)
	    XY_img[k*image_size+j] = src[k][i][j] ;
	
	/* Canny */
	canny_core(CannyGuassianSigma, image_size, z_size, XY_img, XY_derivative_mag, XY_magnitude, XY_orient);
	
	for(k=0;k<z_size;k++)
	  for(j=0;j<image_size;j++)
	    {
	      if( outp[k][i][j]<XY_magnitude[k*image_size+j] )
		outp[k][i][j] = XY_magnitude[k*image_size+j] ;
	      XY_magnitude[k*image_size+j] = 0 ; /* reset */
	    }
      }
    /* 2D: free */
     free(XY_img) ;
     free(XY_orient) ;
     free(XY_derivative_mag) ;
     free(XY_magnitude) ;

  /* Y dimension: */
   printf("\n\nY dimension\n") ;
   /* Canny part--memory applciation*/
    XY_img           = (unsigned char *)malloc(z_size*image_size*sizeof(unsigned char));
    XY_magnitude     = (unsigned char *)calloc(z_size*image_size,1) ;
    XY_orient        = (unsigned char *)calloc(z_size*image_size,1) ;
    XY_derivative_mag= (unsigned char *)calloc(z_size*image_size,1) ;

    /* edge detection*/
    for(j=0;j<image_size;j++)
      {
	if(j%10==0)
	  printf("%dth slice\n", j) ;
	
	for(k=0;k<z_size;k++)
	  for(i=0;i<image_size;i++)
	    XY_img[k*image_size+i] = src[k][i][j] ;
	
	/* Canny */
	canny_core(CannyGuassianSigma, image_size, z_size, XY_img, XY_derivative_mag, XY_magnitude, XY_orient);
	
	for(k=0;k<z_size;k++)
	  for(i=0;i<image_size;i++)
	    {
	      	      if( outp[k][i][j]<XY_magnitude[k*image_size+i] )
		outp[k][i][j] = XY_magnitude[k*image_size+i] ;
	      XY_magnitude[k*image_size+i] = 0 ; /* reset */
	    }
      }
    /* 2D: free */
     free(XY_img) ;
     free(XY_orient) ;
     free(XY_derivative_mag) ;
     free(XY_magnitude) ;



  /* write the edge map */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);

 
  /* 3D: free */
   UCfree3d(outp, z_size, image_size) ;
   UCfree3d(src,  z_size, image_size) ;
}




void Reflect3DImageByHalfY(int argc, char *argv[], int fold)
{
  unsigned char ***src, ***outp ;
  int  i,j,k, image_size, z_size ;
  FILE *fp;
  int ii,jj,kk ;
  int GreyValue ;


  printf("Reflect 3D image by y=image_size/2.\n") ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);


  /* relfect image by y=image_size/2. */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  /* refelction point */
	  ii = i ;
	  jj = 2*fold-j; /*jj = (image_size-1) - j ;*/ /*Dec 18 1999*/
	  kk = k ;

	  if(jj>=0 && jj<image_size) 
	    GreyValue = (src[k][i][j]+src[kk][ii][jj])/2 ;
	  else 
	    GreyValue = src[k][i][j] ;


	  outp[k][i][j] = GreyValue ;
	}


  /* write the smoothed image */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);


  /* free */
  UCfree3d(outp, z_size, image_size) ;
  UCfree3d(src,  z_size, image_size) ;
}


void Dig3DpartANDscaleIt(int argc, char *argv[], float zoomm)
{
  unsigned char ***src, ***outp ;
  int  i,j,k, image_size, z_size ;
  FILE *fp;
  Fvector3d center,sum;
  float total, ii,jj,kk, SumGrey ;
  int ni,nj,nk ;
  float x, y, z ;
  char choice ;

  x=135. ;
  y=122. ;
  z=60. ;

  printf("Will extract a 3D part cented at (%f, %f, %f),\n", x,y,z) ;
  printf("and then zoomed by %f.\n", zoomm) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);

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

  printf("Input the center manually or automatically detect the image center (y/n):") ;
  scanf("%c", &choice) ;

  if(choice=='y')
    {
      printf("Input the center manually, i.e. x=%f y=%f z=%f\n", x,y,z) ;
      scanf("%f%f%f", &x, &y, &z) ;
 
      /* center */
      center.x=x;
      center.y=y;
      center.z=z;

      printf("center.x=%f center.y=%f center.z=%f\n",center.x,center.y,center.z);
    }
  else
    {
      /* center of the input object */
      sum.x=sum.y=sum.z=total=0.0;
      SumGrey = 0 ;
      
      for(k=0;k<z_size;k++)
	for(i=0;i<image_size;i++)
	  for(j=0;j<image_size;j++)
	    {
	      if(src[k][i][j]>10)
		{
		  sum.x+=((float) (i));
		  sum.y+=((float) (j));
		  sum.z+=((float) (k));
		  total+=1.0;

		  SumGrey += src[k][i][j] ;
		}
	    }
      center.x=sum.x/total;
      center.y=sum.y/total;
      center.z=sum.z/total;
      SumGrey /= (total*255.) ;
      printf("center.x=%f center.y=%f center.z=%f   ratio=%f\n",center.x,center.y,center.z, SumGrey);
    }

  /* scaling by linear interpolation ... */
printf("\nscaling in z direction:%f\n", Z_THICKINK) ;
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  ii = (i-image_size/2)/zoomm + center.x ;
	  jj = (j-image_size/2)/zoomm + center.y;
	  kk = (k-z_size/2)/zoomm/Z_THICKINK + center.z ; /*1.*/

	  ni = (int)ii ;
	  nj = (int)jj ;
	  nk = (int)kk ; 

      /*

        In the original code which has been copied from the HAMMER copy that was previously included
        with the DTIDROID repository, the tri-linear interpolation never took effect as the
        outp[k]i][j] was always set to src[nk][ni][nj]. Therefore, comment this block, which is
        a slightly (only syntactically) modified version of the original.

	  niP1 = ni+1 ;
	  njP1 = nj+1 ;
	  nkP1 = nk+1 ;

      if(nk>=0 && nk<z_size-1 && ni>=0 && ni<image_size-1 && nj>=0 && nj<image_size-1)
	    {
	      a = kk-nk ;        a1 = 1.-a ;
	      b = ii-ni ;        b1 = 1.-b ;
	      c = jj-nj ;        c1 = 1.-c ;
	      
	      outp[k][i][j] = static_cast<unsigned char>((src[nk][ni][nj]*(a1*b1*c1) + src[nk][niP1][nj]*(a1*b*c1) + src[nk][ni][njP1]*(a1*b1*c) + src[nk][niP1][njP1]*(a1*b*c) + src[nkP1][ni][nj]*(a*b1*c1) + src[nkP1][niP1][nj]*(a*b*c1) + src[nkP1][ni][njP1]*(a*b1*c) + src[nkP1][niP1][njP1]*(a*b*c))/((a1*b1*c1)+(a1*b*c1)+(a1*b1*c)+(a1*b*c)+(a*b1*c1)+(a*b*c1)+(a*b1*c)+(a*b*c))) ; 
	    }
      */

      outp[k][i][j] = src[nk][ni][nj];
	}


  /* write the smoothed image */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);


  /* free */
  UCfree3d(outp, z_size, image_size) ;
  UCfree3d(src,  z_size, image_size) ;
}



void Scale3DImageByLinearInterpolation(int argc, char *argv[], float enlarge)
{
  unsigned char ***src, ***outp ;
  int  i,j,k, image_size, z_size ;
  FILE *fp;
  Fvector3d center,sum;
  float total ;

  /* for linear interpolation */
  float ii,jj,kk, b,c,d, b1,c1,d1;
  int   ni,nj,nk, niP1,njP1,nkP1 ;


  printf("enlarge rate=%f\n", enlarge) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);

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

/* center of the input object */
  sum.x=sum.y=sum.z=total=0.0;

  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if(src[k][i][j]>10)
	    {
	      sum.x+=((float) (i));
	      sum.y+=((float) (j));
	      sum.z+=((float) (k));
	      total+=1.0;
	    }
	}
  center.x=sum.x/total;
  center.y=sum.y/total;
  center.z=sum.z/total;

  printf("center.x=%f center.y=%f center.z=%f\n",center.x,center.y,center.z);

  printf("Z_THICKINK=%f", Z_THICKINK) ;

  /* scaling ... */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  ii = (i-image_size/2)/enlarge + center.x ;
	  jj = (j-image_size/2)/enlarge + center.y;
	  kk = (k-z_size/2)/enlarge/Z_THICKINK + center.z ; /*1.*/

	  ni = (int)ii ;
	  nj = (int)jj ;
	  nk = (int)kk ;

	  niP1 = ni+1 ;
	  njP1 = nj+1 ;
	  nkP1 = nk+1 ;

	  if(ni>=0 && ni<image_size-1  &&  nj>=0 && nj<image_size-1  &&  nk>=0 && nk<z_size-1 )
	    {
	      b = ii-ni ;        b1 = 1.-b ;
	      c = jj-nj ;        c1 = 1.-c ;
	      d = kk-nk ;        d1 = 1.-d ;

	      outp[k][i][j] = static_cast<unsigned char>(
                  (d1 * (src[nk  ][ni][nj]*(b1*c1) + src[nk  ][niP1][nj]*(b*c1) + src[nk  ][ni][njP1]*(b1*c) + src[nk  ][niP1][njP1]*(b*c))
                 + d  * (src[nkP1][ni][nj]*(b1*c1) + src[nkP1][niP1][nj]*(b*c1) + src[nkP1][ni][njP1]*(b1*c) + src[nkP1][niP1][njP1]*(b*c)))
                 / (d1*((b1*c1)+(b*c1)+(b1*c)+(b*c)) + d*((b1*c1)+(b*c1)+(b1*c)+(b*c)))) ;
	    }
      else if((ni==image_size-1 && nj>=0 && nj<image_size-1 && nk>=0 && nk<z_size-1) || (ni>=0 && ni<image_size-1 && nj==image_size-1 && nk>=0 && nk<z_size-1) || (ni>=0 && ni<image_size-1 && nj>=0 && nj<image_size-1 && nk==z_size-1))
        {
	      outp[k][i][j] = src[nk][ni][nj] ; 
        }
      else
        {
          outp[k][i][j] = 0 ;
        }
	}


  /* write the smoothed image */
  printf("Write sampled file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);


  /* free */
  UCfree3d(outp, z_size, image_size) ;
  UCfree3d(src,  z_size, image_size) ;
}



void Scale3DsegmentedImage(int argc, char *argv[], float enlarge)
{
  unsigned char ***src, ***outp ;
  int  i,j,k,ii,jj,kk, image_size, z_size ;
  FILE *fp;
  Fvector3d center,sum;
  float total ;


  printf("enlarge rate=%f, Z_THICKINK=%f\n", enlarge, Z_THICKINK) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  outp = 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);

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

/* center of the input object */
  sum.x=sum.y=sum.z=total=0.0;

  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if(src[k][i][j]>10)
	    {
	      sum.x+=((float) (i));
	      sum.y+=((float) (j));
	      sum.z+=((float) (k));
	      total+=1.0;
	    }
	}
  center.x=sum.x/total;
  center.y=sum.y/total;
  center.z=sum.z/total;
  /*
    center.x = 134.624252 ;
    center.y = 126.082588 ;
    center.z = 68.208946 ;
  */
/* Added for SYYU

      printf("Input the center manually, i.e. x=%f y=%f z=%f\n", center.x,center.y,center.z) ;
      scanf("%f%f%f", &center.x, &center.y, &center.z) ;
*/

  printf("center.x=%f center.y=%f center.z=%f\n",center.x,center.y,center.z);


  /* scaling ... */
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  ii = static_cast<int>((i - image_size / 2) / enlarge              + center.x);
	  jj = static_cast<int>((j - image_size / 2) / enlarge              + center.y);
	  kk = static_cast<int>((k - z_size     / 2) / enlarge / Z_THICKINK + center.z);

	  if(kk>=0 && kk<z_size && ii>=0 && ii<image_size && jj>=0 && jj<image_size) {
	    outp[k][i][j] = src[kk][ii][jj] ;
      }
	}


  /* write the smoothed image */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(outp[k][i],1,image_size,fp);
  fclose(fp);


  /* free */
  UCfree3d(outp, z_size, image_size) ;
  UCfree3d(src,  z_size, image_size) ;
}

void Denoising3DedgeImage(int argc, char *argv[], int winWID)
{
  unsigned char ***Edge, ***outEdge ;
  int  i,j,k,l,m,n, image_size, z_size, number, cube, Greythreshold ;
  FILE *fp;


  printf("Window Width=%d\n", winWID) ;


  /* image size in xy plane */
  image_size=imageSz;

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

  /* apply for memory */ 
  Edge    = UCalloc3d(image_size,image_size,z_size);
  outEdge = UCalloc3d(image_size,image_size,z_size);

  /* read data */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fread(Edge[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++ )
	outEdge[k][i][j] = 0 ;
	

  /* filtering ... */
  cube = (2*winWID+1)*(2*winWID+1)*winWID ;
  Greythreshold = 100 ;
  for( i=winWID; i<image_size-winWID; i++ )
    for( j=winWID; j<image_size-winWID; j++ )
      for( k=winWID; k<z_size-winWID; k++ )
	{
	  if(Edge[k][i][j]>Greythreshold)
	    {
	      number = 0 ;
	      for( l=-winWID; l<=winWID; l++ )
		for( m=-winWID; m<=winWID; m++ )
		  for( n=-winWID; n<=winWID; n++ )
		    if( Edge[k+l][i+m][j+n]> Greythreshold) number++ ;
	      
	      if( number>cube )  outEdge[k][i][j] = Edge[k][i][j];
	    }
	}

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


  /* free */
  UCfree3d(outEdge, z_size, image_size) ;
  UCfree3d(Edge, z_size, image_size) ;
}




void CannyLikedEdgeDetectionWithGuassianSmoothing(int argc, char *argv[], float Sigma)
{
  unsigned char ***src, ***Edge ;
  int  i,j,k,l, image_size, z_size ;
  FILE *fp;

  /* variables for Guassian filetering */
  int    filter_width;               /* length of 1-D gaussian mask */
  float  ***gsmooth_x, ***gsmooth_y, ***gsmooth_z;
  float  ***derivative_x, ***derivative_y, ***derivative_z;
  int    n;            /* counters */
  float a,b,c,d;        /* mask generation intermediate vars*/
  float ti,tj,tk,total,weights;
  float G[20],dG[20],D2G[20];    /*gaussian & derivative filter masks*/
  int efficient ;


  printf("Sigma=%f\n", Sigma) ;


  /* image size in xy plane */
  image_size=imageSz;

  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);
  Edge=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);

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

	 
  /* calc coeffs for 1-dimensional G, dG/dn and Delta-squared G filters */
  for(n=0; n<20; n++)
    {
      a=gaussian((float)n, Sigma);
      if(a>0.005 || n<2)
	{
	  b=gaussian((float)n-0.5, Sigma);
	  c=gaussian((float)n+0.5, Sigma);
	  d=gaussian((float)n, Sigma*0.5);
	  fprintf(stderr,"a,b,c: %lf,%lf,%lf\n",a,b,c);
	  G[n]=(a+b+c)/3/(6.283185*Sigma*Sigma);
	  dG[n]=c-b;
	  D2G[n]=1.6*d-a; /* DOG */
	}
      else break;
    }
  filter_width=n; printf("n=%d\n", n);

for(i=0; i<filter_width; i++)
printf("dG[%d]=%f\n", i, dG[i]) ;

  /* normalize ... */
  weights = G[0] ;
  for(k=1;k<filter_width;++k)
    weights += G[k]*2.0 ; 
  for(k=0;k<filter_width;++k)
    G[k] /= weights ; 

  
  fprintf(stderr,"smoothing pic\n");
  /* allocate space for gaussian smoothing arrays */
  gsmooth_x = Falloc3d(image_size,image_size,z_size);
  gsmooth_y = Falloc3d(image_size,image_size,z_size);
  gsmooth_z = Falloc3d(image_size,image_size,z_size);
  printf("succeeded in memory application!\n") ;

  /* produce x-, y- and z- convolutions with gaussian filter */
  for( i=filter_width-1; i<image_size-(filter_width-1); i++ )
    {
      if(i%30==0) printf("i=%d\n", i) ;

      for( j=filter_width-1; j<image_size-(filter_width-1); j++ )
	for( k=filter_width-1; k<z_size-(filter_width-1); k++ )
	  {
	    ti = src[k][i][j]*G[0];
	    tj = ti;
	    tk = ti ;
	    for(l=1; l<filter_width; l++)
	      {
		tk += G[l]*(src[k+l][i][j]+src[k-l][i][j]);
		ti += G[l]*(src[k][i+l][j]+src[k][i-l][j]);
		tj += G[l]*(src[k][i][j+l]+src[k][i][j-l]);
	      }
	    gsmooth_x[k][i][j] = ti ;
	    gsmooth_y[k][i][j] = tj ;
	    gsmooth_z[k][i][j] = tk ;
	  }
    }


  /* derivative */
  /* x: */
  derivative_x = Falloc3d(image_size,image_size,z_size);
  for( i=filter_width-1; i<image_size-(filter_width-1); i++ )
    {
      if(i%30==0) printf("i=%d\n", i) ;
      
      for( j=filter_width-1; j<image_size-(filter_width-1); j++ )
	for( k=filter_width-1; k<z_size-(filter_width-1); k++ )
	  {
	    ti=0;
	    for(l=1;l<filter_width;l++)
	      ti += dG[l]*(gsmooth_x[k][i-l][j]-
			   gsmooth_x[k][i+l][j]);
	    derivative_x[k][i][j] = ti;
	  }
    }
  Ffree3d(gsmooth_x, z_size, image_size) ;
  printf("pass x\n") ;

  /* y: */
  derivative_y = Falloc3d(image_size,image_size,z_size);
  for( i=filter_width-1; i<image_size-(filter_width-1); i++ )
    {
      if(i%30==0) printf("i=%d\n", i) ;
      for( j=filter_width-1; j<image_size-(filter_width-1); j++ )
	for( k=filter_width-1; k<z_size-(filter_width-1); k++ )
	  {
	    tj=0;
	    for(l=1;l<filter_width;l++)
	      tj += dG[l]*(gsmooth_y[k][i][j-l]-
			   gsmooth_y[k][i][j+l]);
	    derivative_y[k][i][j] = tj;
	  }
    }
  Ffree3d(gsmooth_y, z_size, image_size) ;
  printf("pass y\n") ;

  /*z:*/
  derivative_z = Falloc3d(image_size,image_size,z_size);
  for( i=filter_width-1; i<image_size-(filter_width-1); i++ )
    {
      if(i%30==0) printf("i=%d\n", i) ;
      for( j=filter_width-1; j<image_size-(filter_width-1); j++ )
	for( k=filter_width-1; k<z_size-(filter_width-1); k++ )
	  {
	    tk=0;
	    for(l=1;l<filter_width;l++)
	      tk += dG[l]*(gsmooth_z[k-l][i][j]-
			   gsmooth_z[k+l][i][j]);
	    derivative_z[k][i][j] = tk;
	  }
    }
  Ffree3d(gsmooth_z, z_size, image_size) ;
  printf("pass z\n") ;


  /* merge ... */       
  printf("merging ...\n") ;
  efficient = 2 ;
  for( i=(filter_width-1)*efficient; i<image_size-(filter_width-1)*efficient; i++ )
    {
      if(i%30==0) printf("i=%d\n", i) ;
      
      for( j=(filter_width-1)*efficient; j<image_size-(filter_width-1)*efficient; j++ )
	for( k=(filter_width-1)*efficient; k<z_size-(filter_width-1)*efficient; k++ )
	  {
	    total = sqrt((derivative_x[k][i][j]*derivative_x[k][i][j]+
			  derivative_y[k][i][j]*derivative_y[k][i][j]+
			  derivative_z[k][i][j]*derivative_z[k][i][j])/1.) ;
	    if( total>255. ) total=255. ;
	    
	    if(total<100.0) total=total ; /*0*/
	    Edge[k][i][j] = (unsigned char)total ;
	  }
    }

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


  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(Edge, z_size, image_size) ;
  Ffree3d(derivative_x, z_size, image_size) ;
  Ffree3d(derivative_y, z_size, image_size) ;
  Ffree3d(derivative_z, z_size, image_size) ;
  /*  Ffree3d(gsmooth_x, z_size, image_size) ;
  Ffree3d(gsmooth_y, z_size, image_size) ; 
  Ffree3d(gsmooth_z, z_size, image_size) ;*/
 }

void ExtractInterestingPart(int argc, char *argv[], int threshold)
{
  unsigned char ***src,***seg;
  int i,j,k, image_size,z_size ;
  FILE *fp;


  /* image size in xy plane */
  image_size=imageSz;

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


  fp=myopen(argv[1],"r");  printf("%s\n",argv[1]) ;
  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);
  seg=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);

  
  /* extraction */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if(src[k][i][j]==threshold)
	    seg[k][i][j] = 255 ;
	  else 
	    seg[k][i][j] = 0 ;
	}

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

  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(seg, z_size, image_size) ;
}


void ExtractInterestingPartWithVariation(int argc, char *argv[], int threshold, int variation)
{
  unsigned char ***src,***seg;
  int i,j,k, image_size,z_size ;
  FILE *fp;


  /* image size in xy plane */
  image_size=imageSz;

  printf("threshold=%d\nvariation=%d\n\n\n", threshold, variation) ;


  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);
  seg=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);

  
  /* extraction */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if(src[k][i][j]>=threshold-variation && src[k][i][j]<=threshold+variation)
	    seg[k][i][j] = 255 ;
	  else 
	    seg[k][i][j] = 0 ;
	}

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

  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(seg, z_size, image_size) ;
}


void GetEdgeFrom3DimageWithoutSmoothing(int argc, char *argv[], int threshold)
{
  unsigned char ***src,***seg;
  int i,j,k, image_size,z_size, temp, tempk, tempi, tempj, l,m ;
  FILE *fp;


  /* image size in xy plane */
  image_size=imageSz;


  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);
  seg=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);

  
  /* edge detection */ 
  for(k=1;k<z_size-1;k++)
    for(i=1;i<image_size-1;i++)
      for(j=1;j<image_size-1;j++)
	{
	  tempk = 0 ;
	  for(l=-1; l<=1; l++)
	    for(m=-1; m<=1; m++)
	      {
	      tempk += abs(src[k+1][i+l][j+m]  - src[k-0][i+l][j+m]) ;
	      /*   if(tempk!=0) 
		printf("k=%d i=%d j=%d    %d\n", k,i,j, tempk) ;	*/       
	      }
	  tempk /= 9 ;

	  tempi = 0 ;
	  for(l=-1; l<=1; l++)
	    for(m=-1; m<=1; m++)
	      tempi += abs(src[k+l][i+1][j+m]  - src[k+l][i-0][j+m]) ;
	  tempi /= 9 ;

	  tempj = 0 ;
	  for(l=-1; l<=1; l++)
	    for(m=-1; m<=1; m++)
	      tempj += abs(src[k+l][i+m][j+1]  - src[k+l][i+m][j-0]) ;
	  tempj /= 9 ;


	  if( tempk>tempi && tempk>tempj )  temp = tempk ;
	  if( tempi>tempk && tempi>tempj )  temp = tempi ;
	  if( tempj>tempk && tempj>tempi )  temp = tempj ;
	  if( tempj==tempk==tempj==0 )  temp = tempj ;

	  if( (tempi+tempj+tempk)/2>temp ) temp=(tempi+tempj+tempk)/2 ;

	  if(temp>threshold)
	    seg[k][i][j] = 255 ;
	  else 
	    seg[k][i][j] = 0 ;
	}

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

  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(seg, z_size, image_size) ;
}


void BlackReversedToWhiteInImg(int argc, char *argv[])
{
  unsigned char ***src;
  int  i,j,k, image_size,z_size, temp ;
  float  ns ;
  FILE *fp;
  long idum ;


  /* image size in xy plane */
  image_size=imageSz;


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

  
  /* half */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	if( src[k][i][j]<10 )  src[k][i][j] = 255 ;

  /* write the result */
  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(src,     z_size, image_size) ;
}


void HalfGreyLevelOfImg(int argc, char *argv[])
{
  unsigned char ***src;
  int  i,j,k, image_size,z_size, temp ;
  float  ns ;
  FILE *fp;
  long idum ;


  /* image size in xy plane */
  image_size=imageSz;


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

  
  /* half */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	src[k][i][j] /= 2 ;

  /* write the result */
  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(src,     z_size, image_size) ;
}



void DoubleGreyLevelOfImg(int argc, char *argv[])
{
  unsigned char ***src;
  int  i,j,k, image_size,z_size, temp ;
  float  ns ;
  FILE *fp;
  long idum, intensity ;


  /* image size in xy plane */
  image_size=imageSz;


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

  
  /* half */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  intensity = src[k][i][j]*2 ;

	  if( intensity>255 )
	    src[k][i][j] = 255 ;
	  else
	    src[k][i][j] = intensity ;
	}

  /* write the result */
  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(src,     z_size, image_size) ;
}



void AddNoiseTo3Dimage(int argc, char *argv[], float Variance)
{
  unsigned char ***src,***corrupt;
  int  i,j,k, image_size,z_size, temp ;
  float  ns ;
  FILE *fp;
  long idum ;


  /* image size in xy plane */
  image_size=imageSz;


  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);
  corrupt=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);

  
  /* edge detection */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  temp = (int)(src[k][i][j]) + (int)(gasdev(&idum)*Variance) ;

	  if(temp>255) temp = 255;
	  if(temp<0)   temp = 0;
	  corrupt[k][i][j] = (unsigned char) temp ;
	}

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

  /* free */
  UCfree3d(src,     z_size, image_size) ;
  UCfree3d(corrupt, z_size, image_size) ;
}


void MappingIntensitiesOfImage(int argc,char *argv[])
{
  FILE          *fp;
  extern char   *optarg;
  char          filename[80] ;
  unsigned char ***Orig, ***Mapped;
  int           num, c,i,j,k ;
  int           gm, wm, vn ;
  int           image_size,z_size;

  /* open images */
   /* original image */
   image_size = imageSz ;
   
   fp=myopen(argv[1],"r"); printf("%s\n", filename) ;
   fseek(fp,0,SEEK_END);
   z_size=ftell(fp)/(image_size*image_size);
   rewind(fp);
   
   Orig   = UCalloc3d(image_size,image_size,z_size);
   Mapped = UCalloc3d(image_size,image_size,z_size);
   
   for(k=0;k<z_size;k++)
     for(i=0;i<image_size;i++)
       fread(Orig[k][i],1,image_size,fp);
   fclose(fp);

  /* setting averages */
  gm = 82 ;
  wm = 112;
  vn = 35 ;
  for( i=0; i<image_size; i++ )
    for( j=0; j<image_size; j++ )
      for( k=0; k<z_size; k++ )
	{
	  if( Orig[k][i][j]>0 &&  Orig[k][i][j]<=90 )
	    Mapped[k][i][j]=vn ;
	  if( Orig[k][i][j]>90 &&  Orig[k][i][j]<=150 )
	    Mapped[k][i][j]=gm ;
	  if( Orig[k][i][j]>150)
	    Mapped[k][i][j]=wm ;
	}

  /* write image */
  printf("Write file!\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(Mapped[k][i],1,image_size,fp);
  fclose(fp);
}

void  CropPartOfImgAroundCenterOfImg(int argc,char *argv[], int sX, int sY, int sZ)
{
  unsigned char ***src,***seg;
  int i,j,k, image_size,z_size,x,y,z ;
  FILE *fp;
  Fvector3d center,sum;
  float total ;
  int  z_min, z_max ;
  int  x_min, x_max ;
  int  y_min, y_max ;

  /* image size in xy plane */
  image_size=imageSz;

  printf("sX=%d sY=%d sZ=%d\n\n\n", sX, sY, sZ) ;


  fp=myopen(argv[1],"r"); printf("%s\n",argv[1]) ;
  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);
  seg=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);

  /* z range */
  z_min = 2000 ;
  z_max = -100 ;
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	if( src[k][i][j]>10 ) 
	  {
	    if(k<z_min) z_min=k ;
	    if(k>z_max) z_max=k ;
	  }
  printf("z_min=%d z_max=%d dif=%d\n", z_min, z_max, z_max-z_min) ;


  /* x, y range */
  x_min = 2000 ;  x_max = -100 ;
  y_min = 2000 ;  y_max = -100 ;
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	if( src[k][i][j]>10 ) 
	  {
	    if(i<x_min) x_min=i ;
	    if(i>x_max) x_max=i ;

	    if(j<y_min) y_min=j ;
	    if(j>y_max) y_max=j ;
	  }
  printf("x: (%d->%d)\ny: (%d->%d)\n", x_min, x_max, y_min, y_max );



  /* center of the input object */
  sum.x=sum.y=sum.z=total=0.0;
      
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if(src[k][i][j]>10)
	    {
	      sum.x+=((float) (i));
	      sum.y+=((float) (j));
	      sum.z+=((float) (k));
	      total+=1.0;
	    }
	}
  center.x=sum.x/total;
  center.y=sum.y/total;
  center.z=sum.z/total;
  printf("center.x=%f center.y=%f center.z=%f\n",center.x,center.y,center.z);

  
  /* extraction */ 
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  x = static_cast<int>(i-center.x);
	  y = static_cast<int>(j-center.y);
	  z = static_cast<int>(k-center.z);
	  if( fabs(x)<sX && fabs(y)<sY && (z<sZ-3 && z>-sZ))
	    seg[k][i][j] = src[k][i][j] ;
	  else 
	    seg[k][i][j] = 0 ;
	}

  /* write the extraction results */
  printf("write cropping result ...\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(seg[k][i],1,image_size,fp);
  fclose(fp);

  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(seg, z_size, image_size) ;
}

void  EroseImageFrom0Background(int argc, char *argv[], int Erose_times)
{
  unsigned char ***src, ***res, ***temp ;
  int i,j,k, image_size,z_size,x,y,z,t, sx,sy,sz ;
  FILE *fp;
  float total ;
  int  state ;
  Fvector3d center,sum;


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

  /* image size in xy plane */
  image_size=imageSz;

  fp=myopen(argv[1],"r");  printf("%s\n",argv[1]) ;
  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);
  res=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);
  

  /* center of the input object */
  sum.x=sum.y=sum.z=total=0.0;
      
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if(src[k][i][j]>10)
	    {
	      sum.x+=((float) (i));
	      sum.y+=((float) (j));
	      sum.z+=((float) (k));
	      total+=1.0;
	    }
	}
  center.x=sum.x/total;
  center.y=sum.y/total;
  center.z=sum.z/total;
  printf("center.x=%f center.y=%f center.z=%f\n",center.x,center.y,center.z);


  /* extraction */ 
  for(t=0;t<Erose_times;t++)
    {
      for(k=0;k<z_size;k++)
	for(i=0;i<image_size;i++)
	  for(j=0;j<image_size;j++)
	    {
	      res[k][i][j] = src[k][i][j] ;

	      if( src[k][i][j]>0 && src[k][i][j]<100 ) /* not white and grey matter */
		{
		  state = 0 ;
		  for(sz=k-1;sz<=k+1;sz++)
		    for(sx=i-1;sx<=i+1;sx++)
		      for(sy=j-1;sy<=j+1;sy++)
			if(sz>=0 && sz<z_size && sx>=0 && sx<image_size && sy>=0 && sy<image_size && src[sz][sx][sy]==0) /*beackground*/
			  state ++ ;
		  if( state>=1 )       
		    {
		      res[k][i][j]=0 ;

		      if( k>center.z+13 || i>center.x+55 || 
			  (fabs(j-center.y)<=2 && k>center.z+2) ||
			  ((k>center.z+10 || i>center.x+40) && fabs(j-center.y)<=6) )   
			src[k][i][j]=0 ;
		    }
		}   
	    }

      /*exchange pointers*/
      temp = src ;
      src = res ;
      res = temp ;
    }
  
  
  /* write the extraction results */
  printf("write erose result ...\n") ;
  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(src, z_size, image_size) ;
  UCfree3d(res, z_size, image_size) ;
}

void  UrethraSmooth(int argc, char *argv[], int size)
{
  unsigned char ***src, ***res ;
  int i,j,k, kk,image_size,z_size,x,y,z,t, sx,sy,sz ;
  FILE *fp;
  float total ;
  int  state, num ;
  Fvector3d center,sum;


  printf("Urethra smoothing\n") ;

  /* image size in xy plane */
  image_size=imageSz;

  fp=myopen(argv[1],"r");  printf("%s\n",argv[1]) ;
  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);
  res=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);
  

  for(k=0;k<z_size;k++)
    {
      printf("k=%d\n", k) ;
      for(i=0;i<image_size;i++)
	for(j=0;j<image_size;j++)
	  {
	    total = 0 ; num=0;
	    for(kk=k-size;kk<=k+size;kk++)
	      if(kk>=0 && kk<z_size) 
		{
		  num ++ ;
		  total += src[kk][i][j] ;
		}
	    res[k][i][j] = static_cast<unsigned char>(total/num);
	  }
    }

  /* write the extraction results */
  printf("write smooth result ...\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      fwrite(res[k][i],1,image_size,fp);
  fclose(fp);

  /* free */
  UCfree3d(src, z_size, image_size) ;
  UCfree3d(res, z_size, image_size) ;
}



void SampleDownImg(int argc, char *argv[])
{
  unsigned char ***src, ***outp;
  int  i,j,k, image_size,z_size, temp ;
  float  ns ;
  FILE *fp;
  long idum, intensity ;
  float Scale_image, Scale_z ;
  int out_image_size, out_z_size ;

  sscanf(optarg, "%f,%f,%d", &Scale_image, &Scale_z, &image_size) ;

  /* image size in xy plane */
  /*  image_size=imageSz;*/


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

  out_image_size = static_cast<int>(image_size * Scale_image);
  out_z_size     = static_cast<int>(z_size * Scale_z);

  printf("Input  image size: %d, %d, %d\n", image_size, image_size, z_size) ;
  printf("Output image size: %d, %d, %d\n", out_image_size, out_image_size, out_z_size) ;

  /* apply for memory */ 
  src    =UCalloc3d(image_size,image_size,z_size);
  outp   =UCalloc3d(out_image_size,out_image_size,out_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);

  
  /* sample ...  */ 
  for(k=0;k<out_z_size;k++)
    for(i=0;i<out_image_size;i++)
      for(j=0;j<out_image_size;j++)
	outp[k][i][j] = src[(int)(k/Scale_z)][(int)(i/Scale_image)][(int)(j/Scale_image)];


  /* write the result */
  fp=myopen(argv[2],"w");
  for(k=0;k<out_z_size;k++)
    for(i=0;i<out_image_size;i++)
      fwrite(outp[k][i],1,out_image_size,fp);
  fclose(fp);

  /* free */
  UCfree3d(src,     z_size, image_size) ;
  UCfree3d(outp,    out_z_size, out_image_size) ;
}


void AverageOfSamples(int argc, char *argv[])
{
  unsigned char ***src, ***res ;
  int   ***AverageImg ;
  int   max; 
  int   i,j,k,x,y,z, nk ;
  FILE  *fp;
  char  filename[80];
  int   image_size,z_size ;
  int   NumOfSamples, MemoryApplicationStatus=FAIL ;
  int   MaxSliceNum=120 ;
  UCvector3d ***TypeOnEachPixel ;

  /* image size */
  image_size=imageSz;
  /* number of samples */
  sscanf(optarg, "%d", &NumOfSamples) ;  printf("NumOfSamples=%d\n", NumOfSamples) ;
  /* sample file name */
  const char* SampleFile = argv[1];
  printf("%s\n", SampleFile) ;

  /* average space */
  AverageImg = Ialloc3d(image_size,image_size, MaxSliceNum);
  res        = UCalloc3d(image_size,image_size, MaxSliceNum);
  TypeOnEachPixel=UCvector3dalloc3d(image_size,image_size, 40);

  MaxSliceNum = 0 ;
  for(nk=0; nk<NumOfSamples; nk++) 
    {
      /* read data*/
      sprintf(filename, "%s%d", SampleFile, nk+1); printf("%s\n", filename) ;
      fp=myopen(filename,"r");
      fseek(fp,0,SEEK_END);
      z_size=ftell(fp)/(image_size*image_size);
      rewind(fp);

      if(z_size>MaxSliceNum) MaxSliceNum=z_size ;

      if(MemoryApplicationStatus==FAIL)
	{
	  src = UCalloc3d(image_size,image_size,z_size);
	  MemoryApplicationStatus=OOKK ;
	}

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

      for(z=0;z<z_size;z++)
	for(x=0;x<image_size;x++)
	  for(y=0;y<image_size;y++)
	    {
	      AverageImg[z][x][y] += src[z][x][y] ;

	      /* other average choice */
	      if( src[z][x][y]==0 ) /* background */
		TypeOnEachPixel[z][x][y].x ++ ;
	      if( src[z][x][y]==255 ) /* white stuff */
		TypeOnEachPixel[z][x][y].y ++ ;
	      if( src[z][x][y]==204 ) /* grey stuff */
		TypeOnEachPixel[z][x][y].z ++ ;
	    }
    }

  /* average */
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  /*res[z][x][y] = AverageImg[z][x][y]/NumOfSamples ;*/

	  /* other choise */
	  if( TypeOnEachPixel[z][x][y].x>TypeOnEachPixel[z][x][y].y && TypeOnEachPixel[z][x][y].x>TypeOnEachPixel[z][x][y].z) 
	    res[z][x][y] = 0 ;
	  if( TypeOnEachPixel[z][x][y].y>=TypeOnEachPixel[z][x][y].x-1 && TypeOnEachPixel[z][x][y].y>=TypeOnEachPixel[z][x][y].z-2) 
	    res[z][x][y] = 255 ;
	  if( TypeOnEachPixel[z][x][y].z>=TypeOnEachPixel[z][x][y].x-1 && TypeOnEachPixel[z][x][y].z>TypeOnEachPixel[z][x][y].y+2) 
	    res[z][x][y] = 204 ;
	}
	       
  /* write the extraction results */
  printf("write average image ...\n") ;
  fp=myopen(argv[2],"w");
  for(k=0;k<MaxSliceNum;k++)
    for(i=0;i<image_size;i++)
      fwrite(res[k][i],1,image_size,fp);
  fclose(fp);
}

void TestInReadingFloat(int argc, char *argv[])
{
  Fvector3d ***src ;
  int   i,j,k,x,y,z ;
  int   image_size,z_size ;
  FILE  *fp;

  /* image size */
  image_size=128;
  z_size=61 ;
  /* number of samples */
  sscanf(optarg, "%d,%d,%d", &i, &j, &k) ;  printf("position (%d,%d,%d)\n", i,j,k) ;


  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename);
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  rewind(fp);
  
  src = Fvector3dalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],sizeof(Fvector3d),image_size,fp);
  fclose(fp);

  /* printf()*/
  printf("def[%d][%d][%d]=(%f,%f,%f)\n", i,j,k, src[k][i][j].x, src[k][i][j].y, src[k][i][j].z) ;
}

void PicKBrainRegionSaveIt(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  Fvector3d position ;


  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename);
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);


  /* save WM, GM, CSF regions  */
  printf("write brain coordinates ...\n") ;
  fp=myopen(argv[2],"w");
   for(z=0;z<z_size;z++)
     for(x=0;x<image_size;x++)
       for(y=0;y<image_size;y++)
	 {
	   if( src[z][x][y]==WM || src[z][x][y]==GM || src[z][x][y]==VN )
	     {
	       position.x = x ;
	       position.y = y ;
	       position.z = z ;
	       fwrite(&position,sizeof(Fvector3d),1,fp);
	     }
	 }
   fclose(fp);
}

void SelectCertainRangeIntensities_KeepOriginalIntensities(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  int   Increased, current, max, min ;
  float gamma, rate ;


  sscanf(optarg, "%d,%d", &min, &max) ;  printf("min=%d max=%d\n", min, max) ;


  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename);
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);


  /* change intensity */
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	if( src[z][x][y]<=max && src[z][x][y]>=min )
	  src[z][x][y] = src[z][x][y] ;
	else 
	  src[z][x][y] = 0 ;

  /* write the extraction results */
  printf("write  binary image ...\n") ;
  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);
}


void SelectCertainRangeIntensities(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  int   Increased, current, max, min ;
  float gamma, rate ;


  sscanf(optarg, "%d,%d", &min, &max) ;  printf("min=%d max=%d\n", min, max) ;


  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename);
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);


  /* change intensity */
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	if( src[z][x][y]<=max && src[z][x][y]>=min )
	  src[z][x][y] = 255 ;
	else 
	  src[z][x][y] = 0 ;

  /* write the extraction results */
  printf("write  binary image ...\n") ;
  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);
}


void IncreaseIntensityAtEachVoxel(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  int   Increased, current, max, min ;
  float gamma, rate ;

  /* image size */
  image_size=imageSz;
  /* number of samples */
  sscanf(optarg, "%d, %f", &Increased, &gamma) ;  printf("SliceNo=%d  gamma=%f\n", Increased, gamma) ;


  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);


  /* change intensity */
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  if( src[z][x][y]!=0 )
	    {
	      current = src[z][x][y]+ Increased ;
	      
	      if( current>255 ) current=255 ;
	      if( current<0 )   current=0;

	      src[z][x][y] = current ;
	    }
	}


  printf("gamma=%f\n", gamma) ;
  /* mapping */
    for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  /*current = 255*pow((float)src[z][x][y]/255.0, 1.0/gamma) ;*/
	  current = static_cast<int>(src[z][x][y]*1.5); /*2.5*/
	  if( current>255 ) current=255 ;
	  if( current<0 )   current=0;
	  src[z][x][y] = static_cast<unsigned char>(current);
	}
  
  /* looking for min and max */
  /*  min = 2000 ;
  max = -1000 ;
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  if( src[z][x][y]<min ) min=src[z][x][y] ;
	  if( src[z][x][y]>max ) max=src[z][x][y] ; 
	}
  rate = 255.0/(max-min) ;
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  current = (src[z][x][y]-min)*rate ;
	  
	  if( current>255 ) current=255 ;
	  if( current<0 )   current=0;
	  
	  src[z][x][y] = current ;
	}
  */

  /* write the extraction results */
  printf("write image ...\n") ;
  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);
}


void MappingIntensityAtEachVoxel(int argc, char *argv[], int Raise)
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  int   Increased, current, max, min ;
  float gamma, rate ;

  /* image size */
  image_size=imageSz;
  /* number of samples */
  sscanf(optarg, "%d, %f", &Increased, &gamma) ;  printf("SliceNo=%d  gamma=%f\n", Increased, gamma) ;


  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);


  printf("joint=%d\n", Raise) ;
  /* mapping */
    for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  if( src[z][x][y]>0 )
	    {
	      current = Raise + src[z][x][y] ;
	      if( current>255 ) current=255 ;
	  
	      src[z][x][y] = current ;
	    }
	}

  /* write the extraction results */
  printf("write image ...\n") ;
  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);
}


void RemoveCertainSliceOfImg(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  int   SliceNo ;

  /* image size */
  image_size=imageSz;
  /* number of samples */
  sscanf(optarg, "%d", &SliceNo) ;  printf("SliceNo=%d\n", SliceNo) ;


  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

  /* remove */
  for(x=0;x<image_size;x++)
    for(y=0;y<image_size;y++)
      src[SliceNo][x][y] = 0 ;
	       
  /* write the extraction results */
  printf("write image ...\n") ;
  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);
}


void MappingIntensitiesBasedOnNearestDistance(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;
  int   SliceNo ;
  float distanceVN, distanceGM, distanceWM, distanceBG, distanceCSF;

  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

  /* mapping  */
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      for(y=0;y<image_size;y++)
	{
	  distanceBG = (src[z][x][y]-BG)*(src[z][x][y]-BG) ;
	  distanceCSF= (src[z][x][y]-CSF)*(src[z][x][y]-CSF) ;
	  distanceVN = (src[z][x][y]-VN)*(src[z][x][y]-VN) ;
	  distanceGM = (src[z][x][y]-GM)*(src[z][x][y]-GM) ;
	  distanceWM = (src[z][x][y]-WM)*(src[z][x][y]-WM) ;

	  if( distanceVN<distanceGM && distanceVN<distanceWM )
	    src[z][x][y] = VN ;
	  if( distanceGM<distanceVN && distanceGM<distanceWM )
	    src[z][x][y] = GM ;
	  if( distanceWM<distanceVN && distanceWM<distanceGM )
	    src[z][x][y] = WM ;
	  if( sqrt(distanceBG)<5 )
	    src[z][x][y] = BG ;
	  if( sqrt(distanceCSF)<5 )
	    src[z][x][y] = CSF ;
	}
	       
  /* write the extraction results */
  printf("write image ...\n") ;
  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);
}

void CropCertainSliceAndZoomItOutputIt(int argc, char *argv[], char WhichDimension, int Slice, float zoomm)
{
  unsigned char ***src ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size, out_image_size, out_z_size ;
  unsigned char **SliceImg ;
  float xres, yres, zres ;

  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

  /* zoomed size */
  xres = zoomm ;
  yres = zoomm ;
  zres = zoomm ; /**1.5/0.9375/1.5 ;*/
  out_z_size     = static_cast<int>(zres*z_size);
  out_image_size = static_cast<int>(xres*image_size);

  if( WhichDimension=='x' || WhichDimension=='X' )
    {
      SliceImg = UCalloc2d(out_z_size, out_image_size);
      x = Slice ; printf("Slice=%d\n", Slice) ;

      printf("zoomed size: %d %d\n\n", out_z_size, out_image_size ) ; 

	for(k=0;k<out_z_size;k++)
	  for(j=0;j<out_image_size;j++)
	  {
	    z = (int)( k/zres ) ;
	    y = (int)( j/yres ) ;
	    
	    if( z<z_size && y<image_size )
	      SliceImg[k][j] = src[z][x][y];
	    else 
	      SliceImg[k][j] = 0 ;
	  }

	/* write the extraction results */
	printf("write image ...\n") ;
	fp=myopen(argv[2],"w");
	for(k=0;k<out_z_size;k++)
	    fwrite(SliceImg[k],1,out_image_size,fp);
	fclose(fp);
    }

  if( WhichDimension=='y' || WhichDimension=='Y' )
    {
	SliceImg = UCalloc2d(out_z_size, out_image_size);
	y = Slice ;
	printf("zoomed size: %d %d\n\n", out_z_size, out_image_size ) ; 

	for(k=0;k<out_z_size;k++)
	  for(i=0;i<out_image_size;i++)
	  {
	    z = (int)( k/zres ) ;
	    x = (int)( i/yres ) ;
	    
	    if( z<z_size && x<image_size )
	      SliceImg[k][i] = src[z][x][y];
	    else 
	      SliceImg[k][i] = 0 ;
	  }

	/* write the extraction results */
	printf("write image ...\n") ;
	fp=myopen(argv[2],"w");
	for(k=0;k<out_z_size;k++)
	    fwrite(SliceImg[k],1,out_image_size,fp);
	fclose(fp);
    }

  if( WhichDimension=='z' || WhichDimension=='Z' )
    {
	SliceImg = UCalloc2d(out_image_size, out_image_size);
	z = Slice ;
	printf("zoomed size: %d %d\n\n", out_image_size, out_image_size ) ; 

	for(i=0;i<out_image_size;i++)
	  for(j=0;j<out_image_size;j++)
	  {
	    x = (int)( i/yres ) ;
	    y = (int)( j/yres ) ;
	    
	    if( y<image_size && x<image_size )
	      SliceImg[i][j] = src[z][x][y];
	    else 
	      SliceImg[i][j] = 0 ;
	  }

	/* write the extraction results */
	printf("write image ...\n") ;
	fp=myopen(argv[2],"w");
	for(i=0;i<out_image_size;i++)
	    fwrite(SliceImg[i],1,out_image_size,fp);
	fclose(fp);
    }
}


void AddBlankSlice(int RequiredSliceNo, int argc, char *argv[])
{
  unsigned char ***src, ***blank ;
  int   i,j,k,x,y,z ;
  FILE  *fp;
  int   image_size,z_size ;

  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

 	       
  /* write the extraction results */
  printf("write image ...\n") ;
  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);
  if( RequiredSliceNo>z_size )
    {
      blank = UCalloc3d(image_size,image_size,RequiredSliceNo-z_size);
      for(k=z_size;k<RequiredSliceNo;k++)
	for(i=0;i<image_size;i++)
	  fwrite(blank[k-z_size][i],1,image_size,fp);
    }
  fclose(fp);
}

void FindMaxMinFromImg(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z;
  FILE  *fp;
  int   image_size,z_size ;
  int max, min ;

  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

 	       
  /* find max and min */
  max = 0 ;
  min = 1000 ;
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if( src[k][i][j]>max ) 
	    max = src[k][i][j] ;

	  if( src[k][i][j]<min && src[k][i][j]>0 ) 
	    min = src[k][i][j] ;
	}
  printf("max=%d min=%d\n", max, min) ;
}


void TaotalValuesInImage(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z;
  FILE  *fp;
  int   image_size,z_size ;
  double TotalIntensities = 0 ;
  int   num ;

  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

 	       
  /* total intensities  */
  TotalIntensities = 0 ;
  num = 0 ;
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  TotalIntensities += src[k][i][j] ;
	  if( src[k][i][j]>0 )
	    num ++ ;
	}

  printf("TotalIntensities=%f\n", TotalIntensities) ;
  printf("average Intensities=%f\n", TotalIntensities/num) ;
  printf("total number in the regions: %d\n", num );
}

void GetVolumeOfVoxelsWithCertainLabel(int Labels, int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z;
  FILE  *fp;
  int   image_size,z_size ;
  int   volume, size ;
  double average ;


  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

 	       
  /* histogram */
  volume = 0 ;

  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	if( src[k][i][j]==Labels ) 
	  volume ++ ;
  printf("                                 Label=%d   ->  Volume=%d\n", Labels, volume) ;
}


void Histogram_Distribution(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z;
  FILE  *fp;
  int   image_size,z_size ;
  double TotalIntensities = 0 ;
  int   num, Hist[256], size ;
  double average ;


  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

 	       
  /* histogram */
  average = 0;
  num = 0 ;

  for(i=0; i<256; i++)
    Hist[i] = 0 ;

  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if( src[k][i][j]>0 ) 
	    {
	      Hist[src[k][i][j]] += 1 ;

	      num ++ ;
	      average += src[k][i][j] ;
	    }
	}
  printf("average uintensity in the skull-stripped brain=%f\n", average/num) ;
  for(i=0;i<256;i++)
    printf("Number of intensity %d: %d\n", i, Hist[i]) ;


  /* output */
  filename = argv[2];
  printf("%s\n", filename) ;
  fp=myopen(filename,"w");
  for(i=0;i<256;i++)
    fprintf(fp, "%d:      %f\n", i, (float)Hist[i]/(float)num) ;
    /*fprintf(fp, "%f      %f\n", i/25.0, (float)Hist[i]/(float)num) ;*/
  fclose(fp);

  /* save the result */
  fp=myopen("LabelVol_Results.txt","w") ;
  fprintf(fp, "labels              %s\n", argv[1]) ;
  for(i=0;i<256;i++)
    if( Hist[i]!=0 )
      fprintf(fp, "%3d %22d \n", i, Hist[i]) ;
  fclose(fp);
  

  printf("total number in the regions: %d\n", num );
}


void FindTheSizeOfHead(int argc, char *argv[])
{
  unsigned char ***src ;
  int   i,j,k,x,y,z;
  FILE  *fp;
  int   image_size,z_size ;
  int   GMsize, WMsize ;

  /* image size */
  image_size=imageSz;

  /* read data*/
  const char* filename = argv[1];
  printf("%s\n", filename) ;
  fp=myopen(filename,"r");
  fseek(fp,0,SEEK_END);
  z_size=ftell(fp)/(image_size*image_size);
  rewind(fp);
  
  src = UCalloc3d(image_size,image_size,z_size);
  
  /* read data */ 
  for(z=0;z<z_size;z++)
    for(x=0;x<image_size;x++)
      fread(src[z][x],1,image_size,fp);
  fclose(fp);

 	       
  /* find size */
  GMsize = 0 ; 
  WMsize = 0 ;
  for(k=0;k<z_size;k++)
    for(i=0;i<image_size;i++)
      for(j=0;j<image_size;j++)
	{
	  if( src[k][i][j]==WM ) 
	    WMsize ++ ;
	  if( src[k][i][j]==GM ) 
	    GMsize ++ ;
	}
  printf("%s: Total WM and GM voxel number=%d\n", filename, WMsize+GMsize) ;
  printf("WM voxel number=%d\n", WMsize) ;
  printf("GM voxel number=%d\n", GMsize) ;

  /* free */
  UCfree3d(src, z_size, image_size) ;
}



// ===========================================================================
// main
// ===========================================================================

int main(int argc,char *argv[])
{
  int   c, num, Erose_times ;
  int   level, variation, threshold, winWID, choice=-1, fold=128 ;
  float Sigma, enlarge, zoomm, CannyGuassianSigma, Smooth_Sigma, Variance, SampleSize ;
  int   BinaryThreshold;
  int   sX, sY, sZ, Ssize ;
  int   Slice;
  char  WhichDimension ;
  FILE  *fp;
  int   RequiredSliceNo, SmoothTimes, FocusIntensity ;
  int   Labels, Raise ;



  num=3;
  if(argc<num) {
    show_usage();
    exit(1);
  }


  /* functions and their parameters */

  level = 50 ;
  variation = 10 ;
  threshold = 10 ;
  Variance = 1. ;
  Sigma = 2.0 ;
  winWID = 1;
  enlarge = 1.0 ;
  zoomm = 1.0 ;
  CannyGuassianSigma =2.0 ;
  BinaryThreshold = 100 ;
  Smooth_Sigma = 1.0 ;
  SmoothTimes = 1 ;
  imageSz = 256 ;

  while((c=getopt(argc-2,argv+2,"s:l:g:n:f:w:e:d:x:y:z:r:c:k:t:b:p:aihmC:E:U:S:A:R:NZ:L:T:P:Q:BI:MHVDW:F:x:")) != -1)
    {
      switch(c)
	{
	case 'x':
	  imageSz=atoi(optarg);
	  break;

	case 'a':
	  HalfGreyLevelOfImg(argc, argv) ;
	  break;

	case 'i':
	    DoubleGreyLevelOfImg(argc, argv) ;
	  break;

	case 'h':
	  BlackReversedToWhiteInImg(argc, argv) ;
	  break;

	case 's':
	  level=atoi(optarg);
	  ExtractInterestingPart(argc, argv, level) ;
	  break;

	case 'l':
	  sscanf(optarg, "%d,%d", &level, &variation) ;
	  ExtractInterestingPartWithVariation(argc, argv, level, variation) ;
	  break;

	case 'g':
	  threshold=atoi(optarg);
	  GetEdgeFrom3DimageWithoutSmoothing(argc, argv, threshold) ;
	  break;

	case 'n':
	  Variance=atof(optarg);
	  AddNoiseTo3Dimage(argc, argv, Variance) ;
	  break;

	case 'f':
	  Sigma=atof(optarg);
	  CannyLikedEdgeDetectionWithGuassianSmoothing(argc, argv, Sigma) ;
	  break;

	case 'w':
	  winWID=atoi(optarg);
	  Denoising3DedgeImage(argc, argv, winWID) ;
	  break;

	case 'e':
	  enlarge=atof(optarg);
	  Scale3DsegmentedImage(argc, argv, enlarge) ;
	  /*Scale3DImageByLinearInterpolation(argc, argv, enlarge) ;*/
	  break;

	case 'r':
	  fold=atoi(optarg);
	  Reflect3DImageByHalfY(argc, argv, fold) ;
	  break;

	case 'd':
	  zoomm=atof(optarg);
	  choice = 1 ;
	  break;

	case 'c':
	  CannyGuassianSigma=atof(optarg);
	  Cubic_edgeDetectionBy2DcannyMethod(argc, argv, CannyGuassianSigma) ;
	  break;

	case 'k':
	  CannyGuassianSigma=atof(optarg);
	  Cubic_edgeDetectionBythree2DcannyMethod(argc, argv, CannyGuassianSigma) ;
	  break;


	case 't':
	  BinaryThreshold=atoi(optarg);
	  ThresholdingInput3DImage(argc, argv, BinaryThreshold) ;
	  break;

	case 'b':
	  Smooth_Sigma=atof(optarg);
	  Smoothing3DimageWithGuassian(argc, argv, Smooth_Sigma) ;
	  break;

	case 'p':
	  /*Smooth_Sigma=atof(optarg);*/
	  sscanf(optarg, "%f,%d,%d", &Smooth_Sigma, &SmoothTimes, &FocusIntensity) ;
	  printf("SmoothTimes=%d  Smooth_Sigma=%f FocusIntensity=%d\n", SmoothTimes, Smooth_Sigma, FocusIntensity) ;
	  Smoothing3DimageWithShapeBasedGuassian(argc, argv, Smooth_Sigma, SmoothTimes, FocusIntensity) ;
	  break;

	case 'm':
	  MappingIntensitiesOfImage(argc, argv) ;
	  break;

	case 'C':
	  sscanf(optarg, "%d,%d,%d", &sX, &sY, &sZ) ;
	  CropPartOfImgAroundCenterOfImg(argc, argv, sX, sY, sZ) ;
	  break;

	case 'E':
	  Erose_times=atoi(optarg);
	  EroseImageFrom0Background(argc, argv, Erose_times) ;
	  break;

	case 'U':
	  Ssize = atoi(optarg);
	  UrethraSmooth(argc, argv, Ssize) ;
	  break;

	case 'S':
	  /*SampleSize = atof(optarg);*/
	  SampleDownImg(argc, argv) ;
	  break;

	case 'A':
	  AverageOfSamples(argc, argv) ;
	  break;

	case 'R':
	  RemoveCertainSliceOfImg(argc, argv) ;
	  break;

	case 'N':
	  MappingIntensitiesBasedOnNearestDistance(argc, argv) ;
	  break;

	case 'Z':
	  sscanf(optarg, "%c,%d,%f", &WhichDimension, &Slice, &zoomm) ;
	  CropCertainSliceAndZoomItOutputIt(argc, argv, WhichDimension, Slice, zoomm) ;
	  break;

	case 'L':
	  IncreaseIntensityAtEachVoxel(argc, argv) ;
	  break;

	case 'T':
	  TestInReadingFloat(argc, argv) ;
	  break;

	case 'P':
	  SelectCertainRangeIntensities(argc, argv) ;
	  break;

	case 'Q':
	  SelectCertainRangeIntensities_KeepOriginalIntensities(argc, argv) ;
	  break;

	case 'B':
	  PicKBrainRegionSaveIt(argc, argv) ;
	  break;

	case 'I':
	  RequiredSliceNo=atoi(optarg); 
	  AddBlankSlice(RequiredSliceNo, argc, argv) ;
	  break;

	case 'M':
	  FindMaxMinFromImg(argc, argv) ;
	  break;

	case 'H':
	  FindTheSizeOfHead(argc, argv) ;
	  break;

	case 'V':
	  TaotalValuesInImage(argc, argv) ;
	  break;

	case 'D':
	  Histogram_Distribution(argc, argv) ;
	  break;

	case 'W':
	  Labels=atoi(optarg); 
	  GetVolumeOfVoxelsWithCertainLabel(Labels, argc, argv) ;
	  break;

	case 'F':
	  Raise = atoi(optarg) ;
	  MappingIntensityAtEachVoxel(argc, argv, Raise) ;
	  break;

	default:
	  break;
	}
    }


  if(choice==1)
    Dig3DpartANDscaleIt(argc, argv,  zoomm) ;

  return 0;
}

