=Paper=
{{Paper
|id=Vol-2133/rjcia-paper6
|storemode=property
|title=Inférence bayésienne pour l’estimation de déformations larges par champs gaussien : application au recalage d’images multi-modales(Bayesian inference for estimation of large deformations by Gaussian random fields: application to multimodal image registration)
|pdfUrl=https://ceur-ws.org/Vol-2133/rjcia-paper6.pdf
|volume=Vol-2133
|dblpUrl=https://dblp.org/rec/conf/rjcia/DeregnaucourtSY18
}}
==Inférence bayésienne pour l’estimation de déformations larges par champs gaussien : application au recalage d’images multi-modales(Bayesian inference for estimation of large deformations by Gaussian random fields: application to multimodal image registration)==
https://ceur-ws.org/Vol-2133/rjcia-paper6.pdf
Inférence bayésienne pour l’estimation de déformations larges par champs
gaussien: application au recalage d’images multi-modales
Thomas Deregnaucourt1 Chafik Samir1 Anne-Françoise Yao2
1
LIMOS, CNRS UMR 6620, Université Clermont Auvergne, France
2
LMBP, CNRS UMR 6158, Université Clermont Auvergne, France
thomas.deregnaucourt@uca.fr - chafik.samir@uca.fr - anne.yao@uca.fr
Résumé 1 Introduction
Le problème de recalage d’images consiste à estimer la Le recalage d’images est une méthode qui vise à estimer
déformation globale entre une image source I1 et une im- la transformation, soumise à certaines contraintes, d’une
age cible I2 . Dans ce contexte, nous nous intéressons image source I1 vers une image cible I2 , afin de fusion-
à l’estimation d’un champ de déformation U sur le do- ner leurs informations complémentaires. Cette méthode
maine Ω = [0, 1]2 de I1 sachant U sur un ensemble fini est utilisée dans de nombreux domaines d’applications
de courbes β ∈ Ω. Pour ce faire, nous proposons une [11, 15, 2]. En imagerie médicale on utilise le recalage
nouvelle méthode basée sur des modèles gaussiens pour d’images pour détecter des maladies, valider un traite-
recaler des images multimodales. La méthode proposée ment, comparer les données du patient avec des atlas
commence par résoudre le problème de correspondance anatomiques, etc. [12]. L’estimation de cette déformation
entre les courbes βs puis estime le déplacement sur tout est basée soit sur les intensités, soit sur des caractéristiques
Ω. La solution optimale est calculée à l’aide du maximum géométriques, soit sur les deux [12, 10]. Dans le premier
de vraisemblance et de l’inférence bayésienne. D’après les cas, on cherche une transformation concervant la corre-
résultats obtenus sur des données réelles et simulées, la spondance entre les niveaux de gris, et dans le second la
déformation résultante a l’avantage d’être exacte sur les correspondance entre des points, des courbes, etc.
observations et d’être lisse sur Ω. Pour cet article, nous nous sommes intéressés au problème
de l’endométriose. Cette maladie est provoquée par
Mots Clef l’apparition de muqueuse utérine, aussi appelé endomètre,
en dehors de la cavité utérine. Elle touche approxima-
Recalage d’images, Statistique spatiale, Processus tivement 10% des femmes en âge de procréer, et peut
gaussiens, Inférence bayésienne. provoquer divers symptômes tels que des douleurs pelvi-
ennes chroniques, une dysenterie sévère, une infertilité,
Abstract etc. [3]. Il est alors nécessaire de pouvoir détecter si
une patiente à l’endométriose afin de la traiter efficace-
Image registration aims to estimate the global deformation ment, que ce soit par des antalgiques, des traitements hor-
between a target image I1 and a reference image I2 . In monaux, ou par chirurgie dans les cas les plus sévères.
this context, we will focus on estimating a random field U L’endomètre pouvant pénétrer d’autres tissus et organes,
on the I1 domain Ω = [0, 1]2 based on observations of U les méthodes de détection de l’endométriose se basent sur
on a finite set of curves β ∈ Ω. Indeed, we present a new plusieurs modalités d’images, donnant des informations
multimodal image registration method based on Gaussian complémentaires. Plus précisément, l’échographie permet
random fields. The proposed method first find the optimal d’avoir une estimation de l’infiltration de l’endomètre dans
correspondences between curves βs then estimate the de- d’autres tissus, et l’imagerie par résonance magnétique
formation vector field on Ω. The optimal solution is com- (IRM) une position précise des kystes [3]. La fusion des
puted using Maximum Likelihood and Bayesian inference. données IRM/échographie permet alors d’avoir un diag-
Based on results using both real and simulated data, the re- nostic précis, mais nécessite un recalage entre ces deux
sulting deformation has the advantage of being exact on the modalités.
observations as being sufficiently smooth over the whole Ω. Dans notre cas, les deux images I1 et I2 représentent re-
spectivement l’échographie et l’IRM d’un même organe.
Keywords Cependant, comme le montre la Figure 1, les modalités
d’images ont des distributions d’intensités différentes, ce
Image registration, Spatial statistics, Gaussian process, qui rend inefficace les méthodes de recalage basées sur les
Bayesian inference. intensités. D’autre part, un spécialiste peut extraire le con-
tour des organes présents dans les deux images. Nous al- toute l’image, et la mise à l’échelle n’est pas une nui-
lons alors utiliser ces dernières pour effectuer le recalage. sance. Ainsi, nous cherchons seulement l’invariance à la
Pour ce faire, nous estimerons le champ de déformation re-paramétrisation. Par brièveté, nous ne décrivons le pro-
entre les deux images, à l’aide des champs gaussiens. cessus que pour les courbes ouvertes, mais celui-ci peut
être étendu simplement à des courbes fermés [7].
Soit β : [0, 1] → R2 une courbe ouverte paramétrisée.
On utilise par la suite la représentation square-root velocity
function (SRVF) q de β, défini par:
q : [0, 1] → R2
˙
β(t)
t 7→ q
||β̇(t)||2
(a) Echographie (b) IRM Le mapping β ⇐⇒ (β(0), q) étant une bijection, il est pos-
sible de revenir aux courbes originales en stockant le pre-
Figure 1: Exemple de courbes extraites manuellement is- mier point de ces dernières. On note C l’espace de SRVFs:
sues d’images multimodales représentant les mêmes or- Z 1
ganes (en rouge et en jaune respectivement) : échographie
C = q ∈ L2 ([0, 1], R2 ) | ||q(t)||22 dt = 1
à gauche et IRM à droite. 0
Comme on recherche une représentation des courbes in-
Le reste du papier est organisé de la manière suivante. Dans variantes aux re-paramétrisations, nous allons utiliser des
la section 2 nous formaliserons notre problème, et expli- classes d’équivalence. Nous définissons d’abord le groupe
querons en détail notre méthode de résolution. Puis nous des re-paramétrisations Γ:
présenterons nos résultats, sur des données synthétiques et
réelles, dans la section 3. Enfin, la section 4 conclura ce Γ = {γ : [0, 1] → [0, 1] | γ(0) = 0, γ(1) = 1, 0 < γ̇ < ∞}
papier.
La re-paramétrisation d’une courbe β par γ ∈ Γ est
2 Méthodologie donnée par β ◦ γ, et la SRVF √ de cette courbe re-
paramétrisée est alors (q ◦ γ) γ̇. Ainsi, pour unifier tous
2.1 Formulation du problème les éléments de C représentant la même courbe, √ on défini
nos classes d’équivalence par [q] = {(q ◦ γ) γ̇ | γ ∈ Γ}.
Soient Ω un domaine fermé de R2 , β1 ⊂ I1 les courbes
On note S = C/Γ = {[q], q ∈ C} l’ensemble des
de l’échographie et β2 ⊂ I2 celles de l’IRM. On cherche
classes d’équivalence. Afin de comparer deux courbes,
à estimer un champ de déformation Ψ sur Ω transformant
on impose la métrique L2 à S. Sous la représentation
I1 en I2 . Ce champ doit déformer β1 en β2 , tout en étant
SRVF, la métrique L2 correspond à une métrique élastique
lisse. Pour résoudre ce problème, il nous faut:
sur l’espace original des courbes [8], ce qui permet de
1. Trouver une correspondance optimale entre β1 et β2 . déformer les courbes pour avoir la correspondance opti-
male. Cette déformation optimale entre deux points sur S
2. Estimer une déformation Ψ induite par un champ de est obtenue par le chemin géodésique, et la distance entre
déplacement U , c’est-à-dire telle que: elles est définie par la longueur du chemin.
Afin de voir comment ceci peut résoudre notre problème,
Ψ:Ω → Ω on note q1 et q2 la représentation SRVF respective des deux
X 7→ Ψ(X) = X + U (X) courbes β1 et β2 . Afin de calculer la géodésique entre leurs
classes d’équivalence [q1 ] et [q2 ], on fixe q1 , et on cherche
et vérifiant la contrainte Ψ(β1 ) = β2 la re-paramétrisation √ optimale de q2 en résolvant γ̂ =
arg inf ||q1 − (q2 ◦ γ) γ̇||22 . La re-paramétrisation γ̂ don-
γ∈Γ
Nous allons tout d’abord nous intéresser au premier
nera alors la correspondance optimale entre les courbes.
problème.
Par la suite, on note {Xi , i = 1, · · · , N } l’ensemble de
2.2 Correspondance optimale entre courbes N points représentant la discrétisation de β1 et {Ui =
U (Xi ), i = 1, · · · , N } les déplacements correspondants
Afin de trouver une correspondance optimale entre les
donnés par β2 ◦ γ̂. Afin de résoudre le second problème,
courbes, nous adaptons les travaux de Srivastava et al. [13].
nous supposons que U est un champ gaussien.
Dans ce papier, les auteurs s’intéressaient à l’analyse des
formes, et cherchaient une invariance aux transformations 2.3 Champ de déformation gaussien
préservant la forme, c’est-à-dire à la translation, la rota- On rappelle qu’un champ gaussien U est défini par:
tion, la mise à l’échelle et la re-paramétrisation. Dans
notre cas, la translation et la rotation sont déjà fixées pour U (X) ∼ N (µ(X), C(X))
où µ(X) et C(X) sont respectivement la moyenne et la où τ Vα,ν = Σθ . Il n’existe cependant pas de solu-
variance de U (X). Une loi gaussienne étant entièrement tion analytique à cette équation, et devons alors utiliser
décrite par sa moyenne et sa variance, il nous suffit de trou- des méthodes numériques pour déterminer θ̂. Pour ce
ver ces derniers pour définir notre champ. Pour ce faire, faire, nous avons choisi de comparer les méthodes de
nous supposons tout d’abord que U est un champ station- Nelder-Mead [9], la descente du gradient et de Newton.
naire, c’est-à-dire que µ(X) = µ, ∀X ∈ Ω. On a alors: L’optimisation sur ν étant difficile dans notre cas, nous
avons estimé ce paramètre par validation croisée.
N (µ(X), C(X)) = µ + N (0, C(X))
Estimation par inférence bayésienne. Malgré
ce qui implique que µ est une translation sur l’image I1 . l’utilisation de méthodes itératives, notre estimation
Nous pouvons alors supposer que µ = 0. Nous estimons de l’estimateur du maximum de vraisemblance peut
ensuite C(X) par une méthode paramétrique. Il existe un converger vers des maximums locaux. Afin d’éviter cela,
grand choix de fonctions de covariance candidates dans la nous avons choisi d’utiliser d’autres estimateurs basés sur
littérature [14, 1]. Dans ce travail, nous avons choisi C l’inférence bayésienne. On va alors trouver des estimateurs
comme fonction de covariance de Matérn : à partir de la loi a posteriori de nos paramètres f (θ | UX ).
21−ν Cette dernière est construite, avec une loi a priori π(θ) sur
C(h) = τ (hα)ν Kν (hα) (1)
Γ(ν) nos paramètres, à l’aide de la règle de Bayes:
où h est la corrélation spatiale et K la fonction de Bessel f (UX | θ)π(θ)
modifiée de seconde espèce. Généralement, τ > 0 est f (θ | UX ) = ∝ L(θ | UX )π(θ)
π(UX )
appelé le paramètre de variance (marginale), α > 0 le
paramètre d’échelle et ν > 0 le paramètre de lissage. Si Cependant, la loi de la densité a posteriori de nos
ν = 21 + k, k ∈ N, l’équation 1 se réduit au produit d’une paramètres n’étant pas calculable, nous allons alors
exponentielle et d’un polynôme[5]: échantillonner cette distribution. Pour ce faire nous util-
k isons une méthode de Monte-Carlo par chaı̂ne de Markov
X (k + l)! k (MCMC): l’algorithme de Metropolis-Hastings [6]. Ce
C(h) = τ e−hα (2hα)k−l
(2k)! l dernier construit, à partir d’une loi de proposition q(. | θ),
l=0
une chaine de Markov de la manière suivante.
Lorsque ν = 12 , C est la fonction de covariance expo-
nentielle, et elle devient gaussienne pour ν = +∞. De • Choisir θ1 ∼ π(θ)
manière plus générale, pour ν = 12 + k, le champ U sera
• Créer θt+1 à l’aide de θt
de classe C k . Par conséquent, définir notre champ re-
vient à estimer l’hyperparamètre θ = (τ, α, ν) de la fonc- 1. Générer θ∗ ∼ q(. | θt )
tion de covariance. On note Cθ la fonction de covariance
2. Calculer la probabilité d’acceptation p:
d’hyperparamètre θ, et Σθ la matrice de covariance as-
sociée aux points {Xi , i = 1, · · · , N }, donnée par: π(θt+1 )L(θt+1 | UX )q(θt | θt+1 )
p = min 1,
Cθ (||X1 − X1 ||) · · · Cθ (||X1 − XN ||)
π(θt )L(θt | UX )q(θt+1 | θt )
Σθ =
.. .. ..
. . . 3. Choisir θt+1 = θ∗ avec probabilité p, sinon
Cθ (||XN − X1 ||) · · · Cθ (||XN − XN ||) choisir θt+1 = θt
Afin de construire notre chaine de Markov, nous devons
2.4 Estimation des paramètres de la fonction définir les lois de π(θ) et (q(θ | .). Pour ce faire, nous
de covariance supposons tout d’abord que ces lois sont séparables, c’est-
Estimation par maximum de vraisemblance. Le pre- à-dire que q(θ | .) = q(τ | .)q(α | .)q(ν | .) et π(θ) =
mier estimateur considéré est obtenu par maximum de π(τ )π(α)π(ν). N’ayant que peu d’informations a priori
vraisemblance. Dans notre modèle, la fonction de vraisem- sur nos paramètres, nous choisissons de mettre des lois peu
blance est définie par: informatives, dont un résumé est présenté en Tableau 1.
T −1
1 U Σ
X θ
UX
L(θ | UX ) = f (UX | θ) = N/2 1/2
e− 2 Lois a priori π(.) Lois de proposition q(.|θ̃)
(2π) |Σθ |
α ∼ U[0, 1] α ∼ U[α̃ − 0.05, α̃ + 0.05]
où UX = U1 · · · UN . Pour ce faire, nous de-
vons minimiser la fonction de log-vraisemblance négative, τ ∼ U[0, 500] τ ∼ U[τ̃ − 50, τ̃ + 50]
donnée par: ν ∼ U[[ 12 , · · · , 11 ν ∼ U[[ν̃ − 1, · · · , ν̃ + 1]]
2 ]]
N N ln |Vα,ν |
− log(L(θ | UX )) = ln(2π) + ln(τ ) + Tableau 1: Distribution a priori et de proposition sur le
2 2 2
UXT −1
Vα,ν UX paramètre θ.
+
2τ
Une fois la chaine MCMC construite, nous avons choisi Méthode RMSE DF MaxLap
d’estimer θ̂ à partir de la loi a posteriori f (α | UX ), plutôt
que de f (θ | UX ). Pour ce faire, on estime tout d’abord α̂, Nelder-Mead 0.0713 0.5813 0.2902
Gradient 0.0699 0.5541 0.2889
puis on sélectionne l’hyperparamètre θ̂ correspondant. Il
Newton 0.0698 0.5477 0.2878
existe trois estimateurs bayésiens, que nous utiliserons par
la suite: le maximum a posteriori (MAP), la moyenne et la MAP 0.0718 0.5909 0.2982
médiane. Moyenne 0.0721 0.5947 0.3003
Mediane 0.0718 0.5918 0.2990
2.5 Interpolation sur une nouvelle position
Une fois l’hyperparamètre θ̂ estimé, la fonction de covari- Tableau 2: Evaluation quantitative des différentes
ance Ĉ est connue, et ainsi la loi de U également. Il méthodes d’estimation des paramètres sur un exemple de
nous reste alors à interpoler le déplacement U (X ∗ ) sur données synthétique présenté en Figure 2.
une nouvelle position X ∗ . Pour ce faire, nous allons cal-
culer l’espérance conditionnelle. Comme U est un proces-
sus gaussien, ceci revient à effectuer un krigeage simple[4],
3.3 Application sur des données réelles
défini par:
Contrairement aux données synthétiques, plusieurs
Ĉ(||X ∗ − X1 ||)
courbes peuvent être en correspondance pour le même
U (X ∗ ) = U1 · · · UN Σ̂−1
··· recalage d’images IRM/échographie. La Figure 3 montre
∗
Ĉ(||X − XN || un exemple de ce type de données, dont les résultats
obtenus sont résumés dans le Tableau 3. Sur cet exemple,
les méthodes d’optimisation de Nelder-Mead et de Newton
3 Applications ne donnent pas de bons résultats. Pour la descente du
Pour chaque exemple, on discrétise les courbes en 300 gradient et l’inférence bayésienne, le recalage est lisse et
points. Avant de présenter nos résultats, nous présentons donne de faibles erreurs d’interpolation.
les critères de qualité du recalage. Afin d’étudier la stabilité de notre méthode, nous l’avons
appliqué sur 7 données réelles. Les résultats obtenus
3.1 Critères sur la qualité du recalage sont résumés sur la Figure 4. En moyenne, l’erreur
Un bon recalage doit avoir une faible erreur d’interpolation est faible pour tous les estimateurs, avec un
d’interpolation, et doit être lisse. Pour évaluer la champ de déformation relativement lisse. Les estimateurs
qualité d’interpolation de notre méthode , on utilise 200 bayésiens sont plus performants car ils fournissent un
points pour effectuer le recalage, et les 100 restants pour champ de déformation plus lisse, tout en gardant une petite
l’évaluation. La qualité est alors estimée entre les points erreur d’interpolation.
d’évaluation à l’aide de la racine carré des erreurs en
moyenne quadratique (RMSE). Comme on utilise des Méthode RMSE DF MaxLap
courbes pour effectuer notre recalage, on peut aussi utiliser
ces dernières pour estimer la qualité d’interpolation. Nous Nelder Mead 0.1272 0.1796 0.7550
avons choisi d’utiliser la distance de Fréchet (FD) définie Gradient 0.1353 0.1132 0.4091
par: Newton 0.1268 0.1755 0.7058
MAP 0.1383 0.1200 0.3667
dF (F1 , F2 ) = inf max {d(F1 ◦ γ1 , F2 ◦ γ2 )} Moyenne 0.1380 0.1190 0.3680
γ1 ,γ2 ∈Γ
Mediane 0.1381 0.1201 0.3673
où γ1 et γ2 sont des reparamétrisations. Concernant la
régularité du champ de déformation, nous utilisons la carte Tableau 3: Évaluation quantitative des différentes
de la norme du laplacien. Un maximum de la norme du méthodes d’estimation des paramètres sur un exemple de
laplacien (MaxLap) faible signifie alors que le champ de données réelle présenté en Figure 3.
déformation est lisse.
3.2 Application à des données synthétiques
Afin d’évaluer la performance de notre approche, nous
nous intéressons à des données synthétiques présentant
4 Conclusion
différents degrés de déformation. Un exemple de résultat Nous avons construit un outil de recalage d’images, basé
obtenu est présenté en Figure 2. En ce basant sur la courbe sur les champs gaussiens. Cette méthode est efficace
de la norme du laplacien, on remarque que la déformation car le champ estimé est lisse et admet de faibles erreurs
est lisse et locale. De plus, d’après le Tableau 2, qui résume d’interpolation.
l’évaluation de chaque méthode sur cet exemple, le champ Afin d’améliorer l’estimation des paramètres, et donc la
de déformation admet de faibles erreurs d’interpolation. qualité du recalage, plusieurs approches sont envisagées
(a) (b) (c)
(d) (e) (f)
Figure 2: Exemple de champ de déformation estimé sur des données synthétiques: (a) courbes sources, (b) courbes cible,
(c) courbes sources (en noir) et cible (en bleu) superposées, (d) le champ de déformation connu (rouge) et estimé (bleu), (e)
la grille uniforme déformé par le champ estimé, et (f) la carte de la norme du laplacien. De haut en bas, les résultats sont
obtenus après utilisation de la méthode d’optimisation suivante: Nelder-Mead, descente du gradient, Newton, Map.
(a) (b) (c)
(d) (e) (f)
Figure 3: Exemple de champ de déformation estimé sur un recalage IRM/échographie: (a) courbes sources, (b) courbes cible,
(c) courbes sources (en noir) et cible (en bleu) superposées, (d) le champ de déformation connu (rouge) et estimé (bleu), (e)
la grille uniforme déformé par le champ estimé, et (f) la carte de la norme du laplacien. De haut en bas, les résultats sont
obtenus après utilisation de la méthode d’optimisation suivante: Nelder-Mead, descente du gradient, Newton, Map.
(a) RMSE (b) FD (c) MaxLap
Figure 4: Évaluation quantitative des différentes méthodes d’estimation des paramètres sur 7 données réelles: (1) Nelder-
Mead, (2) descente du gradient, (3) Newton, (4) MAP, (5) moyenne, et (6) médiane. Pour chaque méthode, nous présentons
un boxplot de (a) la RMSE,(b) la distance de Fréchet, et (c) le maximum de la norme du laplacien.
et feront l’objet de futurs travaux. Tout d’abord, pour les [8] W. Mio, A. Srivastava, and S. Joshi. On shape of
méthodes de la descente du gradient et de Newton, une plane elastic curves. International Journal of Com-
approximation des dérivées partielles semble nécessaire puter Vision, 73(3):307–324, 2007.
pour éviter les erreurs numériques. De plus, un pas adap-
tatif pourrait améliorer la convergence de ces algorithmes. [9] J. Nelder and R. Mead. A simplex method for func-
Dans un second temps, on pourrait utiliser des lois a priori tion minimization. The Computer Journal, 7, issue.
plus informatives sur nos paramètres, ainsi que des lois de 4:308–313, 1965.
propositions plus restreintes pour améliorer la convergence [10] K. Rohr. Landmark-based image analysis: Using ge-
des chaines MCMC. ometric and intensity models. Kluwer Academic Pub-
lishing, 2001.
Références
[1] R. J. Adler and J. E. Taylor. Random Fields and Ge- [11] J. E. Roos, D. Weishaupt, S. Wildermuth, J. K. Will-
ometry. Springer Monographs in Mathematics, 2007. mann, B. Marincek, and P. R. Hilfiker. Experience
of 4 years with open mr defecography: pictorial re-
[2] O. Arandjelović, D.-S. Pham, and S. Venkatesh. Ef- view of anorectal anatomy and disease. Radiograph-
ficient and accurate set-based registration of time- ics, 22(4):817–832, 2002.
separated aerial images. Pattern Recognition,
48(11):3466–3476, 2015. [12] A. Sotiras, C. Davatazikos, and N. Paragios. De-
formable medical image registration : A survey. IN-
[3] L. P. Chamié, R. Blasbalg, R. M. A. Pereira, RIA Report, september 2012.
G. Warmbrand, and P. C. Serafini. Findings of pelvic
endometriosis at transvaginal us, mr imaging, and la- [13] A. Srivastava, E. Klassen, S. Joshi, and I. Jermyn.
paroscopy. Radiographics, 31(4):E77–E100, 2011. Shape analysis of elastic curves in Euclidean spaces.
IEEE Transactions on Pattern Analysis and Machine
[4] N. Cressie. Statistics for Spatial Data, Revised Edi- Intelligence, 33:1415–1428, 2011.
tion. Wiley, 1993.
[14] M. L. Stein. Interpolation of Spatial Data. Springer
[5] T. Gneiting, W. Kleiber, and M. Schlather. Matérn Series in Statistics, 1999.
cross-covariance functions for multivariate random
[15] Z. Tu, W. Xie, J. Cao, C. Van Gemeren, R. Poppe,
fields. Journal of the American Statistical Associa-
and R. C. Veltkamp. Variational method for joint op-
tion, 105:1167–1177, 2010.
tical flow estimation and edge-aware image restora-
[6] W. K. Hastings. Monte carlo sampling methods us- tion. Pattern Recognition, 65:11–25, 2017.
ing 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 re-
lated features. Journal of the American Statistical
Association, 107(499):1152–1165, 2012.