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. */gradmax(img_array, xdir, ydir, edges, vert_size, horiz_size)     float **img_array;     float **xdir;     float **ydir;     char **edges;     int vert_size, 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. */gsmooth(img_array, sigma, vert_size, horiz_size)     float **img_array;     float sigma;     int vert_size, 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);}}/*===========================================================================*/void main(argc, argv)     int argc;     char *argv[];{  float h_threshold;  int bytes;  int vert_size, horiz_size;  float **picture, **xdirection, **ydirection, **hgrm_working_space;  char *im, **edges, **w1space, **w2space;  char infile[MAX_FILENAME_CHAR];  char line_outfile[MAX_FILENAME_CHAR];  char mag_outfile[MAX_FILENAME_CHAR];  /* This is the default for smoothing (must be > 0) */  float sigma = CANNY_SIGMA_DEF;  /* This is the default for 'growing' (must be >= 0) */  int passes = CANNY_PASSES_DEF;  /* This is the default for thresholding (must be 0 < hp < 100) */  float h_percent = CANNY_H_PERC_DEF;  /* This is the default for thresholding (must be 0 < tf < 1 */  float thr_factor = CANNY_THR_FACTOR;    parse_options(argc, argv, &sigma, &h_percent, &thr_factor, &passes,		infile, line_outfile, mag_outfile);  /* read input file */  read_image(infile, &im, &horiz_size, &vert_size, &bytes);    if (bytes != 1) {    fprintf(stderr,"Sorry, only reads char files.\n");    exit(-1); }  /* allocate memory */  picture = (float **) matrix(0, vert_size-1, 0, horiz_size-1);  xdirection = (float **) matrix(0, vert_size-1, 0, horiz_size-1);  ydirection = (float **) matrix(0, vert_size-1, 0, horiz_size-1);  hgrm_working_space = (float **) matrix(0, vert_size-1, 0, horiz_size-1);  edges = (char **) cmatrix(0, vert_size-1, 0, horiz_size-1);  w1space = (char **) cmatrix(0, vert_size-1, 0, horiz_size-1);  w2space = (char **) cmatrix(0, vert_size-1, 0, horiz_size-1);    convert_1duc_2df(im, picture, vert_size, horiz_size);	  /* gsmooth smooths an array with a gaussian filter */  gsmooth(picture, sigma, vert_size, horiz_size);  /* gradmax takes the smoothed picture and replaces it     with the gradient magnitude (a mutation). It also     returns the x and y derivatives, and it returns the      map of local maximums. */    gradmax(picture, xdirection, ydirection, edges, vert_size, horiz_size);    /* find_thres takes the x derivative and y derivative     arrays, a histogram percentile value, and an array     (as working space) to estimate the high threshold     for thresholding */    h_threshold = find_thres(xdirection, ydirection, h_percent,			   hgrm_working_space,vert_size, horiz_size);  /* threshold takes arguments: gradient magnitude array,     map of local maximums, a value for high threshold,     the number of passes for 'growing', and two arrays     for work space. It returns the final edge map in     the array of local maximums (a mutation). */    threshold(picture, edges, h_threshold, thr_factor, passes,	    w1space, w2space, vert_size, horiz_size);    /* write edge map and magnitude of gradient at edges */    convert_2dc_1dc(edges, im, vert_size, horiz_size);  write_image(line_outfile, im, horiz_size, vert_size, 1);  magedges(picture, edges, vert_size, horiz_size);  convert_2df_1duc(picture, im, vert_size, horiz_size);  write_image(mag_outfile, im, horiz_size, vert_size, 1);    /* free memory */  free_matrix(picture, 0, vert_size-1, 0, horiz_size-1);  free_matrix(xdirection, 0, vert_size-1, 0, horiz_size-1);  free_matrix(ydirection, 0, vert_size-1, 0, horiz_size-1);  free_matrix(hgrm_working_space, 0, vert_size-1, 0, horiz_size-1);  free_cmatrix(edges, 0, vert_size-1, 0, horiz_size-1);  free_cmatrix(w1space, 0, vert_size-1, 0, horiz_size-1);  free_cmatrix(w2space, 0, vert_size-1, 0, horiz_size-1);}