Brevet FR2977696A1 — publié le 20130111

La présente invention concerne le domaine général du traitement d'images et leur exploitation à des fins d'analyses notamment scientifiques et/ou médicale. L'invention concerne plus particulièrement un procédé de traitement pour déterminer les mouvements d'organes dans des processus d'imagerie médicale tels que l'imagerie par résonance magnétique (IRM) ou l'échographie. Le procédé de l'invention vise tout particulièrement l'application à l'estimation du mouvement cardiaque en imagerie par résonance magnétique cardiaque. L'imagerie par résonance magnétique cardiaque (IRM cardiaque) est utilisée en routine clinique pour quantifier les pathologies cardiaques et leur impact sur la fonction globale (volumes diastoliques et systoliques (ml), fraction d'éjection (%), masse ventriculaire (g)) ou régionale du muscle cardiaque. Les séquences d'imagerie à la base de tout examen d'IRM cardiaque sont dites "ciné-IRM". Ces séquences d'images permettent de différencier spontanément la paroi myocardique, de la cavité sanguine et des tissus avoisinants (poumons, médiastin, foie...). Le myocarde présente visuellement un aspect homogène et seules les variations cycliques de sa forme permettent une évaluation de la fonction mécanique cardiaque. En routine clinique, la fonction régionale est évaluée qualitativement à partir de l'épaississement régional de la paroi ou peut être estimée à partir de la mesure de l'épaisseur en diastole (phase de repos) et systole (phase de contraction maximum) du coeur, permettant de calculer le pourcentage d'épaississement segmentaire. Cette approche reste néanmoins peu reproductible en raison des erreurs de segmentation des interfaces endocardiques et épicardiques. L'exploration plus fine de la mécanique cardiaque est possible grâce à l'utilisation de séquences d'imagerie spécifiques, réalisées en complément des séquences cine-IRM, avec pour objectif déclaré de quantifier la fonction mécanique régionale et de s'assurer que la mesure soit réalisée dans les mêmes régions au cours du temps. Deux techniques ont pour se faire été développées : les séquences dites de "tagging" ou de marquage tissulaire d'une part et les séquences de contraste de phase avec codage des déplacements par écho stimulés d'autre part. Cette dernière technique est aussi connue sous l'acronyme "DENSE", pour "displacement encoding with stimulated echoes" en anglais. Les séquences de marquage tissulaire sont construites sur la base d'une séquence de ciné-IRM conventionnelle mais modifiées par l'ajout d'une préparation de l'aimantation spatialement localisée immédiatement après le pic R d'un électrocardiogramme R [1 ,2]. Il en résulte une séquence dynamique d'images avec superposition d'un motif géométrique (lignes, grille) ou tatouage tissulaire ayant la propriété de se déformer avec la déformation du tissu lui-même. La technique d'imagerie en contraste de phase permet d'obtenir des informations quantitatives sur le mouvement myocardique en se basant sur le fait que la phase du signal image varie en IRM en fonction du mouvement [3,4]. Longtemps limité par la puissance des gradients des imageurs, cette technique a pris un nouvel essor avec les imageurs récents dotés de gradients à hautes performances. Les séquences d'images dynamiques obtenues par marquage tissulaire et par contraste de phase constituent des séquences spécifiques, additionnelles et spécialisées, nécessitant obligatoirement des algorithmes de traitement spécifiques. On connaît dans la littérature d'autres méthodes d'estimation du mouvement myocardique dans des images ciné-IRM à partir de l'intensité des pixels des images (amplitude de l'image complexe) ou encore, dans les séquences d'images IRM par marquage tissulaire, en utilisant la phase monodimensionnelle linéarisée des images. Les méthodes basées sur l'exploitation de l'intensité des pixels estiment un champ dense ([5] [6] [7]) ou épars dans le cas du suivi des contours [8] avec une finalité de segmentation. D'autres modalités sont utilisées pour l'estimation du champ dense, comme par exemple en imagerie tomographique CT [9] ou en échocardiographie [13]. Les méthodes utilisant la phase monodimensionnelle des images IRM marquées (HARP [10], SinMod [11]) utilisent des signaux analytiques monodimensionnels afin de calculer la phase des images IRM marquées dans une seule direction, et estiment la composante du vecteur de déplacement dans cette direction. Ainsi, afin de trouver les deux composantes du vecteur, ces approches sont séparément appliquées deux fois. La phase monodimensionnelle est rendue linéaire en utilisant un filtrage préalable des images afin d'obtenir une image harmonique. Cela permet, une fois les images complexes calculées, de retrouver des phases linéaires par morceaux. Les séquences en ciné-IRM bien que très largement acquises ne sont généralement pas exploitées pour une quantification de la fonction myocardique locale principalement du fait de l'absence de disponibilité de méthodes performantes. A notre connaissance, très peu de logiciels sont disponibles pour réaliser cette étude. Le but de la présente invention est de proposer une technique alternative d'estimation du mouvement dans des séquences dynamiques d'images telles que des séquences ciné-IRM. L'invention vise tout particulièrement à proposer une technique d'estimation de mouvement qui ne nécessite pas l'acquisition de séquences spécifiques d'images pour sa mise en oeuvre comme pour les techniques de marquage tissulaire. Un but de l'invention est également de proposer une technique d'estimation de mouvement, notamment pour l'analyse de séquences ciné-IRM du myocarde, qui permette un post-traitement de séquences d'images ciné-IRM déjà acquises, afin de permettre un suivi longitudinal de patients par exemple. Ces différents buts sont atteints par un procédé de traitement d'une séquence dynamique d'images présentant une texture ou motif persistant au cours du temps dans la séquence, consistant à estimer et visualiser un champ de mouvement dense entre les images de la séquence dynamique en suivant les étapes suivantes : a) on calcule à partir d'une séquence de N images dynamiques d'intensité {I(x,k),xe [1,dimXjx[1,dimX2],kE [1,N]} une séquence d'images de phases 0(x, k) à deux composantes 01(x, k ), 0z (x, k ) , où x = (x,, xz ) représente les coordonnées spatiales des images dynamiques, u = (u,,uz) leurs coordonnées fréquentielles dans le domaine de Fourier et k =1,2,...,N, le temps discret dans une séquence composée de N images, et b) On réalise un maillage rectangulaire régulier des images de la séquence d'images de phase, chaque maille dudit maillage étant délimitée par quatre noeuds, et on calcule un vecteur déplacement en chaque noeud du maillage entre tout couple d'images successives de la séquence d'images de phase, c) On interpole le déplacement en tout point de chaque maille pour tout couple d'images de phase successives dans la séquence à deux instants successifs k=n et k=n+1, puis d) On définit des régions d'intérêt rectangulaires de largeur LI et de hauteur Lz dont le centre est situé en un noeud du maillage de l'image à l'instant k=n, et e) On calcule le déplacement des coins C de chaque région d'intérêt par différences des composantes de phase A01, 002 entre deux images successives de la séquence d'images de phases aux instants k=n et k=n+1, dans un repère local dont l'origine est située au centre de chaque région d'intérêt à l'instant k=n, et f) On calcule les 8 paramètres a d'une transformation bilinéaire locale du vecteur déplacement au noeud courant à l'instant k=n, g) On calcule le champ dense de mouvement par interpolation bilinéaire en tout point des mailles à partir du déplacement estimé aux 4 noeuds à l'étape f) et h) On affiche sur les images de la séquence dynamique d'intensité les champs denses de mouvement estimés. Le procédé de traitement de l'invention est basé sur l'utilisation de la phase 30 spatiale bidimensionnelle des images d'une séquence dynamique d'images, sans filtrage préalable des images. Ce procédé est ainsi applicable sur les images de toute séquence d'images ciné-IRM que comporte tout examen IRM cardiaque. De plus, par rapport aux méthodes existantes, le procédé de l'invention estime des mouvements locaux non rigides, contrôlés par un modèle paramétrique bilinéaire. Le procédé de l'invention se distingue notamment des méthodes existantes telles que les méthodes HARP et SinMod précitées par les points suivants : - il permet un calcul direct de la phase spatiale bidimensionnelle des images qui permet d'estimer le vecteur déplacement en chaque point des images, - la phase spatiale calculée n'est pas linéarisée mais modélisée par un polynôme à deux variables de degré supérieur à 1, ce qui permet de s'affranchir d'un filtrage initial des images de la séquence dynamique d'images traitée, et - l'estimation locale du mouvement est obtenue par un maillage déformable contrôlé par un modèle de déplacement bilinéaire. Selon une variante de réalisation du procédé de l'invention, on déforme la région d'intérêt courante en appliquant le modèle bilinéaire calculé et on réitère les étapes d) à f) pour l'ensemble des régions d'intérêt de l'image n. Le procédé de l'invention est alors un procédé itératif pour lequel deux itérations suffisent. Par ailleurs, selon le procédé de l'invention, le calcul de la séquence d'images de phase 0(x, k) à l'étape a) comporte les opérations suivantes : ai) on calcule la transformée de Fourier 2D directe S(u,k) de l'image d'intensité I(x, k) de la séquence dynamique, a2) on calcule ensuite la transformée de Fourier 2D inverse du premier quadrant de Fourier h (x, k) de g(u, k) , et du deuxième quadrant de Fourier Îz (x, k) de 3(u,k) pour obtenir une séquence d'images dynamiques complexes à deux composantes h (x, k) et î2 (x, k ) a3) on calcule la séquence d'images de phase 0(x, k) à deux composantes (x, k ) , 02 (x, k) telle que : 0, (x, k) - arg (TI (x, k) ) 02(x,k)- arg (22(x,k)) Par ailleurs selon l'invention, à l'étape b) on balaie le maillage de noeud en noeud pour calculer le vecteur déplacement en chaque noeud, la valeur initiale du vecteur déplacement en un des noeuds étant obtenue à partir du déplacement estimé sur le noeud précédent. De façon avantageuse, à l'étape e) on estime le déplacement des coins C des régions d'intérêt rectangulaires de largeur LI, et de hauteur Lz en fonction de la différence entre les images de phase çi(x',n) et 0(x',n+l) exprimées dans un repère local x' centré sur le noeud courant Nc, tel que : 0 (x',n) =~i (x',n+l)-~i (x',n) A02 (x',n) =02 (x',n+1)-02 (x',n) L L avec x'= [x'l ,x', ] pour x'1 E - 2 2 et x'2 E Toujours selon l'invention à l'étape e), on détermine, après avoir estimé les différences de phase A01, A02, les déplacements locaux de chaque coin C de chaque région d'intérêt en moyennant les différences d'images de phase appelées et 002, et on calcule les déplacements des coins Ci (i =1 à 4) appelés du (Ci ) et dy (Ci) tels que : LZ LZ 2'2 [du(Ci) _ d,(C;) 0' 0' - iXl lx2 0' 0' 2X, 2X2 - za9 où X' a, , , i =1, 2 et _, x~ représentent les dérivées spatiales tard (x» moyennes (sur x') des phases 01(x, k ) et 02 (x, k ) A noter que ces dérivées spatiales sont calculées par des différences finies dans le cas des images Ciné-IRM. Pour des images IRM marquées, elles sont remplacées par les expressions suivantes : ~lx = 27Gf1, ~ix = 27Gf2, (î)2x = -27Gf1, (î)2xz = 27Gf2 où fl et fz sont les fréquences fondamentales du marquage tissulaire supposées connues selon les deux directions x'1 et x'z du repère local x'. Diverses autres caractéristiques ressortent de la description faite ci-dessous en référence aux dessins annexés qui explicite le procédé de traitement d'images de l'invention et des exemples d'applications pratiques de celui-ci en imagerie par résonance magnétique. Sur les figures: - La figure 1 exprime le calcul des images de phase pour une séquence de N images à k fixé (k =1,...,N ) ; - la figure 2 représente un schéma d'estimation du mouvement dans une séquence de Nimages de phase à 2 composantes selon l'approche eulérienne ; - la figure 3 représente un schéma synoptique du procédé de traitement d'image par estimation du mouvement selon la présente invention ; - les figures 4A et 4B représentent le maillage de deux images de phase successives «x,n) et 0(x,n+l) et l'estimation du mouvement élastique sur ces deux images de phase consécutives à deux instants k=n et k=n+l autour d'un noeud courant Nc dans l'image 0(x,n) ; - la figure 5 représente schématiquement l'estimation à deux itérations du modèle bilinéaire pour une région d'intérêt donnée suivant le procédé de l'invention, - la figure 6 décrit le principe de calcul du champ dense de déplacement par interpolation bilinéaire ; - les figures 7A à 7C représentent respectivement une image par résonance magnétique cardiaque de marquage tissulaire, notée I(x,k) sur la figure 1, une paire d'images de phase 01(x,k) et 0(x,k) obtenues à partir de l'image I(x,k) conformément au calcul présenté sur la figure 1 ; - la figure 8 représente l'application du processus d'estimation de mouvement décrit sur la figure 3 à une séquence dynamique d'images par résonance magnétique cardiaque de marquage tissulaire analogue à celle de la figure 7a ; - les figures 9A à 9C représentent respectivement une image par résonance magnétique cardiaque ciné, notée I(x,k) sur la figure 1, une paire d'images de phase 01(x,k) et 0(x,k) obtenues à partir de l'image I(x,k) conformément au calcul présenté sur la figure 1 ; - la figure 10 représente l'application du processus d'estimation de mouvement décrit sur la figure 3 à une séquence dynamique d'images par résonance magnétique cardiaque ciné analogue à celle de la figure 9a ; - les figures 11A et 11B représentent respectivement une image par RM marquée extraite de la séquence d'images utilisée pour estimer le mouvement cardiaque et les trajectoires de mouvement estimées des quatre blocs représentés sur la figure 11A ; - La figure 12 présente un mode d'exploitation logiciel par cartographie des champs de déplacement estimés avec la méthode de l'invention dans les images IRM de la figure 11A. La présente invention propose un nouveau procédé de traitement de séquences d'images dynamiques par résonance magnétique présentant une texture ou motif persistant au cours du temps dans la séquence, ci-après notées « images RM » dans la suite de la présente description, en particulier d'images RM cardiaques, pour déterminer et quantifier les mouvements du muscle cardiaque et ses parois en particulier pour aider au diagnostic de certaines pathologies. Le procédé de l'invention, contrairement aux procédés connus et employés à ce jour consiste essentiellement à exploiter la phase spatiale de ces images RM en combinaison avec une technique de mise en correspondance de blocs déformables définis dans les images RM, pour déterminer précisément et dans toute séquence d'images RM des champs de mouvement représentant fidèlement les mouvements de l'organe observé sur les images RM. En effet, les techniques principalement employées à ce jour pour évaluer les mouvements d'organes cardiaques sur des images RM consistent essentiellement en l'exploitation des intensités lumineuses des pixels des images pour évaluer les variations d'épaisseurs des parois de l'organe puis reconstituer par calcul les volumes de l'organe et suivre leur évolution dans le temps. Cependant ces méthodes comportent un fort risque d'erreur et reposent grandement sur les capacités propres d'analyse des médecins étudiant les images. Le procédé de l'invention vise à proposer une méthode fiable et reproductible de traitement de séquences d'images dynamiques RM permettant de retranscrire et évaluer localement et le plus fidèlement possible les mouvements de l'organe observé pour permettre un diagnostic amélioré par tous les praticiens. Ce procédé va être décrit ci-après en détail en référence aux figures 1 à 6 annexées. En particulier, le procédé de l'invention est présenté sous forme de synoptique à la figure 3. Selon l'invention, le procédé de traitement d'image proposé consiste de façon générale essentiellement à : - calculer la phase spatiale 2D d'images RM d'une séquence dynamique d'images RM pour estimer le vecteur déplacement, - modéliser, sans linéarisation, la phase spatiale 2D par un polynôme à deux variables de degré supérieur à 1, et - estimer localement le mouvement au sein des images par un maillage déformable contrôlé par un modèle de déplacement bilinéaire. Dans la description qui suit, on considère pour expliciter le procédé de l'invention une séquence S de N images RM, N étant un nombre entier supérieur ou égal à 1, d'intensité notée mathématiquement {I(x, k), xE [1, dim XI]x[1, dim X2], k e [1, N]} . Pour cette séquence S, on note x=(xl, x2) les coordonnées spatiales de chaque image constituant ladite séquence S, et u=(ul, uz) leurs coordonnées fréquentielles dans le domaine de Fourier, k=1, 2, 3,...,N, étant le temps discret dans une séquence composée de N images. La première étape du procédé de traitement d'images de l'invention consiste, comme représenté sur la figure 1 et la figure 3, à calculer, à partir des N images RM de la séquence dynamique S, une séquence d'images de phase, notée O(x, k), à deux composantes 01(x, k), 02(x, k) selon le principe décrit sur la figure 1. Ainsi, on projette dans un premier temps la séquence S d'images RM d'intensité I(x,k) dans le domaine fréquentiel de Fourier par calcul de la transformée de Fourier 2D, notée 3(u, k), de l'intensité I(x, k) de la séquence S. La projection de la séquence S dans le domaine de Fourier 3(u, k) comporte alors 5 deux composantes I(u,k), 32(u,k), correspondant aux coordonnées fréquentielles de ladite projection dans les premier et second quadrants du plan de Fourier. Ces composantes I(u,k), 32(u,k) sont ensuite soumises chacune au calcul de la transformée inverse de Fourier 2D pour obtenir deux séquences d'images complexes II (x,k), ÎZ (x,k), dont on calcule ensuite l'argument pour obtenir les 10 composantes d'images de phase 01(x, k), 02(x, k) de la séquence d'images de phase O(x, k)= (01(x, k) 02(x, k)). A partir de cette séquence d'image de phase O(x, k), représentée schématiquement sur la figure 2, on va ensuite déterminer le champ dense de mouvement, noté (u(x,n) v(x,n)) aux instants k=n et k=n+I. 15 La seconde étape du procédé de l'invention consiste ensuite, comme indiqué dans le synoptique de la figure 3, à traiter cette séquence d'images de phase O(x, k)= (01(x,k) 0z(x,k)) à deux composantes pour les subdiviser suivant un maillage rectangulaire comprenant aux quatre coins de chacune des mailles un noeud Ainsi, comme représenté à la figure 4, aux instants k=n les images de phase 20 01(x,n) et 0z(x,n) sont subdivisées selon un maillage rectangulaire, dont les paramètres de maille sont déterminés en fonction de la dimension des structures dans l'image. Aussi, pour chaque noeud du maillage, il convient de calculer le vecteur déplacement obtenu après déformation du maillage sur les images de phase 01(x,n) et 0z(x,n) et identification au modèle paramétrique du mouvement 25 comme cela va être décrit ci-dessous. Une fois les images de phase subdivisées selon le maillage appliqué, ledit maillage est ensuite balayé de noeud en noeud suivant un chemin de lecture, de préférence de gauche à droite et de haut en bas, pour estimer le déplacement en chaque noeud du maillage selon le processus décrit figure 3. Ainsi, considérant un 30 noeud quelconque du maillage, que l'on appelle noeud courant Nc, la valeur initiale des coordonnées spatiales de ce noeud est obtenue en retenant la partie entière du déplacement estimé sur le noeud précédent. On calcule ensuite un vecteur déplacement en chaque noeud du maillage entre tout couple d'images successives de la séquence d'images de phase puis on interpole le déplacement en tout point de chaque maille pour tout couple d'images de phase successives dans la séquence aux deux instants successifs k=n et k=n+1. On définit ensuite dans une troisième étape, comme représenté aux figures 3 et 5, des régions d'intérêt rectangulaires de largeur LI et de hauteur Lz autour de chaque noeud N du maillage de l'image à l'instant k=n, et on calcule ensuite, dans une quatrième étape le déplacement des coins Ci (i=1 à 4) de chaque région d'intérêt par différences des composantes de phase AOI, A02 entre deux images successives de la séquence d'images de phase aux instants k=n et k=n+1, dans un repère local dont l'origine est située au centre de chaque région d'intérêt à l'instant k=n, c'est-à-dire en chacun des coins Ci desdites régions d'intérêt. Cette démarche permet une estimation du mouvement élastique entre les deux images de phase successives aux instants k=n et k=n+1 autour du noeud courant Nc. Dans la représentation donnée aux figures 4 et 5, le mouvement est estimé en recherchant les correspondants des coins CI à C4 d'un bloc rectangulaire de largeur LI et longueur Lz de l'image 0(x, n) de référence autour d'un noeud Nc dans l'image 0(x, n+1) cible. L'estimation du mouvement est effectuée sur chacun des noeuds de la grille de maillage des images discrète selon un pas régulier horizontal G1 et vertical G2. On obtient ainsi la déformation locale bilinéaire de chaque région d'intérêt. L'estimation des translations des coins Ci se fait par combinaison de la 25 différence des images de phase. Ainsi, pour chacun des 4 coins Ci (i=1 à 4) d'une région d'intérêt autour d'un noeud Nc du maillage de l'image de phase le déplacement est estimé par bloc de taille LixL2 dans un repère x' attaché à chaque noeud courant Nc, et définit selon le calcul suivant: 30 4O1(x',n)=O1(x',n+1)-01(x',n) etA02(x',n)=02(x',n+1)-02(x',n) L avec x'=[x'1, x,21 pour x'1 e - z , zL L et pour x'2 E [-L z , z et, en considérant, sans perte de généralité, que LI et L2 sont des entiers pairs. En moyennant les différences de phase 4O1 (x', n) et 402 (x', n) et après avoir éliminé les sauts de phase on obtient deux équations, avec deux inconnues et, les déplacements à estimer pour chacun des coins Ci. Les deux différences de phase moyennes, notées 4O1 (x', n) et 4~2 (x', n), permettent de remonter aux estimations analytiques des déplacements du (Ci) et dy (Ci) en chaque coin Ci (i =1 à 4) autour du noeud courant par les expressions suivantes : Oi' Xl iX2 02xi ' ; 02'x 2 0 Oi A 02 '':)i où ~P,; _ , i =1,2 et j = x~, xz représentent les dérivées x. -1 tard (x» remplacées par les expressions i = 2nf1, ~lx = 2nf2, 0' . = -2nf1,= 2nf2 où fl et f2 sont les fréquences fondamentales du marquage tissulaire supposé connu sur les deux directions x'1 et x'2 du repère local x'. Une fois les déplacements des 4 coins connus, l'étape suivante du procédé de 20 l'invention consiste à déterminer les paramètres de la transformée bilinéaire locale pour la région d'intérêt considérée autour du noeud courant Nc étudié. Ces spatiales moyennes (sur x') des phases 0,(x,k) et 02(x,k) A noter que ces dérivées spatiales sont calculées par des différences finies dans le cas des images ciné-IRM. Pour les images IRM marquées, elles sont suivantes paramètres sont notés a = (au bu cu ~av bv cvd \ T u v On considère pour cela le système de coordonnées locales x'=[x'I, x'21, avec l'origine dans le noeud courant Nc, tel que précédemment défini et dans lequel les coordonnées spatiales des images de phase peuvent être définies par : Iu(x'rx'2)=au.x'I+bu.x'2+cu.x'I.x'2+du v(x'I,x'z)=a,, -x'I+b,, -x'z+c,, -x'I-x'z+dv 5 En prenant par la suite la notation d = rdu(C,) du(C2) du(C3) du(C4)~T \.dU(Cl) dv(Ca) dv(C3) dv(C4)i pour définir les translations estimées des coins C (i=1 à 4) de la région d'intérêt autour du noeud courant Ncconsidéré les paramètres a de la transformée bilinéaire locale correspondant sont alors calculés par la relation suivante : 1 1 1 1 _ LI LI LI _ LI 1 1 1 1 L2 L2 L2 L2 d 2 2 2 2 LI L2 LI L2 LI L2 LI L2 1 1 1 1 10 2 2 2 2 i où LI et L2 sont la largeur et la hauteur de la région d'intérêt étudiée. Une fois les paramètres a de la transformée bilinéaire locale calculés pour une région d'intérêt autour d'un noeud considéré du maillage d'une image de phase, on bénéficie alors d'un modèle de transformation bilinéaire applicable pour déformer 15 toute zone de l'image de phase de façon itérative en deux itérations, la deuxième itération ayant pour but d'affiner la détermination des valeurs des paramètres du modèle bilinéaire par déformation de l'image de phase de référence. On procède donc selon une boucle B1 représentée à la figure 3 à la nouvelle itération de confirmation du modèle bilinéaire déterminé ; 20 Une fois la dernière itération de déformation réalisée on reproduit ensuite en chaque noeud du maillage des images de phase l'estimation de vecteurs déplacement des coins d'une région d'intérêt centrée sur chaque noeud puis le calcul de la transformée bilinéaire correspondante. Cette étape est réalisée selon une boucle B2, telle que représentée sur la figure 3. On peut ensuite calculer, une fois tous les noeuds du maillage testés, le champ dense du mouvement obtenu en considérant les déplacements de tous points A d'une image et de tous les noeuds. Pour les noeuds, le déplacement est obtenu en utilisant le modèle bilinéaire estimé dans la région d'intérêt correspondante et au point de coordonnées (0,0) dans le repère x'=(x'I, x'z) précédemment défini. /a (N.) b (N.) Ç (Ni) d .(Ni)' On note alors a~ (N~) = ay (N~) by (N~) cy (N~) dy (N~) les paramètres du modèle bilinéaire calculés au noeud Ni. Pour les points A, leur déplacement (u(xA,n),v(xA,n)) est obtenu en moyennant les estimations du champ déplacement des 4 noeuds voisins, Ni (i=1 à 4) selon l'expression suivante : (u(xA, n), v(xA, n))T =1 / 4r. aT (Ni)g(Ni) Où les valeurs des distances de chacun des 4 noeuds au point A, notées g(Ni) et déterminées d'après la figure 6, donnent : g(NI) = 1/2(- Gl g(N2)=1/2(Gi g(N3) =1 / 2(Gi 20 g(N4) =1/2(-GI Ainsi, en prenant pour exemple une image de phase quelconque d'une séquence d'images de phase telle que représentée sur la figure 6, dans laquelle on considère un point A central et quatre noeuds NI, N2i N3, N4 délimitant une maille de côtés G1, G2 autour de ce point A et en considérant 4 régions d'intérêt de 25 dimensions LIxL2 contenant chacun des noeuds respectivement on obtient pour chacun des noeuds les valeurs de g(N) représentant la distance de chacun des 4 noeuds au point A considéré. -G2 GIG2/2 2Y -G2 -GIG2/2 2)T G2 GIG2 /2 2Y, et Gz -GIG2/2 2Y. Cette façon de calculer le champ dense du mouvement offre d'une part l'avantage d'obtenir des champs de mouvement réguliers, et permet d'autre part d'éliminer des éventuelles estimations aberrantes. On reproduit ensuite pour chaque image de la séquence d'image de phase l'ensemble du traitement appliqué à une image tel que précédemment décrit, suivant une boucle B3 telle que représentée à la figure 3, jusqu'à ce que toutes les images de la séquence aient été traitées, ce après quoi on affiche les résultats sur un moniteur ou écran connecté à un dispositif de calcul et de traitement informatique utilisé de façon préféré pour mettre en oeuvre le procédé de l'invention. Les figures 7 à 12 représentent différents résultats de mise en oeuvre concrète du procédé de l'invention sur des séquences d'images RM cardiaques. La figure 7A représente une image RM cardiaque de marquage tissulaire, notée I(x,k) sur la figure 1, extraite d'une séquence d'images RM cardiaque particulière. Lorsqu'on applique la procédé de l'invention sur cette image on obtient dans une premier temps, comme représenté sur la figure 1, une paire d'images de phase 01(x,k) et 0(x,k) obtenues à partir de l'image I(x,k), représentées aux figures 7B et 7C. Comme cela ressort des figures 7B et 7C, les images de phases 01(x,k) et 0(x,k) présentent un rendu texturé, conférant une impression de relief des images. Par application du processus d'estimation de mouvement décrit sur la figure 3 à une séquence dynamique d'images par résonance magnétique cardiaque de marquage tissulaire comportant six images RM marquées analogues à celles de la figure 7A on peut déterminer le mouvement au sein de l'organe imagé, comme représenté à la figure 8. On peut ainsi constater sur la figure 8 qu'après application du procédé d'estimation de l'invention les images RM de la séquence considérée sont dotées de flèches représentant le champ de mouvement au sein des images, ce champ de mouvement variant dans chaque image de la séquence, permettant ainsi une observation du mouvement au sein du coeur observé dans l'image, ce mouvement étant indépendant d'une déformation du marquage. De façon analogue, les figures 9A à 9C représentent respectivement une image cardiaque ciné-IRM, et une paire d'images de phase 01(x,k) et 0(x,k) obtenues à partir de l'image I(x,k) conformément au calcul présenté sur la figure 1. Après application du processus d'estimation de mouvement décrit sur la figure 3 à une séquence dynamique d'images cardiaques ciné-IRM comportant 6 images on obtient une séquence de 6 images dotées d'une représentation interne du champ de mouvement au sein des images, telle que représentée à la figure 10. Comme pour la figure 8, les images RM de la séquence considérée sont dotées de flèches représentant le champ de mouvement au sein des images, ce champ de mouvement variant dans chaque image de la séquence, permettant ainsi une observation du mouvement au sein de l'organe observé dans l'image. Le procédé de l'invention peut s'avérer particulièrement avantageux avec des séquences d'images ciné-IRM dans la mesure où il peut être appliqué à des séquences d'images anciennes, conservées dans le dossier médical d'un patient, sans nécessiter un quelconque retraitement. On peut ainsi dresser un historique de la pathologie du patient et observer son évolution dans le temps. Les figures 11A et 11B donnent une représentation de l'utilité et la performance du procédé de l'invention pour déterminer le mouvement au sein d'images RM et permettre au praticien d'en tirer un diagnostic par son interprétation des images obtenues. La figure 11A représente une image RM marquée extraite de la séquence d'images utilisée pour estimer le mouvement cardiaque et les trajectoires de mouvement estimées de quatre blocs définis au niveau de deux zones à risque du coeur à observer (blocs 1, 2) ainsi que dans deux zones normales (blocs 3, 4). La figure 11B représente les trajectoires de mouvement calculées pour chacun des blocs définis dans les images de la figure 11A après mise en oeuvre du procédé de l'invention. On peut ainsi observer que l'amplitude de mouvement des blocs 1, 2 couvrant les zones à risque du coeur est deux à trois fois plus faible sur le cycle cardiaque que pour les deux blocs 3, 4 couvrant les régions normales. Il est ainsi possible très rapidement de déterminer l'existence d'une pathologie cardiaque en observant les trajectoires de mouvement obtenues après traitement de séquences d'images cardiaques d'un patient. Enfin, la figure 12 présente une façon d'exploiter les champs de déplacement estimés avec la méthode de l'invention pour que le médecin puisse caractériser la fonction mécanique du myocarde d'un patient. Des cartographies de paramètres peuvent être construites et présentées sous forme d'image couleur comme présentée (en noir et blanc sur la figure 12). Ces images de cartographies de paramètres peuvent notamment être obtenues par exemple à l'aide du logiciel inTag, (http://www.creatis.insa-lyon.fr/inTag/). La figure 12 présente ainsi une cartographie d'amplitude du mouvement qui donne en chaque point du myocarde l'amplitude du déplacement (correspondant à la longueur des vecteurs déplacement comme représenté figure 8 ou figure 10). Ceci procure au médecin une aide déterminante dans l'établissement d'un diagnostic de pathologie cardiaque. Bib/ioaravhie [1] E. A. Zerhouni, D. M. Parish, W. J. Rogers, A. Yang, and E. P. Shapiro, "Human heart: tagging with MR imaging - A method for noninvasive assessment of myocardial motion," Radiology, vol. 169, pp. 59-63, 1988. [2] L. Axai and L. Dougherty, "MR imaging of motion with spatial modulation of magnetization," Radio%ogy, vol. 171, pp. 841-5, Jun 1989. [3] A. H. Aletras, S. Ding, R. S. Balaban, and H. Wen, "DENSE: displacement encoding with stimulated echoes in cardiac functional MRI," J Magn Reson, vol. 137, pp. 247-52, Mar 1999. [4] A. H. Aletras, R. S. Balaban, and H. Wen, "High-resolution strain analysis of the human heart with fast-DENSE," J Magn Reson, vol. 140, pp. 41-57, Sep 1999. [5] S. C. Amartu and H. J. Vesselle, "A new approach to study cardiac motion: The optical flow of cine MR images," Magnetic Resonance in Medicine, vol. 29, pp. 59-67, 1993. [6] A. C. Planat and I. E. Magnin, "Estimation de mouvement par maillage actif avec prise en compte de discontinuités," Traitement du Signal, vol. 18, pp. 583-401, 2001. [7] X. Papademetris, A. J. Sinusas, D. P. Dione, R. T. Constable, and J. S. Duncan, "Estimation of 3-D left ventricular deformation from medical images using biomechanical models," IEEE Transactions on Medical Imaging, vol. 21, pp. 786-800, Jul 2002. [8] P. Shi, A. Sinusas, R. T. Constable, E. Ritman, and J. S. Duncan, "Pointtracked quantitative analysis of LV motion from 3D image sequences," IEEE Transactions on Medical Imaging, vol. 19, pp. 36-50, 2000. [9] J. M. Gorce, D. Friboulet, and I. E. Magnin, "Estimation of threedimensional cardiac velocity fields: Assessment of a differential method and application to 3D CT data," Medical Image Analysis, vol. 1, pp. 245-261, 1997. [10] N.F. Osman, J.L. Prince, Visualizing myocardial function using HARP MRI, Phys. Med. Biol., 2000, Vol. 45, p. 1665-1682. [11] T. Arts, F.W. Prinzen, T. Delhaas, J. Milles, P. Clarysse, Mapping displacement and deformation of the heart with local sine wave modelling, IEEE Trans. on Medical Imaging, vol. 29, n° 5, pp. 1114-1123, may 2010. [12] A. Basarab, Estimation du mouvement dans des séquences d'images 5 échographiques : application à l'élastographie de la thyroide, Thèse de Doctorat, INSA Lyon, 2008, 135p. [13] J. Crosby, B. H. Amundsen, T. Hergum, E. W. Remme, S. Langeland, and H. T. H, "3-D speckle tracking for assessment of regional left ventricular function," U/trasound in Medicine & Bio%gy, vol. 35, pp. 458-471, 2009. 10 REVENDICATIONS1 - Procédé de traitement d'une séquence dynamique d'images présentant une texture ou motif persistant au cours du temps dans la séquence, consistant à estimer et visualiser un champ de mouvement dense entre les images de la séquence dynamique en suivant les étapes suivantes : a) on calcule à partir d'une séquence de N images dynamiques d'intensité {I(x, k), x e [1, dim X~ ]x [1, dim x k E [1, N]} une séquence d'images de phases 0(x, k) à deux composantes 0, (x, k), 02 (x, k), où x = (xl, x2) sont les coordonnées spatiales des images dynamiques, u = () leurs coordonnées fréquentielles dans le domaine de Fourier et k =1,2,...,N, le temps discret dans une séquence composée de N images, et b) On réalise un maillage rectangulaire régulier des images de la séquence d'image de phase, chaque maille dudit maillage étant délimitée par quatre noeuds, et on calcule un vecteur déplacement en chaque noeud du maillage entre tout couple d'images successives de la séquence d'images de phase, c) On interpole le déplacement en tout point de chaque maille pour tout couple d'images de phase successives dans la séquence à deux instants successifs k=n et k=n+1, puis d) On définit des régions d'intérêt rectangulaires de largeur LI et de hauteur Lz 20 dont le centre est situé en un noeud du maillage de l'image à l'instant k=n, et e) On calcule le déplacement des coins C de chaque région d'intérêt par différences des composantes de phase 41, 002 entre deux images successives de la séquence d'images de phases aux instants k=n et k=n+1, 25 dans un repère local dont l'origine est située au centre de chaque région d'intérêt à l'instant k=n, et f) On calcule les 8 paramètres a d'une transformation bilinéaire locale du vecteur déplacement au noeud courant à l'instant k=n, g) On calcule le champ dense de mouvement par interpolation bilinéaire en tout 30 point des mailles à partir du déplacement estimé aux 4 noeuds à l'étape f) eth) On affiche sur les images de la séquence dynamique d'intensité les champs denses de mouvement estimés. 2 - Procédé de traitement selon la revendication 1, caractérisé en ce que l'on déforme la région d'intérêt courante en appliquant le modèle bilinéaire calculé et on réitère les étapes d) à f) pour l'ensemble des régions d'intérêt de l'image n. 3 - Procédé de traitement selon la revendication 1 ou 2, caractérisé en ce qu'à l'étape a), le calcul de la séquence d'images de phase 0(x, k) comporte les opérations suivantes : ai) on calcule la transformée de Fourier 2D directe 3(u,k) de l'image d'intensité I(x, k) de la séquence dynamique, a2) on calcule ensuite la transformée de Fourier 2D inverse du premier quadrant de Fourier h(x,k) de 3(u,k) , et du deuxième quadrant de Fourier Îz(x,k) de 3(u,k) pour obtenir une séquence d'images dynamiques complexe à deux composantes h(x,k) et Îz(x,k) a3) on calcule la séquence d'images de phase 0(x,k) à deux composantes (x, k ) , 02 (x, k) telle que : 0,(x,k)= arg (h(x,k)) 02(x,k)= arg (Îz(x,k) 4 - Procédé selon la revendication 3, caractérisé en ce qu'à l'étape b) on balaie le maillage de noeud en noeud pour calculer le vecteur déplacement en chaque noeud, la valeur initiale du vecteur déplacement en un des noeud étant obtenue à partir du déplacement estimé sur le noeud précédent. 5 - Procédé selon les revendications 1 à 3, caractérisé en ce qu'à l'étape e) on estime le déplacement des coins C des régions d'intérêt rectangulaires de largeur LI, et de hauteur Lz en fonction de la différence entre les images de phase çi(x',n) et çi(x',n+l) exprimées dans un repère local x' centré sur le noeud courant Nc, tel que : 401(x',n)=01(x',n+l)-01(x',n) 4 02(x', n) _ 02(x', n+ l )- 02(x', n)avec x'=[x'l,x'2] pour x'I E L- 21 , 211 et x'2 E L- 22 , 22J 6 - Procédé selon la revendication 5, caractérisé en ce que à l'étape e) on détermine, après avoir estimé les différences de phase 4(te1 (x', n) et 4(I)2 (x', n), moyenné les différences d'images de phase, deux équations avec deux inconnues et, les déplacements du (ci) et d, (ci) à estimer en chaque coin Ci (i =1 à 4) autour du noeud courant par les expressions suivantes : 1x, 02x1 _ a; où 0Z; =tard (x') ,1 =1,2 et j = x1, x2 représentent les dérivées spatiales moyennes (sur x') des phases 01(x, k) et 02(x, k) 7 - Procédé selon les revendications 1 à 5, caractérisé en ce que le déplacement de tout point de la région d'intérêt centrée au noeud courant Nc et régi par un modèle bilinéaire tel que : Iu(x'I,x'2)=au x'I+bu x'2+cu x'I x'2+du v( x'I , x'2 ) = a' - x'l + b' - x'2 + cv - x'I - x'2 + 02x'2. _ où x'1 ,x'2 sont les coordonnées dans le repère local et a = bu cu T ~a b d cv v v v les paramètres du modèle, estimés à partir des déplacements des 4 coins de la région d'intérêt selon :1 1 1 1 LI LI LI LI 1 1 1 1 L2 L2 L2 L2 d 2 2 2 2 LI L2 LI L2 LI L2 LI L2 1 1 1 1 2 2 2 2 , les déplacements estimés des ~dv(C1) dv(Cz) dv(C3) dv(C4)i 4 coins de la région d'intérêt. (du(CI) du(C2) du(C3) du(C4)~T avec d =5