[Image] Filtre de Malik et Perona

**millie** · 03/05/2007, 12h09

Le filtre de malik est perona est un des premiers filtres basé sur des équations différentielles permettant une diffusion anisotropique.

Au début, certains filtres avaient été proposés (basé sur l'équation de la chaleur), mais avait l'inconvénient de ne pas garder les contours.

Malik et Perona propose de déterminer l'état stationnaire de l'équation :

du/dt = div(g(|grad u|) grad u) où u est l'image.
g est une fonction qui peut prendre plusieurs forme.

Celle que j'ai prise fonctionne bien et a été proposé par malik et perona :
g : x -> lambda²/(lambda² + x²)

J'ai trouvé sur le papier ici :http://www.cs.berkeley.edu/~malik/papers/MP-aniso.pdf assez efficace, simple et rapide.

Déjà la simple fonction qui calcule g (on note lambda2 = lambda²) :

Code C++ :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
  float PeronaMalikFlow::computeDiffusion(float v)
  {
    return (lambda2)/((lambda2) + Math::square(v));
  }

(Math::square calcule le nombre au carré)

Voici le code qui permet de calculer le flux en un point (x,y,canal) précis :

Code C++ :

Sélectionner tout - Visualiser dans une fenêtre à part

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
float PeronaMalikFlow::compute(const Image & image, int x, int y, int canal)
  {
    float current = image.getPixel(x,y, canal);
    int width = image.getWidth();
    int height = image.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;
 
    float ixp = image.getPixel(px, y, canal);
    float ixn = image.getPixel(nx, y, canal);
    float iyp = image.getPixel(x, ny, canal);
    float iyn = image.getPixel(x, py, canal);
 
    float diffxn = computeDiffusion(current - ixn);
    float diffxp = computeDiffusion(current - ixp);
    float diffyn = computeDiffusion(current - iyn);
    float diffyp = computeDiffusion(current - ixp);
 
 
    float delta = (diffxn * (ixn - current))
                  + (diffxp * (ixp - current))
                  + (diffyp * (iyp - current))
                  + (diffyn * (iyn - current));
 
 
    return delta;
  }

Et voici la résolution :

Code C++ :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
 
  /**
   * Résout une équation de type perona et malik suivant un diffuseur particulier
   *
   * @param out l'image de sortie
   * @param in l'image d'entrée
   * @param dt la discrétisation spatiale
   * @param nbiter le nombre d'itération pour la résolution
 
   *
   * @throw Millie::IllegalArgument si le nombre de composantes de out et
   *   de in ne sont pas pareils
   * @throw Millie::IllegalArgument si dt<=0
   * @throw Millie::IllegalArgument si nbiter<0
   */
void Millie::peronaMalikOperator(Image & out, const Image & in,
                                float dt,
                                 int nbiter)
{
  if(out.getNumComponents() != in.getNumComponents())
    throw IllegalArgument("peronaMalikOperator");
 
  if(dt<0)
    throw IllegalArgument("peronaMalikOperator negative value");
 
  if(nbiter<=0)
    throw IllegalArgument("peronaMalikOperator : <= iteration");
 
  /*image temporaire*/
  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 value = dt * ( compute(out, i, j, canal)) //calcul de la diffusion
                                    + out.getPixel(i,j, canal);
          itempo.setPixel(i, j, canal, value);
        }
    out = itempo;
 
  }
}

A noter que dans la fonction : compute, on ne prend en compte que les pixels au nord, est, sud et ouest.
Il existe une version plus sympatique qui prend en compte les 4 autres points (NE, SE, NO, SO) :

Code C++ :

Sélectionner tout - Visualiser dans une fenêtre à part

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
float NeighbourPeronaMalikFlow::compute(const Image& image,
                     int x,
                     int y,
                     int canal)
{
    float current = image.getPixel(x,y, canal);
    int width = image.getWidth();
    int height = image.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;
 
    float ixp = image.getPixel(px, y, canal);
    float ixn = image.getPixel(nx, y, canal);
    float iyp = image.getPixel(x, ny, canal);
    float iyn = image.getPixel(x, py, canal);
 
    float diffxn = computeDiffusion(current- ixn);
    float diffxp = computeDiffusion(current- ixp);
    float diffyn = computeDiffusion(current- iyn);
    float diffyp = computeDiffusion(current- ixp);
 
    float iNE = image.getPixel(nx, py, canal);
    float iSW = image.getPixel(px, ny, canal);
    float iNW = image.getPixel(px, py, canal);
    float iSE = image.getPixel(nx ,ny, canal);
    float diffNE = computeDiffusion(current  - iNE);
    float diffSW = computeDiffusion(current- iSW);
    float diffNW = computeDiffusion(current- iNW);
    float diffSE = computeDiffusion(current- iSE);
 
 
    float delta =   diffxn * (ixn - current)
                  + diffxp * (ixp - current)
                  + diffyp * (iyp - current)
                  + diffyn * (iyn - current)
                  +
                  0.5f * (
                 diffNE * (iNE- current)
              +  diffSW * (iSW - current)
              +  diffSE * (iSE - current)
              +  diffNW * (iNW -current)
             );
 
 
    return delta;
 
}

Voici par exemple ce que ça donne :
dt = 0.1, lambda = 5, nbiter = 20

Image bruitée :

Image corrigée :

**pseudocode** · 03/05/2007, 16h12

Celle que j'ai prise fonctionne bien et a été proposé par malik et perona :
g : x -> lambda²/(lambda² + x²)

Les fonctions de controle de continuité de malik et perona, c'est pas censé être:

f1(x) = Math.Exp( - (x/k)² )

f2(x) = 1 / (1 + (x/k)²)

non ? Me trompe-je ?

**millie** · 03/05/2007, 16h25

Envoyé par pseudocode

Les fonctions de controle de continuité de malik et perona, c'est pas censé être:

f1(x) = Math.Exp( - (x/k)² )

f2(x) = 1 / (1 + (x/k)²)

non ? Me trompe-je ?

f2, c'est pareil que la mienne

(suffit de multiplier par k²/k²) (k = lambda)

f1, c'est effectivement une autre qui a été indiqué dans le papier, mais je trouvais qu'elle donnait de moins bon résultat... (enfin, ça prend une ligne à changer cela dit)

**pseudocode** · 03/05/2007, 17h05

Envoyé par millie

f2, c'est pareil que la mienne

(suffit de multiplier par k²/k²) (k = lambda)

f1, c'est effectivement une autre qui a été indiqué dans le papier, mais je trouvais qu'elle donnait de moins bon résultat... (enfin, ça prend une ligne à changer cela dit)

The function (f1) privileges wider regions over smaller ones, while the function (f2) privileges high-contrast edges over low-contrast ones.

Sinon je te propose celle la:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
 
 
                               k * exp( - ||Gradient(Ixy)|| / M )
f( Gradient(Ixy) ) = ----------------------------------------------------
                       max( exp( Laplacian(Ixy) ) , 1+||Gradient(Ixy)|| )
 
avec,
 
M = max ( ||Gradient(Ixy)|| )  
    Vxy
 
Vxy = voisinage du pixel xy

la source --> http://www.ee.columbia.edu/~qibin/pa...004_icip_2.pdf

**rouliane** · 20/02/2008, 21h08

Bonsoir,

Je regardais un peu l'algo de Millie, et je me demandais ce qu'était l'entier canal ?

**millie** · 20/02/2008, 21h10

Envoyé par rouliane

Bonsoir,

Je regardais un peu l'algo de Millie, et je me demandais ce qu'était l'entier canal ?

C'est juste qu'une image à plusieurs canaux, par exemple, R, G et B. Pour avoir un algorithme qui fonctionne avec tout type d'image (NB, RGB), j'avais utilisé un canal et son numéro

**rouliane** · 20/02/2008, 21h26

D'accord, merci.

**rouliane** · 26/02/2008, 18h07

Bonjour,

J'essaye de programmer ce filtre avec Matlab, mais je trouve des résultats sensiblement différents ( j'ai besoin de plus d'itérations pour converger )

Je me demande déjà pourquoi, dans le calcul de g(|gradu|) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

float diffxn = computeDiffusion(current - ixn);

Vous ne faites jamais intervenir le gradient justement.
Pour ma part j'ai discrétisé le gradient, ce qui fait que le coefficient diffxn ( qui est le coefficient cN dans le document de Berkeley ) est :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

cN=1/( 1+(normgrad(sortie, i-1,j)/K).^2 );

où normgrad est donnée par :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
function sortie = normgrad(image,i,j)
 
Dx=(image(i+1,j)-image(i-1,j))/2;
Dy=(image(i,j+1)-image(i,j-1))/2;
sortie =sqrt( Dx^2+Dy^2);

Pourquoi n'avez vous pas calculé le gradient ?

Voici ce que j'obtiens pour 40 itérations :

**pseudocode** · 26/02/2008, 18h35

Envoyé par rouliane

Pour ma part j'ai discrétisé le gradient, ce qui fait que le coefficient diffxn ( qui est le coefficient cN dans le document de Berkeley ) est :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

cN=1/( 1+(normgrad(sortie, i-1,j)/K).^2 );

Pourquoi n'avez vous pas calculé le gradient ?

C'est une simplification usuelle dans la discretisation du gradient.

Au lieu de calculer:

Cn * Dn = g( || gradient(i+0.5,j) || ) * ( I(i-1,j) - I(i,j) )

on utilise:

Cn * Dn = g( |I(i-1,j) - I(i,j)| ) * ( I(i-1,j) - I(i,j) )

C'est vrai que c'est moins précis mais ça donne des suffisamment bons résultats:

**rouliane** · 26/02/2008, 18h48

Ok, merci beaucoup, je ne connaissais pas cette approximation.

**pseudocode** · 26/02/2008, 18h53

Envoyé par rouliane

Ok, merci beaucoup, je ne connaissais pas cette approximation.

Elle n'a rien d'officielle (au sens mathématique). C'est plus pour accélerer les calculs qu'autre chose.

**Lost in** · 31/05/2008, 19h16

Salut,

d'après ce que j'ai compris le filtrage anisotrope n'affecte pas les contours "francs" de l'image.

J'ai remarqué que plus le nombre d'itérations augmente les contours "francs" commencent à être diffusés

Je n'ai pas trouvé d'explication alors pourriez vous m'aider SVP

**pseudocode** · 31/05/2008, 20h38

Envoyé par Lost in

Salut,

d'après ce que j'ai compris le filtrage anisotrope n'affecte pas les contours "francs" de l'image.

J'ai remarqué que plus le nombre d'itérations augmente les contours "francs" commencent à être diffusés

Je n'ai pas trouvé d'explication alors pourriez vous m'aider SVP

C'est parce que la fonction de diffusion n'est jamais nulle, même pour les fortes différences d'intensité. Il y a donc toujours un mélange entre les pixels, et il devient visible avec beaucoup d'iterations.

Un moyen d'eviter cela c'est de forcer la fonction de diffusion à zéro au dela d'un seuil.

**Lost in** · 31/05/2008, 22h23

Mille mercis pseudocode pour ton aide

**Lost in** · 16/06/2008, 20h56

Salut,

Décidément mes questions ne s'arrêtent jamais

Je suis à la recherche de critère d'arrêt pour le filtrage anisotrope, dans la littérature on en trouve mais j'aimerais bien savoir si vous avez déjà testé un et est ce qu'il a donné de bon résultats

**pseudocode** · 16/06/2008, 22h27

Je n'ai jamais trouvé de critère universel, mais j'aime bien l'estimateur local de ROAD (couplé à un seuil pour savoir si on continue ou pas les iterations)

**pseudocode** · 16/06/2008, 23h44

Envoyé par Lost in

puisque tu l'utilises je suis sûre que c'est le meilleur.

Peux tu s'il te plait m'expliquer le principe ?

Ca consiste a estimer l'homogénéité de la zone autour d'un pixel.

On prend un voisinage du pixel P à estimer
On calcule D_i = abs(P_i - P) pour chaque P_i dans le voisinage
on trie la liste des D_i
on somme les "n" premiers termes de la liste triée (généralement "n" = moitié de la liste)

la valeur obtenue est l'estimateur de ROAD du pixel P. Plus cette valeur est faible, plus il est probable que "P" fasse partie d'une zone homogène du voisinage. Au contraire, une valeur haute indique que P est probablement un outlier.

**Lost in** · 17/06/2008, 23h37

Bonsoir,

Merci pour l'explication pseudocode

Une autre question

:

J'ai remarqué que tous les critères d'arrêt utilisent un seuil, par exemple quand on utilise la variance il faut s'arrêter lorsque le rapport
varience(itération n)/variance(image initiale) est inférieur à un seuil r

Comment choisir ce seuil r ?

**Lost in** · 13/07/2008, 01h40

Excusez le nombre important de mes questions, mais je n'y peux rien car
le sujet m'intéresse beaucoup

Envoyé par pseudocode

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
 
 
                               k * exp( - ||Gradient(Ixy)|| / M )
f( Gradient(Ixy) ) = ----------------------------------------------------
                       max( exp( Laplacian(Ixy) ) , 1+||Gradient(Ixy)|| )
 
avec,
 
M = max ( ||Gradient(Ixy)|| )  
    Vxy
 
Vxy = voisinage du pixel xy

Il y a un petit problème lorsque M=0 (division par zéro) dans ce cas qu'est ce qu'on doit faire

Dans l'attente de vous lire je vous dis à bientôt !!

**pseudocode** · 14/07/2008, 13h22

Envoyé par Lost in

Il y a un petit problème lorsque M=0 (division par zéro) dans ce cas qu'est ce qu'on doit faire

on doit réflechir...

Si M=0, ca veut dire que le plus grand gradient (en norme) du voisinage est nul.

Donc ca veut dire que tous les gradients (en norme) du voisinage sont nuls.

Ce qui signifie que le pixel considéré est au milieu d'une zone completement uniforme (tous les pixels ont la meme valeur).

Donc quelque soit la valeur du coefficient de diffusion, la valeur du pixel filtré restera inchangée (barycentre de n valeurs identiques = la valeur).

On peut donc forcer la valeur de f() à n'importe quelle valeur.

En particulier à 0, ce qui revient a ne pas faire de diffusion pour ce pixel. En terme d'implémentation, ca veut dire qu'on peut faire un court-circuit et passer au pixel suivant.

[Image] Filtre de Malik et Perona

Contribuez

Discussions similaires

Partager

Partager