/**************************************************************************

              	       D A N I E L S O N  algorithm 
edm.c

***************************************************************************/
/*****************************************************************************
FILENAME   : edm.c
PURPOSE    : Medial Axis Functions 
*****************************************************************************/

#include </sys/h/types.h>
#include </usr/include/curses.h>
#include <stdio.h>
#include <math.h>
#include "isthdefs.h"
#include "piece.h"
#include </usr/users/header/arrkey.h>
#include </sys/h/time.h>

#define EUCLID  1      /*   */

extern int no_circle;
extern radius;

extern int image[512][512];
extern int temp_image[512][512];
extern int border_x[512][512];
extern int border_y[512][512];
extern int deltax, deltay;
extern char string[60];
extern WINDOW *funcwin;

int Lij, Lij1, Lij2, Lij3, Lij4;
int total, dd, ee;
struct timeval f;
int oldsec;
int total_time;
int max_elev = 0;

/*****************************************************************************
                   EDM  
*****************************************************************************/
edm (minx,miny,maxx,maxy)
int minx,miny,maxx,maxy;
{
   int r, j,i,e, z;
   int y, looky;
   printit("Euclidean distance map danielson -> "); 

  deltax = (maxx - minx) + 5;
  deltay = (maxy - miny)/* + 5*/; 
  y = 0;
  for (j=miny; j<=maxy; j++)
  {
    VIS$READ_HORIZONTAL_LINE (minx, j, deltax+10, image[y]);
    y++;
  }
  if (!no_circle)
     euclid_lut (radius-1, 11); 
  printit("starting daneilson  -> "); 
  for (j=1; j<=deltay; j++)
  {
   for (i=1; i<=deltax-1; i++)
   {
       if (image[j][i] == 1) /* object */
       {
         border_x[j][i] =  512;
         border_y[j][i] =  512;
       }
       else   /* background */
       {
         border_x[j][i] =  0;
         border_y[j][i] =  0;
       }
   } /* endfor i */
  } /* end for loop */
  
  printit("done copy to L(i,j) format "); 
#ifdef ROGER
  r = gettimeofday(&f);
  oldsec = f.tv_sec;
#endif

  for (j=1; j<=deltay; j++)
  {
   for (i=1; i<=deltax-1; i++)
     do_5mins (i,j, i-1,j-1, i,j-1, i+1,j-1, i-1,j);
   for (i=1; i<=deltax-1; i++)
     do_2mins (i,j, i-1,j);
   
   for (i=deltax-2; i>=1; i--)
     do_2mins (i,j, i+1,j);
  } /* end for j loop */
 
  /*printit ("   starting bottom ");   */
  
  for (j=deltay-2; j >= 1; j--)
  {
   for (i=1; i<=deltax-1; i++)
     do_5mins (i,j, i-1,j+1, i,j+1, i+1,j+1, i-1,j);
   for (i=1; i<=deltax-1; i++)
     do_2mins (i,j, i-1,j);
   
   for (i=deltax-2; i>=1; i--)
     do_2mins (i,j, i+1,j);

  } /* end for j loop */
  
  show_distance (deltax, deltay, minx, maxx, miny, maxy);
#ifdef ROGER
  r = gettimeofday(&f);
  total_time = f.tv_sec - oldsec;
  sprintf (string, "total time is %d in seconds only\n", total_time);
  printit(string); 
#endif
  lut_roger(max_elev);  /* fix the code to get max_elev*/
  VIS$SELECT_LUT(1,14);
} /* end function */

do_5mins (i,j, i1,j1, i2,j2, i3,j3, i4,j4)
int i, j, i1,j1, i2,j2, i3,j3;
{
        Lij = (border_x[j][i] * border_x[j][i])
	       + (border_y[j][i] * border_y[j][i]);
         
        Lij1 = (   (border_x[j1][i1] + 1) * (border_x[j1][i1] + 1) )
	       + ( (border_y[j1][i1] + 1) * (border_y[j1][i1] + 1) );
	 
        Lij2 = (   (border_x[j2][i2] + 0) * (border_x[j2][i2] + 0) )
	       + ( (border_y[j2][i2] + 1) * (border_y[j2][i2] + 1) );
	 
        Lij3 = (   (border_x[j3][i3] + 1) * (border_x[j3][i3] + 1) )
	       + ( (border_y[j3][i3] + 1) * (border_y[j3][i3] + 1) );
        
	Lij4 = (   (border_x[j4][i4] + 1) * (border_x[j4][i4] + 1) )
	       + ( (border_y[j4][i4] + 0) * (border_y[j4][i4] + 0) );
    dd = min( Lij, Lij1);
    ee = min( Lij2, Lij3);
    total = min(dd, ee);
    total = min(total, Lij4);
    
    if (total == Lij)
    {
      /* printf ("itself is min \n");  */;
    }
    else 
    if (total == Lij1)
    {
        border_x[j][i] = border_x[j1][i1] + 1;
        border_y[j][i] = border_y[j1][i1] + 1;
    }
    else 
    if (total == Lij2)
    {
        border_x[j][i] = border_x[j2][i2] + 0;
        border_y[j][i] = border_y[j2][i2] + 1;
    }
    else
    if (total == Lij3)
    {
        border_x[j][i] = border_x[j3][i3] + 1;
        border_y[j][i] = border_y[j3][i3] + 1;
    }
    else
    if (total == Lij4)
    {
        border_x[j][i] = border_x[j4][i4] + 1;
        border_y[j][i] = border_y[j4][i4] + 0;
    }
} /* end function */

do_2mins (i, j, i1,j1)
int i, j, i1,j1;
{
        Lij =    (border_x[j][i] * border_x[j][i])
	       + (border_y[j][i] * border_y[j][i]);
         
        Lij1 = (   (border_x[j1][i1] + 1) * (border_x[j1][i1] + 1) )
	       + ( (border_y[j1][i1] + 0) * (border_y[j1][i1] + 0) );
	 
    dd = min( Lij, Lij1);
    if (dd == Lij1)
    {
        border_x[j][i] = border_x[j1][i1] + 1;
        border_y[j][i] = border_y[j1][i1] + 0;
    }
} /* end function */

show_edm (deltax, minx, maxx, miny, maxy)
 int deltax, minx, maxx, miny, maxy;
{
  int y, j;
  y = 0;
  for (j=miny; j<=maxy; j++)
  {
    VIS$WRITE_HORIZONTAL_LINE (maxx+50, j, deltax, border_x[y]);
    y++;
  }
  y = 0;
  for (j=miny; j<=maxy; j++)
  {
    VIS$WRITE_HORIZONTAL_LINE (maxx+150, j, deltax, border_y[y]);
    y++;
  }
}

show_distance (deltax, deltay, minx, maxx, miny, maxy)
 int deltax, deltay, minx, maxx, miny, maxy;
{
  int y, i, j;
  int square;
  y = 0;
  
  for (j=deltay-2; j >= 1; j--)
  {
   for (i=1; i<=deltax-1; i++)
   {
     square  = (border_x[j][i] * border_x[j][i])
	         + (border_y[j][i] * border_y[j][i]);
     image[j][i] = ( round ( sqrt((float)square) ) );  
     if (image[j][i] > max_elev) 
       max_elev = image[j][i];
   }
  }
  
  for (j=miny; j<=maxy; j++)
  {
    VIS$WRITE_HORIZONTAL_LINE (minx, j, deltax, image[y]);
    y++;
  }
}
/*****************************************************************************
                                END      
*****************************************************************************/