canedge.c

feret人脸图象数据库处理代码

  return(hthreshold);
}

/*===========================================================================*/

/* This function takes as arguments an array to operate on, and three
arrays to fill with results.  The function differentiates the source
array 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 stored
in the third array (character array).
Note: Realize that the last row and column of the x derivative, y
derivative, 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) and
a float value representing the sigma for a gaussian filter. The function
mutates 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);
}