pms_eigsearch.c

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

/*----------------------------------------------------------------------
PROGRAM: pms_eigsearch.c
DATE:    11/4/94
AUTHOR:  Baback Moghaddam, [e-mail]
------------------------------------------------------------------------


  Multiple-Scale Local feature search by eigentemplates with a
  Gaussian-Mixture probabilistic option


  This routine looks for 2 (3) BF files in the datapath (defined below):

  1. features.bf  which is an N-by-2 matrix defining the size (row,col)
     of each of the N templates.

  2. template[n].bf  where n=1...N is the eigenvector ROW matrix where
     the 1st row is the mean feature, and the remaining rows
     are the principal eigenvectors 1:M. 

     NOTE: The eigenvectors are stored in column-order (as in MATLAB)

  3. variances.bf which is an N-by-M matrix, where each row represents
     the rank-ordered eigenvalues associated with that feature


     NOTE: The set of scales over which the search is to be performed
           is provided in an ASCII file of the format:

	   s1
	   s2
	   .
	   .
	   .
	   sN

----------------------------------------------------------------------

   NOTE: In this version if unitmap'ing is selected it will use
         the Mahalanobis mean/stddev in addition to the standard
	 distance computation.

------------------------------------------------------------------------

     Configure this code using the following defines:
     
     DO_VAR:     will set each pixel whose surrond's variance
                 is < VAR_THRESHOLD to VAR_MAXVALUE
	       
---------------------------------------------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "util.h"
#include "io.h"
#include "matrix.h"
#include "affine.h"

struct Points {
  int n;
  int *x;
  int *y;
  float *f;
};

struct Point {
  int x;
  int y;
  float f;
};


#define  PI  3.1415926535

/* ----------- CONFIGURES ---------------- */

#define DO_VAR                  0     /*  variance check     */
#define VAR_THRESHOLD           100   /*  threshold for low-var patches */
#define VAR_MAXVALUE            1e7   /*  replacement in distance map   */

#define MAX_NUM_TEMPLATES       5
#define MAX_NUM_SCALES          32
#define MAX_NUM_MIXTURE_COMPONENTS 32

#define MAX_NUM_EIGENVECTORS    128   /* for static storage in dffs() */
#define MAX_VECTOR_LENGTH       128    /* for static storage in mahalanobis_distance() */
#define MAX_CHARS               512

/* ----------- Command-Line Parsing Stuff ------- */

extern int optind;
extern char *optarg;
char *progname;          /* used to store the name of the program  */
char comline[MAX_CHARS];       /* used to store the entire command line  */

#define OPTIONS "i:l:s:d:o:a:b:k:mvguecpfz"

char *usage = "\t-i indir -l list -s scalefile [-d datadir]\n\
\t\t\t[-o outfile] [-a first_ev] [-b last_ev] [-k feature]\n\
\t\t\t[-g] [-u] [-m] [-p] [-c] [-e] [-f] [-z]\n";

char *help = "\
Multiscale Eigentemplate Search\n\n\
-i indir   \t input directory\n\
-l listfile\t ASCII file listing indir files to process (one per line)\n\
-s scalefile\t ASCII file of scales\n\
-d datadir \t data dir holding eigentemplate BF files (default = ./)\n\
-o outfile \t output file (default = ./ms_eigsearch.out)\n\
-a first_ev\t first eigenvector (default = 1)\n\
-b last_ev \t last eigenvector  (default = 5)\n\
-k feature \t search for the k-th feature only\n\
-g         \t pre-process with graymap\n\
-u         \t pre-process with unitmap\n\
-m         \t use Mahalanobis distance instead of DFFS\n\
-v         \t use normalized correlation (distance = 1 - nc)\n\
-p         \t compute gaussian-mixture inverse log-likelihood\n\
-c         \t use spatial priors\n\
-e         \t output detection maps (in indir/Maps)\n\
-f         \t output ASCII minima files to outdir\n\
-z         \t use template masking (requires files: template(i)_mask.bf)\n";


/* --------- data structures ------------- */

struct Mixture {
  int m;
  int n;
  float *P;
  float **U;
  float **S[MAX_NUM_MIXTURE_COMPONENTS];
  float **InvS[MAX_NUM_MIXTURE_COMPONENTS];
  float *Denom;
};


/* --------- Function Prototypes ---------- */

float dffs(float *patch, int N, float **eigvectors, int first, int last);
float mahalanobis(float *patch, int N, float **eigvectors, float *eigvalues,
		  float lambda_star, int first, int last);
float graymap(float *patch, int N);
float unitmap(float *patch, int N);
void statistics(float *p, int N, float *mean, float *sigma);
float umap_musigma_distance(float mean, float sigma, float **umap_musigma);
float mahalanobis_distance(float *x, int N, float *Mu, float **inv_Sigma);
float mixture_probability(float *patch, int N, float **eigvectors,
			  float *eigvalues, int first, int last,
			  float lambda_star, struct Mixture *mixture_p);
void minima(float **A, int imin, int imax, int jmin, int jmax,
            struct Points *Ps);
float norm_corr(float *patch, int N,
		float *template, float t_mean, float t_sigma);


/* -------- Globals ---------------- */

float **template[MAX_NUM_TEMPLATES];
float **template_mask[MAX_NUM_TEMPLATES];
float **errormap[MAX_NUM_TEMPLATES];
float **variances;
float **umap_musigma[MAX_NUM_TEMPLATES];
float umap_musigma_denom[MAX_NUM_TEMPLATES];
float **spatial_musigma[MAX_NUM_TEMPLATES];
float spatial_musigma_denom[MAX_NUM_TEMPLATES];

/* ------- Command Line Defaults ------------- */

int do_dmdumps = 0; /* default is NOT to do errormap dumps  */
int do_graymap = 0; /* default is NOT to graymap each patch */
int do_unitmap = 0; /* default is NOT to unitmap each patch */
int do_mixture = 0; /* default is NOT to compute probabilities */
int do_mahalanobis = 0; /* default is NOT to do this computation */
int do_spatial = 0; /* default is NOT to add [x,y] Mahalanobis dist */
int do_minima  = 0;  /* default is not to output minima files */
int do_masking = 0;  /* default is not to mask input templates */
int do_correlation = 0; /* default is not to do normalized correlation */
int eig_first = 1;  /* search using first 5 eigenvectors by default */
int eig_last  = 5;




/*----------------------------------------------------------------------*/
/* ---------------------------- MAIN ---------------------------------- */

main(int argc, char *argv[])
{

  register int i,j,k,l,ii,jj;
  int f,s,c,feature=1,nframe,nfeatures,sets, bytes_pixel;
  int nrow, ncol, max_N, xc, yc;
  char command[MAX_CHARS],indir[MAX_CHARS],infile[MAX_CHARS];
  char listfile[MAX_CHARS], scalefile[MAX_CHARS], line[MAX_CHARS];
  char datapath[MAX_CHARS],outdir[MAX_CHARS],filename[MAX_CHARS];
  char outfile[MAX_CHARS];

  float **image, **image_orig;
  float fval1, fval2, fval, mean, sigma;
  float maxerror;
  FILE *fp, *fp1, *fp2;         /* for output values dump */

  int rowl[MAX_NUM_TEMPLATES],rowr[MAX_NUM_TEMPLATES];
  int coll[MAX_NUM_TEMPLATES],colr[MAX_NUM_TEMPLATES];
  int imin[MAX_NUM_TEMPLATES],imax[MAX_NUM_TEMPLATES];
  int jmin[MAX_NUM_TEMPLATES],jmax[MAX_NUM_TEMPLATES];
  int N_dim[MAX_NUM_TEMPLATES]; /* the dimensionality of each template */
  int M_dim[MAX_NUM_TEMPLATES]; /* the # of eigenvectors for each template */

  int numscales;
  float scales[MAX_NUM_SCALES];
  int rowm[MAX_NUM_SCALES][MAX_NUM_TEMPLATES];
  int colm[MAX_NUM_SCALES][MAX_NUM_TEMPLATES];
  float error[MAX_NUM_SCALES][MAX_NUM_TEMPLATES];
  float **Warp;
  float *my_vector;
  int xdim, ydim;
  int   best_scale[MAX_NUM_TEMPLATES];
  float best_error[MAX_NUM_TEMPLATES];
  int   best_rowm[MAX_NUM_TEMPLATES];
  int   best_colm[MAX_NUM_TEMPLATES];

  
  int **templatesize;
  float *patch;

  struct Points myPoints_array[MAX_NUM_TEMPLATES];
  int N_minima;

  struct Mixture mixture[MAX_NUM_MIXTURE_COMPONENTS];
  
  float lambda_star[MAX_NUM_TEMPLATES];  /* optimal variance for DFFS */
  float template_mean[MAX_NUM_TEMPLATES];   /* statistics for norm_corr */
  float template_sigma[MAX_NUM_TEMPLATES];  /* statistics for norm_corr */

  /* required input flags */
  
  int errflag   = 0;
  int inflag    = 0;
  int scaleflag = 0;
  int listflag  = 0;
  

  /* command line defaults */

  int singlefeature = 0;  
  int firstframe = 0;
  int lastframe  = 0;
  strcpy(datapath,".");               /* default datapath directory to cwd   */


  /* setup program name and command line strings */
  
  progname = argv[0];
  for (i=0; i<argc; i++)
    strcat(comline, argv[i]),strcat(comline, " ");


  /* set up default outfile name */

  sprintf(outfile,"%s.out",progname); /* default output file is progname.out */  

  /* --- DEBUG quadratic ---- 
  {
    float *x, *y, **C, val;
    
    x = vector(1, 10);
    y = vector(1, 10);
    C = matrix(1, 10, 1, 10);

    for (i=1; i<=10; i++) {
      x[i] = i;
      y[i] = 11-i;
      for (j=1; j<=10; j++) 
	C[i][j] = 1.0*i/j;
    }
    val = matrix_vector_quadratic(x, C, y, 10);
    printf("quad = %f\n",val);
    exit(1);
  }
  ------------------------------ */

	  
  /* ----------------------  Command Line Parse ------------------------ */
  
  while ((c = getopt(argc, argv, OPTIONS)) != EOF)
    switch (c) {
      
    case 'i':
      strcpy(indir, optarg);
      inflag = 1;
      break;
      
    case 'l':
      strcpy(listfile, optarg);
      listflag = 1;
      break;

    case 'd':
      strcpy(datapath, optarg);
      break;

    case 's':
      strcpy(scalefile, optarg);
      scaleflag = 1;
      break;

    case 'k':
      feature = atoi(optarg);
      singlefeature = 1;
      break;

    case 'o':
      strcpy(outfile, optarg);
      break;

    case 'a':
      eig_first  = atoi(optarg);
      break;

    case 'b':
      eig_last = atoi(optarg);
      break;

    case 'g':
      do_graymap = 1;
      break;

    case 'u':
      do_unitmap = 1;
      break;

    case 'e':
      do_dmdumps = 1;
      break;

    case 'f':
      do_minima = 1;
      break;

    case 'm':
      do_mahalanobis = 1;
      break;

    case 'p':
      do_mixture = 1;
      break;

    case 'v':
      do_correlation = 1;
      break;

    case 'c':
      do_spatial = 1;
      break;

    case 'z':
      do_masking = 1;
      break;

    case '?':
      errflag = 1;
      break;
      
    }
  
  
  /* command line error check */
  
  if (errflag || !inflag || !scaleflag || !listflag) {
    fprintf(stderr,"\nUSAGE: %s %s\n%s\n", progname, usage, help);
    exit(1);
  }

  sprintf(outdir,"%s/Maps",indir);   /* default output dir is indir/Maps    */

  /* override mixture likelihood with Mahalanobis */
  if (do_mixture==1) 
    do_mahalanobis = 0;


  /* -------- Load eigentemplate data --------------- */

  sprintf(filename,"%s/features.bf",datapath);
 
  { 
    float **m;
    max_N = 0;
    m = read_BIN(filename, &nfeatures, &ncol);
    templatesize = imatrix(1, nfeatures, 1, 2);
    
    for (i=1; i<=nfeatures; i++) {
      for (j=1; j<=ncol; j++)
	templatesize[i][j] = (int) m[i][j];
      N_dim[i] = templatesize[i][1]*templatesize[i][2];
      if (N_dim[i]>max_N) 
	max_N = N_dim[i];
      free_matrix(m, 1, nfeatures, 1, ncol);
    }
  }

  
  for (i=1; i<=nfeatures; i++) {
    sprintf(filename,"%s/template%d.bf",datapath,i);
    template[i] = read_BIN(filename, &nrow, &ncol);
    if (do_correlation) {
      statistics(template[i][1], ncol, &template_mean[i], &template_sigma[i]);
      fprintf(stdout,"Template-%d \t mean = %f \t sigma = %f\n\n",
	      i, template_mean[i], template_sigma[i]);
    }
    M_dim[i] = nrow;
    if (ncol!=N_dim[i])
      myerror("Template sizes in BF don't match those is definition file");

    if (do_masking) {
      sprintf(filename,"%s/template%d_mask.bf",datapath,i);
      template_mask[i] = read_BIN(filename, &nrow, &ncol);
      if ((nrow != templatesize[i][1]) || (ncol != templatesize[i][2])) {
	fprintf(stderr,
		"ERROR: template%d_mask.bf dimensions (%d,%d) are invalid\n\n",
		i, nrow, ncol);
	exit(1);
      }
    }
  }


  if (do_mahalanobis || do_mixture) {
    sprintf(filename,"%s/variances.bf",datapath);
    variances = read_BIN(filename, &nrow, &ncol);
    if (nrow!=nfeatures) {
      fprintf(stderr,
	      "ERROR: variances.bf has %d rows! (must have %d)\n\n",
	      nrow, nfeatures);
      exit(1);
    }
    if (eig_last>ncol) {
      fprintf(stderr,
	      "ERROR: variances.bf has %d cols! (must have atleast %d)\n\n",
	      ncol, eig_last);
      exit(1);
    }

    /* now estimate lambda_star */
    for (i=1; i<=nfeatures; i++) {
      lambda_star[i] = 0.0;
      for (j=eig_last+1; j<=ncol-1; j++)  /* sum to next to last in case <=0 */
	lambda_star[i] += variances[i][j];
      lambda_star[i] += ((N_dim[i]-ncol+1) * variances[i][ncol-1]);