db_eigsearch.c

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

  
  if (singleframe) {
    nframe = 1;
    facelist = ivector(1, nframe);
    facelist[1] = frame; }
  else {
    sprintf(filename,"%s",infile);
    if ((fp = fopen(filename, "r")) == NULL)
      myerror("Could not open facelist file.");
    fscanf(fp,"%d",&nframe);
    facelist = ivector(1, nframe);
    for (i=1; i<=nframe; i++)
      fscanf(fp,"%d",facelist+i);
    fclose(fp);
  }



  /* ----  open output data file in outdir ---- */

  sprintf(filename,"%s/db_eigsearch.out",outdir);
  if ((fp = fopen(filename, "w")) == NULL) {
    fprintf(stderr,"ERROR Could not open output file %s \n\n", filename);
    exit(1);
  }
    

  /* ------- DEBUG -------- */
  fprintf(stdout,"%s\n\n",comline);     /* echo command line */
  fprintf(fp,"%s\n\n",comline); 

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

  for (f=1; f <= nframe; f++) {
    
    /* ------- 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;
    
    /* ---- read current frame ----- */
        
    sprintf(filename,"%s/%4d",DBASE_PATH,facelist[f]);
    read_RAW(filename, image, nrow, ncol);
    

    /* ----------------- Process current frame ------------------ */
    
    for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); 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+coll[k]; jj++)
	    for (ii=i-rowl[k]; ii<=i+rowr[k]; ii++) 
	      patch[++c] = image[ii][jj];

	  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);
	  errormap[k][i][j] = fval1;
	  
	  if (fval1<fval) {
	    fval = fval1;
	    rowm[k] = i;
	    colm[k] = j;
	    error[k] = fval;
	  }
	  
	}

      fprintf(fp,"%4d \t %1d \t %3d %3d \t %1.6e \n",
	      facelist[f], k, rowm[k], colm[k], error[k]);
      fflush(fp);
      fprintf(stdout,"%4d \t %1d \t %3d %3d \t %1.6e \n",
	      facelist[f], k, rowm[k], colm[k], error[k]);

      if (do_dmdumps) {
	sprintf(filename,"%s/dmf%dt%d.bf",outdir,facelist[f],k);
	write_BIN(filename, errormap[k], nrow, ncol);
      }

    }
    
    
  } /* end of frame sequence loop */
  
  
  if (writeerrormap>0) {
    write_descriptor("Error", nframe, ncol, nrow, bytes_pixel, comline);
    write_descriptor(outdir, nframe, ncol, nrow, bytes_pixel, comline);
  }
  
  
  /* --- free up allocated matrices ----- */
  
  free_vector(patch, 1, max_N);
  free_imatrix(templatesize, 1, nfeatures, 1, 2);
  free_matrix(image, 1, nrow, 1, ncol);
  free_ivector(facelist, 1, nframe);
  fclose(fp); 
  
  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/N;             /* return MEAN square-error (ie. per pixel) */
}

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

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,std,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, &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;
      val = SQR(C[k]);
      Ec += val;
      Em += val/eigvalues[k];
    }
  }


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

  error = Em + (E - Ec)/eigvalues[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;
}