Bruits et nombres aléatoires cohérents

Introduction du cours

Bonjour à tous !

Si vous avez déjà essayé de générer un terrain, une texture de nuage ou de marbre, de simuler l'oscillation des herbes sous l'effet du vent, vous vous êtes vite rendu compte qu'il faut intégrer un facteur aléatoire à l'algorithme.

Prenons par exemple la génération d'un terrain, sous forme d'une heightmap (donc un tableau bidimensionnel stockant l'altitude de chaque point de la carte). L'idée « naïve » qui vient souvent est de faire une double boucle et de remplir le tableau avec une fonction aléatoire rand(). Comme on s'en rend vite compte, les résultats ne sont pas très glorieux, on obtient de la neige comme sur les vieux postes de télévision. Il va donc falloir creuser un peu plus afin d'avoir un effet à la fois aléatoire et cohérent.

Pour avoir un effet cohérent, il faut préserver une certaine régularité dans les valeurs générées. Par là, on entend qu'il ne faut pas avoir de trop grands écarts entre deux valeurs successives. Le problème, c'est que plus on augmente la cohérence, plus on diminue le côté « aléatoire »...

La plupart des langages actuels fournissent des outils facilitant la génération de nombres aléatoires, mais si on veut quelque chose de cohérent, il faudra les retravailler soi-même. Pour cela, on utilise plusieurs fois le même jeu de valeurs, en les combinant tout en accordant plus d'importance à celles qui sont les plus cohérentes.

Voici quelques exemples expliqués dans ce tutoriel :

A propos de l'appellation "bruit de Perlin"
(Cette note est destinée à lever une ambiguïté sur l’appellation "perlin noise", si à ce stade vous n'y comprenez rien, passez juste sans vous en préoccuper).

On réfère couramment sur internet et même dans la littérature à l'algorithme que je vais introduire comme l'algorithme de Perlin, ou le bruit de Perlin (Perlin noise). Cette nomination est pourtant erronée.

En fait la technique ici employée de sommation de bruits, bien que souvent utilisée avec des bruits de Perlin, n'apparaît même pas dans l'article original de Ken Perlin. La génération d'une fonction de bruit est aussi totalement différente. Le bruit de Perlin est un bruit dit "gradient" (gradient noise), alors que les bruits utilisés ici sont des "bruits de valeurs" (value noise), c'est-à-dire que l'on interpole des valeurs aléatoires discrètes.

L'origine de l'erreur vient sans doute de ce site web, qui est le premier résultat lors d'une recherche google, et qui introduit un bruit semblable à celui que nous allons voir ensemble. Le fait est que, lorsque l'on parle de génération procédurale, et en particulier de bruit, il est impossible de passer à côté de Ken Perlin. Ses travaux dans le domaine sont précurseurs de bon nombre de techniques couramment utilisées aujourd'hui dans le monde de l'image de synthèse. Que l'algorithme que je présente ici soit de Perlin ou non, je n'en sais rien, bien que l'idée semble moins récente, mais dans tous les cas, il ne s'agit pas du bruit "de Perlin".

Alors il est vrai qu'on retrouve aujourd'hui couramment cette erreur, et on pourrait même dire que, pour beaucoup de gens qui ne sont pas spécialistes du sujet, le bruit de Perlin fait référence aussi bien à l'un qu'à l'autre. Cela étant, il me semble que conserver la même appellation pour deux choses différentes n'est pas une bonne idée, et je bannirai donc cette appellation de ce tutoriel.

Notez que j'ai moi-même fait l'erreur dans les précédentes versions de ce tutoriel.

Introduction aux bruits

En algorithmique, on utilise le terme de bruit pour parler d'un jeu de valeurs aléatoires erratiques.
Voici un exemple de bruit (sur une seule dimension) :

Petits mots de vocabulaire : l'amplitude du bruit est l'écart entre la plus grande et la plus petite valeur du bruit. On parle aussi d'effectif total pour désigner le nombre de valeurs du bruit.

La première étape consiste donc à générer un bruit. Il y a plusieurs méthodes possibles, et pour rester général, nous allons donc créer une fonction double bruit(int i) qui renvoie la i-ème valeur de notre bruit.
Pour générer ce bruit, il y a plusieurs méthodes :

• on génère un tableau préalablement à partir d'une fonction rand(), comme il y en a dans tous les langages ;

• on définit explicitement un jeu de valeurs dans le code source. C'est notamment nécessaire lorsqu'il est difficile de faire persister les données entre les différentes exécutions de l'algorithme (dans les shaders, par exemple) ;

• on peut aussi écrire sa propre fonction pseudo-aléatoire.

  Avec cette fonction vous aurez toujours, avec le même argument, le même bruit généré. Ajouter une constante, qui peut être modifiée entre deux générations successives de bruit, permet d'obtenir une grande variété de bruits. Notons que cette constante est généralement appelée graine (seed en anglais).

Si vous comptez utiliser une des deux premières méthodes, vous êtes maintenant en droit de vous demander quel est l'effectif total de notre bruit. J'en reparlerai plus tard, étant donné qu'il dépend d'un paramètre que nous n'avons pas encore abordé.

Par mesure de simplicité, disons que notre fonction bruit() renvoie un nombre entre 0 et 1.

Et si l'on voulait faire un bruit bidimensionnel ?

Dans ce cas, nous aurons besoin d'un bruit en deux dimension, c'est-à-dire d'un jeu de valeurs bivarié. On aurait donc une fonction similaire double bruit2D(int i, int j).

Maintenant on a notre bruit. Mais bon, ce n'est qu'un jeu de points discrets. L'étape suivante va être de construire une fonction assez régulière à partir de ça. À cette fonction, nous donnerons le nom de fonction de bruit.

Pour construire cette fonction, on associe à certaines de ses valeurs celles du bruit, et on tente de les relier entre eux. Pour associer les valeurs, on les dispose à intervalle régulier (on appelle la distance entre deux points le pas).
En clair, notre fonction est telle que fonction_bruit(i * pas) == bruit(i) (pour tout entier positif i plus petit que notre effectif total, évidemment).

Reprenons notre image précédente, et voici un résultat qu'on pourrait obtenir, avec un pas unitaire.

Maintenant, il ne reste plus qu'à calculer les points entre nos valeurs. Pour cela, nous allons introduire la notion d'interpolation.

Interpolation des valeurs discrètes

Ce procédé consiste à construire une fonction continue à partir d'un nombre fini de points. Il existe pas mal d'interpolations, mais nous n'en verrons ici que trois. Comme vous l'aurez compris, on cherche à relier les points entre eux. Pour cela, la stratégie habituelle consiste à utiliser une petite fonction différente pour chaque intervalle entre deux points (c'est le cas de toutes les interpolations utilisées ici, mais ce n'est pas une généralité — les interpolations polynomiales, par exemple, ne sont pas construites ainsi). Évidemment, toutes ces petites fonctions se ressembleront fortement, c'est-à-dire qu'elles seront de même nature, mais leurs coefficients varieront.

Lorsque ces « petites » fonctions sont des polynômes de degré fixé, on parle de spline. Par exemple, l'interpolation linéaire porte aussi le nom de spline linéaire, et c'est aussi le cas de la cubique.

Interpolation linéaire

Cette interpolation est simplissime, elle relie simplement deux points par un segment.

La fonction suivante interpole les points a et b, avec x le facteur qui varie entre 0 et 1 (0 correspond évidemment à a et 1 à b).

double interpolation_lineaire(double a, double b, double x) { return a * (1 - x) + b * x; }

Le seul avantage que l'on puisse lui trouver est la rapidité d'exécution, mais le rendu final aura vite fait de nous en dégoûter.

Pour faire une interpolation linéaire en deux dimensions des points formant le carré $abcd$, il suffit d'interpoler a et b, et c et d sur x, puis les résultats obtenus sur y.

double interpolation_lineaire2D(double a, double b, double c, double d, double x, double y) { i1 = interpolation_lineaire(a, b, x); i2 = interpolation_lineaire(c, d, x); return interpolation_lineaire(i1, i2, y); } L'interpolation cosinusoïdale

Le principe est le même que pour la première : on a un point a, un point b, et on relie les deux. Sauf que, au lieu de les relier linéairement, on va utiliser une courbe en forme de S. Jetons donc un œil à la fonction $\frac{(1 - \cos(\pi x))}2$.

Comme on le voit sur le graphe de la fonction, elle varie entre 0 et 1 lorsque x est entre 0 et 1. Ainsi, on peut considérer $\frac{(1 - \cos(\pi x))}2$ plutôt que x comme facteur d'interpolation.

double interpolation_cos(double a, double b, double x) { double k = (1 - cos(x * PI)) / 2; return interpolation_lineaire(a, b, k); }

Pour la faire en 2D, le principe est le même qu'au-dessus :

double interpolation_cos2D(double a, double b, double c, double d, double x, double y) { double x1 = interpolation_cos(a, b, x); double x2 = interpolation_cos(c, d, x); return interpolation_cos(x1, x2, y); } L'interpolation cubique

L'interpolation cubique consiste à relier deux points en calculant une cubique les reliant.

Une cubique est un polynôme du troisième degré, c'est-à-dire une fonction du type $f(x) = ax^3 + bx^2 + cx + d$ où a, b, c, d sont des constantes (avec d non nul de préférence).

Bon comme vous l'avez peut-être remarqué, il existe une infinité de cubiques passant par ces deux points, et elles ne donnent pas toutes des résultats très glorieux. Pour réduire le nombre de cubiques possibles, nous allons aussi regarder le point précédent ainsi que le suivant. Il y a d'autres critères concernant la continuité des dérivées premières et secondes, mais cela sort du cadre de ce tutoriel.

double interpolation_cubique(double y0, double y1, double y2, double y3, double x) { a = y3 - y2 - y0 + y1; b = y0 - y1 - a; c = y2 - y0; d = y1; return a *x * x * x + b * x * x + c * x + d; }

Eh bien, passons à la 2D, nous verrons... On va faire comme pour les autres alors.
Alors le hic, c'est qu'on a 4 points, non plus 2.
Alors on n'aura pas 3 interpolations de 2 points, mais bien 5 de 4 points.

double interpolation_cubique2D(double y00, double y01, double y02, double y03, double y10, double y11, double y12, double y13, double y20, double y21, double y22, double y23, double y30, double y31, double y32, double y33, double x, double y) v0 = interpolation_cubique(y00,y01,y02,y03,x); v1 = interpolation_cubique(y10,y11,y12,y13,x); v2 = interpolation_cubique(y20,y21,y22,y23,x); v3 = interpolation_cubique(y30,y31,y32,y33,x); return interpolation_cubique(v0,v1,v2,v3,y);

Voilà qui est fini pour les interpolations. Il reste maintenant à construire notre fonction.

Pour les matheux en herbe :

Il est intéressant ici de remarquer le sens mathématique de ce que l'on appelait tantôt la « régularité » de la fonction. En fait, une fonction est plus régulière plus elle est continument dérivable un grand nombre de fois (donc plus le $n$ tel qu'elle appartienne à $C_n$ est grand). Notons que comme on n'a rien de plus méchant que des polynômes et des cosinus, et que ces fonctions sont infiniment continument dérivables, le seul problème qu'on pourrait avoir se situe aux « joints » entre les différentes fonctions.

Regardons la régularité des interpolations proposées :

• interpolation linéaire $\in C_0$ ;

• interpolation cosinusoïdale $\in C_1$ ;

• interpolation cubique $\in C_2$.

Cela appuie donc notre intuition.

Notons que ce n'est plus nécessairement vrai en 2D...

Interpolons notre bruit

Nous avons maintenant une fonction qui nous permet de relier deux points successifs. L'étape suivante consiste donc à obtenir une fonction qui, pour une valeur donnée, trouve les deux points entre lesquels se situe la valeur, ainsi que la position entre les deux points (sous la forme d'une pourcentage, entre 0 et 1, où 0 est le point à gauche et 1 le point à droite) et les relies grâce à la fonction précédemment définie.

Comme les différents points sont espacés d'une distance pas, le point à gauche de notre valeur (appelons la x) à pour indice un entier i tel que i * pas <= x < (i+1) * pas (et donc i+1 est l'indice du point à droite).

On trouve donc l'indice du point à gauche en calculant le plus grand entier plus petit ou égal à x / pas (en utilisant floor(), ou en faisant un cast en entier).

La position entre les deux points se calcule en prenant distance entre la valeur et le point à gauche moins la distance entre les deux points, c'est-à-dire (x - i*pas) / pas. Cela dit, comme c'est équivalant à x / pas - i, et que i est la partie entière de x/pas (c'est comme ça qu'on a calculé i), on peut aussi prendre la partie décimale de x / pas (donc on fait modulo 1).

Évidement, dans le cas d'une interpolation cubique, il faut trouver plus de points, mais ce n'est guère différent. Voilà le code de la fonction :

// int pas = ... // int effectif = ... double fonction_bruit(double x) { int i = (int) (x / pas); return interpolation_cos(bruit(i), bruit(i + 1), (x / pas) % 1); // On peut aussi choisir une des deux autres interpolations : // return interpolation_lineaire(bruit(i), bruit(i + 1), (((double) x) / pas) % 1); // return interpolation_cubique(bruit(max(0, i - 1)), bruit(i), bruit(i + 1), bruit(min(i + 2, effectif)), (x*f)%T); }

Remarquez pour l'interpolation cubique, l'utilisation de min() et de max() pour éviter les erreurs de débordement si vous utilisez des tableaux.

Et en 2D, rien de bien étonnant :

// int pas = ... double fonction_bruit2D(double x, double y) { int i = (int) (x / pas); int j = (int) (y / pas); return interpolation_cos2D(bruit2D(i, j), bruit2D(i + 1, j), bruit2D(i, j + 1), bruit2D(i + 1, j + 1), (x / pas) % 1, (y / pas) % 1); // Évidemment, on...