Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

Inliers.cpp

Go to the documentation of this file.
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();//memory allocation
00033         Initialise_IndexInliers();//usefull when the outliers are considered, SHOULD BE REMOVED!!
00034 
00035         Assigne_data_Euclidian();//assignement of the data points according to their euclidian distance to the Gaussian Means.
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)//while there are not enough components in the mixture...iterate the procedure
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;// varaible to check how good the algorithm has converged 
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)) //precision inversely proportional to the number of components
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         //Update_outliers();
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     // no singularity in the covariance matrix
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     //singularity : independence intra bande simplification
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     // no singularity in the covariance matrix
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     //singularity : independence intra bande simplification
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                         //max=Proba_Assigne_Gauss(i,firstNearestA[i])-Proba_Assigne_Gauss_PostRemoval(i,secondNearestA[i],firstNearestA[i]);
00458                         //if (max<0){
00459                         //cout<<"ecart "<<max<<endl;
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;}//convergence finale!!!
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         //Covariance Matrix
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                         //INVERSE of COV MATRIX
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                         //DETERMINANT OF THE INVERSE COV MATRIX
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         //COVARIANCE POSTREMOVAL ESTIMATION : 1- New inliers contribution
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         //COVARIANCE POSTREMOVAL ESTIMATION : 2- Old inliers contribution
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                                 //if (((float)effectif_cluster[l]+(float)cross_effectif_cluster[k][l])==0){cout<<"putain\n"<<effectif_cluster[l];}
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         //cout<<"Number singularities : "<<nb1<<endl;
00828         //cout<<"Number singularities PostRemoval : "<<nb2<<endl;
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         //INVERSE COVARIANCE POSTREMOVAL ESTIMATION
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                                 //INVERSE of COV MATRIX
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                                 //DETERMINANT OF THE INVERSE COV MATRIX
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          //if (ratio<0){
00955          //cout<<Proba_Assigne_Gauss(i,firstNearestA[i])-Proba_Assigne_Gauss_PostRemoval(i,secondNearestA[i],firstNearestA[i])<<endl;
00956          //}
00957          delete [] auxligne;
00958          return ratio;
00959 }
00960 /*--------------------------------------------------------------------------------------------*/
00961 float Inliers::Error_FreeInliers_Singularity(unsigned long int i)
00962 {
00963         //independence between bands assumption for simplication!!!
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                         //bands without a singularity
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                       //bands with a singularity (except case "singularity->singularity" : zero cost exclusion)
00979                       else
00980                       {
00981                         if (covariance_cluster_PostRemoval[secondNearestA[i]][firstNearestA[i]][j][j]!=0)//singularity->no singularity 
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)// no singularity->singularity 
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         //variation error brougth by the new free inliers
01055         for (i=0;i<nb_inliers;i++)
01056         {
01057                 ratio=0;
01058                 
01059                 //cost error free inlier when no singularity
01060                 if ((singularity[firstNearestA[i]]==0) && (Singularity_PostRemoval[secondNearestA[i]][firstNearestA[i]]==0))
01061                 {       
01062                         ratio=Error_FreeInliers(i);                                                             
01063                 }
01064                 // cost  error free inlier when singularities
01065                 else 
01066                 {       //cout<<singularity[firstNearestA[i]]<<" && "<<Singularity_PostRemoval[secondNearestA[i]][firstNearestA[i]]<<endl;
01067                         ratio=Error_FreeInliers_Singularity(i);
01068                 
01069                 }
01070                 
01071                 SumErrorRatio[firstNearestA[i]]=SumErrorRatio[firstNearestA[i]]+ratio;
01072                 
01073         }
01074         
01075         //variation error brougth by the others inliers du to the change of their variance
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;// factor for calculating the variation of the encoding length ofteh Gaussian paramters. THIS IS A VERY ROUGH APPROXIMATION AND SHOULD BE REPLACED BY THE TERM DEFINED IN THEORITICAL NOTE  
01099                 Inliers_Variation_Index_term_cost2(Indexcost2);
01100 
01101                 //search of the greatest decrease
01102                 maxidecrease=Inf;
01103                 indexmaxidecrease=0;//default
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                 //suppression of the greatest decrease                                          
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 }                                

Generated on Thu Feb 17 11:01:08 2005 for Gaussian Mixture Modeling by the MDL principle by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002