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
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
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
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
00143
00144 float eigenthreshold=0.000001;
00145 int nb_rand_departure=5;
00146 float cbegin=tan(80*M_PI/180);
00147 float eps=1e-3;
00148 float factor=2;
00149 float factor_cons=1;
00150 int nb_step, nb_max_step=50;
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
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
00266 cout<<"\n";
00267 printf("Reading images\n");
00268 PIX a[dimx][dimy] ;
00269 float **clust= new float *[nb_dim];
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);
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
00298 cout<<"\nDecorralating procedure"<<endl;
00299
00300 float **Y= new float *[nb_dim];
00301 for(m=0;m<nb_dim;++m)
00302 { Y[m]=new float[dim];}
00303 float **Y1= new float *[nb_dim];
00304 for(m=0;m<nb_dim;++m)
00305 { Y1[m]=new float[dim];}
00306 float *sauvACP= new float [dim];
00307 float *moy= new float [nb_dim];
00308 float *M= new float [nb_dim];
00309 float *eigenvalues= new float [nb_dim+1];
00310 float **cov= new float *[nb_dim+1];
00311 float **eigenvectors= new float *[nb_dim+1];
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
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
00332 for (l=0;l<nb_dim;l++){
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
00341 }
00342 }
00343
00344 jacobi(cov,nb_dim,eigenvalues,eigenvectors,nrot);
00345 eigsrt(eigenvalues,eigenvectors,nb_dim);
00346
00347
00348
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;
00366 cout<<"Number of component representing "<<keepEnergie<<" of the signal energie : "<<nb_dim<<endl;
00367
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
00393 delete [] eigenvectors;
00394
00395
00396 delete [] eigenvalues;
00397 delete [] sauvACP;
00398
00399
00400
00401
00402 cout<<"\nInitialisation of optimisation procedure"<<endl;
00403 int nb_proj=0;
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];
00408 float **new_base = new float*[nb_dim];
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];
00422
00423 int sub,skip=0;
00424 p_i=1000;
00425 p_i_opt=p_i;
00426
00427
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++)
00435 {
00436
00437 if (skip==0){
00438
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;}
00445
00446 tmp=0;
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]);
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;}
00456 }
00457
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;
00472
00473
00474
00475
00476
00477 for(i=0;i<4;i++){c[i]=cbegin;}
00478 while (c[1]>eps)
00479 {
00480
00481
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;}
00489 }
00490
00491
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
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
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 }
00601
00602
00603
00604
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 }
00645 skip=0;
00646 }
00647 printf("%s %e\n","Final projection index=",p_i_opt);
00648
00649
00650
00651
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 }
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
00666 }
00667
00668
00669
00670
00671
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];
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];
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
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) ;
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)
00759 {
00760
00761 for (l=0;l<m;l++)
00762 {
00763 b[l+1][1]=0;
00764 for (k=0;k<nb_dim;k++)
00765 {
00766 a[l+1][k+1]=(float)new_base[l][k];
00767 }
00768 }
00769 for (l=m;l<nb_dim;l++)
00770 {
00771 b[l+1][1]=(float) new_base[m][l];
00772 for (k=0;k<nb_dim;k++)
00773 {
00774 a[l+1][k+1]=0;
00775 }
00776 a[l+1][l+1]=1;
00777 }
00778
00779 gaussj(a,nb_dim,b,1);
00780 for (k=0;k<nb_dim;k++){ new_base[m][k]=(float)b[k+1][1];}
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;}
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
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
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);
00837
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