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
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
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
00060 pippo.Read_pathfiles_in(path_in_Timeless,PathInTimeless,2);
00061
00062
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
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
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
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
00114 Temp_classfile.Allocate(nx,ny,1,nb_images);
00115 Timeless_classfile.Allocate(nx,ny,1,1);
00116
00117
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
00123 supClass=new SuperposeClass[nb_images];
00124
00125
00126 }
00127
00128 void GraphDynaCluster::Aquire_data()
00129 {
00130
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
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
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
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];
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
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);
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
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
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
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
00292 }
00293
00294
00295 }
00296
00297
00298 }
00299
00300 }
00301 }
00302
00303 void GraphDynaCluster::Out_Evolutions(int format, bool intensiveCalculus)
00304 {
00305
00306 if (format==2){
00307 Create_TxtGraphOutPutVersion2(intensiveCalculus);
00308 ReadWriteGraphVersion2JAVA("./grapheV2",1);
00309 }
00310 if (format==3){
00311 Create_TxtGraphOutPutVersion2(intensiveCalculus);
00312 ReadWriteGraphVersion2JAVA("./grapheV2",0);
00313 }
00314 else{
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
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," ");
00339 fprintf(ifp,"%d %s",(int) nb_dim," ");
00340 fprintf(ifp,"%d %s",(int) nb_images," ");
00341
00342 for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00343 cout<<"MT class "<<kMT<<endl;
00344 fprintf(ifp,"%d %s",kMT," ");
00345
00346 for(int t=0;t<nb_images;t++){
00347 fprintf(ifp,"%d %s",t," ");
00348 Index_CurrentImage=t;
00349 Load_MonoCurrentPMT(kMT);
00350 length=(int)TmpClass[kMT].AssignedTimelessClass[t].size();
00351 fprintf(ifp,"%d %s",length," ");
00352
00353 for (int kTL=0;kTL<(int)length;kTL++){
00354 fprintf(ifp,"%d %s", (int)TmpClass[kMT].AssignedTimelessClass[t][kTL]," ");
00355 fprintf(ifp,"%u %s", (unsigned int)TmpClass[kMT].Weight[t][kTL]," ");
00356 fprintf(ifp,"%e %s", (float)TmpClass[kMT].DivergenceValue[t][kTL]*weight," ");
00357
00358 if (t>0){
00359
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]," ");
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," ");
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){
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.);
00374 for (int jaux=0;jaux<nb_dim;jaux++){resolutionAux[jaux]=resolutionAux[jaux]*0.5;}
00375 }
00376 if (DivResult<=0){DivResult=divObject.MeanDifference();}
00377
00378 fprintf(ifp,"%e %s",DivResult," ");
00379
00380 }
00381 else{
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]," ");}
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
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," ");
00414 fprintf(ifp,"%d %s",(int) nb_dim," ");
00415 fprintf(ifp,"%d %s",(int) nb_images," ");
00416
00417 for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00418 cout<<"MT class "<<kMT<<endl;
00419 fprintf(ifp,"%d %s",kMT," ");
00420 weight=(float)supClass[0].MTPopulation(kMT);
00421 fprintf(ifp,"%u %s",(unsigned int)weight," ");
00422
00423 for(int t=0;t<nb_images;t++){
00424 fprintf(ifp,"%d %s",t," ");
00425 Index_CurrentImage=t;
00426 Load_MonoCurrentPMT(kMT);
00427
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," ");
00436
00437 for (int kTL=0;kTL<(int)length;kTL++){
00438 fprintf(ifp,"%d %s", (int)TmpClass[kMT].AssignedTimelessClass[t][kTL]," ");
00439 fprintf(ifp,"%u %s", (unsigned int)TmpClass[kMT].Weight[t][kTL]," ");
00440 fprintf(ifp,"%e %s", (float)TmpClass[kMT].DivergenceValue[t][kTL]*weight," ");
00441
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
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]," ");
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," ");
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.);
00463 for (int jaux=0;jaux<nb_dim;jaux++){resolutionAux[jaux]=resolutionAux[jaux]*0.5;}
00464 }
00465 if (DivResult<=0){DivResult=divObject.MeanDifference();}
00466
00467 fprintf(ifp,"%e %s",DivResult," ");
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]," ");}
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;
00503 fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00504 nb_dim=auxRead;
00505 fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00506 nb_images=auxRead;
00507
00508
00509
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
00523 TmpClass=new DynaClass[nb_TimelessClasses];
00524 for (int kMT=0;kMT<nb_TimelessClasses;kMT++){TmpClass[kMT].InitRead(nb_images,nb_dim);}
00525
00526
00527 for (int kMT=0;kMT<nb_TimelessClasses;kMT++){
00528 fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00529 fscanf(ifp,"%f",&weightMT);fprintf(ifpw,"%u %s",(unsigned int) weightMT," ");
00530
00531 for(int t=0;t<nb_images;t++){
00532 fscanf(ifp,"%d",&auxRead);fprintf(ifpw,"%d %s",(int) auxRead," ");
00533 TmpClass[kMT].MTAlloc(t);
00534 TmpClass[kMT].MTWeight[t][0]=(unsigned long int)weightMT;
00535
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;
00548 TmpClass[kMT].TLAlloc(t,length);
00549 if (t>0){TmpClass[kMT].flowDivAlloc(t,(int)TmpClass[kMT].AssignedTimelessClass[t-1].size(),length);}
00550
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;
00554 fscanf(ifp,"%f",&auxReadF);fprintf(ifpw,"%u %s",(unsigned int) auxReadF," ");
00555 TmpClass[kMT].Weight[t][kTL]=(unsigned long int)auxReadF;
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];
00558
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
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," ");}
00574 fscanf(ifp,"%f",&auxReadF);
00575 if (edge==1){fprintf(ifpw,"%u %s",(unsigned int) auxReadF," ");}
00576 TmpClass[kMT].flow[t][kTLpre][kTL]=auxReadF;
00577 fscanf(ifp,"%f",&auxReadF);
00578 if (edge==1){fprintf(ifpw,"%e %s",(double) auxReadF," ");}
00579 TmpClass[kMT].div[t][kTLpre][kTL]=auxReadF;
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;
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
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