FR2797978A1 - Procede de recalage automatique d'images - Google Patents

Abstract

Procédé de recalage automatique d'images d'angiographie bi- et tridimensionnelles par comparaison de l'image d'angiographie numérique soustraite bidimensionnelle avec des données relatives à une image tridimensionnelle reconstruite à partir de séquences d'angiographie rotationnelle, dans lequel on estime un champ de distorsions dans l'image, on estime une matrice de projection conique et on effectue une approximation d'une transformation rigide dans l'espace égale à la différence entre un recalage initial basé sur le champ de distorsions et sur la matrice de projection conique et un recalage parfait.

Description

Procédé de recalage automatique d'images.

La présente invention concerne le domaine du traitement

d'images, en particulier d'images radiologiques bi- et tridimensionnelles.

De façon connue, les appareils de radiologie comprennent un moyen d'émission d'un faisceau de rayons X tel qu'un tube à rayons X et un moyen de réception dudit faisceau tel qu'un détecteur à l'état solide ou

encore un scintillateur et une caméra vidéo, par exemple de type CCD.

Le moyen d'émission et le moyen de réception de rayons X sont généralement supportés par un système mobile à un ou plusieurs axes, permettant des prises de vue selon différentes incidences. Le moyen de réception du faisceau de rayons X est relié à des moyens de traitement d'images permettant la génération d'images tridimensionnelles à partir

d'une série d'images bidimensionnelles prises par le moyen de réception.

Il s'avère important d'obtenir une bonne correspondance entre une image tridimensionnelle et une image bidimensionnelle prise lors d'une étape au cours de laquelle le patient subit une action particulière telle que la pose d'un cathéter dans le domaine de l'angiographie, pour pouvoir suivre le déplacement du cathéter dans l'image bidimensionnelle

mais également dans l'image tridimensionnelle.

On connaît, par l'article "Fully automatic 3D/2D Subtracted Angiography Registration" de Messieurs KERRIEN, LAUNAY, BERGER, MAURINCOMME, VAILLANT et PICARD, publié dans CARS' 99, Proceedings of the 13th International Congress and Exhibition, Paris, 23-26 juin 1999, un procédé de recalage automatique d'images d'angiographie soustraites. Cet article indique que des images tridimensionnelles reconstruites de la vascularisation cérébrale peuvent être obtenues à partir d'un grand nombre d'images d'angiographie soustraites et qu'en neurologie interventionnelle, il est intéressant de recaler le volume tridimensionnel reconstruit sur les images

bidimensionnelles instantanées.

La présente invention a pour objet un procédé de recalage perfectionné. La présente invention a également pour objet un procédé de recalage de précision millimétrique ou submillimétrique à temps de calcul faible. Le procédé de recalage automatique, selon l'invention, est destiné à des images d'angiographie bi- et tridimensionnelles par comparaison d'une image d'angiographie numérique soustraite bidimensionnelle avec des données relatives à une image tridimensionnelle reconstruite à partir de séquences d'angiographie rotationnelle, dans lequel on estime un champ de distorsions dans l'image, on estime une matrice de projection conique et on effectue une approximation d'une transformation rigide dans l'espace égale à la différence entre un recalage initial basé sur le champ de distorsions et sur la matrice de projection conique et un recalage parfait. La transformation rigide comprenant une partie translation et une partie rotation, on effectue une approximation de la partie translation en considérant que la partie rotation est connue, la transformation optimale étant atteinte lorsque la corrélation de l'image d'angiographie numérique soustraite et des données

relatives à l'image tridimensionnelle reconstruite est maximale.

De préférence, on suppose que l'erreur résiduelle de recalage est due à une erreur de positionnement en translation et en rotation de l'image tridimensionnelle reconstruite dans le repère d'un moyen de prise d'images. De préférence, la différence entre les dites données relatives à l'image tridimensionnelle reconstruite et l'image angiographique numérique soustraite bidimensionnelle est considérée comme étant due à un petit déplacement suivant une transformation rigide, on détermine la

dite transformation rigide par une technique de flux optique modifiée.

De préférence, les dites données relatives à l'image tridimensionnelle reconstruite sont une image bidimensionnelle obtenue par projection conique d'intensité maximale de l'image tridimensionnelle reconstruite. Ce type de projection associe un unique voxel à chaque pixel. On peut utiliser pour chaque pixel, lors de la projection, la profondeur du voxel élu comme valeur d'une des coordonnées spatiales à un instant donné. Avantageusement, lors de chaque déplacement de l'image d'angiographie numérique soustraite, la nouvelle image bidimensionnelle est calculée à partir de l'ancienne image bidimensionnelle et ce, sans

nouveau calcul de projection de l'image tridimensionnelle reconstruite.

On réduit considérablement la durée de calcul nécessaire, de plusieurs

heures à une ou deux minutes.

Dans un mode de réalisation de l'invention, le recalage initial est

donné par un ensemble marqueur de position tridimensionnelle connue.

Dans un autre mode de réalisation de l'invention, le recalage

initial est donné par étalonnage de la machine d'angiographie.

Ainsi à partir d'un recalage initial, on réduit l'erreur initiale en translation par maximisation de la corrélation entre l'image d'angiographie numérique soustraite et les données relatives à l'image tridimensionnelle reconstruite, et par affinement du recalage par une technique de flux optique modifié, on réussit à obtenir un recalage de précision millimétrique ou infra millimétrique parfaitement suffisant pour une utilisation réelle au cours d'une opération, notamment de neurochirurgie ou de radiologie interventionnelle, au bout d'une durée de

calcul limitée à quelques minutes, voire à quelques dizaines de secondes.

L'invention s'applique particulièrement au cas de la radiologie avec des images tridimensionnelles o seuls les vaisseaux sanguins sont

visibles. Le contenu informatif de telles images est assez faible.

La présente invention sera mieux comprise à l'étude de la

description détaillée d'un mode de réalisation pris à titre d'exemple

nullement limitatif.

En radiologie interventionnelle, il est souhaitable que le praticien puisse savoir à tout instant o son cathéter se trouve dans le corps du patient avec la meilleure précision possible. Il peut à l'heure actuelle déduire cette information d'images d'angiographie numériques soustraites dites DSA, qu'il relie à des images préopératoires d'imagerie par résonance magnétique grâce à sa connaissance anatomique. Les images DSA offrent entre autres une imagerie temps réel à très haute résolution spatiale. Depuis peu, on se 'sert également de reconstructions tridimensionnelles des vaisseaux sanguins dites "3DXA" à partir de séquences d'angiographie rotationnelles effectuées par rotation rapide du tube à rayons X et du moyen de prises d'images sur un demi-tour et prise d'une cinquantaine d'images DSA qui sont les projections en entrée d'un algorithme de tomographie qui produit en sortie l'image 3DXA. Pour plus de précision sur cette technique, le lecteur est invité à se reporter à la thèse de Monsieur LAUNAY "Localisation et reconstruction 3D à partir d'angiographies stéréotaxiques " thèse de doctorat, Institut National

Polytechnique de Lorraine, Nancy, France, 1996.

Ces reconstructions permettent une très bonne appréciation de l'angioarchitecture. De plus, ces images tridimensionnelles peuvent être utilisées en temps réel selon plusieurs types de visualisation, tels que la projection d'intensité maximale, l'isosurface, la fonte de volume, l'endoscopie virtuelle ou encore la coupe reformatée, et apportent une

aide complémentaire aux diagnostics des praticiens.

L'invention permet une mise en correspondance entre les images

DSA et les images 3DXA.

On entend ici par "image DSA", l'image de maximum d'opacification jusqu'au rang N dans la séquence acquise, c'est-à-dire que chaque pixel de l'image résultat prend la plus petite valeur rencontrée sur les N premières images de la séquence, ou l'image de rang N dans la séquence acquise. Le rang N de l'image est soit choisi par l'utilisateur, soit

fixé par rapport à la cadence d'acquisition.

Par ailleurs, le volume 3DXA est projeté en vue dite "MIP" (projection d'intensité maximum) conique. Chaque voxel du volume 3DXA peut être considéré comme une sphère de rayon égal à la taille du voxel, et est projeté dans le plan image. La valeur finale de chaque pixel est la valeur maximale projetée en ce lieu. L'image ainsi obtenue est désignée par le terme "image MIP". En variante, on peut utiliser une

projection conique par addition des valeurs de voxels en chaque pixel.

Le recalage entre les deux images s'effectue par estimation d'un champ de distorsions', de l'image, technique connue par l'article "Machine precision assessment for 3D/2D digital subtracted angiography images registration" de Messieurs KERRIEN, VAILLANT, LAUNAY, BERGER, MAURINCOMME et PICARD, In SPIE medical Imaging, volume 3338, pages 39-49, 1998, auquel le lecteur est invité à se reporter pour plus de

précisions sur cette technique.

Après l'estimation du champ de distorsions dans l'image, on effectue une estimation d'une matrice de projection conique. Un étalonnage permet de connaître les paramètres intrinsèques de cette matrice, les paramètres extrinsèques étant déterminés comme expliqué

plus loin.

La machine de radiologie, une fois étalonnée, fournit un recalage initial qui diffère du recalage recherché ou recalage parfait par une transformation rigide (rotation + translation) dans l'espace tridimensionnel. Cette transformation rigide possède une partie rotation de faible amplitude alors que sa partie translation peut être considérée

comme inconnue.

Dans un premier temps, on retrouve la partie translation en considérant que la partie rotation de la matrice de projection est parfaitement connue. La position optimale est atteinte quand la corrélation normalisée centrée de l'image DSA avec l'image MIP est maximale. On suppose que l'erreur résiduelle de recalage est due à une petite erreur de positionnement, à la fois en rotation et en translation du volume 3DXA, dans le repère de la caméra correspondant au recalage obtenu. On considère alors que l'image MIP est prise à un instant t1, puis que le volume 3DXA bouge suivant une petite transformation rigide. Une image est acquise à l'instant t2 et on obtient l'image DSA. Le mouvement apparent que l'on observe entre l'image MIP actuelle et l'image DSA donne des informations concernant la petite transformation rigide recherchée. Ces informations sont extraites par une technique de flux

optique modifiée.

La corrélation normalisée centrée entre l'image DSA, Id, et l'image MIP, Im, est définie par (x,y)V(Id(x, y)- Id)(Im(x, Y) -)(O) = V//E(xV)D(d(x, y) - Id)2 Z(xy)D(-m(x Y)- -)2 o D est le domaine commune aux deux images et Id, Im sont les moyennes des images Id et Im respectivement. On considère alors que le recalage courant est donné par la matrice de projection M. On cherche une nouvelle matrice de projection M' qui diffère de M par la translation D = (Dx, Dy, Dz). En considérant un point P qui se projette dans le lieu de coordonnée (u, v) dans l'image initiale obtenue par projection du volume selon la matrice M et dans le lieu de coordonnée (u', v') dans l'image finale obtenue après la translation. Nous avons alors M'P = Ms'v' = MP + DD = + Dy (1) On normalise alors la matrice M par la méthode de Toscani ("Système de calibration et perception du mouvement en vision artificielle" G. Toscani, thèse de doctorat, Université de Paris Sud, Orsay, 1987) et on en déduit: (, s Dz, =s+Dz s + Dz (2) s + Dy V -+Dv z Cette équation permet une interprétation de la translation D. Dx et Dy provoquent une translation de l'image initiale alors que le terme Dz effectue un zoom. Cette homothétie a pour centre le coin supérieur gauche de l'image et s'accompagne donc d'une translation des structures dans

l'image. Cette dépendance des paramètres peut nuire à l'optimisation.

Si, maintenant, nous avons une translation (du, dv) parallèle au plan image associé à un grossissement G par rapport au centre l'image (uc vc), nous obtenons: u'= Gu + Gdu + (1 - G)uc (3 v' = Gv + Gdv + (1- G)vc Ces paramètres sont indépendants et se prêtent bien à l'optimisation de la corrélation grâce au fait qu'ils sont plus intuitifs. Ils ne sont cependant pas équivalents au vecteur D. En effet, s dépend des coordonnées X, Y et Z du point P. Nous avons donc considéré que s était constant sur tout le volume, ce qui revient à considérer que ce volume est ponctuel. On peut estimer que la reconstruction du volume 3DXA occupe une sphère, par exemple de 15 cm de diamètre, située à mi-distance entre le point focal et le plan image, dans un système dont la distance focale est par exemple de l'ordre du mètre. L'approximation est suffisante, car les

matrices obtenues indiquent que la variation de s est de l'ordre de 1%.

L'hypothèse est donc valide et on identifie les deux équations ci-dessus

pour déterminer une bijection entre le vecteur D et le triplet (du, dv, G) .

La procédure d'optimisation de ces trois paramètres est la suivante: recherche exhaustive en basse résolution, par exemple 64x64 pixels. Les bornes de variation des paramètres sont soit contraintes par la machine d'angiographie, soit répondant à un critère large qui est que les

deux images se recouvrent sur au moins un quart de la surface.

- recherche en résolution maximale, qui consiste en une recherche exhaustive uniquement suivant les axes du repère de l'espace des paramètres. On ne modifie qu'un seul paramètre à la fois pour générer l'espace de recherche. Dans le cas de trois paramètres, nous avons donc seulement 6 calculs de corrélation au lieu de 26 à effectuer. Il faut enfin noter que le paramètre G étant très différent de du et dv en amplitude et influence, nous alternons entre son optimisation et l'optimisation commune de du et dv. On distingue ainsi les paramètres qui agissent dans le plan axial (à forte résolution) et ceux qui exercent leur influence dans la

direction orthogonale (à plus faible résolution).

Les paramètres du, dv et G qui maximisent la corrélation entre l'image DSA et l'image MIP, étant trouvés, on modifie en conséquence la matrice de projection M et on projette le volume 3DXA pour générer une nouvelle image MIP. Cette nouvelle image est suffisamment proche de l'image DSA pour pouvoir se prêter à un calcul de flux optique. La différence entre l'image MIP et l'image DSA s'explique par une erreur sur la position du volume 3DXA dans le repère de la caméra. Cette erreur s'apparente à un petit déplacement rigide. Considérons un point P = (X, Y, Z) dans l'espace au temps t et le point P' = (X', Y', Z') atteint au temps t' à partir de P par un petit déplacement rigide, composé de la rotation R = RARBRC (A, B et C étant les angles de rotation autour des vecteurs de base du repère tridimensionnel) et la translation T = (U, V, W). Nous avons alors

P' = RP+T

Dans l'hypothèse o le mouvement est petit, nous avons PP'-P =P=f x P+T avec Q= (A,B,C) (4) Dans le repère interne de la caméra, le point P se projette sur le pixel de coordonnée (u, v) de l'image suivant les formules

X Y

u = a- et v = a- (5) z O cx est le rapport de la distance focale sur la taille pixel pour la caméra. En dérivant cette équation par rapport au temps et en la combinant avec la précédente, nous obtenons: = aB - Cv -uv + -- + AB2 +o 2UW 2+5aa -aWv (6) B A v2 C1!/Vv -aA + Cu ±uav - -+ ce a Z o Cr Z Cette équation correspond presque à une transformation dans le plan image si ce n'est pour le terme en 1/Z. Ceci est le problème posé par la résolution du mouvement tridimensionnel à partir du mouvement apparent de l'équation définissant P' - P. Seuls cinq paramètres sur les six peuvent être résolus, car il subsiste une indétermination sur la profondeur de l'objet dans le champ donnée par la variable Z. Un objet en translation parallèle au plan image génèrera le même mouvement apparent qu'un objet deux fois plus petit, situé deux fois plus près du centre optique, et animé d'un mouvement parallèle au premier, mais deux fois moins fort. Nous pouvons cependant tirer parti de la projection de type MIP. En effet, cette projection associe un unique voxel à chaque pixel. Nous pouvons par conséquent associer, pour chaque pixel, lors du calcul de la projection MIP, la profondeur du voxel élu comme valeur de Z. L'équation déterminant les dérivés de u et de v relie alors complètement les coordonnées des pixels dans l'image aux six paramètres du mouvement recherché.

Il reste à déterminer les valeurs de u et v pour chaque pixel.

D'après l'hypothèse du flux optique, l'intensité de l'objet entre l'image DSA et l'image MIP est constante. Cela n'est pas absolument respecté en ce sens que les deux images sont obtenues par des moyens différents et ne représentent pas exactement les mêmes objets, objet réel dans un cas, reconstruction tomographique dans l'autre. En suivant cependant cette idée, on pose: dI = V.t +I=0 (7) dt En combinant cette équation avec l'équation précédente, nous obtenons pour chaque pixel deux équations que doivent respecter les six paramètres du mouvement. Nous aboutissons alors à un système largement surdéterminé à six inconnues, que l'on résout aux moindres

carrés par détermination du pseudo-inverse.

Les deux phases expliquées ci-dessus reposent sur un grand nombre de projections MIP du volume. Ceci est évident dans la procédure

d'optimisation de la corrélation.

Dans le cas du flux optique modifié, les résultats sont qualitativement bons, mais quantitativement mauvais, l'amplitude du mouvement étant très sous-estimée bien que sa direction soit correcte. On met en oeuvre alors une résolution itérative du petit mouvement résiduel, ce qui demande un certain nombre de projections de type MIP, à savoir une

par itération.

Or, la génération d'une projection de type MIP conique d'un volume 3DXA prend de l'ordre de 1 seconde. A titre d'exemple, dans le cas o l'image à basse résolution est de taille 64*64, le nombre de projections requises lors du calcul exhaustif de la corrélation à basse résolution est le produit du nombre de valeurs des paramètres du, dv et G, respectivement égal à 32, 32 et 35. On doit donc effectuer 35840 calculs de corrélation et par conséquent autant de projections de type MIP. A raison d'une seconde par calcul de projection, la durée devient excessive. Cependant, l'équation définissant u' et v' et celle définissant u et v, définissent chacune une transformation dans le plan image, si nous connaissons les paramètres du mouvement associé (du, dv et G pour la corrélation; A, B, C, U, V et W pour le flux optique). Ainsi, étant donné un mouvement de translation ou de petit déplacement rigide, nous pouvons calculer la nouvelle image MIP à partir de l'ancienne sans passer par une reprojection du volume 3DXA. Le temps de calcul en est

considérablement amélioré.

Dans le cas du flux optique modifié, la transformation dans le plan image ne peut être inversée à cause du terme en 1/Z. Il s'ensuit que les images subissent un artéfact. On peut toutefois localiser les zones dégradées, c'est-à-dire les pixels non atteints par la transformation, et les

éviter dans l'établissement du système d'équations.

L'invention permet donc un recalage entièrement automatique entre l'image DSA et le volume 3DXA. Un essai sur une station UltraSparc cadencée à 200 MHz, a été effectué en 1,5 minute. Le

recalage de référence peut être donné en utilisant un cadre de stéréotaxie.

Etant donné deux images DSA prises sensiblement perpendiculairement et recalées avec un volume 3DXA, on peut pointer dans les images DSA des structures particulières telles que des embranchements d'artères, des courbes prononcées, etc., afin de reconstruire la position du curseur ainsi

définie dans le volume 3DXA.

En d'autres termes, le recalage des images DSA avec des images volumiques des vaisseaux sanguins s'effectue par étalonnage de la machine d'angiographie ou par utilisation de marqueurs, réduction de l'erreur initiale pae estimation d'une translation tridimensionnelle par maximisation de la corrélation entre l'image DSA et l'image en projection conique du volume, par exemple MIP, et utilisation d'une technique de

flux optique modifiée pour affiner le recalage.

Grâce à l'invention, on peut localiser un cathéter de façon précise dans un volume 3DXA, visualiser un flux sanguin dans le volume 3DXA, estimer le mouvement relatif du patient par rapport à la machine d'angiographie durant une acquisition d'image et parvenir à une suppression de l'utilisation du cadre de stéréotaxie lors du traitement d'un patient. On évite également l'utilisation de la présegmentation des images, qui est difficile à effectuer et qui consomme de grandes capacités de calculs et qui nécessite une étape rrianuelle de vérification du résultat par l'utilisateur. De plus, la segmentation perturbe les résultats de recalage à cause d'erreurs, de superpositions d'artères dans l'image DSA, de bifurcations, de tangences des artères, etc., et qui peut donc nuire à la précision. Grâce à l'invention, on atteint une haute précision en utilisant

des variables calculées à partir des données de base de l'image.

Claims (7)

REVENDICATIONS

1. Procédé de recalage automatique d'images d'angiographie bi-

et tridimensionnelles par comparaison de l'image d'angiographie numérique soustraite bidimensionnelle avec des données relatives à une image tridimensionnelle reconstruite à partir de séquences d'angiographie rotationnelle, dans lequel on estime un champ de distorsions dans l'image, on estime une matrice de projection conique et on effectue une approximation d'une transformation rigide dans l'espace égale à la différence entre un recalage initial basé sur le champ de distorsions et sur la matrice de projection conique et un recalage parfait, la transformation rigide comprenant une partie translation et une partie rotation, on effectue une approximation de la partie translation en considérant que la partie rotation est connue, la transformation optimale étant atteinte lorsque la corrélation de l'image d'angiographie numérique soustraite et des données relatives à l'image tridimensionnelle

reconstruite est maximale.

2. Procédé selon la revendication 1, dans lequel on suppose que l'erreur résiduelle de recalage est due à une erreur de positionnement en translation et en rotation de l'image tridimensionnelle reconstruite dans le

repère d'un moyen de prise d'images.

3. Procédé selon la revendication 2, dans lequel la différence entre les dites données relatives à l'image tridimensionnelle reconstruite et l'image angiographique numérique soustraite bidimensionnelle est considérée comme étant due à un petit déplacement suivant une transformation rigide, et on détermine la dite transformation rigide par

une technique de flux optique modifiée.

4. Procédé selon l'une quelconque des revendications

précédentes, dans lequel les dites données relatives à l'image tridimensionnelle reconstruite sont une image bidimensionnelle obtenue par projection conique d'intensité maximale ou par addition des valeurs de

voxels de l'image tridimensionnelle reconstruite.

5. Procédé selon la revendication 4, dans lequel, lors de chaque déplacement de l'image d'angiographie numérique soustraite, la nouvelle image bidimensionnelle est calculée à partir de l'ancienne image bidimensionnelle et ce, sans nouveau calcul de projection de l'image

tridimensionnelle reconstruite.

6. Procédé selon l'une quelconque des revendications

précédentes, dans lequel le recalage initial est donné par un ensemble marqueur de position tridimensionnelle connue.

7. Procédé selon l'une quelconque des revendications 1 à 5, dans

lequel le recalage initial est donné' par étalonnage de la machine d'angiographie.