Octrees optimisés grâce au code de Morton

Les octrees (arbres octaires) sont des structures de données utilisées dans les applications 3D pour accélérer une série d'opérations dans l'espace, comme la recherche des voisins d'un point ou la détection de collisions.

Une primitive pour bon nombre de ces applications est la recherche dichotomique. Elle peut être accélérée en employant une représentation particulière des coordonnées, le code de Morton. Il peut également servir pour itérer rapidement sur les éléments de l'arbre.

16 commentaires Donner une note à l'article (5)

Article lu   fois.

L'auteur

Liens sociaux

Viadeo Twitter Facebook Share on Google+   

I. Introduction

Les octrees (arbres octaires) sont des arbres de recherches dichotomiques tridimensionnels travaillant en partitionnant l'espace. Ces arbres sont indispensables pour travailler en 3D, autant dans le calcul scientifique et l'ingénierie que dans les jeux vidéo. Sans eux, la complexité (performance) des programmes devient vite ingérable.

Image non disponible
Un octree généré et coloré aléatoirement.

Image non disponible
Vue en coupe d'une sphère approchée par octree.

Un octree est une structure de données de type arbre dans laquelle chaque nœud compte soit zéro, soit huit fils. Chaque nœud occupe l'espace d'un cube et peut donc être subdivisé en huit parties égales.

Image non disponible
Illustration: un octree et son arborescence

Le code de Morton est une façon de représenter les coordonnées tridimensionnelles en une seule coordonnée. L'avantage de ce code est d'offrir un ordre sur un volume. Il devient alors beaucoup plus simple de travailler sur un octree.

Nous allons aussi créer une classe qui enregistrera le chemin (la liste des pères) du cube précédemment recherché afin d'améliorer la recherche suivante. Avec cette modification, en moyenne, la complexité de la recherche deviendra constante ; sans elle, elle est forcément linéaire en fonction de la profondeur.

II. Notion d'arbre octaire

II-A. Utilisations courantes des octrees

Les octrees sont très utilisés, car ils permettent une recherche de voisins efficace et une subdivision adaptative. Il en résulte un grand nombre d'applications :

  • la détection de collision en trois dimensions dans les logiciels de CAO ou les moteurs physiques ;
  • le sparse voxel octree pour réaliser du ray-casting ou ray-tracing ;
  • l'élimination des objets hors du cône de vue dans le cadre d'un rendu 3D ;
  • les maillages, d'une manière générale, notamment dans le cadre de la méthode des éléments finis ;
  • le problème à n corps, simulation de Barnes-Hut utilisée en astrophysique pour simuler les interactions gravitationnelles entre étoiles au sein d'une galaxie.

II-B. Vocabulaire

Avant de commencer, établissons le vocabulaire utilisé. Ici, « arbre » peut être remplacé par « octree » :

  • un nœud est un élément de l'arbre. Dans un octree, un nœud correspond à un cube dans l'espace 3D à partitionner ;
  • la racine d'un arbre est l'unique nœud ne possédant pas de parent, il permet d'accéder à tous les éléments de l'arbre ;
  • un nœud interne est un élément qui a des fils (huit dans notre cas) ;
  • une feuille de l'arbre est un nœud qui n'a pas de fils ;
  • la profondeur d'un nœud est la distance c'est-à-dire le nombre d'arêtes, de la racine au nœud ;
  • l'étiquette d'un nœud est l'information stockée dans le nœud.

Dans notre cas, nous allons introduire un autre terme au vocabulaire : l'exposant. Il s'agit du logarithmique de la taille du cube représenté par le nœud. Il donne donc à la fois des informations sur la profondeur et sur la taille du cube. Surtout, il donne l'emplacement des bits qui contiennent l'information que nous utilisons pour la recherche et l'insertion.

Il vérifie les propriétés suivantes :

  • Image non disponible avecImage non disponiblela profondeur etImage non disponiblel'exposant du cube racine. L'exposant décroit avec la profondeur de telle sorte que l'exposant est maximum à la racine et nul aux feuilles les plus basses ;
  • Image non disponibleavec Image non disponible la largeur du cube. L'exposant est lié à la taille des cubes qui lui sont associés.

II-C. Recherche et insertion

La première chose que nous voulons pouvoir faire sur un octree est d'ajouter des éléments à une certaine position dans l'espace et profondeur dans l'arbre et accéder à ces éléments. L'insertion est équivalente à l'accès, puisqu'insérer revient à rechercher l'emplacement où nous allons insérer l'élément. La complexité d'insertion de recherche d'un élément sera logarithmique en fonction du nombre d'éléments.

Pour cela, il faut choisir récursivement le bon fils parmi les huit au sein des cubes. Une première solution serait de regarder si le point recherché est en haut ou en bas, à gauche ou à droite, derrière ou devant le centre du cube, puis de choisir le fils correspondant. Pour cela, nous assemblons le résultat des trois tests qui forment trois bits en un nombre décimal entre 0 et 7 qui nous permet d'indexer les huit fils. Pour réaliser l'opération, il faut cependant connaître les coordonnées du centre du cube.

Voici un code permettant de choisir le bon fils d'un nœud avec cette première solution. La structure Octree doit posséder un membre fils qui est un pointeur sur un tableau de huit Octree ainsi que trois autres membres x, y et z qui sont les coordonnées du centre du cube qui représente le nœud. Ces coordonnées doivent préalablement être rentrées dans le nœud.

Choix du fils d'un nœud
Sélectionnez
Octree &oct2 = oct1.fils[ (posx > oct1.x) + ((posy > oct1.y) << 1) + ((posz > oct1.z) << 2) ];
// oct2 est une référence sur le fils de oct1 dans lequel se trouve le point de coordonnée (posx, posy, posz)

Cependant, il existe un moyen plus simple et équivalent de choisir le bon fils sans connaître les coordonnées du centre du cube : il suffit d'utiliser la décomposition binaire des coordonnées des points que nous voulons insérer. En effet, pour savoir si, dans un cube de largeur 2^k, un point est du côté droit ou gauche du cube, nous regardons le bit numéro k-1 de la coordonnée en x du point : si celui-ci est à 1, c'est qu'il est du côté droit du cube, sinon, du côté gauche. Il en va de même pour la hauteur et la profondeur. En assemblant les trois bits des coordonnées d'un même niveau (un pour chaque coordonnée), nous formons aussi un nombre entre 0 et 7 que nous utilisons pour indexer les fils d'un nœud.

Finalement, l'accès à un élément à la position :

(x,y,z) = (5,9,1)10 = (0101,1001,0001)2

serait réalisé par le code :

oct.fils[010].fils[001].fils[000].fils[111] ;

Plus besoin d'enregistrer le centre des cubes et de gérer les huit cas possibles séparément, il suffit d'entrelacer les bits des coordonnées du point que nous voulons insérer et le tour est joué.

Voilà à quoi ressemble une fonction de recherche au sein de l'octree. Nous avons besoin de la position cartésienne entière, de l'exposant du cube que nous voulons trouver et, évidemment, de la racine de l'octree. Nous voyons aussi les opérations bit à bit permettant d'entrelacer les bits des coordonnées pour obtenir l'indice du fils. La structure Octree doit posséder un membre fils qui est un pointeur sur un tableau de huit Octree. Si la fonction ne trouve pas l'octree, elle retourne la feuille à laquelle elle s'est arrêtée.

Recherche d'un élément
Sélectionnez
Octree* getcube(Octree* oct, int posx, int posy, int posz, int e){
    int e0=30; // Conventionnellement, dans ce code, la largeur du cube racine est de 2^30
    while((*oct).fils!=NULL && e0 > e){// tant qu'il y a des fils et que la profondeur n'est pas atteinte
        e0--;// on descend d'un niveau dans la profondeur de l'arbre
        oct=&(*oct).fils[ ((posx >> e0) & 1) + (((posy >> e0) & 1) << 1) +  (((posz >> e0) & 1) << 2) ];
    }
    return oct;
}

Image non disponible
Quadtree avec les points A(50,41) et B(46,22).

Cette manière de faire est bien correcte. En effet, la position du cube que nous choisissons est égale aux déplacements successifs générés par le choix des fils, et ces déplacements sont les puissances de 2. Comme les coordonnées du cube sont l'arrondi des coordonnées du point que nous recherchons, nous retrouvons ces puissances de deux dans les bits de coordonnées du point. Ainsi, le choix des fils lors d'une recherche dichotomique au sein d'un octree est directement déterminé par les bits des coordonnées du point.

II-D. Suppression et simplification

Il faut aussi pouvoir supprimer des éléments. Par exemple, lorsqu'un corps se déplace et change de « cube » au sein de l'octree, il faut supprimer l'élément du cube s'il n'y est plus et le replacer dans le cube qui lui correspond. Pour effectuer cette étape efficacement, il vaut mieux que les étiquettes des nœuds de l'octree soient des pointeurs ou de petits objets, pour éviter les allocations importantes de mémoire.

Une fois les éléments supprimés, il faut simplifier l'octree en supprimant les branches vides. Cette étape s'effectue de manière récursive : pour chaque nœud qui a huit feuilles vides, supprimer le nœud et tester le père.

II-E. Recherche de voisins et zones

La recherche de voisins est très utile pour détecter des collisions ou pour réaliser des interactions qui ont lieu à courte portée.

Il y a deux cas possibles :

  • soit nous cherchons les cubes voisins à un cube, il y en a alors 26 ;
  • soit nous cherchons une zone autour d'un point contenant au moins tous les éléments situés dans un certain rayon autour du point. Nous pouvons définir une telle zone avec huit cubes seulement.
Image non disponible


Les cubes autour du cube central, il y a donc 27 cubes au total (en comptant le cube central).

Image non disponible


Le point contenu dans les huit cubes est à une distance inférieure à la moitié d'un gros cube des bords de l'ensemble.

Si nous réalisons la détection de collision ou le calcul des interactions sur tout l'octree à la fois, nous pouvons ne prendre que la moitié des 26 cubes entourant le cube (soit 13 cubes). L'autre moitié sera traitée par les voisins.

Pour trouver les huit cubes qui contiennent la sphère de rayon r autour du point p, nous cherchons les cubes de côté supérieur ou égal à 2r qui contiennent les points p+(-r,-r,-r) ; p+(+r,-r,-r) ; p+(-r,+r,-r)…

En fin de compte, nous avons avec les huit cubes une zone plus grossière qu'avec les 27 cubes.

La zone de huit cubes peut aussi stocker des données qui ont une portée sur une certaine zone : nous pouvons ainsi stocker dans les nœuds internes toutes les variables qui ont un intérêt pour les fils du nœud.

III. Itérateur à octree

Nous avons créé un format de données et nous voulons pouvoir le parcourir facilement. Nous allons donc créer un itérateur. Pour cela, nous allons utiliser le code de Morton. Son avantage est d'offrir un ordre sur un volume qui est adapté aux octrees. Cet ordre correspond au parcours récursif d'un octree. Ce code est un outil pratique pour transformer les fonctions récursives sur l'octree en fonctions itératives.

III-A. Le code de Morton

Le code de Morton est réalisé en entrelaçant les bits de coordonnées en un seul nombre :

(x,y,z) = (5,9,1)10 = (0101,1001,0001)2

Le code de Morton de cette position est donné par : 0100010001112 = 109510

D'après ce que nous avons vu précédemment, par exemple pour récupérer la feuille à la position (5,9,1)10 dont le code de Morton est« 0100010001112 », il suffit de regrouper les chiffres de ce dernier trois par trois et de les utiliser pour indexer les fils :

oct.fils[010].fils[001].fils[000].fils[111] ;

Ce parcours de l'arbre peut se mettre très rapidement sous la forme d'un algorithme pour récupérer un nœud de l'octree à partir de son code de Morton :

Accès d'un élément
Sélectionnez
Octree* getcube(Octree* oct, int mcode, int e){
    int e0=30; // Conventionnellement, dans ce code, la largeur du cube racine est de 2^30
    while((*oct).fils!=NULL && e0 > e){// tant qu'il y a des fils et que la profondeur n'est pas atteinte
        e0--;// on descend d'un niveau dans la profondeur de l'arbre
        oct=&(*oct).fils[ (mcode >> 3*e0) & 7 ];
    }
    return oct;
}

En pratique, nous prenons la coordonnée x en poids faible et z en poids fort, car, même si, sur le papier, l'ordre est inversé (1,1,0) → (011), au niveau des indices, c'est le bon ordre :

  • lorsque nous stockons les variables sous forme de vecteur, la composante x est dans la première case (indice 0) et le z est dans la troisième case (indice 2) ;
  • lorsque nous entrelaçons les coordonnées x,y,z dans un même nombre, x est de poids faible, le masque correspondant est (1<<0), z est de poids fort, le masque correspondant est (1<<2).
Image non disponible
Illustration: disposition des fils au sein d'un nœud

III-B. Transformation de la récursion à l'itération

Qu'est-ce que nous effectuons lorsque nous parcourons l'octree de façon récursive ? Pour chaque cube, nous parcourons récursivement les sous-cubes dans l'ordre dans lequel nous les avons rangés dans le tableau des fils (0, 1, 2… 7).

Analysons l'évolution du code de Morton du cube courant pendant le parcours des feuilles.

Lorsque nous venons de parcourir une feuille et que nous passons à la suivante, les trois bits de position du code de Morton relatif au cube courant sont incrémentés au niveau du fils en question. Le code de Morton est donc incrémenté de la puissance de 8 correspondant à la profondeur.

oct.fils[001].fils[110] ; → oct.fils[000].fils[111] ; évolution du code de Morton : (14→ 15)

Si la feuille était la dernière du tableau de fils, nous changeons de père et nous prenons le premier fils du père. Cette opération correspond à mettre à 0 la suite de 1 qui formait la position de la feuille et de mettre un bit de poids plus fort à 1 qui correspond au nouveau père. C'est comme si l'on avait additionné la puissance de 8 correspondant à la profondeur de la dernière feuille au code de Morton : l'addition a généré une retenue qui s'est propagée, en mettant à zéro la suite de 1 jusqu'à atteindre le premier bit à 0 pour le mettre à 1.

oct.fils[001].fils[111] ; → oct.fils[010].fils[000] ; évolution du code de Morton : (15 → 16)

Il suffit donc d'additionner au code de Morton la puissance de 8 correspondant à la profondeur de la feuille actuelle pour avoir le code de Morton de la prochaine feuille.

Sur l'animation ci-dessous, un octree est parcouru en utilisant le code de Morton. À chaque itération, le code de Morton est incrémenté deImage non disponible, ce qui correspond au volume du cube courant.

Image non disponible

Ainsi, une boucle permettant de parcourir le cube sans fonction récursive serait, en utilisant le code de Morton :

Parcours récursif d'un octree
Sélectionnez
int mcode; // On y range le code de Morton
int e = 30; // On y récupérera le niveau des feuilles
for(mcode = 0; mcode < CODE_MAX; mcode += 1 << 3*e){// Addition de la puissance de 8 correspondant à la profondeur
    Octree &oct2 = *getlastcube(oct, mcode, &e);
    // Actions sur la références oct2...
}
/* Ici getlastcube() va chercher le dernier cube de l'octree oct
à la position de code de Morton mcode et donne son exposant dans «e» */

III-C. Réalisation de l'opération d'incrémentation en coordonnées cartésiennes

Le code de Morton oblige à utiliser des valeurs sur 64 bits pour avoir une profondeur maximale de 21, car il faut placer les trois coordonnées cartésiennes dans un seul entier. Il faut aussi convertir le code de Morton en coordonnées cartésiennes et inversement pour pouvoir l'utiliser. Ce n'est pas très pratique. Il y a cependant des méthodes efficaces pour calculer le code de Morton.

Nous allons donc exprimer l'opération en coordonnées cartésiennes qui correspond à l'incrémentation du code de Morton pour travailler exclusivement en coordonnées cartésiennes. L'avantage sera donc l'utilisation d'une seule représentation (donc pas besoin de conversion) et la possibilité de travailler sur des profondeurs allant de 30 bits (en 32 bits en gardant un bit de dépassement et le bit de signe) à 62 bits (en 64 bits) tout en conservant — voire en améliorant — la rapidité de l'algorithme.

L'incrémentation correspond à l'addition d'une puissance de 8 au code de Morton. Cela revient à additionner une puissance de 2 à la coordonnée de poids faible tout en propageant les retenues sur toutes les coordonnées à la fois. D'abord, nous propageons les retenues 3 par 3, puis nous gérons les bits de la profondeur à laquelle la retenue s'arrête. Voilà une solution :

Addition d'une puissance de 2
Sélectionnez
int dif = posx & posy & posz; 
dif = ((dif + (1 << e)) ^ dif); // Tous les bits qui seront modifiés par l'addition de la puissance e.
int p = (dif + (1 << e) ) >> 1; // Puissance d'arrêt de retenue        
int i = (posx & p) | ((posy & p) << 1) | ((posz & p) << 2);
i += p; // Nouveau code pour la profondeur d'arrêt de retenue
posx &= ~dif; posx |= i & p; // Mise à jour des nouvelles positions
posy &= ~dif; posy |= (i >> 1) & p;
posz &= ~dif; posz |= (i >> 2) & p;

Cela fait beaucoup d'opérations, mais l'avantage est de n'avoir aucun branchements (pas de boucles ni de conditions).

IV. Optimisation de la recherche de proximité

Pour améliorer la recherche du prochain élément en tenant compte de la liste des pères du cube précédemment recherché, nous devons calculer la profondeur du père commun le plus proche entre le cube précédent et le cube que l'on veut trouver. Puis, il suffira de partir de ce père pour trouver le nouveau cube, au lieu de partir de la racine. Cette opération va nous permettre d'accélérer les recherches de proximité.

Nous allons trouver la profondeur à partir de laquelle les coordonnées sont égales. Pour cela, nous allons réaliser une différence bit à bit des coordonnées grâce à un « ou exclusif » (qui correspond, en C++, à l'opérateur ^).

Finalement, nous regardons à partir de quelle profondeur, lorsque l'on décale vers la gauche les bits (pour supprimer les bits de poids faible) de notre entier de différence, celui-ci est nul (ce qui correspond à aucune différence de coordonnées).

calcul du niveau commun
Sélectionnez
unsigned int difference = (posx1 ^ posx2) | (posy1 ^ posy2) | (posz2 ^ posz2); 
while(difference >> e != 0) { e++; } // rend l'exposant commun dans «e»
// en supposant que «e» était d'exposant inférieur à l'exposant commun entre les deux positions 1 et 2

IV-A. Calcul rapide du logarithme

L'opération précédente revient en fait à calculer la partie entière d'un logarithme en base 2. Nous pouvons ainsi remplacer la boucle par une conversion int → float de la différence. En effet, en transformant l'entier en nombre réel à virgule flottante, le processeur réalise au passage un logarithme en base 2 pour placer l'exposant du nombre dans sa forme flottante normalisée (voir norme IEEE 754).

Il ne reste plus qu'à récupérer les bits de l'exposant pour avoir la profondeur commune des pères. Cependant, puisque nous utiliserons souvent le traceur pour aller chercher des cubes proches, il ne suffira que d'une itération de la boucle pour trouver le père commun. Ainsi, la complexité moyenne en utilisant la boucle ou la conversion est à peu près équivalente.

Voilà une fonction qui renvoie la partie entière du logarithme en base 2 d'un nombre.

extraction de l'exposant
Sélectionnez
union int_float{//pour faire des opérations bit à bit sur des float
    int i;
    float f;
};
inline int getexposant(float valeur){
    union int_float dif;
    dif.f = valeur; // On rentre le nombre sous forme d'un flottant dans un type union
    return ((dif.i >> 23) - 127); // On opère sur le flottant comme un entier
}

V. Implémentation du traceur à octree

Maintenant, nous pouvons créer une classe pour itérer sur les octrees et aller chercher des octrees de proximité avec une complexité moyenne constante (et au pire logarithmique).

Cette classe enregistrera la liste des pères du cube courant. Puisque l'on travaillera avec une profondeur de 30 nœuds maximum, la classe aura un tableau de pointeurs de 32 octrees afin de stocker les pères, la racine et un pointeur sur la racine. Ainsi, à tout moment, nous pourrons avoir accès aux pères des octrees que nous considérons, ce qui est pratique lorsque nous voulons stocker des données qui ont des portées variables et y accéder directement.

Nous intégrerons aussi à notre traceur une échelle et un vecteur de translation initial, pour pouvoir utiliser des coordonnées réelles (et non simplement entières) et ne pas être contraint par la taille du cube racine qui doit être une puissance de deux. Évidemment, nous ne serons pas obligé d'utiliser l'échelle et pourrons utiliser les coordonnées cartésiennes entières.

classe Traceur
Sélectionnez
template<class T> 
class Traceur{ 
private: 
    int posxyz[3]; // les coordonnées cartésiennes
    int e; // l'exposant du cube courant
    Vect3f translation; 
    float scale,_scale;// échelle et son inverse
    Octree<T>* tab[E_MAX+2]; // le tableau de pères (32 pointeurs)
public:
   // les fonctions associées :
   // incrémentation de la position, recherche d'élément, ajout d'élément
};

Elle travaillera sur une classe Octree toute simple :

classe Octree
Sélectionnez
template<class T1> 
class Octree{ 
public: 
    Octree<T1> *ptr; 
    T1 node;
    
    // les fonctions relatives à l'octree :
    // constructeurs, destructeur, opérateur = de copie récursive...
};

VI. Les sources

octree.h : implémentation des classes Octree et Traceur.

La classe Traceur a des fonctions pour rechercher un cube à une certaine position, pour insérer et créer un cube à une certaine position et profondeur. Ces fonctions peuvent être appelées avec des coordonnées entières ou flottantes. Elles sont aussi disponibles en mode protégé avec un test d'appartenance à l'octree. Une autre fonction parcourt tous les éléments ou toutes les feuilles de l'octree. Il y a des fonctions de sélection et de création de zone.

vect3f.h : vecteur de trois nombres à virgule flottante.

test.cpp : exemple d'utilisation.

Il n'y a aucune dépendance.

sources C++

VII. Conclusion

Nous avons réalisé un itérateur à octree qui est très pratique à utiliser et qui offre de très bonnes performances. De plus, cet itérateur à octree enregistre aussi le chemin de l'octree, pour s'en servir et ainsi accélérer les prochaines recherches.

VIII. Remerciements

Merci à dourouc05 pour ses relectures complètes et rigoureuses.

Merci aussi à Littlewhite pour ses relectures.

Merci à f-leb pour ses corrections.

Vous avez aimé ce tutoriel ? Alors partagez-le en cliquant sur les boutons suivants : Viadeo Twitter Facebook Share on Google+   

  

Les sources présentées sur cette page sont libres de droits et vous pouvez les utiliser à votre convenance. Par contre, la page de présentation constitue une œuvre intellectuelle protégée par les droits d'auteur. Copyright © 2016 Pierre Bénet. Aucune reproduction, même partielle, ne peut être faite de ce site ni de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.