Tas de sable

A propos du code

Le code dans cet article a été écrit avec Code::Blocks et la SDL 2.
Vous pouvez trouver ici un guide pour installer ces logiciels.
Bien qu'il soit basé sur la SDL, je n'utilise pas ses fonctions directement. J'ai écrit une petite bibliothèque
avec quelques fonctions basiques pour faciliter la compréhension et la portabilité dans un autre langage.
Vous pouvez en apprendre plus sur cette bibliothèque ici.

Règles des tas de sable

L'idée pour cet article vient d'une vidéo de Numberphile.
Cette vidéo décrit un objet mathématique appelé tas de sable.
La première partie de la vidéo définit des opérations arithmétique sur cet objet.
Mais à la fin elle montre une façon de créer des images fractales que je vais essayer de reproduire.

Alors voyons ce qu'est un tas de sable et comment il évolue.
Considérez une grille carrée. Chaque case de cette grille contient un nombre entre 0 et 3.
Si a un moment une case contient une valeur plus grande que 3, une règle "d'effondrement" s'applique. La valeur de
cette case est réduite de 4, et la valeur de ses 4 voisines directes est augmentée de 1.


Et cette règle est répétée jusqu'à ce qu'il n'y ait plus de case avec une valeur supérieure à 3.
Pour visualiser les tas de sable on peut assigner une couleur à chaque case en fonction de sa valeur.

Maintenant, pour produire une fractale, on a simplement besoin d'une grille assez grande, on met une grande valeur
dans sa case centrale, et on répète la règle d'effondrement tant qu'on le peut.

Premier essai

Comme la grille doit être assez grande, on affichera un pixel par case.
Comme la règle d'effondrement modifie les cases voisines, dans le code on stockera la grille dans 2 tableaux
d'entiers.
L'un contiendra la grille courante qu'on est en train de traiter et l'autre contiendra l'état précédent.
De cette façon on peut appliquer la règle sur chaque case sans avoir à se soucier des effets de bord.
		static int table[2][SCREEN_WIDTH][SCREEN_HEIGHT];
		int currentTable = 0;
				
On initialise un des tableaux avec des 0, sauf pour la case centrale où on met une grande valeur:
		// initialise the table
		for (int y = 0; y < SCREEN_HEIGHT; ++y)
			for (int x = 0; x < SCREEN_WIDTH; ++x)
				table[currentTable][x][y] = 0;

		table[currentTable][SCREEN_WIDTH / 2][SCREEN_HEIGHT / 2] = 400000;
				
Pour afficher les cases on va avoir besoin d'une palette qui assigne une couleur à chaque valeur de case possible.
J'ai essayé de reproduire la palette qui apparait à 22:46 dans la vidéo:
- 0 est bleu
- 1 est cyan
- 2 est jaune
- 3 est rouge
- et j'ai ajouté une couleur blanche pour les valeurs plus grandes que 3
		// init palette
		Color palette[5];
		palette[0] = Color(0, 0, 255);
		palette[1] = Color(0, 255, 255);
		palette[2] = Color(255, 255, 0);
		palette[3] = Color(255, 0, 0);
		palette[4] = Color(255, 255, 255);
				
La partie affichage est simple, on retrouve simplement la bonne couleur dans la palette pour chaque case:
		// draw the current table
		for (int y = 0; y < SCREEN_HEIGHT; ++y)
			for (int x = 0; x < SCREEN_WIDTH; ++x)
			{
				int color = table[currentTable][x][y];

				if (color > 4)
					color = 4;
				gfx.setPixel(x, y, palette[color]);
			}

		gfx.render();
		sys.processEvents();
				
La partie traitement est assez intuitive aussi. D'abord on copie l'ancienne grille dans la courante.
		// update the table
		for (int y = 0; y < SCREEN_HEIGHT; ++y)
			for (int x = 0; x < SCREEN_WIDTH; ++x)
			{
				table[1 - currentTable][x][y] = table[currentTable][x][y];
			}
				
Puis on boucle sur chaque case pour appliquer la règle d'effondrement si on trouve une valeur supérieur à 3.
Remarquez qu'on teste aussi pour chaque case "voisine" si elle est en dehors de la grille.
		for (int y = 0; y < SCREEN_HEIGHT; ++y)
			for (int x = 0; x < SCREEN_WIDTH; ++x)
			{
				if (table[currentTable][x][y] >= 4)
				{
					table[1 - currentTable][x][y] -= 4;

					if (x - 1 >= 0)
						table[1 - currentTable][x - 1][y] += 1;

					if (x + 1 < SCREEN_WIDTH)
						table[1 - currentTable][x + 1][y] += 1;

					if (y - 1 >= 0)
						table[1 - currentTable][x][y - 1] += 1;

					if (y + 1 < SCREEN_HEIGHT)
						table[1 - currentTable][x][y + 1] += 1;
				}
			}
				
Enfin on échange les grilles.
		currentTable = 1 - currentTable;

		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Voici le résultat de ce programme:


Le seul problème c'est que sur mon PC actuel avec un i3-4330, ce code prend plus de 30 minutes avant de se terminer.
Alors on va optimiser ça.

Accélération (1)

Dans notre boucle de traitement, on doit tester chaque case voisine pour voir si elle est à l'intérieur de la
grille.
Ces tests coûtent cher mais il y a une méthode simple pour s'en passer.
On va simplement étendre la grille d'une case dans chaque direction et on ne traitera et n'affichera que la partie
centrale.
		static int table[2][SCREEN_WIDTH + 2][SCREEN_HEIGHT + 2];

		[...]

		for (int y = 1; y <= SCREEN_HEIGHT; ++y)
			for (int x = 1; x <= SCREEN_WIDTH; ++x)
			{
				if (table[currentTable][x][y] >= 4)
				{
					table[1 - currentTable][x][y] -= 4;
					table[1 - currentTable][x - 1][y] += 1;
					table[1 - currentTable][x + 1][y] += 1;
					table[1 - currentTable][x][y - 1] += 1;
					table[1 - currentTable][x][y + 1] += 1;
				}
			}
				
Et au lieu d'utiliser une boucle pour copier l'ancienne grille dans la courante avant le traitement, on peut
utiliser la fonction memcpy qui est plus efficace:
		int screenSize = (SCREEN_WIDTH + 2) * (SCREEN_HEIGHT + 2);
		memcpy(table[1 - currentTable], table[currentTable], screenSize * sizeof(int));
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Ce code prend environ 26 minutes et 40s sur mon ordinateur.
C'est plus rapide, mais on peut faire mieux.

Accélération (2)

Redessiner tout l'écran à chaque frame est coûteux parce que ma bibliothèque n'est pas particulièrement optimisée
pour l'affichage.
Et pour ce programme on n'a pas vraiment besoin de voir chaque étape du dessin. Le plus important c'est le résultat
final.
On peut facilement modifier le code pour n'afficher la grille qu'une fois tous les 16 traitements:
		int counter = 0;

		while (sys.isQuitRequested() == false)
		{
			[...]

			counter++;
			if ((counter % 16) == 0)
				gfx.render();
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Ce code s'exécute en 18 minutes et 30s.

Accélération (3)

Maintenant observons les opérations de la règle d'effondrement.
Si une case est supérieure ou égale à 4, on lui soustrait 4 et on propage 1 dans chaque direction.
Mais si après lui avoir soustrait 4 la valeur de la case est toujours supérieure ou égale à 4, on doit attendre la
boucle de traitement suivante pour continuer.
Ça serait plus efficace si on répétait la règle d'effondrement jusqu'à ce que la valeur passe en-dessous de 4.
Notre boucle de traitement ressemblerait à ça:
		for (int y = 1; y <= SCREEN_HEIGHT; ++y)
			for (int x = 1; x <= SCREEN_WIDTH; ++x)
			{
				while (table[currentTable][x][y] >= 4)
				{
					table[1 - currentTable][x][y] -= 4;
					table[1 - currentTable][x - 1][y] += 1;
					table[1 - currentTable][x + 1][y] += 1;
					table[1 - currentTable][x][y - 1] += 1;
					table[1 - currentTable][x][y + 1] += 1;
					table[currentTable][x][y] -= 4;
				}
			}
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Remarquez qu'on utilise la case de l'ancienne grille comme un compteur. On la décrémente et on la teste pour sortir
de la boucle.
Cette version prend seulement 2min et 35s à s'exécuter.

Accélération (4)

Le programme tourne beaucoup plus vite maintenant, mais on peut encore aller plus loin en faisant un petit
compromis.

D'abord on peut supprimer la boucle qu'on a introduit dans la dernière partie.
En effet, si on observe cette boucle, on réduit la valeur de la case de 4 à chaque fois.
Donc on peut prédire combien de fois on va boucler: c'est le résultat de la division par 4 de la valeur initiale.
Et on peut réécrire notre boucle de traitement comme ceci:
		for (int y = 1; y <= SCREEN_HEIGHT; ++y)
			for (int x = 1; x <= SCREEN_WIDTH; ++x)
			{
				int height = table[currentTable][x][y] / 4;
				table[1 - currentTable][x][y] -= 4 * height;
				table[1 - currentTable][x - 1][y] += height;
				table[1 - currentTable][x + 1][y] += height;
				table[1 - currentTable][x][y - 1] += height;
				table[1 - currentTable][x][y + 1] += height;
			}
				
Ensuite, on a vu qu'on utilisait l'ancienne grille comme un compteur.
Alors on pourrait probablement n'utiliser qu'une seule grille au lieu de deux.
Ça éviterait l'opération de memcpy.

Il y a toujours un danger d'avoir un effet de bord quand on fait ce genre de chose, parce que la règle
d'effondrement va modifier la prochaine case qu'on traitera.
Mais cette règle est simple, et même si les opérations ne se font pas dans le même ordre, tout se résoudra avec les
boucles suivantes.

Le principal compromis qu'on fait quand on n'utilise qu'une grille c'est que quand le programme s'exécutera,
l'affichage n'aura pas l'air symétrique.
Alors voici le code complet quand on n'utilise qu'une seule grille:
		#include <stdio.h>
		#include <stdlib.h>
		#include "main.h"
		#include "Graphics.h"
		#include "System.h"

		#define SCREEN_WIDTH    600
		#define SCREEN_HEIGHT   600

		int main(int argc, char* argv[])
		{
			// init the window
			gfx.init("Sandpiles 5", SCREEN_WIDTH, SCREEN_HEIGHT);
			gfx.init2D();
			gfx.clearScreen(Color(0, 0, 0, SDL_ALPHA_OPAQUE));

			static int table[SCREEN_WIDTH + 2][SCREEN_HEIGHT + 2];

			// initialise the table
			for (int y = 1; y <= SCREEN_HEIGHT; ++y)
				for (int x = 1; x <= SCREEN_WIDTH; ++x)
					table[x][y] = 0;

			table[SCREEN_WIDTH / 2][SCREEN_HEIGHT / 2] = 400000;

			// init palette
			Color palette[5];
			palette[0] = Color(0, 0, 255);
			palette[1] = Color(0, 255, 255);
			palette[2] = Color(255, 255, 0);
			palette[3] = Color(255, 0, 0);
			palette[4] = Color(255, 255, 255);

			int counter = 0;

			while (sys.isQuitRequested() == false)
			{
				// draw the current table
				for (int y = 0; y < SCREEN_HEIGHT; ++y)
					for (int x = 0; x < SCREEN_WIDTH; ++x)
					{
						int color = table[x + 1][y + 1];

						if (color > 4)
							color = 4;
						gfx.setPixel(x, y, palette[color]);
					}

				counter++;
				if ((counter % 16) == 0)
					gfx.render();

				sys.processEvents();

				// update the table
				int screenSize = (SCREEN_WIDTH + 2) * (SCREEN_HEIGHT + 2);

				for (int y = 1; y <= SCREEN_HEIGHT; ++y)
					for (int x = 1; x <= SCREEN_WIDTH; ++x)
					{
						int height = table[x][y] / 4;
						table[x][y] -= height * 4;
						table[x - 1][y] += height;
						table[x + 1][y] += height;
						table[x][y - 1] += height;
						table[x][y + 1] += height;
					}
			}

			gfx.quit();

			return EXIT_SUCCESS;
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Et ce programme ne prend que 1min 33s pour dessiner l'image finale. C'est 20 fois plus rapide que le premier
programme.
Et bien sur on obtient la même image qu'avec le premier.

Plusieurs points de départ

Maintenant qu'on a un moyen rapide d'afficher notre fractale, on peut essayer de jouer avec pour voir si on peut
changer sa forme.
D'abord voyons ce qui se passe si on démarre avec 3 cases de grande valeur au lieu d'une.
Est-ce qu'elles vont entrer en collision et changer le forme de la fractale finale ?

Pour faire ça on a seulement besoin de changer l'initialisation de notre programme:
		int cx = SCREEN_WIDTH / 2 + 1;
		int cy = SCREEN_HEIGHT / 2 + 1;
		table[cx][cy-70] = 200000;
		table[cx-50][cy+50] = 150000;
		table[cx+50][cy+50] = 150000;

		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Voilà une image de ce qu'on voit au début de l'exécution:


Et voilà ce qu'on obtient à la fin.


Même si la partie chaotique au centre n'a pas la même forme, l'anneau "ordonné" autour n'a pas l'air très différent.

Essayons avec un carré

Que va-t-il se passer si on essaye de dessiner un carré au départ au lieu de quelques points ?
Si on modifie l'initialisation comme ça::
		int cx = SCREEN_WIDTH / 2 + 1;
		int cy = SCREEN_HEIGHT / 2 + 1;
		for (int i = -50; i <= 49; ++i)
		{
			table[cx - 50][cy + i] = 2500;
			table[cx + 49][cy + i] = 2500;
			table[cx + i][cy - 50] = 2500;
			table[cx + i][cy + 49] = 2500;
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Voilà ce qu'on voit au début:


Et ce qu'on obtient à la fin.


Encore une fois la forme ne semble pas très différente.

Un espace initial non vide

Au lieu de démarrer avec une grille vide, essayons d'en remplir la moitié avec des cases à 1.
		// initialise the table
		for (int y = 1; y <= SCREEN_HEIGHT; ++y)
			for (int x = 1; x <= SCREEN_WIDTH; ++x)
				table[x][y] = (y > x ? 1 : 0);

		int cx = SCREEN_WIDTH / 2;
		int cy = SCREEN_HEIGHT / 2;
		table[cx][cy] = 400000;

		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Voici ce que ça donne:


Ça c'est intéressant. La partie inférieure a l'air étirée comme si on l'observait à travers une loupe.
Vous pouvez aussi remarquer que sa forme n'est plus vraiment circulaire. Elle ressemble plus à un octogone.

Motif

Maintenant qu'est-ce qu'on obtient si la moitié de la grille n'est plus remplie par une valeur unie, mais plutôt
par un motif ?
Essayons de l'initialiser avec un motif de 1 et de 2 alternés.
		// initialise the table
		for (int y = 1; y <= SCREEN_HEIGHT; ++y)
			for (int x = 1; x <= SCREEN_WIDTH; ++x)
				table[x][y] = (y > x ? (((y+x)&1) ? 2 : 1) : 0);

		int cx = SCREEN_WIDTH / 2 + 1;
		int cy = SCREEN_HEIGHT / 2 + 1;
		table[cx][cy] = 400000;

		Télécharger le code source
		Télécharger l'exécutable pour Windows
				
Voici l'image finale:


Le motif a l'air d'avoir "survécu" à l'opération et les triangles qu'on obtient on l'air d'être remplis avec
d'autres types de motifs.

La vidéo d'origine montrait aussi des images créées avec d'autres règles. Par exemples ils autorisaient les cases à
prendre des valeurs plus grandes que 3.
Ça ne devrait pas être difficile de modifier ce code pour s'adapter à ces nouvelles règles.
J'espère que vous vous amuserez à explorer ce monde de fractales cools.

Links

Vidéo de quelques programmes de cet article

La vidéo originale de Numberphile

Vidéo: Tas de sable en processing