Infe´rence baye´sienne pour l'estimation de de´formations larges par champs gaussien: application au recalage d'images multi-modales

-

Image registration aims to estimate the global deformation between a target image I1 and a reference image I2. In this context, we will focus on estimating a random field U on the I1 domain Ω = [ 0, 1 ]2 based on observations of U on a finite set of curves β ∈ Ω. Indeed, we present a new multimodal image registration method based on Gaussian random fields. The proposed method first find the optimal correspondences between curves βs then estimate the deformation vector field on Ω. The optimal solution is computed using Maximum Likelihood and Bayesian inference. Based on results using both real and simulated data, the resulting deformation has the advantage of being exact on the observations as being sufficiently smooth over the whole Ω. Le recalage d’images est une me´thode qui vise a` estimer la transformation, soumise a` certaines contraintes, d’une image source I1 vers une image cible I2, afin de fusionner leurs informations comple´mentaires. Cette me´thode est utilise´e dans de nombreux domaines d’applications [ 11, 15, 2 ]. En imagerie me´dicale on utilise le recalage d’images pour de´tecter des maladies, valider un traitement, comparer les donne´es du patient avec des atlas anatomiques, etc. [ 12 ]. L’estimation de cette de´formation est base´e soit sur les intensite´s, soit sur des caracte´ristiques ge´ome´triques, soit sur les deux [ 12, 10 ]. Dans le premier cas, on cherche une transformation concervant la correspondance entre les niveaux de gris, et dans le second la correspondance entre des points, des courbes, etc. Pour cet article, nous nous sommes inte´resse´s au proble`me de l’endome´triose. Cette maladie est provoque´e par l’apparition de muqueuse ute´rine, aussi appele´ endome`tre, en dehors de la cavite´ ute´rine. Elle touche approximativement 10% des femmes en aˆge de procre´er, et peut provoquer divers symptoˆmes tels que des douleurs pelviennes chroniques, une dysenterie se´ve`re, une infertilite´, etc. [ 3 ]. Il est alors ne´cessaire de pouvoir de´tecter si une patiente a` l’endome´triose afin de la traiter efficacement, que ce soit par des antalgiques, des traitements hormonaux, ou par chirurgie dans les cas les plus se´ve`res. L’endome`tre pouvant pe´ne´trer d’autres tissus et organes, les me´thodes de de´tection de l’endome´triose se basent sur plusieurs modalite´s d’images, donnant des informations comple´mentaires. Plus pre´cise´ment, l’e´chographie permet d’avoir une estimation de l’infiltration de l’endome`tre dans d’autres tissus, et l’imagerie par re´sonance magne´tique (IRM) une position pre´cise des kystes [ 3 ]. La fusion des donne´es IRM/e´chographie permet alors d’avoir un diagnostic pre´cis, mais ne´cessite un recalage entre ces deux modalite´s.

Dans notre cas, les deux images I1 et I2 repre´sentent respectivement l’e´chographie et l’IRM d’un meˆme organe. Cependant, comme le montre la Figure 1, les modalite´s d’images ont des distributions d’intensite´s diffe´rentes, ce qui rend inefficace les me´thodes de recalage base´es sur les intensite´s. D’autre part, un spe´cialiste peut extraire le contour des organes pre´sents dans les deux images. Nous allons alors utiliser ces dernie`res pour effectuer le recalage. Pour ce faire, nous estimerons le champ de de´formation entre les deux images, a` l’aide des champs gaussiens. (a) Echographie (b) IRM

Le reste du papier est organise´ de la manie`re suivante. Dans la section 2 nous formaliserons notre proble`me, et expliquerons en de´tail notre me´thode de re´solution. Puis nous pre´senterons nos re´sultats, sur des donne´es synthe´tiques et re´elles, dans la section 3. Enfin, la section 4 conclura ce papier. 1. Trouver une correspondance optimale entre β1 et β2. 2. Estimer une de´formation Ψ induite par un champ de de´placement U , c’est-a`-dire telle que: Ψ : Ω

→ X 7→ Ω Ψ(X) = X + U (X) et ve´rifiant la contrainte Ψ(β1) = β2 Nous allons tout d’abord nous inte´resser au premier proble`me. 2.2

Correspondance optimale entre courbes Afin de trouver une correspondance optimale entre les courbes, nous adaptons les travaux de Srivastava et al. [ 13 ]. Dans ce papier, les auteurs s’inte´ressaient a` l’analyse des formes, et cherchaient une invariance aux transformations pre´servant la forme, c’est-a`-dire a` la translation, la rotation, la mise a` l’e´chelle et la re-parame´trisation. Dans notre cas, la translation et la rotation sont de´ja` fixe´es pour toute l’image, et la mise a` l’e´chelle n’est pas une nuisance. Ainsi, nous cherchons seulement l’invariance a` la re-parame´trisation. Par brie`vete´, nous ne de´crivons le processus que pour les courbes ouvertes, mais celui-ci peut eˆtre e´tendu simplement a` des courbes ferme´s [ 7 ]. Soit β : [ 0, 1 ] → R2 une courbe ouverte parame´trise´e. On utilise par la suite la repre´sentation square-root velocity function (SRVF) q de β, de´fini par: q : [ 0, 1 ]

→ R2 Le mapping β ⇐⇒ (β(0), q) e´tant une bijection, il est possible de revenir aux courbes originales en stockant le premier point de ces dernie`res. On note C l’espace de SRVFs: Z 1

0 C = q ∈ L2([ 0, 1 ], R2) | ||q(t)||22dt = 1 Comme on recherche une repre´sentation des courbes invariantes aux re-parame´trisations, nous allons utiliser des classes d’e´quivalence. Nous de´finissons d’abord le groupe des re-parame´trisations Γ: Γ = {γ : [ 0, 1 ] → [ 0, 1 ] | γ(0) = 0, γ(1) = 1, 0 < γ˙ < ∞} La re-parame´trisation d’une courbe β par γ ∈ Γ est donne´e par β ◦ γ, et la SRVF de cette courbe reparame´trise´e est alors (q ◦ γ)√γ˙ . Ainsi, pour unifier tous les e´le´ments de C repre´sentant la meˆme courbe, on de´fini nos classes d’e´quivalence par [q] = {(q ◦ γ)√γ˙ | γ ∈ Γ}. On note S = C/Γ = {[q], q ∈ C} l’ensemble des classes d’e´quivalence. Afin de comparer deux courbes, on impose la me´trique L2 a` S. Sous la repre´sentation SRVF, la me´trique L2 correspond a` une me´trique e´lastique sur l’espace original des courbes [ 8 ], ce qui permet de de´former les courbes pour avoir la correspondance optimale. Cette de´formation optimale entre deux points sur S est obtenue par le chemin ge´ode´sique, et la distance entre elles est de´finie par la longueur du chemin.

Afin de voir comment ceci peut re´soudre notre proble`me, on note q1 et q2 la repre´sentation SRVF respective des deux courbes β1 et β2. Afin de calculer la ge´ode´sique entre leurs classes d’e´quivalence [q1] et [q2], on fixe q1, et on cherche la re-parame´trisation optimale de q2 en re´solvant γˆ = arg γin∈fΓ ||q1 − (q2 ◦ γ)√γ˙ ||22. La re-parame´trisation γˆ donnera alors la correspondance optimale entre les courbes. Par la suite, on note {Xi, i = 1, · · · , N } l’ensemble de N points repre´sentant la discre´tisation de β1 et {Ui = U (Xi), i = 1, · · · , N } les de´placements correspondants donne´s par β2 ◦ γˆ. Afin de re´soudre le second proble`me, nous supposons que U est un champ gaussien. ou` μ(X) et C(X) sont respectivement la moyenne et la variance de U (X). Une loi gaussienne e´tant entie`rement de´crite par sa moyenne et sa variance, il nous suffit de trouver ces derniers pour de´finir notre champ. Pour ce faire, nous supposons tout d’abord que U est un champ stationnaire, c’est-a`-dire que μ(X) = μ, ∀X ∈ Ω. On a alors:

N (μ(X), C(X)) = μ + N (0, C(X)) ce qui implique que μ est une translation sur l’image I1. Nous pouvons alors supposer que μ = 0. Nous estimons ensuite C(X) par une me´thode parame´trique. Il existe un grand choix de fonctions de covariance candidates dans la litte´rature [ 14, 1 ]. Dans ce travail, nous avons choisi C comme fonction de covariance de Mate´rn : 21−ν Γ(ν) C(h) = τ (hα)ν Kν (hα) (1) ou` h est la corre´lation spatiale et K la fonction de Bessel modifie´e de seconde espe`ce. Ge´ne´ralement, τ > 0 est appele´ le parame`tre de variance (marginale), α > 0 le parame`tre d’e´chelle et ν > 0 le parame`tre de lissage. Si ν = 12 + k, k ∈ N, l’e´quation 1 se re´duit au produit d’une exponentielle et d’un polynoˆme[ 5 ]:

C(h) = τ e−hα Xk (k + l)! k

(2k)! l l=0 (2hα)k−l Lorsque ν = 21 , C est la fonction de covariance exponentielle, et elle devient gaussienne pour ν = +∞. De manie`re plus ge´ne´rale, pour ν = 12 + k, le champ U sera de classe Ck. Par conse´quent, de´finir notre champ revient a` estimer l’hyperparame`tre θ = (τ, α, ν) de la fonction de covariance. On note Cθ la fonction de covariance d’hyperparame`tre θ, et Σθ la matrice de covariance associe´e aux points {Xi, i = 1, · · · , N }, donne´e par: Σθ =    Cθ(||X1 − X1||) . . .

· · · . . .

Cθ(||XN − X1||) · · ·

Cθ(||X1 − XN ||) 

...  Cθ(||XN − XN ||) 2.4

Estimation des parame`tres de la fonction de covariance Estimation par maximum de vraisemblance. Le premier estimateur conside´re´ est obtenu par maximum de vraisemblance. Dans notre mode`le, la fonction de vraisemblance est de´finie par: L(θ | UX ) = f (UX | θ) = (2π)N/2|Σθ|1/2 e− UXT Σθ−21UX 1 ou` UX = U1 · · · UN . Pour ce faire, nous devons minimiser la fonction de log-vraisemblance ne´gative, donne´e par: − log(L(θ | UX )) N 2 ln(τ ) + ln |Vα,ν | 2 = +

N 2

ln(2π) + UXT Vα−,ν1UX ou` τ Vα,ν = Σθ. Il n’existe cependant pas de solution analytique a` cette e´quation, et devons alors utiliser des me´thodes nume´riques pour de´terminer θˆ. Pour ce faire, nous avons choisi de comparer les me´thodes de Nelder-Mead [ 9 ], la descente du gradient et de Newton. L’optimisation sur ν e´tant difficile dans notre cas, nous avons estime´ ce parame`tre par validation croise´e. Estimation par infe´rence baye´sienne. Malgre´ l’utilisation de me´thodes ite´ratives, notre estimation de l’estimateur du maximum de vraisemblance peut converger vers des maximums locaux. Afin d’e´viter cela, nous avons choisi d’utiliser d’autres estimateurs base´s sur l’infe´rence baye´sienne. On va alors trouver des estimateurs a` partir de la loi a posteriori de nos parame`tres f (θ | UX ). Cette dernie`re est construite, avec une loi a priori π(θ) sur nos parame`tres, a` l’aide de la re`gle de Bayes: f (θ | UX ) = f (UX | θ)π(θ) π(UX ) ∝ L(θ | UX )π(θ) Cependant, la loi de la densite´ a posteriori de nos parame`tres n’e´tant pas calculable, nous allons alors e´chantillonner cette distribution. Pour ce faire nous utilisons une me´thode de Monte-Carlo par chaˆıne de Markov (MCMC): l’algorithme de Metropolis-Hastings [ 6 ]. Ce dernier construit, a` partir d’une loi de proposition q(. | θ), une chaine de Markov de la manie`re suivante.

Choisir θ1 ∼ π(θ) Cre´er θt+1 a` l’aide de θt 1. Ge´ne´rer θ∗ ∼ q(. | θt) 2. Calculer la probabilite´ d’acceptation p: p = min 1, π(θt+1)L(θt+1 | UX )q(θt | θt+1)

π(θt)L(θt | UX )q(θt+1 | θt) 3. Choisir θt+1 = θ∗ avec probabilite´ p, sinon choisir θt+1 = θt Afin de construire notre chaine de Markov, nous devons de´finir les lois de π(θ) et (q(θ | .). Pour ce faire, nous supposons tout d’abord que ces lois sont se´parables, c’esta`-dire que q(θ | .) = q(τ | .)q(α | .)q(ν | .) et π(θ) = π(τ )π(α)π(ν). N’ayant que peu d’informations a priori sur nos parame`tres, nous choisissons de mettre des lois peu informatives, dont un re´sume´ est pre´sente´ en Tableau 1.

Lois a priori π(.)

α ∼ U[ 0, 1 ] τ ∼ U[0, 500] ν ∼ U[[ 21 , · · · , 121 ]]

Lois de proposition q(.|θ˜) α ∼ U[α˜ − 0.05, α˜ + 0.05]

τ ∼ U[τ˜ − 50, τ˜ + 50] ν ∼ U[[ν˜ − 1, · · · , ν˜ + 1]]

Une fois la chaine MCMC construite, nous avons choisi d’estimer θˆ a` partir de la loi a posteriori f (α | UX ), plutoˆt que de f (θ | UX ). Pour ce faire, on estime tout d’abord αˆ, puis on se´lectionne l’hyperparame`tre θˆ correspondant. Il existe trois estimateurs baye´siens, que nous utiliserons par la suite: le maximum a posteriori (MAP), la moyenne et la me´diane. Une fois l’hyperparame`tre θˆ estime´, la fonction de covariance Cˆ est connue, et ainsi la loi de U e´galement. Il nous reste alors a` interpoler le de´placement U (X∗) sur une nouvelle position X∗. Pour ce faire, nous allons calculer l’espe´rance conditionnelle. Comme U est un processus gaussien, ceci revient a` effectuer un krigeage simple[ 4 ], de´fini par:

U (X∗) = U1 · · ·

UN Σˆ −1 Cˆ(||X∗ − X1||) 

· · · Cˆ(||X∗ − XN ||  3

Applications Pour chaque exemple, on discre´tise les courbes en 300 points. Avant de pre´senter nos re´sultats, nous pre´sentons les crite`res de qualite´ du recalage. 3.1

Crite`res sur la qualite´ du recalage Un bon recalage doit avoir une faible erreur d’interpolation, et doit eˆtre lisse. Pour e´valuer la qualite´ d’interpolation de notre me´thode , on utilise 200 points pour effectuer le recalage, et les 100 restants pour l’e´valuation. La qualite´ est alors estime´e entre les points d’e´valuation a` l’aide de la racine carre´ des erreurs en moyenne quadratique (RMSE). Comme on utilise des courbes pour effectuer notre recalage, on peut aussi utiliser ces dernie`res pour estimer la qualite´ d’interpolation. Nous avons choisi d’utiliser la distance de Fre´chet (FD) de´finie par: dF (F1, F2) =

inf γ1,γ2∈Γ

max {d(F1 ◦ γ1, F2 ◦ γ2)} ou` γ1 et γ2 sont des reparame´trisations. Concernant la re´gularite´ du champ de de´formation, nous utilisons la carte de la norme du laplacien. Un maximum de la norme du laplacien (MaxLap) faible signifie alors que le champ de de´formation est lisse. 3.2

Application a` des donne´es synthe´tiques Afin d’e´valuer la performance de notre approche, nous nous inte´ressons a` des donne´es synthe´tiques pre´sentant diffe´rents degre´s de de´formation. Un exemple de re´sultat obtenu est pre´sente´ en Figure 2. En ce basant sur la courbe de la norme du laplacien, on remarque que la de´formation est lisse et locale. De plus, d’apre`s le Tableau 2, qui re´sume l’e´valuation de chaque me´thode sur cet exemple, le champ de de´formation admet de faibles erreurs d’interpolation.

Me´thode

Nelder-Mead Gradient Newton MAP

Moyenne Mediane Contrairement aux donne´es synthe´tiques, plusieurs courbes peuvent eˆtre en correspondance pour le meˆme recalage d’images IRM/e´chographie. La Figure 3 montre un exemple de ce type de donne´es, dont les re´sultats obtenus sont re´sume´s dans le Tableau 3. Sur cet exemple, les me´thodes d’optimisation de Nelder-Mead et de Newton ne donnent pas de bons re´sultats. Pour la descente du gradient et l’infe´rence baye´sienne, le recalage est lisse et donne de faibles erreurs d’interpolation.

Afin d’e´tudier la stabilite´ de notre me´thode, nous l’avons applique´ sur 7 donne´es re´elles. Les re´sultats obtenus sont re´sume´s sur la Figure 4. En moyenne, l’erreur d’interpolation est faible pour tous les estimateurs, avec un champ de de´formation relativement lisse. Les estimateurs baye´siens sont plus performants car ils fournissent un champ de de´formation plus lisse, tout en gardant une petite erreur d’interpolation.

Me´thode

Nelder Mead Gradient Newton MAP

Moyenne Mediane

RMSE Nous avons construit un outil de recalage d’images, base´ sur les champs gaussiens. Cette me´thode est efficace car le champ estime´ est lisse et admet de faibles erreurs d’interpolation.

Afin d’ame´liorer l’estimation des parame`tres, et donc la qualite´ du recalage, plusieurs approches sont envisage´es et feront l’objet de futurs travaux. Tout d’abord, pour les me´thodes de la descente du gradient et de Newton, une approximation des de´rive´es partielles semble ne´cessaire pour e´viter les erreurs nume´riques. De plus, un pas adaptatif pourrait ame´liorer la convergence de ces algorithmes. Dans un second temps, on pourrait utiliser des lois a priori plus informatives sur nos parame`tres, ainsi que des lois de propositions plus restreintes pour ame´liorer la convergence des chaines MCMC.

Re´fe´rences

References

[1] R. J. Adler and J. E. Taylor . Random Fields and Geometry . Springer Monographs in Mathematics, 2007 .

[2] O. Arandjelovic´ , D.-S. Pham , and S. Venkatesh . Efficient and accurate set-based registration of timeseparated aerial images . Pattern Recognition , 48 ( 11 ): 3466 - 3476 , 2015 .

[3] L. P. Chamie ´ , R. Blasbalg , R. M. A. Pereira , G. Warmbrand, and P. C. Serafini . Findings of pelvic endometriosis at transvaginal us, mr imaging, and laparoscopy . Radiographics , 31 ( 4 ): E77 - E100 , 2011 .

[4] N. Cressie . Statistics for Spatial Data, Revised Edition . Wiley, 1993 .

[5] T. Gneiting , W. Kleiber , and M. Schlather . Mate´rn cross-covariance functions for multivariate random fields . Journal of the American Statistical Association , 105 : 1167 - 1177 , 2010 .

[6] W. K. Hastings . Monte carlo sampling methods using markov chains and their applications . Biometrika , 57 ( 1 ): 97 - 109 , 1970 .

[7] S. Kurtek , A. Srivastava , E. Klassen, and Z. Ding . Statistical modeling of curves using shapes and related features . Journal of the American Statistical Association , 107 ( 499 ): 1152 - 1165 , 2012 .

[8] W. Mio , A. Srivastava , and S. Joshi . On shape of plane elastic curves . International Journal of Computer Vision , 73 ( 3 ): 307 - 324 , 2007 .

[9] J. Nelder and R. Mead . A simplex method for function minimization . The Computer Journal , 7 , issue. 4 : 308 - 313 , 1965 .

[10] K. Rohr . Landmark-based image analysis: Using geometric and intensity models . Kluwer Academic Publishing, 2001 .

[11] J. E. Roos , D. Weishaupt , S. Wildermuth , J. K. Willmann , B. Marincek , and P. R. Hilfiker . Experience of 4 years with open mr defecography: pictorial review of anorectal anatomy and disease . Radiographics , 22 ( 4 ): 817 - 832 , 2002 .

[12] A. Sotiras , C. Davatazikos , and N. Paragios . Deformable medical image registration : A survey . INRIA Report , september 2012 .

[13] A. Srivastava , E. Klassen , S. Joshi , and I. Jermyn. Shape analysis of elastic curves in Euclidean spaces . IEEE Transactions on Pattern Analysis and Machine Intelligence , 33 : 1415 - 1428 , 2011 .

[14] M. L. Stein . Interpolation of Spatial Data . Springer Series in Statistics, 1999 .

[15] Z. Tu , W. Xie , J. Cao , C. Van Gemeren , R. Poppe , and R. C. Veltkamp . Variational method for joint optical flow estimation and edge-aware image restoration . Pattern Recognition , 65 : 11 - 25 , 2017 .