pms_eigsearch.c

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

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

  if (do_mixture) {
    for (k=1; k<=nfeatures; k++) {
      for (i=1; i<=mixture[k].m; i++) {
	free_matrix(mixture[k].S[i], 1, mixture[k].n, 1, mixture[k].n);
	free_matrix(mixture[k].InvS[i], 1, mixture[k].n, 1, mixture[k].n);
	free_matrix(mixture[k].U, 1, mixture[k].n, 1, mixture[k].m);
	free_vector(mixture[k].P, 1, mixture[k].m);
	free_vector(mixture[k].Denom, 1, mixture[k].m);
      }
    }
  }

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


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

  if (do_masking) {
    for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) 
      free_matrix(template_mask[k], 
		  1, templatesize[k][1],
		  1, templatesize[k][2]);
  }
  
  fclose(fp);
  fclose(fp2); 
  
  return 0;

}

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

float norm_corr(float *patch, int N,
		float *template, float t_mean, float t_sigma)
{
  register int i;
  float sum,mean,sigma;

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

  sum = 0;
  for (i=1; i<=N; i++)
    sum += patch[i]*template[i];

  return (sum/N - t_mean*mean)/(t_sigma*sigma);
}

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

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,
		  float lambda_star, 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)/lambda_star; /*FIX: was variances[eig_last+1];*/
  return error;
}

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

float mixture_probability(float *patch, int N, float **eigvectors, 
			  float *eigvalues, int first, int last,
			  float lambda_star, struct Mixture *mixp)

     /* returns the log-likelihood (mixture probability) of the patch 
	for the given Gaussian-Mixture density (it also incorporates
	the DFFS probability)


	NOTE: this routine normalizes the KL-projection coefficients
	      by dividing by sqrt(lambda_1)
	*/
{
  register int i,j,k;
  float sum,val,E, Ec, dffs;
  int M = 1;
  float Pm,Ps,logPs,*X,*Y;

  X = vector(1, mixp->n);
  Y = vector(1, mixp->n);


  if (first>0)
    M = last - first + 1;
  
  if (first != 1) {
    fprintf(stderr,
	    "mixture_probability(): first = %d, should be 1\n", first);
    exit(first);
  }

  if (M > mixp->n) {
    fprintf(stderr,
	    "mixture_probability(): projection dim (%d) greater than mixture dimension (%d)\n", M, mixp->n);
    exit(M);
  }
    

  /* graymap/unitmap the patch */

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



  /* subtract the mean vector */
  E = 0.0;
  for (i=1; i<=N; i++) {
    patch[i] -= eigvectors[1][i];
    E += SQR(patch[i]);
  }

  /* project onto the appropriate eigenvectors
     and normalize by sqrt(lambda_1) for the mixture    */
  val = sqrt(eigvalues[1]);
  Ec = 0.0;
  for (k = first+1; k<=last+1; k++) {
    sum = 0.0;
    for (i=1; i<=N; i++)
      sum += patch[i] * eigvectors[k][i];
    X[k-1] = sum;
    Ec += SQR(sum);
    X[k-1] = sum/val;
  }
  
  /* compute the in-space probability */
  dffs = (E-Ec);
  logPs = -0.5 * (dffs/lambda_star + (N-M)*(log(2.0*PI) + log(lambda_star)));

  /* compute the mixture probability */
  Pm = 0;
  for (i=1; i<=mixp->m; i++) {

    /* subtract the cluster mean */
    for (j=1; j<=mixp->n; j++)
      Y[j] = X[j] - mixp->U[j][i];

    /* compute the quadratic in the exponential */
    /* val = matrix_vector_quadratic(Y, mixp->InvS[i], Y, mixp->n); */
       val = matrix_vector_quadratic(Y, mixp->InvS[i], Y, M);
    /* printf("d%d = %f\n", i, val); */
    val = mixp->P[i] * exp(-val/2.0) / mixp->Denom[i];
    /* printf("P%d = %f\n", i, val); */
    
    Pm += val;
  }

  free_vector(X, 1, mixp->n);
  free_vector(Y, 1, mixp->n);

  return (log(Pm+FLT_MIN) + logPs);
}


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


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