FR 2623642 A1 19890526 FR 8716156 A 19871123 Procédé de segmentation d'images angiographiques vasculaires par suivi vectoriel d'axes et détection de contours asservie L'invention se rapporte au domaine du traitement d'images médicales, et en particulier aux images numériques angiographiques obtenues par les systèmes de radiologie numérique. La segmentation des images issues d'un système de radiologie, notamment numérique, consiste a extraire de l'image les éléments d'informations utiles. Pour les images angiographiques des réseaux vasculaires, cette segmentation doit permettre de retrouver dans l'image les artères et veines du réseau. Dans ce domaine, plusieurs méthodes ont déjà été développées pour réaliser la segmentation d'images angiographiques. Une première méthode - décrite par NGUYEN T.V. -SKLANSKY J. "A fast skeieton-Finder for coronary Arteries", IEEE Comp. Society, n0 742, vol.1, pages 481-483, October 1986 - utilise, pour rechercher les pixels caractéristiques, deux fenêtres de dimensions limitées, l'une de côtés parallèles aux axes de l'image, l'autre de côtés orientés à 450 de ceux de la première. Ces deux fenêtres sont déplacées sur toute l'image, et les maxima, c'est-à-dire les niveaux de gris les plus élevés, sont recherchés dans ces fenêtres suivant quatre directions (deux directions par fenêtre).Si le nombre de maxima suivant une direction donnée est suffisant, ces maxima sont validés comme faisant partie de la "crête" d'un vaisseau. La crête correspond à ce que l'on peut également appeler l'axe du vaisseau dans l'image. Les contours sont obtenus par application de masques, et la comparaison entre les éléments de contours et les éléments d'axe permet d'éliminer un certain nombre de fausses détections. Une autre méthode - décrite par HOFFMANN K. R et ai dans un article intitulé "Automated tracking of the vascular tree in DSA images using a sqiiare-box region of search, algorlthm", SPIE vol.626, Medicine XIV pages 326-33, 1986 consiste à détecter les points d'un vaisseau par balayage du périmètre d'un carré centré sur un point déjà connu d'un vaisseau. Pour l'initialisation, le premier point est indiqué par l'utilisateur. La dimension de la fenêtre est variable suivant le niveau de gris du point courant. Les bifurcations sont recherchées sur le périmètre du carré, à l'extérieur de la structure en cours d'étude.Les vaisseaux détectés sont marqués par des points espacés de manière variable, et la largeur des vaisseaux est obtenue par comparaison du profil réel de l'image du vaisseau avec le profil de la projection d'un élément cylindrique modélisant le vaisseau D'une manière générale, ces méthodes conduisent à l'extraction de points caractéristiques, considérés comme des points d'axe ou des points de contour, qui n'ont pas obligatoirement de liens les uns avec les autres, si ce n'est lorsque l'on cherche une suite de maxima suivant une direction donnée dans une fenêtre prédéterminée. L'inconvénient de ces systèmes est qu'ils produisent donc des structures discontinues, et qu'il est nécessaire ensuite de les traiter, au moyen de calculs complexes, pour les rendre continus.De plus, les éléments détecteurs de contours sont généralement peu fiables et il est nécessaire, compte tenu du nombre de points, de chercher des cohérences pour aboutir à une définition à peu près correcte des coordonnées des points des vaisseaux, pour reconstituer le réseau vasculaire. L'invention a pour objet un procédé de segmentation d'images angiographiques, qui n'a pas les inconvénients sus-mentionnés et qui,# de par sa conception, conduit nécessairement à une .structure continue et fermée. Cette structure ne comporte pas de trous, même aux bifurcations. Selon l'invention, les principales caractéristiques du procédé de segmentation d'images angiographiques vasculaires sont les suivantes - les axes des vaisseaux sont déterminés par un suivi vectoriel produisant une suite de points connexes - le suivi de contour est directement lié aux axes déterminés dans la phase précédente, en partant de l'hypothèse que, de part et d'autre de l'axe, il existe un contour ; le suivi de contour à partir de cette ligne centrale conduit également à des suites de points connexes - le lien étroit entre les axes et les contours, chaque contour étant affecté à un axe parfaitement identifié, permet dans une phase finale la rectification de la position des axes. Ce procédé permet d'obtenir une très bonne précision dans les données extraites, et permet donc de constituer une base de données fiable, notamment pour la reconstruction tridimensionnelle des vaisseaux à partir d'un petit nombre d'images de projection. De plus, le procédé selon l'invention extrait de manière automatique et précise les vaisseaux à partir d'une radiographie numérique, sans prétraitement de l'image. Bien entendu, le procédé de segmentation selon l'invention met en oeuvre; pour certaines de ses phases; des algorithmes déjà décrits En particulier, le procédé selon l'invention utilise également une fenêtre pour la détection de points de l'axe, mais seulement pour un ensemble de points de départ. Par ailleurs, un algorithme de recherche d'un contour par programmation dynamique a déjà été décrit - par exemple par POPE. D . L et ai, dans un article intitulé "Left ventricular Birder Recognition Using a Dynamic search algorithm, Radiology nb 155, pages 513-518, 1985 -. Le procédé selon l'invention utilise également une programmation dynamique pour la recherche des contours mais qui s'appuie sur les lignes centrales. Le procédé selon l'invention est optimisé, notamment - par affinement du critère du choix de la direction de l'exploration d'un vaisseau, du calcul, et des tests - par l'intégration d'une procédure automatique de détection des points germes sur toute structure connexe à celle étudiée - par la contrainte de maintien du suivi de contour par rapport à la ligne centrale - par la connexité obtenue tant pour les axes des vaisseaux que pour les contours - par l'affectation de tout contour à une branche parfaitement identifiée - et la possibilité de rectification de la position de l'axe qui en découle. Selon l'invention, un procédé de segmentation d'images angiographiques numérisées pour l'extraction des axes et des contours de réseaux vasculaires, est caractérisé en ce qu'il comporte - une première phase au cours de laquelle un squelette du réseau est extrait, par une méthode de suivi vectoriel des lignes centrales des vaisseaux, ou axes, recherchant, à partir de points isolés reconnus comme appartenant au réseau par vérification d'un critère Ex et affectés chacun d'une direction, le ou les points connexes suivants, dans un nombre limité de directions prédéfinies fonction de la direction affectée au point courant, et qui satisfont un critère prédéfini Cd fonction des niveaux de gris d'un ensemble de points prédéfinis, - et une phase de détection asservie des contours s appuyant sur le squelette du réseau obtenu dans la phase précédente pour rechercher, de part et d'autre des axes, un contour également formé de suites de points connexes, correspondant à des points de transition importante du niveau de gris recherchés dans des zones prédéfinies de l'image. L'invention sera mieux comprise et d'autres caractéristiques apparaîtront à l'aide de la description qui suit en référence aux figures annexées. - Les figures 1 et 2 schématisent les différentes étapes du procédé de segmentation d'images selon l'invention; - la figure 3 est un schéma explicatif montrant les directions d'investigation - la figure 4 est un autre schéma explicatif illustrant les données utiles au suivi du vaisseau - la figure 5 illustre le critère de choix de la direction - la figure 6 illustre les vecteurs retenus pour l'établissement d'une matrice de gradients - les figures 7a et 7b illustrent les formes des ^ zones de recherche de contours, autour d'un segment d'axe, selon l'inclinaison du segment - la figure 8 illustre, dans le cas le plus simple le chaînage des éléments de contours - la figure 9 illustre également le chaînage des éléments de contours dans un cas où les zones de recherche ne sont pas jointives - la figure 10 illustre le ~ traitement du premier segment d'une branche. La description qui suit du procédé de segmentation d'images angiographiques selon l'invention donnée à titre d'exemple non limitatif fait référence aux figures 1 et 2 respectivement pour les étapes de traitement relatives à ltextraction du "squelette" du réseau vasculaire, et pour les étapes de traitement relatives au suivi des contours. Une image obtenue par angiographie numérique, résultant du balayage d'une zone analysée par un faisceau de rayons X, est une suite de points de niveau de gris variable, par exemple entre O et 256 si les niveaux de gris sont donnés sur 8 bits. Ces rayons sont en effet absorbés différemment suivant les tissus traversés. Cette image est organisée en lignes, d'indice variable i, et en colonnes, d'indice variable j. D'une manière générale, la segmentation d'images angiographiques vasculaires par suivi vectoriel des axes selon l'invention consiste à extraire de cette image angiographique numérique un "squelette" qui par hypothèse est continu, et qui correspond aux axes des vaisseaux du réseau vasculaire. Puis, en partant de l'hypothèse que, de part et d'autre de chaque axe. il existe un contour du vaisseau correspondant, on recherche les contours associés par une détection asservie. Pour la création du squelette, une phase d'initialisation permet de rechercher des points dans l'image ayant une bonne probabilité d'appartenir au squelette. Du fait de la méthode d'analyse, les points appartenant au squelette sont des points de niveau de gris "maximum", c'est-à-dire dont les points voisins, lors d'un balayage d'une ligne ou d'une colonne, ont des niveaux de gris qui lui sont inférieurs. Dans la description qui suit le niveau de gris associé à un point est soit le niveau absolu de gris du pixel correspondant, soit la moyenne des niveaux de gris des pixels d'un voisinage, par exemple 3x3, centré sur ce pixel; cette deuxième mesure est préférable pour s'affranchir de défauts dus au bruit. Pour l'initialisation, la première phase consiste donc à rechercher des maxima sur une ligne (ou sur une colonne de l'image). Un premier balayage d'une ligne i (ou d'une colonne j) permet de détecter les points de niveau de gris maximum, qui ont donc une certaine probabilité d'appartenir à un vaisseau. La phase suivante est alors l'éventuelle validation de ces maxima détectés après une séquence de test : à partir de la moyenne des niveaux de gris dans une zone peu vascularisée de l'image, le "niveau local" de gris, est calculé pour les différents maxima détectés dans la phase préalable par différence avec cette moyenne, et en tenant compte éventuellement des niveaux des points voisins. La moyenne MOY préalablement calculée peut être par exemple 118 pour une échelle de gris à 256 niveaux. Le niveau local du point détecté peut être défini de la manière suivante M(i,j+k) étant le pixel situé sur la ligne i et sur la colonne j+k et N(i,j+k) étant son niveau absolu de gris ; on tient ainsi compte des points voisins de la même ligne dans les colonnes précédente et suivante. Les points de niveau de gris maximum sont retenus seulement si - NIV1 est croissant sur un nombre prédéterminé de points consécutifs, sur une ligne, - NIV1 devient décroissant ensuite, - il n'y a pas de maximum détecté préalablement dans un voisinage VOIS prédéterminé de ce point. Un maximum retenu n'appartient au squelette du réseau que s'il satisfait alors un critère d'appartenance résultant de l'étude d'un voisinage de ce maximum, ou plutôt de l'étude des voisinages de ce maximum dans les différentes directions possibles à partir de ce point : un maximum est retenu comme appartenant à la ligne d'axe du réseau si la moyenne des niveaux de gris de h points dans une direction d'investigation o est supérieure à un seuil, S; "h " caractérise l'horizon o d'investigation, et S est fonction de la moyenne, MOY, d'un coefficient correcteur s- prédéterminé dans l'image, et de l'écart type g du niveau de gris dans une zone peu vascularisée de l'image A partir d'un maximum retenu dans la phase préalable, M(i,j) de niveau N(i,j), les 8 directions d'investigation possibles sont représentées sur la figure 3 et notées Do à D7 avec leurs coefficients directeurs dans le sens i,m# (lignes) et dans le sens j,nx (colonnes). Le critère d'appartenance au squelette, calculé dans la direction D (x = O à 7) de coefficients directeurs mx et nx x est alors M(i,j) est un maximum appartenant au squelette si au moins l'un des Ex, pour x = O à 7 est supérieur au seuil x S = MOY + s A partir de ces pixels validés, la phase d'initialisation se poursuit pour la détermination de l'orientation du vaisseau auquel appartient ce point validé. Soit DxO(mO,nO) la direction qui maximise E ; il est x logique d'émettre l'hypothèse que le vaisseau suit au voisinage immédiat du pixel M(i,j) l'orientation de la droite dx ayant comme une de ses directions Dxo et comme autre direction D 8) Le suivi du vaisseau se fait alors en deux étapes, suivant les deux directions opposées de la droite dx Ces deux directions sont donc gardées en mémoire, avec le maximum validé comme appartenant au réseau. A partir d'au moins un maximum validé, l'étape suivante de la phase d'initialisation est la création d'un fichier dit de points "germes". En effet pour obtenir le squelette continu du réseau vasculaire, il est nécessaire de disposer d'au moins un point "germe" par branche. Ce fichier de points "germes" est obtenu de la manière suivante : à partir d'un maximum validé, et par une exploration systématique du périmètre d'une fenêtre centrée sur ce maximum, le traitement consiste à rechercher tous les maxima, comme dans la première étape de la phase d'initialisation ; les maxima détectés doivent ensuite être validés, en utilisant le même test par rapport au seuil. Toute structure voisine du maximum d'origine est ainsi marquée par un point "germe". Chacun de ceux-ci est ensuite pris comme centre d'une nouvelle fenêtre, etc . . Cette méthode, utilisant le caractère arborescent du réseau garantit le marquage de toutes les branches du réseau par au moins un point "germe" à partir d'un seul maximum détecté dans la première étape de la phase d'initialisation, lors de l'analyse d'une seule ligne, ou d'une seule colonne. Ce fichier de points "germes" étant créé, l'étape suivante de la phase d'initialisation consiste éventuellement à effectuer un classement hiérarchique de ces points "germes". Le but de ce classement est d'éviter l'arrêt de l'exploration d'un vaisseau important, suite à la rencontre d'une branche secondaire, alors que cette branche secondaire aurait déjà été détectée. De plus il est intéressant de permettre à l'utilisateur de faire apparaître les axes des vaisseaux principaux avant les axes des vaisseaux secondaires.Le critère de classement utilisé pour cette. hiérarchisation des points "germes" est le niveau de gris d'un voisinage centré sur chacun des points "germes". Ce voisinage a une dimension prédéterminée, c'est-à-dire qu'il comporte un certain nombre de points voisins du maximum validé, et les points "germes" sont alors classés dans un fichier en fonction du niveau de gris de ce voisinage. Ainsi, à la fin de cette opération le fichier définitif comporte une liste ordonnée de points classés selon les valeurs décroissantes des niveaux de gris des voisinages. Ces points seront les origines successives utilisées lors de la phase de suivi des vaisseaux, comme il sera expliqué cl-après. La phase d'initialisation étant terminée, la phase de suivi des vaisseaux par leurs axes, destinée à établir le fichier du squelette de l'axe vasculaire peut alors commencer. Pour cette phase, illustrée par la figure 4 on utilise la propriété suivante : l'image d'un élément de vaisseau est approximativement équivalente à la projection d'une portion de cylindre, et elle est donc caractérisée par une ligne de crête des niveaux de gris comprise entre deux flancs de pentes sensiblement égales. Pour le suivi des vaisseaux, c'est-à-dire pour passer d'un point courant au point suivant, l'analyse utilise cette propriété. A chacun des points courants, deux décisions peuvent être prises - soit l'arrêt de l'exploration si on arrive à l'extrémité d'un vaisseau, - soit la poursuite du suivi du vaisseau, ce qui impose un critère de choix de la direction de déplacement pour arriver su point suivant. Soit M(i,3)P 1 le pixel retenu a l'itération précédente, et soit M(1,3)P le pixel retenu a l'itération juste terminée: la direction (Dxo) P 1 a permis de passer du point déterminé à l'étape p-l au point déterminé à l'étape p.L'exploration du voisinage pour la recherche du point suivant du vaisseau n'est effectuée que suivant trois directions ~ la direction Dxo précédente et les directions D et Dxd respectivementà xg xd respectlveme#t à gauche et à droite par rapport à la direction précédente Dxo, et dont les coefficients directeurs sont parfaitement déterminés connaissant Dxo : Dxd(md##dJ et D#g (m g ,n g) Cette limitation dans la recherche tient compte du fait qu'il ne peut y avoir un changement de direction brutal le long d'un vaisseau d'un pixel au pixel suivant. Dans l'étape de suivi du vaisseau, le point M(i,j) ne peut donc être suivi que par ] 'un des trois points voisins situés dans ces trois directions, soit M(i+m0, j+n0), M(i+md,j+n#) et M(l+m ,j+n ), respectivement marqués M1, M2 et g g M3 sur la figure 4. Les critères suivants sont donc utilisés pour déterminer le suivi du vaisseau - critère d'arrêt de l'exploration : pour chacune des trois directions x définies précédemment, Dxo, Dxg et DXd, une amplitude moyenne des niveaux de gris d'un certain nombre de points situés dans cette direction, soit ho points (c'est-à-dire jusqu'à l'horizon ho), est calculée selon le même critère Ex que précédemment, dans l'étape d'initialisation (voir formule 1 ci-dessus). L'arrêt de l'exploration est décidé si les moyennes calculées dans les trois directions sont toutes inférieures au seuil S.L'arrêt de l'exploration est également décidé s'il existe un point déjà exploré et différent du point ## MP 1(i,) dans un voisinage de dimensions 3x3 centré sur le point MP(i,j) c'est-à-dire si on rencontre une bifurcation ou un croisement avec un axe déjà détecté. - choix de la direction : le critère retenu pour le choix de la direction, illustré par la figure 5, tient compte de la propriété des vaisseaux de se projeter suivant une ligne de crête des niveaux de gris compris entre deux flancs de pentes sensiblement égales; le critère de décision est formé d'une composante fonction des niveaux de gris de l'axe du vaisseau et d'une composante fonction de la différence entre les niveaux de gris de part et d'autre de l'axe, sur les bords des vaisseaux. Pour cela le critère est fonction des niveaux de gris des points situés dans la direction d'investigation et dans deux directions qui lui sont parallèles à une distance prédéterminée de l'axe. La figure 5 illustre pour la première direction d'investigation Dx, les deux directions orthogonales D(x+2)mod 8 et D(x+6)mod 8 et les parallèles à une certaine distance de l'axe suppose. Le critère calculé pour le choix de la direction est alors où - h est, comme précédemment, l'horizon d'investigation o dans la direction d'investigation, variable selon le niveau de gris du point courant ou d'un voisinage, c'est-à-dire fonction de l'épaisseur du vaisseau, et lo est un paramètre: 1 = i à ho; o o - m et n sont les coefficients directeurs de la direction pour laquelle le critère est calculé, - i' = i+m'.np et j' = j+n'.np ; i" = i+m".np et jll = i+m".np, où m', n' et m", n" sont respectivement les coefficients directeurs selon les deux directions orthogonales à la direction testée, - h1 est l'horizon fixe d'investigation sur les bords supposés du vaisseau de part et d'autre de l'axe à une distance np de l'axe, et 11 est un paramètre: li 11 à hi. - k est un coefficient de pondération prédéterminé. Cet algorithme est dit algorithme réalisant un suivi vectoriel des axes car il tient compte de pixels situés sur des vecteurs de même direction que la direction testée pour le vaisseau, le long de l'axe supposé et le long des bords. h o l'horizon d'exploration dans le sens de l'sxe est fonction du niveau de gris du point courant, et h1 représentant le module du vecteur sur les bords, pour vérifier la symétrie du vaisseau par rapport à l'axe, est un horizon fixe comportant, par exemple quatre à six pixels. La direction retenue, parmi les 3 possibles, est celle qui maximise le critère Cd. Cet algorithme est peu sensible au bruit et aux irrégularités dues par exemple à la mauvaise répartition du produit de contraste dans les vaisseaux. Il permet de combler des absences d'informstion dans le cas de vaisseaux fins et faiblement contrastés, et il explore la structure sur un horizon de plusieurs points, 4 à 10 points par exemple, avant de décider du déplacement d'un pas sur l'image ; le pixel sélectionné n'a pas nécessairement le niveau de gris le plus grand parmi ceux des trois pixels possibles mais, compte tenu de l'élaborstion du critère, il est celui qui est le plus probablement situé sur l'axe du vaisseau. La procédure est répétée pour la recherche du point suivant sur l'axe du vaisseau. Lorsque l'arrêt de ltexploration sur une branche est décidé, le point "germe" suivant, dans la liste ordonnée des points "germes", est pris comme base pour l'exploration d'une nouvelle branche du réseau vasculaire, si un axe déjà détecté ne se trouve pas dans son voisinage. La direction prise au début du suivi d'une nouvelle branche à partir d'un point germe est l'une de celles (opposées) mises en mémoire avec le point germe, dans la phase d'initialisation. On obtient ainsi, lorsque le fichier des points germes est épuisé, l'ensemble des pixels constituant le squelette, c'est-à-dire les points d'axes du réseau vasculaire dans l'image. Cette base de données étant constituée, la phase suivante du procédé consiste à structurer cette base de façon à utiliser au mieux les données résultant de la structure du réseau, à savoir la connexité et la variation faible de direction sur chaque branche formant le réseau. Une propriété importante du suivi vectoriel des axes est en effet, comme indiqué ci-dessus, de produire une suite de points connexes. Le parcours d'une branche ne nécessite que l'exploration d'un voisinage restreint à une fenêtre 3x3 dans la mesure où par hypothèse les axes n'ont qu'un pixel d'épaisseur et sont constitués par la ligne qui joint une suite de pixels connexes. La structuration de la base de données retenue pour la recherche des contours est une décomposition en branches. Chaque branche se termine à l'extrémité d'un vaisseau, que quele point courant n'a pas de successeur, et débute à une bifurcation, c'est-à-dire où le point courant a deux successeurs. La décomposition est un peu plus complexe que la décomposition en segments (un segment est ici un élément d'axe compris entre deux bifurcations ou croisement) : le choix du chemin à la rencontre d'une bifurcation ou d'un croisement est fonction d'un critère tenant compte de la direction relative du segment suivant et du niveau de gris du début de celui-ci. Mais, les branches étant de longueurs plus importantes que les segments, le suivi de contour sans rupture est -facilité. Le processus de décomposition en branches peut être détaillé de la manière suivante - dans un phase d'initialisation, l'origine du tronc commun est détecté par recherche dans une zone prédéterminée de l'image du squelette à niveau de gris le plus élevé (les angiographies sont prises suivant des angles connus, l'arbre sur une vue a donc toujours une position similaire); - en parcourant la- branche incluant le tronc commun, des bifurcations sont alors détectées : chacune est le début d'une branche - à la fin du parcours de la première branche la première bifurcation est reprise pour le parcours de la branche suivante et chaque nouvelle bifurcation est à nouveau mémorisée - en itérant ce processus, tout l'arbre est décomposé en branches de proche en proche. La phase suivante du procédé de segmentation d'images angiographiques vasculaires selon l'invention est alors la détection asservie des contours, par programmation dynamique, à partir des axes formés par le squelette. L'axe d'un vaisseau fournit en effet une information capitale sur les contours: l'axe se trouve nécessairement entre deux contours. La zone de recherche des contours peut donc être restreinte au voisinage du squelette. Les méthodes classiques c'est-à-dire l'application de masques, de type dérivée première ou dérivée seconde, ne sont pas suffisantes pour éviter des interruptions dans les contours. Par contre l'utilisation de la programmation dynamique garantit la continuité des contours, comme le procédé utilisé dans la phase précédente garantissait la continuité des axes. La programmation dynamique, définie en général, est une méthode algorithmique permettant par approches successives de trouver la solution optimale d'une partie du problème global, l'élargissement de la partie analysée conduisant à un nouveau chemin optimal à partir de la solution précédente. Pour la détection asservie des contours, l'algorithme de programmation dynamique repose sur la recherche d'un parcours optimal dans une matrice de gradients dont la dimension est fonction de la taille d'un segment, et qui comporte des valeurs représentatives des variations des niveaux de gris dans des directions compatibles avec celles de l'axe de la branche analysée.A partir de cette matrice de gradients, on crée une matrice de coûts dont les valeurs résultent du cumul des informations sur des chemins définis dans la matrice de gradients, en passant d'une ligne à la suivante par la voie qui maximise le coût cumulé ; puis le chemin optimal dans la matrice de coûts est détecté, il correspond au contour. La première étape de cette phase de détection de contours est donc la création des matrices de gradients : MGRAD(x,y), à partir des informations stockées pour une branche telle que définie ci-dessus, une matrice de gradients étant définie pour chaque segment de la branche Une branche est découpée en segments élémentaires de longueurs N5 pixels, NS étant variable suivant la courbure de la branche. Les coordonnées de début et de fin de chaque segment sont mémorisées. La "direction" de chaque segment est définie par les coordonnées cartésiennes de ses extrémités. Cette division en segments permet d'obtenir des zones de recherche de contours dans lesquelles la direction de l'axe varie peu et qui ont des dimensions convenables pour les calculs ultérieurs, comme expliqué ci-après. Compte tenu de la structure du réseau, une branche a son origine en un point de bifurcation détecté sur la ligne centrale du réseau. Suivant l'orientation de ce vaisseau dans l'espace à trois dimensions par rapport au plan image, la bifurcation réelle peut en fait se trouver en arrière, en avant, ou sur le côté du vaisseau réel donnant naissance à la branche. En conséquence la détection des bords, ou contours, à la naissance de la branche est particulièrement délicate et sera effectuée comme expliqué ci-après d'une manière particulière. Indépendamment de ce problème d'initialisation, la matrice de gradients comporte les paramètres suivants : du fait que le contour est toujours sensiblement parallèle à l'sue du vaisseau, le gradient est fonction de la variation du niveau de gris d'une part, et de la direction du squelette d'autre part. Pour répondre à la première condition, le calcul du gradient est effectué sur une espace élémentaire ayant trois pixels de côté. La direction de ce gradient peut varier dans des proportions importantes notamment du fait de la présence de bifurcations ou de croisements de vaisseaux. Les valeurs de ce paramètre ne sont alors plus représentatives des variations des niveaux de gris induites par la branche étudiée mais sont en fait une combinaison de deux vecteurs - un vecteur V1 lié à la branche en cours d'étude - un vecteur V2 représentatif du contour du vaisseau extérieur à cette branche. La dimension du domaine utilisé pour le calcul du gradient (3 pixels x 3 pixels) rend illusoire la recherche d'une grande précision sur la direction du gradient obtenu, et donc de la contribution exacte du vecteur lié à la branche étudiée. Il est cependant intéressant de réduire l'influence du vecteur V2 qui tend à faire diverger le contour- par rapport à l'axe. Pour résoudre ce problème, les valeurs de la matrice de gradient sont les valeurs algébriques des projections des gradients sur la perpendiculaire au segment représentatif de l'axe. La figure 6 montre l'axe du vaisseau et un segment AB sur cet axe, le gradient G au point M(l,j), et la projection GP de ce gradient sur la perpendiculaire à l'axe du vaisseau, dont la valeur algébrique est mise en mémoire dans la matrice de gradients. Ce paramètre étant défini, une autre étape préalable consiste à définir la zone de recherche du contour, à partir d'un segment donné. Cette zone de recherche est différente suivant l'inclinaison du segment analysé dans le plan de l'image. En effet cette zone est définie à partir des composantes du segment et à partir de la largeur maximale d'un vaisseau du réseau vasculaire en cours d'étude.Les figures 7a et 7b montrent les deux configurations possibles des espaces de recherche du contour autour du segment élémentaire analysé: sur la figure 7a, iB-3A est supérieur à iB-iA, alors que sur la figure 7b, iB#iA est supérieur 3B-iA. Si EC est l'écart maximal entre l'axe et le contour, c 'est-à-dire la largeur de l'espace de recherche d'un côté du segment, et si LG est la longueur de l'espace de recherche, les coordonnées x et y de la matrice de gradient MGRAD(x,y) sont telles que O La matrice de gradients MGRAD(x,y) contient autant de valeurs qu'il y a de points dans chaque zone de recherche de contours, et les coordonnées (i,j) des pixels centres des espaces élémentaires (3x3) de calcul des gradients dans le plan image sont stockés. Pour chaque côté du segment, on dispose d'une matrice de gradients de LG lignes et de EC colonnes contenant des valeurs fonctions de la probabilité de présence d'un contour au point (x,y) correspondant aux coordonnées (i,3) de l'image. Le problème à résoudre est alors de définir une méthode permettant de déterminer le contour le plus probable. La programmation dynamique permet de décomposer la recherche du meilleur chemin au travers de la matrice de gradients en une série de processus décisionnels élémentaires donnant le pas à effectuer . pour passer d'une ligne de la matrice à la suivante. La première étape de l'algorithme consiste donc å construire une matrice de coûts, MCOUT(x,y), à partir de la matrice de gradient MGRAD(x,y). Cette construction est effectuée de la manière suivante - des chemins élémentaires sont établis entre la première et la seconde ligne de la matrice de gradients avec une déviation latérale maximale autorisée de l pixel - chaque point de la deuxième ligne de la matrice de coûts prend la valeur cumulée correspondant au chemin de plus fort coût entre les trois chemins possibles - le processus est répété de la ligne 2 à la ligne 3 etc.. et ce jusqu'à la dernière ligne de la matrice de gradient de cet espace de recherche. L'étape suivante est alors la construction du chemin optimal d'un élément de contour. Cette construction se fait dans le sens inverse du cumul des valeurs dans la matrice de coûts - à partir de la dernière ligne de la matrice de coûts on cherche le maximum correspondant au chemin optimal - la recherche du chemin jusqu'à la ligne précédente se fait en sélectionnant de proche en proche le maximum parmi les trois points les plus proches - la procédure est renouvelée jusqu la première ligne de la matrice de coûts. Elle est arrêtée si l'on atteint la colonne extérieure de la matrice de coûts - les coordonnées de fin de chemin optimal sont mémorisées. Les opérations décrites ci-dessus permettent de rechercher, dans une zone de recherche prédéfinie, le chemin optimal d'un élément de contours. Il reste alors une étape ultérieure qui consiste à relier entre eux les éléments de contours ainsi détectés. Cette opération est le chaînage des éléments de contours. Plusieurs cas sont à distinguer selon le rang du segment dans la branche. Pour simplifier l'exposé on décrit d'abord ci-après le cas d'un segment intermédiaire de la branche La division d'une branche en segments ordonnés avait pour but de réduire les dimensions du graphe de recherche du chemin optimal et de faciliter le maintien de la contrainte sur les positions relatives des axes et des contours.En effet si les segments sont trop longs, c'est-à-dire si la décomposition est insuffisante, l'hypothèse selon laquelle un élément de contour existe de part et d'autre de l'axe n'est plus vérifiée obligatoirement, notamment lorsque la direction du vaisseau varie notablement ; si le segment est trop long la direction du segment, seulement définie par ses extrémités, serait telle que les deux contours se trouveraient sur une même zone de recherche. Mais, du fait de cette décomposition, il est nécessaire alors de "chaîner" les contours sans créer de discontinuités aux points de raccordement. Du fait de la méthode de détection des axes constituant le squelette du réseau vasculaire, tout segment a un successeur excepté si l'on se trouve à l'extrémité d'un vaisseau. Par conséquent les coordonnées de fin d'un segment correspondent aux coordonnées de début du segment suivant. De ce fait, il devra en être de même pour les contours, un élément de contour ayant nécessairement une coordonnée de fin identique à la coordonnée de début de l'élément de contours correspondant au segment suivant. Le cas le plus simple pour le chaînage des segments et des contours correspondants, est représenté sur la figure 8: les zones de recherche de contours sont jointives du fait que leurs largeurs sont parallèles (sur la figure 8 dans le sens de l'axe des i). Dans ce cas pour les segments successifs Si et le ledébut du contour associé à une zone de recherche d'un côté de Si+1 est la fin du contour associé à la zone de recherche précédente du même côté de Si : DSi+1 = FSI et il en est de même pour les zones de recherche de l'autre côté des segments. Le cumul des valeurs dans la matrice de coûts, MCOUT(x, y) se fait normalement en sens inverse du parcours de la branche, c'est-à-dire en partant de l'extrémité de la branche. Ainsi, sur la figure 8 le cumul des valeurs se fera de FS. à DS. (de la fin au début du segment), le parcours de la I i branche étant effectué de DS. à FS ; il en est de même pour le I i segment Si+1 (cumul de FSi+1 à DSi+i et parcourt de DSi+1 à FS1+1). La dernière ligne de la matrice MGRADi+i liée au segment Si+1 correspond à la verticale passant par DSi+1. Compte tenu de la procédure utilisée, le maximum sur la dernière ligne de la matrice MCOUT. ne se trouve pas nécessairement confondu avec la fin du contour associé au segment précédent, FC.. Le maximum sur la ligne ymax-l de MCOUTi+1 peut également être à plus de 1 point à gauche ou à droite de FC.. Mais la continuité du contour est privilégiée : le chemin dans la matrice de coût caractérisant le contour peut donc pour deux segments successifs être légèrement différent du chemin optimal, la contrainte de connexité des points du chemin étant maintenue sur tout le parcours. Par ailleurs lorsque la direction du segment de varie par rapport à celle du segment précédent Si de telle façon que l'on passe de #i > bj pour Si, à a i #j), les deux zones de recherche passant par les extrémités des segments ne sont plus jointives limites parallèles å j pour Si et parallèles à i pour Si+1. Dans ce cas - le point de fin de contour associé à Si, FCi n'est pas sur la première ligne de la matrice de gradients MGRADi, qui est la ligne passant par FSi, - d'un côté du segment, la position du point de contour associé à Si, FCi, début du contour associé à Si+1, DCi+1, entraîne la réduction du nombre de lignes de MGRADî+i de L LG = (jmax-jmin)+AL avec bL - Par contre, de l'autre côté, la position du point de fin de contour associée à Si, F'C1, nécessite une augmentation du nombre de lignes de MGRAD' d'une valeur fl L', pour que F'C. soit inclus dans la zone de recherche associée au segment suivant Si+1: LG = (jmax-jmln) + ss L avec tL > O. Pour résoudre ce problème de continuité des contours, les perfectionnements suivants sont apportés - Lors de la phase de cumul des valeurs dans MCOUT si, à la ligne FS. , il n'y a pas de point suffisamment consistant avec la présence d'un contour, cette ligne est rejetée. Il en est de même pour les lignes suivantes jusqu'à la rencontre d'une valeur de gradient projeté supérieure à un seuil : toutes les lignes suivantes sont alors prises en considération. - Lors du parcours du chemin, si celui-ci approche l'extrémité d'une ligne (x = xmax - 1) dans MCOUT(x,y), le contour est considéré comme sortant du champ du fait d'une courbure. Le dernier point validé est pris comme point terminal du contour Ci,FCi et comme point initial du contour Ci+1,DCi+l, et ce point définit également la dernière ligne de la matrice MGRAD. (comme représenté sur la figure 8). La description ci-dessus explique le suivi des contours pour des segments quelconques des branches à l'exception des premiers. En effet leurs origines sont des bifurcations dont les positions sont définies de manière imprécise par le suivi des axes lors de l'obtention du squelette ; le début réel des contours peut être soit au bord du vaisseau donnant naissance à une branche, soit superposé à une partie de ce vaisseau, ceci suivant les directions relatives du vaisseau dans ltespace 3D et dans le plan de projection. En conséquence la détection des contours est plus fiable à une certaine distance des bifurcations. Pour cette raison, et seulement pour le premier segment de chaque branche, le procédé ci-dessus est conservé, à l'exception du sens de cumul des valeurs et du sens de parcours du chemin dans la matrice de coût qui sont inversés par rapport à ceux décrits ci-dessus, comme indiqués sur la figure 10 qui représente une branche Bn avec son axe et ses contours et une branche Bm prenant naissance au point de bifurcation Bif sur l'axe de Bn : pour S1 le cumul des valeurs se fait à partir de vers la fin du segment FS i et cette opération ne commence qu'à partir de la première ligne où MGRAD contient au moins une valeur de gradient projeté supérieure à un seuil et donc caractérisant éventuellement un point de contour.Par contre le parcours du chemin dans MCOUT est effectué à partir du début de contour associé à S 2,DC2 même si celui-ci ne coïncide pas exactement avec FC1 (maximum du cumul), pour les raisons exposées ci-dessus. On obtient ainsi du fait de la méthode des contours continus directement liés aux axes définis dans la phase préalable. Une étape finale éventuelle peut alors consister à corriger la position des axes en fonction des contours ainsi détectés, notamment si dans un traitement d'images ultérieur les données relatives aux axes sont utilisées, par exemple pour une reconstruction en trois dimensions du réseau vasculaire à partir d'images de projection, les données utiles de ces images de projection étant limitées aux axes. Pour cette étape finale, il suffit de décaler l'axe d'un nombre de pixels tels qu'il se situe en position médiane par rapport aux deux contours précédemment définis. Un tel procédé réalise une segmentation fiable et très efficace des images radiographiques pour toutes les applications d'analyse de réseaux vasculaires, par hypothèse continus. REVENDICATIONS 1. Procédé de segmentation d'images angiographiques numérisées pour l'extraction des axes et des contours de réseaux vasculaires, caractérisé en ce-qu'il comporte - une première phase au cours de laquelle un squelette du réseau est extrait par une méthode de suivi vectoriel des lignes centrales des vaisseaux, ou axes, recherchant, à partir de points isolés reconnus comme appartenant au réseau par vérification d'un critère Ex et affectés chacun d'une direction, le ou les points connexes suivants, dans un nombre limité de directions prédéfinies fonction de la direction affectée au point courant, et qui satisfont un critère prédéfini Cd fonction des niveaux de gris d'un ensemble de points prédéfinis, - et une phase de détection asservie des contours s'appuyant sur le squelette du réseau obtenu dans la phase précédente pour rechercher, de part et d'autre des axes, un contour également formé de suites de points connexes, correspondant à des points de transition importante du niveau de gris, recherchés dans des zones prédéfinies de l'image. 2. Procédé selon la revendication 1, caractérisé en ce que la phase d'extraction du squelette comporte une phase d'initialisation au cours de laquelle est constitué un ensemble de points germes marquant de façon certaine toutes les branches du réseau, ces points étant reconnus comme appartenant au squelette du réseau si, dans au moins une direction, définie par le point testé et l'un de ses voisins, le critère Ey est où h, horizon d'investigation paramétré par 1, est égal à un nombre prédéfini de pixels dans la direction d'investigation D x de coefficients directeurs mx et nx, i et j sont les coordonnées du point courant M(i,j) dans l'image, et N(i,j) son niveau de gris ; la direction qui maximise Ex étant affectée au point connexe du point courant dans cette- direction reconnu comme appartenant au squelette du réseau. 3. Procédé selon la revendication 2, caractérisé en ce que les points testés pour la constitution de l'ensemble des points germes sont, dans une phase initiale des pixels correspondant à des maxima de niveau de gris le long d'une droite dans l'image, puis les pixels correspondant à des maxima de niveau de gris ~ sur le pourtour de fenêtres successivement centrées sur les points germes précédemment retenus, toutes les branches du réseau étant ainsi marquées au moins d'un point germe. 4. Procédé selon la revendication 3, caractérisé en ce qu'un point courant est retenu comme correspondant à un maximum de niveau de gris le long d'une portion de ligne ou de colonne d'image pour un maximum du niveau local de gris, ce niveau local étant défini comme la moyenne des niveaux absolus, de gris de trois pixels consécutifs sur la ligne ou la colonne, éventuellement diminué du niveau moyen d'une zone peu vascularisée de l'image, ce niveau local étant croissant sur un nombre prédéterminé de pixels consécutifs, puis devenant décroissant après le pixel retenu. 5. Procédé selon l'une des revendications 2 à 4, caractérisé en ce que l'ensemble des points germes est classé par valeurs décroissantes des niveaux de gris de voisinages centrés sur ces points germes, pour former un fichier ordonné des points germes. 6. Procédé selon la revendication 1, caractérisé en ce que le suivi des axes pour l'obtention du squelette est réalisé en sélectionnant, à partir d'un point précédemment retenu, seulement trois directions de recherche possibles : la direction Do afffectée au point courant précédent et les deux directions les plus voisines de part et d'autre de la direction Do, soient Dg et Dd, le point connexe suivant s'il existe étant nécessairement dans l'une de ces trois directions, la suite des points connexes étant interrompue par reconnaissance d'une est inférieur à un seuil S pour les trois directions de recherche possibles, où h est l'horizon d'investigation égal à o un nombre prédéfini de pixels dans la direction d'investigation D de coefficients directeurs mx et n I i et j les coordonnées x du point courant et N(l+mxl,i+nxl) le niveau de gris des pixels successifs de coordonnées (i+m#l,j+n#l). 7. Procédé selon la revendication 6, caractérisé en ce que la poursuite du suivi, lorsque le critère E est supérieur x au seuil S pour au moins l'une des trois directions, conduit à la sélection du point connexe suivant parmi les trois points possibles, en utilisant un critère de choix fonction d'une part de la persistance d'un niveau de gris important caractérisant l'axe d'une crête de vaisseau dans la direction d'investigation et d'autre part d'une symétrie suffisante de part et d'autre de la direction d'investigation, sur les bords des vaisseaux, parallèlement à la direction d'investigation. 8. Procédé selon la revendication 7, caractérisé en ce que le point connexe suivant est retenu dans la direction qui maximise le critère Cd tel que h étant l'horizon d'investigation paramétré par 1 dans la o o direction d'investigation de coefficients directeurs m et n, h étant un horizon fixe paramétré par 11; (m',n') et (m",n") étant les coefficients directeurs des directions orthogonales à la direction d'investigation, et i', j', i", j" étant tels que i' = i+m'.np ; i" = i+m".np ; jt = j+n'.np ; j" = j+n".np où np est l'écart en nombre de pixels entre la ligne d'axe et les lignes des bords, i et j étant les coordonnées du point de départ du test. 9. Procédé selon la revendication 1, caractérisé en ce que pour chacun des points du squelette du réseau, une base de données comporte ses coordonnées, le niveau de gris calculé par la moyenne des niveaux dans un voisinage centré sur ce point, et la direction qui lui est affectée, et en ce que cette base de données est structurée en branches définies entre une extrémité et une origine pour la branche principale et entre une extrémité et un point de bifurcation pour les autres branches, le critère de détermination de la continuation d'une branche à un point de bifurcation ou de croisement tenant compte des directions affectées aux points suivants par rapport à la direction précédente, et des niveaux de gris de ces points suivants. 10. Procédé selon la revendication 9, caractérisé en ce que la détection asservie des contours s'appuyant sur le squelette du réseau met en oeuvre une programmation dynamique pour la recherche de ces contours dans des zones de dimensions variables situées de part et d'autre des axes du réseau. 11. Procédé selon la revendication 10, caractérisé en ce que, pour la mise en oeuvre d'une programmation dynamique, chaque branche est divisée en segments de longueur variable suivant la courbure locale de la branche, et de coordonnées de début et de fin connues, - en ce que les zones de recherche de part et d'autre de l'axe ont comme longueur la longueur des segments et comme largeur une largeur variable en fonction du niveau de gris moyen du segment et supérieure à l'écart maximum entre l'axe et le contour, - en ce qu'une matrice de- gradients représentatifs des varaitions de gris des points, ou de leurs voisinages, relativement à la direction du segment est établie pour chaque zone de recherche, - en ce qu 'une matrice de coûts résultant du cumul des valeurs de gradients, le long des chemins élémentaires d'une ligne de la matrice de gradients à la suivante, selon le chemin de plus fort coût et avec une déviation latérale maximale autorisée de un point est établie, - et en ce que le chemin optimal est suivi à partir du maximum de la dernière ligne de la matrice de coûts jusqu la première ligne, en sélectionnant à chaque étape le maximum parmi les trois points les plus proches de la ligne précédente, la procédure étant arrêtée éventuellement avant la première ligne lorsque le point se trouve appartenir à la limite de la zone de recherche. 12. Procédé selon la revendication 11, caractérisé en ce que pour un segment quelconque d'une branche, à l'exception d'un premier segment de branche ayant pour origine un point de bifurcation, le sens du cumul dans la matrice de gradients pour construire la matrice de coûts part de l'extrémité du segment pour aller vers son origine, le sens de construction du chemin optimal caractérisant le contour partant de la ligne de la matrice associée au début du segment pour arriver à la ligne de la matrice associée à son extrémité, ces sens de cumul et de recherche du chemin optimal étant inversés pour les premiers segments des branches pour tenir compte des incertitudes liées aux bifurcations.