new_ms_eigsearch.c

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

  Warp = matrix(1, 3, 1, 3);



  /* ----  open output data file ---- */

  sprintf(filename,"%s", outfile);
  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_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], lambda_star[k],
				    eig_first, eig_last);
	      }
	      else if (do_correlation) {

		fval1 = 1.0 - norm_corr(patch, N_dim[k], template[k][1],
					template_mean[k], template_sigma[k]);
	      }
	      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);
	  }
	  

	  if (do_minima) {
	    minima(errormap[k], imin[k], imax[k], jmin[k], jmax[k],
		   &myPoints_array[k]);
	    sprintf(filename,"%s/%s_s%d_t%d_minima", outdir, infile, s, 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",
		      (int) (myPoints_array[k].x[i]/scales[s]+0.5),
		      (int) (myPoints_array[k].y[i]/scales[s]+0.5),
		      myPoints_array[k].f[i]);

	    { /* RENA begin --- patch dump block --- */
	      register int iii, iiii, jjjj;
	      int image_i, image_j;
	      int i_orig, j_orig;
	      
	      if (do_patch_dump) {
		
		/* for each minima found, output associated patch to patch_dump file */
		for (iii=1; iii<=myPoints_array[k].n; iii++) {
		  
		  i_orig = myPoints_array[k].x[iii];
		  j_orig = myPoints_array[k].y[iii];
		  
		  for (image_i=i_orig-rowl[k], iiii=1; image_i<=i_orig+rowr[k], iiii<=templatesize[k][1]; image_i++, iiii++) 
		    for (image_j=j_orig-coll[k], jjjj=1; image_j<=j_orig+colr[k], jjjj<=templatesize[k][2]; image_j++, jjjj++)
		      image_patch[k][iiii][jjjj] = (unsigned char) image[image_i][image_j];
		  
		  sprintf(filename,"%s/%s_s%d_t%d_patch_%03d", outdir, infile, s, k, iii);
		  write_RAW(filename, image_patch[k], templatesize[k][1], templatesize[k][2]);
		}
	      }
	    } /* RENA end --- patch dump block --- */
	    fclose(fp1);
	  }
	  
	} /* 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);

  /* RENA begin */
  if (do_patch_dump) {
    for (i=1; i<=nfeatures; i++) {
      free_cmatrix(image_patch[i], 1, templatesize[i][1], 1, templatesize[i][2]);
    }
  }
  /* RENA end */

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