Octrees optimisés grâce au code de Morton

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.

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 :

• avecla profondeur etl'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 ;
• avec la largeur du cube. L'exposant est lié à la taille des cubes qui lui sont associés.

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.

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.

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;
}

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.

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.

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

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.

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)₁₀ = (0101,1001,0001)₂

Le code de Morton de cette position est donné par : 010001000111₂ = 1095₁₀

D'après ce que nous avons vu précédemment, par exemple pour récupérer la feuille à la position (5,9,1)₁₀ dont le code de Morton est« 010001000111₂ », 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 :

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;
}

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é de, ce qui correspond au volume du cube courant.

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

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» */

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 :

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).

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.

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
}