Réaction-Diffusion

Note pour les utilisateurs d'Internet Explorer

Cette page affiche des formules mathématiques en utilisant un script MathML.
Si vous lisez ceci avec Internet Explorer, je vous recommande d'accepter l'exécution de ce script quand IE vous le
demande.

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.

Définition

Un système à réaction-diffusion est une simulation de l'évolution d'une ou plusieurs substances chimiques.
Les systèmes à une seule substance ne sont pas très intéressants. Alors on va en simuler un avec 2 substances.

Alan Turing fut le premier à utiliser ces systèmes pour explique comment les motifs apparaissent sur la fourrure
d'un animal.
Mais de nombreux scientifiques ont poursuivi ses travaux et on découvert des équations plus simples.
Ici nous allons utiliser le modèle de Gray-Scott.
En particulier on va essayer de reproduire la vidéo de Coding Train et on va accélérer le code à un niveau
utilisable.

Donc pour simuler ça, chaque pixel de l'écran va contenir les concentrations de A et de B.
Ces valeurs vont rester entre 0 et 1.

A chaque pas de temps on va mettre à jour les concentration dans chaque pixel en utilisant ces formules:

At+1 = At + ( DA ∇2 ⁡ At - At Bt2 + f ( 1 - At ) ) Δt
Bt+1 = Bt + ( DB ∇2 ⁡ Bt + At Bt2 - ( k + f ) Bt ) Δt
Où:

Maintenant expliquons un peu tout ça. Il y a 3 parties à l'intérieur des parenthèses.
La 3ième partie signifie que la substance A va être introduite dans le système à un taux de "f", pendant que la
substance B sera retirée du système à un taux de "k".
La seconde partie "AB2" dit que 2 "atomes" de B vont réagir avec un "atome" de A et le convertir en un B.
La première partie avec le triangle bizarre est la partie "diffusion". Elle dit à quel vitesse les substances vont
se répandre sur les pixels voisins.
Le laplacien est un genre de dérivée en 2D. On peut l'approximer avec un produit de convolution.
Ca veut dire qu'on va prendre le pixel courant et ses 8 voisins, et on va les additionner avec des poids qui
suivront cette matrice:

Dans cet article, on va utiliser les mêmes valeurs que dans la vidéo:

Dans un autre article on verra ce qu'on peut obtenir si on change ces valeurs.

Programmation

Maintenant programmons ces équations.
On va stocker les 2 concentrations par pixel dans un tableau:
		float table[2][SCREEN_HEIGHT][SCREEN_WIDTH];
On utilisera un autre tableau pour stocker les résultats des calculs pour la frame suivante:
		float nextTable[2][SCREEN_HEIGHT][SCREEN_WIDTH];
Puis on définit nos constantes. 
		float DA = 1.0f;
		float DB = 0.5f;
		float f = 0.055f;
		float k = 0.062f;
On va avoir une fonction pour calculer le laplacien: 
		float Laplacian(int tab, int x, int y)
		{
Voici la matrice des poids: 
			static const float coef[3][3] =
			{
				{0.05f, 0.2f,  0.05f},
				{0.2f,  -1.0f, 0.2f},
				{0.05f, 0.2f,  0.05f}
			};
On démarre 2 boucles qui vont de -1 à 1, parce qu'on veut aller un pixel à gauche, un à droite, et un au-dessus/en
dessous.
			float sum = 0;

			for (int dy = -1; dy <= 1; dy++)
				for (int dx = -1; dx <= 1; dx++)
				{
On calcule les coordonnées du pixel, en le faisant boucler si on sort de l'écran.
					int posX = (x + dx + SCREEN_WIDTH) % SCREEN_WIDTH;
					int posY = (y + dy + SCREEN_HEIGHT) % SCREEN_HEIGHT;
On récupère le coefficient dans le tableau et on ajoute la valeur du pixel à la somme.
					float c = coef[dy + 1][dx + 1];
					sum += c * table[tab][posY][posX];
				}
Et finalement on retourne la somme.
			return sum;
		}
Une fonction de mise à jour va calculer les formules qu'on a vues ci-dessus.
		void update()
		{
On boucle sur chaque pixel:
			for (int y = 0; y < SCREEN_HEIGHT; ++y)
				for (int x = 0; x < SCREEN_WIDTH; ++x)
				{
On récupère les concentrations actuelles de A et B:
					float A = table[0][y][x];
					float B = table[1][y][x];
Et on applique les formules à nextTable:
					nextTable[0][y][x] = A + DA * Laplacian(0, x, y) - A * B * B + f * (1.0f - A);
					nextTable[1][y][x] = B + DB * Laplacian(1, x, y) + A * B * B - (k + f) * B;
				}
Finalement on copie le contenu de nextTable dans table pour être prêts pour la prochaine frame.
			memcpy(table, nextTable, SCREEN_WIDTH * SCREEN_HEIGHT * 2 * sizeof(float));
		}
Dans la fonction main(), on va initialiser la table pleine de A:
		// initialize the table
		for (int y = 0; y < SCREEN_HEIGHT; ++y)
			for (int x = 0; x < SCREEN_WIDTH; ++x)
			{
				table[0][y][x] = 1.0f;
				table[1][y][x] = 0.0f;
			}
Avec un petit carré de B au centre:
		for (int y = 0; y < 10; ++y)
			for (int x = 0; x < 10; ++x)
			{
				int posX = SCREEN_WIDTH  / 2 + x - 5;
				int posY = SCREEN_HEIGHT / 2 + y - 5;
				table[1][posY][posX] = 1.0f;
			}
Puis dans la boucle principale on va afficher le contenu de la table, en utilisant du rouge pour A et du bleu pour B:
		while (sys.isQuitRequested() == false)
		{
			// draw the current table
			for (int y = 0; y < SCREEN_HEIGHT; ++y)
				for (int x = 0; x < SCREEN_WIDTH; ++x)
				{
					Color col(table[0][y][x] * 255,
					          0,
					          table[1][y][x] * 255);
					gfx.setPixel(x, y, col);
				}
On gère les évènements système: 
			gfx.render();
			sys.processEvents();
Et finalement on appelle notre fonction de mise à jour: 
			// update the table
			update();
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Comme vous le verrez, ce programme est plutôt lent.
Si vous êtes patient, une fois qu'il a rempli l'écran vous verrez ça:


Chronométrer

Comme à chaque fois que vous voulez optimiser un programme, vous devez d'abord regarder le temps que prend votre
code.
On va utiliser les mêmes fonction de chronométrage qu'on a utilisées pour les metaballs.
Mais si vous essayez de chronométrer chaque frame séparément, vous verrez que le résultat augmente à chaque fois.
Alors on va plutôt additionner les temps des 2000 premières frames pour avoir une valeur plus stable.

On va utiliser 2 variables. "tim" va stocker la somme des temps et "cnt" sera un simple compteur du nombre de
frames.
		float tim = 0;
		int cnt = 0;
Au début de update() on va démarrer le timer: 
		sys.StartPerfCounter();
Et à la fin de cette fonction, on va incrémenter le compteur, additionner les temps et afficher le résultat à la
2000ième frame:
		cnt++;
		if (cnt < 2000)
			tim += sys.StopPerfCounter();
		else if (cnt == 2000)
			printf("%f\n", tim);

		Télécharger le code source
		Télécharger l'exécutable pour Windows
La valeur qu'on obtient est d'environ 45000 ms.

Utiliser des doubles

Quand j'optimise un code, d'habitude j'essaye de remplacer les floats par des ints.
C'est une habitude que j'ai acquis en travaillant sur de vieux ordinateurs et consoles de jeux.
Dans ce programme particulier on ne peut pas utiliser d'entiers. J'ai essayé avec différentes précisions et à
chaque fois la réaction s'est arrêtée rapidement.

Alors j'ai essayé d'utiliser des double à la place de floats et le résultat m'a étonné.
Voici le code quand vous remplacez tous les floats par des doubles, à part la variable du timer.
		#include <stdio.h>
		#include <stdlib.h>
		#include "main.h"
		#include "Graphics.h"
		#include "System.h"
		#include <string.h>
		#include <stdlib.h>

		#define SCREEN_WIDTH    400
		#define SCREEN_HEIGHT   400

		double table[2][SCREEN_HEIGHT][SCREEN_WIDTH];
		double nextTable[2][SCREEN_HEIGHT][SCREEN_WIDTH];

		double DA = 1.0;
		double DB = 0.5;
		double f = 0.055;
		double k = 0.062;

		double Laplacian(int tab, int x, int y)
		{
			static const double coef[3][3] =
			{
				{0.05, 0.2,  0.05},
				{0.2,  -1.0, 0.2},
				{0.05, 0.2,  0.05}
			};

			double sum = 0;

			for (int dy = -1; dy <= 1; dy++)
				for (int dx = -1; dx <= 1; dx++)
				{
					int posX = (x + dx + SCREEN_WIDTH) % SCREEN_WIDTH;
					int posY = (y + dy + SCREEN_HEIGHT) % SCREEN_HEIGHT;
					double c = coef[dy + 1][dx + 1];
					sum += c * table[tab][posY][posX];
				}

			return sum;
		}

		float tim = 0;
		int cnt = 0;
		void update()
		{
			sys.StartPerfCounter();
			for (int y = 0; y < SCREEN_HEIGHT; ++y)
				for (int x = 0; x < SCREEN_WIDTH; ++x)
				{
					double A = table[0][y][x];
					double B = table[1][y][x];

					nextTable[0][y][x] = A + DA * Laplacian(0, x, y) - A * B * B + f * (1.0 - A);
					nextTable[1][y][x] = B + DB * Laplacian(1, x, y) + A * B * B - (k + f) * B;
				}

			memcpy(table, nextTable, SCREEN_WIDTH * SCREEN_HEIGHT * 2 * sizeof(double));

			cnt++;
			if (cnt < 2000)
				tim += sys.StopPerfCounter();
			else if (cnt == 2000)
				printf("%f\n", tim);
		}

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

			// initialize the table
			for (int y = 0; y < SCREEN_HEIGHT; ++y)
				for (int x = 0; x < SCREEN_WIDTH; ++x)
				{
					table[0][y][x] = 1.0;
					table[1][y][x] = 0.0;
				}

			for (int y = 0; y < 10; ++y)
				for (int x = 0; x < 10; ++x)
				{
					int posX = SCREEN_WIDTH  / 2 + x - 5;
					int posY = SCREEN_HEIGHT / 2 + y - 5;
					table[1][posY][posX] = 1.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)
					{
						Color col(table[0][y][x] * 255,
						          0,
						          table[1][y][x] * 255);
						gfx.setPixel(x, y, col);
					}

				gfx.render();
				sys.processEvents();

				// update the table
				update();
			}

			gfx.quit();

			return EXIT_SUCCESS;
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
J'ai toujours pensé que les float étaient un peu plus rapides que les doubles en prenant en compte les accès
mémoire, mais je me trompais.
Ce code ne prend qu'environ 14000 ms !
Apparemment sur les processeurs PC modernes les doubles sont extrêmement optimisés.
J'ai comparé le code assembleur généré par le compilateur pour les 2 versions, mais je n'ai rien trouvé qui
expliquait ça.

Mais certaines astuces qui marchaient avec les floats ne marchent plus avec les doubles. C'est un tout nouveau
monde d'optimisation.
Maintenant je vais devoir revisiter certains autres programmes de mon site...

On a gagné beaucoup de temps, mais on peut faire encore mieux.

Réécrire le laplacien

La fonction Laplacian() prend probablement la plus grande partie du temps dans ce programme.
Et grâce à ses symétries on peut certainement l'écrire de manière plus efficace.
pour réduire les multiplications, on peut ajouter ensemble tous les pixels qui ont un coefficient de 0.05, puis
tous ceux qui ont un coefficient de 0.2.
Et finalement on va additionner toutes ces sommes partielles en utilisant leurs poids.

On va utiliser des tableaux pour précalculer les coordonnées et rendre les choses plus claires:
		double Laplacian(int tab, int x, int y)
		{
			int xp[3];
			int yp[3];

			for (int i = -1; i <= 1; i++)
			{
			   xp[i + 1] = (x + i +  SCREEN_WIDTH) % SCREEN_WIDTH;
			   yp[i + 1] = (y + i + SCREEN_HEIGHT) % SCREEN_HEIGHT;
			}
Puis on peut calculer nos sommes:
			double sum1 = table[tab][yp[0]][xp[0]] +
			              table[tab][yp[0]][xp[2]] +
			              table[tab][yp[2]][xp[0]] +
			              table[tab][yp[2]][xp[2]];

			double sum2 = table[tab][yp[0]][xp[1]] +
			              table[tab][yp[1]][xp[0]] +
			              table[tab][yp[1]][xp[2]] +
			              table[tab][yp[2]][xp[1]];

			return 0.05 * sum1 + 0.2 * sum2 - table[tab][yp[1]][xp[1]];
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Ce code prend seulement 9000 ms.

Eviter le wrapping

Les calculs qu'on fait pour faire boucler les coordonnées quand on est en dehors de l'écran utilisent un opérateur
modulo "%".
Mais cet opérateur prend autant de temps qu'une division. Et une division c'est très coûteux.

Alors on va supprimer ces calculs. Mais pour éviter de sortir de l'écran quand on calcule le laplacien on va
ajouter une cellule dans chaque direction à nos tables.
Alors leurs largeurs et leurs hauteurs seront augmentées de 2:
		double table[2][SCREEN_HEIGHT+2][SCREEN_WIDTH+2];
		double nextTable[2][SCREEN_HEIGHT+2][SCREEN_WIDTH+2];
Maintenant les coordonnées dans la laplacien sont simplement:
		for (int i = -1; i <= 1; i++)
		{
			xp[i + 1] = x + i;
			yp[i + 1] = y + i;
		}
On doit aussi changer les boucles de la fonction update():
		for (int y = 1; y < SCREEN_HEIGHT+1; ++y)
			for (int x = 1; x < SCREEN_WIDTH+1; ++x)
			{
				double A = table[0][y][x];
Modifier le memcpy():
		memcpy(table, nextTable, (SCREEN_WIDTH+2) * (SCREEN_HEIGHT+2) * 2 * sizeof(double));
Dans l'initialisation on doit aussi remplir les valeurs de nextTable, pour éviter des bordures noires sur l'image.
		// initialize the table
		for (int y = 0; y < SCREEN_HEIGHT+2; ++y)
			for (int x = 0; x < SCREEN_WIDTH+2; ++x)
			{
				table[0][y][x] = 1.0;
				table[1][y][x] = 0.0;
				nextTable[0][y][x] = 1.0;
				nextTable[1][y][x] = 0.0;
			}
Et bien sur la routine d'affichage change aussi:
		// draw the current table
		for (int y = 0; y < SCREEN_HEIGHT; ++y)
			for (int x = 0; x < SCREEN_WIDTH; ++x)
			{
				int r = table[0][y+1][x+1] * 255;
				int b = table[1][y+1][x+1] * 255;
				gfx.setPixel(x, y, Color(r, 0, b));
			}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Ce code n'a besoin que de 6000 ms pour les 2000 premières frames.
Bien sur on n'obtient pas la même image quand l'écran est rempli:


Mais on trouvera un moyen de réimplémenter le wrapping plus tard.

Réécrire le laplacien (2)

Il y a une façon plus intelligente de réécrire la fonction Laplacian().
Un produit de convolution peut être vu comme la multiplication d'une matrice par les pixels à l'écran.

( 0.05 0.2 0.05 0.2 -1 0.2 0.05 0.2 0.05 ) * pixels

Cette matrice a des symétries et elle pourrait presque être écrite comme le produit 2 matrices plus simples.
Essayons.

( 1 4 1 ) * ( 1 4 1 ) = ( 1 4 1 4 16 4 1 4 1 )

Maintenant si on multiplie ça par 0.05:

0.05 * ( 1 4 1 4 16 4 1 4 1 ) = ( 0.05 0.2 0.05 0.2 0.8 0.2 0.05 0.2 0.05 )

On obtient presque la matrice qu'on voulait à part pour l'élément central.
Mais comme on fait la somme de tous ces éléments, ça sera facile de soustraire 1.8 * table[t][y][x] à la fin pour
obtenir le coefficient "-1" qu'on veut.

Alors on peut décomposer le produit de convolution en 2 produits successifs.
Le fait qu'il y a une matrice horizontale et une verticale veut dire qu'une convolution est appliquée dans la
direction x alors que l'autre se fera le long de la direction y.

On va calculer ces opérations pour chaque pixel de l'écran et les stocker dans 2 nouveaux tableaux:
		double prod1Table[2][SCREEN_HEIGHT+2][SCREEN_WIDTH+2];
		double prod2Table[2][SCREEN_HEIGHT+2][SCREEN_WIDTH+2];
Voici la fonction qu'on va utiliser pour remplir ces tables: 
		void product()
		{
			for (int t = 0; t < 2; t++)
			{
				for (int y = 0; y < SCREEN_HEIGHT+2; ++y)
					for (int x = 1; x < SCREEN_WIDTH+1; ++x)
						prod1Table[t][y][x] = table[t][y][x - 1] + 4.0f * table[t][y][x] + table[t][y][x + 1];

				for (int y = 1; y < SCREEN_HEIGHT+1; ++y)
					for (int x = 1; x < SCREEN_WIDTH+1; ++x)
						prod2Table[t][y][x] = prod1Table[t][y - 1][x] + 4.0f * prod1Table[t][y][x] + prod1Table[t][y + 1][x];
			}
		}
Maintenant le fonction Laplacian() utilise simplement le résultat de ces calculs. Ici on soustrait la valeur pour
l'élément central de la matrice.
		double Laplacian(int tab, int x, int y)
		{
			return 0.05 * prod2Table[tab][y][x] - 1.8 * table[tab][y][x];
		}
Et maintenant il suffit d'appeler product() dans la fonction update() avant tout le reste:
		void update()
		{
			sys.StartPerfCounter();

			product();

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Cette version prend 3000 ms.

Réécrire les formules principales

J'ai fait quelques tests sur les formules principales. En essayant de trouver un moyen de les calculer plus
efficacement.
Voici les résultats de mes investigations.
Nos formules actuelles ressemblent à ça: 
		nextTable[0][y][x] = A + DA * Laplacian(0, x, y) - A * B * B + f * (1.0 - A);
		nextTable[1][y][x] = B + DB * Laplacian(1, x, y) + A * B * B - (k + f) * B;
Maintenant que le laplacien est une formule simple, on peut le remplacer dans celles-ci:
		nextTable[0][y][x] = A + DA * (0.05 * prod2Table[0][y][x] - 1.8 * table[0][y][x]) - A*B*B + f*(1.0 - A);
		nextTable[1][y][x] = B + DB * (0.05 * prod2Table[1][y][x] - 1.8 * table[1][y][x]) + A*B*B - (k + f)*B;
On peut remplacer table[0][y][x] et table[1][y][x] par leurs valeurs de A et B:
		nextTable[0][y][x] = A + DA * (0.05 * prod2Table[0][y][x] - 1.8*A) - A*B*B + f*(1.0 - A);
		nextTable[1][y][x] = B + DB * (0.05 * prod2Table[1][y][x] - 1.8*B) + A*B*B - (k + f)*B;
Si on développe les parenthèses:
		nextTable[0][y][x] = A + DA * 0.05 * prod2Table[0][y][x] - 1.8*DA*A - A*B*B + f - f*A;
		nextTable[1][y][x] = B + DB * 0.05 * prod2Table[1][y][x] - 1.8*DB*B + A*B*B - k*B - f*B;
Ensuite, si on regroupe les termes qui contiennent A et B:
		nextTable[0][y][x] = 0.05 * DA * prod2Table[0][y][x] - A*B*B + A - f*A - 1.8*DA*A + f;
		nextTable[1][y][x] = 0.05 * DB * prod2Table[1][y][x] + A*B*B + B - f*B - 1.8*DB*B - k*B;
On factorise A et B:
		nextTable[0][y][x] = 0.05 * DA * prod2Table[0][y][x] - A*B*B + (1.0 - f - 1.8*DA) * A + f;
		nextTable[1][y][x] = 0.05 * DB * prod2Table[1][y][x] + A*B*B + (1.0 - f - 1.8*DB - k) * B;
Ou si on introduit 2 nouvelles variables pour les facteurs de A et B: 
		double ca = 1.0 - f - 1.8 * DA;
		double cb = 1.0 - f - 1.8 * DB - k;
		nextTable[0][y][x] = 0.05 * DA * prod2Table[0][y][x] - A*B*B + ca * A + f;
		nextTable[1][y][x] = 0.05 * DB * prod2Table[1][y][x] + A*B*B + cb * B;
Remarquez que ces 2 variables n'utilisent que des valeurs qui restent constantes quand on parcourt les pixels.
Alors on peut les sortir des boucles.
Maintenant la fonction update() ressemble à ça:
		void update()
		{
			sys.StartPerfCounter();
			double ca = 1.0 - f - 1.8 * DA;
			double cb = 1.0 - f - 1.8 * DB - k;

			product();

			for (int y = 1; y < SCREEN_HEIGHT+1; ++y)
				for (int x = 1; x < SCREEN_WIDTH+1; ++x)
				{
					double A = table[0][y][x];
					double B = table[1][y][x];

					nextTable[0][y][x] = 0.05 * DA * prod2Table[0][y][x] - A * B * B + ca * A + f;
					nextTable[1][y][x] = 0.05 * DB * prod2Table[1][y][x] + A * B * B + cb * B;
				}
				
			[...]
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Ce programme a besoin d'environ 2600 ms.

Eviter le memcpy()

Pour éviter d'avoir à copier les tables, on va se passer de nextTable et ajouter une dimension à table:
		double table[2][2][SCREEN_HEIGHT+2][SCREEN_WIDTH+2];
On va utiliser une variable curT pour dire quelle est la table courante:
		int curT = 0;
Dans la fonction product() une seule ligne va changer: 
		prod1Table[t][y][x] = table[curT][t][y][x - 1] + 4.0f * table[curT][t][y][x] + table[curT][t][y][x + 1];
Et la fonction update() va maintenant ressembler à ça:
		void update()
		{
			[...]

			for (int y = 1; y < SCREEN_HEIGHT+1; ++y)
				for (int x = 1; x < SCREEN_WIDTH+1; ++x)
				{
					double A = table[curT][0][y][x];
					double B = table[curT][1][y][x];

					table[1-curT][0][y][x] = 0.05 * DA * prod2Table[0][y][x] - A * B * B + ca * A + f;
					table[1-curT][1][y][x] = 0.05 * DB * prod2Table[1][y][x] + A * B * B + cb * B;
				}
Le memcpy() va être remplacé par une simple inversion de variable:
		curT = 1 - curT;

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Et bien sur l'initialisation et la partie dessin vont changer aussi, mais je ne vais pas vous ennuyer avec ça.
Ce code prend environ 2200 ms sur mon i7-4770K. C'est plus de 20 fois plus rapide que le premier.

Remettre le wrapping

Si on veut récupérer le bouclage autour de l'écran qu'on a retiré au début, on doit simplement "remplir" les pixels
supplémentaires tout autour de l'écran avec une copie des pixels visibles qui sont à l'opposé.
		void update()
		{
			double ca = 1.0 - f - 1.8 * DA;
			double cb = 1.0 - f - 1.8 * DB - k;

			for (int t = 0; t < 2; ++t)
			{
				for (int x = 0; x < SCREEN_WIDTH + 2; ++x)
				{
					table[curT][t][SCREEN_HEIGHT - 1][x] = table[curT][t][1][x];
					table[curT][t][0][x] = table[curT][t][SCREEN_HEIGHT - 2][x];
				}

				for (int y = 0; y < SCREEN_HEIGHT + 2; ++y)
				{
					table[curT][t][y][SCREEN_WIDTH - 1] = table[curT][t][y][1];
					table[curT][t][y][0] = table[curT][t][y][SCREEN_WIDTH - 2];
				}

				table[curT][t][0][0] = table[curT][t][SCREEN_HEIGHT - 2][SCREEN_WIDTH - 2];
				table[curT][t][0][SCREEN_WIDTH - 1] = table[curT][t][SCREEN_HEIGHT - 2][1];
				table[curT][t][SCREEN_HEIGHT - 1][0] = table[curT][t][1][SCREEN_WIDTH - 2];
				table[curT][t][SCREEN_HEIGHT - 1][SCREEN_WIDTH - 1] = table[curT][t][1][1];
			}

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Bonus: effet cool

Même quand l'écran est rempli, le motif continue d'évoluer.
Mais c'est assez difficile à voir.
Alors on va "éclairer" cette évolution en utilisant du vert.
On doit simplement mettre dans la valeur verte des pixels quelque chose qui dépend des la différence entre la table
courante et la suivante:
		int r = table[0][0][y+1][x+1] * 255;
		int b = table[0][1][y+1][x+1] * 255;
		int g = (table[0][0][y+1][x+1] - table[1][0][y+1][x+1]) * 650000;
		g = ABS(g);

		if (g > 255)
			g = 255;
		gfx.setPixel(x, y, Color(r, g, b));

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Voilà ce qu'on voit quand ça tourne:


Liens

Vidéo des deux derniers programmes de cet article

Vidéo: Réaction Diffusion en JavaScript