Rhythmic diffusion
Réaction-diffusion
Sommaire
Introduction
Un système de réaction-diffusion est un modèle mathématique qui décrit l'évolution des concentrations d'une ou plusieurs substances.
Ces substances sont soumises à deux processus.
- La réaction: processus dans lequel la concentration des substances change.
- La diffusion: processus qui va répartir les substances dans l'espace.
Le modèle de Gray-Scott simule deux espèces de substance, A et B.
Des substances A sont constamment ajoutées dans la simulation à une certaine vitesse, tandis que les substances B en sont retirées également à vitesse définie.
Le processus de réaction fait en sorte que deux substances B convertissent une substance A, en une nouvelle substance B.
Cette simulation va être effectuée dans une grille et chaque élément de la grille sera appelé "cellule".
Toutes les cellules possèderont deux nombres (allant de zéro à un). Ces nombres représenteront la concentration des substances A et B à chaque endroit de la grille.
Équation
Pour faire évoluer la simulation, le programme va appliquer l'équation du modèle de Gray-Scott sur toutes les cellules de la grille.
Les nouvelles concentrations de substances des cellules vont dépendre de leurs voisines ainsi que de la valeur des paramètres de l'équation.
A'
=
A
+(
DA
∇²A
-AB²+
f
(1-A)
)
△t
B'
=
B
+(
DB
∇²B
B
+AB²
-(
k
+
f
)B)
△t
- Nouvelles valeurs des concentrations A et B.
- Valeurs actuelles des concentrations A et B.
- Vitesse de diffusion des substances A et B, appelées "Diffusion Rate".
- Le Laplacien en deux dimensions appliqué à A et B, qui donne la vitesse de variation du gradient de la quantité A ou B entre la cellule actuelle et les cellules environnantes.
- Vitesse d'ajout des substances A, appelée "Feed rate".
- Vitesse de suppression des substances B, appelée "Kill rate".
- "Delta t" correspond à l'écart en temps entre chaque itération, toute l'équation est mise à l'échelle par ce facteur.
A chaque actualisation, ce calcul doit être effectué pour toutes les cellules de la grille.
L'ordre dans lequel les celules vont être mises à jour, n'a pas d'importances. Ce travail est ainsi parfaitement adapté aux compute shader d'OpenGL, qui permettent d'effectuer des opérations de manière parallélisées, en utilisant la carte graphique de l'ordinateur.
Pour ce cas de figure un compute shader est nettement plus efficace que d'utiliser le processeur de l'ordinateur, qui mettrait à jour les cellules une par une.
Le moyen le plus simple d'indiquer à un compute shader quelles sont les données en entrée et où stocker les résultats et d'utiliser des textures. Ces textures seront utilisées pour représenter l'état de la simulation et feront office de grille.
Pour stocker les concentrations des substances A et B dans une texture, chaque cellule (ou autrement dit, chaque pixel de la texture), va stocker dans le canal R (Red ou Rouge) la valeur de la concentration A et dans le canal B (Blue ou Bleu) la concentration B.
Comme vu précédemment, pour calculer le prochain état de la simulation il faut son état actuel.
Il est alors nécessaire d'utiliser deux textures, une servant "d'entrée" au compute shader (celle correspondant à l'état actuelle), et l'autre servant de "sortie" qui accueillera le nouvel état de la simulation.
Après chaque calcul, le statut "d'entrée" et de "sortie" des textures seront échangés pour que la texture ayant l'état le plus avancé de la simulation, soit toujours utilisée comme entrée.
Implémentation en GLSL
Le language utilisé par un compute shader est l'OpenGL Shading Language (GLSL). Ce langage dont la syntaxe est proche de celle du langage C, permet d'exécuter du code en profitant de la puissance de la carte graphique de l'ordinateur.
Utiliser la formule de Gray-Scott ainsi que les textures avec du code peut se traduire ainsi:
1ivec2 pixelCoords = ivec2(gl_GlobalInvocationID.xy);
2vec4 cell = imageLoad(currentSimulationState, pixelCoords);
3
4float a = cell.r;
5float b = cell.b;
6
7vec3 lp = laplacian(pixelCoords);
8float newA = a + (diffusionRateA * lp.r - a*b*b + feedRate * (1.0 - a));
9float newB = b + (diffusionRateB * lp.b + a*b*b - (killRate + feedRate) * b);
10
11imageStore(nextSimulationState, pixelCoords, vec4(newA, 0, newB, 0);Les cinq premières lignes consistent à lire les concentrations de substances A et B depuis la texture correspondant à l'état actuel de la simulation.
Comme vu dans l'explication de l'équation, l'influence des cellules voisines se traduit par l'utilisation de l'opérateur Laplacien.
Le Laplacien analyse le changement d'intensité (ou dans notre cas, de concentration des subtances) de la cellule actuelle et de ses voisins. On parle alors d'un filtre de convolution 3x3.
Le poids attribué aux cellules varie en fonction de leur position (-1 pour la cellule au centre, 0.2 pour les cellules adjacentes et 0.05 pour les cellules en diagonales).
vec3 laplacian(in ivec2 uv)
{
vec3 rg = vec3(0, 0, 0);
rg += imageLoad(currentSimulationState, uv + ivec2(-1, -1)).rgb * 0.05;
rg += imageLoad(currentSimulationState, uv + ivec2(0, -1)).rgb * 0.2;
rg += imageLoad(currentSimulationState, uv + ivec2(1, -1)).rgb * 0.05;
rg += imageLoad(currentSimulationState, uv + ivec2(-1, 0)).rgb * 0.2;
rg += imageLoad(currentSimulationState, uv + ivec2(0, 0)).rgb * -1;
rg += imageLoad(currentSimulationState, uv + ivec2(1, 0)).rgb * 0.2;
rg += imageLoad(currentSimulationState, uv + ivec2(-1, 1)).rgb * 0.05;
rg += imageLoad(currentSimulationState, uv + ivec2(0, 1)).rgb * 0.2;
rg += imageLoad(currentSimulationState, uv + ivec2(1, 1)).rgb * 0.05;
return rg;
}Enfin, une fois que les nouvelles concentrations de la cellule sont calculées, ces valeurs sont écrites dans la texture 'nextSimulationState'.
Initialement la grille est remplie de substances A, la dernière étape avant de constater les résultats de la simulation est de rajouter des substances B dans la grille.
Le moyen le plus simple consiste à mettre manuellement à un, la concentration de substances B de certaines particules.
Le code ci dessous génère une zone circulaire au centre de l'écran rajoutant constamment des substances B à la simulation.
1...
2
3float a = cell.r;
4float b = cell.b;
5
6// add chemicals B at screen center
7const ivec2 screenCenter = ivec2(1920 / 2, 1080 / 2);
8const int radius = 150; // in pixel
9int x = screenCenter.x - pixelCoords.x;
10int y = screenCenter.y - pixelCoords.y;
11if (x*x + y*y < radius * radius)
12 b = 1.0;
13
14vec3 lp = laplacian(pixelCoords);
15...Comme les subtances B sont stockées dans le canal bleu et les substances A, dans le canal rouge. Le résultat n'est autre qu'un cercle bleu sur fond rouge.
Les choses deviennent plus intéressantes dès lors que l'on donne des valeurs aux paramètres de l'équation de réaction-diffusion.
Voici les portées de chaque paramètres, ainsi que quelques exemples de motifs réalisables :
- diffusion rate A: 0 - 1
- diffusion rate B: 0 - 1
- feed rate: 0 - 0.07
- kill rate A: 0 - 0.07
Feed rate = 0.014 | Kill rate = 0.04
Feed rate = 0.003 | Kill rate = 0.022
Feed rate = 0.042 | Kill rate = 0.059
Couleurs
Le résultat est fonctionnel mais ces couleurs ne sont pas agréables à regarder.
Pour avoir plus de contrôle sur le résultat final il faut rajouter une étape entre le calcul de l'équation de réaction-diffusion et l'affichage à l'écran de la texture.
Comme dit précédemment, la texture est composée de bleu et de rouge car c'est dans ces canaux de la texture que sont stockées les concentrations des substances.
Cette nouvelle étape aura pour but d'appliquer un dégradé de couleur en fonction de la concentration de chaque cellule.
Cette nouvelle étape nécessite une texture supplémentaire qui sera utilisée pour accueilir le résultat de la simulation, avec le dégradé de couleurs.
C'est cette texture qui sera désormais affichée à l'écran.
Voici à quoi ressemble maintenant le schéma de l'application :
Il est possible de tout faire dans un seul et même compute shader, cependant, pour rendre le projet plus ordonnée et modulable, un compute shader sera dédié à l'application du dégradé sur l'état le plus récent de la simulation.
Dans le cas d'un dégradé avec deux couleurs l'implémentation est assez direct.
ivec2 pixelCoords = ivec2(gl_GlobalInvocationID.xy);
vec4 existingPixel = imageLoad(inTexture, pixelCoords);
vec4 colorA = vec3(0, 0, 0, 1); // Black
vec4 colorB = vec3(1, 1, 1, 1); // White
float concentration = existingPixel.b;
existingPixel = mix(colorA, colorB, concentration);
imageStore(outTexture, pixelCoords, existingPixel);Le fait que la concentration soit entre zéro et un facilite grandement la tâche.
Il suffit alors d'utiliser la fonction mix du langage GLSL pour effectuer une interpolation linéaire et ainsi appliquer le dégradé.
En revanche pour appliquer un dégradé de plus de deux couleurs, il y a quelques étapes supplémentaires.
Premièrement il faut un moyen de représenter dans le code, toutes les couleurs du dégradé ainsi que leur position.
Dans le code ci-dessous, chaque couleur du dégradé est définie en utilisant un 'vec4'.
Les trois premières valeurs correspondent au rouge, vert et bleu de la couleur et la dernière correspond à sa position dans le dégradé (allant de zéro à un).
1ivec2 pixelCoords = ivec2(gl_GlobalInvocationID.xy);
2vec4 existingPixel = imageLoad(inTexture, pixelCoords);
3float concentration = existingPixel.b;
4
5int colorNumber = 5;
6vec4[colorNumber] gradient = {
7 vec4(0, 0, 0, 0), // Black
8 vec4(1, 0, 0, 0.25), // Red
9 vec4(0, 1, 0, 0.5), // Green
10 vec4(0, 0, 1, 0.75), // Blue
11 vec4(1, 1, 1, 1) // white
12};
13
14// Find color index
15int i = 0;
16while (concentration >= gradient[i].w && i < colorNumber - 1)
17 i++;
18
19float t = (concentration - gradient[i - 1].w) / (gradient[i].w - gradient[i - 1].w);
20vec3 rgb = mix(vec3(gradient[i - 1].xyz), vec3(gradient[i].xyz), t);
21
22existingPixel = vec4(rgb, 1.0); // Setting opacity to 1
23
24imageStore(outTexture, pixelCoords, existingPixel);D'après le dégradé décrit par le code, une forte concentration de substance B sera donc représentée par la couleur blanche, et à l'inverse, une absence de substance B sera indiquée par la couleur noire.
Là où avant une simple instruction 'mix' suffisait, il faut désormais retrouver entre quelles couleurs la concentration est (ligne 15 à 17).
Par exemple si la concentration est de 0.1, il faudra faire une interpolation linéaire avec entre la couleur noire et rouge.
Enfin, il faut retrouver le facteur de l'interpolation linéaire entre les deux couleurs autour de la concentration (ligne 19)
Voici quelques résultats utilisant le dégradé montré dans le code précédent :
Feed rate = 0.0057 | Kill rate = 0.026
Feed rate = 0.016 | Kill rate = 0.0435
Il est assez frappant de voir que les couleurs du dégradé ne sont pas présentes de manières égales. Le rouge est omniprésent, tandis que le bleu n'apparait que sur quelques pixels.
La concentration des substances B dans la simulation est généralement entre les valeurs 0 et 0.5. Avec cette information en tête, il est possible de composer des dégradés produisant des résultats hypnotisant.
Feed rate = 0.002 | Kill rate = 0.014
Feed rate = 0.054 | Kill rate = 0.061
Feed rate = 0.0018 | Kill rate = 0.0125
Feed rate = 0.024 | Kill rate = 0.049
Bruits
La simulation en l'état offre déjà un grand nombre de résultats possibles, cependant ces résultats peuvent être poussées encore plus loin.
Jusqu'à présent la valeur des paramètres de l'équation de réaction-diffusion est identique pour toutes les cellule. Il est cependant possible d'appliquer différentes valeurs de paramètres, en fonction de la position de la cellule.
La première chose nécessaire est de pouvoir stocker pour chaque cellule la valeur des quatres paramètres de l'équation de réaction-diffusion.
Pour cela, le moyen le plus simple est d'utiliser une nouvelle texture.
L'ordre des paramètres et de leur canaux est le suivant :
- Diffusion rate A : Canal rouge
- Diffusion rate B : Canal vert
- Feed rate : Canal bleu
- Kill rate : Canal opacité
La question qui se pose ensuite est comment décider des valeurs qui seront assignées dans cette nouvelle texture ?
Les réponses ici sont infinies, elles peuvent aller d'écrire des valeurs aléatoirement à utiliser une image existante ou encore créer diverses formes.
Ici, c'est l'usage de textures procédurales tel que le bruit de Perlin ou de Voronoï qui a été choisi.
Ces deux types de bruits offrent une transition douce entre les zones sombres et claires, rendant le résultat naturel voir presque organique.
A partir de plusieurs paramètres pour régler l'échelle, l'intensité ou encore un décalage de ces bruits, c'est un nouveau compute shader qui va écrire sur une nouvelle texture les valeurs des paramètres à appliquer à la simulation.
Ainsi, les valeurs des paramètres de la réaction-diffusion ne seront plus issues de valeurs constantes dans le code, mais proviendront de la nouvelle texture qui sera fournie au compute shader.
Voici désormais le schéma de l'état final de la simulation :
Voici plusieurs résultats utilisants les bruits sur différentes propriétés de la simulation :
Conclusion
Cette simulation et tous ses paramètres offrent une large variété de motifs accessibles.
Il y a bien sur moyen d'encore améliorer ceci en rajoutant de nouveaux types de bruits, ou en appliquant différentes couleurs en fonction de l'intensité de certains paramètres.
Les performances sont suffisantes, mais pourrait néanmoins être améliorées.
Avec du recul, l'usage d'un compute shader pour appliquer un dégradé est peut-être excessif et aurait pu être fait dans un fragment shader d'OpenGL.
La structure de la simulation est quand à elle satisfaisante et pourrait sans difficultée accueilir de nouvelles fonctionnalitées.
Celle-ci a notamment été pensée pour que par la suite, l'analyse audio puisse facilement interagir avec elle.
Il sera alors aisé de modifier l'état de la simulation et de ressentir au maximum l'influence de la musique sur les effets visuels.
Sources
Karl Sims | Reaction-Diffusion Tutorial : http://karlsims.com/rd.html
MIT CSAIL | Gray Scott Model of Reaction Diffusion : https://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/
Ciphrd | Reaction diffusion on shader : https://ciphrd.com/2019/08/24/reaction-diffusion-on-shader/
The Coding Train| Reaction Diffusion Algorithm in p5.js: https://www.youtube.com/watch?v=BV9ny785UNc