 /*
  * Khoros: $Id: isef.c,v 1.3 1992/03/20 23:42:44 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: isef.c,v 1.3 1992/03/20 23:42:44 dkhoros Exp $";
#endif

 /*
  * $Log: isef.c,v $
 * Revision 1.3  1992/03/20  23:42:44  dkhoros
 * VirtualPatch5
 *
  */ 


/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>	    file name: isef.c                                 <<<<
   >>>>                                                       <<<<
   >>>>   description: routines used by lvdrf, lvsdef, lvdrg  <<<<
   >>>>                (not documented)                       <<<<
   >>>>      routines: 					      <<<<
   >>>>                                                       <<<<
   >>>> modifications:					      <<<<
   >>>>                                                       <<<<
   >>>> gradadp() was changed so that it would not divide by  <<<<
   >>>  zero.  Mark & Charlie 	(2/15/91)     		      <<<<
   >>>>                                                       <<<<
   >>>>	modified threshold(). Looked like code that had been  <<<<
   >>>> incorrectly converted from fortran to C. see below    <<<<
   >>>> below for more details.  Mark (2/17/91)		      <<<<
   >>>>                                                       <<<<
   >>>>          NOTE: THESE ROUTINES ARE NOT FOR PUBLIC USE! <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

/* Routines used by the optimal Recursif filters */

#include <math.h>
#include "vinclude.h"


/***********************************************************************
*
*  Routine Name: expfilter()
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  
*
* Modifications:
*
***********************************************************************/

void expfilter(a0,mx,my,dim,in,out)
float a0;
int mx,my,dim;
unsigned char *in,*out;
{
int i,j,k,ki,myy,mxxx,v,v1;
myy=my-1; mxxx=mx-2;
/*--------------------------------------------------------------*/
for (i=0; i<dim; i+=mx)
 {
  *(out+i)= *(in+i);
  for (k=i+1; k<i+mx; k++)
   {
    v= *(out+k-1); v1= *(in+k);
    *(out+k)=v+a0*(v1-v)+0.5;
   }
 }
for (i=0; i<dim; i+=mx)
 for (k=i+mxxx; k>=i; k--)
  {
   v= *(out+k+1); v1= *(out+k);
   *(out+k)=v+a0*(v1-v)+0.5;
  }
for (j=0; j<mx; j++)
 for (k=j+mx; k<dim; k+=mx)
  {
   v= *(out+k-mx); v1= *(out+k);
   *(out+k)=v+a0*(v1-v)+0.5;
  }
ki=myy*mx;
for (j=ki; j<ki+mx; j++)
 for (k=j-mx; k>=0; k-=mx)
  {
   v= *(out+k+mx); v1= *(out+k);
   *(out+k)=v+a0*(v1-v)+0.5;
  }
return;
}

/***********************************************************************
*
*  Routine Name: gefextr()
*
*          Date:
*        
*       Purpose: This subprogram detects the step edges by local maximum 
*		 gradient. 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:
*
***********************************************************************/

void gefextr(dx,dy,out,mx,my)
unsigned char *dx,*dy,*out;
int mx,my;
/* *dx,*dy: input image; *out: output images;			        */
/* mx: width of the images; my: height of the image;           		*/
{
int i,j,ki,k;
unsigned char *dix,*diy,*ti1,*ti2,*tii1,*tii2,*t1;


for (i=1; i<my-1; i++)
 {
  ki=i*mx;
  for (j=1; j<mx-1; j++)
   {
    k=ki+j;
    dix=dx+k; diy=dy+k; t1=out+k;
    ti1=dix-1; ti2=dix+1;
    tii1=diy-mx; tii2=diy+mx;
/* gradient direction approaches to the horizontal direction.		*/
    if (*dix> *diy)
     if ((*dix> *ti1) && (*dix>= *ti2))
      *t1= *dix;
      else
       *t1= 0;
     else
/* gradient direction approaches to the vertical direction.		*/
      if ((*diy> *tii1) && (*diy>= *tii2))
       *t1= *diy;
       else
        *t1= 0;
    }
 }
return;
}

     


/***********************************************************************
*
*  Routine Name: gef()
*
*          Date:
*        
*       Purpose: This subprogram calculate 1st derivative. 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:
*
***********************************************************************/

void gef(in,dx,dy,mx,my,dim,a0)
unsigned char *in,*dx,*dy;
int mx,my,dim;
float a0;
/* *in: input image; *dx,*dy: output images;			        */
/* mx: width of the images; my: height of the image; 			*/
/* a0: filter parameter.						*/ 
{
unsigned char *tti,*ti,*t1,*dd;
int i,j,k,ki,mxx,myy,mxxx,myyy,v,v1;


mxx=mx-1; myy=my-1; mxxx=mxx-1; myyy=myy-1;
dd=(unsigned char*)malloc(dim*sizeof(unsigned char));
/*----------------------------------------------------------------------*/
/* 			Horizontal Processing 				*/
/*----------------------------------------------------------------------*/
/*	 vertical filtering		*/
for (j=0; j<mx; j++) /* up to down, input in *in, output in *dx */
 {
  *(dx+j)= *(in+j);
  for (k=j+mx; k<dim; k+=mx)
   {
    v= *(dx+k-mx); v1= *(in+k);
    *(dx+k)=v+a0*(v1-v)+0.5;
   }
 }
ki=myyy*mx;
for (j=ki; j<ki+mx; j++) /* down to up, input in *dx, output in *dx */
 for (k=j; k>=0; k-=mx)
  {
   v= *(dx+k+mx); v1= *(dx+k);
   *(dx+k)=v+a0*(v1-v)+0.5;
  }
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)/* left to right, input in *dx, output in *dd */
 {
  *(dd+i)= *(dx+i);
  for (k=i+1; k<i+mx; k++)
   {
    v= *(dd+k-1); v1= *(dx+k);
    *(dd+k)=v+a0*(v1-v)+0.5;
   }
 }
for (i=0; i<dim; i+=mx)/* right to left, input in *dx, output in *dy */
 {
  *(dy+i+mxx)= *(dx+i+mxx);
  for (k=i+mxxx; k>=i; k--)
   {
    v= *(dy+k+1); v1= *(dx+k);
    *(dy+k)=v+a0*(v1-v)+0.5;
   }
 }
/* calculate 1st derivative in horiz. direction from *dd,*dx and *dy.	*/
/* result in *dx ( 1st derivative ) 			  	     	*/
ti=dd; t1=dx;
for (tti=dy; tti<dy+dim-1; tti++)
  {
   v= *tti; v1= *ti++;
   *t1++ =abs(v-v1);
   }
/*----------------------------------------------------------------------*/
/* 			Vertical Processing 				*/
/*----------------------------------------------------------------------*/
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)/* left to right, input in *in, output in *dy */
 {
  *(dy+i)= *(in+i);
  for (k=i+1; k<i+mx; k++)
   {
    v= *(dy+k-1); v1= *(in+k);
    *(dy+k)=v+a0*(v1-v)+0.5;
   }
 }
for (i=0; i<dim; i+=mx)/* right to left, input in *dy, output in *dy */
 for (k=i+mxxx; k>=i; k--)
  {
   v= *(dy+k+1); v1= *(dy+k);
   *(dy+k)=v+a0*(v1-v)+0.5;
  }
/*	 vertical filtering		*/
for (j=0; j<mx; j++)/* up to down, input in *dy, output in *dd */
 {
  *(dd+j)= *(dy+j); 
  for (k=j+mx; k<dim; k+=mx)
   {
    v= *(dd+k-mx); v1= *(dy+k);
    *(dd+k)=v+a0*(v1-v)+0.5;
   }
 }
ki=myy*mx;
for (j=ki; j<ki+mx; j++)/* down to up, input in *dy, output in *in */
 {
  *(in+j)= *(dy+j);
  for (k=j-mx; k>=0; k-=mx)
   {
    v= *(in+k+mx); v1= *(dy+k);
    *(in+k)=v+a0*(v1-v)+0.5;
   }
 }
/* calculate 1st derivative in vertical direction  from  	*/
/* *dd,*dy and *in. result in *dy ( 1st derivative ) 	        */
ti=dd; t1=dy;
for (tti=in; tti<in+dim-1; tti++)
  {
   v= *tti; v1= *ti++;
   *t1++ =abs(v-v1);
   }
free(dd);
return;
}


/***********************************************************************
*
*  Routine Name: gradadp()
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  
*
* Modifications:
*
***********************************************************************/

void gradadp(mx,my,dim,w,o,l,out)
int mx,my,dim,w;
unsigned char *o,*l,*out;
{
unsigned char *l0,*l1,*l2,*l3,*l4,*l5,*l6,*l7,*l8;
int s,ss,i,ii,j,jj,ki,k,ww;
float ms,mss;

   ww=w/2;
   for (i=1; i<my-1; i++)
   {
      ki=i*mx;

      for (j=1; j<mx-1; j++)
      {
         k=ki+j;
         l0=l+k; 
         if (*l0>0)
         {
            l2=l0-mx; 
            l1=l2-1; 
            l3=l2+1;
            l4=l0-1; 
            l5=l0+1;
            l7=l0+mx; 
            l6=l7-1; 
            l8=l7+1;

            if ( (*l1==0) || (*l2==0) || (*l3==0) || (*l4==0) ||
                 (*l5==0) || (*l6==0) || (*l7==0) || (*l8==0)   ) *(out+k)=1;

            else *(out+k)=0;
         }
         else *(out+k)=0;
      }
   }

   for (i= (ww*mx+ww); i< ((ww+1)*mx-ww); i++)
   {
      for (j=i; j<(i+mx*(my-w)); j+=mx)
      {
         if (*(out+j)>0)
         {
            s=0; ss=0; ms=0; mss=0;

            for (ii= (j-ww*mx-ww); ii< (j-ww*mx+ww); ii++)
            {
               for (jj=ii; jj<(ii+w*mx); jj+=mx)
               {
                  if (*(l+jj)>0)
                  { 
                     ms+= *(o+jj); 
                     s+=1;
                  }
                  else
                  {
           	     mss+= *(o+jj); 
                     ss+=1;
        	  }	

		  /*
		   *  Changed "*(out+j)=fabs(ms/s - mss/ss);" to the following
		   *  code since initially s or ss would equal to zero.
		   *		Mark & Charlie (2/15/91)
		   */
                  if (s == 0)
		     *(out+j)=fabs(mss/ss);
                  else if (ss == 0)
		     *(out+j)=fabs(ms/s);
                  else
		     *(out+j)=fabs(ms/s - mss/ss);
               }
            }
         }
      }
   }
   return;
}

/***********************************************************************
*
*  Routine Name: sdefextr()
*
*          Date:
*        
*       Purpose: This subprogram detects the step edges by combination
*		 of z-c  of  2nd derivative and the sign correspondence
*		 between 1st and 2nd derivatives. 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:
*
***********************************************************************/

void sdefextr(dx,dy,ddx,ddy,out,mx,my,s)
unsigned char *dx,*dy,*ddx,*ddy,*out;
int mx,my;
float s;
/* *dx,*dy,*ddx,*ddy: input image; *out: output images;			*/
/* mx: width of the images; my: height of the image; 			*/
/* s: Threshold for choosing the gradient direction.			*/ 
{
int i,j,ki,k;
unsigned char *dix,*ddix,*diy,*ddiy,*tii,*t1,*ti;
for (i=1; i<my; i++)
 {
  ki=i*mx;
  for (j=1; j<mx; j++)
   {
    k=ki+j;
    ddix=ddx+k; dix=dx+k; ddiy=ddy+k; diy=dy+k; t1=out+k;
    ti=ddix-1;
    tii=ddiy-mx;
    if (*dix>s*(*diy))
/* gradient direction approaches to the horizontal direction.		*/
     if (((*ddix==2) && (*ti<2)) || ((*ti>1) && (*ddix==1)))
      *t1= *dix;
      else
       *t1= 0;
     else
      if (*diy>s*(*dix))
/* gradient direction approaches to the vertical direction.		*/
       if (((*ddiy==2) && (*tii<2)) || ((*tii>1) && (*ddiy==1)))
        *t1= *diy;
        else
         *t1= 0;
       else
        {
/* gradient direction is not in the horizontal direction, either 	*/
/* in the vertical direction.						*/
         if (((*ddix==2) && (*ti<2)) || ((*ti>1) && (*ddix==1)))
          *t1= *dix;
          else
           *t1= 0;
	 if ((*t1==0) || (*diy> *dix))
          if (((*ddiy==2) && (*tii<2)) || ((*tii>1) && (*ddiy==1)))
           *t1= *diy;
        }
    }
 }
return;
}

     
/***********************************************************************
*
*  Routine Name: sdeff()
*
*          Date:
*        
*       Purpose: This subprogram calculate 1st derivative and zero-crossing
*		 of 2nd derivative by the differential operators of the
*		 exponential filter 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:
*
***********************************************************************/

void sdeff(in,dx,dy,ddx,ddy,mx,my,dim,a0)
unsigned char *in,*dx,*dy,*ddx,*ddy;
int mx,my,dim;
float a0;
/* *in: input image; *dx,*dy,*ddx,*ddy: output images;			*/
/* mx: width of the images; my: height of the image; 			*/
/* a0: filter parameter.						*/ 
{
unsigned char *tti,*ti,*t1;
int i,j,k,ki,mxx,myy,mxxx,myyy,v,v1,v2;
mxx=mx-1; myy=my-1; mxxx=mxx-1; myyy=myy-1;
/*----------------------------------------------------------------------*/
/* 			Horizontal Processing 				*/
/*----------------------------------------------------------------------*/
/*	 vertical filtering		*/
for (j=0; j<mx; j++)	  /* up to down, input in *in, output in *dx */
 {
  *(dx+j)= *(in+j);
  for (k=j+mx; k<dim; k+=mx)
   {
    v= *(dx+k-mx); v1= *(in+k);
    *(dx+k)=v+a0*(v1-v)+0.5;
   }
 }
ki=myyy*mx;
for (j=ki; j<ki+mx; j++)/* down to up, input in *dx, output in *dx */
 for (k=j; k>=0; k-=mx)
  {
   v= *(dx+k+mx); v1= *(dx+k);
   *(dx+k)=v+a0*(v1-v)+0.5;
  }
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)/* left to right, input in *dx, output in *ddx */
 {
  *(ddx+i)= *(dx+i);
  for (k=i+1; k<i+mx; k++)
   {
    v= *(ddx+k-1); v1= *(dx+k);
    *(ddx+k)=v+a0*(v1-v)+0.5;
   }
 }
for (i=0; i<dim; i+=mx)/* right to left, input in *dx, output in *dy */
 {
  *(dy+i+mxx)= *(dx+i+mxx);
  for (k=i+mxxx; k>=i; k--)
   {
    v= *(dy+k+1); v1= *(dx+k);
    *(dy+k)=v+a0*(v1-v)+0.5;
   }
 }
/* calculate 1st deriv. and z-c of 2nd deriv. in horiz. direction    */
/* from *ddx,*dx and *dy. result in *dx ( 1st deriv. ) and *ddx ( the*/
/* sign of 1st deriv. and 2nd deriv. ).				     */
ti=ddx; t1=dx;
for (tti=dy; tti<dy+dim-1; tti++)
  {
   v= *tti; v1= *ti; v2= *t1;
   if ((v+v1-v2-v2)>=0)
    if ((v2=v-v1)>=0)
     {
      *ti++ =3; *t1++ =v2;
     }
     else
      {
       *ti++ =2; *t1++ = -v2;
      }
    else
     if ((v2=v-v1)>=0)
      {
       *ti++ =1; *t1++ =v2;
      }
      else
       {
        *ti++ =0; *t1++ = -v2;
       }
   }
/*----------------------------------------------------------------------*/
/* 			Vertical Processing 				*/
/*----------------------------------------------------------------------*/
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)/* left to right, input in *in, output in *dy */
 {
  *(dy+i)= *(in+i);
  for (k=i+1; k<i+mx; k++)
   {
    v= *(dy+k-1); v1= *(in+k);
    *(dy+k)=v+a0*(v1-v)+0.5;
   }
 }
for (i=0; i<dim; i+=mx)/* right to left, input in *dy, output in *dy */
 for (k=i+mxxx; k>=i; k--)
  {
   v= *(dy+k+1); v1= *(dy+k);
   *(dy+k)=v+a0*(v1-v)+0.5;
  }
/*	 vertical filtering		*/
for (j=0; j<mx; j++) /* up to down, input in *dy, output in *ddy */
 {
  *(ddy+j)= *(dy+j); 
  for (k=j+mx; k<dim; k+=mx)
   {
    v= *(ddy+k-mx); v1= *(dy+k);
    *(ddy+k)=v+a0*(v1-v)+0.5;
   }
 }
ki=myy*mx;
for (j=ki; j<ki+mx; j++)/* down to up, input in *dy, output in *in */
 {
  *(in+j)= *(dy+j);
  for (k=j-mx; k>=0; k-=mx)
   {
    v= *(in+k+mx); v1= *(dy+k);
    *(in+k)=v+a0*(v1-v)+0.5;
   }
 }
/* calculate 1st deriv. and z-c of 2nd deriv. in vertical direction    */
/* from *ddy,*dy and *in. result in *dy ( 1st deriv. ) and *ddy ( the  */
/* sign of 1st deriv. and 2nd deriv. ).				       */
ti=ddy; t1=dy;
for (tti=in; tti<in+dim-1; tti++)
  {
   v= *tti; v1= *ti; v2= *t1;
   if ((v+v1-v2-v2)>=0)
    if ((v2=v-v1)>=0)
     {
      *ti++ =3; *t1++ =v2;
     }
     else
      {
       *ti++ =2; *t1++ = -v2;
      }
    else
     if ((v2=v-v1)>=0)
      {
       *ti++ =1; *t1++ =v2;
      }
      else
       {
        *ti++ =0; *t1++ = -v2;
       }
   }
return;
}


/***********************************************************************
*
*  Routine Name: threshold()
*
*          Date:
*        
*       Purpose: This subprogram is used for an image segmentation by the
*		 hysteretic threshold.  
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:  Looked like this routine was converted from fortran code to
*	C code.  threshold() was crashing when passing in row sizes "mx"
*	greater than 513.  Changed fixed static arrays int h[513][3],
*	l[513][5] to use typedef's high and low, and change h & l to
*	"high *h; low *l;".  This allowed the routine to allocate the arrays
*	according to the actual deisred "mx" row size.  
*
*	Also changed all h[][] & l[][] indexes to be one less on the col
*	index.  ie) l[k][5] was changed to l[k][4].  The "h" array was
*	being indexed from h[][1-3] instead of h[][0-2] and te "l" array
*	being indexed l[][1-5] instead of l[][0-4].
*
***********************************************************************/

typedef int high[3];
typedef int low[5];

void threshold(in,out,mx,my,dim,s1,s2,min)
unsigned char *in,*out;
int mx,my,dim,min,s1,s2;
/* *in: input image; *out: output image;				*/
/* mx: image width; my: image height; dim: image dimensin (my x mx).	*/
{
unsigned char *t1,*ti,*ti1;
int c,s;
struct dlink
{
int i,j,w;
struct dlink *next;
};
struct mlink
{
int n,n0,l,k;
struct dlink *dl,*dl1;
};
struct mlink *p,*pp,**ml;
struct dlink *d;
int nh,nl,i,j,jj,k,kk,ki,ml_size;
high *h;
low  *l;
/*-----------------------------------------------------------------------*/
/* detect the connected regions in which every pixel level is superior in*/
/* the threshold s1.							 */
/*-----------------------------------------------------------------------*/	
h = (high *) calloc(mx,sizeof(high));
l = (low *) calloc(mx,sizeof(low));
ml_size = mx;
ml = (struct mlink **) malloc(ml_size * sizeof(struct mlink *));

c=0; nh=0; nl=0;
for (i=1; i<my-1; i++)
 {
  ki=i*mx; *(in+ki)=0; *(in+ki+mx-1)=0;
  for (j=1; j<mx-1; j++)
   {
    ti=in+ki+j;
    if (*ti>s1)
     {
      ti1=ti-1;
      if (*ti1<=s1)
       l[++nl][0]=j;
      ti1=ti+1;
      if (*ti1<=s1)
       l[nl][1]=j;
     };
   };
  if (nl>0)
   {
    if (nh>0)
     {
      s=1;
      for (k=1; k<=nl; k++)
       {
        while ((l[k][1]>=(h[s][0]-1)) && (s<=nh))
         {
          if (h[s][1]>=(l[k][0]-1))
           if (l[k][2]==0)
            l[k][2]=s;
	  s++;
         };
        if (l[k][2]>0)
         l[k][3]= --s;
       };
     };
    l[nl+1][2]=1000000; l[0][3]=1000000;
    for (k=1; k<=nl; k++)
     if (l[k][2]==0)
      {
       l[k][4]= ++c;
       if (c >= ml_size)
       {
          ml_size += mx;
	  ml = (struct mlink **) realloc(ml, ml_size * sizeof(struct mlink *));
	  if (ml == NULL) return;
       }
       p=(struct mlink*)malloc(sizeof(struct mlink));
       ml[c]=p; (*p).k=0; (*p).l=1;
       (*p).n=l[k][1]-l[k][0]+1; (*p).n0=c;
       d=(struct dlink*)malloc(sizeof(struct dlink));
       (*d).i=i; (*d).j=l[k][0]; (*d).w=(*p).n; (*d).next=NULL;
       (*p).dl=d; (*p).dl1=d;
      }
      else
       if (l[k-1][3]==l[k][2])
        {
         p=ml[h[l[k-1][2]][2]];
         while ((*p).k!=0)
          p=ml[(*p).k];
         if (l[k][2]!=l[k][3])
          for (jj=l[k][2]+1; jj<=l[k][3]; jj++)
           {
	    pp=ml[h[jj][2]];
	    while ((*pp).k!=0)
	     pp=ml[(*pp).k];
	    if (pp!=p)
	     {
	      (*((*pp).dl1)).next=(*p).dl;
	      (*p).dl=(*pp).dl;
 	      (*p).n=(*p).n+(*pp).n; (*p).l=(*p).l+(*pp).l; (*pp).k=(*p).n0;
	      (*pp).n=0;
	     };
           };
         if (l[k][3]==l[k+1][2])
          {
           (*p).n=(*p).n+l[k+1][1]-l[k+1][0]+1;
           d=(struct dlink*)malloc(sizeof(struct dlink));
           (*d).i=i; (*d).j=l[k+1][0]; (*d).w=l[k+1][1]-l[k+1][0]+1;
           ((*p).l)++; (*d).next=(*p).dl; (*p).dl=d; l[k+1][4]=(*p).n0;
          };
        }
        else
         {
          p=ml[h[l[k][2]][2]];
          while ((*p).k!=0)
          p=ml[(*p).k];
          (*p).n=(*p).n+l[k][1]-l[k][0]+1;
          d=(struct dlink*)malloc(sizeof(struct dlink));
          (*d).i=i; (*d).j=l[k][0];
          (*d).w=l[k][1]-l[k][0]+1;
          (*d).next=(*p).dl; (*p).dl=d;
          l[k][4]=(*p).n0; ((*p).l)++;
          if (l[k][2]!=l[k][3])
           for (jj=l[k][2]+1; jj<=l[k][3]; jj++)
            {
	     pp=ml[h[jj][2]];
	     while ((*pp).k!=0)
	     pp=ml[(*pp).k];
	     if (pp!=p)
	      {
	       (*((*pp).dl1)).next=(*p).dl;
	       (*p).dl=(*pp).dl;
 	       (*p).n=(*p).n+(*pp).n; (*p).l=(*p).l+(*pp).l; (*pp).k=(*p).n0;
	       (*pp).n=0;
	      };
            };
          if (l[k][3]==l[k+1][2])
           {
	    (*p).n=(*p).n+l[k+1][1]-l[k+1][0]+1;
	    d=(struct dlink*)malloc(sizeof(struct dlink));
	    (*d).i=i; (*d).j=l[k+1][0]; (*d).w=l[k+1][1]-l[k+1][0]+1;
	    ((*p).l)++; (*d).next=(*p).dl; (*p).dl=d; l[k+1][4]=(*p).n0;
           };
         };
    for (k=1; k<=nl; k++)
     {
      h[k][0]=l[k][0]; h[k][1]=l[k][1]; h[k][2]=l[k][4];
      l[k][2]=0; l[k][3]=0;
     };
    l[nl+1][2]=0;
   };
  nh=nl; nl=0;
 };
/*----------------------------------------------------------------*/
/* write the segments to an image, in every segment, at least a	  */
/* pixel level is superior in the threshold s1.			  */
/*----------------------------------------------------------------*/
ti=out;
for (i=0; i<dim-1; i++)
  *ti++ =0;
s=0;
for (i=1; i<=c; i++)
 {
  p=ml[i];
  if ((*p).n>min)
   {
    d=(*p).dl; 
    for (j=1; j<=(*p).l; j++)
     {
      for (k=(*d).j; k<=(*d).j+(*d).w-1; k++)
       {
	t1=in+((*d).i)*mx+k;
        if (*t1>s2)
         {
          s++; d=(*p).dl;
          for (jj=1; jj<=(*p).l; jj++)
           {
	    for (kk=(*d).j; kk<=(*d).j+(*d).w-1; kk++)
      	     {
              ti=out+((*d).i)*mx+kk; *ti=255;
	     }
            d=(*d).next;
	   }
    	  goto nn;
         }
       }
      d=(*d).next;
     }
   }
nn: s=s;
  }
free(h); free(l); free(ml);
return;
}


/***********************************************************************
*
*  Routine Name: 
*
*          Date:
*        
*       Purpose: This subprogram calculate 1st derivative.
*		 (This routine is equivalent to gef(), From what I
*		  understood, it is used for processing images where
*		  the contrast is weak)
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:
*
***********************************************************************/

void geff(in,dx,dy,mx,my,dim,a0)
unsigned char *in,*dx,*dy;
int mx,my,dim;
float a0;
/* *in: input image; *dx,*dy: output images;			        */
/* mx: width of the images; my: height of the image; 			*/
/* a0: filter parameter.						*/ 
{
float a1,a00,c1,c2,u;
unsigned char *t1;
int i,j,k,ki,mxx,myy,mxxx,myyy,mxy,v1,*bu1,*bu2,*b1,*b2;
mxx=mx-1; myy=my-1; mxxx=mxx-1; myyy=myy-1;
a1=1-a0; a00=a0/a1; c1=1/a1; /* c2=2/a1; */
if (mx>my)
 mxy=mx+1;
 else
  mxy=my+1;
bu1=(int*)calloc(mxy,sizeof(int));
bu2=(int*)calloc(mxy,sizeof(int));
/*----------------------------------------------------------------------*/
/* 			Horizontal Processing 				*/
/*----------------------------------------------------------------------*/
/*	 vertical filtering		*/
for (j=0; j<mx; j++)	  /* up to down, input in *in, output in *dx */
 {
  u= *(in+j);
  for (k=j+mx; k<dim; k+=mx)
   {
    *(dx+k-mx)=u; v1= *(in+k);
    u=u+a0*(v1-u);
   }
 }
ki=myyy*mx;
for (j=ki; j<ki+mx; j++)/* down to up, input in *dx, output in *dx */
 {
  u= *(dx+j);
  for (k=j; k>=0; k-=mx)
   {
    *(dx+k+mx)=u; v1= *(dx+k);
    u=u+a0*(v1-u);
   }
 }
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)
 {
  b1=bu1; u= *(dx+i)*c1;/* left to right, input in *dx, output in *bu1 */
  for (k=i+1; k<i+mx; k++)
   {
    *b1++ =u; v1= *(dx+k);
    u=a1*u+a00*v1;
   }
  b2=bu2+mxx; u= *(dx+i+mxx)*c1;/* right to left, input in *dx, output in *bu2 */
  for (k=i+mxxx; k>=i; k--)
   {
    *b2-- =u; v1= *(dx+k);
    u=a1*u+a00*v1;
   }
/* calculate 1st deriv. and z-c of 2nd deriv. in horiz. direction        */
/* from *bu1, *bu2 and *dx. result in *dx ( 1st deriv. ) and *ddx ( the  */
/* sign of 1st deriv. and 2nd deriv. ).				         */
  b1=bu1; b2=bu2;
  for (t1=dx+i; t1<dx+i+mx; t1++)
   *t1=abs(*b2++ - *b1++);
 }
/*----------------------------------------------------------------------*/
/* 			Vertical Processing 				*/
/*----------------------------------------------------------------------*/
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)/* left to right, input in *in, output in *dy */
 {
  u= *(in+i);
  for (k=i+1; k<i+mx; k++)
   {
    *(dy+k-1)=u; v1= *(in+k);
    u=u+a0*(v1-u);
   }
 }
for (i=0; i<dim; i+=mx)/* right to left, input in *dy, output in *dy */
 {
  u= *(dy+i+mxx);
  for (k=i+mxxx; k>=i; k--)
   {
    *(dy+k+1)=u; v1= *(dy+k);
    u=u+a0*(v1-u);
   }
 }
/*	 vertical filtering		*/
ki=myy*mx;
for (j=0; j<mx; j++) 
 {
  b1=bu1; u= *(dy+j)*c1; /* up to down, input in *dy, output in *bu1 */
  for (k=j+mx; k<dim; k+=mx)
   {
    *b1++ =u; v1= *(dy+k);
    u=a1*u+a00*v1;
   }
  b2=bu2+myy; u= *(dy+ki+j)*c1;/* down to up, input in *dy, output in *bu2 */
  for (k=ki+j-mx; k>=0; k-=mx)
   {
    *b2-- =u; v1= *(dy+k);
    u=a1*u+a00*v1;
   }
/* calculate 1st deriv. and z-c of 2nd deriv. in vertical direction      */
/* from *bu1, *bu2 and *dy. result in *dy ( 1st deriv. ) and *ddy ( the  */
/* sign of 1st deriv. and 2nd deriv. ).				         */
  b1=bu1; b2=bu2;
  for (t1=dy+j; t1<dy+dim; t1+=mx)
   *t1=abs(*b2++ - *b1++);
 }
free(bu1);
free(bu2);
return;
}



/***********************************************************************
*
*  Routine Name: sdefff()
*
*          Date:
*        
*       Purpose: This subprogram calculate 1st derivative and zero-crossing
*		 of 2nd derivative by the differential operators of the
*		 exponential filter 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jian Zhao 
*
* Modifications:
*
***********************************************************************/

void sdefff(in,dx,dy,ddx,ddy,mx,my,dim,a0)
unsigned char *in,*dx,*dy,*ddx,*ddy;
int mx,my,dim;
float a0;
/* *in: input image; *dx,*dy,*ddx,*ddy: output images;			*/
/* mx: width of the images; my: height of the image; 			*/
/* a0: filter parameter.						*/ 
{
float a1,a00,c1,c2,u;
unsigned char *ti,*t1;
int i,j,k,ki,mxx,myy,mxxx,myyy,mxy,v,v1,v2,*bu1,*bu2,*b1,*b2;
mxx=mx-1; myy=my-1; mxxx=mxx-1; myyy=myy-1;
a1=1-a0; a00=a0/a1; c1=1/a1; c2=2/a1;
if (mx>my)
 mxy=mx;
 else
  mxy=my;
bu1=(int*)calloc(mxy,sizeof(int));
bu2=(int*)calloc(mxy,sizeof(int));
/*----------------------------------------------------------------------*/
/* 			Horizontal Processing 				*/
/*----------------------------------------------------------------------*/
/*	 vertical filtering		*/
for (j=0; j<mx; j++)	  /* up to down, input in *in, output in *dx */
 {
  u= *(in+j);
  for (k=j+mx; k<dim; k+=mx)
   {
    *(dx+k-mx)=u; v1= *(in+k);
    u=u+a0*(v1-u);
   }
 }
ki=myyy*mx;
for (j=ki; j<ki+mx; j++)/* down to up, input in *dx, output in *dx */
 {
  u= *(dx+j+mx);
  for (k=j; k>=0; k-=mx)
   {
    *(dx+k+mx)=u; v1= *(dx+k);
    u=u+a0*(v1-u);
   }
 }
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)
 {
  b1=bu1; u= *(dx+i)*c1;/* left to right, input in *dx, output in *bu1 */
  for (k=i+1; k<i+mx; k++)
   {
    *b1++ =u; v1= *(dx+k);
    u=a1*u+a00*v1;
   }
  b2=bu2+mxx; u= *(dx+i+mxx)*c1;/* right to left, input in *dx, output in *bu2 */
  for (k=i+mxxx; k>=i; k--)
   {
    *b2-- =u; v1= *(dx+k);
    u=a1*u+a00*v1;
   }
/* calculate 1st deriv. and z-c of 2nd deriv. in horiz. direction        */
/* from *bu1, *bu2 and *dx. result in *dx ( 1st deriv. ) and *ddx ( the  */
/* sign of 1st deriv. and 2nd deriv. ).				         */
  b2=bu2; t1=dx+i; ti=ddx+i;
  for (b1=bu1; b1<bu1+mx; b1++)
   {
    v= *b2++; v1= *b1; v2= *t1*c2;
    if ((v+v1-v2)>=0)
     if ((v2=v-v1)>=0)
      {
       *ti++ =3; *t1++ =v2;
      }
     else
      {
       *ti++ =2; *t1++ = -v2;
      }
    else
     if ((v2=v-v1)>=0)
      {
       *ti++ =1; *t1++ =v2;
      }
      else
       {
        *ti++ =0; *t1++ = -v2;
       }
   }
 }
/*----------------------------------------------------------------------*/
/* 			Vertical Processing 				*/
/*----------------------------------------------------------------------*/
/*	 horizontal filtering		*/
for (i=0; i<dim; i+=mx)/* left to right, input in *in, output in *dy */
 {
  u= *(in+i);
  for (k=i+1; k<i+mx; k++)
   {
    *(dy+k-1)=u; v1= *(in+k);
    u=u+a0*(v1-u);
   }
 }
for (i=0; i<dim; i+=mx)/* right to left, input in *dy, output in *dy */
 {
  u= *(dy+i+mxx);
  for (k=i+mxxx; k>=i; k--)
   {
    *(dy+k+1)=u; v1= *(dy+k);
    u=u+a0*(v1-u);
   }
 }
/*	 vertical filtering		*/
ki=myy*mx;
for (j=0; j<mx; j++) 
 {
  b1=bu1; u= *(dy+j)*c1; /* up to down, input in *dy, output in *bu1 */
  for (k=j+mx; k<dim; k+=mx)
   {
    *b1++ =u; v1= *(dy+k);
    u=a1*u+a00*v1;
   }
  b2=bu2+myy; u= *(dy+ki+j)*c1;/* down to up, input in *dy, output in *bu2 */
  for (k=ki+j-mx; k>=0; k-=mx)
   {
    *b2-- =u; v1= *(dy+k);
    u=a1*u+a00*v1;
   }
/* calculate 1st deriv. and z-c of 2nd deriv. in vertical direction      */
/* from *bu1, *bu2 and *dy. result in *dy ( 1st deriv. ) and *ddy ( the  */
/* sign of 1st deriv. and 2nd deriv. ).				         */
  ti=ddy+j; t1=dy+j; b2=bu2;
  for (b1=bu1; b1<bu1+my; b1++)
   {
    v= *b2++; v1= *b1; v2= *t1*c2;
    if ((v+v1-v2)>=0)
     if ((v2=v-v1)>=0)
      {
       *ti=3; *t1=v2; ti+=mx; t1+=mx;
      }
      else
       {
        *ti=2; *t1= -v2; ti+=mx; t1+=mx;
       }
     else
      if ((v2=v-v1)>=0)
       {
        *ti=1; *t1=v2; ti+=mx; t1+=mx;
       }
       else
        {
         *ti=0; *t1= -v2; ti+=mx; t1+=mx;
        }
   }
 }
free(bu1);
free(bu2);
return;
}
