Mot-clé - Génération procédurale

Fil des billets - Fil des commentaires

05/02/2016

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

04/06/2016

Génération procédurale de terrain: Part. II, Stamp Noise

Avant d'entamer cette seconde publication sur la génération procédurale, je vous propose encore une fois d'introduire quelques nouveaux termes qui seront utilisés ici.

  • Récursivité : La récursivité désigne en algorithmie le procédé par lequel une fonction donnée s'appelle elle-même. Un algorithme récursif s'oppose donc à un algorithme itératif qui, lui, répétera un nombre déterminé de fois un bloc d'instructions. Bien que tout problème résolu avec une approche récursive puisse être implémenté sous une forme itérative, et vice-versa, il est parfois clairement plus commode d'opter pour l'une ou l'autre solution selon le contexte. Dans tous les cas, chacune de ces approches nécessitent une condition d'arrêt.
  • Voxel: Un voxel est dans l'espace tridimensionnel ce qu'est un pixel sur une surface plane. Bien que son utilisation soit rare de nos jours dans le jeu vidéo, on le rencontre souvent en imagerie médicale.
  • Complexité algorithmique: Je vous renvoie à l'article que j'ai écrit sur le sujet.
  • Artefact: Cela désigne plusieurs choses, mais la plupart du temps, l'idée reste la même quel que soit le contexte. En imagerie de synthèse cela fait référence aux défauts visibles sur une image produite par un algorithme. Dans notre cas si l'on cherche à créer des surfaces procédurales avec pour ambition de leur donner un aspect naturel, alors on appelle artefact tout ce qui, sur ces surfaces, trahit le caractère artificiel de l'image. Typiquement, la présence de motifs se répétant, des formes géométriques visibles ou bien encore des distorsions topologiques anormales. Ces défauts peuvent être la conséquence d'un bug dans l'implémentation de l'algorithme ou un défaut dans la conception de l'algorithme lui-même. La précision des nombres à virgule flottante peut également produire ce genre d'effets.

Pour ce second article nous allons traiter d'un algorithme naïf, mais qui fonctionne bien. L'idée générale est tellement simple et facile d'implémentation qu'elle n'a pas vraiment de nom et n'est pas formellement décrite, à ma connaissance en tout cas. En cherchant un peu sur internet le procédé a cependant déjà été évoqué, notamment ici:

Creating Spherical Worlds

En ce qui me concerne, et pour cet article nous parlerons, de Stamp Noise.

Posons un peu le contexte pour bien comprendre et introduire le concept avant de commencer. Dans les années 80 nous n'avions pas encore de machines aussi puissantes qu'en l'an de grâce 2016. Ça tombe sous le sens. Lorsque l'on voulait générer un terrain, on le faisait généralement à partir d'une source de données volumineuse -une texture- qui devait être ensuite transformée en une surface en trois dimensions. Ken Perlin a révolutionné la technique avec son algorithme qui avait notamment la prodigieuse capacité de créer des surfaces à la volé, épargnant substantiellement la bande passante et la mémoire au point que ce genre de surface pouvait éventuellement être générée en temps réel quelques années plus tard. Mais ce n'est pas tout, ce genre de technique permet de s'affranchir, dans une certaine mesure, des problèmes que posent le placage de texture sur un objet 3D complexe: Plutôt que de générer une surface qu'il faut ensuite mapper sur l'objet en 3D, on peut directement travailler sur les vertices ou créer les voxels constituant cet objet quand, par exemple, l'algorithme de Perlin est utilisé en trois dimensions.

Entre temps, la technologie ayant évolué, et n'ayant plus les mêmes contraintes techniques, nous pouvons envisager la génération procédurale sensiblement différemment. Le premier problème que pose l'algorithme de Perlin, c'est qu'il n'est pas spécialement rapide. La faute à l'interpolation et le produit scalaire des vecteurs. Je vous renvoie à mon précédent billet pour vous rafraîchir la mémoire. L'autre problème majeur c'est que les terrains produits avec le Bruit de Perlin, quoique convainquant, sont très neutres. Il n'y a pas d'aspérités ou d'érosions qui confèrent au terrain un caractère particulier.

La stratégie que l'on peut envisager, à la lumière de l'éblouissante puissance de nos ordinateurs, c'est de limiter les calculs en stockant des motifs d'aspérités dans une matrice. On calcule une première fois un motif, et on l'imprime successivement sur notre champ de hauteur (heightmap) vierge, en répétant le procédé à différentes échelles. Ce qui est intéressant dans cette approche, c'est que l'on peut générer ainsi plusieurs tampons. Soit en les générant procéduralement, soit en les chargeant depuis une image. Le bénéfice immédiat c'est un gain substantiel de rapidité puisque l'on effectue les opérations mathématiques, éventuellement, durant la création des tampons et non plus dans la phase critique de la création de la heightmap. Par ailleurs, ces tampons peuvent être stockés durablement dans la mémoire pour être réutilisés, évidemment, à différentes échelles. L'autre bénéfice, c'est que les tampons seront potentiellement très différents et contribueront à donner à notre terrain une certaine personnalité.

Maintenant que nous avons introduit le sujet, nous allons rentrer un peu plus dans le détail de l'algorithme. Dans cet article nous distinguerons deux phases:

  • La création du tampon.
  • L'application du tampon.

J'avais dit tout à l'heure qu'il serait intéressant d'avoir plusieurs tampons. Pour cet article cependant, nous ne détaillerons qu'un seul tampon très simple. Un cône. C'est ce que nous allons voir maintenant.

Le rayon de notre cône aura une longueur de (2^n)/2 soit en fait 2^(n-1), avec n un nombre entier qui déterminera l'aire de notre surface. La technique consiste à parcourir une matrice de 2^(2n) éléments. Pour chaque points de la matrice on calcule le rayon du cercle imaginaire passant par le point courant. Ce cercle est sans surprise centré exactement au milieu de la matrice.

tampon1.png

Pour calculer le rayon, rien de plus simple, on utilise la mythique, ancestrale et légendaire relation de Pythagore.

notons

R = 2^(n-1)

le rayon du cercle inscrit dans le carré que forme le tampon. On a alors la relation suivante

r = √( (x-R)² + (y-R)²)

Une petite optimisation s'impose cependant. Cette opération va être répétée 2^2n fois, mais on peut d'abord pré-calculer dans une table de longueur n les puissances de deux de 0 à n-1. Ce qui allégera substantiellement nos calculs.

radius = sqrt( powersOfTwo[x] + powersOfTwo[y] )

Une fois le rayon calculé en un point donné on peut attribuer en ce point du tampon une valeur telle que celle-ci soit comprise entre 0 et 1.

tampon[x][y] = (halfScale - radius ) / halfScale

Voilà, nous sommes bon pour le tampon! Maintenant, il faut l'appliquer!

Globalement, c'est très simple. Notre algorithme est récursif, et à chaque appel, la fonction de bruit travaille sur un secteur, puis subdivise ce secteur en quatre tout en se rappelant quatre fois elle même à nouveau, et ainsi de suite, tant que le secteur à venir a une taille de deux au moins.

sectors.png

Pour chaque secteurs courants on centre le tampon mit à l'échelle. Comme on travaille avec des puissances de deux il est facile pour une coordonnée donnée de trouver sa correspondance dans le tampon.

inc = realScale / scale

où realScale est la dimension réelle de notre tampon/heightmap et où scale est la dimension du secteur courant. En parcourant dans une double boucle notre heightmap en un secteur donné, à chaque itération en x et en y on incrémente respectivement les coordonnées x et y relatives au tampon par inc.

Notons par ailleurs que notre double boucle ne parcourt pas toute la heightmap. Elle parcourt un SECTEUR de la heightmap. Nous aurons donc quelque chose de la forme

        ...
	for(x=0;x<scale;x++) {
	  stampY = 0;
	    for(y=0;y<scale;y++) {
              heighmap[x+offsetX][y+offsetY] += tampon[stampX][stampY];
	      stampY+=inc;
	    }
	  stampX+=inc;
	}
        ...

scale correspond évidemment ici à la taille du secteur courant.

Le cœur du bruit va être maintenant de déplacer le tampon sur son secteur. On n'y est pas obligé, mais dans notre implémentation, quand le tampon dépasse de la heightmap, on reboucle de l'autre côté de celle-ci. Pour déplacer le tampon on se donne deux valeurs pseudo-aléatoires de déplacement en x et en y.

randX = - halfScale + randomInteger % scale;
randY = - halfScale + randomInteger % scale;

avec scale la dimension de notre secteur courant et halfscale la demi dimension du secteur courant. randomInteger est évidement un nombre aléatoire dont la récupération est laissée à la discrétion du programmeur.

displacement.png

La boucle principale devient donc quelque chose de la forme

        ...
	for(x=0;x<scale;x++) {
	  stampY = 0;
	    for(y=0;y<scale;y++) {
              heighmap[x+offsetX+randX][y+offsetY+randY] += tampon[stampX][stampY];
	      stampY+=inc;
	    }
	  stampX+=inc;
	}
        ...

Dans ce pseudo code qui ressemble à s'y méprendre à du C/C++ on ne gère pas encore le cas où les sommes composants les coordonnées x et y dépassent positivement ou négativement de la heightmap. Il faut donc prévoir le coup et distinguer le cas simple où le tampon est bien contenu à l'intérieur de la heightmap et le cas plus rare où le tampon dépasse sur au moins un axe.

Un autre test optimisant substantiellement le nombre d'instructions à exécuter est celui de la valeur locale du tampon. Si cette valeur vaut zéro alors il est inutile de calculer les coordonnées, les offsets, les dépassements, les tests et autres affectations. En procédant ainsi on effectue approximativement 1.27 fois moins d'opérations dans le cas où le tampon est un cône. C'est propre!

Enfin, pour terminer, à chaque appel récursif on choisit pseudo-aléatoirement le type d'opération sur la heightmap. Est-ce qu'on creuse le terrain? Ou est-ce qu'on le surélève? Pour ce faire la valeur locale du tampon est multipliée par une variable qui vaut soit un ou moins un et qui est déterminée avant la double boucle principale.

// sign, c'est ce que je viens d'expliquer plus haut, à l'instant.
// inc divise le bruit courant, comme on en a déjà parlé, pour diminuer l'influence du bruit en fonction de l'octave.
heightmap[x][y] += sign * currentStampValue / (inc);

Voilà pour les grandes lignes! Je vous propose de jeter un œil à l'implémentation que j'ai faite et même de tester ça chez vous en récupérant les sources complètes ici, sur github: stampNoise

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <math.h>
#include "png.cpp"

// Le générateur de nombres pseudo aléatoires.
unsigned long int getRandom() {
        timespec tStruct;
        clock_gettime(CLOCK_REALTIME, &tStruct);
        return tStruct.tv_nsec;
}

// Le générateur de tampon
void UNIVERSE_STAMP_1(double * matrix, int scale) {
	int x,y,X=0;
	double halfScale = scale / 2 ;
	double radius;

        int * powersOfTwo = (int *) malloc( sizeof(int) * scale);

        // On crée deux tables contenants les valeurs élevées à la puissance de deux.
        // On calcul ainsi n fois ces valeurs au lieu de n².
        for(x=0; x<scale;x++) {
          powersOfTwo[x] = (x-halfScale) * (x-halfScale);
        }

        // Plutôt que d'incrémenter d'une unité x et calculer la position du point courant 
        // dans le tableau en multipliant par scale, on incrémente directement x par scale. La formule
        // pour retrouver le point courant n'est plus 
        // 
        // x*scale +y 
        //
        // mais 
        //
        // x+y
        //
        // Il faut donc également penser à changer la valeur limite de la première boucle
        int limit = scale*scale;
	for(x=0;x<limit;x+=scale) {
		for(y=0;y<scale;y++) {
                        // On calcule le rayon du cercle sur lequel se trouve le point courant.
                        // Opération très TRÈS gourmante en temps CPU
			radius = sqrtl( powersOfTwo[y] + powersOfTwo[X]);
			if ( radius < halfScale ) {
                                // y a plus qu'à dessiner le cône.
				matrix[x+y] = (halfScale - radius) / (halfScale);
			}
			else {
                                // Si on est en dehors du cercle, on se casse pas la tête et on affecte un zero.
				matrix[x+y] = 0;
			}
		}
                //Comme x est incrémenté par scale, et qu'on doit quand même accéder à notre tableau powersOfTwo...
                X++;
	}
        free(powersOfTwo);
}

void UNIVERSE_STAMP_NOISE(double * matrix, double * stamp, int scale, int offsetX, int offsetY, int realScale) {

        // La condition d'arrêt de notre bruit récursif.
        // Selon la granularité que l'on désire, on peut augmenter la valeur limite de scale.
        if (scale == 1) {
		return;
	}

        // Demi dimension courante
        // Comme scale est une puissance de deux, plutôt que de diviser, on opère une rotation binaire.
	int halfScale = scale >> 1;
	int x, y;

        // Deux variables très importantes, ce sont elles qui déterminent ou sera appliqué le tampon.
        // C'est le positionnement aléatoire qui fait toute la "beauté" de la heightmap.
	int randX = - halfScale + getRandom() % scale;
	int randY = - halfScale + getRandom() % scale;

        // À chaque octave il faut diminuer l'influence du bruit.
	// On se sert également de cette variable comme pas d'incrémentation des
        // coordonnées du tampon.
        int inc = realScale / scale;

        // Deux variables incrémentales qui servent à récupérer le pixel locale au tampon, en fonction de l'octave.
        // Elles sont toute les deux incrémentés avec la valeur de inc.
	int stampX=0, stampY=0;

        // Nouvelles coordonnées si dépassement du tampon en dehors de la heightmap
	int wrapX,wrapY;

        // Détermine le signe du tampon.
        // S'il est positif, le terrain se surélève, à l'inverse, il se creuse
	float sign = getRandom() & 2 ? -1.0 : 1.0;

        int tmpCoordX,tmpCoordY;
        double currentStampValue;

	for(x=0;x<scale;x++) {
	  stampY = 0;
	    for(y=0;y<scale;y++) {
              // On économise des calcules fastidieux en stockant cette valeur qui sera solicitée au moins une fois.
              currentStampValue = stamp[stampX*realScale+stampY];

              // Avec ce test le gros bloc d'instructions est répété 1.27 fois moins que s'il n'y avait pas de test.
              if (currentStampValue != 0) {
                // On économise des calculs fastidieux en stockant ces valeurs qui seront solicitées plusieurs fois.
                tmpCoordX = randX + offsetX + x;
                tmpCoordY = randY + offsetY + y;

                // Le cas simple où le tampon ne dépasse pas de la heightmap
	        if ( tmpCoordX < realScale && tmpCoordX >= 0 && tmpCoordY < realScale && tmpCoordY >= 0) {
	          matrix[ (tmpCoordX * realScale) + tmpCoordY] += sign * currentStampValue / (inc);
	        }

	        // Là c'est plus pénible, il faut calculer le décalage à appliquer selon le ou les côtés où le pixel à dépassé.
	        else {
                  // On restore les coordonnées du décalage
                  // Comme il se peut que le pixel ne dépasse que sur un axe, par défaut, le décalage est fixé à zero.
		  wrapX = 0;
	    	  wrapY = 0;
                  // Le pixel dépasse à droite
		  if ( tmpCoordX >= realScale && tmpCoordX < realScale*2 ) {
      		    wrapX = - realScale;
		  }

                  // Le pixel dépasse à gauche
		  else if ( tmpCoordX > -realScale && tmpCoordX < 0) {
		    wrapX = realScale;
	  	  }

                  // Le pixel dépasse en haut
		  if ( tmpCoordY > -realScale && tmpCoordY < 0) {
    		    wrapY = realScale;
		  }

                  // Le pixel dépasse en bas
		  else if ( tmpCoordY < realScale * 2 && tmpCoordY >= realScale) {
		    wrapY = -realScale;
		  }

                  // On peut maintenant repositionner le pixel sur la heightmap.
                  // la coordoonée final dans un tableau simulant une matrice est de la forme:
                  //
                  // (X * hauteur) + Y
                  // Avec X valant la somme du 
                  //  décalage du secteur courant en abcisse, 
                  //  du décalage du tampon courant en abcisse,
                  //  la coordonnée x courante,
                  //  le décalage aléatoire en abcisse
                  //
                  // Avec Y valant la somme du 
                  //  décalage du secteur courant en ordonnée, 
                  //  du décalage du tampon courant en ordonnée,
                  //  la coordonnée y courante,
                  //  le décalage aléatoire en ordonnée
		  matrix[ (  wrapX + tmpCoordX) * realScale + tmpCoordY + wrapY] += sign * currentStampValue / (inc);
      	        }
              }
              // On incrémente à l'échelle la coordonnée locale au tampon
	      stampY+=inc;
	    }
          // On incrémente à l'échelle la coordonnée locale au tampon
	  stampX+=inc;
	}

        // En divisant par deux la dimension courante à chaque récursion, et en modifiant l'offset,
        // on subdivise en permanence la heighmap jusqu'à ce que la dimension ainsi divisée soit égale à un.
        // En procédant ainsi, on travaille récursivement sur différentes
        // portions de la heighmap. Il y a donc quatre portions par secteur et à chaque récursion, chacunes
        // des portions devient elles même un secteur.

        // Portion en haut à gauche du secteur courant
	UNIVERSE_STAMP_NOISE(matrix, stamp, scale/2, offsetX+0, offsetY+0, realScale);

        // Portion en bas à gauche du secteur courant
	UNIVERSE_STAMP_NOISE(matrix, stamp, scale/2, offsetX+0, offsetY+scale/2, realScale);

        // Portion en bas à droite du secteur courant
        UNIVERSE_STAMP_NOISE(matrix, stamp, scale/2, offsetX+scale/2, offsetY+scale/2, realScale);

        // Portion en haut à droite du secteur courant
	UNIVERSE_STAMP_NOISE(matrix, stamp, scale/2, offsetX+scale/2, offsetY+0, realScale);
}

int main(int argc, char ** argv) {
	if (argc == 1) {
		std::cout << "Usage: ./heightMap <power of two>\n";
		return 0;
	}

	int scale = atoi(argv[1]);
	int x,y,i,j;
	double * matrix;
	double * stamp;
	matrix	= (double *) malloc( sizeof( double ) * scale * scale);
	stamp 	= (double *) malloc( sizeof( double ) * scale * scale);
        PIXEL ** png = (PIXEL **) malloc(sizeof(PIXEL *) * scale);

        // On s'assure que la matrice de destination sera bien "vierge"
	for (x=0; x<scale;x++) {
		for(y=0;y<scale;y++){
			matrix[x*scale+y] = 0;
		}
		png[x] = (PIXEL *) malloc(sizeof(PIXEL) * scale);
	}

        // On génère d'abord notre tampon
	UNIVERSE_STAMP_1(stamp, scale);

        // On lance notre bruit récursif.
        // On conmmence la récursion avec l'octave la plus grande.
	UNIVERSE_STAMP_NOISE(matrix, stamp, scale, 0, 0, scale);

        // À partir d'ici, la heightmap est terminé. Il n'y a plus qu'à déterminer les extremums
        // pour normaliser la hauteur.

        long double max=0,min = 65536;
        for (x=0; x<scale;x++) {
                for(y=0;y<scale;y++) {
                        if (matrix[x*scale+y] > max) {
                                max = matrix[x*scale+y];
                        }
                        if (matrix[x*scale+y] < min) {
                                min = matrix[x*scale+y];
                        }
                }
        }

        char color;
        for (x=0; x<scale;x++) {
                for(y=0;y<scale;y++) {
                        color = (((matrix[x*scale+y]) - min ) * 255) / (max - min);
                        png[x][y].Alpha = 0xFF;
                        png[x][y].Red = color;
                        png[x][y].Green = color;
                        png[x][y].Blue = color;
                }
        }

        // On transfére notre heightmap dans un fichier png...
        writePng(png,scale);
	return 0;
}

Et donc, qu'est-ce que ça donne quand on passe la heightmap ainsi générée en postprod?

stampNoise0x03.png

stampNoise0x02.png

stampNoise0x01.png

C'est tout à fait convaincant! Cependant l'algorithme, quand il est utilisé avec un cône uniquement, produit des artefacts... En effet on peut voir apparaître, si on est attentif, des cercles trahissant l'usage d'un objet circulaire pour modeler le terrain.

artefact.gif

Pour le moment, ce n'est pas le drame, mais il peut être intéressant, voir nécessaire, d'utiliser d'autres formes de tampons. Dans le cas du cône une possible solution pourrait être de combiner la formule d'interpolation du bruit de Perlin amélioré avec celle utilisée pour produire le cône:

tampon[x][y] =  interpolation ( (halfScale - radius ) / halfScale)

rappelons que la formule d'interpolation utilisée dans le bruit de Perlin amélioré est

6t⁵- 15t⁴ + 10t³

Cela aurait pour effet d'adoucir les transitions entre la base du cône et son sommet. Cône qui du coup n'en serait plus un et qui deviendrait plutôt une sorte de cloche.

Maintenant qu'en est-il du temps d'exécution et de la complexité de notre algorithme? Toujours sur un core i3, je mesure:

  • 3ms pour 128 pixels²
  • 8 ms pour 256 pixels²
  • 30ms pour 512 pixels²
  • 123ms pour 1024 pixels²
  • 530ms pour 2048 pixels²
  • 2530ms pour 4096 pixels²

C'est beaucoup beaucoup mieux que l'algorithme de Perlin! Pour un n donné avec le Bruit de Perlin on peut générer une surface deux fois plus grande avec approximativement le même temps! En ce qui concerne la complexité du Stamp Noise je la définis ainsi:

C(n) = 2^n + π.(n+1).(2^(2n))/4

Graphiquement ça donne

theoreticalStampNoiseComplexity.png

Et quand on mesure les temps effectivement mesurés en millisecondes, toujours en fonction de n

experimentalStampNoiseComplexity.png

On est dans les clous!

Merci à Théo Fish et Benjamin Loubet pour la relecture!

03/11/2016

Génération procédurale de terrain: Part. I, Le Bruit de Perlin amélioré

Avant de se lancer dans la génération procédurale de terrain nous allons rapidement introduire quelques termes importants dont il sera question dans cet article.

  • Heightmap: Ou carte de hauteur. Il s'agit d'une matrice ou d'une image dont chaque indice ou pixel correspond à une élévation locale de la surface qu'on souhaite représenter. On s'en sert donc notamment pour modéliser des terrains.
  • Bruit: Il s'agit d'un signal aléatoire à n dimensions. Typiquement le son de la mer ou la neige du téléviseur. En informatique créer du bruit c'est en fait chercher à produire des valeurs aléatoires qu'on pourra alors utiliser en cryptographie, par exemple, ou dans notre cas, en imagerie de synthèse.
  • Vertex: En infographie 3D, un vertex (vertices au pluriel) est un point dans l'espace repéré par trois coordonnés X,Y et Z. Ce point dans l'espace, quand il y en a plusieurs est ensuite organisé et assemblé avec d'autres vertices en primitives (des triangles ou des segments) qui eux même formerons des objets complexes en 3D.
  • Displacement Mapping: En plaquant une texture de type heightmap sur une surface en 3D, c'est le procédé par lequel on déplace un vertex en fonction de l'intensité du pixel qui lui est associé dans la heightmap.
  • Bumpmap: Il s'agit d'une technique introduite par James Blinn en 1978 permettant de simuler des aspérités sur une surface sans changer les vertices composant celle-ci. C'est donc très pratique lorsque l'on veut ajouter du réalisme tout en faisant l'économie de vertex.
  • Vecteur: Pour faire très simple, si on considère qu'un point est une position dans l'espace, alors un vecteur est un déplacement dans l'espace. Ce qui implique qu'un vecteur possède une direction et une longueur. Ces propriétés peuvent être calculé à partir des coordonnées du vecteur.
  • Gradient: La notion de gradient est très fortement lié à la notion de dérivée et de dérivée partielle. Dans l'un ou l'autre cas on multiplie le résultat par son vecteur unité respectif. On peut donner en exemple les deux expressions ci-dessous, mais il est inutile de trop s'attarder ici sur la notion de gradient parce que dans l'algorithme du Bruit de Perlin les vecteurs gradients sont en fait des vecteurs de norme 1. grad1.png grad2.png
  • Interpolation: Il s'agit d'un procédé qui consiste à créer une transition continue entre deux valeurs discrètes. Mathématiquement ça se traduit par le fait de joindre deux points par une courbe. Pour cela plusieurs formules existent selon le contexte dans lequel on travaille: Interpolation linéaire, cubique, cosinusoïdale, polynomiale, etc... Le choix du type d'interpolation est crucial quand on le met en œuvre dans une application dont la rapidité d'exécution impacte sur le bon fonctionnement du programme. Une interpolation linéaire sera par exemple plus rapide mais moins "jolie" qu'une interpolation polynomiale beaucoup plus efficace mais plus gourmande en temps d'exécution.

Ce qu'on va donc chercher à faire c'est générer une heightmap. La heightmap peut alors être utilisée d'au moins deux façons simultanées si, par exemple, on souhaite faire une planète. Dans le premier cas la heightmap déterminera l'altitude (displacement mapping) de chaque vertex appartenant à la sphère de sorte de faire apparaître des montagnes et des continents. Dans le second cas, souhaitant économiser des ressources par souci d'optimisation, on peut réduire le nombre de vertex et utiliser la heightmap également en tant que Bumpmap. Quoi qu'il en soit il est très probable que la heightmap soit plus élevée en résolution que le maillage auquel on veut l'appliquer, on a donc tout intérêt à combiner les deux techniques.

Pour ce premier article nous allons présenter et détailler le Bruit de Perlin.

Alors, comme on s'en doute, Il s'agit d'une texture procédurale, c'est à dire d'une technique d'imagerie de synthèse. Ce procédé a été inventé par Ken Perlin lorsqu'il travaillait sur le film TRON en 1981. Cet algorithme permet de créer des surfaces d'apparence réaliste et est très utilisé dans le jeu vidéo ou le cinéma. Dans le titre de cet article je parle de Bruit de Perlin amélioré. Ce n'est pas moi qui l'ai amélioré, mais Ken Perlin lui-même peu de temps après la conception de la première version. Dans cet article nous ne parlerons donc que de cette version améliorée sans s'attarder sur les différences.

La fonction Bruit de Perlin, quand elle est utilisée dans un espace en deux dimensions, est une application de R² dans R. C'est à dire qu'elle prend en paramètres deux coordonnées réelles, et retourne une seule valeur réelle qui correspond à la valeur du bruit locale.

Pour ce faire, l'algorithme va fonctionner en deux temps:

Dans la première partie de l'algorithme, on se donne une grille.

perlinGrid.png

Pour chaque point P à coordonnées réelles on considère les nœuds Q à coordonnées entières de la grille qui encadrent ce point P. Par exemple si P = (2.3, 3.7) alors Q1 = (2,3), Q2 = (3,3), Q3 = (3,4) et Q4 = (2,4). La conséquence immédiate c'est que tant que P se situe entre 2 et 3 en x, et y entre 3 et 4, Q1,Q2,Q3 et Q4 ne changeront pas.

On assigne alors à chacun de ces points Q un vecteur gradient pseudo aléatoire. Pour ce faire on dispose du tableau de gradients suivant dans lequel, justement, nous irons piocher pseudo aléatoirement nos vecteurs.

  • 0,1
  • 0,-1
  • 1,0
  • -1,0
  • 1,1
  • 1,-1
  • -1,-1
  • -1,1

Avant d'aller plus loin, il faut expliquer de quelle manière on récupère nos indices aléatoires. Dans l'algorithme de Perlin, On ne fait à aucun moment appel à la fonction rand(). En fait, on dispose d'un tableau de 256 entrées (certaines implémentations en comportent le double). Ce tableau contient en fait tous les entiers de 0 à 255 dans le désordre. Ce qui signifie que soit ce tableau est fourni comme tel, inscrit en dur dans le programme, soit celui-ci est généré AVANT de lancer l'algorithme. Ce tableau dont on parle est ce qu'on appelle une table de permutations.

La formule pour déterminer l'indice du gradient dans le tableau est donc de cette forme là quand l'algorithme est utilisé dans un espace en deux dimensions:

Gn = P[ ( Qnx + P[ Qny & 255 ] ) & 255 ] & 7
  • Avec Gn l'indice du gradient dans le tableau de vecteurs, associé au point à coordonnées entière Qn.
  • Avec Qnx et Qny, respectivement, les coordonnées entières en x et y de Qn.
  • Avec P qui, ici, représente notre table de permutation à 256 valeurs.

L'usage des masques de bit sert à éviter d'utiliser des opérations modulo gourmandes en temps CPU quand les valeurs Qnx et Qny sont supérieures à 255.

Nous en parlerons plus bas, mais en l'état cette formule n'est pertinente que si l'on utilise notre Bruit de Perlin pour une octave donnée. Elle ne convient donc pas si l'on souhaite faire, par exemple, un terrain procédural. Ou plutôt, ça fonctionnera, mais le Bruit produira des défauts sous la forme de motifs qui se répéteront et qu'un œil pourtant distrait remarquera sans aucun doute. Nous y reviendrons.

À l'aide de la formule donnée nous avons donc quatre vecteurs gradients associés respectivement à chaque points Q.

perlinGrid2.png

L'algorithme nous dit maintenant qu'il faut calculer, pour chaque point Qn, le produit scalaire Gn · (P-Qn). Notons au cas où ce ne serait pas clair qu'ici P ne fait pas référence à la table de permutations mais au point P dont les coordonnées ont été données en paramètre à notre fonction Bruit de Perlin. Rappelons également que le produit scalaire, quand il est effectué dans une base orthonormale est simplifié. C'est le cas ici. On a donc:

Gnx * (x-Qnx) + Gny * (y-Qny);

Avec respectivement, Gnx et Gny correspondant à la coordonnée x et y du gradient Gn

Nous arrivons maintenant à la seconde partie de l'algorithme. Une fois que nous avons calculé nos produits scalaires en chaque point Qn il faut maintenant interpoler ceux-là pour en extraire la valeur au point P. Valeur qui sera, vous l'avez deviné, celle qui sera retournée par notre fonction Bruit de Perlin.

Notons s,t,u et v respectivement le résultat des produits scalaire en Q1,Q2,Q3 et Q4. Nous avons alors à calculer trois interpolations.

  • La première, horizontale, entre s et t dont on notera le résultat iST
  • La seconde, horizontale, entre u et v dont on notera le résultat iUV
  • La troisième, verticale, entre iST et iUV qui sera le résultat retourné par la fonction de bruit.

perlinGrid3.png

Pour interpoler, et bien il nous faut une formule. Nous pourrions utiliser une interpolation linéaire mais le rendu final en souffrirait. On ne détaillera pas ici les raisons exactes de ce choix, mais Ken Perlin propose la formule suivante:

6t⁵- 15t⁴ + 10t³

La formule ci-dessus à la propriété intéressante, si l'on travaille dans l'intervalle {0,1}, de s'annuler en ses extremums, et de le faire également en sa dérivée et dérivée seconde.

Et nous allons justement travailler dans cet intervalle. Rappelons que la coordonnée x du point P varie entre Q1x et Q1x + 1, avec Q1 le point à coordonnée entière qui correspond au sommet en haut à gauche de notre carré qui encadre notre point P. Autrement dit le nombre x - Q1x appartient à l'intervalle {0,1}. Et là, normalement, vous me voyez venir. Ce qu'on va faire, c'est qu'on va injecter dans notre formule d'interpolation la valeur x-Q1x de sorte à récupérer un coefficient dont on se servira pour calculer la valeur interpolée en un point donné. Par exemple, pour calculer l'interpolation horizontale des valeurs s et t, respectivement au point Q1 et Q2 nous ferions:

iST = s +interpolation( x -Q1x ) * (t - s)

Il n'y à plus qu'a appliquer le calcul pour déterminer iUV et recommencer une dernière fois pour interpoler iST et iUV. Et maintenant... Et bien... C'est fini, du moins en théorie.

Nous avons évoqué tout à l'heure, sans plus de détail, le fait de générer du bruit pour une octave donnée. Cette notion est très importante. Dans notre Bruit de Perlin, nous avons mentionné au début de notre article une grille. Mais cette grille est une grille unitaire, pour avoir du sens dans une application concrète il faut que cette grille soit mise à l'échelle. Et c'est précisément l'octave qui détermine l'échelle de notre grille. Selon l'octave que l'on se donne le rendu final du bruit sera plus 'fin' ou 'épais'. Rien ne vaut une image pour illustrer mon propos. Ci-dessous, pour une image de 512x512 dans laquelle on fait varier l'octave d'une puissance de deux à la suivante jusqu'à 512, en partant de deux, pas de un (sinon l'image est complètement noire, l'algorithme ne fonctionne pas). Plus l'octave est grande plus le bruit est épais.

perlinNoiseOctaves.png

Bon mais alors attention, si vous avez suivi, vous vous souvenez sans doute que j'ai dit que la formule permettant de sélectionner pseudo aléatoirement un gradient ne fonctionnait que pour une octave donnée. En effet, lorsque l'on veut générer une heightmap, on additionne pour une coordonnée donnée chaque bruit local à cette coordonnée pour chaque octave que l'on se donne. Mais voilà, on a bien dit que nos gradients étaient choisis pseudo aléatoirement à l'aide d'une formule qui est fonction des coordonnées qu'on lui donne.

Si on s'en tient à cette formule plus haut, on va se retrouver dans la situation fort déplaisante où ce que fera en réalité l'algorithme, pour chaque octave, c'est produire le MÊME bruit à une échelle différente. Et ce n'est pas bon du tout.

perlinHomotetie.png

Nous, ce que nous voulons, c'est produire du bruit différent à chaque octave.

La formule

Gn = P[ ( Qnx + P[ Qny & 255 ] ) & 255 ] & 7

devient donc

Gn = P[ ( Qnx + P[ Qny & 255 ] + octave ) & 255 ] & 7

Pour parachever notre algorithme, il faut assigner à chaque bruit local, en fonction de l'octave, un coefficient. On part de l'octave la plus grande avec un coefficient fixé à un. Et pour chaque itération, on divise par deux l'octave et on divise par deux notre coefficient.

Bon alors c'est super, mais visuellement qu'est-ce que ça donne?

proceduralHeightmap0x0B.png

Et quand on passe ça dans GIMP pour faire ressortir les aspérités avec de la lumière?

proceduralHeightmap0x0C.png

C'est du plus bel effet!

Pour terminer, je vous propose un exemple d'implémentation de mon cru relativement optimisée et abondamment commentée dont vous pouvez trouver les sources complètes ici:

https://github.com/DenisSalem/UNIVERSE/tree/master/PoC/improvedPerlinNoise

#include <stdio.h>
#include <stdlib.h>
#include "png.c"

// On définit une macro d'interpolation, à priori plus rapide qu'un appel de fonction.
// la formule est la suivante:  6t^5 - 15t^4 + 10t^3
#define interpolation(t) (6 * t * t * t * t * t - 15 * t * t * t * t + 10 * t * t * t)

//  Certaines variables et tableaux sont globaaux afin de ne pas être empilés lors de
//  l'appel de la fonction PerlinNoise2D, ce qui nous fera économiser des instructions
//  et du temps CPU.

// G est le tableau de gradient.
// Ici G comporte 8 paires de coordonnées. Dans l'implémentation améliorée G comporte
// 12 triplets de coordonnées pour pouvoir travailler dans un cube, et comme ce n'est pas
// le cas ici... :)
int G[8][2] = {
    {1,1},    {-1,1},   {1,-1},   {-1,-1},
    {1,0},    {-1,0},   {0,1},    {0,-1},
};

// Chaque vecteur gradient est référé par un indice Gn du tableau G
int G1,G2,G3,G4;

// Contient le résultat du produit scalaire de Gn avec P-Q dans une base orthonormale.
double s,t,u,v;

//Contient les coordonnées X et Y des points Qn et du point P.
int Q1[2] ,Q2[2],Q3[2],Q4[2];
double p[2];

// Résultat de la macro d'interpolation en X et Y, tmp permet de stocker P-Q avant d'être
// injecté dans la macro.
double iX,iY,tmp;

// Résultat de l'interpolation horizontal en st et uv
double iST;
double iUV;

// La table de permutations. En l'état, l'algorithme produira toujours le même terrain 
// pseudo-aléatoire. Pour obtenir un vrai terrain pseudo-aléatoire différent à chaque
// lancement du programme il faudrait changer ou désordonner cette table de permutations
// avant d'entrer dans la boucle principale ou est appelée notre fonction de bruit.
int P[256] = {
 235,249,14,239,107,49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,
 127, 4,150,254,138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,
 112,104,218,246,97,228,251,34,242,193,238,210,144,12,191,179,162,241,81,51,145,
 153,101,155,167, 43,172,9,129,22,39,253, 19,98,108,110,79,113,224,232,178,185,
 74,165,71,134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,
 187,208, 89,18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186,3,
 47,16,58,17,182,189,28,42,223,183,170,213,119,248,152, 2,44,154,163, 70,221,
 64,52,217,226,250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,
 203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168, 68,175,
 105,92,41,55,46,245,40,244,102,143,54, 65,25,63,161,1,216,80,73,209,76,132,
 142,8,99,37,240,21,10,23,190, 6,148,247,120,234,75,0,26,197,62,94,252,219,
 151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,
 156,180};


// Les coordonnées x et y doivent être des nombres réels positifs, sinon nous serions en 
// dehors de notre grille. Pour cela on envoie à notre fonction les coordonnées d'un
// point de la matrice, que l'on remet ensuite à l'échelle de la façon suivante
// C = c / scale
// avec C la coordonnée mise à l'échelle et c la coordonné entière d'origine.

// La fonction prend également en paramétre un entier -l'octave-
// positif qui détermine la taille de notre grille. Par exemple si la matrice
// de destination à une résolution de 256² éléments et que notre entier vaut
// 128 alors alors la taille de la grille sera en fait de 2² éléments, soit 3²
// noeuds :
//
//
//    <-- 256 --> taille réelle de la matrix
//
//    0----1----2
//    |    |    |
//    3----1----5   Here stands the Grid. A digital frontiere...
//    |    |    |
//    6----7----8
//
//    <--- 2 ---> taille de la grille mise à l'échelle
    
double PerlinNoise2D(int x, int y, unsigned int scale) {

  // L'initialisation peut sembler alambiquée, mais c'est 
  // pour épargner au processeur des calculs inutiles.

  p[0] = (double) x / scale;
  p[1] = (double) y / scale;

  // Coin supérieur gauche
  Q1[0] = (int) p[0];
  Q1[1] = (int) p[1];
  
  // Coin supérieur droit
  Q2[0] = Q1[0] + 1;
  Q2[1] = Q1[1];

  // Coin inférieur droit
  Q3[0] = Q2[0];
  Q3[1] = Q1[1] + 1;

  // Coin inférieur gauche
  Q4[0] = Q1[0];
  Q4[1] = Q3[1];


  // On récupére pseudo aléatoirement les gradients.
  // Pour éviter la répétition de motifs fractals à chaque octave 
  // l'indice final dépend de l'octave courant
  G1 = P[ (Q1[0] + P[ Q1[1] & 255] + scale) & 255 ] & 7;  // Gradient supérieur gauche
  G2 = P[ (Q2[0] + P[ Q2[1] & 255] + scale) & 255 ] & 7;  // Gradient supérieur droit
  G3 = P[ (Q3[0] + P[ Q3[1] & 255] + scale) & 255 ] & 7;  // Gradient inférieur droit
  G4 = P[ (Q4[0] + P[ Q4[1] & 255] + scale) & 255 ] & 7;  // Gradient inférieur gauche

  // On calcule le produit scalaire Gn . (P-Qn)
  // Avec P faisant référence aux coordonnées du point stocké dans p.
  // (P étant la table de permutations)
  s = G[G1][0] * (p[0]-Q1[0]) + G[G1][1] * (p[1]-Q1[1]);
  t = G[G2][0] * (p[0]-Q2[0]) + G[G2][1] * (p[1]-Q2[1]);
  u = G[G3][0] * (p[0]-Q3[0]) + G[G3][1] * (p[1]-Q3[1]);
  v = G[G4][0] * (p[0]-Q4[0]) + G[G4][1] * (p[1]-Q4[1]);

  
  // On interpole en x sur le segment st et uv
  tmp = p[0]-Q1[0];
  iX = interpolation(tmp);

  iST = s + iX * (t-s);
  iUV = v + iX * (u-v);

  // On interpole y sur le segment iST et iUV
  tmp = p[1]-Q1[1];

  iY = interpolation(tmp);

  // On retourne le résultat normalisé.
  return (1 + iST + iY * (iUV - iST)) * 0.5;
}

int main(int argc, char ** argv) {
  if (argc == 1) {
    printf("usage: ./improvedPerlinNoise <power of two>\n");
    return 0;
  }
  int scale = atoi(argv[1]);


  // Matrice de 256 pixels² simulé avec un tableau de longueur 256²
  // C'est dans ce tableau que nous allons stocker notre heightmap
  double * grid = (double *) malloc(sizeof(double) * scale * scale);
  
  // i et j correspondent respectivement aux axes x et y.
  // k correspond à l'octave courrante.
  int i,j,k;

  float min=2,max=0;

  // Selon le type de texture on peut ne pas utiliser de coef ou
  // l'utiliser différemment. Mais l'idée ici est de diminuer
  // l'influence du bruit à mesure que la fréquence augmente.
  double coef = 1.0; 

  for(j=0; j<scale;j++) {
    for(i=0; i<scale;i++) {
      coef = 1.0;
      grid[ i + j*scale] = 0;
      for(k=scale/2;k>=1; k/=2) {
        grid[ i + j*scale ] += PerlinNoise2D(i,j, k) * coef;
        coef /= 2.0;
      }
      if (min > grid[ i + j*scale]) {
        min = grid[ i + j*scale];
      }
      if (max < grid[ i + j*scale]) {
        max = grid[ i + j*scale];
      }
    }
  }

  // Ici la texture est terminée. Il ne reste plus qu'à la normaliser en
  // vue de l'exploiter.

  double normalizedIntensity;
  PIXEL ** matrix = (PIXEL **) malloc( scale * sizeof(PIXEL * ));
  for(j=0;j<scale;j++) {
    matrix[j] = (PIXEL *) malloc( scale * sizeof(PIXEL));
    for (i=0; i< scale;i++) {
      normalizedIntensity = (grid[ i + j * scale ] - min) * ( 1.0 / (max-min) );
      matrix[j][i].Alpha  = 255;
      matrix[j][i].Red    = (char) (normalizedIntensity * 255);
      matrix[j][i].Blue   = (char) (normalizedIntensity * 255);
      matrix[j][i].Green  = (char) (normalizedIntensity * 255);
    }
  }

  writePng(matrix, scale);
  return 0;
}

Avec cette implémentation, sur un core i3, et en ne gardant que la partie qui génère le bruit, on s'en tire en moyenne avec ces temps d’exécution:

  • 29 ms pour une surface de 256 pixels²
  • 126 ms pour une surface de 512 pixels²
  • 556 ms pour une surface de 1024 pixels²
  • 2249 ms pour une surface de 2048 pixels²

Est-ce qu'on peut mieux faire? Y a-t-il d'autres algorithmes?

La réponse dans le prochain article!

Merci à Benjamin Loubet et Théo Fish pour la relecture orthographique!