/*----------------------------------------------------------------------PROGRAM: klt.cDATE:    10/24/93AUTHOR:  Baback Moghaddam         baback@media.mit.edu------------------------------------------------------------------------Performs the Discrete Karhunen-Loeve Transform to compute the set ofeigenimages for a set of images.The mean float image is written to the output directory as data0and 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;}