Sphère procédurale

Comme nous l'avons vu, nous savons générer des surfaces procédurales. Ces surfaces sont raisonnablement convaincantes. Nous pourrions encore les améliorer, mais il me semble que nous pouvons pour le moment essayer de les utiliser en l'état pour en faire quelque chose, en 3D par exemple !

J'ouvre une parenthèse pour vous signaler que maintenant que nous commençons à faire de la 3D, autant que possible je m'efforcerai de ne pas rentrer dans les détails techniques d'OpenGL pour me concentrer sur l'aspect théorique de la création et des algorithmes. Vous êtes cependant fortement encouragé à apprendre à manipuler cette API si ce n'est pas déjà le cas.

Les terrains procéduraux que nous avons déjà faits jusqu'à présent pourraient être mis bout à bout sur une surface plane infinie. On aurait un paysage, et on pourrait le peupler. Mais ce ne serait pas très intéressant pour notre projet. Ce que nous voudrions plutôt, c'est par exemple créer une planète procédurale, ce qui est très différent techniquement.

Dans le tout premier article à propos de génération procédurale où l'on parlait du bruit de Perlin j'avais rapidement expliqué ce qu'était le displacement mapping. À titre de rappel, c'est le procédé par lequel on fait s'élever un vertex dans une direction donné en fonction de l'intensité du pixel de la heightmap qui lui est associé. Sur une surface plane, c'est très facile. On surélève le vertex selon l'axe z, par exemple.

displacement.png

Dans le cas d'une sphère c'est un petit peu plus compliqué parce qu'on ne travaille pas dans le même système de coordonnées. Tout comme dans un repère cartésien un point a trois coordonnées qui cette fois ne forment plus le triplet (x,y,z) mais (θ,Φ,r) dont chaque symbole Théta majuscule, Phi majusctule et r sont respectivement appelés latitude, longitude et altitude.

La première chose que l'on comprend intuitivement, c'est que l'on peut faire correspondre les trois coordonnées cartésiennes avec leurs équivalents sphériques, en particulier, l'axe z deviens l'altitude r. Ce qui veut dire qu'une fois que nos surfaces procédurales serons projetées sur la sphère, c'est l'altitude r qu'on modifiera en fonction de l'intensité du pixel local.

Mais maintenant se pose un plus gros problème. Un problème majeur même, millénaire pourrait on dire, parce que c'est une affaire de cartographes. Il s'agit de savoir comment projeter une surface plane contenant des informations topologiques sur une sphère. Ce problème n'est pas trivial du tout et il existe de très nombreuses façons de faire selon le contexte dans lequel on travaille. Il va donc falloir poser au moins trois contraintes.

  1. Faible cout en terme de calculs au moment de la génération de la sphère et de la projection de la surface sur celle-ci.
  2. Faible distorsions.
  3. La sphère et la surface projetée doivent être simple à manipuler à postériori.

Je ne vais pas rentrer dans le détail, mais avec ça on peut déjà exclure pas mal d'approches. On va tout de même prendre le temps d'expliquer pourquoi on ne va pas utiliser la plus intuitive d'entre elles: La projection équirectangulaire.

Commençons par voir comment ça marche.

La projection équirectangulaire consiste à plaquer sur une sphère un rectangle de dimensions n en hauteur et 2n en largeur. Pour ce faire on fait correspondre comme on l'a dit les coordonnées x et y aux coordonnées sphériques (θ,Φ). Attention cependant, la latitude θ varie entre -90° et 90° alors que la longitude Φ varie entre -180° et 180° (ou 0 et 360°, selon la convention). Autrement dit il faudra faire varier linéairement les coordonnées sphériques en fonction des coordonnées cartésiennes à l'aide d'un produit en croix.

Voilà pour l'essentiel. Si on procède ainsi cependant la projection produira de très fortes distorsions topologiques et ça ne sera pas joli du tout. Ce qui signifie qu'il faut traiter en amont la heightmap rectangulaire avant de la plaquer sur la sphère. Ce traitement va consister à appliquer des distorsions sur la heightmap qui seront compensées une fois plaquée sur la sphère. Pour vous aider à visualiser le phénomène il suffit d'observer une mappemonde construite sur ce principe.

Sphère

On voit bien que plus on se rapproche des pôles, c'est à dire aux extremums de la latitude, plus les continents semblent écrasés.

On à donc deux calculs à faire, couteux qui plus est. Le premier calcul est celui de la distorsion de la heightmap, le second calcul est celui du placement des vertex de la sphère dont on modifie l'altitude en fonction de la heightmap que l'on vient de transformer.

Sphère

Bon. Cette solution ne convient pas du tout. Elle ne convient pas déjà parce qu'elle est très couteuse en temps de calcul. Par ailleurs, elle à le désavantage majeur de ne pas être commode à manipuler. En effet pour notre planète, si nous voulons la peupler dynamiquement, on devra manipuler en temps réel les informations topologiques contenues dans la heightmap… Hors celle-ci aura été transformée. Ça ne va donc pas l'faire. On pourrait stocker une copie non déformée dont on se servirait comme référence, mais ça impliquerait à chaque fois que l'on pose des éléments sur la heightmap (des villes, des routes, des forêts) de leur appliquer ensuite la distorsion qui convient en fonction de la latitude. Enfin pour le niveau de détail dynamique -fonction de la distance camera/surface- ça va vite devenir très compliqué pour les mêmes raisons.

La projection équirectangulaire n'est donc pas terrible pour ce qu'on se propose de faire.

Non, ce qu'on va faire c'est… Projeter un cube sur une sphère ! La génération du cube n'étant pas compliqué en soi il est très facile ensuite de transformer l'objet en sphère. La distorsion induite par la transformation est minime et chaque face du cube est facilement manipulable à posteriori. Enfin il est possible de plaquer six heightmaps sur chacune des faces du cubes sans traitement préalable particulier.

Cela étant dit on va pouvoir attaquer la génération du cube. Pour ce faire on va déjà voir comment on s'y prend pour générer une face avec un nombre arbitraire de vertex. Nous allons voir que la question est un peu plus compliquée qu'il n'y paraît. Parce qu'il y a plusieurs façon d'organiser les vertex.

La première chose à savoir, c'est que la carte graphique, pour afficher des formes complexes, travaille en fait avec des primitives. C'est à dire des points, des segments, ou des triangles qui seront ensuite assemblés pour former un objet complexe. La plupart du temps on travaille avec des triangles parce qu'ils permettent de créer des surfaces sur lesquelles on peut afficher des textures et des effets.

Du coup, un carré, c'est au moins deux triangles.

squareByTwoTriangles.png

Mais nous on veut avoir avoir un carré, certes, mais composé de nombreux vertex, et donc de nombreux triangles pour pouvoir y appliquer la heightmap. Ce qui nous amène à la problématique d'assemblage des primitives. Dans le cas des triangles, OpenGL propose plusieurs façon de s'organiser. On en retiendra deux.

primitives.png

Dans le premier cas, très simple, il faut effectivement trois vertex par triangle. Ce qui veut dire que pour chaque cellule de notre surface carrée il faut compter six vertex. C'est pas génial. C'est pas génial parce qu'avec le mode GL_TRIANGLE_STRIP on peut drastiquement diminuer le nombre de vertex à envoyer à la carte graphique. Dans ce mode chaque triangle est dessiné en fonction du vertex courant, et des deux derniers vertex placés. Cela fait économiser énormément de vertex. Par exemple sur une ligne complète de carrés successifs. Disons 10 carrés on passe de soixante vertex (GL_TRIANGLES) à vingt vertex (GL_TRIANGLE_STRIP). Pas mal ! C'est donc très pratique, mais l'inconvénient, on s'en doute, c'est que la structure des vertex doit être repensée pour des formes complexes. Parfois, il n'est tout simplement pas possible de représenter une forme géométrique sans triangles dégénérés. Et même l'utilisation de ces triangles dégénérés peut être laborieuse dans certains cas. Heureusement pour générer une grille, c'est relativement simple. C'est ce que nous allons voir maintenant.

Une grille carrée va être composée de lignes de cellules carrées. Il y aura autant de lignes que de colonnes. Dans le cas de la première ligne, quand on arrive au bout on se rend compte que l'on a un problème.

degeneratedTriangles1.png

En effet le 7éme vertex ne permet pas de faire un triangle puisqu'il est parfaitement aligné avec les deux précédents. C'est là qu'interviennent les triangles dégénérés.

Un triangle dégénéré est un triangle aplati, ou même réduit à un point, quand chacun des trois points ont les mêmes coordonnées. Dans notre cas nous aurons besoin de deux triangles dégénérés. La raison est que le premier triangle dégénéré va inverser l'orientation du triangle courant, et de tout ceux qui suivront après. On rajoute donc un second triangle dégénéré (donc un vertex en plus) pour restaurer l'orientation des triangles. On ne va pas rentrer dans les détails mais l'orientation des triangles peut-être problématique plus tard.

degeneratedTriangles2.png

Ainsi le neuvième vertex et le huitième vertex forment le segment qui n'attend plus que le dixième vertex, placé sur le quatrième, clôturant ainsi le premier triangle de la seconde ligne. Il n'y a plus qu'a répéter l'opération pour chaque ligne.

Ces vertex doivent être sauvegardés quelque part, évidemment, et comme on veut éventuellement contrôler la résolution de la planète, on doit déterminer une formule qui indique le nombre de vertex à stocker – en vue d'allouer dynamiquement la mémoire – en fonction d'un paramètre n qui lui détermine la résolution.

La formule pour déterminer le nombre de vertex dans le cas d'une grille dont la résolution est une puissance de deux est donc la suivante :

2^2n * 2 -  2^n - 2

Où ici 2^n est la dimension de la grille.

Cela se démontre assez facilement si on se représente schématiquement la grille et que l'on compte le nombre de vertex en tenant compte des triangles dégénérés.

theGridADigitalFrontier.png

À présent on peut jeter un œil à l'algorithme qui permet de « dessiner » la grille tel que je l'ai implémenté et commenté.

Alors attention, on pourra remarquer qu'en fait je ne place pas exactement les vertex. Je sélectionne des indices qui pointeront sur une grille contenant n² vertex. Ce système est un mécanisme proposé par OpenGL qui permet encore une fois d'économiser de la mémoire et de la bande passante.

L'idée est la suivante : Un vertex est généralement un triplet de nombres à virgule flottante (float). Un vertex, donc, consomme 4 * 3 octets. Si on transfère directement les vertex à la GPU on consommera un nombre d'octets déterminé par la relation suivante

S(n) = (2^2n * 2 -  2^n – 2) * 4 * 3

Avec S, fonction de n, le nombre d'octets nécessaire.

Dans notre cas on crée une grille de n² vertex et on pointe sur chaque vertex à l'aide d'un entier court (short int) codé sur deux octets. L'espace mémoire s(n) consommé par les vertex et les indices devient alors

s(n) = (2^2n * 2 -  2^n – 2) * 2 + (2^n) * 4 * 3

pour un n relativement grand on gagne énormément de place. Pour s'en convaincre je vous propose un petit graphique avec en rouge S(n) et en bleu s(n). La différence est considérable.

memoryConsumption.png

Cela étant dit on peut enfin passer au code :

  // Permet de déterminer la direction du placement des vertex.
  bool reverse = false;

  // Coordonnées initiales 
  int x = 0;
  int y = 0;

  // Indice
  int i = 0;

  // On travaille dans une boucle while, ça me semblait plus commode:
  // On peut contrôler conditionnellement l'incrémentation de i.
  while (true) {

    // Les vertex sont placés par deux, verticalement.
    this->index[ i ]     = this->cubeScale * y     + x; 
    this->index[ i + 1]  = this->cubeScale * (y+1) + x;

    i+=2;

    // Quand on arrive au bout d'une ligne on change de direction
    if (x == this->cubeScale -1 || (x == 0 && i > 2)) {
      // Si en plus on arrive au bout de la dernière ligne, on arrête tout.
      if(y == this->cubeScale -2) {
        break;
      }

      // On change de direction
      reverse = !reverse;
      y++;

      // On Double le vertex courant
      this->index[ i ]     = this->cubeScale * y + x; 
      this->index[ i + 1]  = this->cubeScale * y + x;

      // On place le vertex qui formera le premier segment de la ligne
      this->index[ i + 2]  = this->cubeScale * (y+1)+x;
      i+=3;
    }

    // On se déplace le long de la ligne selon la direction
    if (reverse) {
      x--;
    }
    else {
      x++;
    }
  }

Rajoutons que l'on est pas obligé d'utiliser un tableau d'entiers courts. Nous pourrions utiliser des entiers normaux. La raison est que je me suis fixé une résolution maximale pour la grille que je juge suffisante, soit 7² vertex². On est donc limité par la taille de l'index, mais on gagne en rapidité (deux fois moins de données à transférer). La résolution relativement faible de la grille de vertex pourra être compensé avec d'autres astuces que nous verrons dans un autre article.

Maintenant que nous savons créer une grille (une surface carré), nous savons de fait créer un cube. Il ne reste plus qu'à transformer celui-ci en sphère. Alors attention cependant, dans mon implémentation chaque face du cube est rendu à l'écran séparément. Le cube n'est donc pas un objet indivisible, mais plutôt un objet dont chaque face est indépendante.

cube.png

Maintenant, c'est le moment de faire un peu de géométrie

Pour un point donné sur le cube on sait que la distance entre ce point et le centre du cube est donné par la relation

r = √( x² + y² + z²)

Tout à fait de la même façon qu'un carré inscrit dans un cercle trigonométrique.

trigo.png

En effet à l'intersection du carré sur le cercle, le rayon vaut un sur le cercle unité :

r = √(x² + y²)

Le rayon va varier entre 1 et √(1/2), soit 0,707106 quand on évalue le point sur le carré et non sur le cercle.

Pour projeter le point sur le carré (ou le cube dans notre cas)

On multiplie chaque coordonnée du point par l'inverse de son rayon. Par exemple dans le cas du cube nous aurons :

v = 1.0 / √( x² + y² + z²)

{
  X = x * v,
  Y = y * v,
  Z = z * v
}

Avec X,Y,Z les coordonnées finales. Et voilà ! Avec ça on a notre sphère !

sphere.png

On peut tout de même se fendre d'une petite optimisation en calculant une seule fois tous les coefficients pour une seule face. Ceux là étant calculés et stockés, il n'y a plus qu'a y accéder pour chaque face du cube.

  float * radiusPerVertex = new float[this->vertexSize];
  for (int y=0; y<this->cubeScale; y++) {
    for (int x=0; x<this->cubeScale; x++) {
      v1 = -0.5 + x * step;
      v2 =  0.5 - y * step;
      radiusPerVertex[this->cubeScale*y + x] = 1.0/sqrt( pow(v1, 2) +  pow(v2, 2) +  pow(0.5, 2) );
    }
  }
  
  // create vertices
  for (int y=0; y<this->cubeScale; y++) {
    for (int x=0; x<this->cubeScale; x++) {
      v1 = -0.5 + x * step;
      v2 =  0.5 - y * step;

      // Pour éviter de déréférencer 500 mille fois on stocke le rayon courant une bonne fois pour toute.
      v3 = radiusPerVertex[this->cubeScale*y + x];  

      this->vertex[0][this->cubeScale*y + x].x = v1 * v3;
      this->vertex[0][this->cubeScale*y + x].y = v2 * v3;
      this->vertex[0][this->cubeScale*y + x].z = 0.5 * v3;

      this->vertex[1][this->cubeScale*y + x].x = v1 * v3;
      this->vertex[1][this->cubeScale*y + x].y = v2 * v3;
      this->vertex[1][this->cubeScale*y + x].z = -0.5 * v3;

      this->vertex[2][this->cubeScale*y + x].x = v1 * v3;
      this->vertex[2][this->cubeScale*y + x].y = -0.5 * v3;
      this->vertex[2][this->cubeScale*y + x].z = v2 * v3;

      this->vertex[3][this->cubeScale*y + x].x = v1 * v3;
      this->vertex[3][this->cubeScale*y + x].y = 0.5 * v3;
      this->vertex[3][this->cubeScale*y + x].z = v2 * v3;
      
      this->vertex[4][this->cubeScale*y + x].x = 0.5 * v3;
      this->vertex[4][this->cubeScale*y + x].y = v1 * v3;
      this->vertex[4][this->cubeScale*y + x].z = v2 * v3;

      this->vertex[5][this->cubeScale*y + x].x = -0.5 * v3;
      this->vertex[5][this->cubeScale*y + x].y = v1 * v3;
      this->vertex[5][this->cubeScale*y + x].z = v2 * v3;
    }	
  }

Pour découvrir le PoC entier rendez vous ici : Les parties décrites dans cet article se trouvent dans le répertoire sphere/vertexFactory.cpp

On est bons pour la sphère, à la prochaine pour un nouvel article !

Merci à Benjamin Loubet pour l'aide et la relecture orthographique