/*
 This Work by Takashi OHKUBO is licensed under a Creative 
 Commons Attribution-NonCommercial 2.1 Japan License. 
 Using this work in advertising is allowed. -21st/June/2008-
*/
#include <stdio.h>
#include <stdlib.h>
#include <cv.h>

#define DATA_FILE "data.txt"
#define DATA_NUM 10000
#define DATA_DIM 3
#define PCA_DIM  1

void readData(CvMat *source);
void writeData(CvMat *project, CvMat *result);
void plotData(CvMat *result);
void displayPCA_Products( CvMat *average, 
			  CvMat *eigenValues, CvMat *eigenVectors);

int main(void){
  // CREATE MAT
  CvMat *source,*project,*result;
  CvMat *average,*eigenValues,*eigenVectors;
  source  = cvCreateMat(DATA_NUM,DATA_DIM,CV_32F);  
  result  = cvCreateMat(DATA_NUM,DATA_DIM,CV_32F);
  project = cvCreateMat(DATA_NUM,PCA_DIM,CV_32F);  
  average = cvCreateMat(1,DATA_DIM,CV_32F);
  eigenValues  = cvCreateMat(1,PCA_DIM,CV_32F);
  eigenVectors = cvCreateMat(PCA_DIM,DATA_DIM,CV_32F);
  // READ DATA
  readData(source);  
  // PCA
  cvCalcPCA(source,average,eigenValues,eigenVectors,CV_PCA_DATA_AS_ROW);
  cvProjectPCA(source,average,eigenVectors,project);
  cvBackProjectPCA(project,average,eigenVectors,result);
  // DISPLAY & WRITE LOG
  displayPCA_Products(average,eigenValues,eigenVectors);
  writeData(project,result);
  plotData(result);
  // RELEASE MAT
  cvReleaseMat(&source);
  cvReleaseMat(&project);
  cvReleaseMat(&result);
  cvReleaseMat(&average);
  cvReleaseMat(&eigenValues);
  cvReleaseMat(&eigenVectors);
  
  return 0;
}

void readData(CvMat *source){
  int i,j;
  double x;
  FILE *fp=fopen(DATA_FILE,"rt");
  for(i=0;i<DATA_NUM;i++)
    for(j=0;j<DATA_DIM;j++){
      fscanf(fp,"%lf",&x);
      cvmSet(source,i,j,x);
    }
  fclose(fp);
}

void writeData(CvMat *project, CvMat *result){
  int i,j;
  char fileName[200];
  FILE *fp;
  sprintf(fileName,"project_%s",DATA_FILE);
  fp=fopen(fileName,"w");
  for(i=0;i<DATA_NUM;i++){
    for(j=0;j<PCA_DIM;j++)
      fprintf(fp,"%lf\t",cvmGet(project,i,j));
    fprintf(fp,"\n");
  }
  fclose(fp);
  sprintf(fileName,"result_%s",DATA_FILE);
  fp=fopen(fileName,"w");
  for(i=0;i<DATA_NUM;i++){
    for(j=0;j<DATA_DIM;j++)
      fprintf(fp,"%lf\t",cvmGet(result,i,j));
    fprintf(fp,"\n");
  }
  fclose(fp);
}

void plotData(CvMat *result){
  int i,j;
  FILE *gp=popen("gnuplot -geometry 640x480","w");
  fprintf(gp,"set grid\nset ticslevel 0\n");
  fprintf(gp,"splot '%s' u 1:2:3 t 'source'  w d, '-' t 'ReBuild(PCA(source))'  w d 3\n",DATA_FILE);
  for(i=0;i<DATA_NUM;i++)
    fprintf(gp,"%lf\t%lf\t%lf\n",cvmGet(result,i,0),
	    cvmGet(result,i,1),cvmGet(result,i,2));
  fprintf(gp,"e\n");
  fflush(gp);
  printf("\n> Press any key to quit this program.\n");
  getchar();
  fclose(gp);
}

void displayPCA_Products( CvMat *average, 
		      CvMat *eigenValues, CvMat *eigenVectors){
  int i,j;
  char fileName[200];
  FILE *fp;
  sprintf(fileName,"PCA_products_%s",DATA_FILE);
  fp=fopen(fileName,"w");
  printf("\n> Average\n");
  fprintf(fp,"\n> Average\n");
  for(i=0;i<DATA_DIM;i++){
    printf("DIMENSION[%d]\t%lf\n",i,cvmGet(average,0,i));
    fprintf(fp,"DIMENSION[%d]\t%lf\n",i,cvmGet(average,0,i));
  }
  printf("\n> Eigen Vectors\n");
  fprintf(fp,"\n> Eigen Vectors\n");
  for(i=0;i<PCA_DIM;i++){
    printf("PC[%d]\t",i);
    fprintf(fp,"PC[%d]\t",i);
    for(j=0;j<DATA_DIM;j++){
      printf("%lf\t",cvmGet(eigenVectors,i,j));
      fprintf(fp,"%lf\t",cvmGet(eigenVectors,i,j));
    }
    printf("\n");
    fprintf(fp,"\n");
  }
  printf("\n> Eigen Vectors\n");
  fprintf(fp,"\n> Eigen Vectors\n");
  for(i=0;i<PCA_DIM;i++){
    printf("PC[%d]\t%lf \n",i,cvmGet(eigenValues,0,i));
    fprintf(fp,"PC[%d]\t%lf \n",i,cvmGet(eigenValues,0,i));
  }
}