src/MEDMEM/INTERPOLATION/MEDMEM_dTree.hxx

   1 #ifndef MEDMEM_DTREE_HXX
   2 #define MEDMEM_DTREE_HXX
   3
   4 #include "MEDMEM_dTreeSommet.hxx"
   5
   6 #define DTREE_FANTOME -1
   7 #define DTREE_RACINE 0
   8 #define DTREE_TERMINAL 1
   9 #define DTREE_COURANT 2
  10 // etat_descendance
  11 #define DTREE_NON_CALCULE -1
  12 #define DTREE_AUCUNE 0
  13 #define DTREE_VALIDE 1
  14 // val min nbr noeud
  15 #define DTREE_NBR_MIN_NOEUDS 1
  16 #define DTREE_NBR_MAX_DESC 8
  17 // pour meilleu re lecture
  18 #define _TEMPLATE_ template <class NOEUD,class NUAGENOEUD,int DIMENSION,int NBR_NOEUDS_PAR_CASE>
  19 #define _DTREE_ dTree<NOEUD,NUAGENOEUD,DIMENSION,NBR_NOEUDS_PAR_CASE>
  20
  21 //////////////////////////////////////////////////////////////////
  22 ///                                                            ///
  23 ///                       DECLARATIONS                         ///
  24 ///                                                            ///
  25 //////////////////////////////////////////////////////////////////
  26
  27 /*********************************************************/
  28 /*                                                       */
  29 /*             Calcul Statique de Puissance              */
  30 /*                                                       */
  31 /*********************************************************/
  32
  33 // Permet de Calculer 2^n de façon statique
  34
  35 template < int DIMENSION > struct DeuxPuissance
  36         {
  37         enum { valeur = 2 * DeuxPuissance<DIMENSION-1>::valeur } ;
  38         };
  39
  40 template <> struct DeuxPuissance<0>
  41         {
  42         enum { valeur = 1  } ;
  43         };
  44
  45 /*********************************************************/
  46 /*                                                       */
  47 /*                         DTREE                         */
  48 /*                                                       */
  49 /*********************************************************/
  50
  51 // Classe Utilisateur
  52 // # Bugs connus :
  53 //      -  Problemes avec points sur les faces des hypercubes ?
  54 //      -  Plantage sauvage si le point est plus loin de l'hypercube père que le diametre Linf de celui-ci.
  55 // # Politique de classe :
  56 //      - NOEUD : voir dans MEDMEM_WrapperNodes.hxx la classe Wrapper_Noeud<..>
  57 //      - NUAGENOEUD : typiquement, c'est vector<NOEUD> en rédefinissant NUAGENOEUD<...>::SIZE()=vector<...>::size()
  58 // Remarques :
  59 //      - NBR_NOEUDS_PAR_CASE ne doit pas être modifié sauf peut-être dans le cas où l'utilisateur veut utiliser des d-Tree parallèles
  60
  61 template <class NOEUD,class NUAGENOEUD,int DIMENSION,int NBR_NOEUDS_PAR_CASE=DTREE_NBR_MIN_NOEUDS> class dTree
  62 {
  63 protected :
  64         // types
  65         typedef _DTREE_ * Ptr_dTree;
  66         // Static Const
  67         static const int nbr_descendants=DeuxPuissance<DIMENSION>::valeur;
  68         // champs
  69         NUAGENOEUD * nuage;
  70         Ptr_dTree descendant[nbr_descendants];
  71         vector<int> * noeud_contenu;
  72         int etat;
  73         int niveau;
  74         dTree * pere;
  75         Sommet_dTree<DIMENSION> coord_max;
  76         Sommet_dTree<DIMENSION> coord_min;
  77 public :
  78
  79         void init();
  80         dTree();
  81         dTree(NUAGENOEUD *n);
  82         dTree(const Sommet_dTree<DIMENSION> &A,const Sommet_dTree<DIMENSION> &B,dTree *mypere);
  83         dTree(const dTree &F);
  84         ~dTree();
  85         void Get_Noeuds_Filtre(vector<int> &tmp);
  86         Sommet_dTree<DIMENSION> Get_Max() const;
  87         Sommet_dTree<DIMENSION> Get_Min() const;
  88         int is_in_dTree(NOEUD P) const;
  89         double calcule_distance(NOEUD P) const;
  90         dTree & operator = (const dTree & F);
  91         Sommet_dTree<DIMENSION> dTree::donne_sommet(int selecteur) const;
  92         int a_des_fils() const;
  93         dTree * trouve_dTree_contenant(NOEUD P) const;
  94         int trouve_plus_proche_point_bourrin(NOEUD P) const;
  95         int trouve_plus_proche_point(NOEUD P) const;
  96         int trouve_un_point() const;
  97         int tppp_rec(NOEUD P,double &delta,int &flag) const;
  98         int Localise_Point(NOEUD P,double d) const;
  99         void cree_filiation();
 100         int Get_Nbr_Descendants_Non_Vides() const;
 101         int Get_Nbr_Descendants_Vides() const;
 102         int Get_Profondeur_Max() const;
 103 };
 104
 105
 106 //////////////////////////////////////////////////////////////////
 107 ///                                                            ///
 108 ///                           CODE                             ///
 109 ///                                                            ///
 110 //////////////////////////////////////////////////////////////////
 111
 112 _TEMPLATE_ void _DTREE_::init()
 113         {
 114         int i;
 115         nuage=NULL;
 116         noeud_contenu=NULL;
 117         etat=DTREE_FANTOME;
 118         for (i=0;i<nbr_descendants;i++) descendant[i]=NULL;
 119         niveau=-1;
 120         pere=NULL;
 121         }
 122 _TEMPLATE_ _DTREE_::dTree()
 123         {
 124         init();
 125         coord_max=0.0;
 126         coord_min=0.0;
 127         }
 128 _TEMPLATE_ _DTREE_::dTree(NUAGENOEUD *n)
 129         {
 130
 131         int i,j;
 132         double tmp;
 133
 134         if (n==NULL)
 135                 {
 136                 cerr<<"DTREE : Nuage vide !"<<endl;
 137                 exit(-1);
 138                 }
 139
 140         init();
 141         nuage=n;
 142         etat=DTREE_RACINE;
 143         noeud_contenu=new vector<int>(nuage->size());
 144         niveau=0;
 145
 146         for (i=0;i<DIMENSION;i++)
 147                 {
 148                 coord_max[i]=(*nuage)[0][i];
 149                 coord_min[i]=(*nuage)[0][i];
 150                 }
 151         (*noeud_contenu)[0]=0;
 152         for (i=1;i<nuage->size();i++)
 153                 {
 154                 (*noeud_contenu)[i]=i;
 155                 for (j=0;j<DIMENSION;j++)
 156                         {
 157                         if ((*nuage)[i][j]>coord_max[j]) coord_max[j]=(*nuage)[i][j];
 158                         if ((*nuage)[i][j]<coord_min[j]) coord_min[j]=(*nuage)[i][j];
 159                         }
 160                 }
 161         for (j=0;j<DIMENSION;j++)
 162                 {
 163                 tmp=(coord_max[j]-coord_min[j])*0.01;
 164                 coord_max[j]+=tmp;
 165                 coord_min[j]-=tmp;
 166                 }
 167
 168         cree_filiation();
 169
 170         }
 171 _TEMPLATE_ _DTREE_::dTree(const Sommet_dTree<DIMENSION> &A,const Sommet_dTree<DIMENSION> &B,dTree *mypere)
 172         {
 173         if (mypere!=NULL)
 174                 {
 175                 int i,j;
 176                 double tmp;
 177
 178                 init();
 179
 180                 pere=mypere;
 181                 nuage=pere->nuage;
 182                 niveau=pere->niveau+1;
 183
 184                 etat=DTREE_COURANT;
 185
 186                 noeud_contenu=new vector<int>;
 187                 noeud_contenu->reserve((pere->noeud_contenu->size())/nbr_descendants);
 188
 189                 for (i=0;i<DIMENSION;i++)
 190                         {
 191                         coord_max[i]=(*nuage)[0][i];
 192                         coord_min[i]=(*nuage)[0][i];
 193                         }
 194
 195                 for (i=0;i<DIMENSION;i++)
 196                         {
 197                         if (A[i]>B[i])
 198                                 {
 199                                 coord_max[i]=A[i];
 200                                 coord_min[i]=B[i];
 201                                 }
 202                         else
 203                                 {
 204                                 coord_min[i]=A[i];
 205                                 coord_max[i]=B[i];
 206                                 }
 207                         }
 208                 }
 209         else
 210                 {
 211                 cerr<<"DTREE : Construction de descendance sans père !"<<endl;
 212                 exit(-1);
 213                 }
 214         }
 215 _TEMPLATE_ _DTREE_::dTree(const _DTREE_ &F)
 216         {
 217         // Pas Super Top Moumoute ... Recopie de pointeurs et pas de contenus, merdique
 218         int i,j;
 219         init();
 220         for (i=0;i<nbr_descendants;i++) descendant[i]=F.descendant[i];
 221         noeud_contenu=F.noeud_contenu;
 222         etat=F.etat;
 223         niveau=F.niveau;
 224         pere=F.pere;
 225         coord_max=F.coord_max;
 226         coord_min=F.coord_min;
 227         }
 228 _TEMPLATE_ _DTREE_::~dTree()
 229         {
 230         int i;
 231         nuage=NULL;
 232         pere=NULL;
 233         // cout<<*this<<endl;
 234         if (noeud_contenu)
 235                 {
 236                 delete noeud_contenu;
 237                 }
 238         for (i=0;i<nbr_descendants;i++) if (descendant[i])
 239                 {
 240                 delete descendant[i];
 241                 }
 242         }
 243 _TEMPLATE_ void _DTREE_::Get_Noeuds_Filtre(vector<int> &tmp)
 244         {
 245         int i;
 246         switch (etat)
 247                 {
 248                 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
 249                 case DTREE_COURANT :
 250                         for (i=0;i<nbr_descendants;i++) descendant[i]->Get_Noeuds_Filtre(tmp);
 251                 case DTREE_TERMINAL :
 252                         if (noeud_contenu->size()>0)
 253                                 {
 254                                 for (i=0;i<NBR_NOEUDS_PAR_CASE;i++) tmp.push_back((*noeud_contenu)[i]);
 255                                 }
 256                 default :
 257                         cerr<<"DTREE : dTree Non Valide dans Rempli_Noeuds_Filtre"<<endl;
 258                         exit(-1);
 259                 }
 260         }
 261 _TEMPLATE_ Sommet_dTree<DIMENSION> _DTREE_::Get_Max() const
 262         {
 263         return coord_max;
 264         }
 265 _TEMPLATE_ Sommet_dTree<DIMENSION> _DTREE_::Get_Min() const
 266         {
 267         return coord_min;
 268         }
 269 _TEMPLATE_ int _DTREE_::is_in_dTree(NOEUD P) const
 270         {
 271         int test;
 272         for (int i=0;i<DIMENSION;i++)
 273                 {
 274                 test=((coord_min[i]<=P[i])&&(P[i]<=coord_max[i]));
 275                 if (!test) return 0;
 276                 }
 277         return 1;
 278
 279         }
 280 // calcule la distance Linf d'un point exterieur, si le point est interieur, le resultat sera erroné
 281 _TEMPLATE_ double _DTREE_::calcule_distance(NOEUD P) const
 282         {
 283         double min=fabs(coord_max[0]-P[0]);
 284         double tmp;
 285         int i;
 286         for (i=0;i<DIMENSION;i++)
 287                 {
 288                 tmp=fabs(coord_max[i]-P[i]);
 289                 if (min>tmp) min=tmp;
 290                 tmp=fabs(coord_min[i]-P[i]);
 291                 if (min>tmp) min=tmp;
 292                 }
 293         return min;
 294         }
 295 _TEMPLATE_ _DTREE_ & _DTREE_::operator = (const _DTREE_ & F)
 296         {
 297         // Pas Super Top Moumoute ... Recopie de pointeurs et pas de contenus, merdique
 298         int i,j;
 299         init();
 300         for (i=0;i<nbr_descendans;i++) descendant[i]=F.descendant[i];
 301         noeud_contenu=F.noeud_contenu;
 302         etat=F.etat;
 303         niveau=F.niveau;
 304         pere=F.pere;
 305         coord_max=F.coord_max;
 306         coord_min=F.coord_min;
 307         }
 308 // selecteur doit etre lu comme un nombre en base 2
 309 // il specifie les coordonnes de reference du sommet dont on veut les coordonnees reelles
 310 // ex : en dim 3 : coord_min=0+0*2+0*4=0
 311 //                 coord_max=1+1*2*1*4=7
 312 // donc les unites pour x, les 2aines pour y, les 4aines pour z
 313 _TEMPLATE_ Sommet_dTree<DIMENSION> _DTREE_::donne_sommet(int selecteur) const
 314         {
 315         Sommet_dTree<DIMENSION> p;
 316         int i;
 317         int residu=selecteur;
 318         int reste;
 319         for (i=0;i<DIMENSION;i++)
 320                 {
 321                 reste=residu%2;
 322                 if (reste==0) p[i]=coord_min[i]; else p[i]=coord_max[i];
 323                 residu=(residu-reste)/2;
 324                 }
 325         return p;
 326         }
 327 _TEMPLATE_ int _DTREE_::a_des_fils() const
 328         {
 329         if ((etat==DTREE_COURANT)||(etat=DTREE_RACINE)) return 1;
 330         else return 0;
 331         }
 332 _TEMPLATE_ _DTREE_ * _DTREE_::trouve_dTree_contenant(NOEUD P) const
 333         {
 334         Sommet_dTree<DIMENSION> B(coord_min,coord_max);
 335         int i;
 336         if ((etat==DTREE_RACINE)&&(!is_in_dTree(P))) return NULL; // Le noeud est extérieur a l'dTree
 337         else if (etat==DTREE_TERMINAL) return (dTree *) this;     // Le noeud est dans *this qui est terminal
 338
 339         for (i=0;i<nbr_descendants;i++)
 340                 {
 341                 Sommet_dTree<DIMENSION> A=donne_sommet(i);
 342                 int test,j;
 343                 for (j=0;j<DIMENSION;j++)
 344                         {
 345                         if (A[j]<=B[j]) test=((A[j]<=P[j])&&(P[j]<=B[j]));
 346                         else test=((A[j]>=P[j])&&(P[j]>=B[j]));
 347                         if (!test) break;
 348                         }
 349
 350                 if (test) return descendant[i]->trouve_dTree_contenant(P); // Propagation
 351                 }
 352         }
 353 // si de le dTree n'est pas TERMINAL, scanne tous les points du nuage du pere pour trouver le point le plus proche
 354 // sinon scanne uniquement les points contenus dans le dTree
 355 _TEMPLATE_ int _DTREE_::trouve_plus_proche_point_bourrin(NOEUD P) const
 356         {
 357         int i;
 358         int num_sol;
 359         double min;
 360         double tmp;
 361         if (etat!=DTREE_TERMINAL)
 362                 {
 363                 min=DistanceL2(P,(*nuage)[0]);
 364                 for (i=1;i<nuage->size();i++)
 365                         {
 366                         tmp=DistanceL2(P,(*nuage)[i]);
 367                         if (tmp<min)
 368                                 {
 369                                 num_sol=i;
 370                                 min=tmp;
 371                                 }
 372                         }
 373                 return num_sol;
 374                 }
 375         else
 376                 {
 377                 if (noeud_contenu->size()!=0)
 378                         {
 379                         num_sol=(*noeud_contenu)[0];
 380                         min=DistanceL2(P,(*nuage)[num_sol]);
 381                         for (i=1;i<noeud_contenu->size();i++)
 382                                 {
 383                                 tmp=DistanceL2(P,(*nuage)[(*noeud_contenu)[i]]);
 384                                 if (tmp<min)
 385                                         {
 386                                         num_sol=(*noeud_contenu)[i];
 387                                         min=tmp;
 388                                         }
 389                                 }
 390                         return num_sol;
 391                         }
 392                 else return DTREE_FANTOME;
 393                 }
 394         }
 395 // Fonction de pilotage de la recursion, lance les tests preliminaires et invoque la recursion
 396 _TEMPLATE_ int _DTREE_::trouve_plus_proche_point(NOEUD P) const
 397         {
 398
 399         if (etat==DTREE_FANTOME)
 400                 {
 401                 cerr<<"DTREE : TROUVE PLUS PROCHE POINT : L'octree n'est pas construit !"<<endl;
 402                 exit(-1);
 403                 }
 404         dTree *ref=trouve_dTree_contenant(P);
 405         double delta;
 406
 407         if (ref!=NULL)
 408                 {
 409                 // Le point est intérieur
 410
 411                 if (ref->pere!=NULL)
 412                         {
 413                         delta=DistanceInf((ref->pere->coord_max),(ref->pere->coord_min));
 414                         }
 415                 else
 416                         {
 417                         cerr<<"DTREE : TROUVE PLUS PROCHE POINT : l'octree contenant le noeud n'a pas de pere !"<<endl;
 418                         exit(-1);
 419                         }
 420                 }
 421         else
 422                 {
 423                 // Le point est extérieur
 424
 425                 delta=DistanceL2(P,(*nuage)[0]);
 426                 }
 427
 428         int flag=0;
 429         int retour;
 430
 431         retour=tppp_rec(P,delta,flag);
 432
 433         return retour;
 434         }
 435 _TEMPLATE_ int _DTREE_::trouve_un_point() const
 436         {
 437         if (noeud_contenu!=NULL)
 438                 {
 439                 if (noeud_contenu->size()>0) return (*(noeud_contenu))[0];
 440                 }
 441         else
 442                 {
 443                 int i;
 444                 for (i=0;i<nbr_descendants;i++)
 445                         {
 446                         return descendant[i]->trouve_un_point();
 447                         }
 448                 }
 449         }
 450 // partie recursive de l'algo ci dessus
 451 // change le flag en 1 si un point plus proche est trouvé
 452 // adapte automatiquement la distance delta de filtrage
 453 _TEMPLATE_ int _DTREE_::tppp_rec(NOEUD P,double &delta,int &flag) const
 454         {
 455         if (Localise_Point(P,delta)==0)
 456                 {
 457                 // La distance du point à l'octree est plus grande que delta
 458                 return DTREE_FANTOME;
 459                 }
 460         int num_Ptmp;
 461         double dtmp;
 462         if (etat==DTREE_TERMINAL)
 463                 {
 464                 if (noeud_contenu->size()>0)
 465                         {
 466                         num_Ptmp=trouve_plus_proche_point_bourrin(P);
 467                         dtmp=DistanceL2((*nuage)[num_Ptmp],P);
 468                         if (dtmp<=delta)
 469                                 {
 470                                 // Le point le plus proche minimise delta, c'est un bon candidat, on l'envoie !
 471                                 delta=dtmp;
 472                                 flag=1;
 473                                 return num_Ptmp;
 474                                 }
 475                         // Le point le plus proche ne minimise pas delta
 476                         // ===========> peut etre rajouter exit(-1); ??????????
 477                         return DTREE_FANTOME;
 478                         }
 479                 else
 480                         {
 481                         // L'octree qui contient P ne contient aucun point
 482                         return DTREE_FANTOME;
 483                         }
 484                 }
 485         int i;
 486         int num_sol=DTREE_FANTOME;
 487         for (i=0;i<nbr_descendants;i++)
 488                 {
 489                 num_Ptmp=descendant[i]->tppp_rec(P,delta,flag);
 490                 if ((num_Ptmp!=DTREE_FANTOME)&&(flag==1))
 491                         {
 492                         // On a trouvé un point qui minimise delta dans une branche
 493                         num_sol=num_Ptmp;
 494                         }
 495                 }
 496         // A ce stade, on a trouvé un point qui minimise tous les deltas, c'est donc le point le plus proche
 497         return num_sol;
 498         }
 499 // renvoie 1 si la distance L inf du noeud a l'octree est plus petite que d, 0 sinon
 500 _TEMPLATE_ int _DTREE_::Localise_Point(NOEUD P,double d) const
 501         {
 502         int i;
 503         for (i=0;i<DIMENSION;i++)
 504                 {
 505                 if (P[i]>coord_max[i]+d) return 0;
 506                 if (P[i]<coord_min[i]-d) return 0;
 507                 }
 508         return 1;
 509         }
 510 // Méthode de construction du dTree par propagation
 511 _TEMPLATE_ void _DTREE_::cree_filiation()
 512         {
 513         if (etat!=DTREE_RACINE)
 514                 {
 515                 niveau=pere->niveau+1;
 516                 }
 517         else
 518                 {
 519                 niveau=0;
 520                 }
 521
 522         if (noeud_contenu->size()<=NBR_NOEUDS_PAR_CASE)
 523                 {
 524                 etat=DTREE_TERMINAL;
 525                 }
 526         else
 527                 {
 528                 int i,num_loc,test;
 529
 530                 Sommet_dTree<DIMENSION> centre(coord_max,coord_min);
 531
 532                 for (i=0;i<nbr_descendants;i++)
 533                         {
 534                         descendant[i]=new dTree(centre,donne_sommet(i),this);
 535                         }
 536
 537                 for (num_loc=0;num_loc<noeud_contenu->size();num_loc++)
 538                         {
 539                         int indice=(*noeud_contenu)[num_loc];
 540                         NOEUD & courant=(*nuage)[indice];
 541                         test=1;
 542                         for (i=0;(test)&&(i<nbr_descendants);i++)
 543                                 {
 544                                 if (descendant[i]->is_in_dTree(courant))
 545                                         {
 546                                         descendant[i]->noeud_contenu->push_back(indice);
 547                                         test=0;
 548                                         }
 549                                 }
 550                         }
 551
 552                 delete noeud_contenu;
 553                 noeud_contenu=NULL;
 554
 555                 for (i=0;i<nbr_descendants;i++) descendant[i]->cree_filiation();
 556                 }
 557         }
 558 _TEMPLATE_ int _DTREE_::Get_Nbr_Descendants_Non_Vides() const
 559         {
 560         int i;
 561         int ndnv=0;
 562         switch (etat)
 563                 {
 564                 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
 565                 case DTREE_COURANT :
 566                         for (i=0;i<nbr_descendants;i++) ndnv+=descendant[i]->Get_Nbr_Descendants_Non_Vides();
 567                         return ndnv;
 568                 case DTREE_TERMINAL :
 569                         if (noeud_contenu->size()>0) return 1;
 570                         else return 0;
 571                 default :
 572                         cerr<<"dTree Non Valide dans Get_Nbr_Descendants_Non_Vides"<<endl;
 573                         return -1;
 574                 }
 575         }
 576 _TEMPLATE_ int _DTREE_::Get_Nbr_Descendants_Vides() const
 577         {
 578         int i;
 579         int ndnv=0;
 580         switch (etat)
 581                 {
 582                 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
 583                 case DTREE_COURANT :
 584                         for (i=0;i<nbr_descendants;i++) ndnv+=descendant[i]->Get_Nbr_Descendants_Vides();
 585                         return ndnv;
 586                 case DTREE_TERMINAL :
 587                         if (noeud_contenu->size()==0) return 1;
 588                         else return 0;
 589                 default :
 590                         cerr<<"dTree Non Valide dans Get_Nbr_Descendants_Non_Vides"<<endl;
 591                         return -1;
 592                 }
 593         }
 594 _TEMPLATE_ int _DTREE_::Get_Profondeur_Max() const
 595         {
 596         int i;
 597         int prof=0;
 598         int tmp;
 599         switch (etat)
 600                 {
 601                 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
 602                 case DTREE_COURANT :
 603                         for (i=0;i<nbr_descendants;i++)
 604                                 {
 605                                 tmp=descendant[i]->Get_Profondeur_Max();
 606                                 if (tmp>prof) prof=tmp;
 607                                 }
 608                         return prof;
 609                 case DTREE_TERMINAL : return niveau;
 610                 default :
 611                         cerr<<"dTree Non Valide dans Get_Nbr_Descendants_Non_Vides"<<endl;
 612                         return -1;
 613                 }
 614         }
 615
 616 #undef _TEMPLATE_
 617 #undef _DTREE_
 618
 619
 620 #endif