FR2802002A1 - Procede de recalage automatique d'images tridimensionnelles - Google Patents

Abstract

Procédé de recalage automatique d'images tridimensionnelles permettant une bonne visualisation des vaisseaux sanguins, par comparaison d'une image numérique tridimensionnelle obtenue au moyen d'un appareil de radiologie et d'une image numérique tridimensionnelle obtenue au moyen d'un appareil à résonance magnétique nucléaire, dans lequel à partir d'un point de correspondance entre les deux images tridimensionnelles, on effectue par traitement portant sur des images bidimensionnelles projetées une estimation d'une rotation capable de recaler les deux images tridimensionnelles, puis on effectue par traitement portant sur des images bidimensionnelles projetées une estimation d'une translation capable de recaler les deux images tridimensionnelles, et on recale l'une des deux images tridimensionnelles par rapport à l'autre.

Description

Procédé de recalage automatique d'images tridimensionnelles. La présente invention concerne le domaine du traitement d'images, en particulier d'images radiologiques 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 reception du faisceau de rayons X est relié à des moyens de traitement d'images permettant la génération d'images tridimensionnelles dites 3DXA à partir d'une série d'images bidimensionnelles prises par le moyen réception, ces images bidimensionnelles étant représentatives de l'ensemble des structures traversées par les rayons X. Dans une image tridimensionnelle 3DXA, les voxels sont isotropes et ont une dimension l'ordre de 300 gm. Dans les applications d'angiographie, les images 3DXA permettent de voir les vaisseaux sanguins qui sont injectés produit de contraste mais on distingue peu les autres tissus.

Les appareils à résonance magnétique nucléaire comprennent des moyens de prise d'images en coupe, une image étant représentative la proportion d'eau présente dans les structures observées. A partir d'une serie de telles images dites MR prises selon différents plans de coupe décalés en translation et/ou en rotation, on sait reconstruire une image tridimensionnelle dite 3DMR. Dans une image tridimensionnelle 3DMR voxels sont anisotropes, c'est-à-dire sont susceptibles d'avoir des dimensions différentes selon les axes d'un repère tridimensionnel. La résolution est de l'ordre d'un millimètre. Dans les applications d'angiographie, les images 3DMR permettent de voir vaisseaux sanguins et les autres tissus.

s'avère important d'obtenir une bonne correspondance entre une image 3DXA et une image 3DMR pour affiner la connaissance des structures observées, notamment des vaisseaux sanguins dans leur environnement.

telle correspondance peut être obtenue par marqueurs externes dont l'utilisation est contraignante et crée des risques d'erreurs. présente invention a pour objet un procédé recalage perfectionné.

présente invention a également pour objet un procédé de recalage précision millimétrique ou submillimétrique à temps de calcul faible.

procédé de recalage automatique, selon l'invention, est destiné des images tridimensionnelles permettant une bonne visualisation des vaisseaux sanguins. On compare une image numérique d'angiographie tridimensionnelle obtenue au moyen appareil de radiologie et une image numérique tridimensionnelle obtenue au moyen d'un appareil à résonance magnétique nucléaire. A partir point de correspondance entre les dites deux images tridimensionnelles, on effectue traitement portant sur des images bidimensionnelles projetées une estimation d'une rotation capable de recaler dites deux images tridimensionnelles, puis on recale l'une des dites deux images tridimensionnelles par rapport à l'autre, puis on effectue traitement portant sur des images bidimensionnelles projetées une estimation d'une translation capable de recaler les dites deux images tridimensionnelles, et on recale l'une des dites deux images tridimensionnelles par rapport à l'autre.

Les images bidimensionnelles sont obtenues par projection des images tridimensionnelles.

Dans un mode de réalisation, on effectue à nouveau par traitement portant sur des images bidimensionnelles projetées une estimation d'une rotation capable de recaler l'image tridimensionnelle recalée par rapport à la dite autre image tridimensionnelle, puis recale l'image tridimensionnelle recalée par rapport à l'autre ou l'inverse, puis on effectue par traitement portant sur des images bidimensionnelles projetées une estimation d'une translation capable de recaler l'image tridimensionnelle recalée par rapport à la dite autre image tridimensionnelle, et on recale l'image tridimensionnelle recalée par rapport à l'autre ou l'inverse. Préférablement, on recalera l'image plus faible résolution par rapport à l'image de plus forte résolution.

préférence, le point de correspondance entre les dites deux images tridimensionnelles est choisi, manuellement ou automatiquement, sur un vaisseau sanguin.

Avantageusement, le traitement d'estimation de la rotation comprend étapes de: - sélection des voxels de chaque image tridimensionnelle, compris entre une surface extérieure et une surface intérieure, deux surfaces englobant le point de correspondance, - projection radiale sur une surface sphérique voxels d'intensité maximale parmi les voxels sélectionnés pour chaque image tridimensionnelle, génération d'une image bidimensionnelle pour chaque image tridimensionnelle, par projection sur un plan dans un repère tridimensionnel centré sur le point de correspondance pour aplatir la sphère, calcul de la corrélation entre les images bidimensionnelles projetées décalage angulaire nul, puis avec un décalage angulaire positif puis négatif selon chaque axe du repère tridimensionnel, - détermination du décalage angulaire autour des trois axes du repère de l'espace tridimensionnel présentant la corrélation maximale entre les images bidimensionnelles.

Dans un mode de réalisation, le traitement d'estimation de la rotation comprend des étapes de réitération des deux étapes de calcul de la corrélation et de calcul du décalage angulaire pour un décalage d'un plus petit nombre de pixels.

Dans un mode de réalisation, les surfaces extérieure et intérieure comprennent des parties sphériques concentriques. Le centre des sphères peut le point de correspondance.

Avantageusement, chaque surface extérieure ou intérieure comprend une partie tronconique, le sommet du cône étant le point correspondance. L'intersection du cône et d'une sphère définit un petit cercle la dite sphère qui limite la dite partie sphérique et le dit tronc cône en définissant la base.

Dans un mode de réalisation, la directrice du cône est un cercle par exemple disposé dans un plan perpendiculaire à une droite passant le centre du dit cercle et le point de correspondance.

Avantageusement, le traitement d'estimation de la translation comprend des étapes de: - sélection des voxels de chaque image tridimensionnelle compris dans un parallélépipède de dimensions données, centré sur point correspondance, - projection selon trois axes d'un même repère tridimensionnel centré sur le point de correspondance, des voxels d'intensité maximale parmi les voxels sélectionnés pour chaque image tridimensionnelle générant trois images bidimensionnelles pour chaque image tridimensionnelle, la projection étant, de préférence, selon des droites parallèles, - calcul de la corrélation entre chaque paire d'images bidimensionnelles projetées selon le même axe à décalage nul, puis avec un décalage positif puis négatif d'un nombre donné de pixels selon chaque axe plan de chaque image bidimensionnelle, - calcul de la moyenne des corrélations pour chaque décalage, - calcul du décalage en translation entre les images tridimensionnelles correspondant au dit décalage présentant la moyenne de corrélations maximale entre les images bidimensionnelles.

parallélépipède peut être un cube par exemple de 16 mm de côté. La taille voxel étant différente en 3DXA et en IRM, si le cube a 64 voxels en 3DXA, il englobera moins de voxels en IRM si l'on désire que les deux cubes soient de la même taille.

Dans un mode de réalisation, le traitement d'estimation de la translation comprend des étapes de réitération des deux étapes de calcul de la corrélation et de calcul du décalage en translation selon un pas plus faible pour un décalage moins important.

En d'autres termes, on effectue un recalage de deux images tridimensionnelles au moyen de traitements portant sur images bidimensionnelles résultant de projections des dites images tridimensionnelles. On évite ainsi d'effectuer des traitements directement sur des images tridimensionnelles, qui seraient lents et coûteux. On peut éviter l'utilisation de marqueurs externes. Après obtention recalages bidimensionnels nécessaires, on sait en déduire les recalages tridimensionnels correspondants.

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 et illustré par les dessins annexés, sur lesquels figure 1 est un schéma d'étapes d'un procédé selon mode de réalisation de l'invention; figure 2 est un schéma d'étapes d'un procédé selon un autre mode de réalisation de l'invention; figure 3 est un schéma détaillé de l'étape 2 figures précédentes; figure 4 est un schéma détaillé de l'étape 3 des figures 1 et 2; la figure 5 est une vue en perspective des sphères utilisées pour la projection d'une première image, selon un mode de réalisation de l'invention; la figure 6 est une vue en perspective des sphères utilisées pour la projection d'une deuxième image, selon un mode de réalisation de l'invention* figure 7 est une vue schématique d'un type de projection utilisé ; et la figure 8 est une représentation plane de la couronne tronquée. Depuis peu, on se sert de reconstructions tridimensionnelles des vaisseaux sanguins dites "3DXA" à partir de séquences d'angiographie rotationnelle 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 Lorraine, Nancy, France, 1996.

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

Ces reconstructions permettent une très bonne appréciation 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, fisosurface, 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 3DXA et les images 3DMR.

La machine de radiologie, une fois étalonnée, fournit un recalage initial qui diffère du recalage recherché ou recalage parfait par transformation rigide (rotation + translation) dans l'espace tridimensionnel.

Comme on peut le voir sur la figure 1, le recalage des images tridimensionnelles commence par une étape 1 de choix d'un point de correspondance entre une image numérique tridimensionnelle composée d'une matrice de voxels et obtenue au moyen d'un appareil à rayons X et d'une image numérique tridimensionnelle également composée d'une matrice de voxels et obtenue au moyen d'un appareil d'imagerie à résonance magnétique nucléaire IRM. A ce jour, le choix du point correspondance est effectué par un opérateur. Toutefois, on pourrait envisager d'automatiser cette tâche.

L'opérateur va choisir un point qui lui semble être vu avec précision sur chacune des deux images. En angiographie, la première image de type 3DXA fait ressortir les vaisseaux sanguins. La deuxième image de type 3DMR fait ressortir à la fois les vaisseaux sanguins et autres tissus avoisinants. L'opérateur choisira donc comme point de correspondance un point d'un vaisseau sanguin qui est donc visible de façon précise sur les deux images à la fois. On peut estimer que la précision du choix du point de correspondance est de l'ordre de 1 à 2 mm. L'opérateur pourra effectuer ce choix par déplacement d'un curseur sur l'une puis l'autre des images tridimensionnelles affichées sur un écran, au moyen d'une souris, d'une boule, d'un clavier ou de tout autre moyen apte à commander un curseur.

Cette étape 1 étant achevée, l'opérateur lance le recalage automatique proprement dit qui commence par une étape 2 d'estimation d'une rotation définie par trois angles de rotation autour de trois axes d'un repère tridimensionnel dont l'origine est le point de correspondance choisi à l'étape 1. la fin de l'étape 2, on connaît donc trois angles notés 6, p et (D, permettant un recalage angulaire entre les deux images tridimensionnelles.

On passe ensuite à une étape 3 dans laquelle on applique ledit recalage défini les trois angles 6, p et (D à l'une des deux images tridimensionnelles. Le recalage peut être effectué aussi bien sur l'image 3DXA que sur l'image 3DMR.

A l'étape 4, on effectue une estimation d'une translation capable de recaler les images tridimensionnelles, précédemment recalées angulairement, l'une par rapport à l'autre. La translation est définie par trois coordonnées X, Y et Z sur chacun des trois axes d'un repère tridimensionnel dont l'origine est le point de correspondance. Avantageusement, le repère est le même que celui utilisé à l'étape 2 pour la rotation.

A l'étape 5, on applique le recalage en translation défini par les coordonnées (X, Y, Z) à l'une des deux images tridimensionnelles. On obtient donc ainsi deux images mutuellement recalées qui présentent donc une correspondance améliorée entre leurs voxels et permettent une meilleure appréciation de l'anatomie du patient et, dans le cas de l'angiographie, de la position des vaisseaux sanguins par rapport aux tissus voisins.

Dans le mode de réalisation illustré sur la figure 2, on a ajouté des étapes de procédé 6 à 9 aux étapes 1 à 5 décrites ci-dessus. En effet, l'estimation du décalage, tant en rotation qu'en translation, est effectuée par comparaison de la corrélation calculée à partir d'images bidimensionnelles et issues des images tridimensionnelles 3DXA et 3DMR, et de la corrélation entre les mêmes images bidimensionnelles mais décalées, soit en rotation à l'étape 2, soit en translation à l'étape 4, d'un petit décalage, par exemple 4 pixels en translation, et par choix du décalage conférant la corrélation maximale, sachant qu'à partir du décalage entre les images bidimensionnelles, on peut calculer un décalage tant en rotation qu'en translation entre les images tridimensionnelles 3DXA et 3DMR.

Il est donc particulièrement intéressant de répéter les étapes 2 à 5 illustrées sur la figure 1, avec un décalage proposé plus faible de façon qu'on accroisse la précision du recalage. Ainsi, les étapes 6 à 9 illustrées sur la figure 2 sont identiques aux étapes 2 à 5, à ceci près qu'on estime la rotation à l'étape 6 et la translation à l'étape 8 avec une précision plus grande, par exemple deux fois plus grande.

Selon la précision souhaitée du recalage, on peut encore réitérer une nouvelle fois ce bloc de quatre étapes avec une précision encore plus grande, jusqu'à obtenir un recalage subpixélique sur les images bidimensionnelles issues des images tridimensionnelles 3DXA et 3DMR.

Sur la figure 3, sont illustrées les sous-étapes effectuées à l'étape 2. L'étape 2 commence par une sous-étape 10 de sélection de certains voxels de l'image 3DXA de l'image 3DMR compris dans un volume dont la définition est identique pour lesdites images 3DXA et 3DMR. Ce volume est délimité par une surface extérieure référencée 19 sur la figure 5 et par une surface intérieure référencée 20. La surface extérieure 19 englobe la surface intérieure 20, les deux surfaces 19 et 20 pouvant posséder des parties communes. Le point de correspondance 21 est compris dans ledit volume et peut également être compris dans les surfaces 19 et 20. On entend ici que le volume comprend les points des surfaces 19 et 20 qui le définissent.

Supposons que nous connaissons un couple de points homologues dans les deux modalités: Pr dans la modalité 3DXA et P f dans la modalité 3DMR. Si ces points sont connus avec une précision extrême, la translation est entièrement déterminée (T= Pr - Pf) et seule la rotation R reste inconnue. Or, toute sphère centree sur le point fixé est invariante par rotation. De même, tout ensemble de points situés entre deux sphères centrées sur le point fixe est invariant par rotation. D'un autre côté, la position des points situés proches du centre de rotation est peu sensible à l'amplitude de la rotation. L'idée donc de ne considérer que les points situés entre une distance minimale Rmin et une distance maximale Rniax du point fixé. La distance minimale définissant la surface intérieure 20 assure que les points considérés apportent une information significative pour caractériser la rotation. La distance maximale définissant la surface extérieure 19 limite l'ensemble de points à l'intérieur de la boîte crânienne. L'ensemble des voxels entre les surfaces 19 et 20 est appelé couronne.

A titre d'exemple, sur la figure 5, on voit que les surfaces 19 et 20 comprennent une partie hémisphérique de rayons différents et une partie plane circulaire fermant la demi-sphère qui est centrée sur le point de correspondance 21.

Dans une variante préférée qui sera expliquée plus loin en référence à la figure 7, les surfaces 19 et 20 comprennent une partie sphérique supérieure à la moitié d'une sphère et une partie tronconique, le sommet du cône étant confondu avec le point de correspondance 21, de même que le centre de la sphère, l'intersection du cône et de la sphère définissant un petit cercle de ladite sphère et la base du tronc de cône. Le cône est commun pour les deux surfaces 19 et 20.

En effet, les sinus, situés en arrière du nez, montrent une prise de contraste lors d'un examen d'IRM injectée au gadolinium. Une grande zone en hypersignal est donc présente dans le volume IRM. Cette zone constitue un pôle vers lequel les artères également en hypersignal peuvent être attirées. Nous éliminons cette zone de manière automatique en retirant un cône de la couronne. La couronne est alors dite tronquée, mais cela ne remet nullement en cause le principe de projection MIP.

La couronne tronquée est extraite dans chaque modalité puis projetée sur la surface externe. La rotation peut alors se retrouver en tournant la surface extraite de l'IRM autour de son centre et en évaluant la superposition avec la surface extraite du 3DXA. Ce processus compare deux surfaces, c'est-à-dire deux ensembles à deux dimensions. Nous pouvons donc passer par une représentation dans le plan pour l'effectuer. A titre d'exemple, le cône peut présenter un angle d'ouverture de 60 . Dans d'autres variantes, on pourrait prévoir que les surfaces 19 et 20 soient définies par une fonction non constante des angles que fait un point donné par rapport aux axes d'un repère tridimensionnel dont l'origine est le point de correspondance 21.

A la sous-étape 11, on effectue une projection d'intensité maximale voxels sélectionnés sur la surface extérieure 19. En d'autres termes, en prenant une demi-droite ayant pour origine le point de correspondance 20 et sécante avec la surface extérieure 19, on retient comme valeur d'intensité du voxel situé à l'intersection de ladite demi- droite et de la surface extérieure 19, l'intensité du voxel, parmi les voxels de la couronne tronquée, de plus forte intensité situé sur ladite demi- droite.

On réduit ainsi de manière automatique le nombre de voxels sans perdre d'information importante sur la rotation. Les artères apparaissent dans les deux modalités en hypersignal. Le reste des voxels peut être assimilé à du bruit : il n'existe aucune correspondance dont on pourrait tirer parti quant à l'intensité des voxels en dehors des artères. Une étape supplémentaire peut alors être franchie vers la simplification des données par projection radiale en maximum d'intensité (MIP) vers la surface extérieure 19.

La distance maximale joue alors un rôle important. Lui donner une valeur limitée permet d'éviter que des zones de peau intersectent la couronne. La peau présentant un hypersignal, cela aurait pour conséquence de créer une tache blanche sur la surface de la sphère externe après la projection MIP, tache qui n'aurait pas son équivalent en 3DXA. Une partie des artères serait noyée dans cette tache et par conséquent une partie de l'information pertinente pour le recalage serait perdue.

De même, choisir une valeur de distance minimale trop faible donnerait une place démesurée aux vaisseaux proches du point 21. Comme cas limite, il suffit d'imaginer ce que donnerait la projection MIP d'une couronne dont le Rmin a été fixé à 1 voxel, si le point central est situé à l'intérieur d'une artère : on trouve une valeur en hypersignal dans toutes les directions. L'image projetée sur la surface externe est donc uniformément blanche donc inexploitable.

A la sous-étape on réalise une projection des voxels calculée à l'étape 11 sur un plan dans un repère tridimensionnel dont l'origine est le point de correspondance (voir aussi figure 8).

La surface est donc décrite par deux angles : 0 qui varie de -?C à+ 7t autour de l'axe [0y); (p varie de 0 (alignement avec l'axe [0y) dans le sens négatif) à (p max tel que n - (P n,ax est l'angle au sommet du cône retranché de la surface.

La représentation plane peut être obtenue en utilisant ces deux variables comme coordonnées polaires dans le plan image<B>:</B> 6 est l'angle polaire et (p est le module. Si nous prenions la couronne entière (cp max = 180 degrés), les sinus apparaîtraient comme une bande blanche sur le pourtour. Cette bande est d'aucune utilité et est potentiellement gênante, se trouvant dans le même domaine d'intensité que les artères. Enfin, la troncature réduit taille de la couronne tout en conservant les voxels utiles, améliorant d'autant le temps de traitement.

Sur un point P = (x,y,z) situé sur la surface externe. Nous pouvons nous ramener trivialement au cas de la sphère unitaire en normalisant P. Les angles cylindriques 0 et (p sont reliés aux coordonnées cartésiennes par les formules

L'angle (p étant dans l'intervalle [0,n], les formules précédentes peuvent être inversées sans difficulté

Cette formule force (p à être positif. Le passage dans le plan image est alors aisé. Notons (u,v)les coordonnées en pixels du point correspondant à P (u = v = 0 dans le coin supérieur gauche de l'image). Considérons que l'image est de taille Nx <I>N</I> pixels, la formule de passage de la sphère vers l'image est alors

Dans le cas où (p nul (x = z = 0 et y = - 1), les cosinus et sinus de 0 sont indéfinis mais leur produit par (p = 0 dans la formule précédente place (cc,v) au centre de l'image.

Considérons maintenant un point de l'image (u,v). Il aura un correspondant sur la sphère 'et seulement s'il se trouve dans le disque qui en forme la trace

3P <SEP> sur <SEP> la <SEP> sphère <SEP> tronquée <SEP> <I>p' <SEP> - <SEP> 2u <SEP> - <SEP> 1 <SEP> .,</I> <SEP> <B>+</B> <SEP> <I>G-L' <SEP> < <SEP> 1</I>
<tb> <I>@r <SEP> @r</I>
<tb> (4) où p est positif.

Si cette condition remplie, nous pouvons extraire les angles décrivant la sphère

On retrouve aisément les coordonnées du point P de la sphère externe correspondant au pixel (u,v). Anoterque sipestnul, alors (p ainsi que son sinus le sont également. Le cosinus et le sinus de 0 ne sont pas définis mais cela ne gêne nullement car la formule (1) donne x = z = 0.

Soit maintenant une rotation, donnée par une matrice 3 x 3 R, appliquée à la couronne tronquée. Pour tout point (cc,v) satisfaisant à la condition (4), on peut trouver un point P correspondant sur la sphère par application des formules inverses (5) puis (1). Ce point P est transformé en P' par la rotation R (P'= R(P-P f)+<B>P f).</B> Ce point P' appartient également à la sphère, celle-ci étant globalement invariante par la rotation. L'application des formules directes (2) puis (3) donne la localisation du pixel (ii ;v') correspondant au point P'. La relation entre (u,v) et (u',v') est ainsi déterminée. Nous avons vu que toutes les formules utilisées étaient bijectives. L'effet d'une rotation est donc une bijection du plan de l'image dans lui-même. Nous pouvons ainsi travailler dans un espace à deux dimensions au lieu de trois. Les effets de bord dûs aux points entrant (P') ou sortant (P) de la zone conique peuvent se détecter grâce à la condition (4) et se traitent de manière triviale.

On peut aussi effectuer la projection en considérant que les différents voxels sont projetés, pour un plan donné, parallèlement les uns aux autres et parallèlement à une droite perpendiculaire au plan. Toutefois, cette technique de projection limite le choix des surfaces 18 et 19 à une derni-sphère ou plus généralement à une surface délimitée par un plan comprenant le point de correspondance 20 et parallèle au plan sur lequel on projette.

On pourra utiliser une technique de projection différente qui correspond en quelque sorte à dérouler la surface 18 sur le plan sur lequel elle est projetée. Un exemple simple d'une telle méthode est qu'un voxel de coordonnées sphériques (p, 0, cp) aura pour coordonnées cartésiennes (0,cp) dans le plan sur lequel il est projeté. En d'autres termes, une coordonnée cartésienne d'un pixel projeté est une fonction linéaire des coordonnées sphériques du voxel d'origine, et non une fonction sinusoïdale. On obtient ainsi à la fin de la sous-étape 12 une image bidimensionnelle obtenue par projection dans le repère tridimensionnel et ce pour chacune des images tridimensionnelles 3DXA et 3DMR d'origine.

Cette représentation semble bien adaptée au problème représenter 0 et (p non plus comme coordonnées polaires mais comme coordonnées cartésiennes (l'axe des 0 croissants étant suivant les lignes de l'image et l'axe des (p croissants suivant les colonnes). Les formules de la transformation dans le plan image sont simples (des translations, sauf pour la rotation autour de l'axe [0y) qui est un peu plus complexe). Cependant, une certaine continuité manque à ce type de représentation les pixels situés sur les bords verticaux des images passeraient d'un côté à l'autre sous l'effet d'une rotation suivant 0. Mais l'inconvénient le plus sérieux est que cette représentation déforme les structures vasculaires, suivant leur orientation. Cette représentation s'avère une solution pour accélérer le temps calcul de l'algorithme mais la représentation que nous pouvons par comparaison qualifiée de "polaire"est préférée à la représentation que nous nommerons alors<I>"cartésienne".</I> parce qu'elle déforme peu les structures et produit des images qui sont plus aisément interprétables : tout passe comme si l'observateur était situé au centre des sphères et possédait une vision sur 2 cp max degrés.

A la sous-étape 13, on calcule la corrélation entre les deux images bidimensionnelles projetées à décalage angulaire nul, autrement dit, on calcule la corrélation entre l'image bidimensionnelle provenant de l'image 3DXA et l'image bidimensionnelle provenant de l'image 3DMR. On peut utiliser à cet effet une corrélation croisée obtenue en multipliant les valeurs d'intensité des pixels des deux images bidimensionnelles.

Ensuite, on calcule de la même façon une valeur de corrélation pour les mêmes images bidimensionnelles mais décalées l'une par rapport à l'autre d'un certain angle selon l'un des axes du repère tridimensionnel. On réitère ce même calcul pour un décalage négatif du même angle selon le même axe. Puis on effectue les mêmes calculs pour des décalages selon les deux autres axes. On compare les sept valeurs de corrélation obtenues pour une image bidimensionnelle et on retient le décalage conférant la corrélation la plus forte entre les deux images bidimensionnelles, l'une issue de l'image 3DXA et l'autre issue de l'image 3DMR.

On peut remarquer trois problèmes potentiellement gênants la comparaison des images - le fond de l'image 3DMR n'est pas uniforme, avec en particulier de larges zones en hypersignal dues aux grosses artères et aux sinus (partie inférieure de l'image) et des taches claires ici et là dans l'image ; - le diamètre des artères n'est pas exactement le même : des problèmes de résolution dans le volume IRM apparaissent - toutes les artères visibles en 3DXA ne le sont pas dans l'image 3DMR et surtout des vaisseaux (notamment des veines) viennent s'ajouter artères intéressantes dans l'image 3DMR. Pour rehausser les vaisseaux, on utilise un simple opérateur morphologique dit "chapeau haut de forme"<I>(Top-Hat).</I> Pour rappel, il s'agit simplement d'une ouverture, suivie d'une soustraction de l'image ouverte à l'image originale. L'élément structurant est un disque centré sur le pixel d'intérêt et dont le rayon va dépendre de la taille des artères nous souhaitons garder. De manière formelle, si nous désirons conserver toutes les artères dont le diamètre n'excède pas d millimètres, le diamètre D de l'élément structurant sera déterminé par

où d/Rmin est l'angle maximal sous lequel vue une longueur équivalente à la largeur des artères à conserver et (p max /N est la taille angulaire d'un pixel de l'image de représentation plane. Cet opérateur est appliqué aux deux images (3DXA et 3DMR).

En plus du deuxième problème, l'opérateur morphologique fait apparaître plus fortement les variations d'intensité présentes dans l'IRM instabilité de l'intensité le long d'une artère et intensité plus faible pour les artères de petit calibre. Nous employons un critère d'intercorrélation normalisée et centrée, cette fois, par rapport à la moyenne locale des images (en pratique un voisinage carré de 15 pixels de côté est utilisé).

Enfin le dernier problème est résolu en initialisant les paramètres de rotation par une phase de recherche exhaustive du maximum du critère à basse resolution, ce qui permet d'éviter les maxima locaux.

D'autres critères de similarité classiques sont évidemment envisageables à la place de la corrélation : information mutuelle ou rapport de corrélation par exemple.

A la sous-etape 14,à partir des différentes valeurs de décalage retenues à la sous-étape 13 comme fournissant la corrélation la plus forte entre les images bidimensionnelles, on calcule le décalage angulaire entre les images tridimensionnelles 3DXA et 3DMR correspondant au décalage obtenu à la fin de sous-étape 13.

Supposons à présent que nous connaissons parfaitement la rotation mais le couple de points homologues (Pr,Pf) est approximativement connu. Nous avons toujours une estimation de la translation (T = - Pf) mais celle-ci est désormais approximative et demande à être raffinée.

Il n'est besoin de connaître le volume entier pour déterminer une translation. Un sous-volume correctement choisi suffit, l'effet de la translation étant la même en tout point du volume (contrairement à la rotation, l'amplitude de la transformation ne dépend pas de la zone observée). Imposons donc à ces points appariés d'être au voisinage d'une structure artérielle ne présentant aucune invariance par translation, autrement dit la structure doit s'épanouir dans les trois dimensions : un exemple typique étant une bifurcation au niveau de laquelle les trois artères prennent directions très différentes. La translation peut alors être estimée en mettant en correspondance les voisinages de ces points.

Sur la figure 4, sont illustrées les sous-étapes 15 à 18 de l'étape 4. A la sous etape 15, on sélectionne les voxels de chaque image tridimensionnelle 3DXA, 3DMR compris dans un parallélépipède, par exemple un cube dimensions données, centré sur le point de correspondance 21. On peut choisir un cube de 16 mm de côté, ce qui correspond environ à la taille de 64 voxels en 3DXA. A la sous-étape 16, on effectue trois projections, chacune selon un axe d'un même repère tridimensionnel. Chaque face du cube pourra avantageusement être perpendiculaire un axe du repère tridimensionnel. On effectue chaque projection en conservant, pour une lignée de voxels donnée, la valeur d'intensité du voxel présentant l'intensité maximale. Ainsi, si le voxel de coordonnées cartésiennes x1, y1, z1 présente l'intensité la plus forte de tout le parallélépipède sélectionné à la sous- étape 15, on retrouvera cette même valeur d'intensité dans les images bidimensionnelles projetées aux coordonnées suivantes : x l, y 1 pour l'image perpendiculaire à l'axe Z et dont les voxels ont été projetés selon ledit axe Z, (x1, z1) dans l'image bidimensionnelle obtenue par projection selon l'axe Y, et (y1, z1) dans l'image bidimensionnelle obtenue par projection selon l'axe X.

A la fin de la sous-étape 16, pour chaque image tridimensionnelle 3DXA et 3DMR, on obtient trois images bidimensionnelles projetées 2DXAP 2DMRP.

A la sous-étape 17, on calcule la corrélation entre les deux images 2DXAP et 2DMRP projetées selon l'axe X à décalage nul. Puis on effectue ce même calcul pour les mêmes images décalées d'un nombre donné de pixels, par exemple +4, puis le même calcul pour les mêmes images décalées négativement du même nombre donné de pixels, par exemple -4. On effectue les mêmes calculs avec les images 2DXAP et 2DMRP projetées selon l'axe Y, puis celles projetées selon l'axe Z.

Les images de chaque couple d'images homologues (axiales, sagittales et coronales) sont comparées par le critère de corrélation centré normalisé. La taille du sous-cube étant réduite, ces images présentent toutes un fond à peu près uniforme, même en IRM. Le centrage est donc fait par rapport à la moyenne des pixels situés dans la zone commune aux deux images. La normalisation permet de renforcer la robustesse face à une faible taille de zone commune.

On utilise un critère de corrélation faible issu de la moyenne des trois scores : un score fort dans au moins deux images conduit à un critère fort.

Etant donné la faible résolution de l'IRM par rapport au 3DXA, et donc le fort niveau d'interpolation présent dans les images MIP du sous- cube extrait de l'IRM, ainsi que la possibilité d'artéfacts (par exemple dûs au flux) dans l'IRM, ce critère faible permet de ne pas être pénalisés par ces imprécisions.

A la sous-étape 18, on calcule pour chaque décalage, à savoir décalage nul, positif sur l'axe des X, négatif sur l'axe des X, positif sur l'axe des Y, négatif sur l'axe des Y, positif sur l'axe des Z et négatif sur l'axe Z, la moyenne des corrélations obtenue à la sous-étape 17. On retient ensuite le décalage présentant la moyenne de corrélation maximale entre images 2DXAP et 2DMRP.

On passe ensuite à l'étape 5 de recalage au moyen du décalage calcule en coordonnées cartésiennes à la sous-étape 18. Le fait de limiter les voxels à ceux contenus dans un parallélépipède de dimensions données, permet de réduire fortement les temps de calcul nécessaires.

De préférence, on utilise un seul et même repère tridimensionnel pour toutes les étapes et sous-étapes.

Sur la figure 5, est représenté en perspective un exemple d'image tridimensionnelle de type 3DXA avec un réseau d'artères 25 qui sont bien visibles grâce à l'injection de produit de contraste. De cette image, on retient voxels compris entre les surfaces 19 et 20, la partie sphérique de la surface extérieure 19 présentant un rayon de 4 cm et la partie sphérique de la surface intérieure 20 présentant un rayon de 2 cm. Bien entendu, ces valeurs de rayon sont indicatives. Elles sont adaptées à une image tridimensionnelle du cerveau humain, le grand rayon de 4 cm évitant d'avoir dans l'image les structures extérieures au cerveau notamment la peau, et le petit rayon de 2 cm évitant de prendre des vaisseaux trop près du centre des sphères et qui prendraient trop d'importance lors de la projection. L'algorithme d'optimisation choisi suppose tour à tour la rotation, puis la translation parfaitement connue. Lorsque nous estimons la rotation, on suppose donc que l'erreur en translation petite devant l'erreur en rotation. On peut estimer que ce problème d'équilibre entre les deux erreurs se rencontre surtout à l'initialisation : les precisions finales de l'optimisation pseudo-exhaustive étant suffisamment petites pour négliger ce problème par la suite. On suppose donc que pire erreur en translation que l'on puisse avoir est rencontrée à finitialisation manuelle des points Pr et Pf. Nous l'avons estimée équivalente à résolution de l'IRM, soit 1 mm. L'angle sous lequel une telle erreur est vue depuis la distance Rmin est llRmin radians. L'initialisation de la rotation est faite avec N = 64. La taille angulaire d'un pixel de l'image de la couronne tronquée est donc : 2(p ma.,/N. Il faut donc comparer Rmin à N12 #p rnax 12,22 mm, d'où Rmin = 2 cm.

Les deux sphères sont centrées sur le point de correspondance 21 déterminé avant le lancement du procédé de recalage.

Sur la figure 6, est représenté une image tridimensionnelle de type 3DMR, sur laquelle on voit également le point de correspondance 21 et les surfaces extérieure 19 et intérieure 20. Ainsi, on conserve pour le traitement de recalage les mêmes voxels dans l'image 3DXA et dans l'image 3DMR aux effets de bord près.

Sur la figure 7, est illustrée de façon schématique la projection qui peut être utilisée lors de l'étape d'estimation de la rotation. Les voxels compris entre les surfaces extérieure 19 et intérieure 20 sont projetés tout d'abord sur la surface extérieure 19. A titre d'exemple, on a représenté la demi-droite 22 qui présente un angle (Dl par rapport à l'axe de projection 23. Tous les voxels de coordonnée (Dl se projettent sur la surface extérieure 19 par projection d'intensité maximum. On retient donc la valeur d'intensité la plus forte des voxels de coordonnées < Dl.

Le voxel projeté de coordonnée (Dl se situe donc à l'intersection de la demi-droite 22 et de la surface extérieure 19 et est projeté en rotation sur le plan 24 sur lequel se forme une image bidimensionnelle. La projection permet de dérouler la partie sphérique de la surface extérieure 19. En d'autres termes, la coordonnée cartésienne x1 du voxel projeté correspondant au voxel de coordonnée sphérique < Dl sera telle que x1 est proportionnel à < Dl et ce contrairement à la projection classique où x1 sera proportionnel à une fonction sinusoïdale de (Dl.

Pour la translation, le choix d'un premier décalage proposé de quatre pixels puis d'une réitération des étapes d'estimation de rotation et de translation avec un décalage proposé de deux pixels puis d'un pixel, résulte de la précision du choix du point de correspondance dans les images 3DXA et 3DMR qui est estimé entre 1 et 2 mm et de la résolution de l'image 3DXA qui est de l'ordre de 0,3 mm. Ainsi, le pixel le plus fin d'une image bidimensionnelle projetée a une taille de l'ordre de 0,3 mm. En commençant le recalage avec quatre pixels, soit un décalage de 1,2 mm, puis deux pixels, soit 0,6 mm, puis un pixel, soit 0,3 mm, on parvient à rattraper au maximum une erreur en translation de 2,1 mm, ce qui est largement suffisant par rapport à l'erreur raisonnable que peut commettre un opérateur lors du choix du point de correspondance.

Sur la figure 8, est illustré un type de projection préférée dite "polaire". L'angle (p en coordonnées sphériques qui est compris entre - cPmax + cpmax détermine le module en coordonnées polaires.

Grâce à l'invention, on dispose d'un moyen de recalage de deux images tridimensionnelles qui est rapide et peu coûteux dans la mesure où seules images bidimensionnelles sont traitées, en particulier lors calculs de corrélation. On peut ainsi utiliser plus facilement des images tridimensionnelles de radiologie et de résonance magnétique lors d'interventions chirurgicales, ce qui augmente la sécurité et l'efficacité telles opérations.

Claims (6)

REVENDICATIONS

1. Procédé de recalage automatique d'images tridimensionnelles permettant une bonne visualisation des vaisseaux sanguins, par comparaison d'une image numérique d'angiographie tridimensionnelle obtenue au moyen d'un appareil de radiologie et d'une image numérique tridimensionnelle obtenue au moyen d'un appareil résonance magnétique nucléaire, dans lequel: - à partir d'un point de correspondance entre dites deux images tridimensionnelles, on effectue par traitement portant sur des images bidimensionnelles projetées une estimation d'une rotation capable recaler les dites deux images tridimensionnelles, puis - on recale l'une des dites deux images tridimensionnelles par rapport à l'autre, puis - on effectue par traitement portant sur images bidimensionnelles projetées une estimation d'une translation capable de recaler les dites deux images tridimensionnelles, et - on recale l'une des dites deux images tridimensionnelles par rapport à l'autre.

2. Procédé selon la revendication 1, dans lequel on effectue à nouveau par traitement portant sur des images bidimensionnelles projetées une estimation d'une rotation capable de recaler l'image tridimensionnelle recalée par rapport à la dite autre image tridimensionnelle, puis on recale l'image tridimensionnelle recalée par rapport à l'autre ou l'inverse, puis on effectue par traitement portant sur images bidimensionnelles projetées une estimation d'une translation capable de recaler l'image tridimensionnelle recalée par rapport à la dite autre image tridimensionnelle, et on recale l'image tridimensionnelle recalée par rapport à l'autre ou l'inverse.

3. Procédé selon la revendication 1 ou 2, dans lequel le point de correspondance entre les dites deux images tridimensionnelles est choisi, manuellement ou automatiquement, sur un vaisseau sanguin.

4. Procédé selon l'une quelconque des revendications precédentes, dans lequel, chaque image tridimensionnelle étant composée d'une pluralité de voxels, le traitement d'estimation de la rotation comprend des étapes de: sélection des voxels de chaque image tridimensionnelle, compris entre une surface extérieure et une surface intérieure, les deux surfaces englobant le point de correspondance, projection radiale sur une surface sphérique voxels d'intensité maximale parmi les voxels sélectionnés pour chaque image tridimensionnelle, génération d'une image bidimensionnelle pour chaque image tridimensionnelle, par projection sur un plan dans repère tridimensionnel centré sur le point de correspondance pour aplatir la sphère, - calcul de la corrélation entre les images bidimensionnelles projetées à décalage angulaire nul, puis avec un décalage angulaire positif puis négatif selon chaque axe du repère tridimensionnel, - calcul du décalage angulaire entre les images tridimensionnelles correspondant au dit décalage présentant la corrélation maximale entre les images bidimensionnelles.

5. Procédé selon la revendication 4, dans lequel le traitement d'estimation de la rotation comprend des étapes de réitération deux étapes de calcul de la corrélation et de calcul du décalage angulaire pour un décalage moins important.

6. Procédé selon la revendication 4 ou 5, dans lequel surfaces comprennent des parties sphériques concentriques. Procédé selon la revendication 6, dans lequel chaque surface comprend une partie tronconique, le sommet du cône étant point de correspondance. Procédé selon l'une quelconque des revendications précédentes, dans lequel le traitement d'estimation de la translation comprend des étapes de: sélection des voxels de chaque image tridimensionnelle, compris dans un parallélépipède de dimensions données, centré sur le point de correspondance, projection selon trois axes d'un même repère tridimensionnel centré sur le point de correspondance, des voxels d'intensité maximale parmi les voxels sélectionnés pour chaque image tridimensionnelle, genérant trois images bidimensionnelles pour chaque image tridimensionnelle, - calcul de la corrélation entre chaque paire d'images bidimensionnelles projetées selon le même axe à décalage nul, puis avec décalage positif puis négatif d'un nombre donné de pixels selon chaque axe du plan de chaque image bidimensionnelle, - calcul de la moyenne des corrélations pour chaque décalage, - calcul du décalage en translation entre les images tridimensionnelles correspondant au dit décalage présentant la moyenne de corrélations maximale entre les images bidimensionnelles. 9. Procédé selon la revendication 8, dans lequel le traitement d'estimation de la translation comprend des étapes de reitération des deux étapes de calcul de la corrélation et de calcul du décalage en translation selon un pas plus faible pour un décalage moins important. 10. Procédé selon l'une quelconque des revendications 8 ou 9, dans lequel la projection est effectuée, pour chaque , selon des droites parallèles audit axe.