1 #ifndef MEDMEM_DTREE_HXX
2 #define MEDMEM_DTREE_HXX
4 #include "MEDMEM_dTreeSommet.hxx"
6 #define DTREE_FANTOME -1
8 #define DTREE_TERMINAL 1
9 #define DTREE_COURANT 2
11 #define DTREE_NON_CALCULE -1
12 #define DTREE_AUCUNE 0
13 #define DTREE_VALIDE 1
15 #define DTREE_NBR_MIN_NOEUDS 2
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>
21 //////////////////////////////////////////////////////////////////
25 //////////////////////////////////////////////////////////////////
27 /*********************************************************/
29 /* Calcul Statique de Puissance */
31 /*********************************************************/
33 // Permet de Calculer 2^n de façon statique
35 template < int DIMENSION > struct DeuxPuissance
37 enum { valeur = 2 * DeuxPuissance<DIMENSION-1>::valeur } ;
40 template <> struct DeuxPuissance<0>
45 /*********************************************************/
49 /*********************************************************/
51 // Construit un d-Tree sur un nuage de noeud et pour un noeud donné donne lequel des noeuds du nuage est le plus proche
54 // - Problemes avec points sur les faces des hypercubes ?
55 // - Plantage sauvage si le point est plus loin de l'hypercube père que le diametre Linf de celui-ci.
56 // # Politique de classe :
57 // - NOEUD : voir dans MEDMEM_WrapperNodes.hxx la classe Wrapper_Noeud<..>
58 // - NUAGENOEUD : typiquement, c'est vector<NOEUD> en rédefinissant NUAGENOEUD<...>::SIZE()=vector<...>::size()
60 // - 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
61 // ou utilise des nuages de noeud trop grands
63 template <class NOEUD,class NUAGENOEUD,int DIMENSION,int NBR_NOEUDS_PAR_CASE=DTREE_NBR_MIN_NOEUDS> class dTree
67 typedef _DTREE_ * Ptr_dTree;
69 static const int nbr_descendants=DeuxPuissance<DIMENSION>::valeur;
72 Ptr_dTree descendant[nbr_descendants];
73 // numéro des noeuds contenus
74 vector<int> * noeud_contenu;
78 // extrémités de l'hypercube du dTree
79 Sommet_dTree<DIMENSION> coord_max;
80 Sommet_dTree<DIMENSION> coord_min;
85 // Ce constructeur est le seul constructeur utilisateur, il construit le père puis tous les descendants
87 // Ce Constructeur est utilisé par le précédent, pour un pere donné, avec deux extrémités données,
88 // il instantie un dTree sans en calculer les noeuds contenus
89 dTree(const Sommet_dTree<DIMENSION> &A,const Sommet_dTree<DIMENSION> &B,dTree *mypere);
90 dTree(const dTree &F);
92 // Renvoie les numéros de noeuds contenus dans le dTree
93 void Get_Noeuds_Filtre(vector<int> &tmp);
94 // renvoie les extrémités
95 Sommet_dTree<DIMENSION> Get_Max() const;
96 Sommet_dTree<DIMENSION> Get_Min() const;
97 // renvoie vrai si P est dans le dTree
98 int is_in_dTree(NOEUD P) const;
99 // calcule la distance topologique locale de P au dTree
100 double calcule_distance(NOEUD P) const;
101 dTree & operator = (const dTree & F);
102 // retourne le sommet du dTree codé par l'entier selecteur (voir explications dans le code)
103 Sommet_dTree<DIMENSION> donne_sommet(int selecteur) const;
104 int a_des_fils() const;
105 // renvoi une reference sur le dTree terminal descendant de this contenant P
106 dTree * trouve_dTree_contenant(NOEUD P) const;
107 // renvoie le point le plus proche de P dans this par la méthode exhaustive brutale
108 int trouve_plus_proche_point_bourrin(NOEUD P) const;
109 // renvoie le point le plus proche de p dans this par la méthode dtree
110 int trouve_plus_proche_point(NOEUD P) const;
111 // renvoie un numéro de point contenu dans this
112 int trouve_un_point() const;
113 // méthode récursive utilisée par trouve_plus_proche_point
114 int tppp_rec(NOEUD P,double &delta,int &flag) const;
115 // dit si P est d-proche de l'hypercube de this
116 int Localise_Point(NOEUD P,double d) const;
117 // utilisé par le constructeur pour créer tout la filiation du dTree
118 void cree_filiation();
119 // méthodes de mesure
120 int Get_Nbr_Descendants_Non_Vides() const;
121 int Get_Nbr_Descendants_Vides() const;
122 int Get_Profondeur_Max() const;
126 //////////////////////////////////////////////////////////////////
130 //////////////////////////////////////////////////////////////////
132 _TEMPLATE_ void _DTREE_::init()
138 for (i=0;i<nbr_descendants;i++) descendant[i]=NULL;
142 _TEMPLATE_ _DTREE_::dTree()
148 _TEMPLATE_ _DTREE_::dTree(NUAGENOEUD *n)
155 cerr<<"DTREE : Nuage vide !"<<endl;
162 noeud_contenu=new vector<int>(nuage->size());
165 // calcule les extrémités du dTree pere
166 for (i=0;i<DIMENSION;i++)
168 coord_max[i]=(*nuage)[0][i];
169 coord_min[i]=(*nuage)[0][i];
171 (*noeud_contenu)[0]=0;
172 for (i=1;i<nuage->size();i++)
174 (*noeud_contenu)[i]=i;
175 for (j=0;j<DIMENSION;j++)
177 if ((*nuage)[i][j]>coord_max[j]) coord_max[j]=(*nuage)[i][j];
178 if ((*nuage)[i][j]<coord_min[j]) coord_min[j]=(*nuage)[i][j];
181 for (j=0;j<DIMENSION;j++)
183 tmp=(coord_max[j]-coord_min[j])*0.01;
188 // méthode récursive qui construit la filiation
192 _TEMPLATE_ _DTREE_::dTree(const Sommet_dTree<DIMENSION> &A,const Sommet_dTree<DIMENSION> &B,dTree *mypere)
204 niveau=pere->niveau+1;
208 noeud_contenu=new vector<int>;
209 noeud_contenu->reserve((pere->noeud_contenu->size())/nbr_descendants);
211 for (i=0;i<DIMENSION;i++)
213 coord_max[i]=(*nuage)[0][i];
214 coord_min[i]=(*nuage)[0][i];
217 for (i=0;i<DIMENSION;i++)
233 cerr<<"DTREE : Construction de descendance sans père !"<<endl;
237 _TEMPLATE_ _DTREE_::dTree(const _DTREE_ &F)
239 // Pas Super Top Moumoute ... Recopie de pointeurs et pas de contenus, merdique
242 for (i=0;i<nbr_descendants;i++) descendant[i]=F.descendant[i];
243 noeud_contenu=F.noeud_contenu;
247 coord_max=F.coord_max;
248 coord_min=F.coord_min;
250 _TEMPLATE_ _DTREE_::~dTree()
255 // cout<<*this<<endl;
258 delete noeud_contenu;
260 for (i=0;i<nbr_descendants;i++) if (descendant[i])
262 delete descendant[i];
265 _TEMPLATE_ void _DTREE_::Get_Noeuds_Filtre(vector<int> &tmp)
270 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
272 for (i=0;i<nbr_descendants;i++) descendant[i]->Get_Noeuds_Filtre(tmp);
273 case DTREE_TERMINAL :
274 if (noeud_contenu->size()>0)
276 for (i=0;i<NBR_NOEUDS_PAR_CASE;i++) tmp.push_back((*noeud_contenu)[i]);
279 cerr<<"DTREE : dTree Non Valide dans Rempli_Noeuds_Filtre"<<endl;
283 _TEMPLATE_ Sommet_dTree<DIMENSION> _DTREE_::Get_Max() const
287 _TEMPLATE_ Sommet_dTree<DIMENSION> _DTREE_::Get_Min() const
291 _TEMPLATE_ int _DTREE_::is_in_dTree(NOEUD P) const
294 for (int i=0;i<DIMENSION;i++)
296 test=((coord_min[i]<=P[i])&&(P[i]<=coord_max[i]));
302 // calcule la distance Linf d'un point exterieur, si le point est interieur, le resultat sera erroné
303 _TEMPLATE_ double _DTREE_::calcule_distance(NOEUD P) const
305 double min=fabs(coord_max[0]-P[0]);
308 for (i=0;i<DIMENSION;i++)
310 tmp=fabs(coord_max[i]-P[i]);
311 if (min>tmp) min=tmp;
312 tmp=fabs(coord_min[i]-P[i]);
313 if (min>tmp) min=tmp;
317 _TEMPLATE_ _DTREE_ & _DTREE_::operator = (const _DTREE_ & F)
319 // Pas Super Top Moumoute ... Recopie de pointeurs et pas de contenus, merdique
322 for (i=0;i<nbr_descendans;i++) descendant[i]=F.descendant[i];
323 noeud_contenu=F.noeud_contenu;
327 coord_max=F.coord_max;
328 coord_min=F.coord_min;
330 // selecteur doit etre lu comme un nombre en base 2
331 // il specifie les coordonnes de reference du sommet dont on veut les coordonnees reelles
332 // ex : en dim 3 : coord_min=0+0*2+0*4=0
333 // coord_max=1+1*2*1*4=7
334 // donc les unites pour x, les 2aines pour y, les 4aines pour z
335 _TEMPLATE_ Sommet_dTree<DIMENSION> _DTREE_::donne_sommet(int selecteur) const
337 Sommet_dTree<DIMENSION> p;
339 int residu=selecteur;
341 for (i=0;i<DIMENSION;i++)
344 if (reste==0) p[i]=coord_min[i]; else p[i]=coord_max[i];
345 residu=(residu-reste)/2;
349 _TEMPLATE_ int _DTREE_::a_des_fils() const
351 if ((etat==DTREE_COURANT)||(etat=DTREE_RACINE)) return 1;
354 _TEMPLATE_ _DTREE_ * _DTREE_::trouve_dTree_contenant(NOEUD P) const
356 Sommet_dTree<DIMENSION> B(coord_min,coord_max);
358 if ((etat==DTREE_RACINE)&&(!is_in_dTree(P))) return NULL; // Le noeud est extérieur a l'dTree
359 else if (etat==DTREE_TERMINAL) return (dTree *) this; // Le noeud est dans *this qui est terminal
361 for (i=0;i<nbr_descendants;i++)
363 Sommet_dTree<DIMENSION> A=donne_sommet(i);
365 for (j=0;j<DIMENSION;j++)
367 if (A[j]<=B[j]) test=((A[j]<=P[j])&&(P[j]<=B[j]));
368 else test=((A[j]>=P[j])&&(P[j]>=B[j]));
372 if (test) return descendant[i]->trouve_dTree_contenant(P); // Propagation
375 // si de le dTree n'est pas TERMINAL, scanne tous les points du nuage du pere pour trouver le point le plus proche
376 // sinon scanne uniquement les points contenus dans le dTree
377 _TEMPLATE_ int _DTREE_::trouve_plus_proche_point_bourrin(NOEUD P) const
383 if (etat!=DTREE_TERMINAL)
385 min=DistanceL2(P,(*nuage)[0]);
386 for (i=1;i<nuage->size();i++)
388 tmp=DistanceL2(P,(*nuage)[i]);
399 if (noeud_contenu->size()!=0)
401 num_sol=(*noeud_contenu)[0];
402 min=DistanceL2(P,(*nuage)[num_sol]);
403 for (i=1;i<noeud_contenu->size();i++)
405 tmp=DistanceL2(P,(*nuage)[(*noeud_contenu)[i]]);
408 num_sol=(*noeud_contenu)[i];
414 else return DTREE_FANTOME;
417 // Fonction de pilotage de la recursion, lance les tests preliminaires et invoque la recursion
418 _TEMPLATE_ int _DTREE_::trouve_plus_proche_point(NOEUD P) const
421 if (etat==DTREE_FANTOME)
423 cerr<<"DTREE : TROUVE PLUS PROCHE POINT : L'octree n'est pas construit !"<<endl;
426 dTree *ref=trouve_dTree_contenant(P);
431 // Le point est intérieur
435 delta=DistanceInf((ref->pere->coord_max),(ref->pere->coord_min));
439 cerr<<"DTREE : TROUVE PLUS PROCHE POINT : l'octree contenant le noeud n'a pas de pere !"<<endl;
445 // Le point est extérieur
447 delta=DistanceL2(P,(*nuage)[0]);
453 retour=tppp_rec(P,delta,flag);
457 _TEMPLATE_ int _DTREE_::trouve_un_point() const
459 if (noeud_contenu!=NULL)
461 if (noeud_contenu->size()>0) return (*(noeud_contenu))[0];
466 for (i=0;i<nbr_descendants;i++)
468 return descendant[i]->trouve_un_point();
472 // partie recursive de trouve_plus_proche_point
473 // change le flag en 1 si un point plus proche est trouvé
474 // adapte automatiquement la distance delta de filtrage
475 _TEMPLATE_ int _DTREE_::tppp_rec(NOEUD P,double &delta,int &flag) const
477 if (Localise_Point(P,delta)==0)
480 // La distance du point à l'octree est plus grande que delta
482 return DTREE_FANTOME;
486 if (etat==DTREE_TERMINAL)
488 if (noeud_contenu->size()>0)
490 num_Ptmp=trouve_plus_proche_point_bourrin(P);
491 dtmp=DistanceL2((*nuage)[num_Ptmp],P);
495 // Le point le plus proche minimise delta, c'est un bon candidat, on l'envoie !
502 // Le point le plus proche ne minimise pas delta
504 // ===========> peut etre rajouter exit(-1); ??????????
505 return DTREE_FANTOME;
510 // L'octree qui contient P ne contient aucun point
512 return DTREE_FANTOME;
516 int num_sol=DTREE_FANTOME;
517 for (i=0;i<nbr_descendants;i++)
519 num_Ptmp=descendant[i]->tppp_rec(P,delta,flag);
520 if ((num_Ptmp!=DTREE_FANTOME)&&(flag==1))
523 // On a trouvé un point qui minimise delta dans une branche
528 // A ce stade, on a trouvé un point qui minimise tous les deltas, c'est donc le point le plus proche
529 // REMARQUE : cette affirmation est à la base de l'algorithme par dTree mais est loin d'étre évidente
534 // renvoie 1 si la distance L inf du noeud a l'octree est plus petite que d, 0 sinon
535 _TEMPLATE_ int _DTREE_::Localise_Point(NOEUD P,double d) const
538 for (i=0;i<DIMENSION;i++)
540 if (P[i]>coord_max[i]+d) return 0;
541 if (P[i]<coord_min[i]-d) return 0;
546 // Méthode de construction du dTree par propagation
547 _TEMPLATE_ void _DTREE_::cree_filiation()
549 if (etat!=DTREE_RACINE)
551 niveau=pere->niveau+1;
558 if (noeud_contenu->size()<=NBR_NOEUDS_PAR_CASE)
566 Sommet_dTree<DIMENSION> centre(coord_max,coord_min);
568 for (i=0;i<nbr_descendants;i++)
570 descendant[i]=new dTree(centre,donne_sommet(i),this);
573 for (num_loc=0;num_loc<noeud_contenu->size();num_loc++)
575 int indice=(*noeud_contenu)[num_loc];
576 NOEUD & courant=(*nuage)[indice];
578 for (i=0;(test)&&(i<nbr_descendants);i++)
580 if (descendant[i]->is_in_dTree(courant))
582 descendant[i]->noeud_contenu->push_back(indice);
588 delete noeud_contenu;
591 for (i=0;i<nbr_descendants;i++) descendant[i]->cree_filiation();
594 _TEMPLATE_ int _DTREE_::Get_Nbr_Descendants_Non_Vides() const
600 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
602 for (i=0;i<nbr_descendants;i++) ndnv+=descendant[i]->Get_Nbr_Descendants_Non_Vides();
604 case DTREE_TERMINAL :
605 if (noeud_contenu->size()>0) return 1;
608 cerr<<"dTree Non Valide dans Get_Nbr_Descendants_Non_Vides"<<endl;
612 _TEMPLATE_ int _DTREE_::Get_Nbr_Descendants_Vides() const
618 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
620 for (i=0;i<nbr_descendants;i++) ndnv+=descendant[i]->Get_Nbr_Descendants_Vides();
622 case DTREE_TERMINAL :
623 if (noeud_contenu->size()==0) return 1;
626 cerr<<"dTree Non Valide dans Get_Nbr_Descendants_Non_Vides"<<endl;
630 _TEMPLATE_ int _DTREE_::Get_Profondeur_Max() const
637 case DTREE_RACINE : int pourlefunlecompilolesttropcon;
639 for (i=0;i<nbr_descendants;i++)
641 tmp=descendant[i]->Get_Profondeur_Max();
642 if (tmp>prof) prof=tmp;
645 case DTREE_TERMINAL : return niveau;
647 cerr<<"dTree Non Valide dans Get_Nbr_Descendants_Non_Vides"<<endl;