/*----------------------------------------------------------------------PROGRAM: eigcode.cDATE:    01/15/94AUTHOR:  Baback Moghaddam         baback@media.mit.edu------------------------------------------------------------------------Reads input list of files & eigenvector DAT files and compues the projectioncoefficients. Outputs ASCII vectors in output directoryThe mean float image is data0 in the eigenvector direcotry and the remaining eigenimages should be named data1, data2 .... dataN-1.------------------------------------------------------------------------ */#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#include <float.h>#include "util.h"#include "io.h"/* ----------- 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:e:o:n:dtrs:m:a:"char *usage = "\t -i indir -l list -e eigdir -o outdir\n\\t\t [-n int] [-d] [-t descriptor_dir]\n\\t\t [-r] [-s sigma] [-m N_its] [-a alpha]\n";char *help = "\Computes the Eigenvector Projection Coefficients for a set of images\n\n\-i indir   \t input directory\n\-l list    \t ASCII file listing indir files to be processed (one per line)\n\-e eigdir  \t DAT-file directory of eigenvectors (same size as input files)\n\-o outdir  \t output directory for coefficients (same filenames as input)\n\-n integer \t computes the first n eigenvector coefficients only\n\-d         \t computes the residual error for each image\n\-t desc_dir\t directory containing input descriptor (default = indir)\n\-r         \t select robust projection method\n\-s sigma   \t sigma factor for robust norm (default = 50)\n\-m N_its   \t number of iterations for robust method (default = 50)\n\-a alpha   \t iteration update coefficient (default = 1.0)\n";#define MAX_CHARS 256/* ---------- Globals for Robust Method ----------- */float ** robust_Phi;float * robust_X;float * robust_E;float * robust_W;int     robust_Ndim;float robust_project(float *X, int N, float **Phi, int M,		     float *C, float sigma, int N_its, float alpha);/*----------------------------------------------------------------------*//* ---------------------------- MAIN ---------------------------------- */main(int argc, char *argv[]){  register int i, j, k, l;  int f, c, nframe, neigframe, sets, channels, bytes_pixel, nrow, ncol,       in_nrow, in_ncol, in_sets, in_channels;  char command[MAX_CHARS],infile[MAX_CHARS],eigdir[MAX_CHARS],       outdir[MAX_CHARS],filename[MAX_CHARS],listfile[MAX_CHARS],       line[MAX_CHARS],indir[MAX_CHARS],descdir[MAX_CHARS];  float  fval, **image, **test, error, ***eigimages, *coeffs;  FILE *fp1, *fp2;  /* required input flags */    int errflag   = 0;  int inflag    = 0;  int outflag   = 0;  int eigflag   = 0;  int listflag  = 0;  /* command line defaults */    int neigframe_select = 0;  int neigframe_in;  int residual_select = 0;  int desc_flag = 0;  int robust_select  =  0;  float robust_sigma = 50;  float robust_alpha = 1.0;  int   robust_N_its = 50;  /* 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 'e':      strcpy(eigdir, optarg);      eigflag = 1;      break;    case 'o':      strcpy(outdir, optarg);      outflag = 1;      break;    case 'n':      neigframe_in = atoi(optarg) + 1;      neigframe_select = 1;      break;    case 'd':      residual_select = 1;      break;    case 't':      strcpy(descdir, optarg);      desc_flag = 1;      break;    case 'r':      robust_select = 1;      break;    case 's':      robust_sigma = atof(optarg);      break;    case 'm':      robust_N_its = atoi(optarg);      break;    case 'a':      robust_alpha = atof(optarg);      break;    case '?':      errflag = 1;      break;          }      /* command line error check */    if (errflag || !inflag || !eigflag || !outflag) {    fprintf(stderr,"\nUSAGE: %s %s\n%s\n", progname, usage, help);    exit(1);  }    /* ----  read Eigenimage DAT-file parameters and eigenvectors -------- */    /* read_descriptor(eigdir, &neigframe, &sets, &bytes_pixel, &ncol, &nrow); */  read_descriptor2(eigdir, &neigframe, &sets, &channels, 		   &bytes_pixel, &ncol, &nrow);  if (sets>1)     myerror("Eigenvector files must be single-set DAT files!");  if (neigframe_select==1) {    if (neigframe_in>neigframe)      myerror("Can't use more Eigenimages than available!");    neigframe = neigframe_in;  }  eigimages = f3tensor(0, neigframe-1, 1, nrow*channels, 1, ncol);  for (f=0; f<neigframe; f++) {    if (bytes_pixel==1)      read_DAT_frame(eigimages[f], eigdir, f, neigframe, nrow*channels, ncol);    if (bytes_pixel==4) {      sprintf(filename,"%s/data%d", eigdir, f);      read_RAW_float(filename, eigimages[f], nrow*channels, ncol);    }  }  fprintf(stdout,	  "Read %3d  %d-by-%d  %d-channel %s eigenvector files from %s/\n\n",	  neigframe, nrow, ncol, channels, 	  (bytes_pixel==4 ? "float":"uchar"), eigdir);  coeffs = vector(0, neigframe);  image = matrix(1, nrow*channels, 1, ncol);    /* ----  read indir descriptor file -------- */  /*read_descriptor(indir, &nframe, &sets, &bytes_pixel,&in_ncol,&in_nrow); */  if (desc_flag == 0)    read_descriptor2(indir, &nframe, &in_sets, &in_channels, 		     &bytes_pixel, &in_ncol, &in_nrow);  else    read_descriptor2(descdir, &nframe, &in_sets, &in_channels, 		     &bytes_pixel, &in_ncol, &in_nrow);   if (nrow!=in_nrow || ncol!=in_ncol) {    fprintf(stderr,	    "ERROR: Input dimensions %d-by-%d don't match eigenvector dimensions %d-by-%d.\n",	    in_nrow, in_ncol, nrow, ncol);    exit(1);  }  if (channels != in_channels) {    fprintf(stderr,	    "ERROR: Input channels (%d) don't match eigenvectors' (%d)\n", 	    channels, in_channels);    exit(2);  }  if (in_sets>1)     myerror("Input files must be single-set DAT files!");    if (residual_select>0)    test = matrix(1, nrow*channels, 1, ncol);  if (robust_select>0) {        robust_Ndim   = nrow*ncol*channels;    robust_Phi = matrix(1, robust_Ndim, 1, neigframe-1);    robust_X   = vector(1, robust_Ndim);    robust_E   = vector(1, robust_Ndim);    robust_W   = vector(1, robust_Ndim);    /* create transposed eigenvector matrix */    for (k=1; k<neigframe; k++) {      l = 1;      for (i=1; i<=nrow*channels; i++)	for (j=1; j<=ncol; j++)	  robust_Phi[l++][k] = eigimages[k][i][j];    }  }  /* ---- loop over input list file and warp ------- */  if ((fp1 = 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, fp1)) {    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {            nframe++;            /* ------- get filename ------- */             sscanf(line,"%s",infile);        sprintf(filename,"%s/%s", indir, infile);      if (bytes_pixel==4)	read_RAW_float(filename, image, nrow*channels, ncol);      else	read_RAW(filename, image, nrow*channels, ncol);      fprintf(stdout,"Read %d-by-%d  %d-channel %s input file %s\n",	      nrow, ncol, channels, (bytes_pixel==4 ? "float":"uchar"),	      filename);                  /* -------- compute projection coefficients ------ */            if (robust_select>0) {	      	l = 1;	for (i=1; i<=nrow*channels; i++)	  for (j=1; j<=ncol; j++)	    robust_X[l++] = image[i][j] - eigimages[0][i][j]; 	for (i=0; i<neigframe; i++)	  coeffs[i] = 0;	fval = robust_project(robust_X, robust_Ndim, robust_Phi,			      neigframe-1, coeffs, robust_sigma,			      robust_N_its, robust_alpha);	fprintf(stdout,"Performed %3d robust iterations,  Norm = %f\n",		robust_N_its, fval);       }            else {   /* --- standard projection --- */	for (i=1; i<=nrow*channels; i++) 	  for (j=1; j<=ncol; j++) 	    image[i][j] -= eigimages[0][i][j];		for (k=1; k<neigframe; k++) {	  fval = 0.0;	  for (i=1; i<=nrow*channels; i++)	    for (j=1; j<=ncol; j++) 	      fval += (image[i][j] * eigimages[k][i][j]);	  coeffs[k] = fval;	}      }      fprintf(stdout,"Computed %d coefficients, ", neigframe-1);                  /* --------- DEBUG:  reconstruction -------------- */      if (residual_select>0) {		for (i=1; i<=nrow*channels; i++) 	  for (j=1; j<=ncol; j++) {	    fval = 0.0;	    for (k=1; k<neigframe; k++) 	      fval += (eigimages[k][i][j] * coeffs[k]);	    test[i][j] = fval;	  }	error = 0.0; 	for (i=1; i<=nrow*channels; i++) 	  for (j=1; j<=ncol; j++) { 	    error += SQR(test[i][j] - image[i][j]);	  }	fprintf(stdout,"residual = %g  \n", error/(nrow*channels*ncol));      }      else	fprintf(stdout,"\n");                  /* ---- write results to output directory -------- */            sprintf(filename,"%s/%s", outdir, infile);          if ( (fp2 = fopen(filename, "w")) == NULL)	myerror("could not open output file.");      for (i=1; i<neigframe; i++)	fprintf(fp2, "%f\n", coeffs[i]);      fclose(fp2);      fprintf(stdout,"Wrote %d-by-1 ASCII coefficient file to %s\n\n", 	      neigframe-1, filename);    }  }  /* end of main processing loop */    fclose(fp1);  /* ----  open output data file in outdir ---- */  /* write_descriptor(outdir, nframe, 1, neigframe-1, 0, comline); */  write_descriptor2(outdir, nframe, 1, 1, 1, neigframe-1, 0, comline);    /* --- free up allocated matrices ----- */  free_f3tensor(eigimages, 0, neigframe-1, 1, nrow*channels, 1, ncol);  free_matrix(image, 1, nrow*channels, 1, ncol);  if (residual_select>0)    free_matrix(test, 1, nrow*channels, 1, ncol);  free_vector(coeffs, 0, neigframe);  if (robust_select>0) {    free_matrix(robust_Phi, 1, robust_Ndim, 1, neigframe-1);    free_vector(robust_X, 1, robust_Ndim);    free_vector(robust_E, 1, robust_Ndim);    free_vector(robust_W, 1, robust_Ndim);  }    return 0;  }/* -----------END OF MAIN ----------------------------------------- */float robust_project(float *X, int N, float **Phi, int M,		    float *C, float sigma, int N_its, float alpha){  register int n,i,j,k;  float sum;  float s2 = sigma*sigma;  float J;  for (n=1; n<=N_its; n++)  {    /* compute error:  E = X - Phi*C */        for (i=1; i<=N; i++) {      sum = 0.0;      for (k=1; k<=M; k++)	sum += Phi[i][k] * C[k];      robust_E[i] = X[i] - sum;    }    /* compute cost */        sum = 0;    for (i=1; i<=N; i++)       sum += SQR(robust_E[i]) / (s2 + SQR(robust_E[i]));    J = sum;    fprintf(stderr,"robust_project(): J(%d) = %f\n",n,J);    /* compute weights:  W = 2*E*sigma^2 ./ (sigma^2 + E.^2).^2 */    for (i=1; i<=N; i++)       robust_W[i] = 2*robust_E[i]*s2 / SQR((s2 + SQR(robust_E[i])));    /* update coefficients:  C = C + alpha * Phi' * W    */         for (i=1; i<=M; i++) {      sum = 0.0;      for (j=1; j<=N; j++)	sum += (Phi[j][i] * robust_W[j]);      C[i] += alpha * sum;    }  }  /* -- end of iteration loop -- */   return J;}