#include <stdio.h>

 #include <string.h>

 #include "cv.h"

 #include "cvaux.h"

 #include "highgui.h"

 using namespace cv;

 //globle variables

 int nTrainFaces         = ;    // number of trainning images

 int nEigens             = ;    // number of eigenvalues

 IplImage** faceImgArr         = ;     // array of face images

 CvMat* personNumTruthMat     = ;    // array of person numbers

 IplImage* pAvgTrainImg         = ;    // the average image

 IplImage** eigenVectArr     = ;    // eigenvectors

 CvMat* eigenValMat         = ;    // eigenvalues

 CvMat* projectedTrainFaceMat     = ;    // projected training faces

 //// Function prototypes

 void learn();

 void recognize();

 void doPCA();

 void storeTrainingData();

 int loadTrainingData(CvMat** pTrainPersonNumMat);

 int findNearestNeighbor(float* projectedTestFace);

 int loadFaceImgArray(char* filename);

 void printUsage();

 int main( int argc, char** argv )

 {

     if((argc != ) && (argc != )){

     printUsage();

     return -;

     }

     if( !strcmp(argv[], "train" )){

     learn();

     } else if( !strcmp(argv[], "test") ){

     recognize();

     } else {

     printf("Unknown command: %s\n", argv[]);

     }

     return ;

 }

 void printUsage(){

     printf("Usage: eigenface <command>\n",

     "  Valid commands are\n"

     "    train\n"

     "    test\n"

     );

 }

 void learn(){

     int i;

     // load training data

     nTrainFaces = loadFaceImgArray("train.txt");

     if( nTrainFaces < ){

     fprintf(

         stderr,

         "Need 2 or more training faces\n"

           "Input file contains only %d\n",

         nTrainFaces

     );

     return;

     }

     // do PCA on the training faces

     doPCA();

     // project the training images onto the PCA subspace

     projectedTrainFaceMat = cvCreateMat(nTrainFaces, nEigens, CV_32FC1);

     for(i = ; i < nTrainFaces; i ++){

     cvEigenDecomposite(

         faceImgArr[i],

         nEigens,

         eigenVectArr,

         , ,

         pAvgTrainImg,

         projectedTrainFaceMat->data.fl + i*nEigens

     );

     }

     // store the recognition data as an xml file

     storeTrainingData();

 }

 int loadFaceImgArray(char* filename){

     FILE* imgListFile = ;

     char imgFilename[];

     int iFace, nFaces = ;

     // open the input file

     imgListFile = fopen(filename, "r");

     // count the number of faces

     while( fgets(imgFilename, , imgListFile) ) ++ nFaces;

     rewind(imgListFile);

     // allocate the face-image array and person number matrix

     faceImgArr = (IplImage **)cvAlloc( nFaces*sizeof(IplImage *) );

     personNumTruthMat = cvCreateMat( , nFaces, CV_32SC1 );

     // store the face images in an array

     for(iFace=; iFace<nFaces; iFace++){

     //read person number and name of image file

       fscanf(imgListFile, "%d %s", personNumTruthMat->data.i+iFace, imgFilename);

     // load the face image

     faceImgArr[iFace] = cvLoadImage(imgFilename, CV_LOAD_IMAGE_GRAYSCALE);

     }

     fclose(imgListFile);

     return nFaces;

 }

 void doPCA(){

     int i;

     CvTermCriteria calcLimit;

     CvSize faceImgSize;

     // set the number of eigenvalues to use

     nEigens = nTrainFaces - ;

     // allocate the eigenvector images

     faceImgSize.width = faceImgArr[]->width;

     faceImgSize.height = faceImgArr[]->height;

     eigenVectArr = (IplImage**)cvAlloc(sizeof(IplImage*) * nEigens);

     for(i=; i<nEigens; i++){

     eigenVectArr[i] = cvCreateImage(faceImgSize, IPL_DEPTH_32F, );

     }

     // allocate the eigenvalue array

     eigenValMat = cvCreateMat( , nEigens, CV_32FC1 );

     // allocate the averaged image

     pAvgTrainImg = cvCreateImage(faceImgSize, IPL_DEPTH_32F, );

     // set the PCA termination criterion

     calcLimit = cvTermCriteria( CV_TERMCRIT_ITER, nEigens, );

     // compute average image, eigenvalues, and eigenvectors

     cvCalcEigenObjects(

     nTrainFaces,

     (void*)faceImgArr,

     (void*)eigenVectArr,

     CV_EIGOBJ_NO_CALLBACK,

     ,

     ,

     &calcLimit,

     pAvgTrainImg,

     eigenValMat->data.fl

     );

 }

 void storeTrainingData(){

     CvFileStorage* fileStorage;

     int i;

     // create a file-storage interface

     fileStorage = cvOpenFileStorage( "facedata.xml", , CV_STORAGE_WRITE);

     // store all the data

     cvWriteInt( fileStorage, "nEigens", nEigens);

     cvWriteInt( fileStorage, "nTrainFaces", nTrainFaces );

     cvWrite(fileStorage, "trainPersonNumMat", personNumTruthMat, cvAttrList(, ));

     cvWrite(fileStorage, "eigenValMat", eigenValMat, cvAttrList(,));

     cvWrite(fileStorage, "projectedTrainFaceMat", projectedTrainFaceMat, cvAttrList(,));

     cvWrite(fileStorage, "avgTrainImg", pAvgTrainImg, cvAttrList(,));

     for(i=; i<nEigens; i++){

     char varname[];

     sprintf( varname, "eigenVect_%d", i);

     cvWrite(fileStorage, varname, eigenVectArr[i], cvAttrList(,));

     }

     //release the file-storage interface

     cvReleaseFileStorage( &fileStorage );

 }

 void recognize(){

     int i, nTestFaces = ;        // the number of test images

     CvMat* trainPersonNumMat = ;    // the person numbers during training

     float* projectedTestFace = ;    

     // load test images and ground truth for person number

     nTestFaces = loadFaceImgArray("test.txt");

     printf("%d test faces loaded\n", nTestFaces);

     // load the saved training data

     if( !loadTrainingData( &trainPersonNumMat ) ) return;

     // project the test images onto the PCA subspace

     projectedTestFace = (float*)cvAlloc( nEigens*sizeof(float) );

     for(i=; i<nTestFaces; i++){

     int iNearest, nearest, truth;

      // project the test image onto PCA subspace

     cvEigenDecomposite(

         faceImgArr[i],

         nEigens,

         eigenVectArr,

             , ,

         pAvgTrainImg,

         projectedTestFace

     );

     iNearest = findNearestNeighbor(projectedTestFace);

     truth = personNumTruthMat->data.i[i];

     nearest = trainPersonNumMat->data.i[iNearest];

     printf("nearest = %d, Truth = %d\n", nearest, truth);

     }

 }

 int loadTrainingData(CvMat** pTrainPersonNumMat){

     CvFileStorage* fileStorage;

     int i;

     // create a file-storage interface

     fileStorage = cvOpenFileStorage( "facedata.xml", , CV_STORAGE_READ );

     if( !fileStorage ){

     fprintf(stderr, "Can't open facedata.xml\n");

     return ;

     }

     nEigens = cvReadIntByName(fileStorage, , "nEigens", );

     nTrainFaces = cvReadIntByName(fileStorage, , "nTrainFaces", );

     *pTrainPersonNumMat = (CvMat*)cvReadByName(fileStorage, , "trainPersonNumMat", );

     eigenValMat = (CvMat*)cvReadByName(fileStorage, , "eigenValMat", );

     projectedTrainFaceMat = (CvMat*)cvReadByName(fileStorage, , "projectedTrainFaceMat", );

     pAvgTrainImg = (IplImage*)cvReadByName(fileStorage, , "avgTrainImg", );

     eigenVectArr = (IplImage**)cvAlloc(nTrainFaces*sizeof(IplImage*));

     for(i=; i<nEigens; i++){

     char varname[];

     sprintf( varname, "eigenVect_%d", i );

     eigenVectArr[i] = (IplImage*)cvReadByName(fileStorage, , varname, );

     }

     // release the file-storage interface

     cvReleaseFileStorage( &fileStorage );

     return ;

 }

 int findNearestNeighbor(float* projectedTestFace){

     double leastDistSq = DBL_MAX;

     int i, iTrain, iNearest = ;

     for(iTrain=; iTrain<nTrainFaces; iTrain++){

     double distSq = ;

     for(i=; i<nEigens; i++){

         float d_i = projectedTestFace[i] -

         projectedTrainFaceMat->data.fl[iTrain*nEigens + i];

         distSq += d_i*d_i;

     }

     if(distSq < leastDistSq){

         leastDistSq = distSq;

         iNearest = iTrain;

     }

     }

     return iNearest;

 }