FR2877465A1 - Procede d'estimation du mouvement et/ou de la deformation d'au moins un objet dans un ensemble d'images 3d, et procede de recalage tridimensionnel correspondant - Google Patents

Abstract

L'invention concerne un procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D comprenant au moins deux images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas. Selon l'invention, ce procédé comprend les étapes suivantes : a) obtention (12) d'une surface 3D dudit au moins un objet contenu dans l'une des images 3D dudit ensemble, dite image de départ, par traitement de l'image de départ, ladite surface 3D étant représentée sous la forme d'un maillage comprenant des noeuds ; b) choix (11) de l'image de départ comme image courante ; c) détermination (13) de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante dudit ensemble, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation dudit au moins un objet entre les images courante et suivante.

Description

Procédé d'estimation du mouvement et/ou de la déformation d'au moins un

objet contenu dans un ensemble d'images 3D, et procédé de recalage tridimensionnel correspondant.

1. Domaine de l'invention Le domaine de l'invention est celui des images 3D. Plus précisément, l'invention concerne une technique d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas (modèle préexistant).

La technique de l'invention est par exemple mise en oeuvre en imagerie scanner multi-barrette (MSCT), nouvelle avancée technologique en imagerie scanner et imagerie non invasive, permettant d'acquérir des images de hautes résolutions spatiales et temporelles. Une telle mise en oeuvre de l'invention en imagerie scanner multi-barrette présente de gros avantages par rapport à d'autres modalités de l'imagerie cardiaque, dans le sens où la modalité scanner est une technique d'imagerie non invasive, permettant à moindre coût, dans un temps limité d'examen tout en laissant un confort au patient, de visualiser l'ensemble des structures cardiaques d'intérêt en une seule acquisition.

Il est clair cependant que l'invention peut également être mise en oeuvre avec d'autres modalités d'imagerie, telles que l'IRM[ (Imagerie par Résonance Magnétique) et l'échocardiographie, dès lors que ces techniques (toujours en évolution) peuvent fournir une résolution spatio-temporelle suffisante.

En ce qui concerne l'ensemble d'images 3D, il s'agit par exemple: soit d'un ensemble comprenant deux images 3D acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes; soit d'une séquence temporelle d'images 3D comprenant au moins trois images 3D acquises à des instants successifs au sein d'une même plage temporelle ou de plusieurs plages temporelles non adjacentes, et avec une même modalité d'imagerie ou plusieurs modalités d'imagerie.

Dans le domaine médical, l'invention s'applique notamment, mais non exclusivement, à l'extraction et la caractérisation du mouvement cardiaque en imagerie scanner multi-barrette. Il s'agit par exemple d'estimer le mouvement cardiaque des structures du myocarde, pour l'instant limitées aux cavités, mais pouvant s'étendre à la paroi musculaire avec détection à la fois de l'endocarde et de l'épicarde. Cette application relève de l'estimation de mouvements tridimensionnels (3D) de formes complexes et déformables (cavités cardiaques), soumises à des mouvements de nature non-rigide avec pour objectif d'extraire des informations quantitatives de la dynamique cardiaque (fraction d'éjection, variations de volumes, amplitudes de mouvements par région). L'estimation de mouvement 3D est ici conduite en disposant au départ d'une séquence de volumes scanner reconstruits par post-synchronisation sur l'ECG à différents instants du cycle cardiaque, et d'une segmentation semi-interactive des cavités gauches. L'analyse de ces images scanner est soumise à de nouvelles contraintes telles que, en particulier, la présence d'artefacts de mouvement dans les images, la présence de structures fines accessibles de par la bonne résolution spatiale des images et de gros volumes de données à traiter.

Cette application médicale particulière présente un intérêt important. En effet, en terme de santé publique, les maladies cardio-vasculaires représentent aujourd'hui 40% des examens pathologiques réalisés et constituent la première cause de mortalité. L'athérosclérose, principale maladie artérielle caractérisée par des lésions pariétales artérielles ou plaques d'athérome, est à l'origine de phénomènes ischémiques qui se traduisent par des altérations de la dynamique du myocarde. En imagerie cardio- vasculaire, le diagnostic repose sur l'état de la circulation coronaire mais aussi sur l'analyse de la cinétique ventriculaire globale et locale. Au niveau des modalités d'observation et plus particulièrement, en imagerie scanner, l'apparition récente des systèmes multi-barrette (MSCT) (multi-détection, amélioration des résolutions spatiale et temporelle), permet aujourd'hui l'acquisition de structures en mouvement et ouvre une nouvelle voie en ce qui concerne l'imagerie non invasive cardiaque.

Hors du champ médical, la technique de l'invention peut être appliquée à l'estimation de déformations d'autres structures.

Elle peut aussi répondre aux problèmes de recalage tridimensionnel d'objets qui se posent en imagerie médicale, mais également dans d'autres domaines liés à la vision par ordinateur, quand les observations de scènes sont acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie. En effet, dans ce cas, le recalage tridimensionnel nécessite une estimation préalable du mouvement et/ou de la déformation d'un ou plusieurs objets entre les deux images 3D considérées. Ces dernières sont acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas (modèle préexistant). Dans le dernier cas, on cherche à superposer une image réelle avec des données de type atlas (atlas anatomique par exemple), qui permettent de partir d'une image modèle. Ces objets sont par exemple des formes tridimensionnelles déformables soumises à des mouvements non-rigides. C'est bien sûr le cas de nombreux organes du corps humain, pour lesquels le diagnostic ou la conduite de gestes interventionnels en thérapie ou chirurgie assistée par ordinateur nécessite l'acquisition et le traitement d'images 3D. Mais le domaine de la chirurgie esthétique ou réparation tissulaire est aussi concerné.

La technique de l'invention, du fait qu'elle permet d'obtenir des correspondances entre deux objets inscrits dans des volumes tridimensionnels, peut être généralisée à d'autres applications. L'acquisition tridimensionnelle, si elle est particulièrement présente dans le cadre de l'imagerie médicale, est aussi utilisée dans d'autres domaines comme l'imagerie marine (sonar), l'imagerie satellitaire. D'autres domaines tels que la fabrication de matériaux peuvent aussi bénéficier de la présente invention, pour la mesure de déformations de matières souples ou pour l'évaluation de résistances de matériaux sous l'effet de pression par exemple.

2. Art antérieur Il existe de nombreuses techniques connues d'estimation de mouvement. Dans un souci de simplification, on limite ci-après la discussion aux techniques connues d'estimation de mouvement non rigide, type de mouvement dont sont dotées les structures cardiaques. En effet, la plupart des travaux méthodologiques connus relatifs à l'extraction de mouvement sont appliqués à l'analyse du mouvement cardiaque, du fait que le coeur est soumis à des déformations complexes (translation, torsion, compression, expansion, élongation, raccourcissement) qui relèvent de mouvements non-rigides et en particulier élastiques.

On résume ci-dessous ces travaux, toutes modalités d'acquisition d'images confondues (IRM taggué, IRM à contraste de phase, imagerie par ultrasons,

.), en retenant les avantages et inconvénients de chaque type de méthode, et en les situant par rapport à la technique de l'invention. Il est également fait mention de diverses publications, dont les références précises sont fournies ci-après dans une annexe intitulée bibliographie ...DTD: En vision par ordinateur, trois principales méthodes d'estimation de mouvement élastique existent et ont été développées dans la littérature: - les méthodes basées sur les modèles déformables; - les méthodes de flot optique; et - les méthodes de suivi de forme, ou de mise en correspondance.

Dans les méthodes basées sur les modèles déformables, les modèles déformables, modèles géométriques, paramétriques [Bardinet '96, Chen '94] ou non paramétriques [Choi '01, Lin '00, Benayoun '951, après avoir été initialisés, subissent des déformations jusqu'à obtention d'une correspondance optimale entre modèle et surface. Dans le cas des modèles paramétriques, l'évolution des paramètres au cours du temps permet d'avoir une description condensée du rouvement. Dans le cas des modèles non paramétriques, le plus souvent basés sur des modèles de type masseressort et/ou éléments finis avec résolution de contraintes différentielles, le suivi temporel de points de surface permet d'obtenir l'estimation de leur mouvement. Le principal défaut de ces modèles géométriques reste la difficulté à modéliser des surfaces complexes. Dans le cas des volumes cardiaques acquis par scanner mufti-barrette, la très haute qualité de la résolution spatiale fait que la surface peut être par endroit très découpée. Ces aspérités locales seront difficilement restituées par ce type de modélisation alors qu'elles peuvent être significatives du point de vue anatomique.

Les méthodes de flot optique se basent sur trois contraintes: la conservation de la densité des points lors du mouvement, l'absence de divergence et le lissage du mouvement. Une fois ces trois contraintes respectées, la minimisation d'une fonction de pénalité permet d'obtenir un champ dense de déplacement défini sur tout le volume [Song '91, Coatrieux '94, Gorce '97]. L'application de ce type de méthodes dans le présent cas se heurte à deux difficultés: d'une part, la taille des bases de données utilisées fait que ce type de calcul s'avèrerait extrêmement coûteux en ressources machine. D'autre part, et surtout, la contrainte de conservation de la densité des points n'est pas a priori respectée. En effet, la durée d'acquisition (les données pouvant être acquises sur plusieurs cycles cardiaques) ainsi que la méthode d'acquisition en ellemême (visualisation dynamique de la concentration en produit de contraste, soumis à la diffusion) et que le principe de reconstruction des images par post-synchronisation, font qu'un point peut changer de densité au cours cle l'acquisition. Bien que cette contrainte de conservation d'intensité puisse être quelque peu relâchée dans ce type de méthodes, le calcul d'un champ de déplacement dense sur tout le volume étudié est difficilement envisageable sur un gros volume de données à traiter.

Les méthodes de suivi de forme, ou de mise en correspondance [Amini '92, Kambhamettu '03, Ruan '94, Shi '00], sont basées sur la sélection d'entités et leur suivi à partir d'une mesure de correspondance établie sur plusieurs paramètres descriptifs.

La technique de la présente invention se situe dans cette troisième famille de méthodes d'estimation de mouvement, à savoir celle des méthodes de mise en correspondance d'entités descriptives, entre les différentes images 3D prises à différents instants d'observation de l'objet et/ou avec différentes modalités d'imagerie.

En effet, les méthodes de mise en correspondance offrent des avantages par rapport aux méthodes des deux autres familles précitées.

Tout d'abord, contrairement aux méthodes basées sur les modèles déformables, les méthodes de mise en correspondance permettent, par la sélection d'entités fines, une description précise d'une surface complexe. Cependant, plus le nombre d'entités est important, plus le temps de calcul pour trouver les correspondances optimales est important.

En outre, contrairement aux méthodes de flot optique, les méthodes de mise en correspondance permettent de se limiter aux données significatives (la paroi myocardique dans le cas du coeur). Cependant, elles imposent d'avoir préalablement accès à ces données par des outils de segmentation.

De plus, les méthodes connues de mise en correspondance présentent l'inconvénient majeur de nécessiter la segmentation de toutes les images 3D étudiées.

Or, les outils de segmentation, très nombreux, ne sont pas efficaces dans tous les cas. Ils ne garantissent pas la cohérence temporelle des surfaces extraites, puisque les images 3D sont segmentées indépendamment lies unes des autres, avec un risque réel de modification d'une image à l'autre des conditions d'acquisition, des structures environnantes, de la qualité des images,... Ceci est donc source d'erreur pour l'estimation du mouvement et/ou de la déformation.

D'autre part, il est pénalisant en termes de performances de répéter la phase de segmentation pour chaque image 3D. En effet, la segmentation d'une image 3D demeure une étape de traitement difficile et fastidieuse, surtout comme c'est souvent le cas si elle est menée de manière interactive ou semi-interactive (c'est-à-dire si elle nécessite l'intervention de l'utilisateur pour fixer certains paramètres de l'algorithme par exemple).

3. Objectifs de l'invention L'invention a notamment pour objectif de pallier cet inconvénient majeur des méthodes connues de mise en correspondance.

Plus précisément, l'un des objectifs de la présente invention, dans au moins un mode de réalisation, est de fournir une technique d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D, cette technique ne nécessitant la segmentation que d'une image 3D de cet ensemble.

L'invention a également pour objectif, dans au moins un mode de réalisation, de fournir une telle technique permettant de réduire les temps de calcul.

Un autre objectif de l'invention, dans au moins un mode de réalisation, est de fournir une telle technique qui, non seulement n'est pas pénalisée par une grande complexité topologique des objets étudiés, mais au contraire tire parti de cette complexité.

4. Caractéristiques essentielles de l'invention Ces différents objectifs, ainsi que d'autres qui apparaîtront par la suite, sont atteints selon l'invention à l'aide d'un procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D comprenant au moins deux images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas. Le procédé selon l'invention comprend les étapes suivantes: a) obtention d'une surface 3D dudit au moins un objet contenu dans l'une des images 3D dudit ensemble, dite image de départ, par traitement de l'image de départ, ladite surface 3D étant représentée sous la forme d'un maillage comprenant des noeuds; b) choix de l'image de départ comme image courante; c) détermination de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante dudit ensemble, de façon à estimer un charnp de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation dudit au moins un objet entre les images courante et suivante.

L'invention tient son originalité première du fait qu'elle est basée sur une mise en correspondance d'une surface d'objet (obtenue pour l'image courante, acquise par exemple à l'instant t) et d'un volume (qui est l'image 3D suivante, acquise par exemple à un instant suivant t+1). Cette mise en correspondance surface/volume donne accès à des vecteurs de déplacement et/ou de déformation calculés pour chaque entité descriptive de la surface d'objet.

L'invention se démarque donc des méthodes connues de mise en correspondance par le fait qu'une seule image 3D (l'image 3D de départ) nécessite un traitement de type segmentation, permettant d'extraire la surface 3D de l'objet ou des objets d'intérêt. En d'autres termes, le procédé de l'invention ne nécessite la segmentation que d'une des images 3D (volumes) de l'ensemble des images 3D considérées. Le fait de n'avoir qu'une segmentation à réaliser pour une seule image 3D est ici un atout important. En effet, les méthodes connues de mise en correspondance, ainsi que bon nombre de méthodes connues basées sur un modèle déformable, nécessitent la pré-segmentation de toutes les images 3D (volumes) traitées.

Le procédé développé permet de tirer parti de la complexité topologique des objets étudiés. Dans un contexte où les progrès de l'acquisition tridimensionnelle permettent d'obtenir une résolution spatiale de plus en plus élevée, la complexité des surfaces étudiées ira toujours croissante. Contrairement aux méthodes basées sur les modèles déformables qui posent problème dans le cas de formes trop complexes, la méthode développée s'appuie sur cette complexité afin d'extraire les correspondances avec précision.

Il est à noter que la principale contrainte du procédé de l'invention est liée à l'amplitude des correspondances extraites. Celle-ci ne doit pas être trop importante au risque de compromettre la validité des résultats. Cependant, dans le cas des applications concernées, cette limite a un impact très modéré. Dans le cas de l'extraction de mouvement, les progrès des méthodes d'acquisition permettent d'obtenir des résolutions temporelles de plus en plus élevées, entraînant par là la diminution de l'amplitude des mouvements entre chaque instant considéré. Par exemple, pour l'extraction du mouvement cardiaque, les scanners MSCT permettent d'obtenir une résolution temporelle nettement suffisante pour l'utilisation du procédé selon l'invention.

Dans un premier mode de réalisation avantageux de l'invention, ledit ensemble d'images 3D comprend deux images 3D acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas.

Dans un second mode de réalisation avantageux de l'invention, ledit ensemble d'images 3D est une séquence d'images 3D comprenant au moins trois images 3D acquises à des instants successifs au sein d'une même plage temporelle ou de plusieurs plages temporelles non adjacentes, et avec une même modalité d'imagerie ou plusieurs modalités d'imagerie.

De façon avantageuse, dans le cas du second mode de réalisation précité, ladite image de départ est la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes: d) obtention d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation; e) itération des étapes c) et d) jusqu'à ce que toutes les images 3D de la séquence aient été traitées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante utilisée lors de l'itération précédente.

Ainsi, dans le cas d'une séquence d'images 3D, l'invention conduit à la détection couplée des vecteurs de déplacement et/ou déformation associés à la surface de l'objet entre les deux instants t et t+l (résultat de l'étape c) et de la surface de l'objet à l'instant t+l (résultat de l'étape d). De cette façon, dans le cas d'une séquence d'images 3D, seule l'image 3D de départ nécessite un traitement de type segmentation. Les autres images 3D de la séquence analysée ne subissent pas ce traitement, la recherche de correspondances se faisant ensuite directement entre la surface 3D de l'objet reconstruite à un instant et le volume image entier à l'instant suivant.

Selon une variante avantageuse, ladite image de départ est différente de la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes: d') obtention d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation; e') itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un premier sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la dernière image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le premier sens, utilisée lors de l'itération précédente; e") itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un second sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la première image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le second sens, utilisée lors de l'itération précédente.

Préférentiellement, les entités d'une surface 3D donnée sont des noeuds du maillage de ladite surface 3D donnée.

De façon préférentielle, les entités d'une image 3D donnée sont des voxels de ladite image 3D donnée.

Au lieu des noeuds et des voxels, des entités d'autre nature peuvent également être considérées, telle que la maille ou encore un patch de surface (ou ensemble de mailles) pour la surface à mettre en correspondance avec un ensemble de voxels sur le volume. Les voisinages de ces entités et les énergies locales associées à ces entités sont bien sûr dans ce cas à re-définir, tout en pouvant les construire sur le même type de descripteurs de base (position, courbure, ...).

Avantageusement, ladite étape c) comprend une étape de minimisation d'une énergie globale constituée d'une somme d'énergies locales, chaque énergie locale 2877465 10 permettant de comparer une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante.

Selon une caractéristique avantageuse, chaque énergie locale est composée d'une somme pondérée de: - un terme de correspondance locale E,(f,d), exprimant une mesure de cohérence d'un appariement entre une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante; - un terme de régularisation E2(f), exprimant des contraintes de lissage spatio- temporel du mouvement et/ou de la déformation estimée, lesdites contraintes étant telles que des vecteurs de déplacement associés à des entités voisines de la surface 3D doivent être proches en amplitude et/ou en direction.

Avantageusement, ledit terme de correspondance locale E,(f,d) est fonction d'un terme exprimant une correspondance topologique entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante.

De façon avantageuse, ledit terme de correspondance locale E,(f,d) est composé d'une somme pondérée dudit terme exprimant une correspondance topologique et d'au moins un des termes appartenant au groupe comprenant: - un terme exprimant une probabilité que l'entité de l'image 3D suivante appartienne à un élément de contour, dudit au moins un objet contenu dans ladite image 3D suivante; - un terme exprimant une distance entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante; - un terme lié aux densités de niveaux de gris.

Il est à noter que l'élément de contour peut être soit préalablement détecté, soit recherché de manière ponctuelle (si besoin) au moment du calcul du terme d'énergie exprimant une probabilité et impliquant la connaissance de cet élément de contour.

Dans un mode de réalisation particulier de l'invention, ladite étape de minimisation d'une énergie globale est basée sur une formulation statistique d'un processus d'estimation selon un modèle de Markov, telle que le mouvement et/ou la déformation à estimer est considéré comme une réalisation f = {f, , i =1,...,NS} d'un champ aléatoire F = {F, , i =1,... , NS} regroupant toutes les réalisations possibles, les NS sites du champ aléatoire F étant donnés par toutes les entités (par ex. noeuds) de la surface 3D, obtenue pour l'image 3D courante, représentatives du maillage de ladite surface 3D, des labels desdits sites, aussi appelés étiquettes recherchées, étant les entités (par ex. voxels) de l'image 3D suivante jugées correspondantes desdites entités de la surface 3D.

Cette formalisation du problème d'estimation dans le cadre statistique des champs markoviens 3D permet d'optimiser les calculs.

De façon préférentielle, la surface 3D obtenue pour l'image 3D courante, lors de l'étape a) ou d), est considérée comme une surface 3D de résolution maximale. En outre, 10 l'image 3D suivante comprise dans ledit ensemble d'images 3D est considérée comme une image 3D de résolution maximale. Enfin, ladite étape c) est effectuée dans un mode multirésolution et comprend les étapes successives suivantes: cl) obtention d'une part d'au moins une surface 3D de résolution inférieure pour l'image 3D courante, à partir de la surface 3D de résolution maximale pour l'image 3D courante, et d'autre part d'au moins une image 3D suivante de résolution inférieure, à partir de l'image 3D suivante de résolution maximale; c2) choix de la plus basse des résolutions comme résolution courante, et initialisation des correspondances à la résolution courante; c3) pour la résolution courante, détermination de correspondances optimales entre des entités de la surface 3D possédant ladite résolution courante pour l'image 3D courante, et des entités de l'image 3D suivante possédant également ladite résolution courante; c4) si la résolution courante n'est pas égale à la résolution maximale: * propagation de résultats de l'exécution de l'étape c3), de la résolution courante à une résolution supérieure suivante, de façon à permettre une réinitialisation des correspondances; * retour à l'étape c3), en prenant comme nouvelle résolution courante la résolution supérieure suivante; sinon arrêt de l'étape c).

Ainsi, l'intégration d'un processus de multi-résolution permet d'affiner l'estimation de manière guidée et progressive. On tient compte et on exploite les 20 25 30 différents niveaux de détails topologiques des surfaces complexes, et on réduit le temps de calcul.

Avantageusement, ladite étape a) comprend les étapes suivantes: segmentation de l'image de départ, de façon à extraire des contours et points intérieurs audit au moins un objet contenu dans l'image de départ; - reconstruction d'une surface 3D, représentée sous la forme d'un maillage comprenant des noeuds, à partir des contours et points intérieurs obtenus lors de l'étape de segmentation.

De façon avantageuse, dans l'étape c2), l'initialisation des correspondances à la résolution la plus basse est réalisée par des mouvements nuls, chaque entité de la surface 3D obtenue pour l'image 3D courante étant mise en correspondance avec l'entité la plus proche dans l'image 3D suivante.

Avantageusement, dans l'étape c4), la réinitialisation des correspondances est obtenue par interpolation, sur tout le maillage de la surface 3D possédant ladite résolution supérieure suivante, des correspondances optimales déterminées lors de l'étape c3).

De façon avantageuse, ladite étape a) comprend en outre l'étape suivante: pré-traitement de la surface 3D reconstruite, permettant de modifier, si nécessaire, ledit maillage de façon que la longueur de chacune des arêtes constituant une maille soit égale à une des distances entre voxels des images 3D.

Selon une caractéristique avantageuse, les modalités d'imagerie appartiennent au groupe comprenant: l'imagerie scanner; l'imagerie par résonance magnétique; l'imagerie par ultrasons; l'imagerie nucléaire. Cette liste n'est pas exhaustive.

Le procédé selon l'invention s'applique avantageusement dans un domaine d'imagerie appartenant au groupe comprenant: l'imagerie médicale, l'imagerie marine, l'imagerie satellitaire et l'imagerie des matériaux. Cette liste n'est pas exhaustive.

L'invention concerne également un produit programme d'ordinateur comprenant des instructions de code de programme pour l'exécution des étapes du procédé selon l'invention, lorsque ledit programme est exécutésur un ordinateur.

L'invention concerne aussi un moyen de stockage, éventuellement totalement ou partiellement amovible, lisible par un ordinateur, stockant un jeu d'instructions exécutables par ledit ordinateur pour mettre en oeuvre le procédé selon l'invention.

L'invention concerne encore un procédé de recalage tridimensionnel d'au moins un objet contenu dans une image 3D, comprenant les étapes suivantes: estimation (41) du mouvement et/ou de la déformation subi(s), entre une première et une seconde images 3D acquises à deux instants différents et/ou à partir de deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas, par au moins un objet contenu dans lesdites première et seconde images 3D; application (42) audit au moins un objet contenu dans la seconde image d'une fonction corrigeant le mouvement et/ou la déformation estimé lors de l'étape d'estimation. En outre, ladite étape d'estimation (41) est mise en oeuvre selon le procédé d'estimation de l'invention.

5. Liste des figures D'autres caractéristiques et avantages de l'invention apparaîtront à la lecture de la description suivante d'un mode de réalisation préférentiel de l'invention, donné à titre d'exemple indicatif et non limitatif, et des dessins annexés, dans lesquels: la figure 1 présente un organigramme d'un mode de réalisation particulier du procédé d'estimation de mouvement et/ou de déformation selon l'invention; la figure 2 détaille un exemple de réalisation de l'étape référencée 11 sur la figure 1; la figure 3 détaille un exemple de réalisation de l'étape référencée 13 sur la figure 1; et la figure 4 présente un organigramme d'un mode de réalisation particulier du procédé de recalage tridimensionnel selon l'invention.

6. Description d'un mode de réalisation de l'invention On présente cidessous, en relation avec les figures 1 à 3, un mode de réalisation particulier du procédé selon l'invention d'estimation du mouvement et/ou de la déformation d'un objet contenu dans chaque image 3D d'un ensemble d'images 3D.

Dans ce mode de réalisation particulier, l'ensemble d'images 3D considéré est une séquence temporelle d'images 3D acquises au sein d'une même plage temporelle et avec une même modalité d'imagerie (par exemple en imagerie scanner multi-barrette (MSCT)). Il est clair cependant que l'invention peut également être mise en oeuvre avec des images 3D acquises dans des plages temporelles non adjacentes et/ou avec plusieurs modalités d'imagerie. Par ailleurs, l'invention peut aussi s'appliquer si l'ensemble d'images 3D ne comprend que deux images (voir ci-après la description du procédé de recalage tridimensionnel, en relation avec la figure 4).

Dans ce mode de réalisation particulier, à partir d'images 3D acquises à différents instants (séquence temporelle), le procédé selon l'invention permet la détection spatio/temporelle et le suivi d'un mouvement et/ou d'une déformation d'objet, par un principe de mise en correspondance d'une surface représentée sous forme de maillage et d'une image 3D (aussi appelée volume 3D). Les données de départ sont constituées d'une séquence d'images 3D représentant un objet (ou un ensemble d'objets) dans son (ou leur) environnement, cet objet étant soumis, entre chaque image 3D, à un mouvement et/ou à une déformation de type non-rigide et d'amplitudes modérées.

On rappelle que le principe général de l'invention consiste à extraire les correspondances locales entre le maillage initial (surface 3D) décrivant l'objet considéré et l'image 3D suivante, correspondant à la première déformation de cet objet. L'étape de segmentation, qui est cruciale et difficile en analyse d'images, n'est donc ici nécessaire que pour une (préférentiellement la première) image 3D de la séquence analysée. Une fois extraite, ces correspondances permettent alors de modifier le maillage surfacique initial afin d'obtenir un maillage décrivant de façon satisfaisante l'état de l'objet détecté dans le volume considéré après cette première déformation. En d'autres termes, la surface 3D de l'objet contenu dans l'image 3D suivante est obtenue, non pas en effectuant une segmentation de cette image suivante, mais grâce à l'application de vecteurs de déplacement et/ou de déformation, résultant de l'extraction des correspondances, à la surface 3D de l'objet contenu dans l'image 3D courante. Ce processus est appliqué de manière itérative afin d'extraire les correspondances entre les différentes images 3D de la séquence, de façon à décrire tous les mouvements et/ou déformations de l'objet. Il permet d'extraire, de façon simultanée, la position de l'objet déformé dans le volume analysé (image 3D suivante) d'une part, l'amplitude et la direction des déplacements et/ou déformations de l'objet d'autre part, estimés de façon locale.

Comme illustré sur la figure 1, dans le mode de réalisation particulier précité, le procédé selon l'invention d'estimation de mouvement et/ou de déformation comprend les étapes suivantes: - dans une première étape 11, on choisit l'image de départ (la première image 3D de la séquence dans cet exemple) comme image courante; dans une seconde étape 12, cette image courante (image de départ) est soumise à un traitement permettant d'obtenir une surface 3D de l'objet qu'elle contient. Cette surface 3D est représentée sous la forme d'un maillage comprenant des noeuds; dans une troisième étape 13, on détermine des correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante au sein de la séquence, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation de l'objet entre les images courante et suivante; dans une quatrième étape 14, on obtient une surface 3D de l'objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation; dans une cinquième étape 15, on détermine si l'image suivante précitée est la dernière image de la séquence, c'est-à- dire si toutes les images 3D de la séquence ont été traitées; dans l'affirmative, on passe à une sixième étape de fin 16; dans la négative, on effectue une septième étape 17 consistant à prendre l'image suivante précitée comme nouvelle image courante (avec la surface 3D issue de la quatrième étape 14) et l'image suivante de la suivante précitée comme nouvelle image suivante, avant de revenir à la troisième étape 13.

Dans une variante de réalisation, la deuxième étape 12 est effectuée avec une image de départ qui n'est pas la première image 3D de la séquence. Il convient dans ce 30 cas de traiter: 10 15 20 25 - toutes les images situées avant l'image de départ, au sein de la séquence (l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon un premier sens de parcours de la séquence, utilisée lors de l'itération précédente) ; et - toutes les images situées après l'image de départ, au sein de la séquence (l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon un second sens de parcours de la séquence, utilisée lors de l'itération précédente).

Comme illustré sur la figure 2, la deuxième étape 12 comprend par exemple les étapes suivantes: segmentation (121) de l'image de départ, de façon à obtenir les contours et points intérieurs de l'objet contenu dans l'image de départ. Cette étape de segmentation est par exemple réalisée,en combinant une approche de croissance de régions et des informations contour extraites par un opérateur gradient (voir [Guillaume'03 I). Elle permet de fournir comme résultat un ensemble de points 3D supposés appartenir à l'objet étudié ; reconstruction (122) d'une surface 3D, représentée sous la forme d'un maillage comprenant des noeuds, à partir des points identifiés comme appartenant à l'objet lors de l'étape de segmentation. Cette étape de reconstruction est par exemple réalisée en soumettant l'ensemble des points 3D retenus à l'algorithme classique des Marching Cubes (création de topologies locales, représentées sous forme de noeuds et de liens entre ces noeuds, grâce à l'application de masques caractéristiques) ; pré-traitement (123) de la surface 3D reconstruite, permettant de modifier, si nécessaire, le maillage de façon que la dimension des mailles de la surface 3D soit régulière et proche de la dimension de chacun des voxels des images 3D. L'étape (123) de prétraitement permet de s'assurer que les arêtes constituant une maille (une maille étant délimitée par trois noeuds) sont de longueur égale aux distances inter-voxels (elles prennent leur valeurs dans l'ensemble {1, \/2, 1/3}. Ce pré-traitement est assuré par un passage en valeurs entières des coordonnées de tous les noeuds (approximation à l'entier le plus proche), une fusion des 20 25 30 noeuds qui peuvent être rendus confondus par l'étape précédente, une retriangularisation du maillage (décomposition ou fusion des mailles ne comportant pas trois noeuds) Il est à noter que cette contrainte sur la longueur des arêtes des mailles devant être égale à une distance inter-voxels peut être relâchée tout en restant dans le cadre de la présente invention, mais implique de re-définir le voisinage des entités de base choisies et de re-spécifier également la formulation des énergies.

On décrit maintenant plus en détail la troisième étape 13 de détermination de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités de l'image 3D suivante.

On rappelle que dans la présente description, on utilise indifféremment les termes image 3D et volume Afin d'établir des correspondances optimales entre une surface et un volume complet, la présente invention propose une méthode de mise en correspondance. Ce type de méthode nécessite de considérer trois étapes: - la sélection des entités à mettre en correspondance, - la définition d'une énergie locale permettant de comparer deux entités, - la minimisation d'une énergie globale constituée de la somme des énergies locales.

Une des originalités de la méthode de mise en correspondance selon l'invention est qu'elle établit des correspondances entre des entités n'appartenant pas au même domaine. En effet, la mise en correspondance est réalisée entre des entités d'un maillage pré-extrait (défini par des noeuds et des cellules (ou mailles)) et des entités d'un volume complet (composés de voxels, volumes tridimensionnels élémentaires). II est à noter que ces deux ensembles sont très dissemblables: le maillage apporte plus de connaissance (notamment sur l'emplacement de la frontière à étudier, potentiellement difficile à segmenter) que le volume, "données brutes" issues de l'acquisition. Les entités utilisées pour la description du maillage sont les noeuds qui le constituent. Les entités à faire correspondre avec ces noeuds sont les voxels clu volume.

L'énergie locale doit permettre de comparer deux entités, chacune appartenant à son ensemble respectif (noeuds d'un côté, voxels de l'autre) tout en tenant compte des contraintes établies sur les correspondances du voisinage de chaque entité. Elle est par exemple composée de la somme pondérée de deux termes principaux (équation 1) : d'une part un terme de correspondance locale ou terme d'attache aux données E, (f;, d) (exprimant une distance entre estimation du mouvement (f) associée à une entité i (i=1,...,NE) (NE étant le nombre d'entités considérées) et observation (d) ),, d'autre part un terme de régularisation E2 (f,) (exprimant des contraintes de lissage spatio-temporel des déformations).

E(f,,d)=a,El(f,,d)+a2E2(f,) (1) Le terme d'attache aux données El (f,, d) est directement lié aux données du problème, c'est-à-dire à l'estimation du déplacement d'un noeud (f,.) entre les instants t et t+l, conditionnellement aux observations (volumes images à t et t+l) (d).

Il est donné par une mesure de cohérence d'un appariement et a pour but de quantifier: - la probabilité que le voxel considéré à t+l appartienne à la frontière étudiée, - la correspondance topologique entre le noeud à t et le voxel à t+l considérés, - la distance entre le noeud et le voxel.

La probabilité que le voxel considéré appartienne à la frontière étudiée est évaluée grâce à l'application d'un détecteur de contour (qui peut être par exemple le filtrage de Canny) sur le volume image analysé. Ainsi, un voxel détecté comme appartenant au contour est privilégié à un voxel détecté comme n'appartenant pas à une frontière. D'autres types de détecteurs de contour peuvent être retenus dans cette phase, appliqués globalement sur le volume ou localement.

Un terme de correspondance topologique est défini et permet d'effectuer une évaluation de la compatibilité de configuration entre le noeud et le voxel considérés et leur voisinage respectif (le voisinage d'un nceud du maillage est défini par les noeuds ayant une maille en commun avec ce noeud central, et le voisinage d'un voxel, par les voxels vérifiant une connexité première avec ce dernier). Cette correspondance topologique peut se baser sur différents critères descriptifs de la topologie locale. Ainsi, dans le cas de la surface, les courbures moyenne et gaussienne sont de bonnes caractéristiques descriptives. Afin d'implémenter le calcul de ces courbures dans le cas du volume 3D et d'obtenir des descripteurs qui soient comparables à ceux calculés sur la surface, la recherche de la topologie locale du voxel étudié est effectuée. Celle-ci est basée sur l'information de contour fournie par le détecteur de contour et correspond à la recherche de la surface d'iso-intensité (points appartenant au contour) au voisinage du voxel. Une fois cette surface locale trouvée (modélisée par un maillage), le calcul de courbure peut s'effectuer de la même façon que pour la surface, mais est appliqué de manière très locale.

Une autre technique de description de la topologie locale proposée dans cette invention consiste à transposer la topologie du voisinage du noeud étudié au niveau du voxel candidat. Cette topologie est décrite en termes de différences de coordonnées des voisins par rapport aux coordonnées du noeud étudié. Les valeurs des voxels correspondants à cette topologie (dans le volume filtré avec le détecteur de contours) sont sommées. Ainsi, plus cette somme est élevée, plus les topologies locales des deux entités sont considérées comme proches. A ce niveau, le voisinage défini peut être étendu à un nombre plus important de mailles.

Enfin, afin de privilégier les correspondances représentant des mouvements faibles, la distance euclidienne entre noeud considéré et voxel candidat est calculée. Les correspondances à faible distance sont privilégiées.

Le terme d'attache aux données El (f,, d) peut donc être exprimé pour chaque couple Noeud (N, , i =1, ..., NS, NS: nombre de noeuds sur la surface à t) Voxel (Vk, k =1, ...,NV, NV: nombre de voxels du volume à t+1), à partir de trois termes représentatifs des potentiels d'énergie précédemment décrits. l El (/; > d) (N, , Vk) = a,. Ec,,ntour (N, , Vk)+ a, E,01,01 (N, , Vk)+ a d Ed,, (N, , Vk) A ces trois termes composant la fonction d'énergie locale, d'autres termes peuvent également être ajoutés, introduisant par exemple une contrainte sur l'intensité de niveaux de gris de l'image. En effet, un terme d'énergie pourrait être défini en fonction de l'information densitométrique des images à partir du moment où une correspondance entre intensités de niveau de gris peut être établie entre les deux volumes considérés (let volume à partir duquel est défini la surface 3D et 2ème volume dans lequel est recherché la correspondance).

Le terme de régularisation E2 (f) permet d'éviter l'obtention de correspondances aberrantes (par exemple les mouvements de deux noeuds voisins qui se croisent). Elle est basée sur le fait que les correspondances de noeuds voisins doivent engendrer des mouvements proches, que ce soit en amplitude ou en direction. Différentes implémentations de cette énergie peuvent être envisagées (sachant que l'énergie finalement retenue peut être également une combinaison de différentes formulations).

Par exemple, une première formulation peut correspondre à l'écart entre le déplacement évalué et la moyenne des déplacements des voisins du noeud considéré. D'autres formulations peuvent également convenir. Par exemple, une deuxième formulation est proposée dans le cadre de la formulation Markovienne (voir discussion détaillée ci-après), permettant de maintenir l'homogénéité des appariements.

Afin d'estimer globalement les correspondances à partir des énergies locales formulées précédemment, une formulation markovienne est par exemple utilisée.

Le modèle de Markov permet de formuler un problème d'estimation conditionnellement à des observations, sous la forme d'un modèle global décrit par des spécifications locales (introduisant la définition de voisinage).

Dans sa définition, le mouvement 3D recherché est considéré comme une réalisation f = {t, i =1,..., NS} d'un champ aléatoire F = {F,. , i =1,... , NS} regroupant toutes les réalisations possibles. Les sites du champ considéré sont donnés par tous les noeuds représentatifs du maillage de la surface. Les labels de ces sites ou étiquettes recherchées (estimations f) sont les voxels jugés correspondants de ces noeuds. Le voisinage d'un site (ou d'un noeud) est composé des noeuds ayant une maille en commun avec le noeud étudié.

La réalisation la plus probable du champ des étiquettes est donnée par la maximisation d'une probabilité P (f) qui s'écrit: P (f Z exp E T U (f (3) où T est une valeur constante appelée température et U (f) est une fonction d'énergie globale, Z = Z exp T U (f, Q étant l'ensemble des étiquettes possibles.

La réalisation la plus probable revient alors minimiser la fonction d'énergie globale U (f).

La fonction U (f) ou plus généralement U (f, d) s'exprime par l'équation suivante: U(f,d)=a,Ud(f,,d)+a2U,(f) (4) Où Ud (f, d) modélise l'erreur de l'estimation par rapport aux observations et le terme U, (f) modélise l'énergie interne du champ correspondant au mouvement à restituer. Ce deuxième terme joue un rôle de stabilisation dans le problème de régularisation posé.

Le premier terme s'exprime comme une erreur quadratique de manière que la minimisation de cette composante d'énergie favorise les configurations estimées proches 10 des observations. Il est calculé à partir de la mesure d'appariement locale (ou distance) Une formulation de E2 (f,) est également proposée par l'équation suivante: E2 (f)_ V (f, fi), où r); est le voisinage de première connexité associé au site i (7) où V f, f) est le potentiel d'interaction associé à une clique du voisinage (une 20 clique étant représentée ici par deux sites) avec: V (,, )= l -( f . fi), (8) où 8() est une fonction delta Kronecker qui prend la valeur 1 si les étiquettes sont proches et la valeur 0 sinon.

Plusieurs formulations de E2 (f) sont ici possibles en fonction des cliques 25 définies sur le voisinage et de leur potentiel d'interaction associé.

On rappelle qu'une clique définie sur un voisinage est constituée d'un sous-ensemble d'éléments du voisinage, contenant de 1 à N éléments (le site compris), N étant le nombre d'éléments du voisinage. Ces éléments constituant la clique doivent vérifier une connexité (point ou arête de maille en commun dans le cadre de cette invention).

développée précédemment, sur tous les sites du champ: Ud (f, d)_ (E, (f,, d)f, L étant l'ensemble des sites du champ (noeuds du maillage) (5) Le second terme est exprimé de la façon suivante: U,(f)= E2(f) (6) La minimisation de l'énergie globale U (f,d) permet de trouver les correspondances optimales et donc l'estimation de déformation recherchée pour tout l'objet. Les différentes réalisations f peuvent être obtenues selon plusieurs algorithmes d'optimisation. Une méthode pour résoudre ce problème, dans le contexte des champs de Markov, est la recherche aléatoire combinée avec une technique de recuit simulé. L'algorithme de relaxation stochastique de Metropolis est par exemple choisi, combiné avec un recuit simulé, permettant de mieux assurer la convergence vers un minimum global, comparé à un algorithme de relaxation déterministe. D'autres algorithmes d'optimisation sont ici possibles.

Un exemple d'algorithme d'optimisation s'écrit comme suit: Initialisation des correspondances T:=To Tant que {Critère d'arrêt non vérifié} Nbre_essais_à_T:= 0 Nbre_succès:= 0; Tant que {Nbre_essais_à_T < Nbre_essais_à_T _max} ou {Nbre_succès_à_T < Nbre_ succès_à_T _max} Sélection du n ud étudié (qui peut être choisi aléatoirement) Sélection du voxel candidat Nbre succès à T: = Nbre succès à T + 1 Calcul de AU(f, d) Si {AU(f,d)≤0} Validation de la correspondance candidate Sinon Proba = (exp(- AU(f, d) IT)) Si {Proba > random(0,1)} Validation de la correspondance candidate Nbre_succès:= Nbre_succès+l; Fin 20 Fin Fin T:=T-pas_T (diminution de T) Fin La valeur initiale des paramètres impliqués dans cette procédure d'optimisation (température To, Nbre_essais_à_T _max (nombre d'essais maximum à une température donnée) et Nbre_succès_max (nombre maximum de succès à AU(f,d) positif) a été déterminée empiriquement.

La démarche développée dans ce schéma est la suivante. C'est un procédé itératif dans lequel une température donnée initialement va décroître jusqu'à la rencontre d'un critère d'arrêt, pour lequel la solution de correspondance globale est proposée. La loi de décroissance de la température à chaque itération peut prendre plusieurs formes, par exemple résulter de la décrémentation de la température d'un pas.

Pour chaque itération (à une température donnée), plusieurs noeuds sont choisis (soit de manière aléatoire) et traités, ceci jusqu'à ce qu'un nombre suffisant de noeuds ait été considéré (Nbre_essais_à_T) ou jusqu'à ce qu'un nombre donné de succès (Nbre_succès_à_T) ait été atteint. Pour chacun des noeuds sélectionnés, la possibilité de modifier la correspondance en cours est étudiée et acceptée en fonction de la variation d'énergie globale obtenue par ce changement. Un succès correspond à un changement de correspondance pour un noeud, qui a été accepté avec une variation d'énergie globale positive.

La valeur initiale des paramètres impliqués dans cette procédure d'optimisation a été déterminée empiriquement. Ces paramètres sont: To: la température de départ, Pas_T: le pas de décroissance de la température, Nbre_essais_à_T_max: le nombre maximum d'essais à une température donnée, Nbre_succès_à_T _max: le nombre maximum de succès, à une température donnée, avec AU(f,d) (variation de l'énergie globale) positive (c'est-àdire entraînant une augmentation de l'énergie globale).

Le critère d'arrêt peut être défini de différentes façons: soit la température T doit atteindre un seuil (qui peut être la température 0), soit le nombre de modifications de correspondance à une température donnée doit être inférieur à un seuil.

La sélection du noeud étudié peut-être soit complètement aléatoire, soit guidée par des critères d'exhaustivité (par exemple empêcher de sélectionner deux fois le même noeud pour une température donnée).

La sélection du voxel candidat s'effectue en choisissant (soit aléatoirement, soit suivant des critères d'exhausitivité ou de non répétition) un voxel voisin du voxel correspondant courant. Ce voisinage peut être défini comme étant un voisinage de premier ordre (les deux voxels ont alors une face ou une arrête en commun), ou d'ordre supérieur (par exemple pour l'ordre deux, les deux voxels ont une face ou une arrête ou un sommet en commun).Optionnellement, afin de tenir compte et d'exploiter les différents niveaux de détails topologiques des surfaces complexes (surface 3D de l'objet), de converger plus vite vers la solution optimale et dans le but de réduire le temps de calcul, la troisième étape 13 de détermination de correspondances optimales est appliquée dans un mode multi-résolution.

Dans ce cas, la surface 3D obtenue pour l'image 3D courante, lors de la deuxième étape 12 ou la quatrième étape 14, est considérée comme une surface 3D de résolution maximale. En outre, l'image 3D suivante comprise dans l'ensemble d'images 3D est considérée comme une image 3D de résolution maximale.

Le mode multi-résolution nécessite de pouvoir passer d'une résolution à une résolution inférieure, que ce soit pour les volumes de données originaux (c'est-à-dire les images 3D originales de la séquence source) ou pour les maillages 3D (c'est-à-dire les surfaces 3D).

Pour les volumes originaux, le passage d'une représentation de résolution donnée à une représentation de résolution inférieure est par exemple réalisé par un lissage (qui peut être un filtrage gaussien) suivi d'un sous-échantillonnage appliqués au volume image.

Pour les maillages 3D, le passage d'une représentation de résolution donnée à une représentation de résolution inférieure est par exemple réalisé par un lissage (qui peut être un filtrage gaussien) suivi d'un sous-échantillonnage des coordonnées des noeuds, puis une opération de fusion des noeuds confondus et de re-triangularisation du maillage (afin que toutes les mailles soient constituées de trois noeuds).

Pour traiter le premier intervalle temporel, t - t+l, le changement de résolution à t est appliqué sur le maillage de surface, tandis que le changement de résolution à t+l est effectué sur le volume d'origine (image de niveaux de gris). Pour traiter chaque intervalle temporel suivant, t+k - t+k+l, le changement de résolution à t+k est appliqué sur le maillage de surface (issu de l'estimation de correspondance réalisée aux instants précédents et ayant subi le pré-traitement déjà cité (régularisation et retriangularisation)), tandis que le changement de résolution à t+k+1 est effectué sur le volume d'origine.

Dans ce cas, et comme illustré sur la figure 3, la troisième étape 13 précitée comprend les étapes successives suivantes: - dans une étape référencée 131, obtention d'une part d'au moins une surface 3D de résolution inférieure (pour l'image 3D courante, et d'autre part d'au moins une image 3D suivante de résolution inférieure (voir la discussion détaillée ci-dessus donnant un exemple de passage d'une résolution donnée à une résolution inférieure, que ce soit pour les images 3D originales de la séquence source ou pour les surfaces 3D)... On travaille par exemple avec trois niveaux de résolution, y inclus le niveau de résolution maximal; - dans une étape référencée 132, choix de la plus basse des résolutions comme résolution courante, et initialisation des correspondances à la résolution courante; - dans une étape référencée 133, pour la résolution courante, détermination de correspondances optimales entre des noeuds de la surface 3D (possédant cette résolution courante) contenue dans l'image 3D courante, et des voxels de l'image 3D suivante (possédant cette même résolution courante) de la séquence. Cette étape 133 est par exemple réalisée selon l'exemple d'algorithme d'optimisation décrit ci-dessus (algorithme de Metropolis combiné au recuit simulé), dans le cadre d'une formulation markovienne du problème d'estimation; - dans une étape référencée 134, on détermine si la résolution courante est égaleà la résolution maximale; - dans l'affirmative, la troisième étape 13 est terminée (et on passe à l'étape 14 suivante, voir la figure 1) ; - dans la négative, on retourne à l'étape référencée 133, après exécution d'une étape référencée 135 consistant à : * propager les résultats de l'exécution de l'étape référencée 133 (c'est- à-dire les champs de vecteurs de déplacement et/ou de déformation estimés pour les noeuds de la surface 3D), de la résolution courante à une résolution supérieure suivante, de façon à permettre une réinitialisation des correspondances; * prendre comme nouvelle résolution courante la résolution supérieure suivante.

Dans l'étape référencée 132, l'initialisation des correspondances à la résolution la plus basse est par exemple réalisée par des mouvements nuls, chaque noeud (de la surface 3D obtenue pour l'image courante) étant mis en correspondance avec le voxel le plus proche (au sein de l'image 3D suivante).

Dans l'étape référencée 135, la propagation des résultats, d'une résolution donnée à la résolution supérieure, c'est-à-dire l[a réinitialisation des correspondances est par exemple assurée par interpolation, sur tout le maillage de la surface 3D possédant cette résolution supérieure, des correspondances issues de la résolution inférieure (c'est- à-dire déterminées au cours de l'étape référencée 133).

Pour ce faire, pour chaque noeud à la résolution supérieure, le ou les noeuds correspondants à la résolution inférieure (suivant que le souséchantillonnage donne des coordonnées entières ou non) sont sélectionnés et leur déplacement moyenné. Ce déplacement moyen est transposé à la résolution supérieure (par multiplication par le facteur de souséchantillonnage) et appliqué après passage en valeurs entières (un voxel correspondant est estimé).Quand la troisième étape 13 est terminée, on dispose ainsi d'un champ de vecteurs de déplacement et/ou de déformation estimés en chaque point de la surface 3D obtenue pour l'image courante. Ce champ de vecteurs représente le mouvement et/ou la déformation de l'objet entre l'instant t d'acquisition de l'image 3D courante et l'instant t+l d'acquisition de l'image 3D suivante. Ces vecteurs, localisés dans l'espace 3D, sont caractérisés par une orientation et une amplitude.

Comme déjà indiqué ci-dessus, ces vecteurs estimés permettent également d'extraire les points de contours de l'objet à l'instant t+1 et donc de calculer la nouvelle surface correspondant à l'instant t+l (étape référencée 14 sur la figure 1). Puis, ce processus de mise en correspondance est appliqué de la même manière sur tous les instants successifs de la séquence (entre t+l et t+2, puis t+2 et t+3 etc), à partir de la nouvelle surface calculée et du volume image de l'instant suivant. Il fournit ainsi un champ de vecteurs de déplacement à chaque instant de la séquence, associé aux noeuds des mailles des surfaces successivement estimées.

On présente maintenant, en relation avec la figure 4, un mode de réalisation particulier du procédé de recalage tridimensionnel selon l'invention.

De façon classique, un tel procédé de recalage comprend: - une étape (41) d'estimation du mouvement et/ou de la déformation subi(s) , entre une première et une seconde images 3I) acquises à deux instants différents et/ou à partir de deux modalités d'imagerie différentes et/ou à partir d'une modalité d'imagerie et un atlas (modèle préexistant), par (au moins) un objet contenu dans ces première et seconde images 3D; - une étape (42) d'application à l'objet contenu dans la seconde image d'une fonction corrigeant le mouvement et/ou la déformation estimé lors de l'étape d'estimation.

L'originalité réside ici dans le fait que l'étape (41) comprend elle-même les étapes référencées 11 à 13 du procédé d'estimation discuté ci-dessus en relation avec les figures 1 à 3. Ces étapes 11 à 13 ne sont donc pas décrites à nouveau. On notera simplement que l'ensemble d'images 3D considéré ne comprend que deux images 3D (il ne s'agit donc pas d'une séquence temporelle comprenant plus de deux images).

Bien que l'invention ait été décrite ci-dessus en relation avec un nombre limité de modes de réalisation, l'homme du métier, à la lecture de la présente description, comprendra que d'autres modes de réalisation peuvent être imaginés sans sortir du cadre de la présente invention. En conséquence, la portée de l'invention n'est limitée que par les revendications ci-jointes.

Annexe: bibliographie [Amini '92] Amini, A. A. and Duncan, J. S. (1992). "Bending and stretching models for LV wall motion analysis from curves and surfaces." Image and Vision Computing 10(6): 418-430.

[Bardinet '96] Bardinet, E., Cohen, L. D. et al. (1996). "Tracking and motion analysis of the left ventricle with deformable superquadrics." Medical Image Analysis 1(2): 129-49.

[Benayoun '95] Benayoun S., Nastar C., Ayache N. (1995), "Dense Non-Rigid Motion Estimation in Sequences of 3D images Using Differential Constraints", Computer Vision, Virtual Reality and Robotics in Medicine, pp. 309-318.

[Chen '94] Chen, C. W., Huang, T. S. et al. (1994). "Modeling, analysis and visualization of left ventricle shape and motion by hierarchical decomposition." IEEE Transactions on Pattern Analysis and Machine Intelligence 16(4): 342-356.

[Choi '01] Choi, S.-M. and Kim, M.-H. (2001). "Motion visualization of human left ventricle with a time-varying deformable model for cardiac diagnosis." The journal of visualization and computer animation 12: 55-66.

[Coatrieux '94] Coatrieux, J.L., Garreau M., Ruan S., Mao F. (1994). "On motion analysis in medical imaging". Innov. Technol. Biol. Med., 3: 253267.

[Gorce '97] Gorce J.M., Friboulet D., Magnin I.E. (1997), "Estimation of three- dimensional cardiac velocity fields: assessment of a differential method and application to three-dimensional et data", Medical Image Analysis, 1(3) : 245-61.[Guillaume '03] Guillaume H., Garreau M., Boldak C., Larralde A. (2003). "Segmentation de cavités cardiaques en imagerie scanner multbarettes", 12ème Forum des Jeunes Chercheurs en Génie biologique et Médical, Nantes, 92-93.

[Kambhamettu '031 Kambhamettu, C., Goldgof, D., et al. (2003). "3D nonrigid motion analysis under small deformations." Image and Vision Computing 21(3): 229-245.

[Lin '00] Lin, W. and Robb, R. A. (2000). "Visualization of cardiac dynamics using physics-based deformable mode"l. SPIE Medical Imaging 2000: Image Display and Visualization.

[Ruan '94] Ruan S., Bruno A., Coatrieux J.L. (1994). "Three-dimensional motion and reconstruction of coronary arteries frorn biplane cineangiography", Image and Vision Computing, 12,10: 683-689.

[Shi '00] Shi, P., Sinusas, A. J., et al. (2000). "Point-tracked quantitative analysis 5 of left ventricular motion from 3D image sequences." IEEE Transactions on Medical Imaging 19(1): 36-50.

[Song '91] Song, S. M. and Leahy, R. M. (1991). "Computation of 3-D velocity fields from 3-D cine et images of the human heart." IEEE Transactions on Medical Imaging 10(3): 295-306.

Claims (22)

REVENDICATIONS

1. Procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D comprenant au moins deux images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas, caractérisé en ce qu'il comprend les étapes suivantes: a) obtention (12) d'une surface 3D dudit au moins un objet contenu dans l'une des images 3D dudit ensemble, dite image de départ, par traitement de l'image de départ, ladite surface 3D étant représentée sous la forme d'un maillage comprenant des noeuds; b) choix (11) de l'image de départ comme image courante; c) détermination (13) de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante dudit ensemble, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation dudit au moins un objet entre les images courante et suivante.

2. Procédé selon la revendication 1, caractérisé en ce que ledit ensemble d'images 3D comprend deux images 3D acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas.

3. Procédé selon la revendication 1, caractérisé en ce que ledit ensemble d'images 3D est une séquence d'images 3D comprenant au moins trois images 3D acquises à des instants successifs au sein d'une même plage temporelle ou de plusieurs plages temporelles non adjacentes, et avec une même modalité d'imagerie ou plusieurs modalités d'imagerie.

4. Procédé selon la revendication 3, caractérisé en ce que ladite image de départ est la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes: d) obtention (14) d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation; 15 e) itération des étapes c) et d) jusqu'à ce que toutes les images 3D de la séquence aient été traitées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante utilisée lors de l'itération précédente.

5. Procédé selon la revendication 3, caractérisé en ce que ladite image de départ est différente de la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes: d') obtention d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation; e') itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un premier sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la dernière image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le premier sens, utilisée lors de l'itération précédente; e") itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un second sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la première image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le second sens, utilisée lors de l'itération précédente.

6. Procédé selon l'une quelconque des revendications 1 à 5, caractérisé en ce que les entités d'une surface 3D donnée sont des noeuds du maillage de ladite surface 3D donnée.

7. Procédé selon l'une quelconque des revendications 1 à 6, caractérisé en ce que les entités d'une image 3D donnée sont des voxels de ladite image 3D donnée.

8. Procédé selon l'une quelconque des revendications 1 à 7, caractérisé en ce que ladite étape c) comprend une étape de minimiisation d'une énergie globale constituée d'une somme d'énergies locales, chaque énergie locale permettant de comparer une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante.

9. Procédé selon la revendication 8, caractérisé en ce que chaque énergie locale est composée d'une somme pondérée de: - un terme de correspondance locale E,(f,d), exprimant une mesure de cohérence d'un appariement entre une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante; - un terme de régularisation E2(f), exprimant des contraintes de lissage spatio- temporel du mouvement et/ou de la déformation estimée, lesdites contraintes étant telles que des vecteurs de déplacement associés à des entités voisines de la surface 3D doivent être proches en amplitude et/ou en direction.

10. Procédé selon la revendication 9, caractérisé en ce que ledit terme de correspondance locale E,(f,d) est fonction d'un terme exprimant une correspondance topologique entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante.

11. Procédé selon la revendication 10, caractérisé en ce que ledit terme de correspondance locale E,(f,d) est composé d'une somme pondérée dudit terme exprimant une correspondance topologique et d'au moins un des termes appartenant au groupe comprenant: - un terme exprimant une probabilité que l'entité de l'image 3D suivante appartienne à un élément de contour, dudit au moins un objet contenu dans ladite image 3D suivante; - un terme exprimant une distance entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante; - un terme lié aux densités de niveaux de gris.

12. Procédé selon l'une quelconque des revendications 8 à 11, caractérisé en ce que ladite étape de minimisation d'une énergie globale est basée sur une formulation statistique d'un processus d'estimation selon un modèle de Markov, telle que le mouvement et/ou la déformation à estimer est considéré comme une réalisation f = , i =1,..., NS} d'un champ aléatoire F = {F, , i =1,..., NS} regroupant toutes les réalisations possibles, les NS sites du champ aléatoire F étant donnés par toutes les entités (par ex. noeuds) de la surface 3D, obtenue pour l'image 3D courante, représentatives du maillage de ladite surface 3D, des labels desdits sites, aussi appelés étiquettes recherchées, étant les entités (par ex. voxels) de l'image 3D suivante jugées correspondantes desdites entités de la surface 3D.

13. Procédé selon l'une quelconque des revendications 1 à 12, caractérisé en ce que la surface 3D obtenue pour l'image 3D courante, lors de l'étape a) ou d), est considérée comme une surface 3D de résolution maximale, en ce que l'image 3D suivante comprise dans ledit ensemble d'images 3D est considérée comme une image 3D de résolution maximale, et en ce que ladite étape c) est effectuée dans un mode multi-résolution et comprend les étapes successives suivantes: cl) obtention (131) d'une part d'au moins une surface 3D de résolution inférieure pour l'image 3D courante, à partir de la surface 3D de résolution maximale pour l'image 3D courante, et d'autre part d'au moins une image 3D suivante de résolution inférieure, à partir de l'image 3D suivante de résolution maximale; c2) choix (132) de la plus basse des résolutions comme résolution courante, et initialisation des correspondances à la résolution courante; c3) pour la résolution courante, détermination (133) de correspondances optimales entre des entités de la surface 3D possédant ladite résolution courante pour l'image 3D courante, et des entités de l'image 3D suivante possédant également ladite résolution courante; c4) si la résolution courante n'est pas égale à la résolution maximale: * propagation (135) de résultats de l'exécution de l'étape c3), de la résolution courante à une résolution supérieure suivante, de façon à permettre une réinitialisation des correspondances; * retour à l'étape c3), en prenant comme nouvelle résolution courante la résolution supérieure suivante; sinon arrêt de l'étape c).

14. Procédé selon l'une quelconque des revendications 1 à 13, caractérisé en ce que ladite étape a) (12) comprend les étapes suivantes: segmentation (121) de l'image de départ, de façon à extraire des contours et points intérieurs audit au moins un objet contenu dans l'image de départ; - reconstruction (122) d'une surface 3D, représentée sous la forme d'un maillage comprenant des noeuds, à partir des contours et points intérieurs obtenus lors de l'étape de segmentation.

15. Procédé selon l'une quelconque des revendications 13 et 14, caractérisé en ce que, dans l'étape c2), l'initialisation des correspondances à la résolution la plus basse est réalisée par des mouvements nuls, chaque entité de la surface 3D obtenue pour l'image 3D courante étant mise en correspondance avec l'entité la plus proche dans l'image 3D suivante.

16. Procédé selon l'une quelconque des revendications 13 à 15, caractérisé en ce que, dans l'étape c4), la réinitialisation des correspondances est obtenue par interpolation, sur tout le maillage de la surface 3D possédant ladite résolution supérieure suivante, des correspondances optimales déterminées lors de l'étape c3).

17. Procédé selon l'une quelconque des revendications 14 à 16, caractérisé en ce que ladite étape a) comprend en outre l'étape suivante: - pré-traitement (123) de la surface 3D reconstruite, permettant de modifier, si nécessaire, ledit maillage de façon que la longueur de chacune des arêtes constituant une maille soit égale à une des distances entre voxels des images 3D.

18. Procédé selon l'une quelconque des revendications 1 à 17, caractérisé en ce que lesdites modalités d'imagerie appartiennent au groupe comprenant: - l'imagerie scanner; - l'imagerie par résonance magnétique; l'imagerie par ultrasons; - l'imagerie nucléaire.

19. Application du procédé selon l'une quelconque des revendications 1 à 18 dans un domaine d'imagerie appartenant au groupe comprenant: l'imagerie médicale; - l'imagerie marine; - l'imagerie satellitaire; l'imagerie des matériaux.

20. Produit programme d'ordinateur, caractérisé en ce qu'il comprend des instructions de code de programme pour l'exécution des étapes du procédé selon l'une quelconque des revendications 1 à 18, lorsque ledit programme est exécuté sur un ordinateur.

21. Moyen de stockage, éventuellement totalement ou partiellement amovible, lisible par un ordinateur, stockant un jeu d'instructions exécutables par ledit ordinateur pour mettre en oeuvre le procédé selon l'une quelconque des revendications 1 à 18.

22. Procédé de recalage tridimensionnel d'au moins un objet contenu dans une image 3D, caractérisé en ce qu'il comprend les étapes suivantes: estimation (41) du mouvement et/ou de la déformation subi(s), entre une première et une seconde images 3D acquises à deux instants différents et/ou à partir de deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas, par au moins un objet contenu dans lesdites première et seconde images 3D; - application (42) audit au moins un objet contenu dans la seconde image d'une fonction corrigeant le mouvement et/ou la déformation estimé lors de l'étape d'estimation, caractérisé en ce que ladite étape d'estimation (41) est mise en oeuvre selon le procédé 15 d'estimation de l'une quelconque des revendications 1, 2 et 6 à 18, sauf quand elles dépendent de l'une des revendications 3 à 5.