Géométrie sous-riemannienne en dimension infinie et applications à l’analyse mathématique des formes
Analyse mathématique des formes
L’analyse de formes et de déformations d’images a reçu une attention croissante au cours des dernières années. Il s’agit de donner un cadre géométrique au problème afin de permettre, par exemple, l’analyse statistique d’images et d’objets déformables. Ce point de vue a été utilisé avec succès dans de nombreux domaines, et plus particulièrement en anatomie computationnelle. Cette dernière consiste à comparer des images d’organes (généralement des IRM de cerveaux) de différents patients, dans le but de prévoir à l’avance l’apparition de certaines maladies dégénératives [DGM98, GM98, MTY02, MTY06].
La source du problème est que deux formes (par exemple des courbes planes), peuvent représenter des objets similaires pour un observateur humain, mais ˆetre très différentes au niveau de la liste de points qu’un ordinateur utilisera pour les représenter. Il faut donc, pour les comparer, trouver des méthodes qui prennent en compte leur géométrie. Plusieurs fa¸cons de résoudre ce problème en utilisant des concepts mathématiques rigoureux ont été proposées avec succès. Pour cela, il faut tout d’abord définir un espace de formes S sur lequel on puisse travailler. Il faut ensuite trouver une manière de comparer deux formes q0 et q1 dans S. Une des premières solutions a été proposée par Kendall [Ken84]. Il s’agit de définir l’espace des formes discrètes à n points dans R d comme l’ensemble des n-uplets de R d , pour n ∈ N, quotienté par l’action des rotations et des translations. On peut alors estimer la différence entre deux formes q0 = (x1, . . . , xn) et q1 = (y1, . . . , yn) en calculant le carré de leur distance à isométrie près : inf (Xn i=1 |T(xi) − yi | 2 | T : R d → R d isométrie) . Plus généralement, on peut considérer une attache aux données g, c’est-à-dire une fonction à valeurs réelles sur (R d ) n qui compare la différence entre une forme donnée q et la cible q1 : inf n g(T(q)) | T : R d → R d isométrieo.
On prend généralement g(q) = d(q, q1) 2 o`u d est une distance qui peut ˆetre différente de la distance euclidienne. Ce cadre a été le point de départ de plusieurs méthodes pour l’analyse des déformations [DM98, KBCL99, Sma96]. Toutefois, elles sont très restrictives. On ne peut pas, par exemple, étudier des familles de courbes : il faut d’abord trouver un moyen efficace de les discrétiser. Un problème plus important est que ce modèle ne prend pas en compte toutes les propriétés géométriques des formes discrétisées. Par exemple, deux formes q0 et q1 découlant l’une de l’autre par une déformation lisse (par exemple un difféomorphisme de R d ) devraient ˆetre ”plus semblables” que si q1 est un déplacement plus aléatoire de q2. Transformation par difféomorphismes. Une méthode efficace pour tenir compte de ces considérations est d’étudier des transformations induites par des difféomorphismes.
En effet, les formes étudiées sont généralement des objets plongés dans une variété M (le plus souvent, M = R d ). Le groupe des difféomorphismes de M agit alors naturellement sur l’espace des formes par (ϕ, q) 7→ ϕ · q. Par exemple, si q est une courbe paramétrée de M, un difféomorphisme ϕ agit sur q par composition à gauche, tansformant q en une nouvelle courbe ϕ · q = ϕ ◦ q. On estime alors la différence entre deux formes par inf {F(ϕ) + g(ϕ · q0) | ϕ difféomorphisme de M} , o`u F(ϕ) est un coˆut défini sur le groupe des difféomorphismes, et g une attache aux données estimant une distance entre la forme transformée ϕ · q0 et une forme cible q1. Flots de champs de vecteurs. Vient alors le problème de définir un coˆut F avec un bon comportement. Pour cela, on se restreint à des difféomorphismes engendrés par des flots ∂tϕ X(t, x) = X(t, ϕX(t, x)), (t, x) ∈ [0, 1] × M, de champs de vecteurs X dépendant du temps sur M. Un tel flot ϕ X est une déformation de l’espace ambiant M, et induit une déformation t 7→ q(t) = ϕ(t) · q0 d’une forme initiale q0. Un champ de vecteurs X induit donc une déformation infinitésimale X · q d’une forme q, qui s’intègre en une déformation de q. Notons Γ(TM) l’espace des champs de vecteurs de M.
L’ensemble des déformations infinitésimales autorisées est alors défini comme un ensemble de champs de vecteurs X appartenant à un certain sous-espace pré-hilbertien (V,h·, ·i) de Γ(TM). Les éléments de V sont en général des champs de vecteurs lisses. La différence entre deux formes q0 et q1 est alors mesurée par la minimisation, sur l’ensemble des champs de vecteurs t 7→ X(t) ∈ V , d’une fonctionnelle de la forme J(ϕ X) = J(X) 1 = 1 2 Z 1 0 hX(t), X(t)i dt + g(ϕ X(1) · q0) = A(X) + g(ϕ X(1) · q0). (1) Le premier terme mesure l’action A(X) (ou énergie) de la courbe de difféomorphismes définies par le flot ϕ X de X. Celle-ci induit une déformation de l’état initial q0. Le second terme est une attache aux données qui mesure la différence entre l’état déformé final ϕ(1)·q0 et la forme cible q1. Pour le choix de l’espace (V,h·, ·i), on rencontre deux tendances dans la littérature existante. Méthode 1 : métriques riemanniennes faibles. On prend pour V l’ensemble des champs de vecteurs de classe C∞ à support compact.
On choisit ensuite une norme préhilbertienne (généralement une norme Sobolev) pour h·, ·i. C’est une approche utilisée par exemple dans [BHM11, BHM12, MM07]. 1. Remarquons qu’il y a correspondance bijective entre les champs de vecteurs X et les déformations ϕ X. Un tel point de vue engendre naturellement une métrique riemannienne dite ”faible” (voir [Cla09]), invariante à droite sur le groupe de Lie (modelé sur un espace de Fréchet) des difféomorphismes de classe C∞ à support compact, de distance associée dR. L’action d’une déformation t ∈ [0, 1] 7→ ϕ X(t) engendrée par un champ de vecteurs t 7→ X(t) ∈ V est alors donnée par A(ϕ X) = 1 2 Z 1 0 hX(t), X(t)i dt = 1 2 Z 1 0 D ϕ˙ X(t) ◦ ϕ X(t) −1 , ϕ˙ X(t) ◦ ϕ X(t) −1 E dt. On peut également définir une longueur pour cette courbe, et une distance dR comme sur les variétés riemanniennes de dimension finie. Minimiser la fonctionnelle J de (1) revient alors à minimiser 1 2 dR(IdM, ϕ) 2 + g(ϕ · q0) parmi les difféomorphismes de classe C∞ à support compact. Méthode 2 : espaces de Hilbert à noyau reproduisant. On choisit un espace de Hilbert (V,h·, ·i) engendré par un noyau reproduisant K(x, y) : T ∗ y M → TxM [Tro95]. Un tel noyau permet d’effectuer très rapidement certains calculs nécessaires à la minimisation de J. Ce point de vue permet de travailler sur le groupe Ds (M) des difféomorphismes de classe de Sobolev Hs , qui est une variété de Hilbert et un groupe topologique, plutôt qu’une variété de Fréchet comme dans le cas précédent. Son espace tangent en l’identité est l’espace des champs de vecteurs de classe de Sobolev Hs.
Ce cadre, qui est celui adopté dans cette thèse, offre de nombreux avantages d’un point de vue numérique. En plus de préserver de nombreuses propriétés géométriques des formes étudiées, il a permis le développement d’une famille d’algorithmes efficaces de comparaison de formes : les méthodes LDDMM (Large Deformation Diffeomorphic Metric Matching) [BMTY05, DGM98, GTY, GTY06, JM00, MTY02, MTY06, Tro95]. Toutefois, on ne prend g ́en ́eralement pas V ́égal `à tout l’espace des champs de vecteurs de classe de Sobolev Hs , on prend un espace bien plus petit, souvent constitu ́e de champs de vecteurs analytiques. En conséquence, la g ́eom ́etrie invariante à droite engendre ́ e par V n’est pas riemannienne mais sous-riemannienne. Les auteurs des méthodes LDDMM n’étaient pas familiers avec la géométrie sous-riemannienne, et ont traité le sujet du point riemannien
Géométrie sous-riemannienne en dimension finie et infinie
On donne tout d’abord une définition générale de structure sous-riemannienne sur une variété M, introduite par Agrachev, Boscain et al. dans [ABC+10], permettant de prendre en compte des distributions horizontales à rang variable. Il s’agit d’un triplet (H, ξ, g), o`u (H, g) est un fibré vectoriel riemannien sur M et ξ : H → TM est un morphisme de fibrés vectoriels. On dit que (H, ξ) est un fibré tangent relatif sur M. Les courbes horizontales t 7→ q(t) sont les courbes telles qu’il existe un contrôle t 7→ u(t) ∈ Hq(t) tel que ˙q(t) = ξq(t)u(t). Un tel couple t 7→ (q(t), u(t)) est appelé un système horizontal. On définit ensuite la longueur et l’action d’un tel système respectivement par L(q, u) = Z q gq(t) (u(t), u(t))dt, A(q, u) = 1 2 Z gq(t) (u(t), u(t))dt, et on en déduit une distance sous-riemannienne, ainsi que la notion de géodésique sous riemannienne, comme dans le cas habituel. En dimension finie, on sait qu’il existe deux types de géodésiques, données par le principe du maximum de Pontryagin. On définit les hamiltoniens H1 , H0 : T ∗M ⊕ H → R par Hλ (q, p, u) = pξqu − λ 2 gq(u, u). (2) TABLE DES MATIERES ` 9 En (q, p, u) ∈ T ∗M ⊕ H fixé, si ∂uHλ (q, p, u) = 0, ce hamiltonien possèdent un gradient symplectique partiel ∇ωHλ (q, p, u) en (q, p) pour la forme symplectique canonique ω de T ∗M. Théorème 0.1. (Pontryagin [PBGM64]) Soit q : [0, 1] → M une géodésique sous riemannienne issue d’un système horizontal (q, u). Alors il existe (p, λ) 6= (0, 0), o`u p : [0, 1] → T ∗M est de classe de Sobolev H1 , et λ =0 ou 1 tels que ( 0 = ∂uHλ (q(t), p(t), u(t)), ( ˙q(t), p˙(t)) = ∇ωHλ (q(t), p(t), u(t)). (3) Si λ = 0 la trajectoire q est singulière (ou anormale). Elle est dite géodésique normale si λ = 1.
On peut de plus obtenir toutes les géodésiques normales en intégrant le flot du gradient symplectique ∇ωh du hamiltonien normal h : T ∗M → R défini par h(q, p) = max u∈Hq pξqu − 1 2 gq(u, u) . (4) Grˆace à ce dernier, on peut montrer ([Mon02, AS04]) que réciproquement, si (q(·), p(·)) satisfait (3), alors, quand λ = 1, q(·) est bel et bien une géodésique. Lorsque λ = 0, q(·) est un point singulier de l’application point final E qui à un système horizontal (q(·), u(·)) associe son point final q(1) = (q(·), u(·)). Cette application est étudiée plus en détails en section 1.4.2 Passons au cas de la dimension infinie. La géométrie riemannienne en dimension infinie est un sujet connu et étudié en analyse sur les variétés [KM97, Omo74, Cla09]. Elle se scinde en deux cas : les géométries fortes, dont les topologies induites par la métrique riemannienne et celles induites par la structure de variété co¨ıncident, et qui sont très similaire au cas de la dimension finie, et les géométries faibles, o`u les topologies ne co¨ıncident pas. Le cas faible est plus compliqué, car la connection de Levi-Civita et l’application exponentielle n’existent pas toujours. Toutefois, plusieurs variétés riemanniennes faibles ont été abordées avec succès, avec à la clé de nombreuses applications inspirées de travaux d’Arnold [Arn66], comme l’existence et l’unicité de solutions pour de nombreuses équations aux dérivées partielles conservatives.
C’est le cas de l’équation d’Euler des fluides incompressible sur des variétés compactes [EM70], ou de l’équation KdV (voir [KW09] et ses références pour de nombreux autres exemples). En analyse de formes, ces structures ont été particulièrement étudiées par Michor, Mumford et Bauer [BHM11, BHM12, MM07]. La géométrie sous-riemannienne en dimension infinie, en revanche, est encore méconnue. On a bien un r ́esultat de contrôlabilit ́e approch ́e [DS80] dans le cas g ́en ́eral, et de contrôlabilit ́e exacte sur le groupe des difféomorphismes d’une variété compacte engendrés par des flots de champs de vecteurs appartenant à un certain C∞(M)-module [AC09]. Un preprint [GMV12] présente également un équivalent à l’équation hamiltonienne des géodésiques dans un cas très restreint. Mˆeme la géométrie sous-riemannienne forte de dimension infinie est déjà très différente du cas de la dimension finie.
En effet, l’ensemble des points accessibles est le plus souvent un sous-ensemble dense d’intérieur vide, ce qui implique que l’image de la différentielle de l’application point final E est aussi dense d’intérieure vide. En conséquence, il est impossible d’établir un principe du maximum de Pontryagin tel que le théorème 0.1, c’est-à-dire de donner des conditions nécessaires complètes pour qu’une courbe soit une géodésique en toute généralité. En revanche, on peut se contenter d’une réciproque partielle du principe de Pontryagin qui, elle, reste valide en dimension infinie. Cela permet, à défaut de trouver toutes les géodésiques, d’en produire un grande quantité. On montre en effet le résultat suivant.
Espaces de formes contraintes et algorithmes de minimisation
Ces chapitres résument [ATTY13], reproduit en partie I au chapitre 10. On travaille dans le cadre simplifié suivant : M = R d , S est un espace de formes d’ordre ` ∈ N ∗ dans R d et un ouvert d’un espace de Banach B. Rappelons qu’en analyse des déformations, on veut minimiser une fonctionnelle J(X) = 1 2 Z 1 0 hX(t), X(t)i dt + g(ϕ X(1) · q0) = 1 2 Z 1 0 hX(t), X(t)i dt + g(q(1)), o`u X(t) ∈ V un espace de Hilbert de champs de vecteurs, ˙q(t) = ξq(t)X(t) et q(0) = q0 ∈ S. De nombreuses techniques ont été développées pour minimiser J dans diverses exemples d’espaces de formes. Toutefois, l’idée d’ajouter sur les déformations possibles des contraintes qui varient avec la forme est nouvelle. Le seul papier la mentionnant est [You12], et il s’agit de contraintes très particulières, essentiellement utilisées en dimension finie. On souhaite donc étendre notre cadre à des espaces de formes avec contraintes.
Autrement dit, minimiser J sur les champs de vecteurs t 7→ X(t) tels que Cq(t)X(t) = 0 pour presque tout t, ou Cq ∈ L(B, Y ). On voit tout d’abord que sous des hypothèses très faibles, un minimiseur contraint de J existe toujours. Théorème 0.6. Si V a une inclusion continue dans Hs+1(R d , R d ) avec s > d/2 + `, et si q 7→ Cq est continue, g est semicontinue inférieurement et bornée inférieurement, et q0 est à support compact, alors J possède au moins un minimiseur contraint. On prouve également que dans la plupart des cas, on peut approcher le minimiseur contraint de J par une suite de minimiseurs contraints sur des espaces de formes dimension finie. Le résultat principal de ce chapitre est une version adaptée à la dimension infinie du principe du maximum de Pontryagin avec contraintes. Pour cela, notons KV : V ∗ → V l’isométrie induite par h·, ·i, et Kq = ξqKV ξ ∗ q : T ∗ q M → TqM. L’opérateur Kq se calcule facilement à l’aide du noyau reproduisant de V . On définit le hamiltonien du problème h : T ∗S → R par h(q, p) = inf X∈ker Cq pξqX − 1 2 hX, Xi . Théorème 0.7. Supposons que pour tout q ∈ S, Cq soit surjective. Soit X ∈ L 2 (0, 1; V ) un minimiseur de J le long des contraintes.
Alors il existe p ∈ H1 (0, 1; B∗ ) et λ ∈ L 2 (0, 1; Y ∗ ) tels que p(1) + dg(q(1)) = 0 et ( q˙(t) = ∂ph(q(t), p(t)) = Kq(t)p(t) − ξq(t)KV C ∗ q(t)λ(t), p˙(t) = −∂qh(q(t), p(t)) = − ∂q(Kq(t)X(t))∗ p(t) + ∂q(Cq(t)X(t))∗λ(t) presque partout sur [0, 1], o`u le multiplicateur de Lagrange λ est donné par Cq(t)KV C ∗ q(t)λ(t) = Cq(t)KV ξ ∗ q(t) p(t), et X(t) = KV (ξ ∗ q(t) p(t) − C ∗ q(t)λ(t)). Dans le chapitre 7, on développe et on compare deux types d’algorithmes pour la minimisation de J, à la fois dans le cas avec contraintes et le cas sans contraintes. La première méthode consiste à minimiser directement J sur l’espace total des contrôles L 2 (0, 1; V ). La deuxième, appelée méthode du tir, utilise le principe du maximum de Pontryagin pour se restreindre aux géodésiques partant de q0, paramétrées par T ∗ q0 S. Le premier algorithme requiert le calcul du gradient de J sur l’espace L 2 (0, 1; V ). On obtient ∇J(X)(t) = X(t) − KV ξ ∗ q(t) p(t), X ∈ L 2 (0, 1; V ), o`u ˙q(t) = ξq(t)X(t), et p ∈ H1 (0, 1; B∗ ) tel que p(1)+dg(q(1)) = 0 et ˙p(t) = − 1 2 ∂q(p(t)Kq(t)p(t)) presque partout. Un algorithme du gradient permet alors d’obtenir un minimiseur sans contraintes.
Il s’agit d’une formule connue, utilisée avec succès sur de nombreux exemples dans les méthodes LDDMM [TY11, You10]. L’avantage de notre version est qu’elle est valide pour tous les espaces de formes, et mˆeme pour des problèmes de contrôle optimal plus généraux. S’il y a des contraintes, c’est plus compliqué, et comme les formes contraintes n’ont jamais été étudiées, nous avons appliqué la méthode du lagrangien augmenté [IK08]. Plus précisément, il s’agit minimiser successivement des fonctionnelles de la forme JA(X, λr, µr) = 1 2 Z 1 0 hX(t), X(t)i − λr(t)Cq(t)X(t) + 1 2µr kCq(t)X(t)k 2 Y dt + g(q(1)), o`u λr ∈ L 2 (0, 1; Y ) et µr > 0. On obtient, par exemple en utilisant un algorithme du gradient, un contrôle optimal Xr sans contraintes et une trajectoire qr associée qui minimisent cette fonctionnelle. On met ensuite à jour λr+1 = λr − 1 µr CqrXr et µr+1 < µr. Alors, sous des hypothèses appropriées, la suite Xr converge vers un contrôle optimal contraint X, et la suite λr vers le multiplicateur de Lagrange associé à ces contraintes.
Le second algorithme développé, la méthode du tir, utilise notre version du principe du maximum de Pontryagin [PBGM64, AS04, Tré08], selon lequel un minimiseur contraint de J satisfait ∃p ∈ H1 (0, 1; B ∗ ), ( ˙q(t), p˙(t)) = ∇ωh(q(t), p(t)) = (∂ph(q(t), p(t)), −∂qh(q(t), p(t)). Cette équation différentielle a une unique solution maximale à condition initiale p0 ∈ T ∗ q0 S = B∗ fixée. Si une telle courbe est définie sur [0, 1] entier, alors il existe un contrôle X ∈ L 2 (0, 1; V ) de trajectoire associée q, tel que J(X) = 1 2 h(q0, p0) + g(q(1)) = J˜(p0). Minimiser J revient donc à minimiser J˜ sur B∗ . Il suffit, là encore, de calculer le gradient de J˜, que l’on peut ensuite minimiser. Cela s’effectue par la résolution d’une équation adjointe. Le calcul du gradient de J˜ est plus lent que celui des fonctionnelles JA, mais permet une satisfaction totale des contraintes, plutôt qu’une satisfaction approchée. Des simulations numériques de ces algorithmes sont ensuite fournies, sur des exemples ayant des applications en anatomie computationnelle. L’analyse des déformations avec contraintes promet de très nombreuses applications.
La plus importante semble ˆetre la modélisation de multi-formes : l’étude simultanée d’objets indépendants, chacun étant déformé par un flot de champ de vecteurs distinct, mais possédant des intéractions les uns avec les autres. Par exemple, pour étudier simultanément l’hypothalamus et l’hippocampe de différents cerveaux humains, il préférable de les modéliser comme deux objets distincts. Toutefois, ces objets sont immergés dans un milieu continu (le cerveau lui-mˆeme), ce qui les empˆeche de se toucher. On peut modéliser cette configuration par des contraintes sur des multi-formes, voir section 6.2..
Introduction |
