Recalage d’images, PIV et corrélation d’images — Les bases théoriques

Posté par  . Édité par ZeroHeure, Davy Defaud, Benoît Sibaud, NeoX et palm123. Modéré par ZeroHeure. Licence CC By‑SA.
62
4
juil.
2017
Science

Le recalage d’images est utilisé dans la communauté de l’analyse d’images médicales depuis très longtemps (Barnea & Silverman, 1972, Ledbetter et al. 1979), on peut même remonter plus loin dans le temps avec les travaux de Sidney Bertram, 1963 en fournissant tous les outils via des descriptions de circuits analogiques (vraiment impressionnant pour l’époque).
En commençant à travailler sur des images issues d’IRM, j’ai été amené à faire une bibliographie sur les outils utilisés dans l’analyse d’images médicales. J’ai constaté que leurs techniques de recalage d’images étaient très similaires à celles utilisées dans les différents domaines de la mécanique. Après beaucoup de bibliographie, je suis parvenu à situer ce qui est fait en mécanique dans les cadres proposés en analyse d’images médicales.

Cette dépêche est l’occasion pour moi de présenter ces techniques dans le détail.


La prochaine dépêche présentera les logiciels utilisés dans les différentes communautés scientifiques.

Sommaire

Introduction

Dans la recherche que nous conduisons, il y a beaucoup de techniques, et la méthode qui coordonne ces techniques est tout sauf agnostique. Ce qui nous motive tous est l’extraction de la physique de nos mesures.
Le problème qui unit l’ingénieur et le chercheur réside dans cette image du monde que nous transmettent nos mesures car elle n’est jamais exempte de biais qui peuvent jusqu’à empêcher parfois l’extraction de la moindre information physique. Les moyens permettant de diminuer les biais des observations diffèrent en général en fonction des objectifs :

  • lorsque l’on connaît précisément les phénomènes physiques en jeu, il est possible d’introduire la physique en tant que régularisation de l’identification de la transformation ; cette technique est très utilisée en mécanique des solides ;
  • lorsque l’on cherche à comprendre les mécanismes régissant les phénomènes, il est plus délicat de régulariser en se servant de la physique régissant les phénomènes, car celle‐ci n’est pas toujours connue ; de plus, la notion de matériau peut laisser place à une structure mécanique en fonction de l’échelle avec laquelle on regarde (en mécanique, la notion de structure est relative à des matériaux de type hétérogène, ou d’un ensemble de matériaux assemblés).

Le Laboratoire de mécanique de Lille, comme beaucoup d’autres, s’intéresse particulièrement aux mécanismes, que ce soit en tribologie, turbulence de paroi ou en mécanique des matériaux, le fil rouge de ce qui est présenté ici est cette volonté de se plier aux exigences des mécanismes.

Platon a défini dans le livre VII de La République, une distinction entre le monde du sensible et le monde des idées : le monde du sensible ne serait que la projection du monde des idées, qui nous apparaîtrait déformé par notre propre nature ; ce monde du sensible est par essence ce que nous pouvons imager, sentir.
Je considère que le travail premier d’un scientifique qui met en œuvre de l’expérimental est de savoir analyser les images des phénomènes observés pour ensuite essayer d’accéder au monde des idées. Ce point se heurte toujours au problème de la fidélité de ce que l’on observe par rapport à l’idée qui gouverne le phénomène. Le lecteur pourra se référer à l’allégorie de la caverne afin de bien illustrer ce propos.
Tout l’enjeu du scientifique est donc d’être capable d’avoir des images du monde sensible en connaissant les inexactitudes qui y sont attachées afin de pouvoir se faire une idée de ce qui en est la cause la plus probable.
Le monde de l’imagerie numérique nous donne accès à l’observation de l’hétérogénéité des phénomènes et nous permet de remettre en cause ce qui n’était observable jadis que par la variation durant le temps d’une ou plusieurs quantités ponctuelles. En tant que mécanicien, je me limiterai à l’étude des conséquences de la mise en mouvement de la matière, nous verrons dans une prochaine dépêche que, malgré cette limitation, le champ d’application de l’utilisation de l’imagerie avec le miroir de la mécanique permet d’étudier une grande diversité de sujets et d’accéder à une meilleure compréhension des phénomènes en jeu.

Nous traiterons ici de données qui ne sont pas toutes directement associables à des lois a priori, de part la complexité des systèmes observés. Les outils que nous développons s’inscrivent donc dans une logique kantisée, mais elle se heurte malheureusement à également l’opposition des conatus des scientifiques.
Pour faire simple, les méthodes numériques généralistes présentées ici peuvent être utilisées pour appuyer les idées de scientifiques plutôt d’obédience rationaliste mais peuvent tout aussi bien être utilisées par des expérimentateurs plus portés sur l’empirisme. Cette unité de méthode n’empêche pas que les choix techniques (bien qu’appartenant tous au cadre présenté dans cet article) conduiraient à une épistémologie radicalement différente des analyses scientifiques menées par ces deux groupes sur un même objet d’étude.

Il convient donc de mettre en place une dialectique expérimentale propre à chaque problème traité.

Comment approcher ou établir une loi ou idée le plus justement possible ?

Au‐delà des problèmes de biais associés à toute mesure physique, se pose la question du but même de la mesure. Il y a, d’une part, les rationalistes, qui considèrent que toute idée doit être cause ou conséquence d’un enchaînement d’idées rationnelles ; ils chercheront donc à se servir de l’expérimental pour s’approcher de leur loi au travers d’une paramétrisation ad hoc ; et, d’autre part, les empiristes, qui considèrent que toute connaissance prend source dans une expérience sensible, qui chercheront eux à construire des lois au travers de ces observations. À ce titre, je trouve intéressante la citation de Pierre Duhem : « Une loi de physique n’est, à proprement parler, ni vraie ni fausse, mais approchée. » (Traité d’énergétique, 1911), car elle nous permet de discriminer ces courants de pensée expérimentale.
L’analyse d’image en mécanique ne fait pas exception à cette dualité : il nous arrivera de mettre beaucoup d’_a priori
dans la mesure elle‐même, la contraignant ainsi à se plier à notre vue, mais il sera aussi utile de trouver comment ne pas mettre du tout de physique dans la mesure, pour pouvoir construire nos idées. Ainsi, je présenterai dans cet article des approches que nous avons mises au point permettant d’utiliser la mécanique pour extraire au mieux une physique inaccessible sans cette dernière, telle que la corrélation d’images numériques (Digital Image Correlation, DIC) et la Particule Image Velocimetry (PIV), mais également des méthodes innovantes permettant la réjection de solutions aberrantes, telles que la régularisation par filtre médian dans le processus de DIC.

Recalage d’images et flot optique

Que ce soit en analyse d’image ou en mécanique, les techniques se réclamant du flot optique (ou flux optique) sont légion. Le flot optique recouvre énormément de techniques différentes. Beaucoup d’auteurs ont déjà fait des états des lieux de ces techniques (Galvin et al., 1998_, Corpetti et al., 2006_, Baker & Matthews, 2004_), cet article ne sera donc pas une nouvelle revue des techniques déterminant le flot optique. Il n’est cependant pas possible de comprendre les idées qui ont été développées sans situer le flot optique dans ce qui se fait dans la communauté.
Il me semble cependant plus intéressant d’aborder la présentation de ces démarches au travers du concept de recalage d’images qui est un canevas plus général que celui du flot optique.
Il existe plusieurs méthodes utilisées pour recaler des images, certaines se basent sur l’intensité lumineuse de l’image, d’autres cherchent à détecter l’écart informations de type géométrique (c.‐à‐d. la distance), d’autres encore se situent à mi‐chemin et sont appelées « iconiques ». Une classification a été donnée par Cachier et al., 2003, elle sera utilisée.

Principes théoriques du recalage d’images, PIV et corrélation d’images numériques

Le recalage d’images basé sur l’intensité des niveaux de gris est une technique permettant d’apparier des images entre elles au travers d’une transformation. Cette technique est indispensable dans deux domaines : la création de vues panoramiques à partir de plusieurs prises d’images et l’analyse d’images médicales. Nous allons aborder le second point, plus complexe, car il entraîne une description plus générale. Le recalage d’images médicales a pour but de pouvoir comparer deux images de modalités différentes entre elles ou de comparer plusieurs patients entre eux.

La DIC (Digital Image Correlation) consiste à placer une caméra devant un essai (traction, compression, flexion, multi‐axial) où l’on ajoute une texture aléatoire.
La PIV (Particule Image Velocimetry, vélocimétrie par image de particules) consiste à rajouter des particules dans un écoulement (soufflerie, cavité entraînée…) et à les éclairer avec une nappe laser afin de les filmer.
Ces deux techniques ont pour objectif exclusif d’identifier les transformations.

Comme ces méthodes sont basées sur la comparaison de deux images, il convient tout d’abord de savoir comment effectuer une telle comparaison.

Prenons par exemple le cas d’une image d’un patient prise par un IRM et par un scanner. L’IRM repose sur le contraste de résonance magnétique nucléaire alors que le scanner repose sur un contraste d’atténuation des matériaux que traversent les rayons X. Il y a donc peu de chances que les constituants du corps du patient réagissent de la même manière avec ces deux modalités différentes. Ces images ne peuvent donc pas être recalées via le flot optique, il faudra pour garder l’heuristique du recalage trouver un moyen de comparer ces images, ce qui sera introduit par la suite par la notion de métriques.
Là encore, beaucoup a déjà été écrit à ce sujet mais un point d’entrée intéressant se trouve dans l’article de Wachinger et Navab (2010) qui illustre de manière très pédagogique pourquoi une métrique telle que celle utilisée dans le milieu de la mécanique ne peut pas fonctionner dans ce cas précis. Nous allons nous placer dans un canevas développé par Insight Tool Kit (ITK) introduit dans (Yoo et al., 2002), détaillé dans le manuel d’utilisation (Ibanez et al., 2005), et utilisé dans le logiciel de recalage ElastiX (Klein et al., 2010). Cette démarche est présentée en figure 1 :

Figure 1Figure 1 : Principe du recalage d’image ITK et Elastix

Par convention, on appellera l’image fixe, l’image de référence, et l’image mobile, l’image que l’on souhaite recaler par rapport à cette référence. Le principe est donc le suivant :
L’image mobile est recalée par rapport à l’image fixe au travers d’une transformation. Pour identifier au mieux cette transformation, plusieurs étapes sont nécessaires. Il faut tout d’abord définir une métrique qui permettra de comparer les deux images. Cette métrique a comme variables d’entrée un type de transformation (Transformation) (rigide, affine similitude, B‐Spline…), l’image fixe et l’image mobile. Il faut ensuite choisir une stratégie d’optimisation (Optimiseur) dans l’algorithme du gradient (Cauchy, 1847), comme la Méthode du point col, Newton, gradient conjugué… Il est également possible de ne pas évaluer la fonctionnelle (métrique) sur toute l’image mais sur un certain nombre de points (échantillonnage) aléatoire, sur une grille ou total. Enfin, il faut appliquer la transformation à l’image mobile et ce à chaque itération, ce qui nécessite un interpolateur.
Les techniques de type Besnard et al., 2006, Réthoré et al., 2007a, reposent également sur ce schéma mais ne sont qu’un cas particulier, ce qui sera détaillé par la suite.
Le problème à résoudre se met en équation comme présenté :

Dans cette équation : \mathcal{C} est la métrique, {T}_{\mu} la transformation paramétrée par \mu, \mathcal{I_F} l’image fixe et \mathcal{I_M} l’image mobile.

Rappelons que :

\left \{\Phi(x)\right \} est le vecteur des transformations comme elles seront détaillées par la suite.

Ces méthodes étant basées sur un algorithme d’optimisation de descente suivant le gradient, elles peuvent converger vers un résultat faux car le résultat est un minimum local de la fonction objectif (métrique). Pour régler ce problème, une technique de résolution pyramidale est employée, comme l’illustre la figure 2 :

PyramideFigure 2 : Principe de l’approche pyramidale, une transformation sur une image très moyennée est trouvée, puis interpolée pour être appliquée en initialisation de l’échelle plus résolue suivante, et ce, jusqu’à la dernière échelle.

Nous allons détailler quelques points importants du cadre utilisé ici.

Métrique C

Ce qui est appelé ici « métrique », ne correspond pas à la définition mathématique, car il ne satisfait pas les critères imposés pour la définition d’une distance. Il serait plus exacte de l’appeler fonction objectif, mais la communauté du recalage d’image a choisi ce terme, aussi il sera conservé ici. La littérature en compte tellement qu’il est impossible d’en faire une liste exhaustive. Celles‐ci se classent dans deux catégories distinctes, les méthodes statistiques basées sur des probabilités jointes et les méthodes directes. Dans les méthodes directes il faut être en mesure de postuler d’un lien entre les évolutions des niveaux de gris de l’image fixe et ceux de l’image mobile. Par exemple, si le flot optique est conservé entre deux images alors la métrique naturelle est l’erreur quadratique moyenne. S’il y a une relation affine entre l’évolution des niveaux de gris des deux images, alors la métrique naturelle est l’inter‐corrélation normalisée. S’il n’y a pas de relation triviale entre les évolutions de niveaux de gris, alors il faut construire une relation propre à la paire d’images étudiée au travers d’une métrique de type statistique comme l’information mutuelle. Les paragraphes suivants présentent quelques‐unes des métriques utilisées en recalage d’images médicales et applicables dans tous les domaines.

Erreur quadratique moyenne (C=SSD)

L’erreur quadratique moyenne (ou Sum of Squared Differences, SSD aussi appelée mean squared error) correspond à un problème au moindre carré présenté par l’équation suivante :

\Omega_F est le domaine d’intégration, \left|  \Omega_F \right | le nombre de pi(vo)xels (cardinal) pour une transformation \mu donnée.

Inter‐corrélation normalisée (C=NCC)

L’inter‐corrélation normalisée (ou Normalized Cross Corrélation, NCC) est une technique utilisée dans beaucoup de codes de PIV ou DIC de type blockmatching (Limb et Murphy, 1975). Elle n’est pas sensible au changement de luminosité affine et permet donc de travailler sur des images dont l’éclairement de l’image change fortement entre deux positions. C’est notamment le cas de la PIV, où la nappe laser peut ne pas être parfaitement constante, il peut y avoir de forts reflets et des particules peuvent sortir du plan du laser. C’est également le cas en imagerie médicale où les organes peuvent sortir du plan d’observation. Il est donc impossible dans de telles conditions de respecter le flot optique. Il est de plus intéressant de rappeler qu’usuellement les chercheurs utilisant ces moyens pour la turbulence règlent le temps entre les deux lasers de manière à avoir 10 pixels de déplacement moyen, ce qui peut changer fortement l’illumination des particules. La métrique d’inter‐corrélation normalisée présentée par l’équation suivante est donc théoriquement un bon moyen de résoudre ces contraintes expérimentales :

avec

et

On pourrait également parler de l’information mutuelle, mais je pense que cela nuirait à la concision de l’article.

Échantillonnage

L’échantillonnage est défini comme le domaine (c.‐à.‐d. l’ensemble de pixels) sur lesquels la fonction objectif sera évaluée : toute l’image, un masque, une grille, ou tirée aléatoirement. Likar et Pernuš, 2001, ont démontré que le ré‐échantillonnage d’une fonction objectif, telle que l’information mutuelle, lors de l’application de la transformation pouvait conduire à un bruitage important de \mathcal{C}. Beaucoup de solutions existaient avant ce travail mais elles passaient par l’utilisation de bases d’interpolation d’ordre élevé ce qui entraînait un surcoût en temps de calcul. L’idée de Likar et Pernuš, 2001, a été de n’évaluer la fonctionnelle que sur un nombre réduit de points définis par l’utilisateur et tirés aléatoirement à chaque itération de l’algorithme de descente. Dans la pratique cette stratégie se comporte assez bien sur les deux métriques présentées précédemment mais pas avec tous les optimiseurs. Il est également souvent utilisé une sous‐intégration sur une grille ou une intégration totale. La stratégie à choisir dépend beaucoup de l’objectif que l’on se fixe : si le but est d’avoir une mesure la plus fine possible, alors la sous‐intégration pose problème ; en revanche, si l’objectif est d’obtenir une évaluation rapide des paramètres pour, par exemple, envisager une rétroaction à partir de la transformation mesurée (cas du pilotage d’un essai en déformation homogène), alors l’échantillonnage est certainement un point à ne pas négliger.

Transformation

Le modèle de transformation est un autre point clef de la méthode, pour une quantité donnée d’information, nous souhaitons extraire les paramètres d’une transformation. On peut énoncer plusieurs types de transformations :

Transformations globales

Transformation globaleFigure 3 : Illustration de déformation globale d’une image (rotation, translations, déformations homogènes)

Elle sont souvent linéaires, telles que le mouvement de corps rigide, la transformation affine (figure 3), la similitude, la transformation projective, qui sont souvent utilisées dans l’analyse d’images médicales. Il existe aussi des transformations globales utilisées dans la communauté de la mécanique des matériaux, telles que celles associées à la fissuration (Hamam et al., 2007, Hild et Roux, 2006) ou celles associées aux poutres (Hild et al., 2011).

Transformations élastiques (non rigides)

Transformation elastiquesFigure 4 : Illustration des transformations élastiques d’une image (utilisation d’une base quadrangle bilinéaire)

Il existe également une gamme de transformations appelée « élastique » dans la communauté de l’analyse d’images médicales. Elle comprend toutes les transformations qui assurent la partition de l’unité, comme les transformations basées sur les polynômes de Lagrange souvent appelées éléments finis (Besnard et al., 2006, figure 4), mais aussi les grilles de type B-Splines (Klein et al., 2010, Réthoré et al., 2010) et les méthodes étendues de type X-FEM (Réthoré et al., 2007b). Il est important de noter que ce terme « élastique » n’est pas accepté dans la communauté de la mécanique car une déformation élastique a une signification en mécanique : ainsi, il est possible d’avoir une déformation qui respecte la partition unité, mais qui n’est pas élastique (endommagement, plasticité, viscosité, etc.). Les mécaniciens du solide ne distinguent pas cette catégorie et rangent ces méthodes de transformation dans les méthodes globales (car il faut une résolution globale pour accéder à la transformation). On pourra utiliser de manière moins ambiguë la notion de « transformation non rigide ».

Transformations locales

Transformation localesFigure 5 : Illustration des transformations locales d’une image (seules les translations sont ici identifiées)

Ces transformations, présentées en figure 5 s’appliquent à une sous‐partie de l’image et ne cherchent pas à résoudre un problème global comme les deux autres. Elles n’ont pas la nécessité de respecter la partition de l’unité. La plus connue et ancienne (Barnea et Silverman, 1972) est la recherche des translations sur chaque imagette par inter‐corrélation. Il est à noter que la méthode la plus rapide est incontestablement la technique d’inter‐corrélation, qui n’est pas un problème itératif et peut être accélérée en utilisant la FFT. Une amélioration présentée par Sutton et al., 1983, est de calculer la translation ainsi que les déformations de l’imagette. Elle est utilisée dans la plupart des logiciels commerciaux de PIV et DIC. Les transformations locales diffèrent des autres transformations de part le fait que chaque zone d’intérêt (Zone of Interest, ZOI) est traitée indépendamment. Dans la communauté de l’analyse d’images, ce problème a été résolu très tôt par Lucas et Kanade, en 1981, qui détermine le flot optique sur des ZOI. Cette méthode est de loin la plus intéressante car elle permet d’utiliser les transformations définies précédemment, mais à l’échelle locale. C’est de plus la seule méthode de blockmatching qui trouve sa place dans le canevas présenté.

Interpolateur

L’interpolateur sert à évaluer la valeur de \mathcal{I_M}({T}_{\mu}(x_i)) et {\partial \mathcal{I_M}({T}_{\mu}(x_i))}/{\partial x} lors de l’évaluation de la métrique. En effet, les déplacements sont recherchés à une précision inférieure au pixel, il faut donc aller interpoler les valeurs au voisinage de la position définie. Beaucoup d’interpolateurs sont utilisés dans la littérature mais les B‐splines prédominent fortement. Une excellente revue de l’impact d’un interpolateur sur l’image est faite par Getreuer, 2011.

Figure 6Figure 6 : Impact des étapes successives d’interpolation d’une rotation en fonction de l’interpolateur choisi (Getreuer, 2011).

La figure 6 illustre le fait que si l’on cumule des successions de transformations (ici une rotation), alors le choix de l’interpolateur devient fondamental. Sur cette figure, seule une B‐spline d’ordre 11 parvient à ne pas modifier fortement la structure de l’image. Cela montre aussi qu’il est fondamental de ne pas cumuler les interpolations, aussi la stratégie utilisée ici est de toujours travailler avec \mathcal{I_M}({T}_{\mu}(x_i)) (la transformation totale est appliquée à l’image mobile initiale) et non avec des successions d’images interpolées de manière cumulative (comme présenté figure 6).

Optimiseur

Dérivation de C

Il convient de rappeler que les méthodes d’optimisation utilisées ici sont toutes des méthodes basées sur le calcul du gradient de \mathcal{C}. Le procédé est détaillé sur un exemple en calculant le gradient de l’erreur quadratique moyenne définie précédemment en utilisant la composition des dérivées.

Pour chaque métrique utilisée il faut calculer le gradient de la fonctionnelle.

Algorithme du gradient

Formalisation mathématique

Une fois le calcul du gradient de la fonction objectif effectué, il faut trouver une stratégie pour converger vers son minimum.

Le principe de la minimisation de la fonctionnelle est défini par l’équation suivante :

La taille du pas de l’algorithme de gradient est pilotée par la valeur de \alpha_k qui peut être constante ou remise à jour à chaque itération suivant différentes méthodes.

Soit :

Un point fondamental dans une descente suivant le gradient est de déterminer quelle direction de recherche d_k utiliser. Ce choix se fait au travers de \gamma_k ou itérativement par choix des d_k (gradient conjugué et stochastique), nous présenterons ici les principales méthodes utilisées en optimisation.

  • descente suivant le gradient :
  • Newton avec [H_k] la matrice hessienne, telle que H_{ij}(C) = {\partial^2 \mathcal{C}(\mu,\mathcal{I_F},\mathcal{I_M})}/({\partial \mu_i\partial \mu_j}) :
  • Gauß‐Newton :
  • Levenberg‐Marquardt :
  • quasi‐Newton (BFGS) :
  • gradient conjugué :
  • gradient stochastique :

Pour les transformations de type élastique, ces relations sont vraies mais doivent être appliquées au niveau de l’élément (FEM, X-FEM) ou à l’échelle du patch (B‐splines).

Discussion des méthodes d’optimisation

Il est possible de construire l’opérateur :

en se servant de la moyenne des gradients de l’image fixe et mobile ce qui permet d’accélérer la procédure d’optimisation. Lorsque cette moyenne est appliquée à la méthode de Gauß‐Newton la méthode s’appelle ESM pour Efficient Second order Minimization, elle a été inventée par Benhimane et Malis, 2004, qui prouvèrent qu’elle avait une excellente convergence pour la SSD, et validée sur d’autres métriques par Wachinger et Navab, 2009. Cette méthode a également été utilisée dans notre communauté par Réthoré et al., 2007a. Il est important de noter que certains auteurs choisissent de ne pas mettre à jour ce terme en le calculant à partir de l’image fixe, dans ce cas on identifie la transformation opposée. Cette technique a l’avantage de ne calculer qu’une seule fois la dérivée de l’image et la matrice Hessienne (pour \mathcal{C} =SSD), ce qui permet un gain de calcul important dans le cas de la méthode de (quasi ou Gauß) Newton. Ceci peut être très avantageux pour le gradient conjugué, s’il stocke les vecteurs construits au fur et à mesure jusqu’à ce que la base approchée soit stationnaire. Cette base pourrait être utilisée sur l’image suivante d’une séquence comme initialisation, ce qui crée un gain de temps notable et, à partir d’un certain nombre d’itérations, il ne rajoute plus de vecteurs à la base à chaque nouvelle itération du processus d’optimisation. Le principal problème des méthodes de type Newton ou quasi‐Newton est la sensibilité au choix des paramètres initiaux. Si les paramètres sont éloignés, alors il est probable que l’algorithme converge vers un minimum local. S’il est vrai que la technique pyramidale présentée ici permet d’estimer des paramètres à des échelles plus grossières évitant ainsi le piège local, il arrive que cette approche ne suffise pas. D’un autre côté, les méthodes de descente suivant le gradient permettent un écart plus grand par rapport aux paramètres initiaux mais convergent plus lentement. La méthode de Marquardt, 1963, permet de profiter des avantages des deux méthodes mais le choix du \lambda et son évolution au cours du temps peut entraîner une convergence faible au niveau de l’optimum.

Une des méthodes les plus performantes pour résoudre ce système est la méthode quasi‐Newton de type Broyden‐Fletcher‐Goldfarb‐Shanno où les vecteurs V1_{k−1} et V2_{k−1} sont des vecteurs de rang 1 de base différente. Les détails de la création des vecteurs V1_{k−1} et V2_{k−1} sont disponibles sur Wikipédia. Un point important est qu’il est possible de déterminer explicitement l’inverse de l’approximation de la Hessienne (Sherman et Morrison, 1950). De ce fait, cette méthode peut également être utilisée efficacement conjointement avec l’ESM.
Une attention toute particulière doit être apportée à la méthode de type gradient stochastique. Cette technique d’optimisation a été développée pour résoudre des problèmes de grande taille dans le domaine de l’apprentissage automatique. Lorsque les systèmes deviennent très gros, il est délicat de calculer l’inverse de la hessienne du problème, voire parfois impossible. Il reste alors comme possibilité de l’approximer (Gauß‐Newton, quasi‐Newton). Le coût numérique étant alors assez élevé, il est alors possible de faire une simple descente suivant le gradient, mais les propriétés de convergence sont alors beaucoup moins bonnes qu’avec une méthode de Newton, car la direction de recherche n’étant pas bonne, il faut beaucoup plus d’itérations. L’idée du gradient stochastique est de ne pas évaluer la fonctionnelle sur toutes les données mais seulement sur un échantillonnage tiré aléatoirement. Bottou et Bousquet, en 2008, ont montré que dans ce cas, l’utilisation de la Hessienne n’améliorait pas la vitesse de convergence et que cette méthode converge plus vite (en termes de temps et non d’itérations) pour les problèmes de grande taille que les méthodes qui approximent la matrice Hessienne. Cette méthode permet donc de s’affranchir du calcul de la matrice hessienne, assure une convergence vers un minimum global dans le cas d’une fonctionnelle convexe et vers un minimum local pour une fonctionnelle non convexe. Klein et al., en 2005, présente une accélération d’un facteur 500 par rapport à un échantillonnage complet. Il a de plus été montré dans la sous‐section échantillonnage que le sous‐échantillonnage permet de lisser la fonction objectif. Nous rappelons que l’obtention d’une fonctionnelle convexe est assurée par le processus pyramidal.

Stratégie de régularisation

La stratégie la plus utilisée en ce qui concerne la régularisation est celle mise en œuvre par Tikhonov et Glasko, en 1965, et présentée dans l’équation suivante ; elle consiste à favoriser les solutions dont les normes des inconnues sont petites :

Dans notre cas :
A = \alpha_k en fonction de la stratégie d’optimisation choisie et b est égal à la dérivée de la dérivée de \mathcal{C}. La matrice \Gamma peut quand à elle être l’identité ou n’importe quel opérateur linéaire tel que, par exemple, la matrice de raideur d’un système à élément finis ou alors une matrice représentant un filtre passe‐bas peut être utilisée pour éliminer les variations rapides qui sont non différentiables du bruit de mesure.

Après de longues réflexions, nous avons décidé de ne pas implémenter ce canevas pour la régularisation dans nos outils. Bien qu’apparemment séduisantes, les régularisations de Tikonov ne permettent que d’utiliser des normes L_2, car les méthodes d’optimisations implémentées nécessitent une fonctionnelle convexe. Cela entraîne, par exemple, que les régularisations de type L_1 (qui conservent les discontinuités) ne nous sont pas accessibles. Il est possible de formuler le problème différemment, Cachier et al, 2003 ont montré qu’il est possible de formuler le problème de recalage d’image iconique de type démon en tant que problème de flot optique au travers d’une minimisation alternée.
Ainsi, ils formulent le problème de recalage de la manière suivante :

Les termes \mu et \lambda sont les transformations, c’est‐à‐dire des translations pixel à pixel, il n’y a donc pas de base de transformation ici, la cohérence spatiale étant assurée par les énergies de pénalisation pondérées par les coefficients \sigma et \lambda ; R ( \mu ) représente une énergie de pénalisation telle que le laplacien du déplacement en mécanique des fluides ou div ( \sigma_m ) en mécanique du solide.

La méthode itérative de résolution est la suivante :

  • d’abord résoudre l’équation suivante en utilisant un algorithme de descente au gradient en version ESM :
  • puis, résoudre l’équation suivante :

Cette équation pourrait se résoudre également via un algorithme de descente au gradient, mais il existe un moyen de la résoudre par convolution comme démontré en annexe de Cachier et al., 2003_.

Nous nous sommes inspirés de cette idée, mis à part le fait que nous ayons généralisé cette méthode de régularisation à n’importe quelle méthode de filtrage et non plus à une simple convolution. Après différents tests, nous avons constaté que le filtre médian était l’outil de régularisation le plus efficace dans la plupart des cas.


La prochaine dépêche présentera les logiciels utilisés dans les différentes communautés scientifiques.

Aller plus loin

  • # optim non convexe

    Posté par  (site web personnel, Mastodon) . Évalué à 9.

    Il y a aussi des travaux sur l'emploi de méthodes d'optimisation « globales » pour des gros pb non–convexes. C'est de l'approximation mais ça va vite. On peut voir ça comme une descente de gradient stochastique robuste. Utile aussi pour gérer des modèles de déformation. Mots-clefs: randomized search, black-box optimization, metaheuristics.

    • [^] # Re: optim non convexe

      Posté par  . Évalué à 6.

      Ça à l'air vraiment intéressant les randomized search ! Je ne suis pas vraiment un matheux et du coup j'ai du mal à comprendre l'ensemble. Si j'ai bien compris tu te places dans une sphère de dimension n puis tu cherches une solution sur l'hypersphère de dimension n-1. Si je sais comment me placer sur un hyperplan à partir de la fonction objectif (calcul des dérivées), mais du coup lorsque je mets à jour mon algorithme, je quitte le groupe de départ car je me ballade sur le plan tangent et non sur la sphère de départ. L'idée serait ici de passer non pas sur l'hyperplan mais sur l'hypersphère, je n'ai pas bien compris comment on se place sur une hypersphère. Visiblement c'est via la notion de "normal deviate" mais j'ai encore un peu de mal à bien saisir le principe. C'est ce qui remplace l'opération de gradient, mais sa construction me semble encore obscure.

      Ça à l'air très intéressant mais comment cela se comporte il avec l'augmentation de n ?

      • [^] # Re: optim non convexe

        Posté par  . Évalué à 8. Dernière modification le 04 juillet 2017 à 23:29.

        L'idée serait ici de passer non pas sur l'hyperplan mais sur l'hypersphère, je n'ai pas bien compris comment on se place sur une hypersphère. Visiblement c'est via la notion de "normal deviate" mais j'ai encore un peu de mal à bien saisir le principe. C'est ce qui remplace l'opération de gradient, mais sa construction me semble encore obscure.

        C'est comme se déplacer sur une courbe de même latitude sur la surface terrestre (un cercle sur une sphère). Cela doit convoquer les opérateurs de la géométrie différentielle qui permet de faire du calcul différentiel dans des espaces non euclidiens (comme une sphère) à la manière de la relativité générale de tonton Albert.

        Sinon jolie dépêche ! :-) Je n'en ai pas encore approfondi la lecture, mais comme tu abordes tant Planton et sa caverne, que Kant et son épistémologie dans ton introduction, je me sens dans l'obligation de répondre. Cependant, je différerai mon commentaire sur le sujet (j'ai beaucoup à dire dessus) et espère pouvoir le faire d'ici la fin de semaine. En tout cas, je suis heureux de voir qu'un laboratoire de physique s'y intéresse (j'avais quelque peu perdu espoir sur la question ;-).

        Sapere aude ! Aie le courage de te servir de ton propre entendement. Voilà la devise des Lumières.

  • # Super dépêche

    Posté par  . Évalué à 9.

    Tout est bien expliqué je trouve, et ça donne rapidement une idée générale du champ d'étude.

    En ce moment j'étudie les réseaux de neurones, et ça me fait me rendre compte que toutes ces méthodes d'optimisation sont fondamentales et très communément utilisées.

    Aussi petite faute de frappe dans le paragraphe sur le mean squarred error : cardianal / cardinal.

    Merci pour ce gros travail !

    • [^] # Re: Super dépêche

      Posté par  . Évalué à 7.

      Oui, article de haut niveau. De nature encyclopédique donc digne de Wikipédia.

    • [^] # Re: Super dépêche

      Posté par  (site web personnel) . Évalué à 3.

      Corrigé, merci.

    • [^] # Re: Super dépêche

      Posté par  . Évalué à 10.

      Ça m'a fait plaisir de le faire, c'est à la frontière entre plusieurs communauté, c'est donc difficilement publiable. De plus cela peut toujours être délicat de dire que ce que quelqu'un de connu a trouvé dans un domaine existe depuis longtemps ailleurs, je préfère donc m'éloigner un peu de mon domaine pour publier ce genre de choses.

  • # mécanicien :)

    Posté par  . Évalué à 10.

    Bravo pour ce journal, très intéressant.

    Si je peux me permettre, on sent vraiment que l'auteur travaille dans un environnement de mécaniciens, notamment par les exemples pris et les développements mathématiques présentés. N'y voyez pas de jugement :)

    Si l'on veut aborder d'un autre point de vue le recalage d'image, il faut également mentionner toutes les techniques basées sur des descripteurs (parfois appelés motifs):
    - détection de points d'intérêts (FAST, SIFT, SURF…), comme par exemple des coins.
    - calcul de descripteurs pour ces points d'intérêts (on agrège les données dans le voisinage de ces points).
    - mise en correspondance des descripteurs, élimination des "outliers" (points qui ne sont pas reconnus d'un côté ou de l'autre)
    - calcul d'une transformation

    Si on prend hugin qui est pour moi la référence en terme de logiciel libre de recalage d'images, vous avez la possibilité d'utiliser des outils comme autopano-sift-c pour le calcul des points d'intérêt et la correspondance des descripteurs.

    Pour ce qui est des méthodes basées sur la corrélation et la transformée de Fourier, vous avez dans openCV les outils de "templateMatching" qui le font (dans la doc, vous verrez un certain nombre de métriques). Le principal inconvénient de ces techniques est qu'elles ne sont pas robustes au changement d'échelle ni à la rotation. Les méthodes basées sur les descripteurs permettent de s'affranchir de ces problèmes, moyennant un temps de calcul plus important.

    • [^] # Re: mécanicien :)

      Posté par  . Évalué à 8.

      Tu as tout à fait raison. J'ai restreint mon champs d'investigation à ce qui pouvait être commun entre les usages que je connais dans les différentes sciences.

      Dans la présentation des logiciels, je parlerais rapidement de DeepFlow, implémenté dans OpenCV, qui utilise une métrique basée sur la notion de distance au travers de l'agorithme DeepMatchning en calculant ds convolution sur des imagettes de 4x4 puis les agrège de manière à pouvoir traiter les rotations et éliminer les faux positifs. Cela me semble être un bon point d'entrée.
      Tu peux aller voir la dépêche initiale dans l'espace de rédaction. Les modérateurs m'ont proposés de la scinder en 3, de ce fait je n'ai présenté que l'aspect théorique que je maîtrisais suffisamment pour le présenter et en débattre.

      N'hésites pas à contribuer sur la partie logicielle !

  • # Les ToolKits et logiciels desktop libres du traitement d'images

    Posté par  . Évalué à 9.

    Super dépêche.
    Je rejoins l'idée que l'article étant complet, et auto-contenu, il mériterait sa place dans Wikipédia.

    Juste pour préciser, parce que c'est toujours mieux en l'écrivant clairement, et que ça fait de la pub pour une boite engagée dans le libre, et qui se trouve être la mienne.

    Kitware, dont une partie des équipes est à Lyon, développe et maintient ITK et VTK (et CMake, ParaView, 3D Slicer, KWIVER … et plein d'autres outils libres (BSD/Apache), mais c'est hors sujet!).
    Les toolkits (ITK, VTK), bas-niveaux et écrits en C++, sont utiles aux programmeurs pour faire des outils comme ElastiX, ou scripter en Python, C++, Java, Tk, … grâce aux bindings fournis.

    Les applications Desktop, libres aussi, fournissent des interfaces clés en main, avec des interfaces utilisateur basées Qt. L'utilisateur peut programmer sa chaîne de traitement à la souris, et aussi scripter en Python.

    • ParaView est la "GUI" pour la visualisation et l’analyse de données, construite autour de VTK (Visualization Toolkit). Spécialisé autour du parallélisme de données, pour la gestion de données très massives, avec gestion client-serveur, HPC, ferme de serveurs,…
    • 3D Slicer est la "GUI" pour la visualisation et l’analyse d'images 2D/3D, IRM, CT, ultrason, construite autour de ITK.
    • RTK pour la reconstruction rapide de CT à faisceau conique circulaire. RTK se base sur ITK.
    • MAP-Tk GUI pour la reconstruction 3D a partir d'images ou vidéos 2D, construite autour de KWIVER.

    Bastien Jacquet
    Développeur à Kitware

    • [^] # Re: Les ToolKits et logiciels desktop libres du traitement d'images

      Posté par  . Évalué à 7.

      Bonjour Bastien,

      Merci pour le commentaire sur la dépêche et surtout merci pour les briques ITK/VTK !

      Je trouve dommage que les news autour de vos logiciels et bibliothèques ne soient pas relayées ici. Je pense qu'avec la montée en maturité de SimpleITK vos outils sont aujourd'hui accessible aux scientifiques non-informaticiens (i.e : non gourous C++).

      Paraview est un outil exemplaire qui permet de faire des post-traitements fantastiques, scriptable en python.

      VTK est une brique fondamentale du traitement de surfaces triangulaires tridimensionnelles, un peu compliquée à prendre en main mais avec des outils tels que paraview ça simplifie bien les choses !

      Ça serait bien sympa de vous lire plus souvent ici.

Suivre le flux des commentaires

Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.