Une méthode d'estimation conjointe de forme et de mouvement non rigide à partir de séquences temporelles tri-dimensionnelles est proposée. Basée sur une mise en correspondance surface-volume, elle permet, à partir d'une unique forme segmentée, d'estimer la forme et son mouvement sur toute la séquence. Une modélisation Markovienne combinée à un algorithme de recuit simulé estime les correspondances entre les noeuds du maillage modélisant l'objet à l'instant t et les voxels du volume à l'instant t + 1. La méthode a été appliquée à l'extraction de formes et de mouvements cardiaques en tomodensitométrie multibarrette. Les tests, réalisés sur données simulées et données réelles, ont donné des résultats prometteurs

COATRIEUX, Jean-Louis

GARREAU, Mireille

LE BRETON, Hervé

SIMON, Antoine

I-Revues

Modélisation Markovienne pour l’Estimation Combinée de Forme et
de Mouvement : Application au Coeur en Imagerie Scanner
Multibarrette.
Antoine SIMON1, Mireille GARREAU1, Jean-Louis COATRIEUX1, Hervé LE BRETON1,2
1Laboratoire Traitement du Signal et de l’Image, INSERM U642,
Université de Rennes 1, Campus de Beaulieu, Bât. 22, 35042 Rennes Cedex, France
2Centre Cardio-Pneumologique, CHU Pontchaillou, 35033 Rennes, France
antoine.simon@univ-rennes1.fr, mireille.garreau@univ-rennes1.fr
jean-louis.coatrieux@univ-rennes1.fr, herve.le.breton@chu-rennes.fr
Résumé – Une méthode d’estimation conjointe de forme et de mouvement non rigide à partir de séquences temporelles tri-dimensionnelles est
proposée. Basée sur une mise en correspondance surface-volume, elle permet, à partir d’une unique forme segmentée, d’estimer la forme et son
mouvement sur toute la séquence. Une modélisation Markovienne combinée à un algorithme de recuit simulé estime les correspondances entre
les noeuds du maillage modélisant l’objet à l’instant t et les voxels du volume à l’instant t + 1. La méthode a été appliquée à l’extraction de
formes et de mouvements cardiaques en tomodensitométrie multibarrette. Les tests, réalisés sur données simulées et données réelles, ont donné
des résultats prometteurs.
Abstract – We propose a method for joint surface and non-rigid motion estimation from three-dimensional dynamic sequences. Based on
a surface-volume matching, it provides, from one unique segmented surface, both the surface and its motion along the whole sequence. A
Markovian Model, combined with a simulated annealing process, estimates the correspondences between the nodes of the mesh modeling the
object of interest at t time and the voxels of the volume at t + 1 time. The method has been applied to cardiac surface and motion extraction in
Multislice Computed Tomography. Tests realized on simulated and real data have provided promising results.
1 Contexte
L’estimation du mouvement cardiaque à partir de séquences
d’images tri-dimensionnelles a suscité un grand intérêt de la
part de la communauté scientifique au cours des dernières an-
nées. D’un grand apport diagnostique dans la lutte contre les
cardiopathies, elle a bénéficié de progrès importants réalisés en
imagerie dynamique, avec notamment les développements du
ciné-IRM, de l’IRM taggué et des scanners multibarrettes.
Les méthodes d’estimation de mouvement 3D à partir de sé-
quences tri-dimensionnelles peuvent être classées en trois prin-
cipaux types : les méthodes de flot optique [1, 2], les modèles
déformables [3–5], et les méthodes de mise en correspondance
[6,7]. Les méthodes de flot optique, si elle permettent d’obtenir
un champ de déplacements dense, reposent sur une contrainte
(la conservation de l’intensité au cours du temps) difficile à
valider dans le cas de la post-synchronisation combinée à la
diffusion du produit de contraste. Les méthodes basées sur un
modèle déformable paramétrique permettent de décrire la sur-
face et/ou le mouvement grâce à quelques paramètres. Elles
nécessitent cependant l’apport d’une connaissance a priori qui
permet de contraindre le résultat mais n’est pas sans risque dans
le cas de phénomènes pathologiques importants. Les méthodes
de mise en correspondance, basées sur la mise en correspon-
dance d’entités suivant des paramètres descriptifs, permettent
une description fine à la fois de la surface et du mouvement.
Dans la plupart des méthodes de mise en correspondance in-
troduites, les entités utilisées sont issues d’une segmentation
préalable de tous les volumes de la séquence. La qualité du ré-
sultat est alors directement liée à la cohérence temporelle de
la segmentation, elle-même difficile à assurer. Contrairement
à ces approches, la méthode proposée ne nécessite la segmen-
tation que du premier instant de la séquence. En effet celle-ci
établit les correspondances entre la surface à l’instant t et le
volume à l’instant t+1. Plus précisément, la recherche de cor-
respondances est réalisée entre, à l’instant t, les noeuds d’un
maillage issu d’un processus de segmentation et de reconstruc-
tion de surface et, à l’instant t+1, les voxels du volume original
après application d’un détecteur de contours.
Ce processus de mise en correspondance, appliqué itérati-
vement, permet conjointement de détecter la surface de l’objet
dans les différents volumes de la séquence et d’estimer, entre
instants successifs, le champ de déplacements 3D associé.
Le processus de mise en correspondance est basé sur la défi-
nition d’un modèle de Markov qui, combiné à un algorithme de
recuit simulé, permet d’obtenir les correspondances optimales.
La méthode, testée sur données simulées et données réelles a
donné des résultats permettant de valider la démarche propo-
sée.
2 Méthode
Un modèle de champ de Markov est défini afin d’estimer
les déplacements entre deux instants successifs. Un champ de
déplacements 3D est alors considéré comme une réalisation
noeud central
noeud voisin
voisinage de premier ordre
voisinage de deuxième ordre
voisinage de troisième ordre
FIG. 1: Voisinage d’un noeud
f = {fi/i = 1, ..., NS} du champ de Markov
F = {Fi/i = 1, ..., NS} regroupant l’ensemble des configu-
rations possibles (NS correspondant au nombre de noeuds du
maillage à t).
Les sites du champ de Markov sont les noeuds du maillage
correspondant à l’instant t. Ce maillage est obtenu par l’appli-
cation au volume original d’un processus de segmentation et de
reconstruction de surface. L’ensemble de ces sites est noté S.
Les étiquettes affectées à ces sites sont les voxels du volume
à t + 1 en correspondance avec les noeuds du maillage à t.
Le champ F est un champ de Markov selon un système de
voisinage V , basé sur la connexité des noeuds du maillage (cf.
fig. 2). Ce système de voisinage vérifie :
∀ i ∈ S, i /∈ V (i) (1)
∀ {i, j} ∈ S, j ∈ V (i)⇔ i ∈ V (j) (2)
P (f) étant la probablilité associée à la réalisation f , les pro-
priétés Markoviennes suivantes sont vérifiées :
P (f) > 0 ∀ f ∈ F (3)
et :
P (fi|fS−{i}) = P (fi|fV (i)) . (4)
D’après le théorème d’Hammersley-Clifford, si le champ aléa-
toire F est un champ de Markov, alors :
P (f) =
1
Z
exp(−U(f)) (5)
où Z est une constante de normalisation et U(f) une énergie
globale à minimiser pour obtenir la réalisation la plus probable.
Cette énergie est définie pour une réalisation f en relation avec
l’observation d par :
U(f, d) = α1Ud(f, d) + α2Ur(f) (6)
Ud(f, d) est une énergie de correspondance globale, U r(f) est
l’énergie interne du champ, α1 et α2 sont des facteurs de pon-
dération.
L’énergie de correspondance globale Ud(f, d), ou énergie
d’adéquation entre observation d et réalisation f , est constituée
de la somme pondérée, sur l’ensemble des sites S, de termes
permettant d’estimer la correspondance entre le site i et l’éti-
quette associée fi:
Ud(f, d) =
∑
i∈S
( αd.Edist(i, fi)
+αt.Etopol(i, fi) (7)
+αc.Econtour(fi)) (8)
  
  
  
  




  
  
  
  




voxel n’appartenant pas à un contour
voxel appartenant à un contour
noeud voisin
noeud central
  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




  
  
  
  




Volume (instant t+1)
Noeud (instant t)
E_topo = 0
E_topo = 3
E_topo = 2
FIG. 2: Exemples, en deux dimensions, du calcul de
Etopol(i, fi).
αd, αt et αc permettent de moduler l’influence de chacun
des trois termes.
Edist(i, fi) correspond à la distance euclidienne entre
le noeud considéré et le voxel correspondant et permet de pri-
vilégier les correspondances les plus proches.
Econtour(fi) et Etopol(i, fi) utilisent le résultat de l’appli-
cation d’un détecteur de contours (implémenté par un filtre de
Canny) au volume à t+ 1.
Econtour(fi) permet ainsi de privilégier les voxels apparte-
nant à un contour :
Econtour(fi) =
{
0 si fi /∈ contour
1 si fi ∈ contour
Etopol(i, fi) transpose la topologie locale du noeud consi-
déré au niveau du voxel afin d’estimer la correspondance to-
pologique entre entités (cf. fig. 2 pour une illustration, en deux
dimensions et avec un voisinage limité à deux noeuds, du calcul
de cette énergie).
L’ensemble de toutes les cliques, noté C pour l’ensemble du
maillage, est défini tel que toute clique c ∈ C est composée de
deux sites suivant la propriété suivante :
pour c = {s, t} ∈ S2, t ∈ V (s) (9)
L’énergie de régularisation Ur(f) est alors donnée par:
Ur(f) =
∑
c={s,t}∈C
||−→fs −−→ft ||
dist(s, t)
(10)
avec dist(s, t) la distance euclidienne entre les noeuds s et t et−→
fs (respectivement −→ft ) le vecteur mouvement estimé au site s
(respectivement t).
La réalisation f la plus probable est obtenue par minimisa-
tion de l’énergie U(f) (cf. éq. 5). Cette minimisation est réali-
sée grâce à un algorithme de recuit simulé : à chaque itération,
une nouvelle correspondance est évaluée ; en cas de diminution
de l’énergie, la correspondance testée est acceptée ; dans le cas
inverse, elle est soit acceptée soit rejetée suivant un choix aléa-
toire et un paramètre de température, décroissant dans le temps
et favorisant l’acceptation de ce type de correspondances au dé-
but du processus. L’optimisation se termine lorsque le système
est stabilisé.
(a) (b)
(c) (d)
FIG. 3: Amplitude du mouvement simulé (a,c) et estimé (b,d)
suivant deux incidences ((a,b) : oblique antérieure gauche 45 ◦
caudale 25◦, (c,d) : oblique antérieure droite 30◦) (couleurs :
en blanc : mouvements dirigés vers l’extérieur des cavités (ex-
pansion), en noir : mouvements dirigés vers l’intérieur des ca-
vités (contraction)), échelle donnée en voxels.
3 Résultats
3.1 Données simulées
La méthode proposée a été testée sur plusieurs simulations
numériques. Afin d’appliquer la méthode de mise en corres-
pondance surface-volume, une surface (correspondant à l’ins-
tant t) doit être déformée (pour obtenir la surface à l’instant
t + 1) puis insérée dans un volume. Ce volume doit être traité
au préalable afin de simuler au mieux un volume réel après ap-
plication du détecteur de contours.
Plus précisément, la surface à l’instant t est une surface réelle
obtenue par segmentation puis reconstruction de surface. Elle
est soumise à des déformations locales et globales (translation,
torsion, rotation, expansion/compression) afin d’obtenir la sur-
face à l’instant t + 1.
Pour obtenir le volume dans lequel cette surface sera trans-
posée, le volume original correspondant à l’instant t est filtré
par le détecteur de contours puis le contour correspondant à la
surface préalablement segmentée est supprimé.
Enfin, la surface déformée est insérée dans ce volume en mo-
difiant la valeur des voxels ayant les mêmes coordonnées que
les noeuds du maillage.
Cette technique a permis d’évaluer précisément la méthode
proposée et notamment de contrôler l’impact du choix des dif-
férents paramètres impliqués que ce soit dans le calcul des
énergies ou dans la procédure de minimisation. La figure 3
illustre, suivant deux incidences, un exemple de résultats. En
couleur est représentée l’amplitude des mouvements appliqués
(a,c) et celle des mouvements estimés (b,d). En blanc sont re-
présentés les mouvements dirigés vers l’extérieur du volume et
en noir les mouvements dirigés vers l’intérieur. L’erreur moyen-
ne d’estimation a été calculée et correspond, sur cet exemple, à
un voxel.
(a) (b)
FIG. 4: Coupes extraites d’un même volume original.
(a) (b)
FIG. 5: Surface reconstruite après la segmentation du volume
correspondant au premier instant de la séquence (incidences :
(a) : oblique antérieure gauche 45◦ caudale 25◦, (b) : oblique
antérieure droite 30◦).
3.2 Données réelles
La méthode proposée a été appliquée sur données réelles, ac-
quises sur un scanner multibarrette Siemens Somatom. La base,
constituée de dix volumes représentant un cycle cardiaque, a
une résolution spatiale de 0.3 × 0.3 × 0.6 mm (cf. fig. 4). Le
premier volume a été segmenté par un processus de croissance
de région limité par un gradient [8] puis la surface extraite re-
construite par l’algorithme des Marching Cubes (cf. fig. 5). A
partir de cette première segmentation, la méthode d’estimation
conjointe de la forme et du mouvement a été appliquée à l’en-
docarde ventriculaire gauche, à une résolution inférieure, sur
toute la séquence après application d’un détecteur de contour
aux volumes correspondant aux instants suivants.
Les figures (6(a), 6(d)) et (6(b),6(e)) illustrent, sous deux in-
cidences, les résultats obtenus à deux instants consécutifs cor-
respondant à la fin de la diastole ventriculaire. Des mouvements
dirigés vers l’extérieur de la cavité ventriculaire (représentés en
blanc), caractéristiques de la diastole, sont détectés au niveau
du ventricule, notamment pour le deuxième instant correspon-
dant à la phase de systole auriculaire. Les figures (6(c), 6(f))
illustrent le début de la phase de systole ventriculaire, avec
la présence de déplacements dirigés vers l’intérieur de la ca-
vité ventriculaire, pouvant correspondre à des mouvements de
contraction.
Afin d’évaluer le processus d’estimation de forme au cours
du temps, les résultats ont été comparés avec des surfaces, cor-
respondant aux mêmes instants, mais obtenues avec l’outil de
segmentation utilisé, dans notre approche, pour le premier vo-
(a) (b) (c)
(d) (e) (f)
FIG. 6: Amplitude du mouvement estimé à trois instants consécutifs (de la fin de la diastole ventriculaire ((a,d) et (b,e)) au
début de la systole ventriculaire (c,f)) (couleurs : en blanc : mouvements dirigés vers l’extérieur des cavités (expansion), en noir :
mouvements dirigés vers l’intérieur des cavités (contraction)) (incidences : (a,b,c) : oblique antérieure gauche 45 ◦ caudale 25◦,
(d,e,f) : oblique antérieure droite 30◦).
lume de la séquence. On considère bien sûr que ces surfaces
segmentées ne peuvent être exploitées comme des références
mais plutôt comme des indicateurs. En traitant l’estimation de
surface sur l’ensemble de la séquence, les tests ont montré une
différence moyenne de 1.1 mm entre, sur les différents instants,
surface estimée par notre méthode et surface directement seg-
mentée.
4 Conclusion et perspectives
Une méthode d’estimation combinée de la position d’un ob-
jet et de son mouvement dans une séquence 3D est proposée.
Elle est ici appliquée à la détection de cavités myocardiques en
imagerie scanner multibarrette pour l’évaluation de la cinétique
cardiaque. La méthode repose sur une mise en correspondance
surface-volume formalisée dans un cadre markovien. Des tests
réalisés sur données simulées et sur données cardiaques réelles
ont montré des résultats prometteurs. De futurs travaux por-
teront sur la recherche de moyens de visualisation adaptés aux
besoins cliniques et sur l’extraction de paramètres cliniquement
significatifs. Ils porteront également à la fois sur l’amélioration
de la méthode (notamment par le test d’autres méthodes de mi-
nimisation) et sur son évaluation à partir d’un nombre plus im-
portant de données réelles.
Références
[1] J. M. Gorce, D. Friboulet, and I. E. Magnin, “Estimation
of three-dimensional cardiac velocity fields: assessment of
a differential method and application to three-dimensional
CT data,” Medical Image Analysis, vol. 1, pp. 245–261,
Apr. 1997.
[2] S. M. Song and R. Leahy, “Computation of 3-D velocity
fields from 3-D cine CT images of a human heart,” IEEE
Trans. Med. Imag., vol. 10, pp. 295–306, Sept. 1991.
[3] A. F. Frangi, W. J. Niessen, and M. A. Viergever, “Three-
dimensional modeling for functional analysis of cardiac
images, a review,” IEEE Trans. Med. Imag., vol. 20, pp.
2–25, Jan. 2001.
[4] C. W. Chen, T. S. Huang, and M. Arrott, “Modeling, ana-
lysis, and visualization of left ventricle shape and motion
by hierarchical decomposition,” IEEE Trans. Pattern Anal.
Machine Intell., vol. 16, pp. 342–356, Apr. 1994.
[5] S. Benayoun and N. Ayache, “Dense non-rigid motion es-
timation in sequences of medical images using differential
constraints,” Int. Journal of Computer Vision, vol. 26, pp.
25–40, Jan. 1998.
[6] C. Kambhamettu, D. Goldgof, M. He, and P. Laskov,
“3D nonrigid motion analysis under small deformations,”
Image and Vision Computing, vol. 21, pp. 229–245, Jan.
2003.
[7] P. Shi, A. J. Sinusas, R. T. Constable, E. Ritman, and J. S.
Duncan, “Point-tracked quantitative analysis of left ven-
tricular surface motion from 3-D image sequences,” IEEE
Trans. Med. Imag., vol. 19, pp. 36–50, Jan. 2000.
[8] H. Guillaume and M. Garreau, “Segmentation de cavités
cardiaques en imagerie scanner multi-barettes,” in 12ème
Forum des Jeunes Chercheurs en Génie Biologique et Mé-
dical (2003), Nantes, France, may 2003, pp. 92–93.


Modélisation Markovienne pour l'estimation combinée de forme et de mouvement : Application au coeur en imagerie scanner multibarrette

Abstract

Similar works

Full text

Available Versions

I-Revues