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
modifierUn flot dans le plan est une fonction de dans . 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 petit et on calcule et à partir de et avec les formules :
où et sont les composantes de c'est à dire que
( 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
modifierUn 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
modifierPour 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 - 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 par exemple on obtient une transition douce entre la coloration maximale et l'absence de coloration.
Dans 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
modifierJe 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
modifierLe 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
modifierPour 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
modifierL'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 dans
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
modifierLe 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
modifierLe 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 ()
{
}
}
}