Filtre de choc de Osher et Rudin

[millie]

Bonjour,

Je vais vous donner l'algorithme pour appliquer un filtre de choc de base.
C'est à dire le filtre de choc de Osher et Rudin qui consiste à déterminer l'état stationnaire de l'équation suivante :

dI/dt = - sign( Dnn(I)) * || grad I ||
(ou Dnn est la dérivée directionnelle selon le gradient et grad est le gradient) (et I l'image).

Il y a d'autres filtres de chocs (qui font par exemple intervenir un arctang ou des complexes, ou encore une convolution avec une gaussienne). Ces autres méthodes sont souvent de meilleur qualité car elles peuvent être appliquées à des images ayant de bruit.

Le filtre de choc de base a tendance à ne pas trop aimer le bruit.

La méthode de base qui est souvent choisi pour résoudre et de calculer :

Code :

1
2
3
4
5
/* I est l'image d'entrée*/
I_0 = I

I_(n+1) = I_n - dt * sign( Dnn(I)) * || grad I ||

Il suffit d'itérer cette équation (une 20 ou 100 aines de fois), mais attention la discrétisation de grad et de Dnn demande quelques précautions.

Tout d'abord :

Voici les fonctions permettant de calculer la différentielle en un point dans un canal particulier (je suis désolé d'avoir choisi une notation C++ mais ça me donne moins de boulot à traduire. J'espère en tout cas que vous comprendrez la sémantique des opérateurs :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
/** différentielle par rapport à x**/

float diffDx(const Image& in, int x, int y, int canal)
{
  int width = in.getWidth();

  int px = x-1;
  int nx = x+1;
  if (px<0)
    px=0;
  if (nx>=width)
    nx=width-1;

  return((in.getPixel(nx,y,canal) - in.getPixel(px,y,canal))/2.0f);
}

float diffDxN(const Image & in, int x, int y, int canal)
{
  int width = in.getWidth();
  if (x+1>=width)
    return 0;

  return(in.getPixel(x+1, y, canal) - in.getPixel(x,y,canal));
}

float diffDxP(const Image & in, int x, int y, int canal)
{
  if(x-1<0)
    return 0;
  return(in.getPixel(x, y, canal) - in.getPixel(x-1,y,canal));
}


/** différentiel par rapport à y **/

float diffDy(const Image& in, int x, int y, int canal)
{
  int height = in.getHeight();

  int py = y-1;
  int ny = y+1;

  if (py<0)
    py=0;
  if (ny>=height)
    ny=height-1;

  return (in.getPixel(x,ny,canal) - in.getPixel(x,py,canal))/2.0f;
}

float diffDyN(const Image & in, int x, int y, int canal)
{
  int height = in.getHeight();
  if (y+1>=height)
    return 0;

  return(in.getPixel(x, y+1, canal) - in.getPixel(x,y,canal));
}

float diffDyP(const Image & in, int x, int y, int canal)
{
  if(y-1<0)
    return 0;
  return(in.getPixel(x, y, canal) - in.getPixel(x,y-1,canal));
}

/** norme du gradient **/

/** version minmod **/
float diffGradNminmod(const Image& in, int x, int y, int canal)
{
  float dmoins= diffDyP(in, x, y, canal);
  float dplus= diffDyN(in, x, y, canal);

  fy= Math::minmod( dmoins, dplus );

  dplus= diffDxP(in, x, y, canal);
  dmoins= diffDxN(in,  x, y, canal);

  fx= Math::minmod( dmoins, dplus );

  return sqrt(fx² + fy²);
}

/** gradient simple **/

float diffGradN(const Image& in, int x, int y, int canal)
{
  return sqrt(diffDx(in, x, y, canal)² + diffDy(in, x, y, canal)²);
}


/** différentielle selon xx, xy et yy **/


float diffDxx(const Image& in, int x, int y, int canal)
{
  int width = in.getWidth();

  int px = x-1;
  int nx = x+1;
  if (px<0)
    px=0;
  if (nx>=width)
    nx=width-1;
return (in.getPixel(nx,y,canal) + in.getPixel(px,y,canal) -2.0f * in.getPixel(x,y, canal));
}


float diffDyy(const Image& in, int x, int y, int canal)
{
  int height = in.getHeight();

  int py = y-1;
  int ny = y+1;

  if (py<0)
    py=0;
  if (ny>=height)
    ny=height-1;

  return (in.getPixel(x,ny,canal) + in.getPixel(x,py,canal) -2.0f * in.getPixel(x,y, canal));
}

float diffDxy(const Image& in, int x, int y, int canal)
{
  int width = in.getWidth();
  int height = in.getHeight();

  int px = x-1;
  int nx = x+1;
  int py = y-1;
  int ny = y+1;
  if (px<0)
    px=0;
  if (nx>=width)
    nx=width-1;
  if (py<0)
    py=0;
  if (ny>=height)
    ny=height-1;

  return 0.25 * (in.getPixel(nx, ny, canal) + in.getPixel(px, py, canal)
                - in.getPixel(px, ny, canal) - in.getPixel(nx, py, canal));
}


/** différentielle suivant le gradient **/

float Millie::diffDnn(const Image& in, int x, int y, int canal)
{
  float normecarre = diffGradN(in, x,y, canal)²;

  float dx = diffDx(in, x, y, canal);
  float dy = diffDy(in, x, y, canal);
  float dxx = diffDxx(in, x, y, canal);
  float dyy = diffDyy(in, x, y, canal);
  float dxy = diffDxy(in, x, y, canal);
  float temp = dx² * dxx + dy² * dyy
               + 2.0 * dx * dy * dxy;
  return (temp/(normecarre+0.001));
}

Avec la fonction minmod qui est définie par (il existe d'autres formulations) :

Code :

1
2
3
4
5
6
7
     float minmod(float a, float b)
      {
        if(a*b>0)
           return min(abs(a),abs(b));
        else
        return 0;
}

Et si vous saviez pas, sign est défini par :

Code :

1
2
3
4
5
6
7
8
  float sign(float a)
      {
        if(a<0)
          return -1.0f;
        if(abs(a)==0.00) /*possibilité de mettre <0.01*/
          return 0.0f;
        return 1.0f;
      }

Et l'algorithme s'écrit donc :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
/* @param out l'image de sortie
 * @param in l'image d'entrée
 * @param dt discrétisation temporelle (souvent proche de 0.1
 * @param nbiter le nombre d'itération
 *
 * erreur si dt<0 ou si nbiter<=0 ou si le nombre de composantes de out et de 
 * in ne sont pas pareil
 */
void miniShock(Image& out, const Image& in,  float dt, int nbiter)
{
  if(out.getNumComponents() != in.getNumComponents())
    Erreur
  if(dt<0)
    Erreur
  if(nbiter<=0)
    Erreur

  Image itempo(in.getWidth(), in.getHeight(), in.getNumComponents());

  out = in;
  itempo = in;

  int n;
  int i,j;

  for(n=0; n<nbiter;n++)
  {
    for(int canal = 0; canal< in.getNumComponents(); canal++)
      for(j=0;j<in.getHeight();j++)
        for(i=0;i<in.getWidth();i++)
        {

          float dnn = diffDnn(out, i, j, canal);

          float second_order = abs(diffDx(out, i, j, canal)) + abs(diffDy(out, i, j, canal));

          if(second_order < 1.0f) /*possibilité de mettre ==0*/
            dnn = diffDxx(out, i, j, canal);

          float grad = diffGradNminmod(out, i, j, canal);

          float value =   dt*(-sign(dnn) *  grad ) + out.getPixel(i,j, canal);

          itempo.setPixel(i, j, canal, value);
        }

    out = itempo;
  }


}

Voici avec : nbiter = 100 et dt = 0.1


Lenna :

Correction :




Ca, c'est une version avec l'application d'une gaussienne (rayon 2, sigma = 10) avec Inn (je n'ai pas noté l'algorithme) : c'est la méthode d'Alvarez et de Mazorra

[factis]

Voici la version matlab. (et optimise pour)

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
function i=miniShock_perf1( i, dt, nbiter)
    if(dt<0)
        warning('dt negatif')
    end
    if(nbiter<=0)
        warning('nbiter negatif')
    end
    tx=size(i,1);
    ty=size(i,2);
    maskup=ones(tx,ty);
    iup=ones(tx+2,ty+2);
    h=waitbar(0,'Init');
    for(n=0:nbiter)
        for(c=1:size(i,3))
        
            iup=[[i(1,1,c),i(1,:,c),i(1,end,c)];[i(:,1,c),i(:,:,c),i(:,end,c)];[i(end,1,c),i(end,:,c),i(end,end,c)]];

            ioxoy=iup(2:end-1,2:end-1);
            ipxoy=iup(2:end-1,3:end);
            inxoy=iup(2:end-1,1:end-2);
            ioxpy=iup(3:end,2:end-1);
            ioxny=iup(1:end-2,2:end-1);
            ipxpy=iup(3:end,3:end);
            inxpy=iup(3:end,1:end-2);
            inxny=iup(1:end-2,1:end-2);
            ipxny=iup(1:end-2,3:end);


            dx =(ipxoy-inxoy)/2;
            dxn=(ipxoy-ioxoy);
            dxp=(ioxoy-inxoy);
            dy =(ioxpy-ioxny)/2;
            dyn=(ioxpy-ioxoy);
            dyp=(ioxoy-ioxny);
            dxy=(ipxpy+inxny-ipxny-inxpy)/4;
            dxx=(ipxoy+inxoy-2*ioxoy);
            dyy=(ioxpy+ioxny-2*ioxoy);


            gradn=sqrt(dxx.*dxx+dyy.*dyy);
            dnn=(dx.*dx.*dxx+dy.*dy.*dyy+2*dx.*dy.*dxy)./(gradn.*gradn+0.001);

            minmodx=((dxn.*dxp)>0).*min(dxn,dxp);
            minmody=((dyn.*dyp)>0).*min(dyn,dyp);
            gradminmod=sqrt(minmodx.*minmodx+minmody.*minmody);


            second_order = abs(dx) + abs(dy);
            dnn=(second_order<1).*dxx+(second_order>=1).*dnn;
            i(:,:,c) = dt*(-sign(dnn).*gradminmod)+i(:,:,c);
            i(:,:,c) = max(min(i(:,:,c),255*maskup),0*maskup); %seuil
        end

        waitbar(0,h,['iter:',num2str(n)])
        imshow(i/max(max(max(i))))
        drawnow;
        pause(0.1)
    end
    close (h)
end

Il faut avouer que pour les floues de 1 a 5 pixels ca marche plutot bien, en tout cas mieux que photoshop.

Par contre mon floue est plus de l'ordre de 400pxl ... oui c'est un jolie floue. Mais vu que l'image fait 18Mpxl (et qu'il y a tres peu de bruit (100ISO) : quasiement que de la quntificqtion JPEG), il doit bien avoir un moyen d'avoir qqc de potable avec 1 ou 2 Mpixl?
(L'origine du floue est un defaut de Mise au point)

J'ai cherche sur google les filtres choc, mais j'en ai pas trouve beacoups (et je pense pas que ca puisse s'appliquer å mon cas assez "extreme"). J'ai trouve 2 logiciels payant (specialises dans le defloutage) que j'ai essaye en demo et les resultats sont vraiment navrant (pire que photoshop ou l'algo ci dessus).

D'ailleur c'est etonnant de voir sur le web autant de gens qui se posent ce probleme et si peu de reponses commercial.

[pseudocode]

Envoyé par factis

Par contre mon floue est plus de l'ordre de 400pxl ... oui c'est un jolie floue. Mais vu que l'image fait 18Mpxl (et qu'il y a tres peu de bruit (100ISO) : quasiement que de la quntificqtion JPEG), il doit bien avoir un moyen d'avoir qqc de potable avec 1 ou 2 Mpixl?
(L'origine du floue est un defaut de Mise au point)

Il faut une approche multi-echelle : construire une pyramide d'image a des tailles différentes, faire le défloutage sur l'image la plus petite, reporter les modifs dans l'image de taille juste supérieure... et recommencer.

[factis]

Bonne idée, mais il risque pas d'avoir des gros bloc un peu trop visibles lors de l'assemblage?

Si on décompose , on traite (donc des contours plus net), et on désassembles, on n'a que des contour plus net: on perd beaucoup de détail.

Pour avoir une idée du problème:

(ok, en fait y a l'air d'avoir un peu plus de bruit que prévu, mais y a toujours la masse de pixels)

Si on souhaite appliquer une déconvolution, quelle est la forme de la PSF? (ça doit être un problème classique les floues de Mise au point?)

[pseudocode]

Ah oui, c'est quand meme un peu flou. La je pense que c'est un peu trop pour un filtre de choc. Il faut passer par des algos de déconvolution et donc effectivement estimer la PSF. Usuellement la PSF pour du flou est modélisée par une gaussienne, mais il y a aussi des techniques pour avoir une meilleure estimation:

google OUT OF FOCUS BLUR PSF