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

              	           I S T H M U S . C 

*****************************************************************************/
/*****************************************************************************
FILENAME   : isthmus.c
PURPOSE    : Isthmus Finding functions           
*****************************************************************************/

#include </sys/h/types.h>
#include <stdio.h>
#include <math.h>
#include "isthdefs.h"
#include "statline.h"
#include "piece.h"

extern float euclid_dist();
extern float distance();

extern int border_x[512][512];
extern int border_y[512][512];
extern int temp_image[512][512];
extern int dist_image[512][512];
extern int localx[2500], localy[2500], locallength;

extern int inverse_isthmus;
extern int draw_is_on;
extern char strg[79];
extern int p[10];  /* from thin.c */

float sqrt_of_two;
int totalskel = 0;
int first_skel = 1;
int skellength;
int upcnt = 0;
int downcnt = 0;
int prevdist;
int savex = 0;
int savey =0;
int skelx=0;
int persegment = 0;
int sum, biggest=0, otherx, othery;
int sub;

/*****************************************************************************
                     MAKE_GRAPH 
*****************************************************************************/
make_graph (deltax, deltay, minx, miny, maxx, maxy)
int deltax, deltay, minx, miny, maxx, maxy;
{ 
  int x, y,i,j, z ,length;
  int xx, yy, looky;
  int grey = 1; 
  int h;

  first_skel = 1; 
  if (inverse_isthmus) 
  {  
    h = (maxx - minx) + (maxx + 25);
    looky = (maxy + miny) / 2;  
    /* burning the contour of outter edges of the thinned skeleton */ 
VIS$BORDER_FOLLOW(&length, localx, localy, grey, 7, maxx+25, looky, h, maxy, 1);
    printit ("after burning outside edge");
    for (i=1; i<length; i++)
    { 
       x = localx[i] - maxx - 25;
       y = localy[i] - miny;
       temp_image[y][x] = 7; 
    }
    for (i=1; i<length; i++)
    { 
       x = localx[i] - maxx - 25;
       y = localy[i] - miny;
#ifdef ROGER
       printf ("looking at x: %d y: %d\n", x, y); 
       printf ("looking at localx: %d localy: %d\n", localx[i], localy[i]); 
#endif
       VIS$SET_PIXEL (x + maxx + 25, y + miny, 254);
  
   for (z=1;z<=8;z++)
   {
      get_next_neighbor (z, x, y, &xx, &yy);
      if (temp_image[yy][xx] == 1) 
      {
	 skeleton_follow(yy, xx, minx, miny, maxx);
         temp_image[yy][xx] = 0;  /* added */
       } /* endif temp */
      } /* endfor z loop */
    } /* endfor i loop */
    piece[currpiece].inv_skelx[0] = totalskel;  /* count */
    totalskel = 0; 
  } /* endif */
 else           /* POSITIVE Isthmus */
 {
  looky = (maxy - miny) / 3;  /* 1/3 the way */
  for (j=looky; j< deltay; j++) 
    for (i=35; i<deltax; i++)
    { 
       if (temp_image[j][i] == 1) 
       {
	  if (VIS$GET_PIXEL (i + maxx + 25, j + miny) == 1)
	  {  
	     skeleton_follow(j,i, minx, miny, maxx);
	     piece[currpiece].pos_skelx[0] = totalskel;  /* count */
	     totalskel = 0; 
	     return(1);  
          } 
	}
     } 
 } /* end else */
} /* end function */

/*****************************************************************************
                     SKELETON_FOLLOW 
*****************************************************************************/
skeleton_follow (j,i,minx, miny, maxx)
int j,i, minx, miny, maxx;
{
  int n, x,y,z, t, grey; 
  int xx, yy; 
  int rx, ry;
  int skelcolor = 255;

  VIS$SET_PIXEL (i + maxx + 25,j + miny, skelcolor); 
  if (draw_is_on)
  {
     drawline (border_x[j][i]+minx, border_y[j][i]+miny,
	           i+minx, j+miny, 255);
  }
   
  if (first_skel)
  {
    first_skel = 0; 
    persegment = 0; 
    prevdist = dist_image[j][i];
    if (inverse_isthmus)
       drawpiece(currpiece, 255);
   }/* end if firstskel */
  
  if (!inverse_isthmus)
  {
     totalskel++;
     piece[currpiece].pos_skelx[totalskel] = i + minx;
     piece[currpiece].pos_skely[totalskel] = j + miny;
  }
  else
  {
     totalskel++;
     piece[currpiece].inv_skelx[totalskel] = i + minx;
     piece[currpiece].inv_skely[totalskel] = j + miny;
  }

  for (z=1;z<=8;z++)
  {
    get_next_neighbor (z, i, j, &x, &y);
    if (temp_image[y][x] == 1 && 
	(VIS$GET_PIXEL(x + maxx + 25, y + miny) == 1) )
    { 
      p[2] = temp_image[j-1][i]; 
      p[3] = temp_image[j-1][i+1]; 
      p[4] = temp_image[j][i+1]; 
      p[5] = temp_image[j+1][i+1]; 
      p[6] = temp_image[j+1][i]; 
      p[7] = temp_image[j+1][i-1]; 
      p[8] = temp_image[j][i-1]; 
      p[9] = temp_image[j-1][i-1]; 

      n = p[2]+p[3]+p[4]+p[5]+p[6]+p[7]+p[8]+p[9];
    
    if (!persegment)
    {
      if (upcnt > UPTHRESH && savex !=0) /* SAVE THE ISTHMUS POINTS */
      {
	 for (yy=savey-1; yy<=savey+1; yy++)
           for (xx=savex-1; xx<=savex+1;xx++)
            { 
	      sum = abs (border_x[yy][xx] - border_x[savey][savex]) + 
		    abs(border_y[yy][xx] - border_y[savey][savex]);
              if (sum > biggest)
              {
                 biggest = sum;
	         otherx = border_x[yy][xx] + minx; /* + 2; /* + 2 RWW*/
		 othery = border_y[yy][xx] + miny;
	      }
	    } /* end for xx */ 
	 biggest = 0;
	 find_pts_on_contour (border_x[savey][savex]+minx/*+2*/, 
			      border_y[savey][savex]+miny,
	                      otherx, othery, currpiece);
	 upcnt = downcnt = savex = 0; 
         persegment = 1; 
      }
      if (prevdist < dist_image[j][i]
	  || prevdist == dist_image[j][i]) /* going up or same */
      {  
	if (downcnt > 3)
        {
	  savex = i/*-2*/;
          skelx = i/*-2*/;
	  savey = j;
#ifdef ROGER
	  printf ("saved it at savex: %d savey: %d\n", savex, savey); 
#endif
          if (prevdist < dist_image[j][i])
	  {  upcnt = 1;
	  }
	  downcnt = 0; 
	}
	else
        if (prevdist < dist_image[j][i])
	{ 
	  upcnt++; 
        } 
      }
      else
      if (prevdist > dist_image[j][i]) /* going down */
      {  downcnt++;
         upcnt = 0; 
      }
   } /* end persegment */
      
      if (n == 2 || n == 1)
      { 
	 if (n == 1)
         {
	    prevdist = dist_image[y][x];
	    upcnt = downcnt = 0; 
	    savex = 0;
	    persegment = 0; 
	    VIS$SET_PIXEL (i + maxx + 25,j + miny, skelcolor); 
	    prevdist = dist_image[j][i];
            skeleton_follow(y,x, minx, miny, maxx);
            return(0); 
	 }
         else
	 {
	    prevdist = dist_image[j][i];
            skeleton_follow(y,x, minx, miny, maxx);
	 }
      }
      else
       if (n > 2)
       { 
	 prevdist = dist_image[y][x];
	 VIS$SET_PIXEL (i + maxx + 25,j + miny, skelcolor); 

#ifdef ROGER
    printf ("INTERSECTION POINT\n");
#endif
         if (persegment)
	   if (inverse_isthmus)
	   {
	     sub = piece[currpiece].inv_isthnum; 
	     skellength = (int)euclid_dist(skelx, savey,i/*-2*/,j);  
	     piece[currpiece].inv_isth[sub].skeldist = skellength; 
	     piece[currpiece].inv_isth[sub].branch_index = totalskel;
	     sprintf (strg, "euclidean skel dist: %d", skellength); 
             printit (strg);
#ifdef ROGER
          printf ("euclidean skeldist: %d\n", 
	          piece[currpiece].inv_isth[sub].skeldist);
	  printf ("saved it at skelx: %d savey: %d\n", skelx, savey); 
	  printf ("i-2: %d j: %d\n", i-2, j); 
	  printf ("\n");
	  waitforkey(); 
#endif
	   }
           else
	   {
	     sub = piece[currpiece].isthnum; 
	     skellength = (int)euclid_dist(skelx,savey,i/*-2*/,j);  
	     piece[currpiece].dmtisth[sub].skeldist = skellength; 
	     piece[currpiece].dmtisth[sub].branch_index = totalskel;
	     sprintf (strg, "euclidean skel dist: %d", skellength); 
             printit (strg);
#ifdef ROGER
         printf ("skeldist: %d\n", piece[currpiece].dmtisth[sub].skeldist);
	  printf ("saved it at skelx: %d savey: %d\n", skelx, savey); 
	  printf ("i-2: %d j: %d\n", i-2, j); 
	   printf ("\n");
	   waitforkey(); 
#endif
	   }/* end else inver isthmus */
	 upcnt = downcnt = 0; 
	 savex = 0;
         persegment = 0; 
	 
         skeleton_follow(y,x, minx, miny, maxx);
	 upcnt = downcnt = 0; 
         persegment = 0; 
	 savex = 0;
#ifdef ROGER
    printf ("out of recursion at x:%d y:%d\n", i + maxx + 25,j+miny);
#endif
       } /* endif */
    } /* end if not temp_image== 1 */
  } /* end for */ 

}  /* end skeleton follow */
  

/*****************************************************************************
                   FIND_PTS_ON_CONTOUR
*****************************************************************************/
find_pts_on_contour (x1, y1, x2, y2, pcnum)
int x1, y1, x2, y2, pcnum;
{
  int index1=0, index2=0, temp, numisth;
  int newindex1, newindex2;
  int i;
  int one_diff=512;
  int two_diff=512;
  int one_tdiff, two_tdiff;
#ifdef ROGER
  printf ("starting to look pts for on the CONTOUR\n");
  printf ("looking for x1: %d y1: %d\n", x1, y1);
  printf ("looking for x2: %d y2: %d\n", x2, y2);
#endif
  for(i=1;i<=piece[pcnum].pathcnt;i++)
  {  
     one_tdiff = abs (x1 - (int)piece[pcnum].avgx[i]) + 
                 abs (y1 - (int)piece[pcnum].avgy[i]);
       if (one_tdiff < CLOSE_TO_CONTOUR)
       if (one_tdiff < one_diff)
       {
#ifdef ROGER
 printf ("matched x1,y1 at: x[%d]: %d y[%d]: %d\n", i, 
	      (int)piece[pcnum].avgx[i],i, (int)piece[pcnum].avgy[i]);
#endif
          one_diff = one_tdiff;
          index1 = i; 
       } 
     two_tdiff = abs (x2 - (int)piece[pcnum].avgx[i]) + 
                 abs (y2 - (int)piece[pcnum].avgy[i]);
       if (two_tdiff < CLOSE_TO_CONTOUR)
       if (two_tdiff < two_diff)
       {
#ifdef ROGER
 printf ("matched x2,y2 at: x[%d]: %d y[%d]: %d\n", i, 
	      (int)piece[pcnum].avgx[i],i, (int)piece[pcnum].avgy[i]);
#endif
          two_diff = two_tdiff;
          index2 = i; 
       } 
   } /* end for loop*/
  
   if (index1 != 0 && index2 != 0)
   {

      if (index1 > index2)
      {
        temp = index1;
        index1 = index2;
        index2 = temp;
      }
      if ( (index2 - index1) > 0.33 * piece[pcnum].pathcnt)
	 return(0);
      
      narrowest (pcnum, index1, index2, &newindex1, &newindex2); 
      accenthoriz ((int)piece[pcnum].avgx[newindex1], 
		   (int)piece[pcnum].avgy[newindex1], 1); 
      accenthoriz ((int)piece[pcnum].avgx[newindex2], 
		   (int)piece[pcnum].avgy[newindex2], 1); 
     
      drawline ( (int)piece[pcnum].avgx[newindex1], 
	         (int)piece[pcnum].avgy[newindex1], 
                 (int)piece[pcnum].avgx[newindex2], 
	         (int)piece[pcnum].avgy[newindex2], 255); 
   if (inverse_isthmus) 
   { 
     piece[pcnum].inv_isthnum++;
     numisth = piece[pcnum].inv_isthnum;
     piece[pcnum].inv_isth[numisth].indexbeg = newindex1;
     piece[pcnum].inv_isth[numisth].indexend = newindex2;
     piece[pcnum].inv_isth[numisth].lenseg   = newindex2 - newindex1;
     
     piece[pcnum].inv_isth[numisth].index_dist = 
	(int)euclid_dist( (int)piece[pcnum].avgx[newindex1], 
	                  (int)piece[pcnum].avgy[newindex1], 
                          (int)piece[pcnum].avgx[newindex2], 
	                  (int)piece[pcnum].avgy[newindex2]); 
#ifdef ROGER
     printf ("index1: %d\n", piece[pcnum].inv_isth[numisth].indexbeg);
     printf ("index2: %d\n", piece[pcnum].inv_isth[numisth].indexend);
     printf ("length: %d\n", piece[pcnum].inv_isth[numisth].lenseg);
#endif
   } 
   else
   { 
     piece[pcnum].isthnum++;
     numisth = piece[pcnum].isthnum;
     piece[pcnum].dmtisth[numisth].indexbeg = newindex1;
     piece[pcnum].dmtisth[numisth].indexend = newindex2;
     piece[pcnum].dmtisth[numisth].lenseg   = newindex2 - newindex1;
     
     piece[pcnum].dmtisth[numisth].index_dist = 
	(int)euclid_dist( (int)piece[pcnum].avgx[newindex1], 
	                  (int)piece[pcnum].avgy[newindex1], 
                          (int)piece[pcnum].avgx[newindex2], 
	                  (int)piece[pcnum].avgy[newindex2]); 
#ifdef ROGER
     printf ("index1: %d\n", piece[pcnum].dmtisth[numisth].indexbeg);
     printf ("index2: %d\n", piece[pcnum].dmtisth[numisth].indexend);
     printf ("length: %d\n", piece[pcnum].dmtisth[numisth].lenseg);
#endif
   } 
  
  } /* endif index1 && index2 != 0*/
  else    /* didnt find both indices!!! */
  {
     printf ("didnt find both indices\n");
     printf ("index1: %d\n", newindex1);
     printf ("index2: %d\n", newindex2);
  }/* end else */
} /* end function */

get_ldist (pcnum, index1, index2, len)
int pcnum, index1, index2, *len;
{
  float flen = 0.0;
  int x1, x2, y1, y2, i;
  sqrt_of_two = sqrt(2.0); 
  printf ("sqrt 2: %f\n", sqrt_of_two);

  for (i=index1; i<= index2; i++)
  {
     x1 = (int)piece[pcnum].avgx[i];
     x2 = (int)piece[pcnum].avgx[i+1];
     y1 = (int)piece[pcnum].avgy[i];
     y2 = (int)piece[pcnum].avgy[i+1];
     
     if ( abs(x1 - x2) + abs(y1 - y2) == 1)
     {  flen = flen + 1.0;
        printf ("flen: %f\n", flen); 
     }
     else
     {  flen = flen + sqrt_of_two; /* sqrt of two */
        printf ("sqrt 2 flen: %f\n", flen); 
     }
  }
  *len = (int)flen;

} /* end function */

narrowest (pcnum, index1, index2, newindex1, newindex2) 
int pcnum, index1, index2, *newindex1, *newindex2; 
{
   int x, j, x1, x2, y1, y2;
   int small=100000; 
   int dist;
   x1 = (int)piece[pcnum].avgx[index1]; 
   y1 = (int)piece[pcnum].avgy[index1];
   x2 = (int)piece[pcnum].avgx[index2]; 
   y2 = (int)piece[pcnum].avgy[index2];
   dist = (int)distance (x1, y1, x2, y2);
   /*printf ("old distance: %d\n", dist); */

   for (x=index1-NARROW; x <= index1 + NARROW; x++)
   {
      x1 = (int)piece[pcnum].avgx[x]; 
      y1 = (int)piece[pcnum].avgy[x];
      for (j=index2-NARROW; j <= index2 + NARROW; j++)
      {
        x2 = (int)piece[pcnum].avgx[j]; 
        y2 = (int)piece[pcnum].avgy[j];
        dist = (int)distance (x1, y1, x2, y2);
        if (dist < small)
	{
           small = dist;
           *newindex1 = x;
	   *newindex2 = j;
	}
      } /* endfor j */
   } /* endfor x */
   /*printf ("NEW distance: %d\n", small); */

} /* end function */

/*****************************************************************************
                     END 
*****************************************************************************/