/* ------- 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], eig_first, eig_last);	      }	      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);	  }	  	} /* 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);   fclose(fp);  fclose(fp2);     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;}/* ---------------------------------------------------------------------- */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,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)/eigvalues[last+1];  /* 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;}