IFT6145 - Hiver 2007
TP4 - Flux Optique

Description
(Faites attention, le TP a change par rapport a l'an passe)

1. Implantez une fonction de calcul du flux optique basée sur l'algorithme "CONSTANT_FLOW" (CF) du livre Trucco & Verri. Pour le lissage, je suggere 4 voisins, mais vous pouvez utiliser autre chose (8 voisins par exemple). Pour la constante, choisissez une valeur qui donne de bons resultats. Utilisez cette valeur pour toutes les sequences.
2. Implantez une fonction de calcul du flux optique basee sur l'algorithme de "Horn & Schunk" (HS) (voir dans les notes). Choisissez une taille de voisinage qui donne de bons resulats, mais utilisez la meme pour toutes les sequences.
3. Vous devez utiliser (et comparer) quatre variantes du calcul des gradients spatio-temporels:
    (G) la derivée de la gaussienne (voir dkernel[]) avec un sigma bas (<=1),
    (D) la différence de deux gaussiennes (voir dogkernel[]) avec le meme sigma que (A),
4. Exécuter vos deux algorithmnes (CF et HS) sur les 7 séquences fournies. Les séquences contiennent chacune 9 images. Vous devez calculer le flux optique pour l'image du milieu de la sequence.
3. Montrez vos (2 x 2 x 7) résultats sous la forme de champs de vecteurs. Chaque resultat doit etre clairement identifie par l'etiquette algorithme-filtrage-sequence. Par exemple:  CF-G-1, HS-D-6, etc...
4.  Donnez les choix de constantes de lissage et autres parametres que vous avez choisis. Commentez brievement vos résultats (qualité générale, discontinuités, impact de sigma, differences entre les deux algorithmes,choix des parametres, ...)
Notes:
  - Vous pouvez utiliser le code fourni pour calculer les dérivées spatio-temporelles.
  - Vous pouvez utiliser la fonction ListConvolve de Mathematica
  - Remettre un notebook Mathematica

Initialisation

<< Graphics`Graphics`

<< Graphics`PlotField`

Pour ne pas avoir certains message d'erreur ennuyants....

Off[General::spell1];
Off[General::spell];

Pour ne pas afficher un graphique,utiliser dfi,sinon dfd
Example:Plot[Sin[x],{x,0,Pi},dfi]

dfi = DisplayFunction -> Identity;
dfd = DisplayFunction -> $DisplayFunction;

Options pour l'affichage

SetOptions[ListDensityPlot, PlotRange -> All, 
        Mesh -> False, AspectRatio -> Automatic];

SetOptions[DensityPlot, PlotRange -> All, Mesh -> False, 
        AspectRatio -> Automatic, PlotPoints -> 40];

SetOptions[Plot, PlotRange -> All];

Fonctions de base

Load "pgm"  or "pnm" images (ascii or raw)

\(\(loadPnm[name_] := 
        Block[{f, type, xs, ys, 
            data}, \((f = OpenRead[name]; 
If[
              f == $Failed, Return[]]; 
type = 
              Read[f, String]; 
{xs, ys} = 
              ToExpression[
                StringReplace["{" <> 
                    getString[
                      f] <> "}", " " -> ","]]; 
Read[f, String]; \  (*nb\ grey\ level*) data = 
              Switch[type, 
"P2", 
Partition[
                  ReadList[f, Number, xs*ys, RecordLists -> False], 
                  xs], 
"P3", 
Partition[
                  Partition[
                    ReadList[f, Number, xs*ys*3, RecordLists -> False], 
                    3], xs], "P5", 
Partition[
                  ReadList[f, Byte, xs*ys, RecordLists -> False], 
                  xs], 
"P6", 
Partition[
                  Partition[
                    ReadList[f, Byte, xs*ys*3, RecordLists -> False], 
                    3], xs]]; 
Close[
              f]; 
Reverse[data])\)];\)\)

Lit une string,en sautant les commentaires

\(\(getString[f_] := 
        Module[{}, \((
debut = 
              Read[f, Character]; 
While[\((debut == "#")\), \((Read[f, String]; \ 
                debut = Read[f, Character])\)]; 
\((debut <> 
                Read[f, String])\)
)\)];\)\)

Sauvegarde une image en pgm.
    (Suppose tableau d'entiers entre 0 et 255)

\(\(savePgm[name_, img_] := 
        Block[{f, type, xs, ys, 
            data}, \((f = OpenWrite[name]; 
If[
              f == $Failed, Return[]]; 
WriteString[
              f, "P5
"]; 
WriteString[f, 
              ToString[Length[img[\([1]\)]]]\  <> " " <> 
                ToString[
                  Length[img]]\  <> \ "
255
"]; 
WriteString[f, 
              FromCharacterCode[
                Flatten[Reverse[img]]]]; 
Close[
              f];
)\)];\)\)

scaleImg[img_] := 
        Block[{min, 
            max}, (
{min, max} = {Min[img], 
                Max[img]}; 
Round[(img - min)/(max - 
                        min)*255 + 
                0.4999]
)];

Support

norm[v_] := Sqrt[v . v];

Sous-échantillonne un tableau 2D d'un facteur subr et subc

subSample[v_, subr_, subc_] := 
Table[
          v[[r, c]], {r, 1, Length[v], subr}, {c, 1, Length[v[[1]]], 
            subc}];

Fonctions de calcul du flux optique

Rien pour l'instant...
A vous de definir les fonctions...

Sequence #1 : Carré + flux constant

\(\(nom = \(\(("sequences/carre/" <> #)\) &\)\  /@ \ {"img3.pgm", "img4.pgm", "img5.pgm", "img6.pgm", "img7.pgm", "img8.pgm", "img9.pgm", "img10.pgm", "img11.pgm"};\)\)

seq1 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq1]

ListDensityPlot  /@  seq1;

Sequence #2 : Cercle + flux constant

\(\(nom = \(\(("sequences/cercle/" <> #)\) &\)\  /@ \ {"img3.pgm", "img4.pgm", "img5.pgm", "img6.pgm", "img7.pgm", "img8.pgm", "img9.pgm", "img10.pgm", "img11.pgm"};\)\)

seq2 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq2]

ListDensityPlot  /@  seq2;

Sequence #3 : Blob + flux radial

\(\(nom = \(\(("sequences/blob1/" <> #)\) &\)\  /@ \ {"img1.pgm", "img2.pgm", "img3.pgm", "img4.pgm", "img5.pgm", "img6.pgm", "img7.pgm", "img8.pgm", "img9.pgm"};\)\)

seq3 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq3]

ListDensityPlot  /@  seq3;

Sequence #4 : Blob + flux discontinu

\(\(nom = \(\(("sequences/blob2/" <> #)\) &\)\  /@ \ {"img1.pgm", "img2.pgm", "img3.pgm", "img4.pgm", "img5.pgm", "img6.pgm", "img7.pgm", "img8.pgm", "img9.pgm"};\)\)

seq4 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq4]

ListDensityPlot  /@  seq4;

Sequence #5 : Blob + flux rotationel

\(\(nom = \(\(("sequences/blob3/" <> #)\) &\)\  /@ \ {"img1.pgm", "img2.pgm", "img3.pgm", "img4.pgm", "img5.pgm", "img6.pgm", "img7.pgm", "img8.pgm", "img9.pgm"};\)\)

seq5 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq5]

ListDensityPlot  /@  seq5;

Sequence #6 : Rubic

\(\(nom = \(\(("sequences/rubic/" <> #)\) &\)\  /@ \ {"rubic04.pgm", "rubic05.pgm", "rubic06.pgm", "rubic07.pgm", "rubic08.pgm", "rubic09.pgm", "rubic10.pgm", "rubic11.pgm", "rubic12.pgm"};\)\)

seq6 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq6]

ListDensityPlot  /@  seq6;

Sequence #7 : Taxi

\(\(nom = \(\(("sequences/taxi/" <> #)\) &\)\  /@ \ {"taxi04.pgm", "taxi05.pgm", "taxi06.pgm", "taxi07.pgm", "taxi08.pgm", "taxi09.pgm", "taxi10.pgm", "taxi11.pgm", "taxi12.pgm"};\)\)

seq7 = N[loadPnm[#]] &  /@  nom;

Dimensions[seq7]

ListDensityPlot  /@  seq7;

Masques de convolution

Gaussienne de moyenne 0 et variance v

gauss[x_, v_] := 
        1/Sqrt[2 Pi v] Exp[-x^2/(2 v)];

Plot[{gauss[x, 10], gauss[x, 1]}, {x, -10, 10}];

Dérivée de la gaussienne

dgauss[x_, v_] = D[gauss[x, v], x];

Plot[{dgauss[x, 10], dgauss[x, 1]}, {x, -10, 10}];

Difference de gaussiennes

dog[x_, v_] := (gauss[x + 1, v] - gauss[x - 1, v])/
          2;

Plot[{dog[x, 10], dgauss[x, 1]}, {x, -10, 10}];

Filtres

Filtre de lissage spatio-temporel (1 dimension parce que c'est séparable)
écart type  "sig", la taille total du filtre est 2t+1 (donc 't' éléments de chaque coté de 0)

kernel[sig_, t_] := 
        Table[N[gauss[x, sig^2]], {x, -t, t}];

ListPlot[kernel[1.5, 4], PlotJoined -> True];

La somme des éléments devrait être près de 1.

Plus  @@  kernel[1.5, 4]

ListPlot[kernel[1, 2], PlotJoined -> True];

Plus  @@  kernel[1, 2]

Filtre de dérivée spatio-temporelle (1 dimension parce que séparable)
écart type "sig", la taille total du filtre est 2t+1 (donc 't' éléments de chaque cotés de 0)

dkernel[sig_, t_] := 
        Table[N[dgauss[x, sig^2]], {x, -t, t}];

ListPlot[dkernel[1.5, 4], PlotJoined -> True, 
        PlotRange -> All];

ListPlot[dkernel[1, 4], PlotJoined -> True, 
        PlotRange -> All];

ListPlot[dkernel[0.5, 3], PlotJoined -> True, 
        PlotRange -> All];

Filtre de difference de gaussiennes (1 dimension parce que séparable)
écart type "sig", la taille total du filtre est 2t+1 (donc 't' éléments de chaque cotés de 0)

dogkernel[sig_, t_] := 
        Table[N[dog[x, sig^2]], {x, -t, t}];

ListPlot[dogkernel[1.5, 4], PlotJoined -> True, 
        PlotRange -> All];

ListPlot[dogkernel[1, 4], PlotJoined -> True, 
        PlotRange -> All];

ListPlot[dogkernel[0.5, 3], PlotJoined -> True, 
        PlotRange -> All];

Filtrage Spatio-temporel

Convolution d'une séquence avec un filtre de taille 2t+1 en X, ou Y ou T seulement.
La convolution est effectuée seulement pour les pixels à plus de 't' pixels du bord
La séquence est de trois dimensions : t,y,x dans cet ordre.

Convolution de la séquence spatio-temporelle 's' avec le kernel X, le kernel Y, et le kernel T
(suppose que la convolution est séparable)
(suppose que s contient les dimension [T,Y,X] dans l'ordre)

\(\(convolution3D[s_, kernelX_, kernelY_, kernelT_] := 
        Block[{mx, my, mt, i1, i2, 
            i3}, \((
mx = \((Length[kernelX] + 1)\)/
                2; 
my = \((Length[kernelY] + 1)\)/
                2; 
mt = \((Length[kernelT] + 1)\)/
                2; 
 (*\ Convolution\ en\ X, 
              pour\ tous\ Y\ et\ T\ *) 
i1 = 
              Map[ListConvolve[kernelX, #, {mx, mx}] &, 
                s, {2}]; 
 (*\ Convolution\ en\ Y, 
              pour\ tous\ X\ et\ T\ *) 
i2 = 
              Transpose[
                Map[ListConvolve[kernelY, #, {my, my}] &, 
                  Transpose[i1, {1, 3, 2}], {2}], {1, 3, 
                  2}]; 
 (*\ Convolution\ en\ T, 
              pour\ tous\ X\ et\ Y\ *) 
i3 = 
              Transpose[
                Map[ListConvolve[kernelT, #, {mt, mt}] &, 
                  Transpose[i2, {3, 1, 2}], {2}], {2, 3, 1}]; 
	
            i3
)\)];\)\)

Enleve nx,ny,et nt elements à chaque bout des dimensions x,y, et t
(suppose que s contient les dimensions [t,y,x] dans cet ordre)

reduit3D[s_, nx_, ny_, nt_] := 
        Block[{dx, dy, 
            dt}, (
{dt, dy, dx} = 
              Dimensions[s]; 
Take[
              s, {nt + 1, dt - nt}, {ny + 1, dy - ny}, {nx + 1, 
                dx - nx}]
)];

Dérivées spatio-temporelles

On suppose que l'on sait que la séquence contient 9 images,
et l'image du milieu est donc l'image 5.

seq = seq1;

ListDensityPlot[seq[[1]]];

Ecart type est le 1.5 en x, y, et t

sig = 1.5;

Puisque les fitres sont séparables, on utilise un filtre "gaussien" dans deux dimensions et un filtre "dérivée de gaussien" dans la dimension choisie.

ix = 
        convolution3D[seq, 
dkernel[sig, 4], 
          kernel[sig, 4], kernel[sig, 4]];

iy = 
        convolution3D[seq, 
kernel[sig, 4], 
          dkernel[sig, 4], kernel[sig, 4]];

it = 
        convolution3D[seq, 
kernel[sig, 4], kernel[sig, 4], 
          dkernel[sig, 4]];

ix = reduit3D[ix, 4, 4, 4][[1]];

iy = reduit3D[iy, 4, 4, 4][[1]];

it = reduit3D[it, 4, 4, 4][[1]];

maxv = Max[ix, iy, it, -ix, -iy, -it]

ListDensityPlot[ix, 
        PlotRange -> {-maxv, maxv}];
ListDensityPlot[iy, 
        PlotRange -> {-maxv, maxv}];
ListDensityPlot[it, PlotRange -> All];

Flux Normal

On doit éliminer les pixels dont le gradient spatial est trop petit (ici, <0.1) ...

imask = 
        Map[If[Abs[#] < 0.1, 0, 1] &, (ix^2 + iy^2), {2}];

ListDensityPlot[imask];

Le flux normal {un,vn} est défini comme
{un,vn} = -it/(ix^2+iy^2) * {ix,iy}
Mais ici on ajoute dans l'équation le masque de sélection des pixels (imask)

un = (-it*ix*imask)/(ix^2 + 
              iy^2 +  (1 - imask));
vn = (-it*iy*imask)/(ix^2 + 
              iy^2 +  (1 - imask));

ListDensityPlot[un, 
        PlotRange -> All];
ListDensityPlot[vn, PlotRange -> All];

Affichage d'un champs de vitesse

vecteur[val_, {y_, x_}] := {{x, y}, val};

afficheFlux[u_, v_] := afficheFlux[u, v, 1];

afficheFlux[u_, v_, sf_] := 
        Block[{vec, 
            m}, (
  
vec = 
              Transpose[{u, v}, {3, 1, 2}]; 
  
m = 
              Flatten[MapIndexed[vecteur, subSample[vec, 4, 4], {2}], 
                1]; 
ListPlotVectorField[m, 
              ScaleFactor -> sf]
)];

afficheFlux[un, vn];

Pour grouper l'affichage des graphiques...

s0 = Plot[Sin[x], {x, 0, 2  Pi}];

s1 = Plot[x^2, {x, -4, 4}];

s2 = Plot[Tan[x], {x, 0, 1}];

s3 = Plot[Sin[x]/x, {x, 0, 2  Pi}];

Show[GraphicsArray[{{s0, s1}, {s2, s3}}]]