Brevet FR2736454A1 — publié le 19970110

i PROCEDE DE RECONSTRUCTION D'IMAGES TRIDIMENSIONNELLES SUR UN OBJET MOBILE OU DÉFORMABLE DESCRIPTION L'invention concerne un procédé de reconstruction d'images tridimensionnelles d'un objet mobile ou déformable. L'imagerie tridimensionnelle permet de reconstruire les valeurs d'un paramètre particulier de l'objet étudié après avoir acquis une série de projections bidimensionnelles de cet objet. En tomographie aux rayons X, on recueille une information de densité car on obtient des résultats représentatifs de l'atténuation linéaire d'un rayonnement. En tomographie d'émission, c'est l'activité locale d'émission d'un corps marqueur qui est calculée. Les techniques de reconstruction d'images s'appliquent également dans le cas de transmission de neutrons ou d'émission de rayonnement y ou de positons. Historiquement, on a d'abord obtenu des coupes d'objets à l'aide d'une ligne de détecteur. Actuellement, l'imagerie tridimensionnelle consiste à utiliser les réseaux de détecteurs bidimensionnels pour en déduire le contenu du volume de l'objet. Les problèmes sont beaucoup plus compliqués pour les images tridimensionnelles à cause du volume d'informations plus important à traiter et éventuellement de situations géométriques plus variées à considérer. Cependant, le problème à résoudre peut toujours être représenté par l'équation: m = Hf+b (1) o m désigne l'ensemble des mesures, f l'image à retrouver et H l'opérateur de projection; b est le B 12254.3/Jci bruit de mesure. La reconstruction équivaut à l'inversion de cette équation, qui peut être résolue par deux familles de méthodes. Les méthodes analytiques consistent à exploiter une formule ou un algorithme dans lesquels les mesures sont combinées entre elles pour obtenir les valeurs aux points de l'image. Il existe de nombreux exemples, qui impliquent souvent, par exemple quand les capteurs sont focalisés vers un foyer unique en définissant un faisceau de détection conique, l'emploi d'une transformée de Radon de l'image, c'est-à-dire d'une fonction mathématique par laquelle on fait correspondre à chaque point de mesure la somme de la propriété mesurée sur tous les points d'un plan, appelé plan caractéristique, passant par ce point. Un des inventeurs est l'auteur de plusieurs demandes antérieures qui portent sur ces questions; en particulier, il faut citer le brevet européen 0 292 402, qui constitue l'invention de base dans laquelle les caractéristiques et les propriétés des méthodes analytiques d'inversion et en particulier de la transformée de Radon sont rappelées. Une caractéristique de ces méthodes est que la discrétisation des résultats obtenus, pour pouvoir les associer aux points de l'image effectivement retenus par une discrétisation de l'objet, est produite à la fin des calculs. On peut opposer les méthodes algébriques dans lesquelles la discrétisation de l'image est effectuée dès le début, de sorte qu'on obtient un système d'équations o chaque mesure est exprimée en fonction de la contribution des points de l'image que traverse la direction de scrutation du détecteur responsable de la mesure. Le problème se ramène alors ouvertement à l'inversion de la matrice de projection H. Ces méthodes ne sont guère applicables quand les B 12254.3 JCI points sont nombreux, à cause de la taille des calculs, mais le problème est du moins toujours soluble, ce qui n'est pas toujours possible avec des méthodes analytiques dont les formules d'inversion peuvent être insuffisantes pour obtenir des résultats corrects. Un autre avantage des méthodes algébriques est que des informations connues à l'avance ou arbitraires peuvent être introduites beaucoup plus facilement. Si par exemple le contenu d'une partie de l'image est connu à l'avance, il est facile de l'introduire directement dans les équations, ce qui n'est pas le cas avec les méthodes analytiques. En pratique, l'inversion de la matrice est généralement menée suivant les techniques itératives habituelles dans l'analyse numérique. Un domaine important de l'imagerie tridimensionnelle est la médecine et plus précisément l'examen d'organes vivants. Des mouvements du patient sont inévitables, quand on ne doit pas considérer un organe mobile comme le coeur dont les battements s'accompagnent à la fois de déplacements et de déformations. Le problème consiste à obtenir une image cohérente d'une série de mesures qui, dans les procédés ordinaires, sont accomplies successivement en déplaçant le réseau de capteurs autour de l'objet. Une solution simple consiste à disposer un grand nombre de réseaux autour du patient, immobiles et qui prennent simultanément toutes les images utilisées pour la reconstruction. On élude alors le problème de mobilité de l'organe considéré, au prix d'un enchérissement considérable du dispositif utilisé. Une autre solution fréquemment utilisée consiste à ne prendre des mesures qu'à des instants en phase pendant des cycles différents de l'organe, mais cela n'est possible que pour des phénomènes périodiques et nécessite une certaine précision de synchronisation, B 12254.3 JCI tout en allongeant la durée de prise des mesures. Enfin, on peut négliger purement et simplement les mouvements de l'organe et se contenter d'une image brouillée. Mais il est évident que cette façon de procéder possède comme les autres des inconvénients. Citons en passant le brevet français n 2 701 135 ayant un inventeur commun avec celui-ci et qui fournit un procédé de reconstruction d'image d'un objet évoluant d'une prise de projection à l'autre, dans lequel la loi d'évolution est connue au moins approximativement et est incorporée aux calculs pour donner finalement l'image de l'objet à un instant quelconque. C'est le même résultat qui est obtenu ici, mais dans des conditions différentes puisque la loi d'évolution est supposée inconnue et doit être retrouvée. De plus, le brevet antérieur est relatif aux méthodes algébriques d'inversion, alors que celui-ci ressortit aux méthodes analytiques. Et ici, plusieurs jeux de projections sont pris à chaque état, à des incidences différentes autour de l'objet. L'invention a pour objet d'obtenir des résultats satisfaisants d'imagerie sans complication de l'appareil, par un procédé impliquant, dans la réalisation choisie, des itérations améliorant un modèle de l'image en réduisant l'effet des dégradations liées aux conditions d'acquisition (bruitages, artéfacts, etc.) pour reconstruire une image brute (approchée car fortement bruitée en pratique) associée à chaque jeu de mesures et consistant à déterminer les déplacements et les déformations de l'objet pour ce jeu de mesures afin d'inclure les résultats des différentes images approchées dans le modèle. Pour résumer, il s'agit d'un procédé comprenant une prise d'une série de jeux de mesures d'une propriété caractéristique de l'objet par un B 12254.3 JCI réseau bidimensionnel de capteurs mobiles autour de l'objet, les jeux de mesures consistant en un réseau bidimensionnel de projections de la propriété sur les capteurs, les jeux de mesures étant pris à un groupe d'états différents de l'objet en ce qui concerne des variations de la position ou de forme de l'objet, consistant à reconstruire une image brute de l'objet à chacun des états en combinant leurs jeux de mesures respectifs, caractérisé en ce qu'il consiste ensuite à évaluer les variations de l'objet entre une image de l'objet à un état de référence et les images approchées, et à modifier l'image de l'objet à l'état de référence par une combinaison avec les jeux de mesures et les variations de l'objet associées aux jeux de mesures. Avantageusement, les variations de l'objet sont estimées par un champ de déplacement de points de modélisation de l'objet et par une évolution de la propriété respectant une relation de continuité de la propriété dans le temps; elles peuvent être estimées par des variations de contours de l'objet modélisés par des formes géométriques simple et de nature invariable; de plus, si l'objet est sujet à la fois à des déformations et à des déplacements d'ensemble, on peut décider que les déformations et les déplacements d'ensemble sont évalués séparément avec des modélisations différentes de l'objet. Un autre aspect de l'invention consiste à combiner les jeux de mesure en tenant compte de bruitages des mesures et de bruitages des évaluations des variations de l'objet. Les figures suivantes seront désormais commentées pour une description détaillée d'une réalisation préférée de l'invention: * la figure 1 illustre la prise des mesures, B 12254.3 JCI * la figure 2 résume le procédé, * la figure 3 illustre une première modélisation d'un myocarde, * et les figures 4 et 5 illustrent une seconde modélisation d'un myocarde. Sur la figure 1, on a illustré un objet d'étude fréquent, à savoir un coeur 1 et plus précisément son myocarde, qui a été alimenté par un liquide enrichi d'un traceur radioactif émetteur de photons. Les mesures sont prises par un réseau bidimensionnel 2 de capteurs 3, qui sont collimatés dans la même direction. Ils recueillent donc les photons émis sur des lignes de projection p respectives toutes parallèles entre elles. L'ensemble des projections s'étend dans un faisceau 4 de section constante dans lequel le coeur 1 est plongé. Le réseau bidimensionnel 2 est mobile autour du coeur 1 en glissant sur un rail 5 circulaire ou d'une autre forme, et qui peut d'ailleurs basculer autour de pivots 6 pour prendre une inclinaison quelconque. Les capteurs 3 ainsi placés à des incidences différentes prennent une série de jeux de mesures en étant réceptifs aux photons à certains instants, synchronisés sur un groupe déterminés d'états du coeur 1, ou de phases de son cycle. Les moyens de synchronisation, de commande des capteurs 3, de déplacement du réseau bidimensionnel 2 et de mémorisation des mesures, déjà connus, ne seront pas décrits ici car l'invention n'en propose pas de nouveau. Le faisceau à section constante est usuel avec des procédés de détection de photons émis. L'invention est encore applicable à des procédés caractérisés par des faisceaux différents, coniques ou en éventail: seuls les coefficients de la matrice de projection R varient. B 12254.3 JCI Voici une explication détaillée du procédé utilisé. La série de jeux de mesures successifs correspond à divers états de l'objet dont le numéro d'ordre est représenté par l'indice k (entier positif), o k=0 désigne l'état dans lequel l'image doit être reconstruite. L'objet peut être par exemple discrétisé en M3 points formant un cube de M points de côté. L'image de l'objet à chaque état peut être représenté par un vecteurcolonne à M2 composantes, l'opérateur de déformation de l'objet entre l'état 0 et l'état k par une matrice carrée de dimension M3 notée Dk, et les mesures peuvent être représentées par un vecteur- colonne de projections mk à LxM2 composantes. En effet, on suppose qu'il y a M2 capteurs qui recueillent autant de projections passant par une ligne de M points de modélisation de l'objet et que L jeux de mesures ont pu être pris pour le même état de l'objet, à des incidences différentes, en synchronisant les jeux de mesures sur les cycles de l'objet. Comme les projections sont affectées par un bruit bk et comme on appelle Wk un bruit qui représente l'erreur de l'opérateur de déformation, on peut écrire la paire d'équations suivante: { =Dkf0 + wk {fk DkfO + wk (2) tmk =Rfk + bk '(2 en appelant fk et fo les images réelles de l'objet aux états k et 0, et R la matrice de projection des mesures. Selon l'invention, on fait entrer les jeux de mesures les uns après les autres dans les calculs pour reconstruire à chaque fois une estimation de ^k+l l'image, notée f0+1 quand l'état k+1 a été B 12254.3 JCI considéré, et tenter de l'améliorer et l'approcher de f0. Concrètement, on applique l'équation: -k+l -k f+l = f0+ Kk+lk+1. (3) Cette formulation du problème fait apparaître plusieurs nouvelles grandeurs définies de la façon suivante: Kk+l = PkD+lRtP-1 1+ Sk+l ' avec Pk+1 = Pk(I - Dt Rtp-k+1 RDk+1Pk) (4) k+l F-k+l et Pk+l = RDk+iPkD+iRt + RCwk+lRt + Cbk+1; et Sk+1 = mk+1 - RDk+lfk Pk désigne une matrice d'estimation de l'erreur entre l'estimation d'image obtenue k+1 et l'image réelle fo. Pk calculé par récurrence pour chaque jeu de projections à partir d'une valeur primitive P-1 introduite arbitrairement par l'opérateur et qui est en pratique formée d'une matrice dont les coefficients sont tous très grands. PFEk+l désigne la matrice exprimant le bruit entre l'estimation d'image précédente fk (une fois déformée à l'état k+1 par l'opérateur Dk+1 et projetée) et le jeu de projections mk+1 correspondant à l'état k+l. Les matrices Cbk et Cwk désignent les matrice de covariance des bruits bk sur les projections et wk sur l'opérateur de déformation comme ils ont été définis ci-dessus. Ces bruits ne sont pas connus en pratique et doivent être estimée par l'opérateur d'après son expérience et ses observations. Le bruit bk peut être évalué sur chaque image de projections par une analyse de densité spectrale en considérant les hautes fréquences et en B 12254.3 JCI supposant le bruit blanc. On rappelle enfin que la matrice de covariance d'un phénomène X de moyenne X (exprimé sous forme matricielle) est définie par C(X) = E(X - X)t.(X - X) (5), o E désigne l'espérance statistique. Enfin, la lettre t en indice renvoie à la transposée de matrice. Ce procédé revient à appliquer à chaque récurrence une régularisation au sens de Tikhonov tel que -k+1 minimise la fonction de coût k+l +l k E(f) = A!(RDk+lf, mk+l) + 2 (f') (6) o: k+l(x Y) = (x - y)t(RCwK+Rt + Cbk+l) (x - y) Ak+l(x y) = (x - y)tP1(x - y) (7) 2 k o A1 et A2 sont des normes qui expriment l'écart entre les mesures et la solution trouvée et l'écart entre la solution f0 et l'image reconstruite f. L'accomplissement de tous ces calculs impose de trouver l'opérateur de déformation Dk+l. On utilise pour cela l'estimation d'image précédemment obtenue fk et une image brute îk+1 de l'objet à l'état k+l à l'aide des projections mk+l qui viennent d'être prises en considération, en utilisant un algorithme usuel de reconstruction d'image tridimensionnelle dit de rétroprojection des projections filtrées. On rappelle qu'on dispose de L jeux de M2 mesures pour chaque état, ce qui permet de reconstituer l'image brute îk+l. La comparaison entre l'image brute à l'état k+l fk+l dans laquelle on B 12254.3 JCI discerne la forme et la position de l'objet recherché et l'estimation d'image actuelle fI de l'objet à l'état 0 avec les jeux de mesures des états 0 à k permet de déduire cet opérateur de déformation Dk+l. Cette étape, qui est une des plus caractéristiques de l'invention, est décrite dans les lignes qui viennent. Ensuite, l'estimation d'image suivante de l'objet f0k+1 est reconstruite à partir de l'estimation d'image actuelle f, du jeu de projections mk+l et de l'opérateur de déformation Dk+l en appliquant les équations précédentes. On se rapproche ainsi pas à pas vers l'image réelle fo, pour obtenir une image finale fO représentative de l'état 0. La figure 2 résume le procédé. On va donc examiner plus concrètement comment définir l'opérateur des déformations Dk dans le cas particulier de l'examen du myocarde. Comme on le verra, le choix de l'objet ou de l'organe tend à imposer une schématisation de sa forme pour le modéliser et mener les calculs commodément, mais d'autres modélisations seront immédiatement envisagées pour d'autres organes. Dans le cas d'un procédé (TEMP) de tomographie d'émission monophotonique, généralement employé pour étudier le fonctionnement d'un organe et o un produit radioactif est administré à cet organe, le procédé classique pour l'isoler sur l'image brute, reconstituée directement, ik+l consiste tout d'abord à sélectionner manuellement une portion de l'image qui l'englobe, puis à écarter les parties étrangères de l'image qui subsistent par un critère de seuil d'activité: plus précisément, on écarte les parties dont l'activité est inférieure à une proportion, généralement de 35%, de l'endroit d'activité maximale. Il peut subsister alors, à côté de l'organe lui-même, B 12254.3 JCI des objets résiduels de plus petite taille qu'on écarte ensuite, ce qui élimine les fausses détections associées au bruit et aux artéfacts de reconstruction. On peut ainsi éliminer automatiquement les objets dont la taille, mesurée d'après leur surface sur l'image, est inférieure à 30% de la taille de l'objet principal. Dans le cas o on peut considérer que l'activité d'émission globale par image reste constante au cours du temps, l'équation de continuité: a + V. (a v) = 0, Et fréquemment rencontrée en mécanique des milieux continus et mécanique des fluides, est respectée, a(M,t) et v(M,t) désignant la densité d'activité d'émission et la vitesse de l'objet en un point M et à un instant t. V désigne l'opérateur de divergence. On appelle Vk le champ des mouvements vk(M) des points M de l'objet de l'état d'origine (k=0) à l'état k. Il est possible de décomposer ce mouvement par l'équation: v(M) = v (C) + A C M + D(M) (8) o C est le centre d'inertie de l'objet (le myocarde), v(C) et Q les composantes de translation et de rotation du torseur du champ de mouvement rigide de -4 l'objet et d(M) le mouvement non rigide exprimé dans un référentiel lié à l'objet, c'est-à-dire les déformations proprement dites. Dans le cas du myocarde, ce référentiel a pour origine et pour axes directeurs le centre C et les axes principaux d'un ellipsoïde servant à modéliser la B 12254.3 JCI surface barycentrique du myocarde entre ses surfaces interne et externe (endocarde et péricarde). Cette surface peut être définie en coordonnées sphérique pour la pointe du coeur (appelée apex) et cylindrique pour la base du coeur; ce système de coordonnées est noté (p, p,.) par l'équation: pmaxa(p, G,.)pdp pG(%G,.) = (9) pmaxa(p,%G,.)dp o a(p,4G,.) désigne l'activité d'émission au point de coordonnée (p,4G,.) et pmax désigne la distance maximale de scrutation. Pour toute coordonnée angulaire (G,.), le point de la surface barycentrique a donc comme coordonnée radiale pG. Cette méthode semble plus exacte que d'autres si le coeur présente un défaut de vascularisation et en particulier que les méthodes fondées sur la recherche des moments principaux d'inertie. En effet, les défauts d'activités qui apparaissent alors peuvent plus facilement perturber les résultats des calculs. On a représenté sur la figure un myocarde vu en coupe, l'ellipsoïde apte à modéliser la surface barycentrique et le myocarde pour les besoins d'estimation de son mouvement d'ensemble, et l'échantillonnage retenu, qui est hémisphérique du côté de la pointe du myocarde et cylindrique du côté de sa base. Pour cette partie, les rayons de détermination de la surface barycentrique G sont tous parallèles. Ces rayons le long desquels les points Gi sont déterminés sont tous notés par r. On commence donc par calculer les points Gi de coordonnées xi,Yi et zi de la surface barycentrique B 12254.3 JCI G à l'aide de la formule donnée ci-dessus, puis on calcule les paramètres de l'ellipsoïde 7 le plus proche de cette surface en définissant la fonction d'erreur: N 2 2 2 F(x,y,z)=aox +aly +a2z +a3xy+a4xz+asyz+a6x+a7y+a8z=1 (11). On résout le système obtenu pour obtenir les coefficients ao à a8 jusqu'à obtenir le minimum de cette fonction d'erreur; un nouveau référentiel est associé à l'ellipsoïde obtenu, qui implique d'autres rayons de scrutation et de calcul des coordonnées des points Gi; en pratique, on recommence le procédé par itérations successives jusqu'à ce que les variations de position et d'orientation des axes du référentiel associé à l'ellipsoïde 7, et qui découlent des coefficients ao à a8, deviennent insensibles compte tenu de la résolution de l'image. Pour la première itération l'ellipsoïde 7 et le référentiel associé sont obtenus en résolvant l'équation d'erreur avec tous les points de l'image du myocarde, en l'absence de rayons permettant de calculer une surface barycentrique. Les axes dessinés sur la figure sont alors obtenus. L'axe e3 est toujours dirigé dans le grand axe de l'ellipsoïde vers la pointe du myocarde. L'origine C est le centre de l'ellipsoïde. Une fois que le référentiel a été calculé pour l'image k, le mouvement global de translation B 12254.3 JCI entre les états 0 et k est défini par le vecteur cock et la matrice de rotation est calculée par le produit des composantes des deux vecteurs [ k -k kt r. o e1, e2, e3 j. [2 t e2, 3] Le mouvement non rigide peut être estimé à partir des contours du myocarde et de la surface barycentrique G. Une modélisation moins schématique du myocarde est désormais utilisée: l'endocarde et le péricarde sont modélisés par des surfaces en balle de fusil, c'est-à-dire demi-ellipsoidales de révolution à l'extrémité et cylindriques de révolution à la base. La position de ces surfaces est effectuée en recherchant les paramètres de ces surfaces qui permettent de maximiser l'intégrale sur ces surfaces de la norme quadratique du gradient de la propriété d'émission. Cela provient de ce que, comme seul l'organe émet les photons, la variation d'émission est la plus forte en traversant ses surfaces interne et externe, et on cherche donc les surfaces de forme simple qui vérifient au mieux cette propriété. Les paramètres nécessaires à définir les deux surfaces en balle de fusil sont au nombre de cinq (la hauteur des cylindres et les mesures des ellipsoïdes aux axes principaux). Si on appelle a la demi-mesure des ellipsoïdes selon leurs petits axes (perpendiculairement à e3), c leur demi- mesure selon leur grand axe e3, A l'indice de l'endocarde, B l'indice du péricarde et h la hauteur des cylindres, les calculs sont menés en faisant varier aA, aB, cB, et on établit par convention cA=l, 5aA car on constate que ce paramètre est difficile à calculer. De même, le paramètre h ne peut pas être obtenu facilement directement. Pour cela, l'image brute îk+1 est B 12254.3 JCI décomposée en coupes passant par l'axe e3, o p est constant, et on mesure pour chacune de ces coupes le gradient de l'activité le long de la surface comprise entre les courbes SA et SB modélisant l'endocarde et le péricarde. Le plan Q auquel le cylindre s'arrête correspond au gradient maximal de l'activité, sommée sur toutes les coupes. Tout cela est représenté à la figure 4. Il est alors possible de définir complètement une modélisation SG de la surface barycentrique G qui est une surface en balle de fusil définie par les équations: F(x, y,z) = a0x2 + aly2 + a2z2 si z > 0 SG 2 2. - (12) F[r(x,y,z)= aox2 + alY2 si - h A chaque point M du myocarde, on peut faire correspondre une surface en forme de balle de fusil intermédiaire à SA et SB et définie par de nouveaux paramètres a et c; h restant identique. Après une déformation de dilatation ou de contraction, le point M, de coordonnées x, y et z, devient un point M' de coordonnées x', y' et z' présent sur une surface de même forme et définie par les paramètres a', c' et h' tels que: x'=(a'/a)x, y'=(a'/a)y et z'=[(c'+h')/(c+h)]z. (13) En modélisant -comme illustré à la figure - le myocarde déformé comme une superposition de tranches en forme de balles de fusil dilatées inégalement, les déplacements de déformation du point M entre les états O et k+1 suivant les trois axes peuvent être exprimés par les formules suivantes: B 12254.3 JCI d(x)= - 1 x d(M) d(y) = tb _ 1ly (14) d ____ (14> d(z) = cl +hl 1) l ( c + h) De façon plus précise, les coefficients a, b et c associés au point M obéissent à des relations de proportionnalité avec les coefficients de même nom de la surface barycentrique d'une part, et de la surface modélisant soit l'endocarde, soit le péricarde d'autre part. Les coefficients a', b' et c' associés au pont M' obéissent aux mêmes relations avec les coefficients des mêmes surfaces de modélisation à l'état déformé. Il faut ajouter à ces déformations les déplacements d'ensemble pour obtenir le mouvement complet en chaque point du myocarde. L'étape suivant consiste à définir l'opération des déformations Dk+1 entre l'état 0 et l'état k+1. Aucun état intermédiaire n'intervient dans le calcul, ce qui réduit les erreurs numériques. L'activité à l'état k+1 du point M(x,y,z), soit a(M,k+1), est égale à la somme des activités des points Mi à l'état 0, a(Mi,0) qui arrivent au voisinage de M après un mouvement complet de coordonnées vx(M), vy(M) et vz(M) de M, avec une pondération selon la formule: a(M,k + 1) = aia(Mi,0) et ai = (1 - Six)(l - giy)(l - siz), six = x - (Xi + Vx (Mi)si|x -(Xi + vx(Mi)I 1; O sinon O iy = y - (Yi + vy (Mi) si y - (Yi + vy (Mi) B 12254.3 JCI Il faut remarquer que l'activité globale de l'objet est pratiquement identique pour chaque état. En effet, même si elle diminue à chaque mesure à cause de l'épuisement de la radioactivité, les états sont tous mesurés plusieurs fois à des instants successifs, et les mesures des différents états sont entremêlées si bien que certaines sont précoces et d'autres sont tardives pour chacun d'eux. C'est en ce sens que l'équation de continuité est respectée dans ce procédé. D'autres modes d'obtention des mouvements des points de l'objet sont concevables, bien que celui- ci paraisse actuellement le meilleur. Un autre est décrit dans l'article "Computation of 3-D velocity fields from 3-D cine CT images of a human heart" par Song et Leahy, paru dans IEEE transactions on medical imaging, vol. 10, n 3, septembre 1991, pp. 295-306. Le calcul de l'image f0k+1 implique l'inversion de la matrice Pzk+l, conformément à la formule rencontrée plus haut qui donnait Kk+l. Cette matrice de très grande taille ne peut pas être inversée directement en général. Un autre procédé est alors utilisé. On montre que, sous les hypothèses précédentes, (fo- fk+1) et ck+1 sont des vecteurs o aléatoires de moyenne nulle et dont les matrices de covariance sont Pk+l et PEk+l. Ces matrices de covariance peuvent être approchées par des matrices- blocs circulantes qui peuvent être diagonalisées selon les formules suivantes: PCk+l = WmPEk+1Wm (16) Pk+l = VOk+lW10 Vm est une matrice de dimensions LM2xLM2 telle que: B 12254.3 JCI Wm... Vm. : (17) 0...Wm Wm est une matrice carrée de dimensions M2xM2 telle que si X représente une image bidimensionnelle à M2 points, le produit W-1X est la m transformée de Fourier discrète bidimensionnelle de X; fk+l est une matrice-blocs diagonale dont chaque bloc est une matrice diagonale de dimensions M2xM2 dont les éléments diagonaux sont les densités spectrales du bruit entre les projections et le calcul de celles-ci à partir des déformations et de l'estimation d'image, de valeur (mk+1 - RDk+lf); Wo est une matrice de dimensions M xM telle que si Xreprésente une image tridimensionnelle de M3 points, W01X est la transformée de Fourier discrète tridimensionnelle de Xi; et Pk+l est une matrice diagonale dont les éléments diagonaux sont les densités spectrales de l'erreur d'estimation (fo - fk+1) Dans ces conditions, le produit Kk+lsk+l t0à W IlDt RtWmP-1 WlEk + 1 (18). est à WOPkW0k + 1RWmk + 1mk + 1 Comme le myocarde est à peu près incompressible, le produit DkDk est proche de la matrice d'identité. On obtient finalement, pour un faisceau de projections parallèles, l'équation suivante k+l _ k à partir de f0k+ = f + Kk+1Sk+l, en exploitant les propriétés de l'opération de projection R et des opérateurs Wm et Wo: fk + 1 = + DkL+ 1RtDk + l(mk + 1 - RDk + lf) (19) o Rt et Dk+l désignent respectivement une rétroprojection et un banc de L filtres B 12254.3 JCI monodimensionnels dont chacun est défini par une fonction de transfert indépendante de l'angle de projection noté 1; cette fonction est définie par: l + (V) r Fk (v) (20) k + 1 Fk + 1 (v)' o rk(v) et rFSk+l(V) désignent les densités spectrales de Pk et ftk+l; de plus, rk+l(v) = rk(v){l- kt (v) (21) IvISsk+l (v" 2) et rFEk+l(v) = k(v) + wk+(v) + rbk+1(v), Fwk+1 et Fbk+1 désignant les densités spectrales de Wk+1 et bk+1. En pratique et compte tenu des propriétés de l'opération de projection R et de symétrie de révolution de (D + (v), l'équation précédente revient k+, à ok+1 = (1 - 'k+l)fk + Dk-1Rt k+lmk+1, k+1/v + x y) avec qyk+l(Vx,vy,vz) v + v (22); x Y x y et le procédé effectif de reconstruction est accompli en résolvant cette dernière équation. En résumé, à chaque itération k+l1, les mesures mk+l sont filtrées par le filtre de fonction de transfert (i + traduisant la confiance en ces données B 12254.3 JCI et en l'estimation des déformations, puis rétroprojetées et déformées par l'opérateur Dk 1 afin d'être recalées sur l'état 0 dont l'image est restaurée. Cette image f + 1 est alors la somme des mesures filtrées reconstruites et déformées et du modèle filtré par le filtre de fonction de transfert (1-wk+l) opérant une pondération entre les mesures et les informations a priori. Les filtres À et y sont actualisés en fonction des bruits des jeux de mesures et des modélisations des déformations. Le processus est initialisé avec f1 = 0 et P1 = a2I o a-1 est un très grand nombre, assimilable à l'infini. L'état de référence correspondant à la solution fo est librement choisi. Il peut correspondre à un état effectif de prises de mesures ou à un état conventionnellement choisi tel qu'un état primitif défini par une modélisation et inaccessible désormais. Il est d'ailleurs possible de changer d'avis et de remplacer finalement l'image obtenue définitivement fO en se servant encore de l'opérateur de déformations Dk. Les essais préliminaires ont démontré que la combinaison d'informations d'images de l'objet à des états différents, caractéristique de l'invention, donnait des images beaucoup plus nettes que les reconstructions classiques o seules des mesures synchronisées sur un état unique de l'objet étaient prises; en d'autres termes, l'image f est meilleure que fo (ou fo), d'autant plus que k est plus grand. Mais le nombre d'états peut être quelconque; l'invention peut même être employée avec profit pour un seul état de prises de mesures, à condition de disposer d'informations de départ associées à un état de référence; l'image reconstruite, enrichie de ces informations, peut B 12254.3 JCI avantageusement être rendue dans l'état de prise de mesures, car elle est normalement plus précise que l'image brute à. Une application différente de la médecine, à laquelle s'appliquent volontiers certaines des remarques précédentes, concerne les contrôles non destructifs de l'industrie o les états de prises de mesures sont des états déformés et l'état de référence un état défini par des plans ou une modélisation, indemne de déformation. B 12254.3 JCI REVENDICATIONS 1. Procédé de reconstruction d'images tridimensionnelles d'un objet comprenant une prise d'une série de jeux de mesures (m) d'une propriété caractéristique de l'objet (1) par un réseau bidimensionnel (2) de capteurs (3) mobile autour de l'objet, les jeux de mesures (m) consistant en un réseau bidimensionnel de projections (p) de la propriété sur les capteurs, les jeux de mesures étant pris à un groupe d'états (k) différents de l'objet en ce qui concerne des variations de position ou de forme de l'objet, consistant à reconstruire une image brute (f) de l'objet à chacun des états en combinant leurs jeux de mesures respectifs, caractérisé en ce qu'il consiste ensuite à évaluer (Dk) les variations de l'objet entre une image de l'objet à un état de référence et les images brutes, et à modifier l'image de l'objet à l'état de référence (fO) par une combinaison avec les jeux de mesures (mk) et les variations (Dk) de l'objet associées aux jeux de mesures. 2. Procédé de reconstruction d'images tridimensionnelles d'un objet selon la revendication 1, caractérisé en ce que en ce que les variations de l'objet sont estimées par un champ de déplacement de points de modélisation (M) de l'objet et par une évolution de la propriété respectant une relation de continuité de la propriété dans le temps. 3. Procédé de reconstruction d'images tridimensionnelles d'un objet selon l'une quelconque des revendications 1 ou 2, caractérisé en ce que les variations de l'objet sont estimées par des variations de contours de l'objet modélisés par des formes géométriques simples (SA, SB) et de nature invariable. B 12254.3 JCI 4. Procédé de reconstruction d'images tridimensionnelles d'un objet selon la revendication 3, caractérisé en ce que l'objet est modélisé sous forme de points (M) présents sur des formes géométriques simples (SM) et de nature invariable et semblable à celle des formes géométriques de modélisation des contours. 5. Procédé de reconstruction d'images tridimensionnelles d'un objet selon l'une quelconque des revendications 3 ou 4, caractérisé en ce que l'objet est un myocarde (1) modélisé par une demi- sphère prolongée par un cylindre à l'état de référence et par une empilement de tranches de modélisation composées chacun d'un demiellipsoïde prolongé par un cylindre aux états des images brutes. 6. Procédé de reconstruction d'images tridimensionnelles d'un objet selon l'une quelconque des revendications 1 à 5, caractérisé en ce que l'objet est sujet à la fois à des déformations et à des déplacements d'ensemble, et en ce que les déformations et les déplacements d'ensemble sont évalués séparément avec des modélisations différentes de l'objet. 7. Procédé de reconstruction d'images tridimensionnelles d'un objet selon l'une quelconque des revendications 1 à 6, caractérisé en ce que les jeux de mesure sont combinées successivement, pour chacun des états, avec l'image de l'objet à l'état de référence (f0). 8. Procédé de reconstruction d'images tridimensionnelles d'un objet selon l'une quelconque des revendications 1 à 7, caractérisé en ce que les jeux de mesure sont combiné en tenant compte de bruitages (b, w) des mesures et de bruitages des évaluations des variations de l'objet. B 12254.3 JCI 9. Procédé de reconstruction d'images tridimensionnelles d'un objet selon l'une quelconque des revendications 1 à 8, caractérisé en ce que l'objet est supposé incompressible.5 B 12254.3 JCI