Main Page   File List   File Members  

PCA_PP.cpp

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*/
00015 /*---------------------------------------------------------------------------*/
00016 #include <iostream>
00017 #include <string>  
00018 #include <vector>
00019 #include <fstream>
00020 
00021 #include <cstdio>
00022 #include <cstdlib>
00023 #include <cmath>
00024 #include <cassert>
00025 #include <ctime>
00026 
00027 //format of the input images
00028 #define PIX float
00029 #define inf 1e20
00030 
00031 using namespace std;
00032 
00035 void usage(string argv0) {
00036 cout <<"*************************************************************************************************************"<<endl;
00037     cout << "\nusage: " << ">"
00038               << argv0 
00039               << " input_image_file_name output_PPimage_file_name output_PCAimage_file_name nb_lignes nb_colonnes PCAEnergiePercentage PPquantile" << endl; 
00040     cout << "\nOptions for cutting: " 
00041          << "offset_colonnes offset_lignes nb_colonnes_original"<< endl;
00042          
00043     cout << "\nOptional parameters for Projection Pursuit : nb rand_departure, nb_maximum_of_iterations, epsilon,  factor_decreasing" <<endl;   
00044    
00045     cout << "\n Default parameters:\n nb rand_departure = 5\n nb_maximum_of_iterations = 50\n epsilon =1e-3\n factor_decreasing = 2\n quantile = 0.7\n PCAEnergiePercentage= 90\n"<<endl;
00046  
00047  cout <<"*************************************************************************************************************"<<endl;
00048   
00049 }
00050 
00051 //auxiliary functions
00056 float Evaluate_PI(float **x, unsigned long int lenght);
00057 
00062 void base_ortho_plan_finder(float **plan_init,float **new_base,unsigned long int nb_dim);
00063 
00069 void projec_on_base(float **Y,float **newbase,unsigned long int dim, unsigned long int nb_dim);
00070 
00076 void projec_on_oldbase(float **Y,float **newbase,unsigned long int dim, unsigned long int nb_dim);
00077 
00081 void destructuration(float **Y,unsigned long int dim);
00082 
00086 double Total(float *a, unsigned long int lenght);
00087 
00088 //numerical recepies
00090 void eigsrt(float d[], float **v, int n);
00092 void jacobi(float **a, int n, float d[], float **v, int *nrot);
00094 float gasdev(long *idum); 
00096 void gaussj(float **a, int n, float **b, int m);
00097         
00127 int main(int argc, char *argv []){   
00128   
00129 srand((unsigned)time(NULL));
00130 
00131 string ifn1;
00132 string ifn2;
00133 string ifn3;
00134 
00135 unsigned long int           dimx           = 0;
00136 unsigned long int           dimy           = 0;
00137 int                         ofx            = 0;
00138 int                         ofy            = 0;
00139 int                         Dy             = 0;
00140 float                       quantile       = 0.7;
00141 int                         keepEnergie    = 90;
00142 //Default parameters....
00143 
00144 float eigenthreshold=0.000001;
00145 int nb_rand_departure=5;//nb of random departure 
00146 float cbegin=tan(80*M_PI/180);//c initial value
00147 float eps=1e-3;//limit of c
00148 float factor=2;//factor which devides c[1,2,3] for each iteration
00149 float factor_cons=1;//factor which devides c[0] only once
00150 int nb_step, nb_max_step=50;//nb max of iteration before c=c\2
00151 
00152 switch (argc-1) {
00153 
00154          case 14:  factor               = atof(argv[14]);
00155          case 13:  eps                  = atof(argv[13]);
00156          case 12:  nb_max_step          = atoi(argv[12]);
00157          case 11:  nb_rand_departure    = atoi(argv[11]);
00158          case 10:  Dy                   = atoi(argv[10]);
00159          case 9:  ofy                   = atoi(argv[9]);
00160          case 8:  ofx                   = atoi(argv[8]);
00161          case 7:  quantile              = atof(argv[7]);
00162          case 6:  keepEnergie           = atoi(argv[6]);
00163          case 5:  dimy                  = atoi(argv[5]);
00164          case 4:  dimx                  = atoi(argv[4]);
00165          case 3:  ifn3                  = string(argv[3]);
00166          case 2:  ifn2                  = string(argv[2]);
00167          case 1:  ifn1                  = string(argv[1]); break;
00168          case 0:
00169         default: usage(argv[0]); return -1;
00170     }
00171     if ((argc-1) < 5) { usage(argv[0]); return -1;}
00172     if ((argc-1) == 5) {Dy=dimy;} 
00173      cout<<"\n*************************************************************************************************************"<<endl;cout<<"Multidimensionnal analysis and dimension reduction by Principal Component Analysis and Projection Pursuit\n
00174      Author : Patrick HEAS  2005"<<endl;cout<<"*************************************************************************************************************\n"<<endl;
00175     cout << "> " << argv[0] 
00176               << " " << ifn1 
00177               << " " << ifn2
00178               << " " << ifn3 
00179               << " " << dimx 
00180               << " " << dimy 
00181               << " " << keepEnergie 
00182               << " " << quantile
00183               << " " << ofx
00184               << " " << ofy  
00185               << " " << Dy      
00186               << " " << nb_rand_departure 
00187               << " " << nb_max_step 
00188               << " " << eps
00189               << " " << factor << endl;
00190         
00191 
00192 
00193 /*************************************************************************************************************/
00194 
00195 time_t tdeb, t1,t2;
00196 time(&tdeb);
00197 t1=tdeb;
00198 
00199 
00200 unsigned long  int i=0,j=0,k=0,l=0,m=0,n=0;;
00201 unsigned long  int dim=dimx*dimy;
00202  
00203 unsigned long int nb_bands=1;
00204 
00205 string line;
00206 ifstream inaux;
00207 inaux.open(ifn1.c_str(),ifstream::in);
00208 ifstream in;
00209 in.open(ifn1.c_str(),ifstream::in);
00210         
00211 n=0;
00212 while(getline(inaux,line)) {n++;}       
00213 unsigned long  int nb_images=n, nb_dim=n;
00214 cout << "number of images : " << nb_images << std::endl;
00215 
00216 char ***path= new char **[nb_images];
00217 for(m=0;m<nb_images;++m){path[m] =new char*[1];} 
00218 for(m=0;m<nb_images;++m){for(k=0;k<1;++k){path[m][k]=new char[100];}}
00219 char ***save_path= new char **[(int)nb_images/2];
00220 for(m=0;m<(unsigned long int)nb_images/2;++m){save_path[m] =new char*[2];} 
00221 for(m=0;m<(unsigned long int)nb_images/2;++m){for(k=0;(int)k<2;++k){save_path[m][k]=new char[100];}}
00222 char ***save_path2= new char **[nb_images];
00223 for(m=0;m<nb_images;++m){save_path2[m] =new char*[1];} 
00224 for(m=0;m<nb_images;++m){for(k=0;k<1;++k){save_path2[m][k]=new char[100];}}
00225 
00226 /*******************************Reading text file names*********************************************/
00227 
00228 for (i=0;i<(unsigned long int)nb_images;i++){
00229         for (j=0;j<(unsigned long int)nb_bands;j++){
00230                         getline(in,line);
00231                         strcpy(path[i][j],(char *)line.c_str());
00232         }       
00233 }
00234 
00235 string line1;
00236 ifstream inaux1;
00237 inaux1.open(ifn2.c_str(),ifstream::in);
00238 ifstream in1;
00239 in1.open(ifn2.c_str(),ifstream::in);
00240         
00241 cout <<" PP Save Files :"<<endl;                                
00242 for (i=0;i<(unsigned long int)nb_images/2;i++){
00243         for (j=0;j<2;j++){
00244                         getline(in1,line1);
00245                         strcpy(save_path[i][j],(char *)line1.c_str());
00246                         std::cout<<save_path[i][j]<<std::endl;
00247         }       
00248 }
00249 string line2;
00250 ifstream inaux2;
00251 inaux2.open(ifn3.c_str(),ifstream::in);
00252 ifstream in2;
00253 in2.open(ifn3.c_str(),ifstream::in);
00254 
00255 cout<<" PCA Save Files :"<<endl;                                
00256 for (i=0;i<(unsigned long int)nb_images;i++){
00257         for (j=0;j<(unsigned long int)nb_bands;j++){
00258                         getline(in2,line2);
00259                         strcpy(save_path2[i][j],(char *)line2.c_str());
00260                         std::cout<<save_path2[i][j]<<std::endl;
00261         }       
00262 }
00263 
00264 
00265 /*******************************Reading the images*********************************************/
00266 cout<<"\n";
00267 printf("Reading images\n");
00268  PIX a[dimx][dimy]    ;                                            
00269        float **clust= new float *[nb_dim]; //vector clust contains the data
00270         for(m=0;m<nb_dim;++m)
00271              clust[m] =new float[dim];
00272 
00273          char * ifn;
00274          FILE * ifp;
00275          m=0;n=0;
00276          for (l = 0 ;l <nb_dim; l++){
00277            ifn=path[m][n];printf("%s\n",ifn); // Reading  dimensions R, G and B
00278            ifp=NULL;
00279            assert (ifp=fopen(ifn,"r"));
00280            assert (0==fseek(ifp,(sizeof(PIX)*Dy*ofy+sizeof(PIX)*ofx), SEEK_SET));
00281            assert (dimy==fread(&a[0][0],sizeof (PIX),dimy,ifp));    
00282            for (i=1;i<dimx;i++){
00283              assert (0==fseek(ifp,(sizeof(PIX)*(Dy-dimy)), SEEK_CUR));
00284              assert (dimy==fread(&a[i][0],sizeof (PIX),dimy,ifp));
00285            }
00286            fclose(ifp); 
00287            if (n==nb_bands-1){n=0;m++;}
00288            else{n++;}
00289            for (i = 0 ;i < (dimx); i++){ k=(dimy)*i;          
00290              for (j = 0; j <(dimy); j++) { 
00291                clust[l][k]=a[i][j];
00292                k=k+1;
00293              }
00294            }
00295          }  
00296         
00297 /*********************************Spherisation of the data*************************************************/
00298 cout<<"\nDecorralating procedure"<<endl;        
00299 
00300 float **Y= new float *[nb_dim]; //    Sphered data 
00301 for(m=0;m<nb_dim;++m)
00302 { Y[m]=new float[dim];}
00303 float **Y1= new float *[nb_dim]; //    Sphered data bis
00304 for(m=0;m<nb_dim;++m)
00305 { Y1[m]=new float[dim];}
00306  float *sauvACP= new float [dim];//vector to save
00307 float *moy= new float [nb_dim]; //    mean vector 
00308 float *M= new float [nb_dim];  //   Normalizing diag Matrix 
00309 float *eigenvalues= new float [nb_dim+1];//eigen values (index belong to [1;nb_dim])
00310 float **cov= new float *[nb_dim+1];//    covariance Matrix (index belong to [1;nb_dim][1;nb_dim])
00311 float **eigenvectors= new float *[nb_dim+1];//eigen vectors(index belong to [1;nb_dim][1;nb_dim])
00312 for(m=0;m<nb_dim+1;++m)
00313 {cov[m]=new float[nb_dim+1]; 
00314  eigenvectors[m]=new float[nb_dim+1]; }
00315 
00316 int *nrot= new int [1];
00317 double truc;
00318 
00319 
00320 //printf("Mean vector\n");
00321 for (l=0;l<nb_dim;l++)
00322 {  truc=Total(clust[l],dim);
00323    moy[l]=(float)(truc/dim);
00324    M[l]=0;   
00325    for (i=0;i<dim;i++){
00326       M[l]=((clust[l][i]-moy[l])*(clust[l][i]-moy[l]))+M[l];
00327    }
00328    M[l]= sqrt(1/((double)dim)*M[l]);
00329 } 
00330 
00331 //printf("Covariance Matrix");
00332 for (l=0;l<nb_dim;l++){//printf("\n");
00333   for (k=l;k<nb_dim;k++){
00334     cov[l+1][k+1]=0;
00335     for (i=0;i<dim;i++){
00336       cov[l+1][k+1]=(clust[l][i]-moy[l])*(clust[k][i]-moy[k])+cov[l+1][k+1];
00337     }
00338     cov[l+1][k+1]=cov[l+1][k+1]/(dim*M[l]*M[k]);
00339     cov[k+1][l+1]=cov[l+1][k+1];
00340     //printf("%f\t",cov[l+1][k+1]);
00341   }
00342 } 
00343 
00344 jacobi(cov,nb_dim,eigenvalues,eigenvectors,nrot);//compute eigensysteme
00345 eigsrt(eigenvalues,eigenvectors,nb_dim);//sort eigen values and vectors
00346 
00347  
00348 /*Keeping eigen values representing more than "keepEnergie" pourcent of the signal energie*/
00349 float TotalEnergie=0;
00350 for (l=1;l<nb_dim+1;l++){
00351         TotalEnergie+=eigenvalues[l];
00352 }
00353 float Energie=0;
00354 for (l=1;l<nb_dim+1;l++){
00355         Energie+=eigenvalues[l];
00356         if (Energie/TotalEnergie*100>keepEnergie){eigenthreshold=eigenvalues[l];break;}
00357 }
00358 
00359 int dim_less=0;
00360 for (l=1;l<nb_dim+1;l++)
00361   {
00362   if (eigenvalues[l]<=eigenthreshold){ dim_less++; }
00363   }
00364 
00365 nb_dim=nb_dim-dim_less;// Component not very descriminent which have a eigenvalue under eigenthreshold
00366 cout<<"Number of component representing "<<keepEnergie<<" of the signal energie : "<<nb_dim<<endl; 
00367                       // These dimensions are removed and we work now with b_dim-dim_less dimensions
00368 
00369 for (i=0;i<dim;i++){
00370   for (l=0;l<nb_dim;l++){ Y[l][i]=0;
00371     for (k=0;k<nb_dim;k++){
00372       Y[l][i]= Y[l][i]+(clust[k][i]-moy[k])/M[k]*eigenvectors[k+1][l+1];
00373     }
00374     Y[l][i]=Y[l][i]/(sqrt(eigenvalues[l+1]));Y1[l][i]=Y[l][i];
00375   }
00376 }
00377 
00378 printf("saving PCA results\n");
00379 for (j=0;j<nb_dim;j++){
00380   for (i=0;i<dim;i++){sauvACP[i]=(float)Y[j][i];} 
00381   ifp=NULL; 
00382   ifn=save_path2[j][0];
00383   assert (ifp=fopen(ifn,"wb")); 
00384   assert ( dim==fwrite(&sauvACP[0], sizeof(float),dim,ifp)); 
00385   assert (EOF!=fclose(ifp));
00386 }
00387 
00388 for(m=0;m<nb_dim;++m)
00389             { delete [] cov[m];
00390               delete [] eigenvectors[m] ;
00391             }
00392 //delete [] cov ;
00393 delete [] eigenvectors;
00394 //delete [] moy;
00395 //delete [] M;
00396 delete [] eigenvalues;
00397 delete [] sauvACP;
00398 
00399 /*********************************************************************************************************/
00400 /*********************************Global Optimisation **************************************************/
00401 
00402 cout<<"\nInitialisation of optimisation procedure"<<endl;
00403 int nb_proj=0;//number of random departure
00404 int rand_dep;
00405 float *temp1 = new float[4];float *temp2 = new float[4];
00406 float tmp,aux1,aux2, p_i, p_i_old,p_i_opt;
00407 float *c = new float[4];//factor fixing the size of the neighberhood visited
00408 float **new_base = new float*[nb_dim];//orthonormal base ortho gonal to the plan for structure removal
00409 for (m=0;m<nb_dim;++m){new_base[m]= new float[nb_dim];}
00410 float **v_random = new float*[4];
00411 for (m=0;m<4;++m){v_random[m]= new float[nb_dim];}
00412 float **plan_init= new float*[2];
00413 for (m=0;m<2;++m){ plan_init[m]= new float[nb_dim];} 
00414 float **ab1= new float*[8];    
00415 float **ab2= new float*[8];
00416 for (m=0;m<8;++m){  ab1[m]= new float[nb_dim]; ab2[m]= new float[nb_dim];}
00417 float **projec_plan= new float*[8];
00418 for (m=0;m<8;++m){ projec_plan[m]= new float[dim]; }
00419 float **arg= new float*[2];
00420 for (m=0;m<2;++m){ arg[m]= new float[dim]; }
00421 float *sauv= new float [2*dim];//vector to save
00422 
00423 int sub,skip=0;
00424 p_i=1000;
00425 p_i_opt=p_i;
00426 
00427 /*****************************************New projection research*********************************/
00428 printf("Optimisation procedure\n");
00429 while ((p_i_opt>=quantile) && (nb_proj<((int)(nb_images/2)))){
00430 
00431 nb_proj++ ;
00432 printf("%s %d\n","index of current significant projection",nb_proj);
00433 p_i_opt=-1;nb_step=0;
00434 for (rand_dep=0;rand_dep<nb_rand_departure;rand_dep++) // New random departure.........
00435 {
00436 //printf("new random departure\n");
00437  if (skip==0){
00438 //initialisation of a random plan (alpha, beta)
00439 tmp=0;
00440  for (l=0;l<nb_dim;l++){plan_init[0][l]=(((float)rand())/RAND_MAX-0.5);}
00441  for (l=0;l<nb_dim;l++){tmp=plan_init[0][l]*plan_init[0][l]+tmp;}
00442  tmp=sqrt(tmp);
00443  for (l=0;l<nb_dim;l++){
00444  plan_init[0][l]=plan_init[0][l]/tmp;}//unitary condition 
00445 
00446 tmp=0;//beta
00447 int last=(int)((((float)rand())/RAND_MAX)*nb_dim-0.5);
00448 for (l=0;l<nb_dim;l++){plan_init[1][l]=(((float)rand())/RAND_MAX-0.5);} 
00449 for (l=0;l<nb_dim;l++){tmp=plan_init[0][l]*plan_init[1][l]+tmp;}
00450 tmp=tmp-plan_init[0][last]*plan_init[1][last];
00451 plan_init[1][last]=-1.*tmp/(plan_init[0][last]);//orthogonal condition
00452 tmp=0;
00453 for (l=0;l<nb_dim;l++){tmp=plan_init[1][l]*plan_init[1][l]+tmp;}
00454 tmp=sqrt(tmp);
00455 for (l=0;l<nb_dim;l++){plan_init[1][l]=plan_init[1][l]/tmp;}//unitary condition
00456  }
00457 //projection on the plan of the initial solution
00458 
00459 
00460 for (i=0;i<dim;i++) 
00461 {
00462   arg[0][i]=0;arg[1][i]=0;
00463   for (l=0;l<nb_dim;l++)
00464   {   
00465     arg[0][i]=Y[l][i]*plan_init[0][l]+arg[0][i]; 
00466     arg[1][i]=Y[l][i]*plan_init[1][l]+arg[1][i];
00467    }
00468 }
00469 
00470 
00471 p_i=Evaluate_PI(arg,dim);p_i_old=p_i;// chi2
00472 //printf("%s %e\n","Projection indice initial=",p_i);
00473 
00474 //searching the maxima.........................................................................
00475 
00476 // 4 parallele sub_routines.....
00477 for(i=0;i<4;i++){c[i]=cbegin;}//initaialisating at tan(80)=5.67
00478 while (c[1]>eps)
00479   {// printf("%e\n",c[1]);
00480 
00481   //initialisation of 4 random vector v
00482   for(sub=0;sub<4;sub++)
00483   {
00484   tmp=0;
00485   for (l=0;l<nb_dim;l++){ v_random[sub][l]=(((float)rand())/RAND_MAX);}
00486   for (l=0;l<nb_dim;l++){tmp= v_random[sub][l]*v_random[sub][l]+tmp;}
00487   tmp=sqrt(tmp);
00488   for (l=0;l<nb_dim;l++){ v_random[sub][l]= v_random[sub][l]/tmp;}//unitary condition
00489   }
00490 
00491   //calculation of 2 couple of  vectors for the 2 different plans for EACH sub-routine
00492   
00493   temp1[0]=0;temp1[1]=0;temp1[2]=0;temp1[3]=0;
00494   temp2[0]=0;temp2[1]=0;temp2[2]=0;temp2[3]=0;
00495   for (i=0;i<nb_dim;i++)
00496     {
00497     for (sub=0;sub<4;sub++)
00498       {
00499       ab1[sub*2][i]=plan_init[1][i]+c[sub]*v_random[sub][i]; 
00500       ab2[sub*2][i]=plan_init[1][i]-c[sub]*v_random[sub][i];
00501       temp1[sub]=temp1[sub]+ab1[sub*2][i]*ab1[sub*2][i];
00502       temp2[sub]=temp2[sub]+ab2[sub*2][i]*ab2[sub*2][i];
00503       }
00504     }
00505   for (sub=0;sub<4;sub++)  
00506     {
00507       temp1[sub]=sqrt(temp1[sub]);temp2[sub]=sqrt(temp2[sub]); 
00508       for (i=0;i<nb_dim;i++) 
00509         { ab1[2*sub][i]=ab1[2*sub][i]/temp1[sub]; 
00510         ab2[2*sub][i]=ab2[2*sub][i]/temp2[sub];
00511         }
00512     }
00513   
00514    
00515    for (sub=0;sub<4;sub++) 
00516     {
00517     temp1[0]=0;temp1[1]=0;temp1[2]=0;temp1[3]=0;
00518     temp2[0]=0;temp2[1]=0;temp2[2]=0;temp2[3]=0; 
00519     for (i=0;i<nb_dim;i++)
00520       {temp1[sub]=temp1[sub]+ab1[2*sub][i]*plan_init[0][i];
00521       temp2[sub]=temp2[sub]+ab2[2*sub][i]*plan_init[0][i];
00522       }
00523     aux1=0;aux2=0;
00524     for (i=0;i<nb_dim;i++)
00525       {
00526         ab1[2*sub+1][i]=plan_init[0][i]-temp1[sub]*ab1[2*sub][i];
00527         ab2[2*sub+1][i]=plan_init[0][i]-temp2[sub]*ab2[2*sub][i];
00528         aux1=aux1+ab1[2*sub+1][i]*ab1[2*sub+1][i];
00529         aux2=aux2+ab2[2*sub+1][i]*ab2[2*sub+1][i];
00530       }
00531     aux1=sqrt(aux1);
00532     aux2=sqrt(aux2);
00533     for (i=0;i<nb_dim;i++) { ab1[2*sub+1][i]= ab1[2*sub+1][i]/aux1; ab2[2*sub+1][i]= ab2[2*sub+1][i]/aux2;
00534 
00535     }  
00536     }
00537 
00538   //projection on the plan  defined by the first candidate for EACH sub-routine
00539 
00540       for (i=0;i<dim;i++)
00541         {
00542           for (sub=0;sub<4;sub++)  
00543             {
00544               projec_plan[2*sub][i]=0;projec_plan[2*sub+1][i]=0;
00545               for (l=0;l<nb_dim;l++)
00546                 {
00547                   projec_plan[2*sub][i]=Y[l][i]*ab1[2*sub+1][l]+projec_plan[2*sub][i]; 
00548                   projec_plan[2*sub+1][i]=Y[l][i]*ab1[2*sub][l]+projec_plan[2*sub+1][i]; 
00549                 }
00550             }
00551         }
00552     
00553      tmp=Evaluate_PI(projec_plan,dim);
00554      if (tmp>p_i)
00555        {p_i=tmp;p_i_old=p_i;nb_step=0; 
00556        for (l=0;l<nb_dim;l++){plan_init[0][l]=ab1[0][l];plan_init[1][l]=ab1[1][l];}
00557        }
00558      for (sub=1;sub<4;sub++)  
00559        {
00560          for (i=0;i<dim;i++){arg[0][i]= projec_plan[2*sub][i];arg[1][i]= projec_plan[2*sub+1][i];}
00561          tmp=Evaluate_PI(arg,dim);
00562          if (tmp>p_i)
00563            {p_i=tmp;p_i_old=p_i;nb_step=0;
00564            for (l=0;l<nb_dim;l++){plan_init[0][l]=ab1[2*sub][l];plan_init[1][l]=ab1[2*sub+1][l];}
00565            }
00566        }
00567   //projection on the plan defined by the second candidate for EACH sub-routine
00568 
00569       for (i=0;i<dim;i++)
00570         {
00571           for (sub=0;sub<4;sub++)  
00572             {
00573               projec_plan[2*sub][i]=0;projec_plan[2*sub+1][i]=0;
00574               for (l=0;l<nb_dim;l++)
00575                 {
00576                   projec_plan[2*sub][i]=Y[l][i]*ab2[2*sub+1][l]+projec_plan[2*sub][i]; 
00577                   projec_plan[2*sub+1][i]=Y[l][i]*ab2[2*sub][l]+projec_plan[2*sub+1][i]; 
00578                 }
00579             }
00580         }
00581      tmp=Evaluate_PI(projec_plan,dim);
00582      if (tmp>p_i)
00583        {p_i=tmp;p_i_old=p_i;nb_step=0;
00584         for (l=0;l<nb_dim;l++){plan_init[0][l]=ab2[0][l];plan_init[1][l]=ab2[1][l];}
00585        }
00586      for (sub=1;sub<4;sub++)  
00587        {
00588          for (i=0;i<dim;i++){arg[0][i]= projec_plan[2*sub][i];arg[1][i]= projec_plan[2*sub+1][i];}
00589          tmp=Evaluate_PI(arg,dim);
00590          if (tmp>p_i)
00591            {p_i=tmp;p_i_old=p_i;nb_step=0; 
00592            for (l=0;l<nb_dim;l++){plan_init[0][l]=ab2[2*sub][l];plan_init[1][l]=ab2[2*sub+1][l];}
00593            }
00594        }
00595      if ((p_i==p_i_old)&&(nb_step<=nb_max_step)){nb_step++;}
00596      if ((p_i==p_i_old)&&(nb_step>nb_max_step)) {c[0]=cbegin/factor_cons;c[1]=c[1]/factor;c[2]=c[2]/factor;c[3]=c[3]/factor;nb_step=0;}
00597 
00598     
00599 
00600 } //end while loop
00601 /*********************************************************************************************************/
00602 
00603 /*****************************saving best projection of the random departures********************************/
00604 //projection on the plan of the final solution
00605 printf("%s %e\n","Intermediate projection index=",p_i);
00606 if (p_i>p_i_opt) {
00607 p_i_opt=p_i;
00608  
00609 
00610 
00611 for (i=0;i<dim;i++)
00612 {
00613    arg[0][i]=0;arg[1][i]=0;
00614    for (l=0;l<nb_dim;l++)
00615    {
00616      arg[0][i]=Y1[l][i]*plan_init[0][l]+arg[0][i];
00617      arg[1][i]=Y1[l][i]*plan_init[1][l]+arg[1][i]; 
00618    }
00619 }   
00620 
00621 if (nb_proj==((int)(nb_images/2-1))){
00622         for (i=0;i<dim;i++)
00623         {
00624         arg[0][i]=Y1[0][i];
00625         arg[1][i]=Y1[1][i];   
00626         }   
00627 }
00628 for (i=0;i<dim;i++){
00629   for (l=0;l<2;l++){sauv[l*nb_dim+i]=(float)arg[l][i];}}
00630 ifp=NULL; 
00631 ifn=save_path[nb_proj-1][0];
00632 assert (ifp=fopen(ifn,"wb"));
00633 assert ( dim==fwrite(&sauv[0], sizeof(float),dim,ifp));
00634 assert (EOF!=fclose(ifp));
00635 
00636 for (i=0;i<dim;i++){
00637   for (l=0;l<2;l++){sauv[l*nb_dim+i]=(float)arg[1-l][i];}}
00638 
00639 ifp=NULL; 
00640 ifn=save_path[nb_proj-1][1];
00641 assert (ifp=fopen(ifn,"wb"));
00642 assert ( dim==fwrite(&sauv[0], sizeof(float),dim,ifp));
00643 assert (EOF!=fclose(ifp));
00644 }// end overwrite  
00645 skip=0;
00646 } //end for loop of new random departure
00647 printf("%s %e\n","Final projection index=",p_i_opt);
00648 
00649 /*********************************************************************************************************/
00650 
00651 /*************************Structure removal or dimension removal*******************/
00652 base_ortho_plan_finder(plan_init,new_base,nb_dim);
00653 projec_on_base(Y,new_base,dim,nb_dim);
00654 destructuration(Y,dim);
00655 projec_on_oldbase(Y,new_base,dim,nb_dim);
00656 
00657 
00658 skip=1;
00659 /**********************************************************************************************************/
00660 }//end of while p_i>=quantile loop
00661 printf("All the interesting projections have been found...\n");
00662 time(&t2);
00663 cout<< "Processing time duration : "<<t2-t1<<" seconds."<<endl;
00664 return 0;
00665 /************************************End of projection research******************************************/
00666 }
00667 
00668 /**************************************Auxiliary Fonction **************************************************/
00669 
00670 
00671 /******************************CHI2 index************************************/
00678 float Evaluate_PI(float **x, unsigned long int dim)
00679 {
00680 float res=0, p_i=0;
00681 unsigned long int *box=new unsigned long int[433];//index from 1 to 432
00682 unsigned long int i,j,k,ir,aux;
00683 float pi_help,xt1,xt2,r_i,theta_i,help,R;
00684 R=1.25;
00685 float *Ftheo= new float [6];// density in the boxes of the theorical normal distribution (normalized)
00686 Ftheo[0]=0.00387949;
00687 Ftheo[1]=0.010849;
00688 Ftheo[2]=0.0159186;
00689 Ftheo[3]=0.0184859;
00690 Ftheo[4]=0.0186418;
00691 Ftheo[5]=0.0572252;
00692 
00693 for(k=0;k<433;k++)
00694   {box[k]=0;}
00695 for(i=0;i<dim;i++)
00696   {
00697   xt1=x[0][i];xt2=x[1][i];
00698   r_i=sqrt(xt1*xt1+xt2*xt2);
00699   aux=(int)(5*r_i/R);
00700   if (aux<5) {j=72*aux+1;}
00701   else {j=72*5+1;}
00702   theta_i=atan(xt2/xt1);
00703   if (xt1<0){theta_i=theta_i+M_PI;}
00704   else if (xt2<0){theta_i=theta_i+2*M_PI;}
00705   k=j +(int)(theta_i*36/M_PI);
00706   box[k]=box[k]+1;
00707   }
00708 
00709 for (i=1;i<7;i++)
00710   {
00711     ir=(i-1)*72;pi_help=0;
00712     help=box[ir+1]+box[ir+2]+box[ir+3]+box[ir+4]+box[ir+5]+box[ir+6]+box[ir+7]+box[ir+8]+box[ir+9];
00713     pi_help=pi_help+help*help;
00714     for (j=1;j<64;j++) 
00715       {
00716        help=help-box[ir+j]+box[ir+j+9];
00717        pi_help=pi_help+ help*help;
00718       }
00719     for (j=1;j<9;j++) 
00720 
00721       {
00722        help=help-box[ir+j+63]+box[ir+j];
00723        pi_help=pi_help+ help*help;
00724       }
00725     p_i=p_i+pi_help/(Ftheo[i-1]);
00726   }
00727 p_i=p_i/(9*(float)dim*(float)dim)-1;
00728 res=p_i; 
00729 delete [] box;
00730 delete [] Ftheo;
00731 return res;
00732 }
00733 
00734 /*********************************************************************************************************/
00741 void base_ortho_plan_finder(float **plan_init,float **new_base,unsigned long int nb_dim)
00742 {
00743 unsigned long int l,m,k;
00744 float tmp;
00745 //using Gauss Jordan linear equation solution
00746 for (l=0;l<nb_dim;l++){new_base[0][l]=plan_init[0][l];new_base[1][l]=plan_init[1][l];}
00747 for (m=2;m<nb_dim;m++)
00748   {
00749   tmp=0;
00750   for (l=m;l<nb_dim;l++)
00751     {
00752     new_base[m][l]=(((float)rand())/RAND_MAX) ;//fixing degree of freedom 
00753     }
00754   }
00755 float **a= new float*[nb_dim+1];
00756 for (m=0;m<nb_dim+1;++m){a[m]= new float[nb_dim+1];} 
00757  float **b= new float*[nb_dim+1];for (m=0;m<nb_dim+1;++m){b[m]=new float[2];}
00758 for (m=2;m<nb_dim;++m)//nb_dim-2 equation systeme to be solved
00759 {
00760   //defining the systeme to be solved
00761   for (l=0;l<m;l++)
00762     {
00763     b[l+1][1]=0;//zero block
00764     for (k=0;k<nb_dim;k++)
00765       {
00766         a[l+1][k+1]=(float)new_base[l][k];//already fixed vectors
00767       }
00768     }
00769   for (l=m;l<nb_dim;l++)
00770     {
00771     b[l+1][1]=(float) new_base[m][l];//block fixed randomly
00772     for (k=0;k<nb_dim;k++)
00773       {
00774       a[l+1][k+1]=0;//zero block
00775       }
00776     a[l+1][l+1]=1;//unity block
00777     }
00778   //solving the system
00779   gaussj(a,nb_dim,b,1);
00780   for (k=0;k<nb_dim;k++){ new_base[m][k]=(float)b[k+1][1];}//filling the base with one vector more
00781   tmp=0;
00782   for (l=0;l<nb_dim;l++){tmp=new_base[m][l]*new_base[m][l]+tmp;}
00783   tmp=sqrt(tmp);
00784   for (l=0;l<nb_dim;l++){new_base[m][l]=new_base[m][l]/tmp;}//normalizing the new vector
00785 }
00786 
00787 for (m=0;m<nb_dim+1;++m){delete [] a[m];delete [] b[m];}
00788 delete [] a;delete [] b;
00789 }
00790 /*********************************************************************************************************/
00797 void projec_on_base(float **Y,float **newbase,unsigned long int dim,unsigned long int nb_dim)
00798 {
00799 unsigned long int i,m,l;
00800 float **arg= new float *[nb_dim]; 
00801 for(m=0;m<nb_dim;++m)
00802 { arg[m]=new float[dim];}
00803 for (i=0;i<dim;i++)
00804   {
00805     for (m=0;m<nb_dim;m++)
00806       {  
00807         arg[m][i]=0;
00808         for (l=0;l<nb_dim;l++){arg[m][i]=Y[l][i]*newbase[m][l]+arg[m][i];}        
00809       } 
00810   }
00811 for (i=0;i<dim;i++) {for (m=0;m<nb_dim;m++){Y[m][i]=arg[m][i];}}
00812 
00813 for (m=0;m<nb_dim;++m){delete [] arg[m];}
00814 delete [] arg;
00815 }
00816 /*********************************************************************************************************/
00824 void projec_on_oldbase(float **Y,float **newbase,unsigned long int dim,unsigned long int nb_dim)
00825 {
00826 //defining old base in new base coordinates
00827 unsigned long int i,m,l;
00828 float **arg= new float *[nb_dim]; 
00829 for(m=0;m<nb_dim;++m){ arg[m]=new float[dim];}
00830 float **a= new float*[nb_dim+1];float **b= new float*[nb_dim+1];
00831 for (m=0;m<nb_dim+1;++m){a[m]= new float[nb_dim+1];b[m]= new float[nb_dim+1];} 
00832 
00833 //matrix inverse
00834 for(m=0;m<nb_dim;m++){ for(l=0;l<nb_dim;l++){ a[m+1][l+1]=newbase[m][l];b[m+1][l+1]=0;}}
00835 for(m=0;m<nb_dim;m++){b[m+1][m+1]=1;}
00836  gaussj(a,nb_dim,b,nb_dim);//inversion
00837 //multiplication
00838 for (i=0;i<dim;i++)
00839   {
00840     for (m=0;m<nb_dim;m++)
00841       {  
00842         arg[m][i]=0;
00843         for (l=0;l<nb_dim;l++){arg[m][i]=Y[l][i]*a[m+1][l+1]+arg[m][i];}          
00844       } 
00845   }
00846 for (i=0;i<dim;i++) {for (m=0;m<nb_dim;m++){Y[m][i]=arg[m][i];}}
00847 for (m=0;m<nb_dim;++m){delete [] arg[m];}
00848 delete [] arg;
00849 for (m=0;m<nb_dim+1;++m){delete [] a[m];delete [] b[m];}
00850 delete [] a;delete [] b;
00851 }
00852 /*********************************************************************************************************/
00858 void destructuration(float **Y,unsigned long int dim)
00859 {
00860 unsigned long int i;
00861 long int *aux=new long int[1];
00862 aux[0]=-(int)(((float)rand())/RAND_MAX*10000);
00863 for (i=0;i<dim;i++)
00864   {
00865    Y[0][i]=gasdev(aux);Y[1][i]=gasdev(aux);
00866   }
00867 }
00868 /*********************************************************************************************************/
00870 double Total(float *a, unsigned long int lenght)
00871 
00872 {
00873 unsigned long int i;
00874 double res=0;
00875  for (i=0;i<lenght;i++){ res= res+(double)a[i];}
00876 return res;
00877 }
00878 
00879 /*********************************************************************************************************/

Generated on Thu Feb 17 10:58:49 2005 for Dimension Reduction by Principal Component Analysis and Projection Pursuit by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002