klt.c

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

/*----------------------------------------------------------------------
PROGRAM: klt.c
DATE:    10/24/93
AUTHOR:  Baback Moghaddam
         [e-mail]
------------------------------------------------------------------------

Performs the Discrete Karhunen-Loeve Transform to compute the set of
eigenimages for a set of images.

The mean float image is written to the output directory as data0
and the remaining eigenimages are named data1, data2 .... dataN.


------------------------------------------------------------------------ */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "util.h"
#include "io.h"


/* ----------- CONFIGURES ---------------- */



/* ----------- Command-Line Parsing Stuff ------- */

extern int optind;
extern char *optarg;
char *progname;          /* used to store the name of the program  */
char comline[170];       /* used to store the entire command line  */

#define OPTIONS "i:l:o:n:t:"

char *usage = "-i indir -l list -o outdir [-n #_of_eigvecs] [-t descriptor_dir]\n";

char *help="\
Computes the Discrete KLT for a set of images\n\n\
-i indir  \t input directory\n\
-l list   \t ASCII file listing indir files to be used (one per line)\n\
-o outdir \t output directory\n\
-n integer\t computes (and writes) the first n eigenvectors only\n\
-t desc_dir\t directory containing descriptor file (default = indir)\n";

/* --------- Function Prototypes ---------- */

int compute_covariance(float ***data, int N, int nrow, int ncol, 
		       float **mean, float **Cov);

int compute_eigenvectors(float **Cov, int N, float *D, float **V); 

int compute_eigenimages(float ***data, int N, int nframe,
			int nrow, int ncol, float **V);

/* -------- Globals ---------------- */

#define MAX_CHARS 256

/* ------- Command Line Defaults ------------- */



/*----------------------------------------------------------------------*/
/* ---------------------------- MAIN ---------------------------------- */

main(int argc, char *argv[])
{

  int i, j, f, c, nframe, sets, channels, bytes_pixel;
  int nrow, ncol, nrotations;
  char command[MAX_CHARS],indir[MAX_CHARS],infile[MAX_CHARS],
       listfile[MAX_CHARS],outdir[MAX_CHARS],descdir[MAX_CHARS],
       line[MAX_CHARS],filename[MAX_CHARS];

  float ***images, **meanimage;
  float *D, **V, **Cov;

  FILE *fp;

  /* required input flags */
  
  int errflag   = 0;
  int inflag    = 0;
  int outflag   = 0;
  int listflag  = 0;

  /* command line defaults */
  
  int nframe_out_select = 0;
  int nframe_out;
  int desc_flag = 0;

  /* setup program name and command line strings */
  
  progname = argv[0];
  for (i=0; i<argc; i++)
    strcat(comline, argv[i]),strcat(comline, " ");
  
  
  /* ----------------------  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 'o':
      strcpy(outdir, optarg);
      outflag = 1;
      break;

    case 'n':
      nframe_out = atoi(optarg) + 1;   /* this way optarg is # of ev's */
      nframe_out_select = 1;
      break;

    case 't':
     strcpy(descdir, optarg);
      desc_flag = 1;
      break;

    case '?':
      errflag = 1;
      break;
      
    }
  
  
  /* command line error check */
  
  if (errflag || !inflag || !outflag || !listflag) {
    fprintf(stderr,"\nUSAGE: %s %s\n%s\n", progname, usage, help);
    exit(1);
  }
  


  /* ----  read DAT-file parameters and Images -------- */
  
  /* read_descriptor(indir, &nframe, &sets, &bytes_pixel, &ncol, &nrow); */
  if (desc_flag == 0)
    read_descriptor2(indir, &nframe, &sets, &channels, 
		     &bytes_pixel, &ncol, &nrow);
  else
    read_descriptor2(descdir, &nframe, &sets, &channels, 
		     &bytes_pixel, &ncol, &nrow);
  
  if (sets>1) 
    myerror("Input files must be single-set DAT files!");
  

  /* ---  read listfile ------- */
  
  if ((fp = 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, fp))
    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) 
      nframe++;
  rewind(fp);
  if (nframe_out_select==0)
    nframe_out = nframe;
  if (nframe_out>nframe)
    myerror("Can't print more Eigenimages than used in Computation!");


  /* ---- read in images ------- */

  images = f3tensor(0, nframe-1, 1, nrow*channels, 1, ncol);
  f = 0;
  while (fgets(line, MAX_CHARS, fp)) {
    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {
      sscanf(line, "%s", infile);
      sprintf(filename, "%s/%s", indir, infile);
      if (bytes_pixel==1)
	read_RAW(filename, images[f], nrow*channels, ncol);
      if (bytes_pixel==4) 
	read_RAW_float(filename, images[f], nrow*channels, ncol);
      f++;
    }
  }
  fprintf(stdout,"Read %d  %d-by-%d  %d-channel %s files from %s/\n",
	  f, nrow, ncol, channels, (bytes_pixel==4 ? "float":"uchar"), indir);

  /* ----  write output descriptor file in outdir ---- */

  /* write_descriptor(outdir, nframe_out, ncol, nrow, 4, comline); */
  write_descriptor2(outdir, nframe_out, 1, channels, ncol, nrow, 4, comline); 
  


  /* ---- compute covariance matrix and mean images ---- */

  meanimage = matrix(1, nrow*channels, 1, ncol);
  Cov = matrix(1, nframe, 1, nframe);
  compute_covariance(images, nframe, nrow*channels, ncol, meanimage, Cov);
  fprintf(stdout,"Computed %d-by-%d covariance matrix\n",
	  nframe, nframe);
  


  /* ---- compute sorted eigenvectors of Cov ----- */

  D = vector(1, nframe);
  V = matrix(1, nframe, 1, nframe);
  nrotations = compute_eigenvectors(Cov, nframe, D, V); 
  fprintf(stdout,"Computed eigenvectors:  %d jacobi() rotations\n",
	  nrotations);



  /* ---- compute eigenimages from the eigenvectors ---- */
  
  compute_eigenimages(images, nframe_out, nframe, nrow*channels, ncol, V);
  fprintf(stdout,"Computed %d eigenimages\n", nframe_out-1);


  /* ---- write results to output directory -------- */

  sprintf(filename,"%s/data0",outdir);
  write_RAW_float(filename, meanimage, nrow*channels, ncol);
  for (f=0; f<nframe_out-1; f++) {
    sprintf(filename,"%s/data%d",outdir,f+1);
    write_RAW_float(filename, images[f], nrow*channels, ncol);
  }
  fprintf(stdout,"Wrote %d  %d-by-%d  %d-channel float files to %s\n", 
	  nframe_out, nrow, ncol, channels, outdir);


  /* ----- write eigenvalues to eigenvalues ------- */

  sprintf(filename,"%s/eigenvalues",outdir);
  if ((fp = fopen(filename,"w")) == NULL) {
    fprintf(stderr,"ERROR: Couldn't open eigenvalue file: %s\n", filename);
    exit(1);
  }

  /* ---- FIX: set last eigenvalue to zero (cuz it's junk!) --- */
  D[nframe] = 0.0;

  for (i=1; i<=nframe; i++) 
    fprintf(fp,"%e\n", D[i]);
  fprintf(stdout,"Wrote eigenvalues as %d-by-1 ASCII file to %s\n",
	  nframe, filename);
  fclose(fp);


  /* --- free up allocated matrices ----- */

  free_vector(D, 1, nframe);
  free_matrix(V, 1, nframe, 1, nframe);
  free_matrix(Cov, 1, nframe, 1, nframe);
  free_f3tensor(images, 0, nframe-1, 1, nrow*channels, 1, ncol);
  free_matrix(meanimage, 1, nrow*channels, 1, ncol);

  return 0;

}



/*------------------------------------------------------------------- */
/* -------------------- end of main() ------------------------------- */
/*------------------------------------------------------------------- */

int compute_covariance(float ***data, int N, int nrow, int ncol, 
		       float **mean, float **Cov)

/* NOTE: This routine will subtract the mean from the orginal sequence */

{
  register int i,j,k,l;
  float sum = 0.0;

  for (i=1; i<=nrow; i++)
    for (j=1; j<=ncol; j++) {
      for (k=0, sum=0.0; k<N; k++)
	sum += data[k][i][j];
      mean[i][j] = sum/N;
      for (k=0; k<N; k++)
	data[k][i][j] -= mean[i][j];
    }

  for (i=1; i<=N; i++)
    for (j=i; j<=N; j++) {

      sum = 0.0;
      for (k=1; k<=nrow; k++)
	for (l=1; l<=ncol; l++) 
	  sum += (data[i-1][k][l] * data[j-1][k][l]);

      Cov[i][j] = Cov[j][i] = sum/N;
    }

  return 0;
}


/*----------------------------------------------------------------------*/


int compute_eigenvectors(float **Cov, int N, float *D, float **V)

/* computes sorted eigenvectors for the covariance matrix Cov using
   the Numerical Recipes routines jacobi() and eigsrt()               */

{
  int nrotations;
  
  jacobi(Cov, N, D, V, &nrotations);
  eigsrt(D, V, N);

  return nrotations;
}

/*----------------------------------------------------------------------*/

int compute_eigenimages(float ***data, int N, int nframe, 
			int nrow, int ncol, float **V)

/* computes a set of N eigenimages using the original data (mean-subtracted)
   and the eigenvectors V of the covariance matrix                           

   NOTE: This routine overwrites the data tensor with the eigenimages  */

{
  register int i,j,k,l;
  float sum, *v;
  
  v = vector(0, N-1);

  for (i=1; i<=nrow; i++)
    for (j=1; j<=ncol; j++) { 
      for (k=0; k<N; k++) {
	for (l=0, sum=0.0; l<nframe; l++)
	  sum += (V[l+1][k+1] * data[l][i][j]);
	v[k] = sum;
      }
      for (k=0; k<N; k++)
	data[k][i][j] = v[k];
    }


  /* now normalize them */
  
  for (k=0; k<N; k++) {
    for (i=1, sum=0.0; i<=nrow; i++)
      for (j=1; j<=ncol; j++) 
	sum += SQR(data[k][i][j]);
    sum = sqrt(sum);
    if (sum>(10*FLT_MIN)) {
      for (i=1; i<=nrow; i++)
	for (j=1; j<=ncol; j++) 
	  data[k][i][j] /= sum;
    }
  }
 

  free_vector(v, 0, N-1);
  return 0;
}