canny.c

FERET人脸库的处理代码。内函预处理

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 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. */

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) and
a float value representing the sigma for a gaussian filter. The function
mutates 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);}
}