Brevet FR2662524A1 — publié le 19911129

Procédé de visualisation tridimensionnelle par lancer de rayons avec interpolation locale. L'invention concerne tous les domaines qui requièrent des moyens de visualisation tridimensionnelle (3-D) et dans lesquels les bases de données sont décrites à l'acquisition comme un empilement de coupes bidimensionnelles ou directement en 3-D mais avec un échantillonnage anisotrope dans les trois dimensions. Le procédé intéresse donc, entre autres, toutes les sources d'imagerie médicales fournissant ce type de données : Il peut être implanté sur leur informatique hôte mais aussi sur tout calculateur spécialisé ou d'usage général. Les techniques de visualisation tridimensionnelle actuellement peuvent être identifiées par le type de représentation des objets contenus dans la base de données qu'elles adoptent. Il s'agit - soit de représentations surfaciques (description de la surface de l'objet par des approximations d'ordre 1 ou plus, ou encore par des voxels -volume élémentaire généralement cubique), - soit de représentations volumiques, qui contrairement aux précéden tes conservent l'ensemble des informations présente dans la base de données initiale. Dans chacun de ces grands cadres, de nombreuses méthodes ont été proposées. L'invention concerne la représentation par volume et plus particulièrement les techniques par lancer de rayons (par opposition à des techniques de type projection ou de type codage par arbre octal ou "octree"). Les approches basées sur le lancer de rayons sont apparues dans les années 70 en synthèse d'images : la recherche d'effets spéciaux (réflexion multiple, réfraction, ...) a cependant eu comme conséquence de rendre complexes leurs applications. Leur utilisation pour le rendu de volumes ne remonte qu'à trois ou quatre ans. La visualisation de volumes apporte pourtant une simplification majeure, dés lors que des perceptions réalistes des bases de données seront privilégiées par rapport à la recherche d'effets spéciaux. Seuls alors les rayons primaires seront à calculer. Ces techniques sont maintenant parfaitement éprouvées et dans le domaine médical on pourra se référer aux articles de Drebin (DREBIN R., CARPENTER L. et HANRAHAN P., Volume rendering, ACM Com.Graph., 12, 4, Août 88, pp 65-74), de Levoy (LEVOY M., Display of surfaces from volume data, Thesis, University of North Carolina at Chapel Hill, 1989), de Robb (ROBB R. and BARRILLOT C., Interactive 3-D display and analysis, SPIE vol 939 Hybrid Image and Signal Processing, 1988, pp 173-202) et Trousset (TROUSSET. Y., Visualisation d'objets médicaux tridimensionnels, Thèse, ENST Paris, sept 1987 et brevet français FR 2.613.509). Tous ces travaux ont apporté des améliorations au schéma de base du lancer de rayons, soit sur le plan de la qualité des représentations, soit en terme de temps de calcul. La méthode, décrite par la suite, a l'avantage de réduire les capacités mémoires requises pour la représentation 3-D et ceci sans limiter les autres fonctionnalités (transformations géométriques, découpes par des plans quelconques, détection de sous objets,qualités de l'image, ...). En visualisation d'image 3-D, un compromis doit être trouvé entre trois facteurs principaux : (a) la qualité des images et des fonctions offertes à l'utilisateur, (b) les temps de réponse du système en sachant que l'interactivité homme-machine (même pour des tailles élevées de bases de données) doit être assurée, (c) les coûts des stations d'imagerie. L'invention proposées permet d'apporter un gain significatif sur le troisième point (c) en diminuant notamment la taille de la mémoire nécessaires aux systèmes de traitement. La particularité de la plupart des bases de données 3-D constituées à partir de coupes (ou plans) 2-D est de présenter une anisotropie axiale. C'est-à-dire que la résolution -l'échantillonnagedans le plan est meilleure que dans la direction orthogonale. Tous les travaux, publiés à ce jours, rétablissent avant toute visualisation ou manipulation de la base de données l'isotropie du volume en interpolant entre les coupes. La taille totale de la base de données est ainsi multipliée par un facteur N, avec N#2. Il apparaît avantageux, pour éviter cette augmentation du volume de données, de ne pas pratiquer une interpolation globale préalable et de visualiser directement le contenu de la base anisotrope. Une seule publication (TROUSSET. Y., Visualisation d'objets médicaux tridimensionnels, Thèse, ENST Paris, sept 1987) propose un schéma permettant de faire l'économie d'une telle augmentation mais en contraignant le problème de telle sorte que la condition (a) n'est plus totalement respectée. Sa méthode de visualisation d'une base de données tridimensionnelle formées de coupes empilées est caractérisé par le fait que son incidence de vision est choisie de manière à ce que l'écran est placé perpendiculairement aux coupes et que les rayons suivent une géométrie parallèle. Ainsi, à chaque ligne de l'écran correspond l'image de la tranche d'une coupe.Après un lancer de rayons sur la base anisotrope, il lui suffit d'interpoler des ligne entre les lignes de l'écran. Cette méthode n'est généralisable ni dans le cas où l'incidence de vision est quelconque, ni dans celui d'une géométrie conique. La présente invention vise à fournir un procédé de visualisation par lancer de rayons ne nécessitant pas d'interpolation préalable de la base de données et s'affranchissant des restrictions citées préalablement. C'est-à-dire que - la base de données anisotrope initiale reste inchangée, - seules les parties visibles (surfaces ou volumes internes) à l'observateur et indépendamment du point de vue choisi devront être interpolées lors de la visualisation. Dans ce but, l'invention propose une méthode caractérisée par une séquence d'opérations à réaliser pour chaque rayon et qui se décompose en différente phases (1) définition du rayon (trajectoire) dans la base anisotrope, (2) estimation de la position du point d'intersection du rayon avec le ou les sujets à visualiser dans la base de donnés, (3) interpolation sur le voisinage de ce point estimé (le choix de l'interpolateur est libre, il faut seulement qu'il soit local), (4) recherche sur le rayon dans la zone interpolée du point d'intersection, (5) calcul éventuel de la normale, (6) spécification de la valeur au plan image (ombrage). L'invention sera mieux comprise à la lecture de la description qui suit, d'un mode particulier de mise en oeuvre de l'invention, donné à titre d'exemple non limitatif. Cet exemple reprend les principes de base du lancer de rayons et définit les notations indispensables. Dans un but de simplification, l'exemple est décrit en géométrie parallèle. Son principe s'étend cependant aux projections perspectives. La description se réfère aux schémas qui l'accompagnent dans lesquels - la figure 1 est un schéma qui présente les deux repères (écran et objets) et les relations géométriques entre eux, - la figure 2 est un schéma qui représente le lancer de rayons en géométrie parallèle, - la figure 3 est un graphe destiné à illustrer la suite d'opéra tions du lancer de rayons avec interpolation locale, - la figure 4 est un schéma qui pose le problème de l'intersec tion entre le rayon et l'objet décrit dans la base de données, - la figure 5 est un schéma qui illustre une des méthodes d'esti mation approchée de la position de l'intersection rayon/objet, - la figure 6 est un schéma qui illustre une deuxième méthode d'estimation approchée de cette position, - la figure 7 présente des histogrammes permettant l'évaluation des méthodes proposées. Un rayon dans une base isotrope est définit de la manière suivante La scène 3-D possède son propre repère (1) (voir figure 1) appelé repère objet Rob. Celui-ci dépend de la disposition des coupes tomographiques. Traditionnellement ces dernières sont considérées parallèles au plan (xobXyob) et empilées selon l'axe zob. L'autre système de coordonnées (2) -le repère écran ou repère image tYc est lié à la visualisation. L'écran représente le plan (xec,yec). La scène est placée du côté des Zec positifs. Ces deux repères sont liés par des relations de rotation et de translation. Ces relations dépendent des paramètres de vision. La direction de l'observation est défini e par ses coordonnées sphériques dans le repère objet (1) qui forment une matrice de rotation [R] entre les deux repères. Les origines des repères sont liées par un vecteur de translation. Lors du lancer d'un rayon, la position relative d'un point de ce dernier par rapport à l'objet doit être connue en permanence. Ceci impose constamment de passer d'un repère à l'autre. Toutefois, pour visualiser une base de données formée de voxels, il est plus efficace de tourner l'écran autour de l'objet que de tourner l'objet lui-même. Ce sont donc les transformations écran/objet qui interviendront le plus souvent. Ceci revient à appliquer la formule suivante [Xob,Yob,Zob] = [ R ] [Xec,Yec,Zec]+T Le principe du lancer de rayons est le suivant (voir figure 2). A partir du point de départ (3) P0(x0,y0,z0), le rayon (4) est parcouru de l'écran (5) vers la scène (6). A chaque incrément, une transformation de repère écran/objet est effectuée. Sur le résultat de celle-ci -c'est-à-dire les coordonnées de la position du rayon dans le repère objet- sont vérifiées les conditions d'intersection (par une fonction que nous appellerons CONDITION~D~INTERSECTION). Ces dernières peuvent être fondées sur la densité, le gradient local ou bien d'autres attributs du voxel. Une fois l'intersection constatée, la couleur du pixel (xo,y0) est calculée. Le rayon est une droite dans le repère écran. Les paramètres de ce rayon sont - son point de départ Po (x0,y0,z0), - son vecteur directeur d [ dx, dy, dz ], - une variable t représentant la position le long du rayon. Son équation s'écrit x(t) = x0 + doit y(t) = y0 + z(t) = z0 + dz.t Si les rayons sont parallèles entre eux et à l'axe Zec t vaut [ 0,0,1 ] . L'équation du rayon se ramène à x(t) = xo y(t) = y0 z(t) = z0 + t c'est à dire : [xec(t),yec(t),zec(t)] = [ x0,y0,z0+t ] où x0 et y0 sont les coordonnées du pixel de départ, Les coordonnées du rayon dans le repère objet, pour un t donné, sont [xob(t),yob(t),zob(t)] = [R]. [xec(t),yec(t),zec(t)] + T = [R].[x0,y0,z0+t] + T Cette opération se décompose en trois termes xob(t) = x0.r11 + y0.r21 + (z0+t).r31 + tx yob(t) = x0.r12 + y0.r22 + (z0+t).r32 + ty Zob(t) = x0.r13 + y0 r23 + (z0+t).r33 + tz avec rjj les composante de la matrice [R] et ti les valeurs du vecteur T. La transformation des coordonnées à chaque pas du rayon ralenti le processus. Nous incrémenterons donc directement le rayon dans le repère objet. Cet incrément vaut Xob(ti) - xob(tj-1) = (ti-ti-1).r31 Yob(ti) - yob(ti-1) = (ti-ti-1).r32 zob(t;) - zob(ti-1) = (ti-ti-1).r33 Si nous considérons # #t = ti-ti-1 (pas du rayon), alors nous avons: [xob(ti ) ,yob(ti) X zob(ti) ] = [ xo#(ti.i) ,yob(ti-1 ) , Zob(ti-1) + [ r31,r32,r33 ] . En rappel et en référence avec le schéma synoptique de la figure 3 voici les différentes opérations du lancer de rayons sur une base anisotrope (7) (cette dernière, dans cet exemple, sera composée de coupes présentant une seule anisotropie selon Zob): (8) lancer du rayon dans la base anisotrope, (9) estimation de la position du point d'intersection du rayon avec le ou les sujets à visualiser dans la base de donnés, (10) interpolation sur le voisinage de ce point estimé, (11) détection fine de la surface dans cette zone interpolée, (12) ombrage en fonction des paramètres contenus dans cette zone. Le rayon dans la base anisotrope est définit de la manière suivante: L'échantillonnage du repère Rob n'est plus isotrope. Si D représente la distance inter-coupes et p la taille d'un pixel des coupes originales, le rapport de compression (noté e) de Rob selon l'axe zob vaut : e=D/p. Le ré-échantillonnage isotrope, par interpolation de e-l coupes fictives entre deux coupes originales, définit un nouveau repère appelé Riso. La transformation entre R. et Rob est soumis à une matrice de mise à l'échelle E. 1 0 0 E O 1 0 0 O l/e Le lancer de rayons, décrit précédemment est effectué du repère sur sur Rjso. Les transformations géométriques entre les différents repères sont [xiso,yiso,ziso] = [ R ] . [xob,yob,zob] = [ E ] . [ XisosYisoXZiso ] L'estimation de la position de l'intersection du rayon avec le ou les sujets à visualiser dans la base de données peut se dérouler ainsi: Voir figure 4; si pour les plans originaux (13) les rayons (14) touchent bien l'objet (15), pour les plans "intermédiaires" une surface encore inexistante doit être détectée. l'interpolation sur un voisinage volumique est gourmande en temps de calcul. La réduction de ce voisinage impose une estimation de la position de la surface la plus précise possible. Le module de détection est le suivant : le point d'intersection estimé est un voxel du repère Riso qui appartient à la trajectoire du rayon partant de (xec,yec) et qui répond à un critère appelé ESTIMATION. C'est-à-dire Pa(x,y,Z) : point de détection approchée de la surface viso(xiso,yiso,ziso) : voxel de Rjso de position (xiso,yiso,ziso) Rayon(xec,yec) : rayon lancé du pixel (xe#,yec) ESTIMATION(x,y,z) : module d'estimation de la position de la surface Plusieurs modèles de la fonction ESTIMATION sont étudiés, qui font intervenir des critères géométriques et les propriétés physiques des structures. ler modèle Ce modèle, qui se réfère à la figure 5, est fondée sur la discrétisations originale de la base de données. Un lancer de rayons sur une base anisotrope équivaut, en fait, à un lancer de rayons (16) sur une base isotrope pour laquelle chaque coupe (17) retrouve son épaisseur originale (18). Tout ce passe comme si dans le repère Riso les voxels des plans intermédiaires prenaient la valeur des plans originaux. Nous pouvons donc détecter les voxels (18) de la surface de manière classique c'est-à-dire dès qu'un rayon (16) les pénètrent (19). La condition de détection approchée s'écrit CONDITION~D~INTERSECTION est la fonction citée précédemment qui vérifie si le voxel appartient ou non à l'objet. Le voxel candidat a pour adresse et la partie entière de Cette translation de 0,5 permet de répartir l'épaisseur de part et d'autre de la coupe originale. - 2ème modèle Ce modèle se réfère à la figure 6. Le rayon (20), calculé en nombre réel, est entouré d'autres rayons formés par les coordonnées entières du rayon initiale. Ces rayons, que nous pouvons appeler rayons auxiliaires, touchent obligatoirement les zones connues de l'objet. Les coordonnées de la position approchée de l'intersection du rayon sont déterminées par interpolation des positions d'intersections des rayons auxiliaires avec l'objet. Seule l'interpolation selon Zob est déterminante. Les rayons auxiliaires peuvent être réduits aux nombre de deux : celui qui se trouve au-dessus (z entier supérieur) et celui qui se trouve endessous (z entier inférieur). L'estimation se déroule de la façon suivante : les intersections des rayons auxiliaires supérieurs (21) et inférieurs (22) sont recherchées par CONDITIONP' INTERSECTION(x,b, y,b, z,b+l) pour le rayon auxiliaire supérieur et par CONDITION~D'INTERSECTION(xOb,yOb,zOb) pour le rayon auxiliaire inférieur. Le point de la surface estimé est représenté par l'intersection (23) du rayon (20) et de la ligne (24) définie par les deux points d'intersection trouvés (21) et (22). - 36 modèle Le premier modèle peut être affiné. Le paramètre physique (densité, ...) associé au point détecté est calculé par interpolation linéaire des valeurs des coupes environnantes. Si cette valeur répond aux conditions d'intersection, ce point est retenu, sinon le rayon est incrémenté. avec e : rapport de compression d : distance entre coupes originales selon z Le module de l'interpolation locale est le suivant Toutes les techniques d'interpolation à partir de données locales sont utilisables qu'elles soient du premier ordre (interpolation linéaire) ou plus élaborées (fonctions splines, ... ). Cette liste n'est pas limitative, la seule restriction est le caractère local de l'interpolation. La taille et la forme de la zone environnant le point d'intersection estimé et dans laquelle sera menée l'interpolation peut être définie arbitrairement. Différent critères permettent toutefois un ajustement de ces paramètres (cette liste est non limitative) - la précision du modèle de détection, - le rapport e (si e est faible -2 ou 3- la position estimée est proche de la surface "réelle"), - la technique d'interpolation retenue. Il s'avère que les points estimés par l'opérateur de détection sont proches de la surface de l'objet après interpolation linéaire. Si la technique d'interpola tion locale est linéaire, la région peut se résumer aux proches voisins (connexes) du voxel détecté. Pour les techniques d'ordre supérieur cette région doit être étendue, - la direction de vision.La région doit être étendue selon la direc Lion de islon (celle du rayon), la méthode de calcul de l'ombrage (des voisins qui y sont nécessaires). La détection fine de la surface et l'ombrage se déroulent ainsi Les différents voxels, issus du module d'interpolation, sont stockés dans une petite base auxiliaire. Le rayon, transposé dans cette base auxiliaire et isotrope, procède de manière classique (recherche de la surface, ombrage -en fonction du gradient de la surface, de l'éloignement de l'écran, des sources de lumières, de la nature de la surface, .. ...). La validation de la méthode à été effectuée sur deux plans - l'appréciation subjective des différences entre images par comparaison visuelle, - l'analyse quantitative des valeurs issues des procédures proposées. Cette dernière pouvait être réalisée directement dans l'espace 3-D interpolé après détection. La valeur locale des comparaisons et la complexité de mise en oeuvre ont conduit à préférer une comparaison des images projetées sous différents points de vue. La procédure appliquée est la suivante : 1) acquisition de bases de données 3-D anisotrope; 2) interpolation globale et visualisation sous un angle d'incidence particulier; 3) pour le même point de vue, calcul de l'image par interpolation locale sur le rayon; 4) soustraction point par point des images obtenues en (2) et (3) et statistiques sur les différences. Cette évaluation a été faite pour des nombres de coupes interpolées variables, en utilisant les méthodes de détection décrites précédemment et différents schémas d'interpolation. Les résultats après reparcours local du rayon ne présentant pas de différence significatives, seuls les points de détection approchée ont été utilisés, ce qui à comme intérêt de simplifier la méthode. Pour l'une des bases de données ayant fait l'objet de test, base de résolution 128 X 128 X 136 voxels (volume isotrope), la suppression d'une coupe sur deux a été effectuée. La base de données, rendue ainsi anisotrope, a été visualisée soit après interpolation globale (linéaire), soit par lancer de rayons avec interpolation locale. La figure 7 nous montre les histogrammes des différences entre l'image originale (base de données isotrope de départ) et les images produites après interpolation globale (25) ou interpolation locale (ler modèle de détection (26), 2ème (27) et 3ème (28)). Elle montre que le 3ème modèle de détection donne les meilleurs résultats et que ceux ci sont extrêmement proche de ceux de l'interpolation globale. Naturellement, la présente invention n'est pas limitée à la description qui vient d'être faite. En particulier elle peut s'appliquer à toute les techniques d'accélération du processus de lancer de rayons. L'avantage de ce procédé est qu'il peut être généralisable ou s'insérer à toute séquence d'opérateurs le long du rayon exploitants les distribution locales des paramètres physiques ou simplement des critères géométriques. Revendications 1 - Procédé de visualisation tridimensionnelle du contenu d'une base de données volumique présentant une ou plusieurs anisotropies axiales (échantillonnage différent selon les trois axes), selon une incidence quelconque, utilisant le principe de lancer de rayons et caractérisé en ce qu'il comprend les opérations suivantes - définition du rayon (trajectoire) dans la base anisotrope, - estimation de la position du point d'intersection du rayon avec le ou les sujets à visualiser dans la base de données, - interpolation sur le voisinage de ce point estimé (le choix de l'interpolation est libre, il faut seulement qu'il soit local), - recherche sur le rayon dans la zone interpolée du point d'intersection, - calcul éventuel de la normale en ce point et spécification de la valeur au plan image (ombrage). 2 - Procédé selon la revendication 1, caractérisé en ce que la position du point d'intersection du rayon avec le ou les sujets contenus dans la base anisotrope est estimée lorsque le rayon pénètre un voxel (anisotrope) de cette base, voxel qui répond aux conditions d'intersection souhaitées. 3 - Procédé selon la revendication 1, caractérisé en ce que la position du point d'intersection du rayon avec le ou les sujets contenus dans la base anisotrope est estimée par une méthode qui se sert du principe suivant : le rayon, calculé en nombre réel, est entouré d'autres rayons formés par les coordonnées entières du rayon initiale, que ces rayons, que nous pouvons appeler rayons auxiliaires, touchent obligatoirement les zones connues de l'objet; ainsi les coordonnées de la position approchée de. l'intersection du rayon sont déterminées par interpolation des positions d'intersections des rayons auxiliaires avec l'objet. 4 - Procédé selon la revendication 3, caractérisé en ce que la position du point d'intersection du rayon avec le ou les sujets contenus dans la base anisotrope est estimée par une méthode qui se sert du fait que si la base de données ne présente qu'une seule anisotropie axiale seule l'interpolation selon cette direction est déterminante; ainsi les rayons auxiliaires peuvent être réduits au nombre de deux : ceux qui se trouvent de part et d'autre du rayon selon cette direction; l'estimation se déroule le principe suivant les intersections des deux rayons auxiliaires sont recherchées, le point de la surface estimé est représenté par l'intersection (23) du rayon (20) et de la ligne (24) définie par les deux points d'intersection trouvés (21) et (22). 5 - Procédé selon la revendication 2, selon lequel le paramètre physique associé au point trouvé selon le procédé de la revendication 2 est calculé par interpolation linéaire des valeurs des coupes environnantes; la position du point d'intersection du rayon avec le ou les sujets contenus dans la base anisotrope est estimée si la valeur du paramètre physique répond aux conditions d'intersection. 6 - Procédé selon l'une quelconque des revendications 1 à 5, caractérisé en ce que la position estimée de l'intersection entre le rayon et le ou les sujets contenus dans la base de données est considéré comme position exacte de cette intersection. 7 - Procédé selon l'une quelconque des revendications 1 à 6, caractérisé en ce que la taille et la forme de la zone environnant le point d'intersection estimé et dans laquelle sera menée l'interpolation sont définies soit arbitrairement soit selon un ou plusieurs critères de la liste suivante - la précision du modèle de détection, - le ou les rapports caractérisant l'anisotropie, - la technique d'interpolation retenue, - la direction de vision, - la méthode de calcul de l'ombrage (des voisins qui y sont nécessaires). 8 - Procédé selon l'une quelconque des revendications 1 à 7, caractérisé en ce que les différents voxels, issus du module d'interpolation, sont stockés dans une petite base auxiliaire temporaire et que le rayon, transposé dans cette base auxiliaire et isotrope, procède de manière classique pour trouver le ou les objets à visualiser et pour spécifier la valeur du plan image (ombrage).