nextupprevious
suivant: Mise en oeuvre de l'algorithmemonter:Restauration d'images par minimisation précédent: Enjeu: débruitage et déconvolution d'images
 

2  A propos de la variation totale


Sous-sections

2.1   Fonctions à variation bornée

La variation totale mesure l'amplitude totale des oscillations d'un signal. En traitement de l'image, elle dépend de la longueur des lignes de niveaux (formule de la co-aire). La variation totale d'une fonction f est ainsi définie par

displaymath24

Si f n'est pas différentiable, on considère cette dérivée au sens général des distributions. On peut alors considérer la variation totale d'une fonction discontinue. On dit que f est à variation totale bornée si

displaymath30

La variation totale dans le cas d'un signal discretisé à la fréquence 1/T permet de mieux comprendre sa signification.

Soit tex2html_wrap_inline35 le signal obtenu par un échantillonnage uniforme de période T. La variation totale discrète se calcule par une différence finie (semblable à une dérivée d'un signal continu):

eqnarray15

où (tex2html_wrap_inline34) désigne les abscisses des extrema locaux de f . Le signal discret sera alors à variation totale bornée si tex2html_wrap_inline41 est majorée par une constante indépendante de la résolution 1/T.

2.2   Un algorithme de minimisation de la variation totale

L'objectif est de minimiser la fonctionnelle suivante:
\begin{displaymath}F(u) = J_1(u) + \beta{\int}_{\Omega}\vert \nabla u \vert \end{displaymath}

où $J_1(u)$ représente le terme d'attache aux données comme dans la section précédente. L'article de Fatemi Rudin Osher, propose une mise en \oeuvre basée sur une descente de gradient. Le calcul du gradient du terme de variation totale donne $div(\frac{\nabla u}{\vert\nabla u\vert})$. Pour une mise en \oeuvre numérique il faut disposer d'un schéma discret de l'algorithme. Dans l'article originel, les auteurs établissent une approximation de $div(\frac{\nabla u}{\vert\nabla u\vert})$ par éléments finis.

Nous présentons dans ce document une variante de cet algorithme proposée par Casas, Kunisch et Pola [1]. Ces derniers proposent une analyse exacte de l'opérateur en introduisant l'ensemble des fonctions discrètes. Toutefois on notera que les deux schémas aboutissent finalement à des résultats identiques

Référence

[1] Casas, Kunisch, Pola, ``Regularization by Functions of Bounded Variation and Applications to Image Enhancement'', Universidad de Cantabria, num 5-1996, 1996  

2.3   Vers un algorithme numérique

2.3.1  Images sur un domaine discret

Une image continue est une correspondance $u:\Omega \rightarrow\mathbb{R}$. Premièrement, on décompose le domaine $\Omega$ en sous-domaines. On pose
\begin{displaymath}\begin{array}{c}a_k = min (x_k : x = (x_j)_{j=1}^n \in \o......ax (x_k: x = (x_j)_{j=1}^n \in \overline{\Omega})\end{array} \end{displaymath}
\begin{displaymath}D = \prod_{k=1}^n (a_k, b_k)\end{displaymath}
Cela définit un plus petit pavé $D$ qui contient $\Omega$.

Ensuite chaque intervalle $[a_k, b_k]$ est sous-divisé en $m_k$ intervalles, partitionant ainsi $D$ en $m$ sous-pavés $(D_j)_{1\leq j \leq m}$, c'est à dire en $m$ pixels lorsqu'il s'agit d'une image:

\begin{displaymath}D_j = \prod_{k=1}^n (a_{k_j}, b_{k_j}) \end{displaymath}
Maintenant on peut encore définir une partition de $\Omega$ en posant
\begin{displaymath}\Omega_j = D_j \cap \Omega\end{displaymath}
Enfin on associe à cette partition l'espace des fonctions constantes par morceaux sur les pixels ${\Omega_j}$. C'est l'ensemble
\begin{displaymath}V_m = \{ u_m = \sum_{j=1}^m u^j \chi_{\Omega_j}: u_j\in \mathbb{R}\} \end{displaymath}

2.3.2   Expression de la variation totale

Le premier résultat important est l'introduction des ``bords discrets'' $ \partial \Omega_j $ dans l'image, si bien que l'on a la propriété en dimension N=2  
\begin{displaymath}\int_{\Omega} \vert\nabla u_m\vert = \sum_{i < j} \vert u^i......tial\Omega_j \vert _{\mathcal{H}_1} \mbox{ où } u_m \inV_m\end{displaymath} (2)


Le calcul de la variation totale d'une image rectangulaire pxn peut se faire en faisant intervenir un vecteur de dimension $r =(n-1)p + n(p-1)$ qui se calcule à l'aide d'une application linéaire $\mathbb{R}^m \rightarrow \mathbb{R}^r $ représentée par une matrice A. Les éléments du vecteur $Au_m$ sont en fait les différences $ (u^i - u^j)$ - on calcule ainsi des différences de pixels sur les bords horizontaux et verticaux. La norme $l^1$ du vecteur $Au_m$ donne le terme de régularisation souhaité. 
 

2.3.3  Note sur la mise en \oeuvre du calcul

Explicitons la forme de la matrice A. On suppose que l'image comporte p lignes et n colonnes et qu'elle est représentée par un vecteur de dimension nxp qui contient les valeurs des pixels de l'image lorsqu'elle est balayée de la gauche vers la droite et de haut en bas. 

\begin{displaymath}A = \left(\begin{array}{ccc}I_p & \otimes & M_n \\M_p & \otimes & I_n \\\end{array} \right)\end{displaymath}

où $I_p$ est la matrice identité p x p, $I_n$ la matrice identité n x n, $\otimes$ représente le produit de Kronecker, $M_n$ est une matrice (n-1)xn 

\begin{displaymath}M_n = \left(\begin{array}{ccccc}1 & -1 & 0 & \cdots & 0......ots & 0 \\0 & \cdots & 0 & 1 & -1 \\\end{array} \right)\end{displaymath}

et $M_p$ est une matrice (p-1) x p, définie de manière identique à $M_n$. 

Coder A sous forme d'une matrice ordinaire, c'est-à-dire indexée par ses lignes et ses colonnes conduirait à une occupation mémoire inutile. En effet la dimension de A est très grande, typiquement de l'ordre de $m^2 = n^4$. Mais il suffit de remarquer que sur une ligne de cette matrice, il n'y a que deux éléments non nuls: (+1) et (-1). La matrice A peut donc être codée comme un tableau de lignes contenant la position de ces deux éléments. 
 

2.4   Méthode de descente de gradient adoptée

On a vu plus haut qu'il fallait prendre le gradient du terme de régularisation. D'après les résultats précédents, on a en fait: 
\begin{displaymath}J_2(u_m) = \varphi(Au_m) \mbox{ où } \varphi \mbox { estsimplement la norme } l^1\end{displaymath}

On peut vérifier que le calcul du gradient de $J_2(u_m)$ donne  

\begin{displaymath}\nabla J_2(u_m) = A^Tsgn(Au_m) \mbox{ où } sgn(.): \mathbb{R} \rightarrow \{-1, +1\} \mbox{ est lafonction signe.}\end{displaymath} (3)


Ceci aboutit donc à un schéma simple de l'algorithme de descente du gradient. 
 
nextupprevious
suivant: Mise en oeuvre de l'algorithmemonter:Restauration d'images par minimisation précédent: Enjeu: débruitage et déconvolution d'images