Physique et mathématiques visuelles/Comment filmer des flots 2D

Ce chapitre donne et explique un programme, écrit en C sharp, avec lequel il est possible de filmer des flots 2D.

Le calcul des trajectoires

Un flot $\mathbf {f}$ dans le plan est une fonction de $\mathbb {R} ^{2}$ dans $\mathbb {R} ^{2}$. En chaque point du plan le vecteur vitesse du flot est défini.

Quand le flot est connu, le calcul d'une trajectoire à partir d'un point initial peut être très simple.

On choisit un intervalle de temps $dt$ petit et on calcule $x(t+dt)$ et $y(t+dt)$ à partir de $x(t)$ et $y(t)$ avec les formules :

$x(t+dt)=x(t)+f_{1}(x(t),y(t))dt$

$y(t+dt)=y(t)+f_{2}(x(t),y(t))dt$

où $f_{1}$ et $f_{2}$ sont les composantes de $\mathbf {f}$ c'est à dire que

$\mathbf {f} (x,y)=(f_{1}(x,y),f_{2}(x,y))$

( Mais pour améliorer la précision du calcul, et donc sa puissance, on a intérêt à se servir d'un meilleur algorithme. Celui de Runge-Kutta (présenté dans Blanchard, Devaney et Hall 2011) est particulièrement performant et a été adopté par ce programme dans la fonction R2 Go( int nstep, double t) dans le fichier R2.cs. )

Le nuage de particules entraînées par le flot

Un flot transparent et vu de l'intérieur est aussi invisible que l'air. Il faut des nuages ou des poussières pour voir les écoulements.

Si on veut voir une seule ligne de courant, il suffit de créer des particules aléatoirement dans le temps mais toujours à la même place - on a supposé que le flot est stationnaire, c'est à dire qu'il ne varie pas dans le temps. Tout se passe comme si on déposait en permanence un peu d'encre en un point du flot qui l'entraîne.

Dans ce programme, les lignes de courant bleues sont issues d'une même ligne source (ou de plusieurs) qui peut être, mais pas nécessairement, le grand cercle qui entoure l'écran.

Pour visualiser le flot entre les lignes de courant, on crée également des particules noires réparties aléatoirement sur la ligne source.

Pour avoir un film périodique les particules sont recréées à la même place à intervalles de temps réguliers.

Le dessin en haute définition

Pour dessiner des mouvements bien réguliers, il ne faut pas représenter les particules par des pixels ou des groupes de pixels, parce qu'on les verrait sauter de pixel en pixel.

Une particule est représentée par sa position - un vecteur dans $\mathbb {R} ^{2}$ - et on calcule la façon dont elle colore un pixel à partir de sa distance au centre du pixel. Il faut choisir une fonction qui est maximale pour une distance nulle et nulle pour une distance supérieure au rayon de la particule. En choisissant une forme en $\cos ^{4}$ par exemple on obtient une transition douce entre la coloration maximale et l'absence de coloration.

Dans $\mathbb {R} ^{2}$ un pixel est représenté non par un point mais par une surface carrée. Le pixel 0,0 par exemple, c'est à dire le pixel en haut à gauche, est le carré entre les points (0,0) et (1,1), il est donc centré en (0.5,0.5).

Configuration requise

Je compile ce programme avec MonoDevelop dans Ubuntu (Linux) et et je me sers d'Image Magick pour l'assemblage d'une série d'images en un seul film.

Avec quelques adaptations ce programme devrait tourner dans n'importe quel environnement. Il suffit de créer des images Bitmap, de pouvoir les enregistrer en .gif, d'avoir les fonctions SetPixel, qui colore les pixels, et GetPixel, qui donne leur couleur présente, et quelques fonctions mathématiques habituelles : Sin, Cos, Sqrt...

Comment assembler toutes les images en un seul film .gif

Le programme crée une série de 50 fichiers image, numérotés de Im100.gif à Im149.gif et les place dans un dossier Film qu'il faut avoir créé au préalable.

Image Magick permet d'assembler ces 50 images en un seul fichier de film .gif.

Si vous avez Ubuntu, Image Magick est déjà installé sur votre ordinateur.

Pour l'utiliser, il faut ouvrir un terminal (si vous ne savez pas faire, essayez Terminal avec l'onglet Chercher sur votre ordinateur en haut à gauche) puis se déplacer dans l'arborescence des dossiers pour trouver le dossier Film (essayez cette adresse si vous ne savez pas faire). Il faut alors écrire les trois lignes de code suivantes, chacune suivie d'une frappe de la touche Enter :

convert -delay 4 \
Im*.gif \
-loop 0 ./Nom_du_film.gif

4 veut dire qu'il y a 4/100 de seconde entre chaque image, donc 25 images par seconde.

-loop 0 veut dire que le film tourne en boucle infinie. Pour n passages, il suffit de choisir -loop n.

Tous les fichiers du programme

Pour compiler ce programme, vous devez rassembler tous les fichiers Program.cs, Film.cs, Calculate.cs, R2.cs, Draw.cs, Cloud.cs, Line.cs et LineSource.cs dans un seul projet (un console project avec Monodevelop).

.cs veut dire C sharp ou C#.

Le fichier de définition d'un film

L'en-tête donne les noms des paramètres qui servent à calculer un film.

La fonction R2 Flow(R2 v) est un flot 2D, c'est à dire une fonction de $\mathbb {R} ^{2}$ dans $\mathbb {R} ^{2}$

Ici elle est définie par Flow(x,y)=(x+y,y)

R2 est la classe des couples de réels.

Un nombre réel est un double pour la machine. Cela veut dire double précision, par opposition à la simple précision (le type float, nommé ainsi pour sa virgule flottante) qui n'est plus très utilisée.

int veut dire nombre entier (integer).

Le fichier Film.cs :

using System; // Technical
using System.Drawing; // Technical : external library for the class Color

namespace See // Technical
{
	public class Film // There will be only one element in this class, the film to be made.
	{
		// List of names of parameters (Technical)
			// Parameters of the film
		public Color screenColor; 
		public int wsc, hsc;// w : width / h : height / sc : screen
		public int nim;     // n : number of / im : images

		public double x1, y1, x2, y2;// rectangle of observation of the flow
		public double T;             // duration of observation

		public int nper;// number of periods of creation of particles
		public double t;// duration of a period

		public int nl;// number of line sources  

		public LineSource[] ls;// list of line sources

			// Technical parameters
		public int nsteptot; // total number of steps of calculus
		public int nstep;    // number of steps in the calculus of one period of creation of particles
		public double inw;   // internal width of a particle (real number of pixels)
		public double exw;   // external width of a particle
		public double antial;// real number of pixels of anti-aliasing

		//
		// The flow to be seen : Flow(x,y)=(x+y,y)
		//
		public static R2 Flow(R2 v) // Flow is a function from R2 to R2
		{
			R2 w= new R2 (); // w is in R2 
			w.x = v.x + v.y ;
			w.y = v.y ;
			return w; // w is the value of the function at the point v : w = Flow(v)
		}

		// Initialization of the parameters of the film
		public Film ()
		{
			Color black = Color.FromArgb (0, 0, 0);
			Color blue = Color.FromArgb (0, 80, 255);
			// Color.FromArgb(r,g,b) means an intensity r of red, g of green and b of blue, 
			// against a black background. 
			screenColor=Color.FromArgb (255, 255, 255);// White screen

			wsc = 1920;// width of screen (in pixels)
			hsc = 1080;// height of screen
			nim =   50;// number of images

			x1 = -16.0;// rectangle of observation of the flow
			y1 = +9.0;
			x2 = +16.0;
			y2 = -9.0;
			T = -10.0; // duration of observation; Negative means that time flows backward.

			nper = 25; // number of periods of creation of particles

			nl = 1;    // number of line sources  

			ls  = new LineSource [nl]; // (technical)

			ls [0]   = new LineSource (black,blue); // black is the color of randomly created paticles
													// blue is the color of streamlines
			ls [0].l = new Line (2, 0.0, 0.0, Calculate.ObservationRadius(x1,y1,x2,y2));
			// 2 means circular line. The other three parameters are
			// the coordinates x0,y0 of its center and its radius r
			ls [0].a = 0.0; 	 // angle from the horizontal of the first point of this circular arc
			ls [0].b = 2*Math.PI;// angle of the last point
			ls [0].npartrnd = 200000;// number of particles randomly created on the line
			ls [0].nstr     = 40;    // number of streamlines originated on the line
			ls [0].npartstr = 3000;  // number of particles on each streamline

			// For a straight line, try this code :
			/*
			ls [0] = new LineSource (black,blue);
			ls [0].l = new Line (1, -8.0, -8.0, 8.0, 8.0);
			// 1 means straight line. The other four parameters are 
			// the coordinates x1,y1,x2,y2 of its extremities.  
			ls [0].npartrnd = 200000;// number of particles randomly created on the line
			ls [0].nstr = 20;        // number of streamlines originated on the line
			ls [0].npartstr = 3000;  // number of particles on each streamline
			*/

			// Technical parameters
			nsteptot=100;// total number of steps of calculus
			inw = 0.0;   // internal width of a particle (real number of pixels)
			exw = 3.0;   // external width of a particle
			antial = 3.0;// real number of pixels of antialiasing

			// Calculated parameters
			t = T/nper; // duration of a period of creation of particles
			nstep = nsteptot / nper; // number of steps in the calculus of one period of creation of particles
		}
	}
}

Le programme Main

Le fichier Program.cs :

using System;

namespace See
{
	class MainClass
	{
		public static void Main (string[] args)
		{
			Film f= new Film(); // f is the film

			Cloud cloud = Calculate.InitialConditions (f);
			// cloud is the cloud of particles driven by the current

			Draw.FilmCloud(f, cloud);

			Console.WriteLine ("The end");
		}
	}
}

Les autres fichiers

Le fichier Calculate.cs :

using System;
using System.Drawing;

namespace See
{
	public class Calculate
	{
		public static double ObservationRadius(double x1, double y1, double x2, double y2)
		{
			double dx = 0.5 * (x2 - x1);
			double dy = 0.5 * (y2 - y1);
			return Math.Sqrt (dx * dx + dy * dy);
		}

		public static Cloud InitialConditions(Film f)
		{
			// Calculation of the initial clouds
			int nn=2*f.nl;// number of sources
			Cloud[] cloud = new Cloud[nn];


			int npt = 0;//nombre de points au total, à incrémenter
			int i;

			// Line sources
						// Random
			for (i = 0; i < f.nl; i++) {
				Console.Write ("    Line source. Random ");
				Console.Write (i);
				Console.Write ("  ");
				cloud [i] = new Cloud (f.ls[i].npartrnd);
				cloud [i].InitLineSource (f, f.ls[i].l, f.ls[i].a, f.ls[i].b);
				cloud [i].InitBichromatic(f.ls[i].col1, f.ls[i].col2, f.ls[i].mix);
				npt += cloud [i].n;
			}
						// Streamlines (dotted line source)
			for (i = f.nl; i < 2*f.nl; i++) {
				Console.Write ("    Line source. Streamlines ");
				Console.Write (i);
				Console.Write ("  ");
				cloud [i] = new Cloud (f.ls[i-f.nl].nstr*f.ls[i-f.nl].npartstr);
				cloud [i].InitDottedLineSource (f, f.ls[i-f.nl].nstr, f.ls[i-f.nl].l, f.ls [i - f.nl].a,
					f.ls [i - f.nl].b, f.ls[i-f.nl].wst);
				cloud [i].InitBichromatic(f.ls[i-f.nl].colstr1, f.ls[i-f.nl].colstr2, f.ls[i-f.nl].mix);
				npt += cloud [i].n;
			}
			return new Cloud (cloud, nn);
		}
	}
}

Le fichier R2.cs :

using System;

namespace See
{
	public class R2
	{
		public double x,y;

		public double Norm()
		{
			return Math.Sqrt(x*x+y*y);
		}
		public R2 Plus(R2 v)
		{
			R2 sum = new R2 ();
			sum.x = x + v.x;
			sum.y = y + v.y;
			return sum;
		}
		public R2 Times(double lambda)
		{
			R2 mult= new R2();
			mult.x= lambda*x;
			mult.y= lambda*y;
			return mult;
		}
		public R2 Go(int nstep, double t)
		{
			int i; 
			double dt;
			dt = t / nstep;
			R2 via1 = new R2 ();
			R2 via2 = new R2 ();
			R2 via3 = new R2 ();
			R2 via4 = new R2 ();
			R2 dvia1 = new R2 ();
			R2 dvia2 = new R2 ();
			R2 dvia3 = new R2 ();
			R2 dvia4 = new R2 ();

			via1.x = x;
			via1.y = y;

			for (i = 0; i < nstep; i++) {
				// Runge-Kutta method
				dvia1 = Film.Flow (via1);
				via2 = via1.Plus (dvia1.Times (0.5 * dt));

				dvia2 = Film.Flow (via2);
				via3 = via1.Plus (dvia2.Times (0.5 * dt));

				dvia3 = Film.Flow (via3);
				via4 = via1.Plus (dvia3.Times (dt));

				dvia4 = Film.Flow (via3);

				via1 = via1.Plus (dvia1.Plus (dvia2.Times (2.0).Plus (dvia3.Times (2.0).Plus (dvia4))).Times (1.0 / 6).Times (dt));
			}

			return via1;
		}
	}
}

Le fichier Draw.cs :

using System;
using System.Drawing;
using System.Drawing.Imaging;

namespace See
{
	class Draw
	{
		public static void Screen(Bitmap image, Color color)
		{
			int i, j;
			for (i = 0; i < image.Width; i++) {
				for (j = 0; j < image.Height; j++) {
					image.SetPixel (i, j, color);
				}
			}
		}
		public static void Point(Film f, Bitmap im, double x1, double y1, Color colorp)
		{
			double xx1, yy1;
			xx1 = f.wsc * (x1 - f.x1) / (f.x2 - f.x1);
			yy1 = f.hsc * (y1 - f.y1) / (f.y2 - f.y1);
			int n1 = (int) Math.Floor(xx1-f.exw-0.5);
			int n2 = (int) Math.Floor(xx1+f.exw+1.5);
			int m1 = (int) Math.Floor(yy1-f.exw-0.5);
			int m2 = (int) Math.Floor(yy1+f.exw+1.5);
			if (n1<0){n1 = 0;};
			if (m1<0){m1 = 0;};
			if (n2>im.Width){n2 = im.Width;};
			if (m2>im.Height){m2 = im.Height;};
			double posx, posy, varcolor;
			double dist; //distance du centre du pixel à la ligne;
			Color coulp, colorf;
			int i, j;
			for (i = n1; i < n2; i++) {
				posx = i + 0.5;
				for (j = m1; j < m2; j++) {
					posy = j + 0.5;
					dist = Math.Sqrt ((xx1 - posx) * (xx1 - posx) + (yy1 - posy) * (yy1 - posy));
					if (dist < f.exw) {
						varcolor = RoundedStep (dist, f.inw, f.exw);
						colorf = im.GetPixel (i, j);
						coulp= Color.FromArgb (colorf.R + (int)(varcolor * (colorp.R - colorf.R)), 
							colorf.G + (int)(varcolor * (colorp.G - colorf.G)), 
							colorf.B + (int)(varcolor * (colorp.B - colorf.B)));
						if (colorf == f.screenColor) {
							im.SetPixel (i, j, coulp);
						} else {
							im.SetPixel (i,j, Color.FromArgb ( (int)(0.5*(coulp.R + colorf.R)), 
								(int)(0.5*(coulp.G + colorf.G)), 
								(int)(0.5*(coulp.B + colorf.B))));
						}
					}
				}
			}
		}
		public static void Cloud(Film f, Bitmap im, Cloud cloud)
		{
			int i;
			for (i = 0; i < cloud.n; i++) 
			{
				Draw.Point(f, im, cloud.pt [i].x, cloud.pt [i].y, cloud.col[i]);
			}

		}
		public static void FilmCloud(Film f, Cloud cloud)
		{
			Bitmap im = new Bitmap (f.wsc, f.hsc);
			int i;
			for (i = 0; i < f.nim; i++) {
				Draw.Screen (im, f.screenColor);
				Draw.Cloud (f, im, cloud);
				im.Save ("Film/Im"+ Convert.ToString(i+100)+".gif", ImageFormat.Gif);
				cloud.Anim (f.nstep, f.t/f.nim);
			}
		}
		// RoundedStep is for anti-aliasing
		public static double RoundedStep(double x, double win, double wout)
		{
			if (x > wout) {
				return 0.0;
			}
			if (x < win) {
				return 1.0;
			}
			return Math.Pow(Math.Cos(0.5*Math.PI* (x - win) / (wout - win)),4);
		}
	}
}

Le fichier Cloud.cs :

using System;
using System.Drawing;

namespace See
{
	public class Cloud // A cloud of particles
	{
		public int n;      // number of particles
		public R2[] pt;    // positions of particles
		public Color[] col;// colors of particles

		public Cloud (int npt)
		{
			n = npt;
			pt = new R2[n];
			col = new Color[n];
		}
		public Cloud ( Cloud[] cloud, int nn) // Union of nn clouds in a single one
		{
			int i,j,k;
			n = 0;
			for (i = 0; i < nn; i++) {
				n += cloud [i].n;
			}
			pt = new R2[n];
			col = new Color[n];
			k = 0;
			for (i = 0; i < nn; i++) {
				for (j = 0; j < cloud [i].n; j++) {
					pt[k] = new R2 ();
					pt [k].x = cloud [i].pt [j].x;
					pt [k].y = cloud [i].pt [j].y;
					col [k] = cloud [i].col [j];
					k++;
				}
			}
		}
		public void InitBichromatic(Color color1, Color color2, double mix) 
		// Bichromatic clouds are often nice but difficult for the GIF format which is limited to 256 colors 
		// (there are many shades between the two colors and with the color of the screen)
		{
			int i;
			Random hasard;
			hasard = new Random ();
			for (i = 0; i < n; i++) {
				if (hasard.NextDouble () < mix) {
					col [i] = color1;
				} else {
					col [i] = color2;
				}
			}
		}	
		public void InitLineSource(Film f, Line l , double a, double b)
		// t is the duration of a period of creation of particles.
		// Particles are created at each period at the same places. This is the trick to eternal return.
		{
			int i,j,n2,m;
			n2 = (int)(n / f.nper);
			m = 0;
			Random hasard = new Random ();
			for (i = 0; i < n2; i++) {
				pt[i*f.nper] = Line.Def (l, a + (b-a)*hasard.NextDouble ());
				pt[i*f.nper] = pt[i*f.nper].Go (f.nstep, f.t * hasard.NextDouble ());
				m++;
				if (f.nper>1){
					pt[i*f.nper+1] = pt[i*f.nper].Go (f.nstep, -f.t );
					m++;
					if (f.nper > 2) {
						pt [i * f.nper + 2] = pt[i*f.nper].Go (f.nstep, f.t);
						m++;
						for (j = 3; j < f.nper; j++) {
							pt [i * f.nper + j] = pt [i * f.nper + j - 1].Go (f.nstep, f.t);
							m++;
						}
					}
				}
				if (m > 999) {
					m = 0;
					Console.Write (i*f.nper);
					Console.Write (" ");
				}
			}
		}
		public void InitDottedLineSource(Film f, int n1, Line l, double a, double b, double dx)
		{
			int i,j,k,m,n2;
			m = 0;
			n2 = (int)(n / (n1*f.nper));
			double stepx;
			stepx = (b - a) / n1;
			Random hasard = new Random ();
			for (i = 0; i < n1; i++) {
				for (j = 0; j < n2; j++) {
					pt[(i*n2+j)*f.nper] = Line.Def(l, a + stepx* i + hasard.NextDouble () * dx);
					pt[(i*n2+j)*f.nper] = pt[(i*n2+j)*f.nper].Go (f.nstep, f.t * hasard.NextDouble ());
					m++;
					if (f.nper > 1) {
						pt [(i*n2+j)*f.nper + 1] = pt[(i*n2+j)*f.nper].Go (f.nstep, -f.t);
						m++;
						if (f.nper > 2) {
							pt [(i*n2+j)*f.nper + 2] = pt[(i*n2+j)*f.nper].Go (f.nstep, f.t);
							m++;
							for (k = 3; k < f.nper; k++) {
								pt [(i*n2+j)*f.nper + k] = pt [(i*n2+j)*f.nper + k - 1].Go (f.nstep, f.t);
								m++;
							}
						}
					}
					if (m > 999) {
						m = 0;
						Console.Write ((i*n2+j)*f.nper);
						Console.Write (" ");
					}
				}
			}
		}
		public void Anim(int nstep, double t)
		{
			int i;
			R2 end =new R2();
			for (i = 0; i < n; i++) {
				end = pt [i].Go (nstep, t);
				pt[i].x=end.x;
				pt[i].y=end.y;
			}
		}
	}
}

Le fichier Line.cs :

using System;

namespace See
{
	public class Line // Parameterized line
			  // Def(l,x) gives the position of the point with parameter x on the line l
	{
		public int num; // number for the type of line
				// 1 : straight line
				// 2 : circle
		public double x1,x2,x3,x4; // parameters for the definition of the line
					   // straight line from x1,x2 to x3,x4
					   // circle centered on x1,x2 and radius x3

		public Line(int pnum, double px1, double px2, double px3)
		{
			num = pnum;
			x1 = px1;
			x2 = px2;
			x3 = px3;
			x4 = 0.0; // unused parameter for a circle
		}

		public Line(int pnum, double px1, double px2, double px3, double px4)
		{
			num = pnum;
			x1 = px1;
			x2 = px2;
			x3 = px3;
			x4 = px4;
		}

		public static R2 Def(Line l, double x)
		{
			R2 f= new R2 ();
			f.x = 0.0;
			f.y = 0.0;
			if (l.num==1){// straight line from x1,x2 to x3,x4
				f.x = l.x1+x*(l.x3-l.x1);
				f.y = l.x2+x*(l.x4-l.x2);
			}
			if (l.num==2){// circle / center x1,x2 / radius x3 
				f.x = l.x1+l.x3*Math.Cos(x);
				f.y = l.x2+l.x3*Math.Sin(x);
			}
			// other types of lines can be added here
			return f;
		}
	}
}

Le fichier LineSource.cs :

using System;
using System.Drawing;

namespace See
{
	public class LineSource // A line which acts as a source of particles
	{
		public Line l;// parameterized line identifier

		public double a,b; // parameters of the beginning and of the end of the line source on the parameterized line

		public int npartrnd;//number of particles randomly distributed on the line
		public Color col1, col2;// colors of particles
		public double mix;      // coefficient of random color mixing 
		// 0.0 : col1 alone / 1.0 : col2 alone

		public int nstr;// number of streamlines  (particles periodically distributed on the line)
		public int npartstr; // number of particles on a streamline
		public double wst; // width of streamline
		public Color colstr1, colstr2;// colors of particles 
		public double mixstr;         // coefficient of random color mixing 
		// 0.0 : colstr1 alone / 1.0 : colstr2 alone

		public LineSource (Color col, Color colstr)
		{
			// Initialization
			col1 = col;
			col2 = col;
			colstr1 = colstr;
			colstr2 = colstr;
			mix = 0.0;   // monocromatic cloud
			mixstr = 0.0;// monocromatic cloud
			a = 0.0;// parameter of the beginning for a straight line
			b = 1.0;// parameter of the end for a straight line
			wst = 0.0;// infinitely thin streamline
		}

		public LineSource ()
		{
		}
	}
}