Les metaballs 3d ou blobs
                         -------------------------

                                dake/calodox


 Introduction
 ------------

 Les metaballs 3d peuvent être génerés de différentes manières. Celle
 qui est le plus souvent utilisée fait intervenir l'algorithme
 des marching-cubes. Cette algorithme va créer le mesh des metaballs
 (la surface) à partir d'une grille à 3 dimensions et qui contient
 des valeurs. Les metaballs sont des isosurfaces, en gros des
 surfaces génerées à partir de valeurs contenues dans une grille.

 Calcul de la grille
 -------------------

 Les metaballs sont le plus souvent contenus dans un grand cube.
 Ce cube est lui-même subdivisé en de plus petits cubes. C'est
 la grille. Dans la plupart des cas, on utilise une grille de 32x32x32
 mais on peut très bien se contenter d'une grille de 16x16x16.
 Attention ! Augmenter la taille de la grille induit des temps
 de calcul qui n'augmentent pas linéairement.
 Chaque point de la grille contient une valeur (isovaleur),
 ces valeurs proviennent de formules. Dans le cas des metaballs,
 nous allons remplir la grille avec des "spheres", il faut donc
 utiliser la formule de la sphere avec quelques modifications.
 Chaque metaball possède un centre, nous allons calculer les
 isovaleurs pour chaque centre et les additionner. C'est un peu
 comme pour des metaballs en 2D sauf que nous travaillons en 3D.
 Le fait d'additionner les valeurs va créer un "pont" entre
 les spheres suivant leur distance.

 La formule pour des metas est :

                               1.0/(x^2+y^2+z^2)

                               avec :
                               x=(centre_blob.x - position_grille.x)
                               y=(centre_blob.y - position_grille.y)
                               z=(centre_blob.z - position_grille.z)

 =========================================================================
 Code (avec une grille de 32x32x32 et 4 metaballs) :

 #define GRILLE_X     32
 #define GRILLE_Y     32
 #define GRILLE_Z     32
 #define NB_METABALLS 4

 float   metagrille[GRILLE_X][GRILLE_Y][GRILLE_Z];      //defini la grille.

 for (int gz=0;gz<GRILLE_Z;gz++)                
      {
      for (int gy=0;gy<GRILLE_Y;gy++)
           {
           for (int gx=0;gx<GRILLE_X;gx++)       
                {
                 metagrille[gx][gy][gz]=0.0;            //efface la grille 
                   for (int n=0;n<NB_METABALLS;n++)
                       {
                       diffx=gx-metaball[n].centrex; 
                       diffy=gy-metaball[n].centrey;
                       diffz=gz-metaball[n].centrez;
                       distance=1.0/(diffx^2 + diffy^2 + diffz^2);
                       metagrille[gx][gy][gz]+=distance;
                       }
                }
           }
      }
 =========================================================================

 En gros pour chaque point de la grille : 
 1. on efface la grille en remettant l'isovaleur à 0
 2. pour chaque metaball :
                         - on calcule la distance depuis le centre du
                           metaball jusqu'à la position ou nous sommes
                           dans la grille
                         - on rajoute cette valeur à l'isovaleur
                           ou nous sommes dans la grille


 Les marching-cubes, étape I
 ---------------------------

 C'est ici qu'intervient l'algo des marching-cubes (AMC). Celui-ci va
 construire l'objet 3d à partir de cette grille. Le principe est assez simple,
 notre grille forme de petits cubes. Un cube contient 8 isovaleurs (pour
 chaque sommet), le cube est aussi formé de 12 arêtes. Les vertexs du mesh
 se trouveront _sur_ ces arêtes, pas forcément au centre de l'arête, mais
 ils ne se trouveront jamais à l'intérieur du cube. L'AMC va scanner
 chaque cube et va créer les vertexs/triangles correspondant.
 Il y'a un maximum de 5 triangles par cube. Il peut y'en avoir aucun aussi.

 Maintenant, nous allons devoir définir un "isolevel", il s'agit en fait
 d'une "limite", la frontière entre 2 parties de la grille. Prenons
 l'exemple d'une boule de glace qui flotterait dans l'air. L'air serait
 a 20°, la glace à sa surface serait de 5° et commencerait à fondre
 avec un centre à 0°. Notre isolevel serait à 5 car c'est justement
 la température présente à l'endroit le plus proche de la surface.
 Nous ne garderons que les endroits ou la température sera égale
 à 5°, ce qui nous donnera en gros, la surface de la boule et
 seulement la surface, non pas ce qui est à l'intérieur. C'est
 exactement la même chose avec les metaballs, nous allons diviser
 la grille en 2 parties, celle qui représente "l'intérieur" des metaballs,
 et l'autre qui est la partie qui les entoure. La frontière entre
 ces 2 parties formera notre objet en 3d.

 Nous allons à présent scanner chaque cube de la grille (ce qui
 ferait 31^3 pour une grille de 32x32x32, le cube n°32 n'existant
 pas). L'AMC est basé sur des tables de combinaisons qui vont permettre
 de recréer les triangles. Ces tables sont indexées grâce à un offset
 de 8 bits. Le cube se présente comme-ceci :

          4-----------5       numeros de sommets 
         /|          /|
       7'----------6' |
       |  |        |  |
       |  |        |  |
       |  0--------|--1
       | /         | /
       3'----------2'

                
          .-----4-----.       numeros d'arêtes
         7|          5|
       .'-----6----.' |
       |  8        |  9
      11  |       10  |
       |  .----0---|--.
       | 3         | 1
       .'----2-----''
             

 Pour chaque coin, nous allons voir si l'isovaleur est inférieur à
 l'isolevel. Si c'est le cas, nous mettons le bit n°X (X=le numero
 du coin) à 1. Cela nous donne un index de 8 bits, code correspondant :

 =========================================================================
 cindex=0;
 if (cube[0]<isolevel) cindex #=1;
 if (cube[1]<isolevel) cindex #=2;
 if (cube[2]<isolevel) cindex #=4;
 if (cube[3]<isolevel) cindex #=8;
 if (cube[4]<isolevel) cindex #=16;
 if (cube[5]<isolevel) cindex #=32;
 if (cube[6]<isolevel) cindex #=64;
 if (cube[7]<isolevel) cindex #=128;
 =========================================================================

 Si cindex vaut 0 ou 255 après avoir effectué les opérations ci-dessus,
 cela veut dire que nous sommes en présence d'un cube qui n'a pas
 d'intersection(s) avec la surface, nous pouvons passer au cube suivant.
 Dans le cas contraire, nous allons utiliser cindex pour adresser
 dans la table des arêtes. Cette table possède 256 valeurs de 12 bits.
 On les stocke en 16 bits pour faciliter et accélérer les calculs.
 Elle va permettre de savoir quelles arêtes du cube sont touchées
 par les metaballs. Chaque bit de ce code de 12bits représente
 un index. Si le bit n°X est à 1, cela veut dire que l'arête
 n°X est coupée par l'isosurface.

 Les marching-cubes, étape II
 ----------------------------

 Maintenant calculons la position des vertexs sur les arêtes
 touchées par l'isosurface. Ce calcul se fait par le biais d'une
 interpolation linéaire. Voici 1 arête et les 2 isovaleurs correspondantes.
 Si nous fixons l'isolevel à 15, nous aurons notre vertex en plein centre :

          isolevel=15
 10 ----------*---------- 20

 Si nous fixons l'isolevel à 12, le vertex se trouvera plus vers la gauche :

   isolevel=12
 10 ---*----------------- 20

 Etant donnée que chaque sommet du cube possède des coordonnées (X,Y,Z),
 nous allons pouvoir trouver la position (X,Y,Z) du nouveau vertex !
 Il est conseillé de travailler en float pour les coordonnées 3D, afin
 de garder un maximum de précision. Si l'on veut rester en fixed point,
 il suffit de multiplier la position X,Y,Z de chaque sommet du cube
 par un facteur assez grand (10,20,32...).

 La formule consacrée pour l'interpolation est :

 NP=POS1+(POS2-POS1)*(isolevel-VAL1)/(VAL2-VAL1);

 NP=nouvelle coordonnée   POS1=coordonnée 1er point POS2=coordonnée 2eme point
 VAL1=isovaleur 1er point VAL2=isovaleur 2eme point

 On va faire ca pour chaque coordonnée (X Y et Z) :

 =========================================================================
 float temp=(isolevel-val1)/(val2-val1)  //attention au div by 0 !
                                         //faire un petit test au préalable
                                         //dans ce cas, retourner un vertex
                                         //qui se trouve au milieu de l'arête

 vertex.x=p1.x + temp * (p2.x - p1.x);
 vertex.y=p1.y + temp * (p2.y - p1.y);
 vertex.z=p1.z + temp * (p2.z - p1.z);
 =========================================================================

 Maintenant que nous avons calculés tous les vertexs des arêtes
 touchées par la surface. Nous allons créer les triangles. C'est
 extrêmement simple, il s'agit d'un simple accès à une table qui
 nous donne les numéros des vertexs pour chaque triangle. Cette
 table est une liste de 256 combinaisons. Chaque combinaison
 fait 16 _bytes_. Nous avons juste à décaler notre index de
 4 (mul par 16) : cindex<<=4; et nous pouvons accéder aux listes
 de triangles correspondant à notre cas de cube. Une combin. se
 présente comme ceci :

 {3,11,2,4,6,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}

 Le cas précedent veut dire :
 1 triangle constitué des vertexs se trouvant sur les arêtes 3, 11
 et 2 et 1 triangle avec pour vertexs 4,6 et 2. Il suffit de scanner
 les bytes jusqu'à ce qu'on arrive à -1 qui indique qu'il n'y pas
 d'autres triangles présents pour ce cube.

 Et voila ! Si on fait ca pour chaque cube, on obtient un mesh
 que l'on peut afficher comme n'importe quel autre objet 3d.

 Optimisations
 -------------

 Calculer tous les cubes de la grille peut s'avérer couteux
 en temps machine. Il y'a plusieurs moyens pour accélerer le tout,
 du moins, pour les metaballs calculés avec la formule citée plus
 haut dans cette doc.

 La première consiste à ne pas choisir un isolevel trop grand, en
 effet, plus cette valeur est grande, plus l'isosurface possedera
 de triangles.

 La deuxième est qu'il faut éviter les calculs inutiles, c-a-d
 de calculer des cubes qui ne sont pas touchés par la surface. Pour
 cela, on peut recourir à 2 moyens. La récursivité ou les boites
 englobantes.

 La récursivité
 --------------

 Cette méthode qui m'a été proposée par Mrz/Astroidea, consiste à ne
 scanner que les cubes se trouvant à la surface des metaballs.
 On commence au centre d'un des metaball, on continue dans n'importe
 quelle direction jusqu'à ce qu'on touche la surface (on peut savoir cela
 en regardant la valeur de l'index binaire).
 Ensuite, on commence l'algo de recursivité. On check le cube
 du haut, du bas, de gauche, de droite, si l'un d'entre eux est
 touché par la surface, on va dans sa direction et on refait les
 comparaisons. C'est une sorte de flood-filling mais en 3d. Si
 on se trouve dans une impasse, on revient en arrière jusqu'à ce
 qu'on puisse partir dans une autre direction. Il faut faire
 ca pour chaque metaball. Ne pas oublier de mettre des flags pour
 chaque cube, si le flag d'un cube est déjà mis à 1, on n'a pas besoin
 de le calculer. (ne pas oublier aussi d'effacer la table des flags à
 chaque mouvement des metaballs).
 On peut coder cela en sauvant les coordonnées dans la pile ou en utilisant
 des trees (en insérant des branches pour chaque nouveau cube et en
 remontant dans l'arbre dès qu'on ne peut plus avancer..)

 Les boîtes englobantes
 ----------------------

 Plus simple à coder, il suffit juste de créer une boîte qui entoure
 chaque metaball. Ensuite, on ne scanne que les cubes se trouvant
 _dedans_ la boîte. La aussi, il faut veiller à mettre des flags
 si un cube a déjà été calculé. Les boîtes permettent aussi d'accélerer
 le calcul des isovaleurs de la grille. On ne calcule que les isovaleurs
 se trouvant dans chaque boîte lors du remplissage de la grille. On peut
 réduire le nombre de cubes à calculer par un facteur de 2 jusqu'à 4.
 Pour une grille de 24^3, avec une boite englobante de 12x12x12, on arrive à
 environ 5000 cubes à calculer pour 4-5 métas (cela depend de leur
 position), contre 12167 cubes (23^3).

 Les normales
 ------------

 P. Doom/Carrots a écrit un article sur les normales pour les metaballs.
 Disponible dans le diskmag Pain 08.98 à www.space.ch/scene/pain
 Ces "normales" sont suffisantes pour un calcul de luminosité mais
 s'avèrent trop imprécises pour une élimination de faces cachées.

 Autres isosurfaces
 ------------------

 On peut facilement créer d'autres types d'objets en modifiant la
 formule qui remplit la grille ou avec d'autres trucs. L'un d'eux
 permet de faire des "gouttes", des metaballs qui laissent des
 trainées derrière eux. Il s'agit juste d'un "motion blur" de la grille.
 Au lieu de la remettre à 0 à chaque frame, on divise chaque isovaleur
 par un facteur + ou - grand. Cela va laisser une trace derrière
 chaque meta. (cf. yadot/byterapers et kill the clone/fudge)

 Morphing
 --------

 Très simples à mettre en oeuvre, on a juste à interpoler les isovaleurs
 de 2 grilles différentes et à appliquer les marching-cubes sur cette
 grille interpolée. (utilisé dans Stash/TBL, lorsque une boule
 créee des sortes de branches de chaque côté). 

 Conclusion
 ----------

 J'espère que cette doc aura été utile :) En ce qui concerne les
 tables nécessaires aux marching-cubes, elles sont disponibles
 dans le zip, avec la doc en anglais de Paul Bourke (qui contient
 aussi un source en C).

 -dake/calodox-
 -www.space.ch/scene/calodox-