Texture et bruit de Voronoï

Alors, Texture de Voronoï... Qu'est-ce que c'est?

En général ça se présente sous la forme d'une heightmap construite sur la base d'un diagramme de Voronoï.

Il s’agit d’un objet mathématique découpant l’espace en cellules de manière particulière. On parle de pavage de Voronoï, ou de Voronoï tesselation en anglais.

C

Cet objet mathématique est définit ainsi :

Soit le plan affine R² et soit S, un ensemble fini de points du plan.

  • - On appelle germes, graines (seed, en anglais) ou bien encore sites les points de cet ensemble S.
  • - On appelle région ou cellule (cell, en anglais) de Voronoï, associée à un élément p de S, l’ensemble de TOUS les points qui sont plus proches de p que de tout autre point de S. Ici, les éléments de S sont des points de références et on évalue les distances entre ceux là et tous les points qui existent sur le plan. Parmi tous les points d’une cellule de Voronoï V, il ne peut y avoir qu’un seul élément de S : Celui qui a servi à construire V.

À titre de rappel un plan affine est un plan dans lequel il n’y aurait pas de repère. En effet dans ce plan nous ne nous intéressons pas aux positions des points mais à la distance de ceux là entre eux.

Mathématiquement, cela se traduit ainsi :

mathVoronoi.png

Ce qui ce lit L'ensemble des points x appartenant à R² tel que pour tout point q appartenant à S la distance entre le point x et p soit inférieure ou égale à la distance entre le point x et q.

  • - Avec x, un point quelconque du plan R².
  • - Avec p, le point de référence ∈ S.
  • - Avec q, un point quelconque ∈ S.

Enfin, on peut ajouter qu'un point donné du plan qui est à égale distance de deux germes appartient à une frontière de Voronoï si ce point donné n'est pas à une distance moindre d'un autre germe.

Cela étant dit, en considérant les propriétés du diagramme de Voronoï, il est possible de construire une heightmap où la hauteur de chaque pixel serait déterminée selon sa distance avec le site le plus proche. Ce qu'on va faire pour commencer c'est se munir d'une graine suffisamment grande pour contenir tous nos sites. En effet, comme pour le bruit de perlin, ou le stamp noise, nous allons générer plusieurs octaves de notre bruit. Chacune des octaves aurait son propre ensemble de sites.

Tous nos sites sont stockés dans un tableau arbitrairement long de 4096 octects.

  	auto time = std::chrono::high_resolution_clock::now();
        auto count = time.time_since_epoch().count();

	std::default_random_engine generator;
	std::uniform_int_distribution<int> distribution(0,4294967295);
        generator.seed(std::default_random_engine::result_type(count));
		
	auto dice = std::bind ( distribution, generator );
		
	unsigned char * seed = new unsigned char[4096];
	int limit = 4096 / sizeof(unsigned int)
	for (int i = 0; i < limit; i++) {
		((unsigned int*) seed)[i] = (unsigned int) dice();
	}

On utilise les librairies standard chrono et random pour générer nos graines. À partir de là les choses restent relativement simples. On parcourt notre surface en x et en y

	surface[x + y * scale] += persistence * Voronoi(x, y, seeds);

et pour chaque pixel, on incrémente avec la valeur de la fonction Voronoi multipliée par la persistance.À titre de rappel, la persistance est en fait un coefficient définit par l'octave courante. Plus le grain du bruit est fin, plus la persistance est faible. Dans la pratique, et c'était déjà le cas dans le bruit de Perlin, ça permet de contrôler l'intensité maximale de chaque octave.

La fonction Voronoi est implémentée rigoureusement telle que décrit par la définition mathématique:

GLfloat VoronoiCell(int x, int y, unsigned char * seeds) {
	GLfloat tmp, dx, dy;
	
	GLfloat minDist = 16581375.0;
	for (int i = 0; i < seedsNumber; i++) {
		dx = (x - (((int*)seeds)[(seedOffset+i*2) & limit  ] & modScale));
		dy = (y - (((int*)seeds)[(seedOffset+i*2+1) & limit] & modScale));
		tmp = dx*dx+dy*dy;
		if (tmp < minDist) {
			minDist = tmp;
		}
	}

	return sqrt(minDist);
}

Pour un point donné on cherche le site le plus proche. Une fois ce site trouvé, on renvoie la distance associée.

Dans mon implémentation très simple seedOffset désigne la position de départ de l'ensemble des sites stockés à l’intérieur de seeds, le fameux tableau de 4096 float.

La variable limit est en fait une constante fixée a 4095. Pour éviter de faire un modulo sur une puissance de deux (4096), on peut utiliser l'opérateur "ET". On dirait pas comme ça, mais ça fait gagner approximativement 120 millisecondes!

Vous remarquerez que tout ça nous oblige à utiliser la fonction sqrt qui n'est pas connue pour être particulièrement rapide. En fait le programme est très lent. Il n'est pas particulièrement compliqué, mais il est très très lent. Quand il est implémenté côté cpu par exemple, pour une octave contenant 32 sites, la génération de la surface s'exécute en moyenne en 70 millisecondes sur un core i7. C'est pas terrible. C'est même très mauvais si on veut faire du temps réel. Néanmoins on peut substantiellement augmenter les performances en implémentant l'algorithme côté GPU. En effet il est très facile de paralléliser cet algorithme dans un compute shader:

#version 430

layout (local_size_x = 1, local_size_y = 1) in;
layout (rgba32f, binding = 0) uniform image2D voronoi;

uniform uint brushScale;
uniform vec4 voronoiSeeds[64];\n"

const uint area = brushScale*brushScale;

void main() {
	vec2 pixel = vec2(gl_GlobalInvocationID.xy);
	float minDistance = 16581375.0;
	vec2 d
	float tmp;
	for (int i=0; i < 64; i++) {
		d.x = pixel.x - voronoiSeeds[i].x;
		d.y = pixel.y - voronoiSeeds[i].y;
		tmp = d.x*d.x + d.y * d.y;

		if (tmp < minDistance) {
			minDistance = tmp;
		}
	}

	imageStore(voronoi, ivec2(gl_GlobalInvocationID.xy), vec4(vec3(sqrt(minDistance)),1.0));
}

Si vous ne connaissez pas OpenGL et en particulier le fonctionnement de ce qu'on appelle les shaders, il peut-être bon d'expliquer un peu comment ça marche. Un shader est un programme informatique s'exécutant sur la carte graphique de manière parallèle. Le langage de programmation utilisé pour écrire un shader OpenGL est GLSL, il est très proche du langage C.

Dans notre shader ci-dessus, nous avons une fonction main qui sera donc invoquée un certain nombre de fois en parallèle. Ce nombre est déterminé soit par l'utilisateur, soit par la pipeline de rendu. Ce qui est intéressant, c'est que la carte graphique peut manipuler plusieurs centaines d'invocations de shader en même temps. Inutile de vous dire que ça va très très vite et que c'est très pratique pour certains type d'algorithme. Dans notre cas, ça marche vachement bien parce que le programme est facilement parallélisable. On obtient de jolies performances, en effet le shader s'exécute en environ 15 à 20 millisecondes pour un diagramme de Voronoï de 32 sites pour une octave.

Dès qu'on augmente la profondeur, et donc le nombre d'octave, ça devient tout de suite plus lent, évidement, même sur GPU. En pratique, pour U N I V E R S E je ne pense pas utiliser directement ce type d'algorithme. Par contre il peut-être intéressant de combiner cette approche avec le Stamp Noise. C'est à dire que l'on générerait une fois une texture de Voronoï, pour une seule octave, qu'on réutiliserait comme on le faisait avec le cône de base utilisé dans l'implémentation original du Stamp noise.

Bon du coup, visuellement ça donne ça avec la fonction de Voronoi original.

voronoi32seeds1octave.png

voronoi32seeds1octaveBumped.png

Il est possible de sensiblement changer la fonction de Voronoi pour obtenir un résultat très différent, ainsi la fonction devient:

GLfloat VoronoiCrack(int x, int y, unsigned char * seeds) {
	GLfloat tmp, dx, dy;
	
	GLfloat minDist1 = 16581375.0;
	GLfloat minDist2 = 16581375.0;
	for (int i = 0; i < seedsNumber; i++) {
		dx = (x - (((int*)seeds)[(seedOffset+i*2) & limit  ] & modScale));
		dy = (y - (((int*)seeds)[(seedOffset+i*2+1) & limit] & modScale));
		tmp = dx*dx+dy*dy;
		if (tmp < minDist2) {
			if (tmp < minDist1) {
			  minDist2 = minDist1;
			  minDist1 = tmp;
			}
			else {
			  minDist2 = tmp;
			}
		}
	}

	return - sqrt(minDist2) + sqrt(minDist1);
}

Ici, ce qu'on fait, c'est qu'on cherche le second site le plus proche du point courant. On renvoie ensuite la différence. Voilà le résultat, toujours pour une octave et 32 sites:

voronoiCrack32seeds1octave.png

voronoiCrack32seeds1octaveBumped.png

Alors attention là les performances sont un peu réduites. Le temps d'exécution monte à 90 millisecondes sur un core i7.

Une autre variante consiste cette fois à ne plus faire la différence, mais à retourner le produit des deux distances, le temps d'exécution est toujours très élevé est reste approximativement le même que précédemment.

GLfloat VoronoiJelly(int x, int y, unsigned char * seeds) {
	GLfloat tmp, dx, dy;
	
	GLfloat minDist1 = 16581375.0;
	GLfloat minDist2 = 16581375.0;
	for (int i = 0; i < seedsNumber; i++) {
		dx = (x - (((int*)seeds)[(seedOffset+i*2) & limit  ] & modScale));
		dy = (y - (((int*)seeds)[(seedOffset+i*2+1) & limit] & modScale));
		tmp = dx*dx+dy*dy;
		if (tmp < minDist2) {
			if (tmp < minDist1) {
			  minDist2 = minDist1;
			  minDist1 = tmp;
			}
			else {
			  minDist2 = tmp;
			}
		}
	}

	return -sqrt(minDist2) * sqrt(minDist1);
}

voronoiJelly32seeds1octave.png

voronoiJelly32seeds1octaveBumped.png

Comme je le disais, en terme de performance, ce n'est pas folichon, mais visuellement c'est intéressant. Si on utilise des diagrammes de Voronoï avec le stamp noise, il n'y a pas besoin de générer plusieurs octaves sur ce type de bruit. La preuve de concept associée à cette article permet cependant de le faire. Puisqu'on est là, autant voir ce que ça donne!

La boucle de génération devient donc

	GLfloat persistence;

	for (int octave = 1; octave <= depth; octave++) {
	  	seedsNumber *= 4;
		persistence = (1.0 / (GLfloat) (1 << octave));
		for (int x = 0; x < scale; x++) {
			for (int y = 0; y < scale; y++) {
			  	#ifdef VORONOI_CRACK
					surface[x + y * scale] += persistence * VoronoiCrack(x, y, seeds);
				#else
			  		#ifdef VORONOI_JELLY
					surface[x + y * scale] += persistence * VoronoiJelly(x, y, seeds);
					#else
					surface[x + y * scale] += persistence * VoronoiCell(x, y, seeds);
					#endif
				#endif
				if(surface[x + y * scale] < min) {
					min = surface[x + y * scale];
				}
				if(surface[x + y * scale] > max) {
					max = surface[x + y * scale];
				}
			}
		}
		seedOffset += seedsNumber;
	}

Comme les temps de génération augmente vraiment VRAIMENT vite on se limite à quatre octaves. Également, pour un rendu plus intéressant me semble-t-il, j'augmente à chaque octave d'un facteur quatre le nombre de sites. Avec ça on arrive approximativement à 5 secondes sur un core i7 pour chacune des variations de la fonction de voronoï. Visuellement on obtient ceci:

voronoiCell32seeds4octave.png

voronoiCell32seeds4octaveBumped.png

voronoiCrack32seeds4octave.png

voronoiCrack32seeds4octaveBumped.png

voronoiJelly32seeds4octave.png

voronoiJelly32seeds4octaveBumped.png

C'est très esthétique, la combinaison de cette technique avec le Stamp Noise est très prometeuse, nous verrons ça dans un prochain article!

Pour essayer chez vous le programme, rendez-vous sur la page github du Poc.

Un grand merci à mes collègues de l'INRIA pour les conseils et les optimisations du code! Merci également à Laura pour sa précieuse relecture!