float convert(char i){  if (i>=0) return (float)i;  else return((float)(i+256));}float fconv(register float *x, register float *y, register int n){  float sum=0;    if (n>0) {    do {      sum+= *x++ * *y;      y++;    } while (--n>0);  }  return(sum);}/* This function takes a pointer to a float array as argument and   scales the array to be less than 255 */void scaler(float **f_array, int vert_size, int horiz_size){  float scale, max=0;  register row, col;    for(row=0; row < vert_size; row++)    for(col=0; col < horiz_size; col++) {      /*	printf("%e\n",f_array[row][col]);  */      if (f_array[row][col] > max)	max = f_array[row][col];    }    scale = max / 255.0;  if (scale == 0) {    fprintf(stderr, "Div by zero\n");    exit(1); }    for(row=0; row < vert_size; row++)    for(col=0; col < horiz_size; col++)       f_array[row][col] /= scale;}void magedges(float **imgarray, char **edgearray, int vert_size, int horiz_size){  register int row, col;    for(row = 0; row < vert_size; row++)    for(col = 0; col < horiz_size; col++)      if (edgearray[row][col] == 0)	imgarray[row][col] = 0;    scaler(imgarray, vert_size, horiz_size);    /* scale to 255 */}/*===========================================================================*//* Arguments: gradient magnitude array (float), map of local maximums   (char), value for high threshold (float), number of passes for   'growing', and two arrays for workspace.    The function first creates an array of points above the low threshold    (thresfactor times the high threshold).   Then the map of local maximums is run through the high threshold -THIS   IS A MUTATION. Next, for every point above the high threshold, the   eight nearest neighbors are checked to see if they are above the low   threshold. If so, they are marked (and the edge array is mutated again).   Thus, the edges are 'grown'. The number of passes defines how many   times this growing routine is run. */void threshold(float **mag, char **edges, 	  float hthreshold, float thresfactor, int passes, 	  char **low, char **wk_space,	  int vert_size, int horiz_size){  register int row, col;  float lthreshold;  int n;    lthreshold = hthreshold * thresfactor;    /* create an array of points above low threshold */    for(row = 0; row < vert_size; row++)    for(col = 0; col < horiz_size; col++) {            if ((edges[row][col]) && (mag[row][col] > lthreshold)) 	low[row][col] = 255; /* edges[row][col] */      /* indicates local maximum   */	      else 	low[row][col] = 0;    }	    /* now I must create an array of pts above the high threshold */  /* do this by changing edges */  for(row = 0; row < vert_size; row++)    for(col = 0; col < horiz_size; col++) {            if ((edges[row][col]) && (mag[row][col] > hthreshold))	edges[row][col] = 255;            else	edges[row][col] = 0;    }    /* now I have low and high pts */  for(n=0; n < passes; n++)     grow(edges, low, wk_space, vert_size, horiz_size);}void grow(char **r_array, char **chk_array, char **w_array,     int vert_size, int horiz_size){  register int row, col, i, j;    /* copy edge array of high points to workspace */    for(row=0; row < vert_size; row++)    for(col=0; col < horiz_size; col++)      w_array[row][col] = r_array[row][col];    /* now to add more points appropriately */    for(row=1; row < vert_size - 1; row++)    for(col=1; col < horiz_size - 1; col++)             if(w_array[row][col])  /** if high point **/		for(i = -1; i <= 1; i++)	  for(j = -1; j <= 1; j++)	    	    r_array[(row + i)][(col + j)] =	      chk_array[(row + i)][(col + j)];	}/*===========================================================================*//* This function takes as arguments: x and y derivative arrays (float),   histogram percentile, and a work space array (float). The function returns   an estimate of the maximum noise level in the picture.*/float find_thres(float **xdir, float **ydir, float percentile,		 float **h_val, int vert_size, int horiz_size){  int histogram[HISTOGRAM_SIZE];  float x2dx, y2dy, division_size, hthreshold, max;  int bucket, sumpts, total_points, ppt;  register int row, col;    /* initialize histogram */  for(bucket = 0; bucket < HISTOGRAM_SIZE; bucket++)    histogram[bucket] = 0;    /* run 2nd derivative operator in x and y directions */  /* also note the maximum value */  /* remember that the edges are garbage */  for(row = 1, max = 0; row < vert_size - 2; row ++)    for(col = 1; col < horiz_size - 2; col++) {      x2dx = xdir[row][col] * (-2) + xdir[row][(col - 1)]	+ xdir[row][(col + 1)];      y2dy = ydir[row][col] * (-2) + ydir[(row - 1)][col]	+ ydir[(row + 1)][col];            /* square 2nd derivate and store in work space */            h_val[row][col] = x2dx * x2dx + y2dy * y2dy;            if(h_val[row][col] > max)	max = h_val[row][col];    }    /* fill buckets */  for(row = 1; row < vert_size - 2 ; row++)    for(col = 1; col < horiz_size - 2 ; col++) {            if(h_val[row][col] == max)	bucket = HISTOGRAM_SIZE - 1;      else	bucket = (int)(floor((h_val[row][col] / max) *			     HISTOGRAM_SIZE));      (histogram[bucket])++;    }    total_points = (horiz_size - 3) * (vert_size - 3);    division_size = max / HISTOGRAM_SIZE;    ppt = (percentile / 100) * total_points; /* percent to pts  */    for(bucket = 0, sumpts = 0; sumpts < ppt; bucket++)    sumpts += histogram[bucket];    /* calculate the value at percentile and take sqrt */    hthreshold = ((bucket + bucket - 1) / 2) * division_size;    hthreshold = sqrt(hthreshold) * MAGIC_SCALE_FACTOR;    return(hthreshold);}/*===========================================================================*//* This function takes as arguments an array to operate on, and threearrays to fill with results.  The function differentiates the sourcearray in the x and y directions (over a four element window) and stores these results in two of the destination arrays (float).The function also computes the magnitude of the gradient at each point and returns this result in the source array. NOTE: THIS MUTATES THE SOURCE ARRAY.The map of the local maximums of the gradient magnitude array is storedin the third array (character array).Note: Realize that the last row and column of the x derivative, yderivative, and magnitude arrays contain garbage. Sizes of arrays are defined in the include file gdef.h. */void gradmax(float **img_array, float **xdir, float **ydir, char **edges,	int vert_size, int horiz_size){  register int row, col;  int x0y0, x0y1, x1y0, x1y1;  float front, back;	/* temporary storage for max check */  float mag2;    /* start gradient computation */    for(row=0; row < vert_size-1; row++) /* rt and bt edge */    for(col=0; col < horiz_size-1; col++) {            x0y1 = img_array[row][col];      x1y1 = img_array[row][(col + 1)];      x0y0 = img_array[(row + 1)][col];      x1y0 = img_array[(row + 1)][(col + 1)];            xdir[row][col]= x1y1 + x1y0 - x0y0 - x0y1;      ydir[row][col]= x0y1 + x1y1 - x0y0 - x1y0;            /*   compute magnitude and store in original */            mag2 = xdir[row][col] * xdir[row][col] +	ydir[row][col] * ydir[row][col];            img_array[row][col] = sqrt(mag2);		          }      /* Points in the gradient magnitude array are tested to see if they are     local maximums.     Definitions:        front--nearest pixel in the direction defined by the gradient        back---nearest pixel in the direction opposite that defined by               the gradient     If the gradient magnitude at a point is greater than that at the front     and greater than that at the back, it is marked as local maximum.      If the gradient does not point directly at a pixel, linear interpolation     is performed. Note: this is called 'non-maximum suppression' in the      literature. */  /* non-max suppression */  /* remember that edges of arrays are garbage */  for(row = 1; row < vert_size - 2; row++)    for(col = 1; col < horiz_size - 2; col++) {            if(xd==0 && yd==0)	edges[row][col]=0;            /* quadrants I and III */            else if ((xd >= 0 && yd >= 0) || (xd <= 0 && yd <= 0)) {		if ( fabs(yd) >= fabs(xd) ) {	  front = v3 + (xd / yd) * (v2 - v3); /* octant 2 */	  back = v7 + (xd / yd) * (v6 - v7);  /* octant 6 */	  if (V > front && V > back)	    edges[row][col] = 255; /* mark edge */	  else	    edges[row][col] = 0; /* no edge */ }		else if( fabs(yd) <= fabs(xd) ) {	  front = v1 + (yd / xd) * (v2 - v1);  /* octant 1 */	  back = v5 + (yd / xd) * (v6 - v5);   /* octant 5 */	  if (V > front && V > back)	    edges[row][col] = 255; /* mark edge */	  else	    edges[row][col] = 0; /* no edge */ }		else	  printf("Error: nms\n");      }            /* quadrants II and IV */            else if ((xd >= 0 && yd <= 0) || (xd <= 0 && yd >= 0)) {		if ( fabs(yd) >= fabs(xd) ) {	  front = v3 - (xd / yd) * (v4 - v3); /* octant 3 */	  back = v7 - (xd / yd) * (v8 - v7);  /* octant 7 */	  if (V > front && V > back)	    edges[row][col] = 255; /* mark edge */	  else	    edges[row][col] = 0; /* no edge */ }		else if( fabs(yd) <= fabs(xd) ) {	  front = v1 - (yd / xd) * (v8 - v1);  /* octant 8 */	  back = v5 - (yd / xd) * (v4 - v5);   /* octant 4 */	  if (V > front && V > back)	    edges[row][col] = 255; /* mark edge */	  else	    edges[row][col] = 0; /* no edge */ }		else 	  printf("Error: nms\n"); }            else 	printf("Error: non-maximum suppression\n");    }	    /* fill non-valid parts */    for(row = 0; row < vert_size; row++) {  /* fill columns */    edges[row][0] = 0;    edges[row][horiz_size-2] = 0;    edges[row][horiz_size-1] = 0;  }    for(col = 0; col < horiz_size; col++) {  /* fill rows */    edges[0][col] = 0;    edges[vert_size-2][col] = 0;    edges[vert_size-1][col] = 0;  }}/*===========================================================================*//* This function takes as arguments a pointer to a 2-D array of float (the horizontal dimension is defined in the include file gdef.h) anda float value representing the sigma for a gaussian filter. The functionmutates the array by convolving it with the gaussian filter. */void gsmooth(float **img_array, float sigma, int vert_size, int horiz_size){  float gauss[MAX_GAUSS];    /* for filter coefficients */  float worksp[WORKSIZE];    /* for workspace */    register int row, col;  float scale;  int n, gwidth, i, index;  int count, wksprpt, wksplpt, wksptpt, wkspbpt;    /* determine width */  n = (int)(ceil(sqrt(-log(CUTOFF) * 2) * sigma)); /* max n */  gwidth = 2 * n + 1; 	/* gaussian width */    wksplpt= n;	  		/* left end of workspace */  wksprpt= (horiz_size -1) + n;	/* right end of workspace */  wksptpt= n;			/* top of workspace */  wkspbpt= (vert_size - 1) + n;	/* bottom of workspace */  /* calculate scale factor */    scale = 0;  		/* initialize */  for(i = -n; i <= n ; i++) {    scale += exp(-(i*i)/(2 * sigma * sigma));}  /* fill gaussian */    for(i = -n; i <= n; i++) {    index = i + n;    gauss[index] = (exp(-(i*i)/(2 * sigma * sigma)))/scale;}  /* convolve rows */    for(row=0; row < vert_size; row++) {        for(col=0; col < horiz_size; col++)      worksp[(col+n)]=img_array[row][col];        for(count=1; count<=n ; count++) {      worksp[(wksplpt-count)]=worksp[(wksplpt+count)];      worksp[(wksprpt+count)]=worksp[(wksprpt-count)]; }        for(col=0; col < horiz_size; col++)      img_array[row][col] = fconv(gauss,				  &worksp[col],gwidth); }    /* convolve columns */    for(col=0; col < horiz_size ; col++) {        for(row=0; row < vert_size; row++)      worksp[(row+n)]=img_array[row][col];        for(count=1; count<=n; count++) {      worksp[(wksptpt-count)]=worksp[(wksptpt+count)];      worksp[(wkspbpt+count)]=worksp[(wkspbpt-count)]; }    for(row=0; row < vert_size ; row++)      img_array[row][col]= fconv(gauss,					 &worksp[row],gwidth);}}