ms_eigsearch.c

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

  /* ------- echo command line  -------- */

  fprintf(stdout,"%s\n\n",comline);
  fprintf(fp2,"# %s\n\n",comline); 


  /* ------------------------------------------------------ */
  /* -------------- MAIN Sequence Loop -------------------- */
  /* ------------------------------------------------------ */
    

  if ((fp = fopen(listfile, "r")) == NULL) {
    fprintf(stderr,"ERROR: Could not open input file %s \n\n", listfile);
    exit(1);
  }

  nframe = 0;
  while (fgets(line, MAX_CHARS, fp)) {

    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {

      nframe++;
      
      /* ---- read current frame ----- */

      sscanf(line, "%s", infile);
      sprintf(filename,"%s/%s", indir, infile);
      if (bytes_pixel == 4)
	read_RAW_float(filename, image_orig, nrow, ncol);
      else
	read_RAW(filename, image_orig, nrow, ncol);
      
      
      for (s=1; s<=numscales; s++) {
	
	/* ------- initialize error map ----- */

	for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++)
	  for (i=1; i<=nrow; i++)
	    for (j=1; j<=ncol; j++)
	      errormap[k][i][j] = 0.0;
	
	/* ------ apply scaling to image_orig ------ */

	for (i=1; i<=3; i++)
	  for (j=1; j<=3; j++)
	    Warp[i][j] = 0.0;
	Warp[1][1] = Warp[2][2] = scales[s];
	Warp[3][3] = 1.0;
	affine_warp(image_orig, image, nrow, ncol, Warp, 0.0);
	xdim = (int) (scales[s]*nrow + 0.5);
	ydim = (int) (scales[s]*ncol + 0.5);
	
	
	/* ----------------- Process current frame/scale --------------- */
	
	for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) {
	  
	  imax[k] = xdim - rowr[k];
	  jmax[k] = ydim - colr[k];
	  
	  fval = FLT_MAX;
	  for (i=imin[k]; i<=imax[k]; i++)
	    for (j=jmin[k]; j<=jmax[k]; j++) {
	      
	      c = 0;
	      for (jj=j-coll[k]; jj<=j+colr[k]; jj++)
		for (ii=i-rowl[k]; ii<=i+rowr[k]; ii++) 
		  patch[++c] = image[ii][jj];
	      
	      fval2 = 0;
	      if (do_unitmap) {
		statistics(patch, N_dim[k], &mean, &sigma);
		fval2 = umap_musigma_distance(mean, sigma, umap_musigma[k]); 
	      }

	      if (do_mahalanobis) {
		fval1 = mahalanobis(patch, N_dim[k], template[k], 
				    variances[k], eig_first, eig_last);
	      }
	      else
		fval1 = dffs(patch,N_dim[k],template[k],eig_first,eig_last);

	      if (do_unitmap) 
		fval1 += fval2;

	      
	      if (do_spatial) {
		my_vector[1] = (int) (i/scales[s]+0.5);
		my_vector[2] = (int) (j/scales[s]+0.5);
		fval1 +=
		  mahalanobis_distance(my_vector, 2,
				       spatial_musigma[k][1],
				       &spatial_musigma[k][1]);
	      }
	      
	      errormap[k][i][j] = fval1;
	      if (fval1<fval) {
		fval = fval1;
		rowm[s][k] = (int) (i/scales[s]+0.5);
		colm[s][k] = (int) (j/scales[s]+0.5);
		error[s][k] = fval;
	      }
	      
	    }
	  

	  { /* ----- best minimum block ------- */
	    float n1,n2,n3,n4,n5,n6,n7,n8,n9;

	    fval = FLT_MAX;
	    for (i=imin[k]+1; i<=imax[k]-1; i++) {
	      for (j=jmin[k]+1; j<=jmax[k]-1; j++) {
		fval1 = errormap[k][i][j];
		if (fval1<fval) {
		  n1 = errormap[k][i-1][j-1];
		  n2 = errormap[k][i-1][j];
		  n3 = errormap[k][i-1][j+1];
		  n4 = errormap[k][i][j-1];
		  n5 = errormap[k][i][j];
		  n6 = errormap[k][i][j+1];
		  n7 = errormap[k][i+1][j-1];
		  n8 = errormap[k][i+1][j];
		  n9 = errormap[k][i+1][j+1];
		  if (n5<n1 & n5<n2 & n5<n3 & n5<n4 & 
		      n5<n6 & n5<n7 & n5<n8 & n5<n9) {
		    rowm[s][k] = (int) (i/scales[s]+0.5);
		    colm[s][k] = (int) (j/scales[s]+0.5);
		    error[s][k] = fval1;
		    fval = fval1;
		  }
		}
	      }
	    }
	  } /* ---- end of best minimum block ---- */


	  fprintf(stdout,"%s\t %1.3f \t %1d \t %3d %3d \t %1.6e \n",
		  infile, scales[s], k, 
		  rowm[s][k], colm[s][k], error[s][k]);
	  
	  if (do_dmdumps) {
	    sprintf(filename,"%s/%s_s%dt%d.bf", outdir, infile, s, k);
	    write_BIN(filename, errormap[k], xdim, ydim);
	  }
	  
	} /* end of feature sequence loop */
      
	
      }  /* end of scale sequence loop */

      
      for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) { 
	best_error[k] = FLT_MAX;
	for (s=1; s<=numscales; s++) 
	  if (error[s][k]<best_error[k]) {
	    best_error[k] = error[s][k];
	    best_scale[k] = s;
	    best_rowm[k] = rowm[s][k];
	    best_colm[k] = colm[s][k];
	  }
	fprintf(fp2,"%s \t %1.3f  %1d   %3d %3d   %1.6e \n",
		infile, scales[best_scale[k]], k, best_rowm[k], 
		best_colm[k], best_error[k]);
	fflush(fp2);
      }    

    } /* end of valid listfile line if clause */
  } /* end of frame sequence loop */
  

  
  
  /* --- free up allocated matrices ----- */

  if (do_unitmap) {
    for (k=1; k<=nfeatures; k++)
      free_matrix(umap_musigma[k], 1, 3, 1, 2);
  }

  free_vector(my_vector, 1, MAX_VECTOR_LENGTH);
  free_vector(patch, 1, max_N);
  free_imatrix(templatesize, 1, nfeatures, 1, 2);
  free_matrix(image, 1, nrow, 1, ncol);
  free_matrix(image_orig, 1, nrow, 1, ncol);
  free_matrix(Warp, 1, 3, 1, 3);
 
  fclose(fp);
  fclose(fp2); 
  
  return 0;

}



/*------------------------------------------------------------------- */
/* -------------------- end of main() ------------------------------- */
/*------------------------------------------------------------------- */


float dffs(float *patch, int N, float **eigvectors, int first, int last)
{
  register int i,j,k;
  float error,sum,mean,var=1.0,std,E,Ec,Em;
  int M = 1;
  
  static float C[MAX_NUM_EIGENVECTORS];


  if (first>0)
    M = last - first + 1;
  
  /* compute the first 2 moments of image patch */
 
  if (DO_VAR) {
    statistics(patch, N, &mean, &var);
    var = SQR(var);
  }


  /* graymap/unitmap the patch */

  if (do_graymap)  
    graymap(patch, N);
  else if (do_unitmap)
    unitmap(patch, N);

  /* subtract the mean from the patch (note mean is first row)
     and compute the Energy of the patch                       */

  E = 0.0;
  for (i=1; i<=N; i++) {
    patch[i] -= eigvectors[1][i];
    E += SQR(patch[i]);
  }


  /* project onto the appropriate eigenvectors */
  
  Ec = Em = 0.0;   /* if first=0 then this is just mean-template SSD */
  if (first>0) {   
    for (k = first+1; k<=last+1; k++) {
      sum = 0.0;
      for (i=1; i<=N; i++)
	sum += patch[i] * eigvectors[k][i];
      C[k] = sum;
      Ec += SQR(C[k]);
    }
  }

  error = (E - Ec);
  if (DO_VAR)
    if (var<VAR_THRESHOLD)
      error = VAR_MAXVALUE;


  return error;
}

/* ---------------------------------------------------------------------- */

float mahalanobis(float *patch, int N, float **eigvectors, float *eigvalues,
		  int first, int last)
{
  register int i,j,k;
  float error,sum,mean,var=1.0,sigma,E,Ec,Em,val;
  int M = 1;
  
  static float C[MAX_NUM_EIGENVECTORS];


  if (first>0)
    M = last - first + 1;
  
  /* compute the first 2 moments of image patch */
  
  if (DO_VAR) {
    statistics(patch, N, &mean, &sigma);
    var = SQR(sigma);
  }


  /* graymap/unitmap the patch */

  if (do_graymap)  
    graymap(patch, N);
  else if (do_unitmap)
    unitmap(patch, N);

  /* subtract the mean from the patch (note mean is first row)
     and compute the Energy of the patch                       */

  E = 0.0;
  for (i=1; i<=N; i++) {
    patch[i] -= eigvectors[1][i];
    E += SQR(patch[i]);
  }


  /* project onto the appropriate eigenvectors */
  
  Ec = Em = 0.0;   /* if first=0 then this is just mean-template SSD */
  if (first>0) {   
    for (k = first+1; k<=last+1; k++) {
      sum = 0.0;
      for (i=1; i<=N; i++)
	sum += patch[i] * eigvectors[k][i];
      C[k] = sum;
      val = SQR(C[k]);
      Ec += val;
      Em += val/eigvalues[k-1];   /* FIX: used to be [k] */
    }
  }


  if (DO_VAR)
    if (var<VAR_THRESHOLD)
      error = VAR_MAXVALUE;

  error = Em + (E - Ec)/eigvalues[last+1];  /* FIX: used to be [last+2] */
  return error;
}


/*---------------------------------------------------------------------- */


float graymap(float *p, int N)

/* Maps the N-dimensional float vector p from its min/max to 0/255 */

{

  register int i;
  float vmin,vmax,range,scale;

  vmin = FLT_MAX;
  vmax = -vmin;

  for (i=1; i<=N; i++) {
    if (p[i]<vmin)
      vmin = p[i];
    if (p[i]>vmax)
      vmax = p[i];
  }

  range = vmax-vmin;
  if (range==0.0)         /* if no variation don't modify array */
    return 0;

  scale = 255.0/range;

  for (i=1; i<=N; i++)
    p[i] = (int) (0.5 + scale*(p[i]-vmin));

  return range;
}



/*---------------------------------------------------------------------- */


float unitmap(float *p, int N)

/* Maps the N-dimensional float vector p to zero-mean/unit-std dev. */
/* It returns the standard deviation (or zero if sigma<FLT_MIN).    */

{

  register int i;
  float mean = 0, sigma = 0;

  statistics(p, N, &mean, &sigma);

  if (sigma>FLT_MIN) {
    for (i=1; i<=N; i++)
      p[i] = (p[i]-mean)/sigma;
    return sigma;
  }
  else {
    for (i=1; i<=N; i++)
      p[i] -= mean;
    return 0;
  }
}


/* ---------------------------------------------------------------------- */

void statistics(float *p, int N, float *mean, float *sigma)

/* computes the mean and standard deviation of the vector p */

{
  register int i;
  float val, mymean = 0, mysigma = 0; 

  for (i=1; i<=N; i++)
    mymean += p[i];
  mymean /= N;

  for (i=1; i<=N; i++) {
    val = p[i] - mymean;
    mysigma += val*val;
  }
  mysigma = sqrt(mysigma/(N-1));

  *mean = mymean;
  *sigma = mysigma;
}


/* ---------------------------------------------------------------------- */

float umap_musigma_distance(float mean, float sigma, float **umap_musigma)

/* This function takes the mean and sigma of the patch and returns
   the Mahalanobis distance according to the 2nd-order statistics
   in the matrix umap_musigma (whose first row is the mean, and 2nd-3rd
   the INVERSE covariance
*/

{
  register int i,j;
  float sum=0, distance=0;
  float x[3],y[3];

  x[1] = mean  - umap_musigma[1][1];
  x[2] = sigma - umap_musigma[1][2];

  for (i=1; i<=2; i++) {
    sum = 0;
    for (j=1; j<=2; j++)
      sum += umap_musigma[i+1][j]*x[j];
    y[i] = sum;
  }

  distance = 0;
  for (i=1; i<=2; i++)
    distance += x[i]*y[i];

  return distance;
}

/* ---------------------------------------------------------------------- */

float mahalanobis_distance(float *x, int N, float *Mu, float **inv_Sigma)

/*
  This function returns the Mahalanobis distance for a feature vector
*/

{
  register int i,j;
  float sum=0, distance=0;
  float x_n[MAX_VECTOR_LENGTH];

  for (i=1; i<=N; i++)
    x_n[i] = x[i] - Mu[i];

  for (i=1; i<=N; i++) {
    sum = 0;
    for (j=1; j<=N; j++)
      sum += inv_Sigma[i][j]*x_n[j];
    distance += x_n[i]*sum;
  }

  return distance;
}