Filtrage de Perona et Malik

Une image peut être modélisée comme une fonction de classe $\mathcal{C}^1$ de $\R^2$ à valeur dans $\R$. Une image numérique est alors une version discrétisée de cette fonction sur un certain rectangle. Le filtrage par équation de la chaleur consiste à faire évoluer l’image pendant une durée déterminée suivant l’EDP de la chaleur. On suppose donc que l’image est une fonction qui varie au cours du temps, notée alors $u(x,y,t)$. L’équation de la chaleur se présente sous la forme :

$$ \left\lbrace \begin{array}{l} \dfrac{\partial u}{\partial t}(x,y,t) - \Delta u(x,y,t) = 0 \text{, pour } t>0 \\\ \\\ u(0,x,y) = u_0(x,y) \end{array} \right. $$

pour $u_0$ une image initiale donnée de taille $M \times N$. $t$ représente le temps. La solution de cette équation en un temps donné correspond à la convolution avec une gaussienne.

Pour résoudre effectivement ce problème, on considère que l’image d’entrée est une version discrétisée de la fonction $u$ au temps $t=0$. Le laplacien $\Delta$ est calculé, à temps $t$ fixé, via des opérateurs discrets de divergence et de gradient : $\Delta u = \mathrm{div} (\nabla u)$. À temps $t$ fixé, le gradient est un champ de vecteur de taille $M \times N$, et chaque vecteur a deux dimensions.

$$ \nabla u_{k,l} = \left( \begin{array}{l} \left\lbrace \begin{array}{rl} u_{k+1,l,1}-u_{k,l,1} & \text{ si } k < N \text{,}\\\ 0 & \text{ si } k = N \end{array} \right. \\\ \\\ \left\lbrace \begin{array}{rl} u_{k,l+1,1}-u_{k,l,1} & \text{ si } l < M \text{,} \\\ 0 & \text{ si } l = M \end{array} \right. \end{array} \right) $$

La divergence est alors l’opérateur conjugué. Il prend un champ de vecteur à deux dimensions en entrée (le résultat du gradient par exemple), et est défini comme suit :

$$ \mathrm{div}(p)_{k,l} = \left\lbrace \begin{array}{rl} p^1_{k,l} - p^1_{k-1,l} & \text{ si } 1<k<N \text{,}\\\ p^1_{k,l} & \text{ si } k = 1 \text{,}\\\ -p^1_{k-1,l} & \text{ si } k = N \end{array} \right. + \left\lbrace \begin{array}{rl} p^2_{k,l}- p^2_{k,l-1} & \text{ si } 1<l<M \text{,}\\\ p^2_{k,l} & \text{ si } l=1 \text{,}\\\ -p^2_{k,l-1} & \text{ si } l=M \end{array} \right. $$

avec $p = (p^1_{k,l}, p^2_{k,l})_{k,l}$ le champ dont on cherche la divergence. $p^1_{k,l}$ (respectivement $p^2_{k,l}$) représente la première (resp. deuxième) coordonnée du champ de vecteur au point de coordonnées $(k,l)$.

Image originale Equation de la chaleur Filtre de Perona-Malik
Original Equation de la chaleur Filtre de Perona-Malik

La méthode de Perona et Malik évite de flouter l’image en préservant les contours. Celle-ci repose sur l’EDP suivante :

$$ \left\lbrace \begin{array}{l} \dfrac{\partial u}{\partial t}(x,y,t) = \mathrm{div} ( f(x,y) . \nabla u (x,y,t)) \text{, pour } t>0\\\ \\\ u(0,x,y) = u_0(x,y) \end{array} \right. $$

avec $f$ une fonction de $\R^2$ à valeur dans $[0,1]$ a valeur faible près des contours, c’est à dire avec une valeur faible près des zones de l’image présentant un gradient élevé.

La fonction $f$ se calcule à partir de la fonction initiale et permet de limiter l’action du filtrage (floutage) au niveau des zones de contours. En pratique, c’est un tableau de coefficients de la taille de l’image. L’équation de la chaleur correspond à celle de Perona-Malik où $f \equiv 1$.

Afin de ne pas flouter uniformément l’image, on réduit l’action de l’opérateur laplacien au niveau des zones de contours en multipliant le gradient par une fonction ayant une valeur proche de 1 aux endroits où il n’y a pas de contour (norme du gradient faible) et ayant une valeur proche de 0 aux endroits où il y a des contours (norme du gradient élevée). Une telle fonction peut être choisie de la forme $f(x,y) = \exp ( - \Vert \nabla F(x,y) \Vert^2 )$ avec $F$ une convolution de l’image initiale avec une gaussienne. Cette convolution vise à élargir la zone de coefficients faibles afin de mieux préserver les contours.

  • Questions :
    1. Implémenter les gradients et divergence de l’image.
    2. Implémenter le schéma d’Euler sur l’équation de la chaleur ($U_{n+1} = U_n + \delta_t \cdot \Delta U_n$, prendre un pas de temps $\delta_t = 0.1$)
    3. Implémenter la convolution de l’image avec une gaussienne en prenant soin de considérer que les bords de l’image soient symétriques (utiliser scipy.signal.convolve2d). Construire la fonction $f$ et la visualiser.
    4. Implémenter le schéma d’Euler correspondant à l’équation de Perona et Malik avec un pas de temps de 0.2 et 200 itérations.