pms_eigsearch.c

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

      lambda_star[i] /= (N_dim[i] - eig_last);
    }

  }

    
  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 );
     

      /* extract 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];

      /* compute the denominator of umap probability */
      umap_musigma_denom[k] = 
	2.0*PI * sqrt(matrix_determinant(C, 2));

      /* invert the covariance matrix and put it back in same place */
      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];
    }  
    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 );
     
      /* extract 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];
 
      /* compute the denominator of spatial probability */
      spatial_musigma_denom[k] = 
	2.0*PI * sqrt(matrix_determinant(C, 2));

      /* invert the covariance matrix and put it back in same array */
      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];
    }  
    free_matrix(C,1,2,1,2);
    free_matrix(invC,1,2,1,2);
  }


   
  if (do_mixture) {

    int m,n;
    float **vec, **BigS;

    for (k=1; k<=nfeatures; k++) {

      sprintf(filename,"%s/mixture%d_priors.bf",datapath,k);
      vec = read_BIN(filename, &n, &m );
      if (n!=1) {
	fprintf(stderr,"ERROR: mixture%d_priors.bf should be a row vector\n",
		k);
	exit(k);
      }
      mixture[k].m = m; /* number of components */
      mixture[k].P = vector(1, mixture[k].m);
      mixture[k].Denom = vector(1, mixture[k].m);
      for (i=1; i<=mixture[k].m; i++)
	mixture[k].P[i] = vec[1][i];

      sprintf(filename,"%s/mixture%d_means.bf",datapath,k);
      mixture[k].U = read_BIN(filename, &n, &m);
      if (m!=mixture[k].m) {
	fprintf(stderr,
		"ERROR: mixture%d_means.bf should be have %d columns (components)\n",k, mixture[k].m);
	exit(n);
      }
      mixture[k].n = n; /* dimensionality of subspace */

      sprintf(filename,"%s/mixture%d_sigmas.bf",datapath,k);
      BigS = read_BIN(filename, &n, &m);
      if (n!=mixture[k].n) {
	fprintf(stderr,
		"ERROR: mixture%d_sigmas.bf should be have %d rows\n",k, mixture[k].n);
	exit(n);
      }
      if (m!=(mixture[k].n*mixture[k].m)) {
	fprintf(stderr,
		"ERROR: mixture%d_sigmas.bf should be have %d columns\n",k, mixture[k].n*mixture[k].m);
	exit(n);
      }

      for (m=1; m<=mixture[k].m; m++) {
	mixture[k].S[m]    = matrix(1, mixture[k].n, 1, mixture[k].n);
	mixture[k].InvS[m] = matrix(1, mixture[k].n, 1, mixture[k].n);
	n = (m-1)*mixture[k].n;
	for (i=1; i<=mixture[k].n; i++) 
	  for (j=1; j<=mixture[k].n; j++) 
	    mixture[k].S[m][i][j] = BigS[i][n+j];
	/* matrix_inverse(mixture[k].S[m], mixture[k].n, mixture[k].InvS[m]);*/
	matrix_inverse(mixture[k].S[m], eig_last, mixture[k].InvS[m]);
	/* mixture[k].Denom[m] =
	  pow(2.0*PI,mixture[k].n/2.0) * 
	    sqrt(matrix_determinant(mixture[k].S[m], mixture[k].n));*/
	mixture[k].Denom[m] =
	  pow(2.0*PI,eig_last/2.0) * 
	    sqrt(matrix_determinant(mixture[k].S[m], eig_last));
      }
      free_matrix(vec, 1, 1, 1, mixture[k].m);
      free_matrix(BigS, 1, mixture[k].n, 1, mixture[k].n*mixture[k].m);
    }
  }
  
  
  /* ---- 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);
  Warp = matrix(1, 3, 1, 3);



  /* ----  open output data file ---- */

  sprintf(filename,"%s", outfile);
  if ((fp2 = fopen(filename, "w")) == NULL) {
    fprintf(stderr,"ERROR Could not open output file %s \n\n", filename);
    exit(1);
  }
    

  /* ------- 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];
	      
	      if (do_masking) {
		c = 0;
		for (ii=1; ii<=templatesize[k][1]; ii++)
		  for (jj=1; jj<=templatesize[k][2]; jj++)
		    patch[++c] *= template_mask[k][ii][jj];
	      }

	      if (do_correlation) {
		fval1 = 1.0 - norm_corr(patch, N_dim[k], template[k][1],
					template_mean[k], template_sigma[k]);
	      }

	      else if (do_mahalanobis) {
		
		fval1 = mahalanobis(patch, N_dim[k], template[k], 
				    variances[k], lambda_star[k],
				    eig_first, eig_last);
		
		if (do_unitmap) {
		  statistics(patch, N_dim[k], &mean, &sigma);
		  fval1 += umap_musigma_distance(mean, sigma,
						 umap_musigma[k]);
		}
		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]);
		}
	      }
	      
	      else if (do_mixture) {

		fval1 = mixture_probability(patch, N_dim[k], template[k],
					    variances[k], eig_first, eig_last,
					    lambda_star[k], &mixture[k]);
		fval1 = -(fval1);
		if (do_unitmap) {
		  statistics(patch, N_dim[k], &mean, &sigma);
		  fval2 = umap_musigma_distance(mean, sigma,
						umap_musigma[k]);
		  fval2 = exp(-fval2/2.0) / umap_musigma_denom[k];
		  fval1 += -log(fval2+FLT_MIN);
		  
		}
		if (do_spatial) {
		  my_vector[1] = (int) (i/scales[s]+0.5);
		  my_vector[2] = (int) (j/scales[s]+0.5);
		  fval2 =
		    mahalanobis_distance(my_vector, 2,
					 spatial_musigma[k][1],
					 &spatial_musigma[k][1]);
		  fval2 = exp(-fval2/2.0) / spatial_musigma_denom[k];
		  fval1 += -log(fval2+FLT_MIN); 
		}
	      }

	      else  /* just do dffs */ {

		fval1 = dffs(patch,N_dim[k],template[k],eig_first,eig_last);
		
		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);
	  }

	  if (do_minima) {
	    minima(errormap[k], imin[k], imax[k], jmin[k], jmax[k],
		   &myPoints_array[k]);
	    sprintf(filename,"%s/%s_s%d_t%d_minima", outdir, infile, s, k);
	    if ((fp1 = fopen(filename, "w")) == NULL) {
	      fprintf(stderr,"ERROR Could not open minima file %s \n\n", filename);
	      exit(1);
	    }
	    for (i=1; i<=myPoints_array[k].n; i++) 
	      fprintf(fp1, "%3d\t%3d\t%1.6e\n",
		      (int) (myPoints_array[k].x[i]/scales[s]+0.5),
		      (int) (myPoints_array[k].y[i]/scales[s]+0.5),
		      myPoints_array[k].f[i]);
	    fclose(fp1);
	  }

	  
	  /* -------------- DEBUG ------------
	  c = 0;
	  j = colm[s][k];
	  i = rowm[s][k];
	  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];
	  ------------------------------------ */
	  
	} /* 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",