#include <stdio.h>#include <math.h>#include <string.h>#include "canedge.h"#include "nr.h"#include "dat.c"extern int optind;extern char *optarg;#define OPTIONS "s:h:l:i:"/*===========================================================================*/void print_documentation(s, h, l, i)     float s, h, l;     int i;{  char *message = "\Canedge is a simplified version of an edge detector developed by \\nJohn Canny at MIT's Artificial Intelligence laboratory. The detector\\nfirst smooths an image using a gaussian shaped low pass filter. \\nThen the gradient magnitude of the smoothed image is non-maximum\\nsuppressed in the direction of the gradient.  These maxima are then\\nthresholded at a high threshold, and again at a low threshold (which\\nis some fraction of the high threshold).  Every maxima above the\\nhigh threshold is an edge.  The next step is called edge 'growing'.\\nEvery point above the low threshold (but not above the high threshold)\\nis marked as an edge if it is next to a point above the high threshold.\\nThis process is repeated several times and thus, the edges are\\n'grown'.\n";  printf("\nusage: canedge [-s float] [-h float] [-l float] [-i int] in_im line_im mag_im\n\n");  printf("-s float : smoothing sigma\n");  printf("           Default: %f\n",s);  printf("-h float : high threshold\n");  printf("           Default: %f\n",h);  printf("-l float : low threshold\n");  printf("           Default: %f\n",l);  printf("-i int   : growing iterations\n");  printf("           Default: %d\n",i);    printf("\n%s\n",message);}/*===========================================================================*/void parse_options(argc, argv, s, h, l, i, infile, line_outfile, mag_outfile)     int argc, *i;     char *argv[], *infile, *line_outfile, *mag_outfile;     float *s, *h, *l;{  int errflag = 0;  int fileloc = 0;  char c; while ((c = getopt(argc, argv, OPTIONS)) != EOF)   switch (c) {   case 's':      *s = atof(optarg);      fileloc += 2;      break;    case 'l':      *l = atof(optarg);      fileloc += 2;      break;    case 'i':      *i = atoi(optarg);      fileloc += 2;      break;    case '?':      errflag = 1;      break;    }  if (argc != fileloc + 4) {    print_documentation(*s, *h, *l, *i);    exit(-1); }    fprintf(stderr,"%d s %s \n", fileloc, argv[fileloc+1]);  strcpy(infile, argv[fileloc+1]);  strcpy(line_outfile, argv[fileloc+2]);  strcpy(mag_outfile, argv[fileloc+3]);}/*===========================================================================*//* These functions are part of my first version of a float gaussian   smoother.  The include file sizes.h contains size information. */float convert(i)     char i;{  if (i>=0) return (float)i;  else return((float)(i+256));}float fconv(x,y,n)     register float *x, *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 */scaler(f_array, vert_size, horiz_size)     float **f_array;     int vert_size, 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;}magedges(imgarray, edgearray, vert_size, horiz_size)     float **imgarray;     char **edgearray;     int vert_size, 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. */threshold(mag, edges, hthreshold, thresfactor, passes, low, wk_space,	  vert_size, horiz_size)     float **mag;     char **edges;     float hthreshold;     float thresfactor;     int passes;     char **low;     char **wk_space;     int vert_size, 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);}grow(r_array, chk_array, w_array, vert_size, horiz_size)     char **r_array;     char **chk_array;     char **w_array;     int vert_size, 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(xdir, ydir, percentile, h_val, vert_size, horiz_size)     float **xdir;     float **ydir;     float percentile;     float **h_val;     int vert_size, 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;