if (singleframe) {    nframe = 1;    facelist = ivector(1, nframe);    facelist[1] = frame; }  else {    sprintf(filename,"%s",infile);    if ((fp = fopen(filename, "r")) == NULL)      myerror("Could not open facelist file.");    fscanf(fp,"%d",&nframe);    facelist = ivector(1, nframe);    for (i=1; i<=nframe; i++)      fscanf(fp,"%d",facelist+i);    fclose(fp);  }  /* ----  open output data file in outdir ---- */  sprintf(filename,"%s/db_eigsearch.out",outdir);  if ((fp = fopen(filename, "w")) == NULL) {    fprintf(stderr,"ERROR Could not open output file %s \n\n", filename);    exit(1);  }      /* ------- DEBUG -------- */  fprintf(stdout,"%s\n\n",comline);     /* echo command line */  fprintf(fp,"%s\n\n",comline);   /* ------------------------------------------------------ */  /* -------------- MAIN Sequence Loop -------------------- */  /* ------------------------------------------------------ */      for (f=1; f <= nframe; f++) {        /* ------- 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;        /* ---- read current frame ----- */            sprintf(filename,"%s/%4d",DBASE_PATH,facelist[f]);    read_RAW(filename, image, nrow, ncol);        /* ----------------- Process current frame ------------------ */        for (k=feature; k<=(singlefeature>0 ? feature:nfeatures); 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+coll[k]; jj++)	    for (ii=i-rowl[k]; ii<=i+rowr[k]; ii++) 	      patch[++c] = image[ii][jj];	  if (do_mahalanobis)	    fval1 = mahalanobis(patch, N_dim[k], template[k], variances[k],				eig_first, eig_last);	  else	    fval1 = dffs(patch,N_dim[k],template[k],eig_first,eig_last);	  errormap[k][i][j] = fval1;	  	  if (fval1<fval) {	    fval = fval1;	    rowm[k] = i;	    colm[k] = j;	    error[k] = fval;	  }	  	}      fprintf(fp,"%4d \t %1d \t %3d %3d \t %1.6e \n",	      facelist[f], k, rowm[k], colm[k], error[k]);      fflush(fp);      fprintf(stdout,"%4d \t %1d \t %3d %3d \t %1.6e \n",	      facelist[f], k, rowm[k], colm[k], error[k]);      if (do_dmdumps) {	sprintf(filename,"%s/dmf%dt%d.bf",outdir,facelist[f],k);	write_BIN(filename, errormap[k], nrow, ncol);      }    }          } /* end of frame sequence loop */      if (writeerrormap>0) {    write_descriptor("Error", nframe, ncol, nrow, bytes_pixel, comline);    write_descriptor(outdir, nframe, ncol, nrow, bytes_pixel, comline);  }      /* --- free up allocated matrices ----- */    free_vector(patch, 1, max_N);  free_imatrix(templatesize, 1, nfeatures, 1, 2);  free_matrix(image, 1, nrow, 1, ncol);  free_ivector(facelist, 1, nframe);  fclose(fp);     return 0;}/*------------------------------------------------------------------- *//* -------------------- end of main() ------------------------------- *//*------------------------------------------------------------------- */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/N;             /* return MEAN square-error (ie. per pixel) */}/* ---------------------------------------------------------------------- */float mahalanobis(float *patch, int N, float **eigvectors, float *eigvalues,		  int first, int last){  register int i,j,k;  float error,sum,mean,var=1.0,std,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, &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;      val = SQR(C[k]);      Ec += val;      Em += val/eigvalues[k];    }  }  if (DO_VAR)    if (var<VAR_THRESHOLD)      error = VAR_MAXVALUE;  error = Em + (E - Ec)/eigvalues[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;}