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!