Détection de contours par recuit simulé

Hanna Martinsson et Mathieu Renaudat< - dea MVA 2003-2004/STRONG>





 

Table des matières



Introduction

La détection de contours donne et a donné lieu à de nombreuses recherches. Certaines approches utilisent la minimisation d'une énergie associée à un contour. Ces méthodes ont d'après les auteurs de l'article [2] le défaut de souvent dépendre de l'initialisation de ces contours. En effet, l'énergie qu'il faut minimiser admet souvent des minima locaux. Ainsi, le recuit simulé paraît être un outil intéressant pour s'affranchir des initialisation laborieuse de contours. Suite à une initialisation plus ou moins correcte, le contour actif évolue vers un minimum local, alors que l'on souhaite trouver le contour qui minimise de manière globale l'énergie. De plus, l'utilisation des B-splines s'avère être un choix intéressant pour modéliser les contours. L'écart entre les points de contrôle est un paramètre que l'on peut fixé a priori. Ainsi le terme d'énergie de régularisation n'a plus lieu d'être. L'utilisation de ces deux méthodes a déjà été adoptée dans le passé. L'aspect nouveau dans l'article [2] vient d'un échantilloneur de Gibbs continu. Ainsi, l'échantillonnage adaptatif, l'algorithme ``max braketing'' ou encore l'utilisation d'un modèle a priori de la densité de probabilité conditionnelle utilise une hypothèse forte qui n'est pas explicitée: les échantillons peuvent être triés. Ainsi, ce type d'algorithme permet par exemple de traiter le problème de la reconnaissance vocale (formant tracking) puisqu'il s'agit de détecter les contours d'un signal temps-fréquence. Dans une première partie, nous présentons l'intérêt des B-splines pour la détection de contours. Puis nous décrivons le principe du recuit simulé, pour ensuite expliquer son utilisation dans le cadre de la détection de contours par les B-splines et enfin, nous présentons des résultats obtenus par les auteurs de l'article [2] notamment pour l'application reconnaissance vocale.

B-Splines

En traitement d'image, on souhaite souvent détecter les contours d'une forme dans l'image. Dans le cas des contours actifs ou snakes, on se donne un potentiel image $ f(x,y)$ qui est minimum localement au niveau des contours. Et on fait évoluer un contour actif au grés de plusieurs forces: La force du potentiel, $ -\nabla f$, ainsi que les forces de régularisation.

Ces deux ``forces'' donnent lieu à deux énergies: une énergie potentielle donnée par l'image (c'est le terme de fidélité aux contours) et une énergie qui fixe l'élasticité et la raideur du contour.

Cette énergie est calculée le long de la courbe Lorsque l'on utilise des B-splines, le terme de régularisation devient implicite via le nombre de points de contrôle que l'on utilise pour construire notre B-spline. En effet si on utilise peu de points de contrôle alors la B-spline est très lisse. Par contre lorsque l'on utilise beaucoup de points de contrôle on peut faire apparaître plus de détails. On a représenté les deux cas d'un contour avec beaucoup de points de contrôle et d'un contour avec peu de points de contrôle sur les figures (1) et (2).

Figure 1: B-spline avec 50 points de contrôle.
\includegraphics[height=6cm]{bspline50.eps}
La figure (1) montre une B-spline avec 50 points de contrôle.
Figure 2: B-spline avec 5 points de contrôle.
\includegraphics[height=6cm]{bspline5.eps}
Ayant dans la figure (2) limité le nombre de points de contrôle à 5, on voit bien que la B-spline qui en résulte comporte moins de détails que la précédante. Choisir peu de points de contrôle s'avère donc être un paramètre de lissage d'un contour.

L'énergie qu'il faut minimiser ne doit donc pas forcément faire apparaître un terme de lissage. Dans l'article [3], une comparaison est faite entre les snakes admettant une énergie interne pour (pour la régularité) et les cubic B-splines ou l'écart entre les points de contrôle est utilisé comme un paramètre pour fixer l'élasticité de la courbe.

Recuit simulé

Dans le domaine du traitment des images, nous sommes souvent amenés à minimiser une énergie. La fonction d'énergie présente en général un grand nombre de minima locaux, ce qui complique la tâche de trouver un minimum global. Prenons par exemple la méthode classique de descente de gradient: Elle consiste à, pour chaque itération, faire un pas dans la direction de la plus grande pente vers les basses valeurs de la fonction. Si, pour un état initial donné, l'algorithme arrive après quelques itérations à un minimum local, il n'y a plus de direction qui fait diminuer l'énergie est l'algorithme se trouve alors piégé.

L'idée de l'algorithme de recuit simulé est de donner la possibilité de sortir des minima locaux en autorisant l'algorithme à prendre des pas dans une mauvaise direction avec une faible probabilité.

Le paramètre de température

Le terme recuit vient de l'opération qui consiste à faire fondre, puis laisser refroidir lentement un métal. Cette procédure amène le métal à son état de plus basse énergie, un état très ordonné, tandis qu'un refroidissement trop rapide pourrait laisser le métal piégé dans un état peu favorable, car trop énergétique.

Le recuit simulé cherche donc à imiter cela en introduisant un paramètre global, $ T$, appelé température. A chaque instant, le prochain pas se fera alors au hasard, selon une distribution de probabilité qui dépend de la température. Nous introduisons la distribution de Gibbs avec paramètre de température, définie de la manière suivante

$\displaystyle P_T(X=x) = {1 \over Z(T)} e^{{-U(x) \over T}} \textrm{, avec } Z(T) = \sum e^{{-U(x) \over T}}$

La distribution de Gibbs a la bonne propriété de se comporter de la même manière que le modèle physique à températures extrèmes. Pour $ T$ très grande, $ e^{{-U(x) \over T}}$ se rapproche de $ 1$ et $ P_T$ converge vers la probabilité uniforme sur le domaine. Toute direction a donc la même probabilité d'être prise. Pour $ T$ proche de 0, on montre que $ P_T$ est uniformément distribuée sur les minima globaux. Les deux extrèmes correspondent à un algorithme qui cherche le minimum complètement au hasard et à l'algorithme de descente de gradient respectivement.

L'idée d'introduire un paramètre de température à été proposée par Kirkpatrick dans [4].

Echantillonneur de Gibbs

Afin de réaliser l'algorithme du recuit simulé, nous avons besoin d'une méthode d'échantillonner selon la loi de Gibbs. L'échanillonneur de Gibbs, proposé par Geman et Geman dans [1], procède en créant itérativement une suite de B-splines (toutes cherchant à modéliser un même contour). Après un nombre d'itérations suffissamment grand le contour créé est une réalisation de la loi de Gibbs. A l'étape $ n$, l'algorithme consiste à effectuer les opérations suivantes:

$ \circ$
On choisit un point de contrôle du spline,
$ \circ$
Pour ce point, on calcule la probabilité locale conditionnelle,
$ \circ$
Les coordonnées du point sont mis à jour par un tirage aléatoire selon la loi trouvée.

Dans la forme originale de l'échantillonneur de Gibbs, on traite des images considérées comme étant des champs de Markov et la probabilité conditionnelle locale d'un pixel ne dépend alors uniquement de la configuration du voisinage. Ici, le calcul de la probabilité conditionnelle locale est plus compliqué et sera traité ultérieurement.

Si le choix du point de contrôle à l'étape $ n$ se fait de manière à ce que tous les points soient balayés un grand nombre de fois, le splines obtenu après une infinité d'itérations est une réalisation de la loi de Gibbs. La convergeance de l'algorithme a été démontré par Winkler (1995). En pratique, nous allons considérer que l'on a convergé lorsque les changements sont petits.

L'algorithme du recuit simulé

Nous disposons maintenant des outils nécessaires pour présenter l'algorithme du recuit simulé proposé par Geman et Geman dans [1]. Il se présente comme ceci:

$ \circ$
On choisit une température initiale $ T_0$,
$ \circ$
On choisit une B-spline initial,
$ \circ$
Pour $ n = 1$ ..
$ \triangleright$
A partir de la configuration précédente $ x(n-1)$, on utilise l'échantillonneur de Gibbs pour obtenir une réalisation de la loi de Gibbs d'énergie $ U(x) / T(n)$,
$ \triangleright$
On diminue la température et incrémente $ n$,
$ \circ$
On arrète lorsque les changements sont petits.

La température initiale ne doit pas être choisi trop faible, pour assurer la convergeance, mais pas non plus trop élevée, pour assurer un temps de calcul raisonnable. Puisque l'algorithme ne reste pas bloqué dans des minima locaux, la B-spline initiale peut être choisie quelconque.

En général, on fait un balayage de tous les points de contrôle à chaque niveau de température.

Si la température décroît suffissamment lentement, on a convergeance de l'algorithme. En théorie, il faut une décroissance logarithmique, mais en pratique on utilise des décroissances géométriques avec des bons résultats.

L'échantillonnage conditionnelle

Dans le cas classique d'un champs de Markov discret, il est facile de calculer l'énergie associée à chaque valeur possible d'un pixel et d'en déduire la probabilité conditionnelle locale. Or, appliqué à l'éstimation des contours par des B-splines, nous avons tellement de valeurs, ou emplacements, possibles pour chaque point de contrôle que l'on peut considérer le problème comme étant continue. Puisque chaque calcul d'énergie coute cher en temps de calcul, nous devons faire un choix rusé d'emplacements pour lesquels les calculs seront faits. Il n'est pas optimal d'utiliser un quadrillage uniforme d'emplacements, car à température basse la distribution de probabilité conditionnelle est tres localisée est la majorité des points retenus ne serviront à rien à cause de leur trop grande distance du centre de la distribution. Ce problème a été traité par Jamieson dans [2], qui a choisi d'utiliser des méthodes adaptatives pour y remédier. Ces méthodes consiste à prendre les énergies calculées précédemment comme le point de départ, puis d'identifier la région de plus grande incertitude pour y placer d'avantage de points.

B-Splines et le recuit simulé

Les auteurs de l'article [2] expliquent que l'efficacité de l'algorithme qu'ils proposent réside dans l'échantillonnage qui est un échatillonnage proche de l'échantillonnage de Gibbs discret qui engendre des temps de calcul importants. Ici l'échantillonneur utilisé est dit continu. Il s'agit de positionner de calculer des probabilités conditionnelles tout comme dans le cas de l'échantilloneur de Gibbs discret mais en utilisant des méthodes pour limiter le nombre d'échantillonnage coûteux en temps de calcul. Pour cela, il est utilisé des techniques dites adaptative ou encore des modèles a priori pour estimer au mieux cette probabilité conditionnelle avec un minimum de ``mesures''. Par exemple, le cas du ``Formant tracking'' permet de situer le problème en problème 1D. Ainsi, on peut estimer la probabilité conditionnelle par une courbe comme dans la figure (3).
Figure 3: Probabilité conditionnelle estimée par différentes techniques.
\includegraphics[height=6cm]{probacond.eps}
La zone grisée montre la probabilité réelle, en le trait entier l'estimation par un échantillonnage régulier, et le trait pointillé une estimation avec un algorithme adaptatif. Avec le même nombre de mesures faites, l'algorithme adaptatif est plus précis.

L'élément clef qui n'est pas dit explicitement dans cette article est qu'il faut pouvoir classer les paramètres, en l'occurrence les points de contrôle de la B-spline, ce qui n'est pas toujours le cas dans une problématique de traitement d'image pur et dur. Il faut dans le cas contraire se ramener à un cas en 1D. C'est cela qui est fait dans tous les exemples proposés.

Dans l'article [2], il est présenté un algorithme qui permet de se figurer en quoi consiste cet échantillonnage: l'algorithme de Max-Bracketing, dont le principe est d'identifier l'endroit où se trouve l'incertitude maximale. On suppose que les mesures au points $ x_i$ ont été faites et on considère alors le bloc qui a l'aire la plus élevée en termes de probabilité conditionnelle,

$\displaystyle \arg_i \max_i {1 \over 2} (x_{i+1} - x_{i-1}) e^{-H(x) / T},$

puis on fait de nouvelles mesures à l'intérieur de ce bloc afin de concentrer les efforts là où ca nous intéresse le plus. Il s'agit de faire une mesure en fonction des mesures qui ont été faites précédemment pour faire un échantillonnage adaptatif.

Une autre approche consiste à supposer un modèle de la densité de probabilité a priori.

Les résultats sont bien meilleurs lorsque l'on utilise l'un de ces échantillonneurs adaptatifs plutôt qu'un échantillonnage régulier et linéaire. Les figures (4) et (5) présentent des résultats de convergence de recuits simulés. Evidemment lorsque l'on arrête l'algorithme lorsque la température est faible comme sur la figure 4, l'algorithme a presque convergé vers le résultat souhaité, c'est à dire une courbe tracée à la main. Par contre, lorsque l'algorithme est arrêté à haute température, la convergence n'est pas finie, et les courbes obtenues sont relativement loin du résultat souhaité.

Figure 4: Convergence des B-splines à basse température.
\includegraphics[height=6cm]{convabasseT.eps}

Le résultat de la convergence des B-splines vers les contours dessinés à la main lorsque la température est basse est satisfaisante.

Figure 5: Convergence des B-splines à haute température.
\includegraphics[height=6cm]{convahauteT.eps}

Le résultat de la convergence des B-splines vers les contours dessinés à la main lorsque la température est haute est moins satisfaisante que lorsque l'on considère les résultats obtenus avec une température basse.

Implémentation

Algorithme

Afin de démonstrer le fonctionnement du méthode, nous en avons fait une implémentation en Matlab. Nous nous sommes inspirés d'une implémentation faite par Jamieson, évoquée dans [2]. Pour simplifier, nous cherchons un contour qui traverse l'image dans toute sa largeur (ou sa longueur). Ainsi, nous pouvons paramétriser la courbe par son abscisse en x (ou y) et la détéction du contour se ramène à la recherche pour tout x de son ordonnée. Cela revient donc à considerer un problème d'échantillonnage en une seule dimension.

L'initialisation du contour consiste à lancer une B-spline horizontal avec des points de contrôle équidistribués suivant l'axe du milieu de l'image. Pour chaque valeur de température donnée, un ordre de visite des points de contrôle est décidée au hasard. Lors de la visite d'un point, l'algorithme lui attribué un certain nombre d'emplacements candidats en y et calcule pour chacun d'entre eux l'énergie qui y est associée. Cette énergie à deux parties, l'énergie externe liée à un potentiel d'énergie définie par un filtre "laplacian of gaussian" appliqué à l'image, et l'énergie interne liée à la régularité de la courbe. Basé sur ces énergies, la probabilité locale conditionnelle est calculée et l'algorithme effectue ensuite un tirage aléatoire selon la loi trouvée. La nouvelle coordonnée du point est alors donnée par le candidat retenu. Pour simplifier l'implémentation, l'algorithme n'est pas adaptatif, mais fait plutôt un échantillonnage linéair dans le voisinage du point de contrôle en question.

Code source

La fonction principale a pour syntax
[controlPoints, contour] = annealer(imageFile,TempVect,sampleNum,spacing)
où imageFile est le nom du fichier contenant l'image, TempVect le vecteur (décroissant) de températures, sampleNum le nombre d'emplacements à essayer dans chaque étape et spacing le variable d'espacement, c'est-à-dire le nombre de pixels entre les points de contrôle. Les variables de sortie contiennent les vecteurs des points de contrôle et du contour respectivement.

Les fonctions auxiliaires appellées pendant la phase d'initialisation sont
U = potentialEnergy(I),
qui calcule le potentiel de l'image, défini par le filtre "laplacian of gaussian", ainsi que
contour = formBspline(controlPoints, spacing),
laquelle, à partir des points de contrôle, calcule la B-spline qui donne le contour estimé. Cette dernière appelle à son tour la fonction
bVal = bSplineFunc(val),
oû est définie la fonction de base de la B-spline.

Après l'initialisation, le programme execute, pour chaque valeur de température donnée et pour chacun des points de contrôle, la fonction
[controlPoints, contour] = gibbsSampler(U,controlPoints,contour,spacing,Temp,iCP,nSamples),
qui effectue l'échantillonnage de Gibbs. En dehors des fonctions mentionnées précécemment, cette fonction appelle aussi
H = computeH(U,controlPoints,spacing,iCP,Cval),
pour le calcul l'énergie du contour pour une valeur particulière d'un point de contrôle donné.

Résultats

Puisque notre programme cherche un contour qui traverse l'image, il est bien adapté à la recherche de l'horizon ou d'une ligne de crête. Nous l'avons ainsi essayé sur une image avec un contour net entre la terre et le ciel. Le résultat est présenté dans la figure (6). Avec un vecteur de températures qui présente une décroissance logarithmique, l'algorithme finit par trouver une bonne estimation du contour, tracée en rouge par dessus de l'image originale.

Figure 6: Contour tracé avec notre implémentation.

Conclusion

En conclusion, les résultats sont prometteurs. Cette méthode permet de s'affranchir des initialisations fastidieuses ce qui est tout l'intérêt du recuit simulé. Elle utilise un outil de représentation très flexible: la B-spline. Cependant ce qui demande du travail est de se replacer dans un contexte bien particulier qui permet de classer les points de contrôle, c'est-à-dire de passer d'un problème 2D à un problème 1D.

Bibliographie

1
S. Geman and D. Geman, Stochastic Relaxation, Gibbs Distribution and the Bayesian Restoration of Images, IEEE Transactions on Pattern Analysis and Machine Intelligence 6 (1984), no. 6, 721-741.

2
M. Jamieson, P. Fieguth and L. J. Lee, Parametric Contour Estimation by Simulated Annealing, ICIP03, vol. 3, 2003, pp. 449-452.

3
P. Brigger, J. Hoeg and M. Unse, A Flexible Tool for Parametric Contour Detection, IEEE Transactions on Image Processing 9 (2000), 1484-1496.

4
S. Kirkpatrick, C. D. Gellatt and M. P. Vecchi, Optimization by Simulated Annealing, Science (1983), no. 220, 671-680.