// ----
// ----  Extract the thinness metric based skeleton. 
// ----  Uses Toriwaki and Saito's DT algorithm. 
// ----  
// ----  Implementation by : Nikhil Gagvani, Vizlab, Rutgers University
// ----
// ----  Input : Binary 3D volume with sizes. 
// ----  Output: ASCII obj file with x,y,z, DT-MNT for all object voxels
// ----
// $Id: euclidskel.c++,v 1.2 2000/09/22 02:21:23 gagvani Exp $

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define MIN(x,y) (((x) < (y))?(x):(y))
#define MAX(x,y) (((x) > (y))?(x):(y))


main(int argc, char *argv[])
{
  FILE *fin, *fout;
  unsigned char *cf;
  int L,M,N;         // Sizes in x,y,z dimensions
  int i,j,k,n;
  float *f, *buff , df, db, d, w, thresh, tot;
  long idx, slsz, sz, neibidx[26];
  int measureTime = 0;
  struct timeval tp1, tp2;
  struct timezone tz1, tz2;

  if (argc < 7)
  {
    printf("Usage: %s <volfile> <xs> <ys> <zs> <outfile> <thresh> [measureTimeFlag].\n",argv[0]);
    printf("Extract all voxels whose DT-MNT is greater than thresh and create an mskel file.\n");
    exit(1);
  }

  if ((fin = fopen(argv[1],"r")) == NULL)
  {
    printf("Cannot open %s\n",argv[1]);
    exit(1);
  }

  L = atoi(argv[2]);
  M = atoi(argv[3]);
  N = atoi(argv[4]);

  thresh = atof(argv[6]);

  if (argc > 7)
   measureTime = 1;
  
  cf = new unsigned char[L*M*N];
  f = new float[L*M*N];
  slsz = L*M;		// slice size
  sz = slsz*N;

  if ( fread(cf, sizeof(unsigned char), sz, fin) < sz)
  {
    printf("File size is not the same as volume size\n");
    delete f;
    exit(1);
  }

  if (measureTime)
    gettimeofday(&tp1, &tz1);

  for (idx = 0; idx < slsz*N; idx++)
  if (cf[idx] > 0)
    f[idx] = 5000;

  delete [] cf;  

  fclose(fin);

  int maxdim = MAX(L,M);
  maxdim = MAX(maxdim,N);
  buff = new float[maxdim+10];

  // Using Algorithm 3 from Appendix 

  // Step 1  forward scan
  
  for (k = 0; k < N; k++)
    for (j = 0; j < M; j++)
    {
       df = L;
       for (i = 0; i < L; i++)
       {
         idx = k*slsz + j*L + i;
         if (f[idx] !=0)
           df = df + 1;
         else
           df = 0;
         f[idx] = df*df;
       }
     }
 
  //  Step 1 backward scan
  
  for (k = 0; k < N; k++)
    for (j = 0; j < M; j++)
    {
      db = L;
      for (i = L-1; i >=0; i--)
      {
        idx = k*slsz + j*L + i;
        if (f[idx] !=0)
          db = db + 1;
        else
          db = 0; 
        f[idx] = MIN(f[idx],db*db);
      }
    }

  // Step 2

  for (k = 0; k < N; k++)
    for (i = 0; i < L; i++)
    {
      for (j =0; j < M; j++)
        buff[j] = f[k*slsz + j*L +i];
    
      for (j = 0; j < M; j++)
      {
        d = buff[j];
        if (d != 0)
        {
          int rmax, rstart, rend;
          rmax = (int) floor(sqrt(d)) + 1;
          rstart = MIN(rmax, (j-1));
          rend = MIN(rmax, (M-j));
          for (n = -rstart; n < rend; n++)
          {
              if (j+n >= 0 && j+n < M)
              {
                w = buff[j+n] + n*n;
                if (w < d)  d = w;
              }
          }
        }
        idx = k*slsz + j*L +i;
        f[idx] = d;
      }
    }

  // Step 3
  for (j = 0; j < M; j++)
    for (i = 0; i < L; i++)
    {
      for (k =0; k < N; k++)
        buff[k] = f[k*slsz + j*L +i];
    
      for (k = 0; k < N; k++)
      {
        d = buff[k];
        if (d != 0)
        {
          int rmax, rstart, rend;
          rmax = (int) floor(sqrt(d)) + 1;
          rstart = MIN(rmax, (k-1));
          rend = MIN(rmax, (N-k));
          for (n = -rstart; n < rend; n++)
          {
              if (k+n >= 0 && k+n < N)
              {
                w = buff[k+n] + n*n;
                if (w < d)  d = w;
              }
          }
        }
        idx = k*slsz + j*L +i;
        f[idx] = d;
      }
    }

  // Output the obj file.

  if ((fout = fopen(argv[5],"w")) == NULL)
  {
    printf("Cannot open %s for writing\n",argv[5]);
    exit(1);
  }

  for (k = 0; k < N; k++)
    for (j = 0; j < M; j++)
      for (i = 0; i < L; i++)
      {
         idx = k*slsz + j*L + i;
         if (f[idx] !=0)
         {
	   // Compute neibs of 3x3x3 cube 

	   // front face
	   neibidx[0] = (k-1)*slsz + (j-1)*L + i-1;
	   neibidx[1] = neibidx[0] + 1;
	   neibidx[2] = neibidx[0] + 2;
  	   neibidx[3] = (k-1)*slsz + j*L + i-1;
  	   neibidx[4] = neibidx[3] + 1;
  	   neibidx[5] = neibidx[3] + 2;
           neibidx[6] = (k-1)*slsz + (j+1)*L + i-1; 
           neibidx[7] = neibidx[6] + 1;
           neibidx[8] = neibidx[6] + 2;

	   // Middle face

	   neibidx[9] = k*slsz + (j-1)*L + i-1;
	   neibidx[10] = neibidx[9] + 1;
	   neibidx[11] = neibidx[9] + 2;
	   neibidx[12] = k*slsz + j*L + i-1;
	   // Skip the current voxel
	   neibidx[13] = neibidx[12]+2;
	   neibidx[14] = k*slsz + (j+1)*L + i-1; 
           neibidx[15] = neibidx[14]+1;
           neibidx[16] = neibidx[14]+2;

	   // back face
	
	   neibidx[17] = (k+1)*slsz + (j-1)*L + i-1;
           neibidx[18] = neibidx[17]+1;
 	   neibidx[19] = neibidx[17]+2;
	   neibidx[20] = (k+1)*slsz + j*L + i-1; 
	   neibidx[21] = neibidx[20]+1;
	   neibidx[22] = neibidx[20]+2;
           neibidx[23] = (k+1)*slsz + (j+1)*L + i-1; 
	   neibidx[24] = neibidx[23] + 1;
	   neibidx[25] = neibidx[23] + 2;


           // Compute the mean neighbor transform.

	   tot = 0;
	   
	   for (int l = 0; l < 26; l++)
              if (neibidx[l] >=0 && neibidx[l] < sz)
                tot+=sqrt(f[neibidx[l]]);
           //tot = (sqrt(f[idx]) - tot/26.0)/sqrt(f[idx]);
 	   tot = sqrt(f[idx]) - tot/26;
	   
           if (tot >= thresh)
	     if ( f[idx] > 1.0005 )
	       fprintf(fout,"%d %d %d %f %f\n",i,j,k,sqrt(f[idx]), tot);
         }
       }

  delete [] f;
  fclose(fout);
 
  if (measureTime)
  {
    gettimeofday(&tp2, &tz2);
    printf("Total Time for DT+skel+writeskel = %f\n",
        (tp2.tv_sec-tp1.tv_sec) + 1e-06*(tp2.tv_usec-tp1.tv_usec));
  }
}