/*----------------------------------------------------------------------PROGRAM: new_ml_detect.cDATE:    2/27/95AUTHOR:  Baback Moghaddam, baback@media.mit.edu------------------------------------------------------------------------  Feature Geometry Validation using Maximum Likelihood Minimization  This routine reads in BF format files: dmf%dt%.bf from the indir  and a feature model matrix (in BF) format as well as a geometry_musigma.bf  file  And dumps the coods of best feature matches in output file	       ---------------------------------------------------------------------- */#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#include <limits.h>#include <float.h>#include "util.h"#include "io.h"struct Points {  int n;  int *x;  int *y;  float *f;};struct Point {  int x;  int y;  float f;};#include "functions.h"/* ----------- CONFIGURES ---------------- */#define MAX_NUM_TEMPLATES         6#define MAX_NUM_MINIMA          100    /* max num of minima for array storage */#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:m:g:l:o:s:t:fn"char *usage = "\t-i indir -l list -m model.bf -g geometry_musigma.bf \n\\t\t\t[-o outfile] [-s mean_scale] [-t top_N_max] [-f] [-n]\n";char *help="\Eigenfeature Detection with Maximum Likelihood Matching\n\n\-i indir    \t input directory where distancemaps (or minima) reside\n\-l list     \t ASCII list of filenames to process (one per line)\n\-m model.bf \t BF-format file containing canonical feature locations\n\-g priors.bf\t geometry_musigma.bf file containing geometry priors\n\-o outfile  \t output file\n\-s mean_scale\t float, expected scale of feature constellation (default = 1.0)\n\-t top_N_max\t integer, # of local minima for each feature (default = 10)\n\-f          \t read minima files instead of distancemaps\n\-n          \t use new Mahalanobis distance for shape (default = old)\n";/*----------------------------------------------------------------------*//* ---------------------------- MAIN ---------------------------------- */main(int argc, char *argv[]){  register int i,j,k,l,ii,jj;  int f, c, nframe, nfeatures, sets, bytes_pixel;  int nrow, ncol, image_nrow, image_ncol, max_N, xc, yc;  char command[MAX_CHARS],indir[MAX_CHARS],listfile[MAX_CHARS],       modelfile[MAX_CHARS],outfile[MAX_CHARS],line[MAX_CHARS],       infile[MAX_CHARS],filename[MAX_CHARS],geometryfile[MAX_CHARS];  float **image;  float fval1, fval2, fval;  float maxerror;  FILE *fp, *fp1, *fp2;         /* for output values dump */  int i_min, i_max, j_min, j_max, N_minima;  struct Points myPoints_array[MAX_NUM_TEMPLATES];  int top_N[MAX_NUM_TEMPLATES];  int **Combination_matrix;  int num_combinations_max;  int num_combinations;  int best_c;  float **Model;  float **Match;  float **Stiffness;  float match_scale, match_energy;  int Geom_Dim;  float *Geom_Mean;  float **Geom_Cov;  float **Geom_ICov;  /* required input flags */    int errflag   = 0;  int inflag    = 0;  int outflag   = 0;  int modelflag = 0;  int listflag  = 0;  int geomflag  = 0;  /* command line defaults */  float mean_scale = 1.0;  int top_N_max = 10;  int do_read_minima = 0;      /* default is to read detection maps */  int new_geometry = 0;        /* default is the old dist/angle method */  /* setup program name and command line strings */    progname = argv[0];  for (i=0; i<argc; i++)    strcat(comline, argv[i]),strcat(comline, " ");  sprintf(outfile,"%s.out", progname);            /* default output file */      /* ----------------------  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 'm':      strcpy(modelfile, optarg);      modelflag = 1;      break;    case 'g':      strcpy(geometryfile, optarg);      geomflag = 1;      break;    case 'o':      strcpy(outfile, optarg);      outflag = 1;      break;    case 's':      mean_scale = atof(optarg);      break;    case 't':      top_N_max = atoi(optarg);      break;    case 'f':      do_read_minima = 1;      break;    case 'n':      new_geometry = 1;      break;    case '?':      errflag = 1;      break;          }      /* command line error check */    if (errflag || !inflag || !modelflag || !listflag || !geomflag) {    fprintf(stderr,"\nUSAGE: %s %s\n%s\n", progname, usage, help);    exit(1);  }    /* ----  read indir descriptor file -------- */    if (do_read_minima == 0) {    read_descriptor(indir, &nframe, &sets, &bytes_pixel, 		    &image_ncol, &image_nrow);    if (sets>1)       myerror("Input files must be single-set DAT files!");    if (bytes_pixel != 4) {      fprintf(stderr,"ERROR: detection maps in %s must be float!\n", indir);      exit(1);    }    image = matrix(1, image_nrow, 1, image_ncol);  }  /* ---- read model file ----- */  Model = read_BIN(modelfile, &nrow, &ncol);  if (nrow<2) {    fprintf(stderr,"ERROR: Model file %s must have 2 rows!\n\n",	    modelfile);    exit(1);  }  nfeatures = ncol;  /* ---- read geometry file ----- */  {     float **M;    M = read_BIN(geometryfile, &nrow, &ncol);    if (new_geometry>0)       Geom_Dim = (nfeatures-1)*2;         /* new stuff has just x,y pairs */    else {      Geom_Dim = 0;      for (i=1; i<= nfeatures; i++)	for (j=i+1; j<=nfeatures; j++)	  Geom_Dim++;      Geom_Dim = 2 * Geom_Dim;             /* inter-dist + inter-angles */    }    if (ncol != Geom_Dim) {      fprintf(stderr,"ERROR: Geometry file %s must have %d cols not %d!\n\n",	      geometryfile, Geom_Dim, ncol);      exit(2);    }    if (nrow != Geom_Dim+1) {     fprintf(stderr,"ERROR: Geometry file %s must have %d rows not %d!\n\n",	      geometryfile, Geom_Dim+1, nrow);      exit(2);    }        Geom_Mean = vector(1, Geom_Dim);    Geom_Cov  = matrix(1, Geom_Dim, 1, Geom_Dim);    Geom_ICov = matrix(1, Geom_Dim, 1, Geom_Dim);        for (i=1; i<=Geom_Dim; i++)      Geom_Mean[i] = M[1][i];    for (i=1; i<=Geom_Dim; i++)      for (j=1; j<=Geom_Dim; j++)	Geom_Cov[i][j] = M[i+1][j];    matrix_inverse(Geom_Cov, Geom_Dim, Geom_ICov);    free_matrix(M, 1, nrow, 1, ncol);  }  /* ---- now determine max # of minima and allocate space ------ */  N_minima = MAX_NUM_MINIMA;  if (do_read_minima == 0) {    i_max = image_nrow;    j_max = image_ncol;    i_min = j_min = 1;    N_minima = (i_max - i_min) * (j_max - j_min) / 9;  }  for (k=1; k<=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);  }  /* ----  set up minimization stuff ---- */  Stiffness = matrix(1, nfeatures, 1, nfeatures);  for (i=1; i<=nfeatures; i++)    for (j=1; j<=nfeatures; j++)      Stiffness[i][j] = 1.0;  num_combinations_max = 1;  for (i=1; i<=nfeatures; i++) {    top_N[i] = top_N_max;    num_combinations_max *= top_N_max;  }  Combination_matrix = imatrix(1, num_combinations_max, 1, nfeatures);  Match = matrix(1, 3, 1, nfeatures);    /* ----  open output data file in outdir ---- */  sprintf(filename,"%s",outfile);  if ((fp = fopen(filename, "w")) == NULL) {    fprintf(stderr,"ERROR Could not open output file %s \n\n", filename);    exit(1);  }      fprintf(stdout,"%s\n\n",comline);     /* echo command line */  fprintf(fp,"# %s\n\n",comline);   /* ------------------------------------------------------ */  /* -------------- MAIN Sequence Loop -------------------- */  /* ------------------------------------------------------ */   /* ---- loop over input list file and warp ------- */  if ((fp2 = fopen(listfile, "r")) == NULL) {    fprintf(stderr,"ERROR: Could not open list file %s \n\n", listfile);    exit(1);  }  nframe = 0;  while (fgets(line, MAX_CHARS, fp2)) {    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {      nframe++;      /* --- read filename  --- */      sscanf(line, "%s", infile);                  /* ---- find/read the local minima in each map -------- */            if (do_read_minima) {	for (k=1; k<=nfeatures; k++) {	  sprintf(filename,"%s/%s_t%d_minima", indir, infile, k);	  if ((fp1 = fopen(filename, "r")) == NULL) {	    fprintf(stderr,		    "ERROR: Could not open minima file %s \n\n", filename);	    exit(1);	  }	  i=1;	  while (fscanf(fp1, "%d %d %f", 			myPoints_array[k].x+i, myPoints_array[k].y+i,			myPoints_array[k].f+i) != EOF  && i<N_minima)	    i++;	  myPoints_array[k].n = i-1;	  fclose(fp1);	}      }      else {	for (k=1; k<=nfeatures; k++) {	  sprintf(filename,"%s/%s_t%d", indir, infile, k);	  read_RAW_float(filename, image, image_nrow, image_ncol);	  minima(image, 1, image_nrow, 1, image_ncol, &myPoints_array[k]);	}      }                        /* ---- minimzation search ------- */                  for (k=1; k<=nfeatures; k++)	if (myPoints_array[k].n<top_N_max)	  top_N[k] = myPoints_array[k].n;      i = 1;      for (k=1; k<=nfeatures; k++)   /* this is superfluous, cuz inside */	i *= top_N[k];               /* enumerate, the same is computed */            num_combinations = enumerate(top_N,				   nfeatures,				   Combination_matrix);      if (new_geometry>0) {		best_c = new_best_ml_match(myPoints_array,				   Combination_matrix, num_combinations, nfeatures,				   Model,				   Match,				   Geom_Mean,				   Geom_ICov,				   Geom_Dim,				   &match_scale,				   &match_energy);       }      else {		best_c = best_ml_match(myPoints_array,			       Combination_matrix, num_combinations, nfeatures,			       Model,			       Match,			       Geom_Mean,			       Geom_ICov,			       Geom_Dim,			       &match_scale,			       &match_energy);      }                  /* ---------- Dump information to stdout --------- */            fprintf(stdout,"FILE : %s\n", infile);      fprintf(stdout,"Minima: ");      for (k=1; k<=nfeatures; k++) 	fprintf(stdout,"%3d ", myPoints_array[k].n);      fprintf(stdout,"\n");      fprintf(stdout,"Best  : ");      for (k=1; k<=nfeatures; k++)	fprintf(stdout,"%3d ", top_N[k]);      fprintf(stdout," -->  %d  combinations\n", num_combinations);      fprintf(stdout,"Rank  : ");      for (k=1; k<=nfeatures; k++)	fprintf(stdout,"%3d ", Combination_matrix[best_c][k]);      fprintf(stdout,"\nScale : %1.6f \nEnergy: %1.6e \n",	      match_scale, match_energy);      for (k=1; k<=nfeatures; k++) {	fprintf(stdout,"%3d %3d\n", 		(int) Match[1][k],  (int) Match[2][k]);      }      fprintf(stdout,"\n");                              /* --------- dump matches to output file ------------- */            fprintf(fp,"%s \t", infile);      for (j=1; j<=2; j++) 	for (k=1; k<=nfeatures; k++) 	  fprintf(fp,"%3d ",  (int) Match[j][k]);      fprintf(fp,"  %1.6f  %1.6e", match_scale, match_energy);      fprintf(fp,"\n");      fflush(fp);    } /* end of valid list filename */  } /* end of frame sequence loop */    fclose(fp);   fclose(fp2);        /* --- free up allocated matrices ----- */  for (k=1; k<=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);  }  free_matrix(Stiffness, 1, nfeatures, 1, nfeatures);   free_imatrix(Combination_matrix, 1, num_combinations_max, 1, nfeatures);  free_matrix(Match, 1, 3, 1, nfeatures);    if (do_read_minima == 0)    free_matrix(image, 1, image_nrow, 1, image_ncol);    free_vector(Geom_Mean, 1, Geom_Dim);  free_matrix(Geom_Cov,  1, Geom_Dim, 1, Geom_Dim);  free_matrix(Geom_ICov, 1, Geom_Dim, 1, Geom_Dim);    return 0;}/*------------------------------------------------------------------- *//* -------------------- end of main() ------------------------------- *//*------------------------------------------------------------------- */