Main Page   Compound List   File List   Compound Members   File Members   Related Pages  

graphCharac.cpp

Go to the documentation of this file.
00001 
00003 #include <iostream>
00004 #include <string>  
00005 #include <vector>
00006 #include <fstream>
00007 
00008 #include <cstdio>
00009 #include <cstdlib>
00010 #include <cmath>
00011 #include <cassert>
00012 
00013 
00014 #include "graphCharac.hh"
00015 
00016 
00017 using namespace std;
00018 
00019 /*****************************************************************************************************/
00020 GraphDynaCluster::GraphDynaCluster()
00021 {
00022 
00023 }
00024 /*****************************************************************************************************/
00025 GraphDynaCluster::~GraphDynaCluster()
00026 {
00027         
00028 }
00029 /*****************************************************************************************************/
00030 void GraphDynaCluster::Initialise(char * path_in_Temporal, char * path_in_Timeless, int dimx_tmp, int dimy_tmp, int nb_bands)
00031 {
00032         nb_dim=nb_bands;
00033         dim=dimx_tmp*dimy_tmp;
00034         nx=dimx_tmp;
00035         ny=dimy_tmp;
00036         
00037         Acquire_path(path_in_Temporal,nb_dim,path_in_Timeless); 
00038         Allocate_variables();
00039         Aquire_data();
00040 } 
00041 /*****************************************************************************************************/ 
00042 void GraphDynaCluster::Acquire_path(char * path_in_Temporal, int nb_bands, char* path_in_Timeless)  
00043 {
00044         
00045         int m,k;
00046         Image<unsigned char> pippo;
00047         nb_images=pippo.NumberImages_pathfiles_in(path_in_Temporal,2);
00048         
00049         //auxiliary paths
00050         char ***PathInTemp= new char **[nb_images];
00051         for(m=0;m<nb_images;++m){PathInTemp[m] =new char*[2];} 
00052         for(m=0;m<nb_images;++m){for(k=0;k<2;++k){PathInTemp[m][k]=new char[128];}}
00053         //read 
00054         pippo.Read_pathfiles_in(path_in_Temporal,PathInTemp,2);
00055         
00056         char ***PathInTimeless= new char **[1];
00057         for(m=0;m<1;++m){PathInTimeless[m] =new char*[2];} 
00058         for(m=0;m<1;++m){for(k=0;k<2;++k){PathInTimeless[m][k]=new char[128];}}
00059         //read
00060         pippo.Read_pathfiles_in(path_in_Timeless,PathInTimeless,2);
00061         
00062         //allocation path
00063         PathInTemp_Classfile= new char **[nb_images];
00064         for(m=0;m<nb_images;++m){PathInTemp_Classfile[m] =new char*[1];} 
00065         for(m=0;m<nb_images;++m){for(k=0;k<1;++k){PathInTemp_Classfile[m][k]=new char[128];}}
00066         
00067         PathInTimeless_Classfile= new char **[1];
00068         for(m=0;m<1;++m){PathInTimeless_Classfile[m] =new char*[1];} 
00069         for(m=0;m<1;++m){for(k=0;k<1;++k){PathInTimeless_Classfile[m][k]=new char[128];}}
00070         
00071         PathInTemp_MeanCov= new char **[nb_images];
00072         for(m=0;m<nb_images;++m){PathInTemp_MeanCov[m] =new char*[1];} 
00073         for(m=0;m<nb_images;++m){for(k=0;k<1;++k){PathInTemp_MeanCov[m][k]=new char[128];}}
00074         
00075         PathInTimeless_MeanCov= new char **[1];
00076         for(m=0;m<1;++m){PathInTimeless_MeanCov[m] =new char*[1];} 
00077         for(m=0;m<1;++m){for(k=0;k<1;++k){PathInTimeless_MeanCov[m][k]=new char[128];}}
00078 
00079         
00080         //copy auxiliary paths in paths
00081         for (m=0;m<nb_images;m++)
00082         {
00083                 strcpy(PathInTemp_Classfile[m][0],PathInTemp[m][0]);
00084                 strcpy(PathInTemp_MeanCov[m][0],PathInTemp[m][1]);
00085         }
00086         strcpy(PathInTimeless_Classfile[0][0],PathInTimeless[0][0]);
00087         strcpy(PathInTimeless_MeanCov[0][0],PathInTimeless[0][1]);
00088         
00089         
00090         
00091         //delete auxiliary paths
00092         for(m=0;m<nb_images;++m){for(k=0;k<2;++k){delete [] PathInTemp[m][k];}}
00093         for(m=0;m<nb_images;++m){delete [] PathInTemp[m];}
00094         for(k=0;k<2;++k){delete [] PathInTimeless[0][k];}
00095         
00096 
00097         
00098 }
00099 /****************************************************************************************************/
00100 void GraphDynaCluster::Allocate_variables()
00101 {
00102         //aux_variables allocation
00103         Current_C_A=new float*[nb_dim]; Current_C_B=new float*[nb_dim];Current_C_AB=new float*[nb_dim*2]; 
00104         Current_M_A=new float[nb_dim]; Current_M_B=new float[nb_dim]; Current_M_AB=new float[nb_dim*2];
00105         
00106         for (int j=0;j<nb_dim;j++)
00107         {
00108         Current_C_A[j]=new float[nb_dim];
00109         Current_C_B[j]=new float[nb_dim];
00110         }
00111         for (int j=0;j<nb_dim*2;j++){Current_C_AB[j]=new float[nb_dim*2];}
00112         
00113         //classfiles allocation
00114         Temp_classfile.Allocate(nx,ny,1,nb_images);
00115         Timeless_classfile.Allocate(nx,ny,1,1);
00116         
00117         //attributes allocation
00118         Temp_MeanCov.Allocate(255,(nb_dim+1)*nb_dim,1,nb_images);
00119         Timeless_MeanCov.Allocate(255,(nb_dim*nb_images+1)*nb_dim*nb_images,1,1);
00120         
00121         
00122         //SuperposeClass alllocation
00123         supClass=new SuperposeClass[nb_images];
00124                 
00125         
00126 }
00127 /*****************************************************************************************************/
00128 void GraphDynaCluster::Aquire_data()
00129 {
00130         //reading attributes 
00131         Temp_MeanCov.Read_Cut_SubSample(PathInTemp_MeanCov,0,0,Temp_MeanCov.ny,1,1);
00132         Timeless_MeanCov.Read_Cut_SubSample(PathInTimeless_MeanCov,0,0,Timeless_MeanCov.ny,1,1);
00133         
00134         //reading classfiles
00135         Temp_classfile.Read_Cut_SubSample(PathInTemp_Classfile,0,0,ny,1,1);
00136         Timeless_classfile.Read_Cut_SubSample(PathInTimeless_Classfile,0,0,ny,1,1);
00137 }
00138 /*****************************************************************************************************/
00139 void GraphDynaCluster::Load_Current(int j1,int j2)
00140 {
00141         int k,l;
00142         unsigned long int dim_MeanCov=Temp_MeanCov.nx*Temp_MeanCov.ny;
00143         for(k=0;k<nb_dim;k++)
00144         {
00145                 Current_M_A[k]=(*(Temp_MeanCov.position+dim_MeanCov*Index_CurrentImage+(nb_dim+1)*nb_dim*(j1+1)+k));
00146                 
00147                 Current_M_B[k]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(j2+1)+Index_CurrentImage*nb_dim+k));
00148                 if (k==0){Current_M_A[k]/=alpha0;Current_M_B[k]/=alpha0;}
00149                 if (k==1){Current_M_A[k]/=alpha1;Current_M_B[k]/=alpha1;}
00150                 if (k==2){Current_M_A[k]/=alpha2;Current_M_B[k]/=alpha2;}
00151                 for (l=0;l<nb_dim;l++)
00152                 {
00153                         Current_C_A[k][l]=(*(Temp_MeanCov.position+dim_MeanCov*Index_CurrentImage+(nb_dim+1)*nb_dim*(j1+1)+(k+1)*nb_dim+l));
00154                         Current_C_B[k][l]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(j2+1)+(k+1+Index_CurrentImage*nb_dim)*nb_dim*nb_images+l+Index_CurrentImage*nb_dim));
00155                         if (k==0){Current_C_A[k][l]/=alpha0;Current_C_B[k][l]/=alpha0;}
00156                         if (k==1){Current_C_A[k][l]/=alpha1;Current_C_B[k][l]/=alpha1;}
00157                         if (k==2){Current_C_A[k][l]/=alpha2;Current_C_B[k][l]/=alpha2;}
00158                         if (l==0){Current_C_A[k][l]/=alpha0;Current_C_B[k][l]/=alpha0;}
00159                         if (l==1){Current_C_A[k][l]/=alpha1;Current_C_B[k][l]/=alpha1;} 
00160                         if (l==2){Current_C_A[k][l]/=alpha2;Current_C_B[k][l]/=alpha2;} 
00161                 }
00162                 
00163         }
00164 
00165 }
00166 /*****************************************************************************************************/
00167 void GraphDynaCluster::Build(float Factor)
00168 {
00169         int j1,j2,j,k;
00170         int nb_CurrentClasses;
00171         PrecisionFactor=Factor;
00172         
00173         Divergence DivObject(nb_dim);
00174         float *min_bound=new float[nb_dim];float *max_bound=new float[nb_dim];
00175         float *resolution=new float[nb_dim];
00176         float DivResult,MinDivResult;
00177         int IndexMinDiv;
00178         
00179         
00180         nb_TimelessClasses=(int)*(Timeless_MeanCov.position);
00181         nb_CurrentClasses=(int)*(Temp_MeanCov.position);
00182         
00183         //Dynamic classe objects allocation (only for nearest TL class)
00184          TmpClass=new DynaClass[nb_TimelessClasses];
00185         for (j1=0;j1<nb_TimelessClasses;j1++)
00186         {       
00187                 TmpClass[j1].Allocate(nb_images,nb_dim);
00188         }
00189         
00190         
00191         //Divergence Calculation
00192         for (Index_CurrentImage=0;Index_CurrentImage<nb_images;Index_CurrentImage++)
00193         {
00194                 supClass[Index_CurrentImage].AcquirePara(Index_CurrentImage,nb_dim,nb_images,nb_TimelessClasses,&Temp_classfile,&Timeless_classfile);
00195                 
00196                 nb_CurrentClasses=(int)*(Temp_MeanCov.position+Temp_MeanCov.nx*Temp_MeanCov.ny*Index_CurrentImage);
00197                 cout<<"Number of classes at time "<<Index_CurrentImage+1<<" : "<<nb_CurrentClasses<<endl;
00198                 float *VectorDiv=new float[255];  //max nb of classes: 255!!
00199                 int *VectorDivIndex=new int[255];
00200                 float boucle=0;
00201                 float auxINF=2.*INF;
00202                 int CalculateDiv;
00203                 float DivWeight;
00204                 for (j2=0;j2<nb_TimelessClasses;j2++)
00205                 {       
00206                 
00207                         for (j1=0;j1<nb_CurrentClasses;j1++){
00208                                 VectorDiv[j1]=auxINF;
00209                         }
00210                         //Searching the divergence 
00211                         MinDivResult=INF;
00212                         IndexMinDiv=-1;
00213                         for (j1=0;j1<nb_CurrentClasses;j1++){
00214                                 Load_Current(j1,j2);
00215                                 DivObject.Initialise (Current_C_A, Current_C_B, Current_M_A, Current_M_B );
00216                                 DivObject.InvertCovMatrix();
00217                                 CalculateDiv=DivObject.BoundResolution_Calculation(min_bound,max_bound,resolution,PrecisionFactor);                                     DivWeight=(float)supClass[Index_CurrentImage].DivWeight(j1,j2)  ;
00218                                 if (CalculateDiv<1||DivWeight!=0){
00219                                    DivResult=-1;
00220                                    boucle=0;
00221                                    while (DivResult<0 && boucle<3){
00222                                      boucle++; 
00223                                      DivResult=DivObject.Divergence2GaussModel(min_bound,max_bound,resolution,DivWeight);// spatial constrain introduce via DivWeight.
00224                                       for (int jaux=0;jaux<nb_dim;jaux++){resolution[jaux]=resolution[jaux]*0.5;}
00225                                    }
00226                                    if (DivResult<=0){
00227                                         DivResult=DivObject.MeanDifference();
00228                                   }
00229                                 }
00230                                 else{DivResult=INF;}
00231                                 if (DivWeight==0){DivResult=INF;}
00232                                 if (DivResult<MinDivResult){MinDivResult=DivResult;IndexMinDiv=j1;}
00233                                 if (DivResult>=auxINF){DivResult=INF;}
00234                                 for (int jj=0;jj<nb_CurrentClasses;jj++){
00235                                         if (DivResult<VectorDiv[jj]){
00236                                                 for (int jjj=nb_CurrentClasses;jjj>jj;jjj--){ 
00237                                                         VectorDiv[jjj]=VectorDiv[jjj-1];
00238                                                         VectorDivIndex[jjj]=VectorDivIndex[jjj-1];
00239                                                         
00240                                                 }
00241                                                 VectorDiv[jj]=DivResult;
00242                                                 VectorDivIndex[jj]=j1;
00243                                                 break;
00244                                         }       
00245                                 }
00246                                 
00247                                         
00248                         }
00249                         
00250                         cout<<"multitemporal class : "<<j2<<" assigned to timelocalized class "<<IndexMinDiv<<"; measurement : "<<MinDivResult<<endl;
00251                         
00252 
00253                         //Dynamic classe objects allocation (exept for nearest TL class)
00254                         for (int jj=1;jj<nb_CurrentClasses;jj++){
00255                                 int *AuxAssignedTimelessClass=new int[1];
00256                                 unsigned long int *AuxWeight=new unsigned long int[1];
00257                                 float  *AuxDivergenceValue=new float[1];
00258                                 float  **AuxCentroids=new float*[1];
00259                                 float ***AuxCovariances=new float**[1];
00260                                 AuxCentroids[0]=new float[nb_dim];
00261                                 AuxCovariances[0]=new float*[ nb_dim];
00262                                 for (j=0;j<nb_dim;j++)
00263                                 {
00264                                         AuxCovariances[0][j]=new float[nb_dim];
00265                                 }
00266                                 
00267                                 TmpClass[j2].AssignedTimelessClass[Index_CurrentImage].push_back(AuxAssignedTimelessClass[0]);
00268                                 TmpClass[j2].DivergenceValue[Index_CurrentImage].push_back(AuxDivergenceValue[0]);
00269                                 TmpClass[j2].Centroids[Index_CurrentImage].push_back(AuxCentroids[0]);
00270                                 TmpClass[j2].Covariances[Index_CurrentImage].push_back(AuxCovariances[0]);
00271                                 TmpClass[j2].Weight[Index_CurrentImage].push_back(AuxWeight[0]);
00272                                 
00273                         }
00274                         
00275                         for (int jj=0;jj<nb_CurrentClasses;jj++){
00276                         
00277                                 TmpClass[j2].AssignedTimelessClass[Index_CurrentImage][jj]=VectorDivIndex[jj];
00278                                 TmpClass[j2].Weight[Index_CurrentImage][jj]=(unsigned long int)supClass[Index_CurrentImage].DivWeight(VectorDivIndex[jj],j2);
00279                                 //cout<<"timelocalized class "<<VectorDivIndex[jj]<<"; Divergence : "<<VectorDiv[jj]<<" "<<TmpClass[j2].Weight[Index_CurrentImage][jj]<<"";"";
00280                                 
00281                                 TmpClass[j2].DivergenceValue[Index_CurrentImage][jj]=VectorDiv[jj];
00282                                 Load_Current(VectorDivIndex[jj],j2);
00283                                 for (j=0;j<nb_dim;j++)
00284                                 {
00285                                         TmpClass[j2].Centroids[Index_CurrentImage][jj][j]=Current_M_A[j];
00286                                         //cout<<" "<<Current_M_A[j]<<"/"<<Current_M_B[j];
00287                                         for (k=0;k<nb_dim;k++)
00288                                         {
00289                                            TmpClass[j2].Covariances[Index_CurrentImage][jj][j][k]=Current_C_A[j][k];
00290                                         }
00291                                         //cout<<","<<Current_C_A[j][j]<<"/"<<Current_C_B[j][j];
00292                                 }
00293                                 //cout<<""<<endl;
00294                                 
00295                         }
00296                         
00297                         
00298                 }
00299         
00300         }
00301 }
00302 /*****************************************************************************************************/ 
00303 void GraphDynaCluster::Out_Evolutions(int format, bool intensiveCalculus)
00304 {
00305     
00306      if (format==2){//format graph where the mean and covariance are stored many times (for the actual implementation of the learning procedure and  of the GUI) 
00307         Create_TxtGraphOutPutVersion2(intensiveCalculus);
00308         ReadWriteGraphVersion2JAVA("./grapheV2",1);
00309      } 
00310      if (format==3){//format graph where the mean and covariance are stored many times (for the actual implementation of the learning procedure and creation of a format without edge for the GUI) 
00311         Create_TxtGraphOutPutVersion2(intensiveCalculus);
00312         ReadWriteGraphVersion2JAVA("./grapheV2",0);
00313      }
00314      else{//Standard graph format output
00315         if (format!=2){
00316                 Create_TxtGraphOutPutVersion1(intensiveCalculus);
00317         }
00318      }
00319 
00320 }
00321 /*****************************************************************************************************/ 
00322 void GraphDynaCluster::Create_TxtGraphOutPutVersion1(bool intensiveCalculus )
00323 {
00324 
00325         char *path=new char[128];
00326         FILE *ifp;
00327         strcpy(path,"./grapheV1");
00328         ifp=NULL; 
00329         assert (ifp=fopen(path,"w"));
00330         int length;
00331         float flow;
00332         float weight;
00333         //float auxInf=INF;
00334         Divergence divObject(nb_dim);
00335         float *min_bound=new float[nb_dim];float *max_bound=new float[nb_dim];
00336         float *resolution=new float[nb_dim];float *resolutionAux=new float[nb_dim];
00337         int calculateDiv=0;float DivResult=0;
00338         fprintf(ifp,"%d %s",(int) nb_TimelessClasses," ");//nb MT clusters
00339         fprintf(ifp,"%d %s",(int) nb_dim," ");//nb dimension FS
00340         fprintf(ifp,"%d %s",(int) nb_images," ");//nb images samples
00341         //for each MT class
00342         for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00343             cout<<"MT class "<<kMT<<endl;
00344             fprintf(ifp,"%d %s",kMT," ");//MT index
00345             //for each TL window
00346             for(int t=0;t<nb_images;t++){
00347                 fprintf(ifp,"%d %s",t," ");//time sample
00348                 Index_CurrentImage=t;
00349                 Load_MonoCurrentPMT(kMT);
00350                 length=(int)TmpClass[kMT].AssignedTimelessClass[t].size();
00351                 fprintf(ifp,"%d %s",length," ");//number of assigned  TLs
00352                 //for each TL_t class
00353                 for (int kTL=0;kTL<(int)length;kTL++){ 
00354                         fprintf(ifp,"%d %s", (int)TmpClass[kMT].AssignedTimelessClass[t][kTL]," ");//TL_t cluster index
00355                         fprintf(ifp,"%u %s", (unsigned int)TmpClass[kMT].Weight[t][kTL]," ");//pixel quantity belonging to classes kTL_t & kMT
00356                         fprintf(ifp,"%e %s", (float)TmpClass[kMT].DivergenceValue[t][kTL]*weight," ");//divergence between classes kTL_t & kMT
00357                                                         
00358                         if (t>0){
00359                             //for each assigne TL_(t-1)class
00360                             for (int kTLpre=0;kTLpre<(int)TmpClass[kMT].AssignedTimelessClass[t-1].size();kTLpre++){
00361                                 fprintf(ifp,"%d %s",(int)TmpClass[kMT].AssignedTimelessClass[t-1][kTLpre]," ");//TL_(t-1) cluster index
00362                                 flow=(unsigned long int)supClass[t].FlowWeight(TmpClass[kMT].AssignedTimelessClass[t][kTL],TmpClass[kMT].AssignedTimelessClass[t-1][kTLpre],kMT);
00363                                 fprintf(ifp,"%u %s",(unsigned int) flow," ");//pixel quantity belonging to classes kTL_t&kTL_(t-1)&kMT : flow
00364                                 LoadpreTLTL(TmpClass[kMT].AssignedTimelessClass[t-1][kTLpre],TmpClass[kMT].AssignedTimelessClass[t][kTL]);
00365                                 
00366                                 divObject.Initialise(Current_C_A, Current_C_B, Current_M_A, Current_M_B);
00367                                 divObject.InvertCovMatrix();
00368                                 calculateDiv=divObject.BoundResolution_Calculation(min_bound,max_bound,resolution,PrecisionFactor);
00369                                 if (intensiveCalculus==true && calculateDiv<1){ //Divergence between consecutive TL classes (intensive calculation)                             
00370                                         DivResult=-1;
00371                                         for (int jaux=0;jaux<nb_dim;jaux++){resolutionAux[jaux]=resolution[jaux];}
00372                                         while (DivResult<0&&(float)resolution[0]/resolutionAux[0]<=(float)1){
00373                                           DivResult=divObject.Divergence2GaussModel(min_bound,max_bound,resolutionAux,1.);//div TL(t-1) TLt
00374                                           for (int jaux=0;jaux<nb_dim;jaux++){resolutionAux[jaux]=resolutionAux[jaux]*0.5;}
00375                                         }
00376                                         if (DivResult<=0){DivResult=divObject.MeanDifference();}
00377                                         //DivResult=sqrt(DivResult*DivResult);
00378                                         fprintf(ifp,"%e %s",DivResult," ");//div TL(t-1) TLt
00379                                   
00380                                 }
00381                                 else{//Mean difference between consecutive TL classes (fast calculation)
00382                                         DivResult=divObject.MeanDifference();
00383                                         fprintf(ifp,"%e %s",DivResult," ");
00384                                 }
00385                             }
00386                         }
00387                 }
00388 
00389                 if (t>0){ fprintf(ifp,"%f %s", (float)TmpClass[kMT].MutualInfo[t-1]," ");}// Mutual information of MT class between t(t-1)
00390             }
00391         }
00392         assert (EOF!=fclose(ifp)); 
00393         delete [] min_bound; delete [] max_bound;
00394         delete [] resolution;
00395                 
00396 }
00397 /*****************************************************************************************************/ 
00398 void GraphDynaCluster::Create_TxtGraphOutPutVersion2(bool intensiveCalculus )
00399 {
00400         char *path=new char[128];
00401         FILE *ifp;
00402         strcpy(path,"./grapheV2");
00403         ifp=NULL; 
00404         assert (ifp=fopen(path,"w"));
00405         int length;
00406         float flow;
00407         float weight;
00408         //float auxInf=INF;
00409         Divergence divObject(nb_dim);
00410         float *min_bound=new float[nb_dim];float *max_bound=new float[nb_dim];
00411         float *resolution=new float[nb_dim];float *resolutionAux=new float[nb_dim];
00412         int calculateDiv=0;float DivResult=0;
00413         fprintf(ifp,"%d %s",(int) nb_TimelessClasses," ");//nb MT clusters
00414         fprintf(ifp,"%d %s",(int) nb_dim," ");//nb dimension FS
00415         fprintf(ifp,"%d %s",(int) nb_images," ");//nb images samples
00416         //for each MT class
00417         for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00418             cout<<"MT class "<<kMT<<endl;
00419             fprintf(ifp,"%d %s",kMT," ");//MT index
00420             weight=(float)supClass[0].MTPopulation(kMT);
00421             fprintf(ifp,"%u %s",(unsigned int)weight," ");//MT weight
00422             //for each TL window
00423             for(int t=0;t<nb_images;t++){
00424                 fprintf(ifp,"%d %s",t," ");//time sample
00425                 Index_CurrentImage=t;
00426                 Load_MonoCurrentPMT(kMT);
00427                 //mean and cov kMT
00428                 for (int j=0;j<nb_dim;j++){fprintf(ifp,"%f %s",Current_M_A[j]," ");}
00429                 for (int j=0;j<nb_dim;j++){
00430                         for (int l=0;l<nb_dim;l++){
00431                                 fprintf(ifp,"%f %s",Current_C_A[j][l]," ");
00432                         }
00433                 }
00434                 length=(int)TmpClass[kMT].AssignedTimelessClass[t].size();
00435                 fprintf(ifp,"%d %s",length," ");//number of assigned  TLs
00436                 //for each TL_t class
00437                 for (int kTL=0;kTL<(int)length;kTL++){ 
00438                         fprintf(ifp,"%d %s", (int)TmpClass[kMT].AssignedTimelessClass[t][kTL]," ");//TL_t cluster index
00439                         fprintf(ifp,"%u %s", (unsigned int)TmpClass[kMT].Weight[t][kTL]," ");//pixel quantity belonging to classes kTL_t & kMT
00440                         fprintf(ifp,"%e %s", (float)TmpClass[kMT].DivergenceValue[t][kTL]*weight," ");//divergence between classes kTL_t & kMT
00441                         //mean and cov TL_t
00442                         for (int i=0;i<nb_dim;i++){fprintf(ifp,"%f %s",TmpClass[kMT].Centroids[t][kTL][i]," ");}
00443                         for (int i=0;i<nb_dim;i++){
00444                                 for (int j=0;j<nb_dim;j++){fprintf(ifp,"%f %s",TmpClass[kMT].Covariances[t][kTL][i][j]," ");}
00445                         }                               
00446                         if (t>0){
00447                             //for each assigne TL_(t-1)class
00448                             for (int kTLpre=0;kTLpre<(int)TmpClass[kMT].AssignedTimelessClass[t-1].size();kTLpre++){
00449                                 fprintf(ifp,"%d %s",(int)TmpClass[kMT].AssignedTimelessClass[t-1][kTLpre]," ");//TL_(t-1) cluster index
00450                                 flow=(unsigned long int)supClass[t].FlowWeight(TmpClass[kMT].AssignedTimelessClass[t][kTL],TmpClass[kMT].AssignedTimelessClass[t-1][kTLpre],kMT);
00451                                 fprintf(ifp,"%u %s",(unsigned int) flow," ");//pixel quantity belonging to classes kTL_t&kTL_(t-1)&kMT : flow
00452                                 LoadpreTLTL(TmpClass[kMT].AssignedTimelessClass[t-1][kTLpre],TmpClass[kMT].AssignedTimelessClass[t][kTL]);
00453                                 
00454                                 divObject.Initialise(Current_C_A, Current_C_B, Current_M_A, Current_M_B);
00455                                 divObject.InvertCovMatrix();
00456                                 calculateDiv=divObject.BoundResolution_Calculation(min_bound,max_bound,resolution,PrecisionFactor);
00457                                 if (intensiveCalculus==true && calculateDiv<1){
00458                                         
00459                                         DivResult=-1;
00460                                         for (int jaux=0;jaux<nb_dim;jaux++){resolutionAux[jaux]=resolution[jaux];}
00461                                         while (DivResult<0&&(float)resolution[0]/resolutionAux[0]<=(float)1){
00462                                           DivResult=divObject.Divergence2GaussModel(min_bound,max_bound,resolutionAux,1.);//div TL(t-1) TLt
00463                                           for (int jaux=0;jaux<nb_dim;jaux++){resolutionAux[jaux]=resolutionAux[jaux]*0.5;}
00464                                         }
00465                                         if (DivResult<=0){DivResult=divObject.MeanDifference();}
00466                                         //DivResult=sqrt(DivResult*DivResult);
00467                                         fprintf(ifp,"%e %s",DivResult," ");//div TL(t-1) TLt
00468                                   
00469                                 }
00470                                 else{
00471                                         DivResult=divObject.MeanDifference();
00472                                         fprintf(ifp,"%e %s",DivResult," ");
00473                                 }
00474                             }
00475                         }
00476                 }
00477 
00478                 if (t>0){ fprintf(ifp,"%f %s", (float)TmpClass[kMT].MutualInfo[t-1]," ");}// Mutual information of MT class between t(t-1)
00479             }
00480         }
00481         assert (EOF!=fclose(ifp)); 
00482         delete [] min_bound; delete [] max_bound;
00483         delete [] resolution;
00484                 
00485 }
00486 /*****************************************************************************************************/ 
00487 void GraphDynaCluster::ReadWriteGraphVersion2JAVA(char * path, int edge)
00488 {
00489         
00490         FILE *ifp;
00491         ifp=NULL; 
00492         assert (ifp=fopen(path,"r"));
00493         int auxRead;int length;float auxReadF;float weightMT;   
00494         
00495         char *pathw=new char[128];
00496         FILE *ifpw;
00497         strcpy(pathw,"./grapheV2JAVA");
00498         ifpw=NULL; 
00499         assert (ifpw=fopen(pathw,"w"));
00500         
00501         fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00502         nb_TimelessClasses=auxRead;//nb MT clusters
00503         fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00504         nb_dim=auxRead;//nb dimension FS
00505         fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00506         nb_images=auxRead;//nb images samples
00507         
00508         
00509         //aux_variables allocation
00510         Current_C_A=new float*[nb_dim]; Current_C_B=new float*[nb_dim];Current_C_AB=new float*[nb_dim*2]; 
00511         Current_M_A=new float[nb_dim]; Current_M_B=new float[nb_dim]; Current_M_AB=new float[nb_dim*2];
00512         for (int j=0;j<nb_dim;j++)
00513         {
00514         Current_C_A[j]=new float[nb_dim];
00515         Current_C_B[j]=new float[nb_dim];
00516         }
00517         for (int j=0;j<nb_dim*2;j++){Current_C_AB[j]=new float[nb_dim*2];}
00518         
00519         float *min_bound=new float[nb_dim];float *max_bound=new float[nb_dim];
00520         float *resolution=new float[nb_dim];
00521                 
00522         //Dynamic classe objects allocation
00523          TmpClass=new DynaClass[nb_TimelessClasses];
00524         for (int kMT=0;kMT<nb_TimelessClasses;kMT++){TmpClass[kMT].InitRead(nb_images,nb_dim);}
00525         
00526         //for each MT class
00527         for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00528             fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");//MT index
00529             fscanf(ifp,"%f",&weightMT);fprintf(ifpw,"%u %s",(unsigned int) weightMT," ");//MT weight   
00530             //for each TL window
00531             for(int t=0;t<nb_images;t++){
00532                 fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");//time sample
00533                 TmpClass[kMT].MTAlloc(t);
00534                 TmpClass[kMT].MTWeight[t][0]=(unsigned long int)weightMT;
00535                 //mean and cov kMT
00536                 for (int j=0;j<nb_dim;j++){
00537                         fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%f %s",(float) auxReadF," ");
00538                         TmpClass[kMT].MTCentroids[t][0][j]=auxReadF;
00539                 }
00540                 for (int j=0;j<nb_dim;j++){
00541                         for (int l=0;l<nb_dim;l++){
00542                                 fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%f %s",(float) auxReadF," ");
00543                                 TmpClass[kMT].MTCovariances[t][0][j][l]=auxReadF;
00544                         }
00545                 } 
00546                 fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00547                 length=auxRead;//number of assigned  TLs
00548                 TmpClass[kMT].TLAlloc(t,length);
00549                 if (t>0){TmpClass[kMT].flowDivAlloc(t,(int)TmpClass[kMT].AssignedTimelessClass[t-1].size(),length);}
00550                 //for each TL_t class
00551                 for (int kTL=0;kTL<(int)length;kTL++){ 
00552                         fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00553                         TmpClass[kMT].AssignedTimelessClass[t][kTL]=auxRead;//TL_t cluster index
00554                         fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%u %s",(unsigned int) auxReadF," ");
00555                         TmpClass[kMT].Weight[t][kTL]=(unsigned long int)auxReadF;//pixel quantity belonging to classes kTL_t & kMT
00556                         fscanf(ifp,"%e",&auxReadF);fprintf(ifpw,"%e %s",(double) auxReadF," ");
00557                         TmpClass[kMT].DivergenceValue[t][kTL]=(float)auxReadF/(float)TmpClass[kMT].MTWeight[t][0];//divergence between classes kTL_t & kMT
00558                         //mean and cov TL_t
00559                         for (int i=0;i<nb_dim;i++){
00560                                 fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%f %s",(float) auxReadF," ");
00561                                 TmpClass[kMT].Centroids[t][kTL][i]=auxReadF;
00562                         }
00563                         for (int i=0;i<nb_dim;i++){
00564                                 for (int j=0;j<nb_dim;j++){
00565                                         fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%f %s",(float) auxReadF," ");
00566                                         TmpClass[kMT].Covariances[t][kTL][i][j]=auxReadF;
00567                                 }
00568                         }
00569                         if (t>0){
00570                            //for each assigne TL_(t-1)class
00571                            for (int kTLpre=0;kTLpre<(int)TmpClass[kMT].AssignedTimelessClass[t-1].size();kTLpre++){
00572                                 fscanf(ifp,"%d",&auxRead);
00573                                 if (edge==1){fprintf(ifpw,"%d %s",(int) auxRead," ");}//TL_(t-1) cluster index
00574                                 fscanf(ifp,"%f",&auxReadF);
00575                                 if (edge==1){fprintf(ifpw,"%u %s",(unsigned int) auxReadF," ");}
00576                                 TmpClass[kMT].flow[t][kTLpre][kTL]=auxReadF;//pixel quantity belonging to classes kTL_t&kTL_(t-1)&kMT : flow
00577                                 fscanf(ifp,"%f",&auxReadF);
00578                                 if (edge==1){fprintf(ifpw,"%e %s",(double) auxReadF," ");}
00579                                 TmpClass[kMT].div[t][kTLpre][kTL]=auxReadF;//div TL(t-1) TLt
00580                            }
00581                            
00582                         }
00583                         fprintf(ifpw,"\n");
00584                 }
00585 
00586                 if (t>0){       
00587                         fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%e %s\n",(double) auxReadF," ");
00588                         TmpClass[kMT].MutualInfo[t-1]=auxReadF;// Mutual information of MT class between t(t-1)
00589                 }
00590             }
00591         }
00592         dim=0;
00593         for (int kMT=0;kMT<nb_TimelessClasses;kMT++){dim+=TmpClass[kMT].MTWeight[0][0];}
00594         
00595         assert (EOF!=fclose(ifp)); 
00596         delete [] min_bound; delete [] max_bound;
00597         delete [] resolution;
00598                 
00599 }
00600 
00601 /*****************************************************************************************************/
00602 void GraphDynaCluster::Load_CurrentPMT(int kMT)
00603 {
00604         for(int k=0;k<nb_dim;k++)
00605         {
00606                 Current_M_A[k]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+(Index_CurrentImage-1)*nb_dim+k));
00607                 Current_M_B[k]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+Index_CurrentImage*nb_dim+k));
00608                 if (k==0){Current_M_A[k]/=alpha0;Current_M_B[k]/=alpha0;}
00609                 if (k==1){Current_M_A[k]/=alpha1;Current_M_B[k]/=alpha1;}
00610                 if (k==2){Current_M_A[k]/=alpha2;Current_M_B[k]/=alpha2;}
00611                 for (int l=0;l<nb_dim;l++)
00612                 {
00613                         Current_C_A[k][l]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+(k+1+(Index_CurrentImage-1)*nb_dim)*nb_dim*nb_images+l+(Index_CurrentImage-1)*nb_dim));
00614                         Current_C_B[k][l]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+(k+1+Index_CurrentImage*nb_dim)*nb_dim*nb_images+l+Index_CurrentImage*nb_dim));
00615                         if (k==0){Current_C_A[k][l]/=alpha0;Current_C_B[k][l]/=alpha0;}
00616                         if (k==1){Current_C_A[k][l]/=alpha1;Current_C_B[k][l]/=alpha1;}
00617                         if (k==2){Current_C_A[k][l]/=alpha2;Current_C_B[k][l]/=alpha2;}
00618                         if (l==0){Current_C_A[k][l]/=alpha0;Current_C_B[k][l]/=alpha0;} 
00619                         if (l==1){Current_C_A[k][l]/=alpha1;Current_C_B[k][l]/=alpha1;}
00620                         if (l==2){Current_C_A[k][l]/=alpha2;Current_C_B[k][l]/=alpha2;} 
00621                         
00622                 }
00623         }
00624         for(int k=0;k<nb_dim;k++)
00625         {
00626                 Current_M_AB[k]=Current_M_A[k];
00627                 Current_M_AB[k+nb_dim]=Current_M_B[k];
00628         }
00629         for(int k=0;k<nb_dim*2;k++){
00630                 for (int l=0;l<nb_dim*2;l++){
00631                 
00632 Current_C_AB[k][l]=*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+(k+1+(Index_CurrentImage-1)*nb_dim)*nb_dim*nb_images+l+(Index_CurrentImage-1)*nb_dim);
00633                         if (k==0||k==3){Current_C_AB[k][l]/=alpha0}
00634                         if (k==1||k==4){Current_C_AB[k][l]/=alpha1}
00635                         if (k==2||k==5){Current_C_AB[k][l]/=alpha2}
00636                         if (l==0||l==3){Current_C_AB[k][l]/=alpha0}     
00637                         if (l==1||l==4){Current_C_AB[k][l]/=alpha1}
00638                         if (l==2||l==5){Current_C_AB[k][l]/=alpha2}     
00639                 }
00640         }
00641 
00642 }
00643 /*****************************************************************************************************/
00644 void GraphDynaCluster::Load_MonoCurrentPMT(int kMT)
00645 {
00646         for(int k=0;k<nb_dim;k++)
00647         {
00648                 Current_M_A[k]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+Index_CurrentImage*nb_dim+k));
00649                 if (k==0){Current_M_A[k]/=alpha0}
00650                 if (k==1){Current_M_A[k]/=alpha1}
00651                 if (k==2){Current_M_A[k]/=alpha2}
00652                 for (int l=0;l<nb_dim;l++)
00653                 {
00654                                                 Current_C_A[k][l]=(*(Timeless_MeanCov.position+(nb_dim*nb_images+1)*nb_dim*nb_images*(kMT+1)+(k+1+Index_CurrentImage*nb_dim)*nb_dim*nb_images+l+Index_CurrentImage*nb_dim));
00655                         if (k==0){Current_C_A[k][l]/=alpha0}
00656                         if (k==1){Current_C_A[k][l]/=alpha1}
00657                         if (k==2){Current_C_A[k][l]/=alpha2}
00658                         if (l==0){Current_C_A[k][l]/=alpha0}    
00659                         if (l==1){Current_C_A[k][l]/=alpha1}    
00660                         if (l==2){Current_C_A[k][l]/=alpha2}    
00661                         
00662                 }
00663         }
00664 
00665 }
00666 /*****************************************************************************************************/ 
00667 void GraphDynaCluster::LoadpreTLTL(int prekTL,int kTL){
00668 
00669         unsigned long int dim_MeanCov=Temp_MeanCov.nx*Temp_MeanCov.ny;
00670         for(int k=0;k<nb_dim;k++)
00671         {
00672                 Current_M_A[k]=(*(Temp_MeanCov.position+dim_MeanCov*(Index_CurrentImage-1)+(nb_dim+1)*nb_dim*(prekTL+1)+k));
00673                 Current_M_B[k]=(*(Temp_MeanCov.position+dim_MeanCov*Index_CurrentImage+(nb_dim+1)*nb_dim*(kTL+1)+k));
00674                 if (k==0){Current_M_A[k]/=alpha0;Current_M_B[k]/=alpha0;}
00675                 if (k==1){Current_M_A[k]/=alpha1;Current_M_B[k]/=alpha1;}
00676                 if (k==2){Current_M_A[k]/=alpha2;Current_M_B[k]/=alpha2;}
00677                 for (int l=0;l<nb_dim;l++)
00678                 {
00679                         Current_C_A[k][l]=(*(Temp_MeanCov.position+dim_MeanCov*(Index_CurrentImage-1)+(nb_dim+1)*nb_dim*(prekTL+1)+(k+1)*nb_dim+l));
00680                         Current_C_B[k][l]=(*(Temp_MeanCov.position+dim_MeanCov*Index_CurrentImage+(nb_dim+1)*nb_dim*(kTL+1)+(k+1)*nb_dim+l));
00681                         if (k==0){Current_C_A[k][l]/=alpha0;Current_C_B[k][l]/=alpha0;}
00682                         if (k==1){Current_C_A[k][l]/=alpha1;Current_C_B[k][l]/=alpha1;}
00683                         if (k==2){Current_C_A[k][l]/=alpha2;Current_C_B[k][l]/=alpha2;}
00684                         if (l==0){Current_C_A[k][l]/=alpha0;Current_C_B[k][l]/=alpha0;} 
00685                         if (l==1){Current_C_A[k][l]/=alpha1;Current_C_B[k][l]/=alpha1;}
00686                         if (l==2){Current_C_A[k][l]/=alpha2;Current_C_B[k][l]/=alpha2;} 
00687                 }
00688         }
00689 
00690 }
00691 
00692 /*****************************************************************************************************/ 
00693 void  GraphDynaCluster::MIEstimation()
00694 {       
00695         
00696         for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00697                 for (int t=1;t<nb_images;t++){  
00698                         Index_CurrentImage=t;
00699                         Load_CurrentPMT(kMT);
00700                         TmpClass[kMT].MTWeight[t-1][0]=1;
00701                         TmpClass[kMT].MTWeight[t][0]=1;
00702                         TmpClass[kMT].JointWeight[t-1][0]=1;
00703                         for (int j=0;j<nb_dim;j++){
00704                                 TmpClass[kMT].MTCentroids[t-1][0][j]=Current_M_A[j];
00705                                 TmpClass[kMT].MTCentroids[t][0][j]=Current_M_B[j];
00706                                 TmpClass[kMT].JointCentroids[t-1][0][j]=Current_M_AB[j];
00707                                 TmpClass[kMT].JointCentroids[t-1][0][j+nb_dim]=Current_M_AB[j+nb_dim];
00708                                 for (int l=0;l<nb_dim;l++){
00709                                         TmpClass[kMT].MTCovariances[t-1][0][j][l]=Current_C_A[j][l];
00710                                         TmpClass[kMT].MTCovariances[t][0][j][l]=Current_C_B[j][l];
00711                                 }
00712                         }
00713                         for (int j=0;j<nb_dim*2;j++){
00714                                 for (int l=0;l<nb_dim*2;l++){
00715                                         TmpClass[kMT].JointCovariances[t-1][0][j][l]=Current_C_AB[j][l];
00716                                 }
00717                         }
00718                         
00719                 }
00720         }
00721 
00722 cout<<"Calculating the mutual Information for each MT Class between each consecutive times"<<endl;
00723 
00724         
00725         for (int t=1;t<nb_images;t++){  
00726             mi.Allocate(nb_dim);
00727             for (int kMT=0;kMT<nb_TimelessClasses;kMT++){       
00728                 //cout<<"t "<<t-1<<t<<" , kMT "<<kMT<<endl;         
00729                 vector<float**> *C1=&TmpClass[kMT].MTCovariances[t-1];
00730                 vector<float*> *M1=&TmpClass[kMT].MTCentroids[t-1];     
00731                 vector<float**> *C2=&TmpClass[kMT].MTCovariances[t];
00732                 vector<float*> *M2=&TmpClass[kMT].MTCentroids[t];               
00733                 vector<float**> *C12=&TmpClass[kMT].JointCovariances[t-1];
00734                 vector<float*> *M12=&TmpClass[kMT].JointCentroids[t-1];
00735                 mi.Initialise1(C1,M1);
00736                 mi.Initialise2(C2,M2);
00737                 mi.Initialise12(C12,M12);
00738                 TmpClass[kMT].MutualInfo[t-1]=mi.MICalculus();
00739                 
00740             }
00741             
00742         }
00743 cout<<"Mutual Information end"<<endl;
00744         
00745 }
00746 /*****************************************************************************************************/ 
00747 
00748 
00749 

Generated on Thu Feb 17 11:01:55 2005 for Inference of a Graph of Dynamic Cluster Trajectories by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002