/*----------------------------------------------------------------------PROGRAM: new_ms_eigsearch.cDATE:    2/24/95AUTHOR:  Baback Moghaddam, baback@media.mit.edu------------------------------------------------------------------------  Multiple-Scale Local feature search by eigentemplates  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.	 Also fixed the unitmap musigma cost 	 Includes patch dump additions from rena------------------------------------------------------------------------     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;};/* ----------- 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_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:gumnvcefp"      /* RENA added p option */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] [-n] [-v] [-e] [-c] [-f] [-p]\n";       /* RENA added [-p] */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\-n         \t use optimal weight rho^* (default = lambda_{last_ev+1})\n\-v         \t use normalized correlation (distance = 1 - nc)\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\-p         \t output image patches to outdir\n";          /* RENA added -p help description *//* --------- 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);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 **errormap[MAX_NUM_TEMPLATES];float **variances;float **umap_musigma[MAX_NUM_TEMPLATES];float **spatial_musigma[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_mahalanobis = 0; /* default is NOT to do this computation */int do_optimal_rho = 0; /* default is the old Mahalanobis 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_correlation = 0; /* default is not to do normalized correlation *//* RENA begin */int do_patch_dump = 0; /* default is not to do patch dumps *//* RENA end */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;  /* RENA begin */  unsigned char **image_patch[MAX_NUM_TEMPLATES];  /* RENA end */  struct Points myPoints_array[MAX_NUM_TEMPLATES];  int N_minima;  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 */    /* ----------------------  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 'n':      do_optimal_rho = 1;      break;    case 'v':      do_correlation = 1;      break;     case 'c':      do_spatial = 1;      break;    /* RENA begin */    case 'p':      do_patch_dump = 1;      break;    /* RENA end */    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    */  /* -------- 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);      /* RENA begin */      if (do_patch_dump) {	image_patch[i] = cmatrix(1, templatesize[i][1], 1, templatesize[i][2]);	/* RENA end */      }    }  }    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_mahalanobis) {    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++) {      if (do_optimal_rho>0) {	/* average eigenvalues in orthogonal subspace */	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]);	lambda_star[i] /= (N_dim[i] - eig_last);      }      else {	lambda_star[i] = variances[i][eig_last+1];      }    }  }    if (do_unitmap) {    float **C, **invC;    int m, n;    C = matrix(1,2,1,2);    invC = matrix(1,2,1,2);    for (k=1; k<=nfeatures; k++) {      sprintf(filename,"%s/umap_musigma%d.bf",datapath,k);      umap_musigma[k] = read_BIN(filename, &m, &n );           /* invert the covariance matrix portion of umap_musigma[k] */      for (i=1; i<=2; i++)	for (j=1; j<=2; j++)	  C[i][j] = umap_musigma[k][i+1][j];      matrix_inverse(C, 2, invC);      for (i=1; i<=2; i++)	for (j=1; j<=2; j++)	  umap_musigma[k][i+1][j] = invC[i][j];      /* ----------- DEBUG -----------      for (i=1; i<=3; i++) {	for (j=1; j<=2; j++)	  printf("%f ",umap_musigma[k][i][j]);	printf("\n");      }      printf("\n\n");      C[1][1] = 10;      C[1][2] = 10;      printf("dist = %f \n", 	     mahalanobis_distance(C[1], 2, umap_musigma[k][1], &umap_musigma[k][1]));      printf("dist = %f \n",		    umap_musigma_distance(10, 10, umap_musigma[k]));      exit(1);      --------------- DEBUG ------------- */    }      free_matrix(C,1,2,1,2);    free_matrix(invC,1,2,1,2);  }     if (do_spatial) {    float **C, **invC;    int m, n;    C = matrix(1,2,1,2);    invC = matrix(1,2,1,2);    for (k=1; k<=nfeatures; k++) {      sprintf(filename,"%s/spatial_musigma%d.bf",datapath,k);      spatial_musigma[k] = read_BIN(filename, &m, &n );           /* invert the covariance matrix portion of umap_musigma[k] */      for (i=1; i<=2; i++)	for (j=1; j<=2; j++)	  C[i][j] = spatial_musigma[k][i+1][j];      matrix_inverse(C, 2, invC);      for (i=1; i<=2; i++)	for (j=1; j<=2; j++)	  spatial_musigma[k][i+1][j] = invC[i][j];      /* ----------- DEBUG -----------      for (i=1; i<=3; i++) {	for (j=1; j<=2; j++)	  printf("%f ",spatial_musigma[k][i][j]);	printf("\n");      }      printf("\n\n");      C[1][1] = 10;      C[1][2] = 10;      printf("dist = %f \n", 	     mahalanobis_distance(C[1], 2, spatial_musigma[k][1], &spatial_musigma[k][1]));      printf("dist = %f \n",		    spatial_musigma_distance(10, 10, spatial_musigma[k]));      exit(1);      --------------- DEBUG ------------- */    }      free_matrix(C,1,2,1,2);    free_matrix(invC,1,2,1,2);  }  /* ---- determine left/right boundaries for errormaps ----- */  for (i=1; i<=nfeatures; i++) {    rowl[i] = (int) templatesize[i][1]/2.0;    rowr[i] = templatesize[i][1] - rowl[i] - 1;    coll[i] = (int) templatesize[i][2]/2.0;    colr[i] = templatesize[i][2] - coll[i] - 1;  }  patch = vector(1, max_N);  my_vector = vector(1, MAX_VECTOR_LENGTH);  /* ----  read indir descriptor -------- */    read_descriptor(indir, &nframe, &sets, &bytes_pixel, &ncol, &nrow);  if (sets>1)     myerror("Input files must be single-set DAT files!");  image = matrix(1, nrow, 1, ncol);  image_orig = matrix(1, nrow, 1, ncol);  /* --- allocate errormap tensor and set window boundaries --- */  for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); k++) {    errormap[k] = matrix(1, nrow, 1, ncol);    imin[k] = rowl[k] + 1;    imax[k] = nrow - rowr[k];    jmin[k] = coll[k] + 1;    jmax[k] = ncol - 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);    }  }  /* -------- Load scalefile -------------- */  sprintf(filename,"%s",scalefile);  if ((fp = fopen(filename, "r")) == NULL) {    fprintf(stderr,"ERROR Could not open scale file %s \n\n", filename);    exit(1);  }  i = 1;  while (fgets(line, MAX_CHARS, fp)) {    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {      sscanf(line, "%f", scales+i);      i++;    }  }  numscales = i-1;  fclose(fp);