eigsearch.c

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

    printf("row/col l/r %d: %3d %3d %3d %3d\n", k, rowl[k], rowr[k],
	   coll[k], colr[k]);
  --------------------------------------------------------------------- */


  if (do_minima) {
    int i1 = imin[feature], i2 = imax[feature];
    int j1 = jmin[feature], j2 = jmax[feature];
    int numpixels = (i2-i1)*(j2-j1), n;
    
    for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) {
      n = (imax[k]-imin[k])*(jmax[k]-jmin[k]);
      if (n>numpixels) {
	numpixels = n;
	i1 = imin[k];
	i2 = imax[k];
	j1 = jmin[k];
	j2 = jmax[k];
      }
    }
    N_minima = (i2-i1)*(j2-j1) / 9;
    for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) {
      myPoints_array[k].n = N_minima;
      myPoints_array[k].x = ivector(1, N_minima);
      myPoints_array[k].y = ivector(1, N_minima);
      myPoints_array[k].f =  vector(1, N_minima);
    }
  }
  
  
  /* --- check that specified template sizes and windows match image --- */

  for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) 
    if ((imin[k]-rowl[k]<1) || (imax[k]+rowr[k]>nrow) || \
	(jmin[k]-coll[k]<1) || (jmax[k]+colr[k]>ncol)) {
      fprintf(stderr,
	      "ERROR: image dimensions %d-by-%d don't match templates\n",
	      nrow, ncol);
      exit(3);
    }    
  
  
  /* ----  open output data file in outdir ---- */
  
  sprintf(filename,"%s/%s.out", outdir, progname);
  if ((fp2 = fopen(filename, "w")) == NULL) {
    fprintf(stderr,"ERROR Could not open output file %s \n\n", filename);
    exit(1);
  }
    

  /* ------- 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, nrow, ncol);
      else
	read_RAW(filename, image, nrow, ncol);
      
      
      /* ------- 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;
      
      /* ------- 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+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] = i;
	      my_vector[2] = j;
	      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[k] = i;
	      colm[k] = j;
	      error[k] = fval;
	    }
	    
	  }
	
	fprintf(stdout,"%s\t %1d \t %3d %3d \t %1.6e \n",
		infile, k, rowm[k], colm[k], error[k]);

	fprintf(fp2,"%s \t %1d   %3d %3d   %1.6e \n",
		infile, k, rowm[k], colm[k], error[k]);
	fflush(fp2);
	
	if (do_dmdumps) {
	  sprintf(filename,"%s/%s_t%d", outdir, infile, k);
	  write_RAW_float(filename, errormap[k], nrow, ncol);
	}


	if (do_minima) {
	  minima(errormap[k], imin[k], imax[k], jmin[k], jmax[k],
		 &myPoints_array[k]);
	  sprintf(filename,"%s/%s_t%d_minima", outdir, infile, k);
	  if ((fp1 = fopen(filename, "w")) == NULL) {
	    fprintf(stderr,"ERROR Could not open minima file %s \n\n", filename);
	    exit(1);
	  }
	  for (i=1; i<=myPoints_array[k].n; i++) 
	    fprintf(fp1, "%3d\t%3d\t%1.6e\n",
		    myPoints_array[k].x[i], myPoints_array[k].y[i],
		    myPoints_array[k].f[i]);
	  fclose(fp1);
	}

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

  } /* end of frame sequence loop */
  
  
  
  /* --- write outdir descriptor --- */

  if (do_dmdumps>0) {
    write_descriptor(outdir, nframe, ncol, nrow, 4, comline);
  }
  

  
  /* --- free up allocated matrices ----- */
  
  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);

  if (do_minima) {
    for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) {
      free_ivector(myPoints_array[k].x, 1, N_minima);
      free_ivector(myPoints_array[k].y, 1, N_minima); 
      free_vector(myPoints_array[k].f, 1, N_minima);   
    }
  }

  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;
}