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!