00001
00008
00009 #include <iostream>
00010 #include <string>
00011 #include <vector>
00012 #include <fstream>
00013
00014 #include <cstdio>
00015 #include <cstdlib>
00016 #include <cmath>
00017 #include <cassert>
00018
00019 #include "./Inliers.hh"
00020
00021 #define Inf 1.0e20
00022 #define TINY 1.0e-20
00023
00024 using namespace std;
00025
00026 Inliers::Inliers(Image<PIXEL>* data,Outliers * O,int nb)
00027 {
00028 M_LOG2=log(2);
00029 Initialise(data,O,nb);
00030 nb_inliers=(im_point->nx)*(im_point->ny);
00031 nb_init=nb;
00032 Allocate_Variable_StdDev();
00033 Initialise_IndexInliers();
00034
00035 Assigne_data_Euclidian();
00036 Update_effectif();
00037 Remove_Unaffected_RefVector(0);
00038 Update_Centers();
00039 Estimation_Stdev_RefVector();
00040 Adaptation();
00041
00042 int count=0;
00043 while (A[0].nb_ref<nb_init)
00044 {
00045 cout<<A[0].nb_ref<<" < "<<nb_init<<endl;
00046
00047 A[0].AddRefVectors(nb_init);
00048 Assigne_data_Euclidian();
00049 Update_effectif();
00050 Remove_Unaffected_RefVector(0);
00051 Update_Centers();
00052 Estimation_Stdev_RefVector();
00053 Adaptation();
00054 count++;
00055 if (count>10){break;}
00056 }
00057 MinMax();
00058 cout<<"End of Reference Vector Initialization"<<endl;
00059 }
00060
00061 void Inliers::Initialise(Image<PIXEL> * data_point,Outliers * O,int nb)
00062 {
00063 im_point=data_point;
00064 O_point=O;
00065 A=new Ref_Vector[1](im_point,nb);
00066 nbMeanAdapt=0;
00067 nbMeanIt=0;
00068
00069 maxi_dim=new float[im_point->nb_bands];
00070 mini_dim=new float[im_point->nb_bands];
00071 MinMax();
00072
00073
00074 }
00075
00076 void Inliers::Adaptation()
00077 {
00078 ChangingCriterion=Inf;
00079 int i=0;
00080 float range=0;
00081 for (int k=0;k<im_point->nb_bands;k++){range+=(maxi_dim[k]-mini_dim[k]);}
00082 range=range/im_point->nb_bands;
00083 while(ChangingCriterion>(im_point->nb_bands)*A[0].nb_ref*range/(1e3))
00084 {
00085 NearestRefVector_StdDevVar();
00086 Update_Centers();
00087 Estimation_Stdev_RefVector();
00088 i++;
00089 if (i>10){break;}
00090 }
00091 nbMeanIt+=i;
00092 nbMeanAdapt++;
00093 }
00094
00095 void Inliers::Selection()
00096 {
00097 Update_RefVectors();
00098 }
00099
00100 void Inliers::Outliers_detection()
00101 {
00102
00103 }
00104
00105 void Inliers::AddRandomRefVectors()
00106 {
00107 A[0].AddRefVectors(nb_init);
00108 ChangingCriterion=Inf;
00109 Assigne_data_Euclidian();
00110 Update_effectif();
00111 Remove_Unaffected_RefVector(0);
00112 float range=0;
00113 for (int k=0;k<im_point->nb_bands;k++){range+=(maxi_dim[k]-mini_dim[k]);}
00114 range=range/im_point->nb_bands;
00115 while(ChangingCriterion>(im_point->nb_bands)*A[0].nb_ref*range/(1e3))
00116 {
00117 Update_Centers();
00118 Estimation_Stdev_RefVector();
00119 Reassigne_inliers_StdDevVariable();
00120 Remove_Unaffected_RefVector(1);
00121 }
00122 int count=0;
00123 while (A[0].nb_ref<nb_init)
00124 {
00125 cout<<A[0].nb_ref<<" <= "<<nb_init<<endl;
00126 ChangingCriterion=Inf;
00127
00128 A[0].AddRefVectors(nb_init);
00129 Assigne_data_Euclidian();
00130 Update_effectif();
00131 Remove_Unaffected_RefVector(0);
00132 while(ChangingCriterion>(im_point->nb_bands)*A[0].nb_ref*range/(1e3))
00133 {
00134 Update_Centers();
00135 Estimation_Stdev_RefVector();
00136 Reassigne_inliers_StdDevVariable();
00137 Remove_Unaffected_RefVector(1);
00138 }
00139 count++;
00140 if (count>10){break;}
00141 }
00142 Estimation_Stdev_RefVector();
00143 MinMax();
00144 cout<<"End of Reference Vector Re-initialization"<<endl;
00145 }
00146
00147 void Inliers::Allocate_Variable_StdDev()
00148 {
00149 int k,l,j;
00150 unsigned long int dim=(im_point->nx)*(im_point->ny);
00151 if (nb_init==-1){nb_init=(int)(log((double)dim))*3.5*log((double)im_point->nb_bands)*2;}
00152 nb_init=(int)nb_init*2;
00153 nb_inliers=dim;
00154
00155 firstNearestA=new int[dim];secondNearestA=new int[dim];
00156 cross_effectif_cluster=new unsigned long int*[nb_init];
00157 effectif_cluster=new unsigned long int[nb_init];
00158 covariance_cluster=new float**[nb_init];
00159 Inv_covariance_cluster=new float**[nb_init];
00160 covariance_cluster_PostRemoval=new float***[nb_init];
00161 Inv_covariance_cluster_PostRemoval=new float***[nb_init];
00162 center_RefVector_PostRemoval=new float**[nb_init];
00163 determinant_PostRemoval=new float*[nb_init];
00164 determinant=new float[nb_init];
00165 singularity=new int[nb_init];
00166 Singularity_PostRemoval= new int*[nb_init];
00167 Index_Inliers=new unsigned long int[dim];
00168
00169 for (k=0;k<nb_init;k++)
00170 {
00171 cross_effectif_cluster[k]=new unsigned long int[nb_init];
00172 determinant_PostRemoval[k]=new float[nb_init];
00173 covariance_cluster[k]=new float*[im_point->nb_bands];
00174 Inv_covariance_cluster[k]=new float*[im_point->nb_bands];
00175 covariance_cluster_PostRemoval[k]=new float**[nb_init];
00176 Inv_covariance_cluster_PostRemoval[k]=new float**[nb_init];
00177 center_RefVector_PostRemoval[k]=new float*[nb_init];
00178 Singularity_PostRemoval[k]=new int[nb_init];
00179 }
00180 for (k=0;k<nb_init;k++)
00181 {
00182 for (l=0;l<nb_init;l++)
00183 {
00184 covariance_cluster_PostRemoval[k][l]=new float*[im_point->nb_bands];
00185 Inv_covariance_cluster_PostRemoval[k][l]=new float*[im_point->nb_bands];
00186 center_RefVector_PostRemoval[k][l]=new float[im_point->nb_bands];
00187
00188 }
00189 for(j=0;j<im_point->nb_bands;j++)
00190 {
00191 covariance_cluster[k][j]=new float[im_point->nb_bands];
00192 Inv_covariance_cluster[k][j]=new float[im_point->nb_bands];
00193 }
00194 }
00195 for (k=0;k<nb_init;k++)
00196 {
00197 for (l=0;l<nb_init;l++)
00198 {
00199 for(j=0;j<im_point->nb_bands;j++)
00200 {
00201 covariance_cluster_PostRemoval[l][k][j]=new float[im_point->nb_bands];
00202 Inv_covariance_cluster_PostRemoval[l][k][j]=new float[im_point->nb_bands];
00203 }
00204 }
00205 }
00206 nb_init=(int)((float)nb_init/2.);
00207
00208 }
00209
00210 void Inliers::MinMax()
00211 {
00212 unsigned long int i;
00213 int j;
00214 for (j=0;j<im_point->nb_bands;j++){maxi_dim[j]=0;mini_dim[j]=Inf;}
00215 unsigned long int dim=(unsigned long int)(im_point->nx)*(unsigned long int)(im_point->ny);
00216 for (i=0;i<dim;i++)
00217 {
00218 for (j=0;j<im_point->nb_bands;j++){
00219 if (maxi_dim[j]<(float)*(im_point->position+i+dim*j))
00220 {maxi_dim[j]=(float)*(im_point->position+i+dim*j);}
00221 if (mini_dim[j]>(float)*(im_point->position+i+dim*j))
00222 {mini_dim[j]=(float)*(im_point->position+i+dim*j);}
00223 }
00224 }
00225
00226 }
00227
00228 void Inliers::Initialise_IndexInliers()
00229 {
00230 unsigned long int dim=(unsigned long int)(im_point->nx)*(unsigned long int)(im_point->ny);
00231 for (unsigned long int i=0; i<dim;i++){Index_Inliers[i]=i;}
00232 }
00233
00234 void Inliers::Assigne_data_Euclidian()
00235 {
00236 int k;
00237 unsigned long int i;
00238 float max,aux;
00239
00240 for (i=0;i<nb_inliers;i++){
00241 max=-Inf;
00242 for (k=0;k<A[0].nb_ref;k++){
00243 aux=Proba_Assigne_Euclidian(i,k);
00244 if (aux>max){
00245 max=aux;
00246 firstNearestA[i]=k;
00247 }
00248 }
00249 }
00250
00251
00252 }
00253
00254
00255 float Inliers::Proba_Assigne_Euclidian(unsigned long int i, int k)
00256 {
00257 float aux=0;
00258 unsigned long int dim=(im_point->nx)*(im_point->ny);
00259 for (int j=0;j<(im_point->nb_bands);j++){
00260 aux=(float)(*(A[0][k]+j)-(float)*(im_point->position+Index_Inliers[i]+dim*j))*(float)(*(A[0][k]+j)-*(im_point->position+Index_Inliers[i]+dim*j))+aux;
00261 }
00262
00263 return 1.0/aux;
00264
00265 }
00266 float Inliers:: Proba_Assigne_Gauss_PostRemoval( unsigned long int i, int k, int l)
00267 {
00268
00269 float aux=1.;
00270 int j,j1,j2;
00271 unsigned long int dim=(im_point->nx)*(im_point->ny);
00272 float auxbis=0;
00273 float *auxligne = new float[im_point->nb_bands];
00274
00275
00276 if (Singularity_PostRemoval[k][l]!=1)
00277 {
00278 aux=sqrt(determinant_PostRemoval[k][l])/pow(2.0*M_PI,(double)(im_point->nb_bands/2.));
00279 for (j1=0;j1<(im_point->nb_bands);j1++){
00280 auxligne[j1]=0;
00281 for (j2=0;j2<(im_point->nb_bands);j2++){
00282 auxligne[j1]=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-(center_RefVector_PostRemoval[k][l][j2]))*Inv_covariance_cluster_PostRemoval[k][l][j1][j2]+auxligne[j1]; }
00283 }
00284 for (j1=0;j1<(im_point->nb_bands);j1++){
00285 auxbis=auxbis+auxligne[j1]*((float)*(im_point->position+Index_Inliers[i]+dim*j1)-(center_RefVector_PostRemoval[k][l][j1]));
00286 }
00287 aux=aux*exp(-0.5*auxbis);
00288 }
00289
00290
00291 else
00292 {
00293 aux=1;
00294 for (j=0;j<(im_point->nb_bands);j++)
00295 {
00296 if (covariance_cluster_PostRemoval[k][l][j][j]!=0){
00297
00298 aux=1.0/(sqrt(covariance_cluster_PostRemoval[k][l][j][j]*2.0*M_PI))*exp(-((center_RefVector_PostRemoval[k][l][j])-(float)*(im_point->position+Index_Inliers[i]+dim*j))*((center_RefVector_PostRemoval[k][l][j])-(float)*(im_point->position+Index_Inliers[i]+dim*j))/(2*covariance_cluster_PostRemoval[k][l][j][j]))*aux;
00299 }
00300 else
00301 {
00302 if((float)*(im_point->position+Index_Inliers[i]+dim*j)!=(center_RefVector_PostRemoval[k][l][j])){aux=0;}
00303 }
00304 }
00305
00306 }
00307
00308 delete [] auxligne;
00309 return aux;
00310 }
00311
00312
00313
00314 float Inliers::Proba_Assigne_Gauss( unsigned long int i, int k)
00315 {
00316
00317 float aux=1.;
00318 int j,j1,j2;
00319 unsigned long int dim=(im_point->nx)*(im_point->ny);
00320 float auxbis=0;
00321 float *auxligne = new float[im_point->nb_bands];
00322
00323
00324 if (singularity[k]!=1)
00325 {
00326 aux=sqrt(determinant[k])/pow(2.0*M_PI,(double)(im_point->nb_bands/2.));
00327 for (j1=0;j1<(im_point->nb_bands);j1++){
00328 auxligne[j1]=0;
00329 for (j2=0;j2<(im_point->nb_bands);j2++){
00330 auxligne[j1]=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-*(A[0][k]+j2))*Inv_covariance_cluster[k][j1][j2]+auxligne[j1]; }
00331 }
00332 for (j1=0;j1<(im_point->nb_bands);j1++){
00333 auxbis=auxbis+auxligne[j1]*((float)*(im_point->position+Index_Inliers[i]+dim*j1)-*(A[0][k]+j1));
00334 }
00335 aux=aux*exp(-0.5*auxbis);
00336 }
00337
00338
00339 else
00340 {
00341 aux=1;
00342 for (j=0;j<(im_point->nb_bands);j++)
00343 {
00344 if (covariance_cluster[k][j][j]!=0){
00345
00346 aux=1.0/(sqrt(covariance_cluster[k][j][j]*2.0*M_PI))*exp(-(*(A[0][k]+j)-(float)*(im_point->position+Index_Inliers[i]+dim*j))*(*(A[0][k]+j)-(float)*(im_point->position+Index_Inliers[i]+dim*j))/(2*covariance_cluster[k][j][j]))*aux;
00347 }
00348 else
00349 {
00350 if((float)*(im_point->position+Index_Inliers[i]+dim*j)!=*(A[0][k]+j)){aux=0;}
00351 }
00352 }
00353
00354 }
00355
00356 delete [] auxligne;
00357 return aux;
00358 }
00359
00360
00361 void Inliers::Reassigne_inliers_StdDevVariable()
00362 {
00363 int k;
00364 unsigned long int i;
00365 float max,aux;
00366
00367 for (i=0; i<nb_inliers;i++){
00368 max=-Inf;
00369 for (k=0;k<A[0].nb_ref;k++){
00370 aux=Proba_Assigne_Gauss(i,k);
00371 if (aux>max){
00372 max=aux;
00373 firstNearestA[i]=k;
00374 }
00375 }
00376 }
00377 Update_effectif();
00378 }
00379
00380 void Inliers::NearestRefVector_StdDevVar()
00381 {
00382 unsigned long int i;
00383 int k;
00384 float firstaux=0, secondaux=0;
00385 float max;
00386 float secondmax;
00387 float *aux=new float[A[0].nb_ref];
00388 for (i=0;i<nb_inliers;i++)
00389 {
00390 firstNearestA[i]=-1;
00391 secondNearestA[i]=-1;
00392 max=-Inf;
00393 secondmax=-Inf;
00394 for (k=0;k<A[0].nb_ref;k++)
00395 aux[k]=Proba_Assigne_Gauss(i,k);
00396
00397 for (k=0;k<A[0].nb_ref;k++)
00398 {
00399 firstaux=aux[k];
00400 if (firstaux> max)
00401 {
00402 firstNearestA[i]=k;
00403 max=firstaux;
00404 }
00405 }
00406 for (k=0;k<A[0].nb_ref;k++)
00407 {
00408 secondaux=aux[k];
00409 if ( (secondaux>secondmax) && firstNearestA[i]!=k )
00410 {
00411 secondNearestA[i]=k;
00412 secondmax=secondaux;
00413 }
00414 }
00415
00416 }
00417 Update_effectif();
00418 Remove_Unaffected_RefVector(1);
00419 delete [] aux;
00420 }
00421
00422
00423 void Inliers::NearestRefVector_StdDevVar_PostRemoval()
00424 {
00425 unsigned long int i;
00426 int k;
00427 float firstaux=0, secondaux=0;
00428 float max;
00429 float secondmax;
00430 for (i=0;i<nb_inliers;i++)
00431 {
00432 firstNearestA[i]=-1;
00433 secondNearestA[i]=-1;
00434 max=-Inf;
00435 secondmax=-Inf;
00436 for (k=0;k<A[0].nb_ref;k++)
00437 {
00438 firstaux=Proba_Assigne_Gauss(i,k);
00439 if (firstaux> max)
00440 {
00441 firstNearestA[i]=k;
00442 max=firstaux;
00443 }
00444 }
00445 for (k=0;k<A[0].nb_ref;k++)
00446 {
00447 if(firstNearestA[i]!=k)
00448 {
00449 secondaux=Proba_Assigne_Gauss_PostRemoval(i,k,firstNearestA[i]);
00450 if (secondaux>secondmax)
00451 {
00452 secondNearestA[i]=k;
00453 secondmax=secondaux;
00454 }
00455 }
00456 }
00457
00458
00459
00460
00461 }
00462 Estimation_CrossEffectif();
00463 Remove_Unaffected_RefVector(1);
00464 }
00465
00466
00467 void Inliers::Update_effectif()
00468 {
00469 for (int k=0;k<A[0].nb_ref;k++){effectif_cluster[k]=0;}
00470 for (unsigned long int i=0;i<nb_inliers;i++)
00471 {
00472 effectif_cluster[firstNearestA[i]]++;
00473 }
00474
00475 }
00476
00477 void Inliers::Remove_Unaffected_RefVector(int AllCov)
00478 {
00479 for (int k=0;k<A[0].nb_ref;k++)
00480 {
00481 if(effectif_cluster[k]==0)
00482 {
00483 Remove_RefVector(k,AllCov);
00484 k--;
00485 }
00486 }
00487 }
00488
00489 void Inliers::Remove_RefVector(int k, int AllCov)
00490 {
00491 unsigned long int i;
00492 int j,j1,j2,kk;
00493 A[0].nb_ref=A[0].nb_ref-1;
00494 for (i=0;i<nb_inliers;i++){
00495 if (firstNearestA[i]>k){firstNearestA[i]--;}
00496 }
00497 for (kk=k;kk<A[0].nb_ref;kk++){
00498 effectif_cluster[kk]=effectif_cluster[kk+1];
00499 for (j=0;j<im_point->nb_bands;j++){
00500 A[0][kk][j]=A[0][kk+1][j];
00501 }
00502 }
00503
00504
00505 vector<float*>::iterator it1;
00506 vector<float*>::iterator it2;
00507
00508 it1=A[0].begin()+A[0].size()-2;
00509 it2=A[0].begin()+A[0].size()-1;
00510 A[0].erase(it1,it2);
00511 nb_RefVector=A[0].nb_ref;
00512 if (AllCov==1)
00513 {
00514 for (kk=k;kk<A[0].nb_ref;kk++){
00515 singularity[kk]=singularity[kk+1];
00516 determinant[kk]=determinant[kk+1];
00517 for (j1=0;j1<im_point->nb_bands;j1++){
00518 for (j2=0;j2<im_point->nb_bands;j2++){
00519 covariance_cluster[kk][j1][j2]=covariance_cluster[kk+1][j1][j2];
00520 Inv_covariance_cluster[kk][j1][j2]=Inv_covariance_cluster[kk+1][j1][j2];
00521 }
00522 }
00523 }
00524 for (i=0;i<nb_inliers;i++){
00525 if (secondNearestA[i]>k){secondNearestA[i]--;}
00526 }
00527
00528 }
00529 }
00530
00531 void Inliers::Update_Centers()
00532 {
00533 unsigned long int i;
00534 unsigned long int dim=(im_point->nx)*(im_point->ny);
00535 int j,k;
00536 float auxi;
00537 float **aux=new float*[A[0].nb_ref];
00538 for (k=0;k<A[0].nb_ref;k++){aux[k]=new float[im_point->nb_bands];}
00539 for (k=0;k<A[0].nb_ref;k++){for (j=0;j<im_point->nb_bands;j++){aux[k][j]=0;}}
00540
00541 ChangingCriterion=0;
00542 for (i=0;i<nb_inliers;i++)
00543 {
00544 for (j=0;j<im_point->nb_bands;j++)
00545 {
00546 aux[firstNearestA[i]][j]= (float)*(im_point->position+Index_Inliers[i]+dim*j)+aux[firstNearestA[i]][j];
00547 }
00548 }
00549 for (k=0;k<A[0].nb_ref;k++)
00550 {
00551 auxi=0;
00552 for (j=0;j<im_point->nb_bands;j++)
00553 {
00554 auxi=(A[0][k][j]-aux[k][j]/(float)effectif_cluster[k])*(A[0][k][j]-aux[k][j]/(float)effectif_cluster[k])+auxi;
00555 }
00556 ChangingCriterion+=auxi;
00557 }
00558 Change=1;
00559 ChangingCriterion=sqrt(ChangingCriterion);
00560 float range=0;
00561 for (int k=0;k<im_point->nb_bands;k++){range+=(maxi_dim[k]-mini_dim[k]);}
00562 range=range/im_point->nb_bands;
00563 if (ChangingCriterion<(im_point->nb_bands)*A[0].nb_ref*range/(1e4)){Change=0;}
00564
00565 for (k=0;k<A[0].nb_ref;k++){
00566 for (j=0;j<im_point->nb_bands;j++)
00567 {
00568 A[0][k][j]=aux[k][j]/(float)effectif_cluster[k];
00569 }
00570 }
00571
00572 delete [] aux;
00573 }
00574
00575 void Inliers::Estimation_CovarianceMatrix()
00576 {
00577 unsigned long int i;
00578 unsigned long int dim=(unsigned long int)(im_point->nx*im_point->ny);
00579 int l,j,k;
00580
00581 for(k=0;k<A[0].nb_ref;k++){
00582 for(l=0;l<im_point->nb_bands;l++){
00583 for(j=0;j<im_point->nb_bands;j++){ covariance_cluster[k][l][j]=0;}}}
00584
00585
00586 for (i=0;i<nb_inliers;i++)
00587 {
00588 for (l=0;l<im_point->nb_bands;l++)
00589 {
00590 for (j=l;j<im_point->nb_bands;j++)
00591 {
00592
00593 covariance_cluster[firstNearestA[i]][l][j]=((float)*(im_point->position+Index_Inliers[i]+dim*j)-A[0][firstNearestA[i]][j])*((float)*(im_point->position+Index_Inliers[i]+dim*l)-A[0][firstNearestA[i]][l])+covariance_cluster[firstNearestA[i]][l][j];
00594 }
00595 }
00596 }
00597 for(k=0;k<A[0].nb_ref;k++)
00598 {
00599 if (effectif_cluster[k]!=0)
00600 {
00601 for (l=0;l<im_point->nb_bands;l++)
00602 {
00603 for(j=l;j<im_point->nb_bands;j++)
00604 {
00605 covariance_cluster[k][l][j]=covariance_cluster[k][l][j]/effectif_cluster[k];
00606 covariance_cluster[k][j][l]=covariance_cluster[k][l][j];
00607 }
00608 }
00609 }
00610 }
00611 }
00612
00613 void Inliers::Singularity_Detection()
00614 {
00615 int j;
00616 for (int k=0;k<A[0].nb_ref;k++)
00617 {
00618 singularity[k]=0;
00619 for (j=0;j<im_point->nb_bands;j++)
00620 {
00621 if (covariance_cluster[k][j][j]==0)
00622 {
00623 singularity[k]=1;
00624 }
00625 }
00626 }
00627
00628 }
00629
00630 void Inliers::Estimation_Stdev_RefVector()
00631 {
00632 Estimation_CovarianceMatrix();
00633 Singularity_Detection();
00634 Invert_CovarianceMatrix();
00635
00636 }
00637
00638 void Inliers::Invert_CovarianceMatrix()
00639 {
00640 int l,j,k;
00641 float **auxCOV=new float*[im_point->nb_bands+1];
00642 for (k=0;k<im_point->nb_bands+1;k++){auxCOV[k]=new float[im_point->nb_bands+1];}
00643 float *col=new float[im_point->nb_bands+1];
00644 auxCOV[0][0]=0;
00645 float d;
00646 int *indx=new int[im_point->nb_bands+1], n=(int)im_point->nb_bands;
00647 void ludcmp(float **a,int n,int *indx,float *d);
00648 void lubksb(float **a,int n,int *indx, float b[]);
00649
00650 for (k=0;k<A[0].nb_ref;k++)
00651 {
00652 for (l=0;l<im_point->nb_bands;l++)
00653 {
00654 for(j=l;j<im_point->nb_bands;j++)
00655 {
00656 auxCOV[l+1][j+1]=covariance_cluster[k][l][j];
00657 auxCOV[j+1][l+1]=auxCOV[l+1][j+1];
00658 }
00659
00660 }
00661 if (singularity[k]!=1)
00662 {
00663
00664 ludcmp(auxCOV,n,indx,&d);
00665 for (j=1;j<im_point->nb_bands+1;j++){
00666 for (l=1;l<im_point->nb_bands+1;l++){col[l]=0.0;}
00667 col[j]=1.0;
00668 lubksb(auxCOV,n,indx,col);
00669 for (l=1;l<im_point->nb_bands+1;l++)
00670 {
00671 Inv_covariance_cluster[k][l-1][j-1]=col[l];
00672 }
00673 }
00674 for (j=1;j<im_point->nb_bands+1;j++)
00675 {
00676 for (l=1;l<im_point->nb_bands+1;l++)
00677 {
00678 auxCOV[j][l]=Inv_covariance_cluster[k][j-1][l-1];
00679 }
00680 }
00681
00682
00683 ludcmp(auxCOV,n,indx,&d);
00684 for (j=1;j<im_point->nb_bands+1;j++){d=auxCOV[j][j]*d;}
00685 determinant[k]=d;
00686
00687 }
00688 else {determinant[k]=0;}
00689
00690 }
00691 for(k=0;k<im_point->nb_bands+1;k++){delete [] auxCOV[k];}
00692 delete [] auxCOV;delete [] col;
00693
00694
00695
00696 }
00697
00698 void Inliers::Estimation_CrossEffectif()
00699 {
00700 int k;
00701 for(int l=0;l<A[0].nb_ref;l++)
00702 {
00703 effectif_cluster[l]=0;
00704 for (k=0;k<A[0].nb_ref;k++)
00705 {
00706 cross_effectif_cluster[l][k]=0;
00707 }
00708 }
00709 for(unsigned long int i=0;i<nb_inliers;i++)
00710 {
00711 effectif_cluster[firstNearestA[i]]++;
00712 cross_effectif_cluster[firstNearestA[i]][secondNearestA[i]]++;
00713
00714 }
00715
00716 }
00717
00718 void Inliers::Estimation_Center_RefVector_PostRemoval()
00719 {
00720 unsigned long int i;
00721 unsigned long int dim=(im_point->nx)*(im_point->ny);
00722 int l,k,j;
00723 Estimation_CrossEffectif();
00724 for(l=0;l<A[0].nb_ref;l++)
00725 {
00726 for (k=0;k<A[0].nb_ref;k++)
00727 {
00728 for(j=0;j<im_point->nb_bands;j++)
00729 {
00730
00731 center_RefVector_PostRemoval[l][k][j]=0;
00732 }
00733 }
00734 }
00735 for(i=0;i<nb_inliers;i++)
00736 {
00737 for (j=0;j<im_point->nb_bands;j++)
00738 {
00739
00740 center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j]=center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j]+((float)*(im_point->position+Index_Inliers[i]+dim*j));
00741 }
00742 }
00743
00744 for(l=0;l<A[0].nb_ref;l++)
00745 {
00746 for (k=0;k<A[0].nb_ref;k++)
00747 {
00748 for(j=0;j<im_point->nb_bands;j++)
00749 { center_RefVector_PostRemoval[l][k][j]=(center_RefVector_PostRemoval[l][k][j]+A[0][l][j]*(float)effectif_cluster[l])/((float)effectif_cluster[l]+(float)cross_effectif_cluster[k][l]);
00750
00751 }
00752 }
00753 }
00754
00755
00756 }
00757
00758 void Inliers::Estimation_CovarianceMatrix_PostRemoval()
00759 {
00760 unsigned long int i;
00761 unsigned long int dim=(im_point->nx)*(im_point->ny);
00762 float aux1,aux2;
00763 int l,k,j1,j2;
00764 for(l=0;l<A[0].nb_ref;l++)
00765 {
00766 for (k=0;k<A[0].nb_ref;k++)
00767 {
00768 for(j1=0;j1<im_point->nb_bands;j1++)
00769 {
00770 for(j2=0;j2<im_point->nb_bands;j2++)
00771 {
00772 covariance_cluster_PostRemoval[l][k][j1][j2]=0;
00773 }
00774 }
00775 }
00776 }
00777
00778
00779 for(i=0;i<nb_inliers;i++)
00780 {
00781 for (j1=0;j1<im_point->nb_bands;j1++)
00782 {
00783 for (j2=j1;j2<im_point->nb_bands;j2++)
00784 { aux1=((float)*(im_point->position+Index_Inliers[i]+dim*j1)-center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j1]);
00785 aux2=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j2]);
00786
00787 covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j1][j2]=covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j1][j2]+aux1*aux2;
00788 }
00789 }
00790 }
00791
00792
00793 for(l=0;l<A[0].nb_ref;l++)
00794 {
00795 for (k=0;k<A[0].nb_ref;k++)
00796 {
00797 for(j1=0;j1<im_point->nb_bands;j1++)
00798 {
00799 for(j2=j1;j2<im_point->nb_bands;j2++)
00800 {
00801
00802 covariance_cluster_PostRemoval[l][k][j1][j2]=(covariance_cluster_PostRemoval[l][k][j1][j2]+covariance_cluster[l][j1][j2]*(float)effectif_cluster[l])/((float)effectif_cluster[l]+(float)cross_effectif_cluster[k][l]);
00803 covariance_cluster_PostRemoval[l][k][j2][j1]=covariance_cluster_PostRemoval[l][k][j1][j2];
00804 }
00805
00806 }
00807 }
00808 }
00809 }
00810
00811 void Inliers::Singularity_Detection_Postremoval()
00812 {
00813 int j1,j2,j,nb1=0,nb2=0;;
00814 for (j1=0;j1<A[0].nb_ref;j1++){ singularity[j1]=0;
00815 for (j2=0;j2<A[0].nb_ref;j2++){Singularity_PostRemoval[j1][j2]=0;
00816 }
00817 }
00818
00819 for (j=0;j<im_point->nb_bands;j++){
00820 for (j1=0;j1<A[0].nb_ref;j1++){
00821 if (covariance_cluster[j1][j][j]==0){singularity[j1]=1;nb1++;}
00822 for (j2=0;j2<A[0].nb_ref;j2++){
00823 if (covariance_cluster_PostRemoval[j1][j2][j][j]==0){Singularity_PostRemoval[j1][j2]=1;nb2++;}
00824 }
00825 }
00826 }
00827
00828
00829 }
00830
00831 void Inliers::Invert_PostCovarianceMatrix()
00832 {
00833 int l,k,k1,k2,j;
00834 float **auxCOV=new float*[im_point->nb_bands+1];
00835 for (k=0;k<im_point->nb_bands+1;k++){auxCOV[k]=new float[im_point->nb_bands+1];}
00836 float *col=new float[im_point->nb_bands+1];
00837 auxCOV[0][0]=0;
00838 float d;
00839 int *indx=new int[im_point->nb_bands+1], n=(int)im_point->nb_bands;
00840
00841 void ludcmp(float **a,int n,int *indx,float *d);
00842 void lubksb(float **a,int n,int *indx, float b[]);
00843
00844
00845
00846 for (k1=0;k1<A[0].nb_ref;k1++)
00847 {
00848 for (k2=0;k2<A[0].nb_ref;k2++)
00849 {
00850 Singularity_PostRemoval[k1][k2]=0;
00851 for (l=0;l<im_point->nb_bands;l++)
00852 {
00853 if (covariance_cluster_PostRemoval[k1][k2][l][l]==0)
00854 {
00855 Singularity_PostRemoval[k1][k2]=1;
00856 }
00857 for(j=l;j<im_point->nb_bands;j++)
00858 {
00859 auxCOV[l+1][j+1]=covariance_cluster_PostRemoval[k1][k2][l][j];
00860 auxCOV[j+1][l+1]=auxCOV[l+1][j+1];
00861 }
00862 }
00863
00864 if (Singularity_PostRemoval[k1][k2]!=1)
00865 {
00866
00867 ludcmp(auxCOV,n,indx,&d);
00868 for (j=1;j<im_point->nb_bands+1;j++){
00869 for (l=1;l<im_point->nb_bands+1;l++){col[l]=0.0;}
00870 col[j]=1.0;
00871 lubksb(auxCOV,n,indx,col);
00872 for (l=1;l<im_point->nb_bands+1;l++)
00873 {
00874 Inv_covariance_cluster_PostRemoval[k1][k2][l-1][j-1]=col[l];
00875 }
00876 }
00877 for (j=1;j<im_point->nb_bands+1;j++){
00878 for (l=1;l<im_point->nb_bands+1;l++){
00879 auxCOV[l][j]=Inv_covariance_cluster_PostRemoval[k1][k2][l-1][j-1];
00880 }
00881 }
00882
00883 ludcmp(auxCOV,n,indx,&d);
00884 for (j=1;j<im_point->nb_bands+1;j++){d*=auxCOV[j][j];}
00885 determinant_PostRemoval[k1][k2]=d;
00886 }
00887 else {determinant_PostRemoval[k1][k2]=0;}
00888
00889 }
00890 }
00891 for(k=0;k<im_point->nb_bands+1;k++){delete [] auxCOV[k];}
00892 delete [] auxCOV;delete [] col;
00893 }
00894
00895 void Inliers::Estimation_Stdev_RefVector_PostRemoval()
00896 {
00897
00898 Estimation_Center_RefVector_PostRemoval();
00899 Estimation_CovarianceMatrix_PostRemoval();
00900 Singularity_Detection_Postremoval();
00901 Invert_PostCovarianceMatrix();
00902 }
00903
00904 void Inliers::Inliers_Variation_Index_term_cost2(float* aux)
00905 {
00906 int l,k;
00907 Estimation_CrossEffectif();
00908 for (k=0;k<A[0].nb_ref;k++)
00909 {
00910 aux[k]=(float)effectif_cluster[k]*log((float)effectif_cluster[k]/(float)nb_inliers)/M_LOG2;
00911 for (l=0;l<A[0].nb_ref;l++)
00912 {
00913 if (l!=k)
00914 {
00915 aux[k]=aux[k]-(float)((effectif_cluster[l]+cross_effectif_cluster[k][l])*log((float)(effectif_cluster[l]+cross_effectif_cluster[k][l])/(float)nb_inliers)/M_LOG2-effectif_cluster[l]*log((float)effectif_cluster[l]/(float)nb_inliers)/M_LOG2);
00916 }
00917 }
00918 }
00919 }
00920
00921 float Inliers::Error_FreeInliers(unsigned long int i)
00922 {
00923 unsigned long int dim=(unsigned long int)(im_point->nx*im_point->ny);
00924 int j1,j2;
00925 float aux1,aux2;
00926 float *auxligne = new float[im_point->nb_bands];
00927 float ratio=0;
00928
00929 aux2=0;
00930 for (j1=0;j1<(im_point->nb_bands);j1++)
00931 {
00932 auxligne[j1]=0;
00933 for (j2=0;j2<(im_point->nb_bands);j2++)
00934 { auxligne[j1]=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-(A[0][firstNearestA[i]][j2]))*Inv_covariance_cluster[firstNearestA[i]][j1][j2]+auxligne[j1];
00935 }
00936 }
00937
00938 for (j1=0;j1<(im_point->nb_bands);j1++){
00939 aux2=aux2+auxligne[j1]*((float)*(im_point->position+Index_Inliers[i]+dim*j1)-(A[0][firstNearestA[i]][j1]));
00940 }
00941
00942 aux1=0;
00943 for (j1=0;j1<(im_point->nb_bands);j1++){
00944 auxligne[j1]=0;
00945 for (j2=0;j2<(im_point->nb_bands);j2++)
00946 { auxligne[j1]=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j2])*Inv_covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j1][j2]+auxligne[j1];
00947 }
00948 }
00949
00950 for (j1=0;j1<(im_point->nb_bands);j1++){
00951 aux1=aux1+auxligne[j1]*((float)*(im_point->position+Index_Inliers[i]+dim*j1)-center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j1]);
00952 }
00953 ratio=(aux1-aux2)/(2.0*M_LOG2)+0.5/M_LOG2*log(determinant[firstNearestA[i]]/determinant_PostRemoval[secondNearestA[i]][firstNearestA[i]]);
00954
00955
00956
00957 delete [] auxligne;
00958 return ratio;
00959 }
00960
00961 float Inliers::Error_FreeInliers_Singularity(unsigned long int i)
00962 {
00963
00964
00965 unsigned long int dim=(unsigned long int)(im_point->nx*im_point->ny);
00966 int j;
00967 float aux1,aux2,ratio=0;
00968 if (singularity[firstNearestA[i]]==1 || Singularity_PostRemoval[secondNearestA[i]][firstNearestA[i]]==1)
00969 {
00970 for (j=0;j<im_point->nb_bands;j++)
00971 {
00972
00973 if (covariance_cluster[firstNearestA[i]][j][j]!=0&&covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j]!=0)
00974 { aux1=((float)*(im_point->position+Index_Inliers[i]+dim*j)-center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j]);
00975 aux2=((float)*(im_point->position+Index_Inliers[i]+dim*j)-A[0][firstNearestA[i]][j]);
00976 ratio=aux1*aux1/(2.*M_LOG2*covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j])-aux2*aux2/(2.*M_LOG2*covariance_cluster[firstNearestA[i]][j][j])+0.5*log(covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j]/covariance_cluster[firstNearestA[i]][j][j])/M_LOG2+ratio;
00977 }
00978
00979 else
00980 {
00981 if (covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j]!=0)
00982 { aux1=((float)*(im_point->position+Index_Inliers[i]+dim*j)-center_RefVector_PostRemoval[secondNearestA[i]][firstNearestA[i]][j]);
00983 aux1=aux1*aux1/(2.*M_LOG2*covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j])+0.5*log(covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j])/M_LOG2;
00984 ratio=aux1;
00985 }
00986 if (covariance_cluster[firstNearestA[i]][j][j]!=0)
00987 { aux1=((float)*(im_point->position+Index_Inliers[i]+dim*j)-A[0][firstNearestA[i]][j]);
00988 aux1=aux1*aux1/(2.*M_LOG2*covariance_cluster[firstNearestA[i]][j][j])+0.5*log(covariance_cluster[firstNearestA[i]][j][j])/M_LOG2;
00989 ratio=-aux1;
00990 }
00991 }
00992
00993 }
00994 }
00995 return ratio;
00996 }
00997
00998 void Inliers::Error_OtherInliers(float * SumErrorRatio)
00999 {
01000 float ratio,aux1,aux2;
01001 float *auxligne = new float[im_point->nb_bands];
01002 unsigned long int dim=(unsigned long int)(im_point->nx*im_point->ny);
01003 unsigned long int i;
01004 int k,j1,j2;
01005 for (k=0;k<A[0].nb_ref;k++)
01006 {
01007 for (i=0;i<nb_inliers;i++)
01008 {
01009 if (singularity[firstNearestA[i]]==0&&Singularity_PostRemoval[firstNearestA[i]][k]==0)
01010 {
01011 if (firstNearestA[i]!=k)
01012 {
01013 aux2=0;
01014 for (j1=0;j1<(im_point->nb_bands);j1++){
01015 auxligne[j1]=0;
01016 for (j2=0;j2<(im_point->nb_bands);j2++)
01017 { auxligne[j1]=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-*(A[0][firstNearestA[i]]+j2))*Inv_covariance_cluster[firstNearestA[i]][j1][j2]+auxligne[j1];
01018 }
01019 }
01020 for (j1=0;j1<(im_point->nb_bands);j1++){
01021 aux2=aux2+auxligne[j1]*((float)*(im_point->position+Index_Inliers[i]+dim*j1)-*(A[0][firstNearestA[i]]+j1));
01022 }
01023
01024 aux1=0;
01025 for (j1=0;j1<(im_point->nb_bands);j1++){
01026 auxligne[j1]=0;
01027 for (j2=0;j2<(im_point->nb_bands);j2++)
01028 { auxligne[j1]=((float)*(im_point->position+Index_Inliers[i]+dim*j2)-center_RefVector_PostRemoval[firstNearestA[i]][k][j2])*Inv_covariance_cluster_PostRemoval[firstNearestA[i]][k][j1][j2]+auxligne[j1];
01029 }
01030 }
01031
01032 for (j1=0;j1<(im_point->nb_bands);j1++){
01033 aux1=aux1+auxligne[j1]*((float)*(im_point->position+Index_Inliers[i]+dim*j1)-center_RefVector_PostRemoval[firstNearestA[i]][k][j1]);
01034 }
01035 ratio=(aux1-aux2)/(2.0*M_LOG2)+0.5/M_LOG2*log(determinant[firstNearestA[i]]/determinant_PostRemoval[firstNearestA[i]][k]);
01036 SumErrorRatio[k]=SumErrorRatio[k]+ratio;
01037 }
01038 }
01039 }
01040 }
01041 delete [] auxligne;
01042 }
01043
01044 void Inliers::Inliers_Variation_Error_term_cost4(float * SumErrorRatio)
01045 {
01046 unsigned long int i;
01047 int k;
01048 for (k=0;k<A[0].nb_ref;k++) {SumErrorRatio[k]=0;}
01049 float ratio=0;
01050 Estimation_Stdev_RefVector();
01051 Estimation_Stdev_RefVector_PostRemoval();
01052 NearestRefVector_StdDevVar_PostRemoval();
01053
01054
01055 for (i=0;i<nb_inliers;i++)
01056 {
01057 ratio=0;
01058
01059
01060 if ((singularity[firstNearestA[i]]==0) && (Singularity_PostRemoval[secondNearestA[i]][firstNearestA[i]]==0))
01061 {
01062 ratio=Error_FreeInliers(i);
01063 }
01064
01065 else
01066 {
01067 ratio=Error_FreeInliers_Singularity(i);
01068
01069 }
01070
01071 SumErrorRatio[firstNearestA[i]]=SumErrorRatio[firstNearestA[i]]+ratio;
01072
01073 }
01074
01075
01076 Error_OtherInliers(SumErrorRatio);
01077 }
01078
01079 void Inliers::Update_RefVectors()
01080 {
01081
01082 unsigned long int i;
01083 int k;
01084 unsigned long int dim=(unsigned long int)(im_point->nx*im_point->ny);
01085 float ChangeCodeLength;
01086 float* SumErrorRatio=new float[dim];
01087 float* Indexcost2=new float[dim];
01088 for (i=0;i<dim;i++){SumErrorRatio[i]=0;}
01089 int nb_ref_previous=0;
01090 float FactorOverCostModel3=0;
01091 float maxidecrease;
01092 int indexmaxidecrease;
01093 cout<<"number of reference Vector after selection step : "<<A[0].nb_ref<<endl;
01094 while ((nb_ref_previous!=A[0].nb_ref) && (A[0].nb_ref>1))
01095 {
01096 nb_ref_previous=A[0].nb_ref;
01097 Inliers_Variation_Error_term_cost4(SumErrorRatio);
01098 FactorOverCostModel3=(im_point->nb_bands)+1;
01099 Inliers_Variation_Index_term_cost2(Indexcost2);
01100
01101
01102 maxidecrease=Inf;
01103 indexmaxidecrease=0;
01104 for (k=0;k<A[0].nb_ref;k++)
01105 {
01106 ChangeCodeLength=(-1.)*(im_point->nb_bands)*sizeof(PIXEL)*8*FactorOverCostModel3+Indexcost2[k]+SumErrorRatio[k];
01107 if (maxidecrease>ChangeCodeLength){maxidecrease=ChangeCodeLength;indexmaxidecrease=k;}
01108
01109 }
01110 k=indexmaxidecrease;
01111
01112
01113 ChangeCodeLength=(-1.)*(im_point->nb_bands)*sizeof(PIXEL)*8*FactorOverCostModel3+Indexcost2[k]+SumErrorRatio[k];
01114
01115 if (ChangeCodeLength<0)
01116 {
01117 Remove_RefVector(k,1);
01118 Adaptation();
01119 }
01120 cout<<"number of reference Vector intermediate : "<<A[0].nb_ref<<endl;
01121
01122 }
01123 nb_RefVector=A[0].nb_ref;
01124
01125 delete [] SumErrorRatio;
01126 delete [] Indexcost2;
01127
01128 }
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01158
01159
01160 void Inliers::Print_Mean_Covariance(int l)
01161 {
01162 int j2,j1;
01163 if (l<0)
01164 {
01165 for(l=0;l<A[0].nb_ref;l++)
01166 { cout<<"\nReference Vector "<<l<<", determinant: "<<determinant[l]<<", effectif: "<<effectif_cluster[l]<<", Singularity : "<<singularity[l]<<"\n";
01167 cout<<"Mean : Covariance : "<<endl;
01168 for(j2=0;j2<im_point->nb_bands;j2++)
01169 {
01170
01171 cout<<A[0][l][j2]<<" ";
01172 for (j1=0;j1<im_point->nb_bands;j1++)
01173 {
01174 cout<<" "<<(covariance_cluster[l][j2][j1]);
01175 }
01176 cout<<"\n";
01177 }
01178 }
01179 cout<<""<<endl;
01180 }
01181 else
01182 {
01183 cout<<"\nReference Vector "<<l<<", determinant: "<<determinant[l]<<", effectif: "<<effectif_cluster[l]<<", Singularity : "<<singularity[l]<<"\n";
01184 cout<<"Mean : Covariance : "<<endl;
01185 for(j2=0;j2<im_point->nb_bands;j2++)
01186 {
01187
01188 cout<<A[0][l][j2]<<" ";
01189 for (j1=0;j1<im_point->nb_bands;j1++)
01190 {
01191 cout<<" "<<(covariance_cluster[l][j2][j1]);
01192 }
01193 cout<<"\n";
01194 }
01195 }
01196 }
01197
01198 void Inliers::Print_Mean_Covariance_PostRemoval(int l)
01199 {
01200 int j2,j1,k;
01201 if (l<0)
01202 {
01203
01204
01205 for(l=0;l<A[0].nb_ref;l++)
01206 { cout<<"\nReference Vector "<<l<<", determinant: "<<determinant[l]<<", effectif: "<<effectif_cluster[l]<<", Singularity : "<<singularity[l]<<"\n";
01207 cout<<"Mean : Inv_covariance_cluster : "<<endl;
01208 for(j2=0;j2<im_point->nb_bands;j2++)
01209 {
01210
01211 cout<<A[0][l][j2]<<" ";
01212 for (j1=0;j1<im_point->nb_bands;j1++)
01213 {
01214 cout<<" "<<(covariance_cluster[l][j2][j1]);
01215 }
01216 cout<<"\n";
01217 }
01218 for(k=0;k<A[0].nb_ref;k++)
01219 {
01220 cout<<"\n\t If Removal Vector "<<k<<", determinant: "<<determinant_PostRemoval[l][k]<<", effectif: "<<cross_effectif_cluster[k][l]+effectif_cluster[l]<<", Singularity : "<<Singularity_PostRemoval[l][k]<<"\n";
01221 cout<<"\tMean : Inv_covariance_cluster : "<<endl;
01222 for(j2=0;j2<im_point->nb_bands;j2++)
01223 {
01224
01225 cout<<"\t"<<center_RefVector_PostRemoval[l][k][j2]<<" ";
01226 for (j1=0;j1<im_point->nb_bands;j1++)
01227 {
01228 cout<<" "<<(covariance_cluster_PostRemoval[l][k][j2][j1]);
01229 }
01230 cout<<"\n";
01231 }
01232 }
01233
01234 }
01235 cout<<""<<endl;
01236 }
01237 else
01238 {
01239 cout<<"\nReference Vector "<<l<<", determinant: "<<determinant[l]<<", effectif: "<<effectif_cluster[l]<<", Singularity : "<<singularity[l]<<"\n";
01240 cout<<"Mean : Inv_covariance_cluster : "<<endl;
01241 for(j2=0;j2<im_point->nb_bands;j2++)
01242 {
01243
01244 cout<<A[0][l][j2]<<" ";
01245 for (j1=0;j1<im_point->nb_bands;j1++)
01246 {
01247 cout<<" "<<(covariance_cluster[l][j2][j1]);
01248 }
01249 cout<<"\n";
01250 }
01251 for(k=0;k<A[0].nb_ref;k++)
01252 {
01253 cout<<"\n\t If Removal Vector "<<k<<", determinant: "<<determinant_PostRemoval[l][k]<<", effectif: "<<cross_effectif_cluster[k][l]+effectif_cluster[l]<<", Singularity : "<<Singularity_PostRemoval[l][k]<<"\n";
01254 cout<<"\tMean : Inv_covariance_cluster : "<<endl;
01255 for(j2=0;j2<im_point->nb_bands;j2++)
01256 {
01257
01258 cout<<"\t"<<center_RefVector_PostRemoval[l][k][j2]<<" ";
01259 for (j1=0;j1<im_point->nb_bands;j1++)
01260 {
01261 cout<<" "<<(covariance_cluster_PostRemoval[l][k][j2][j1]);
01262 }
01263 cout<<"\n";
01264 }
01265 }
01266 }
01267 }