Brevet FR2773247A1 — publié le 19990702

"PROCEDE DE TRAITEMENT D'UNE IMAGE MULTIDIMENSIONNELLE BRUITEE ET APPAREIL D'IMAGERIE MEDICALE METTANT EN OEUVRE CE PROCEDE" Description: L'invention concerne un procédé de traitement d'image pour réduire le bruit en conservant les bords d'objets représentés dans une image d'intensité, comprenant une étape d'acquisition de l'image d'intensité sous forme d'une matrice multidimensionnelle de points. L'invention concerne également un appareil d'imagerie médicale ayant un système pour mettre en oeuvre ce procédé. Par points, dans une matrice multidimensionnelle, on entend des pixels pour une image bidimensionnelle et des voxels pour une image tridimensionnelle. L'invention trouve son application en particulier dans les systèmes médicaux à rayons X fournissant des images multidimensionnelles. Une méthode de réduction du bruit dans une image numérique bruitée est déjà connue de la publication intitulée "Image Processing: Flows under Min/Max Curvature and Mean Curvature" par R. Malladi et alii, dans "GRAPHIC MODELS AND IMAGE PROCESSING, Vol.58, N 2, March, pp127-141, 1996. ARTICLE N OOl1. Ce document décrit une classe d'algorithmes appropriés à réduire le bruit dans une image et à en rehausser le contraste. Le procédé est applicable aux images bruitées en noir et blanc, ou en couleur. Le procédé comprend deux étapes principales. Premièrement, une formulation d'un ensemble de niveaux de gris ou de niveaux d'intensité est utilisée pour engendrer des courbes joignant les pixels de même intensité dans l'image appelées lignes isophotes. Deuxièmement des moyens sont utilisés pour lisser les irrégularités de courbure des lignes isophotes, ce qui diminue le bruit et rehausse les images. L'opération de lissage des lignes isophotes est menée à bien par des réductions des minima et des maxima de courbure des lignes isophotes, appelées commutation de courbures. Les degrés de réduction de courbure sont sélectionnés dans une fenêtre de manière à diminuer les irrégularités du dessin des lignes isophotes dans des limites. Il en résulte que le bruit est diminué dans cette fenêtre de commutation. Ce procédé est préconisé selon le document cité d'abord parce qu'il comprend seulement un paramètre de rehaussement choisi automatiquement dans la plupart des cas; et ensuite par ce qu'il arrête de lisser un point automatiquement en dépendance des dimensions de la fenêtre de commutation. En fait, ce procédé repose sur des étapes itératives. Une étape itérative est relancée en boucle autant de fois qu'il est nécessaire pour obtenir un lissage correct des lignes isophotes. Pour une image complexe, c'est-à-dire non géométrique, et tridimensionnelle, c'est-à-dire par exemple de 1283 voxels, il a été trouvé par des essais faits par ailleurs, la nécessité de réaliser 300 à 500 itérations, ce qui conduit, dans l'état actuel de la technique à des temps de calcul de l'ordre de 10 heures par image. Un but de l'invention est de fournir un procédé pour éliminer le bruit, c'est-à-dire lisser une image complexe multidimensionnelle, avec un temps de calcul beaucoup plus court, c'est-à-dire qui ne nécessite qu'un faible nombre d'itérations. Ce but est atteint au moyen d'un procédé de traitement d'image pour conserver les bords d'objets représentés dans une image d'intensité et pour réduire le bruit, comprenant des étapes de: acquisition de l'image d'intensité sous forme d'une matrice multidimensionnelle de points, détermination du gradient d'intensité en chaque point, défini par sa direction et son module, passage en chaque point d'un filtre spatial récursif à noyau anisotrope provoquant un degré de lissage plus important perpendiculairement à la direction du gradient que dans la direction du gradient d'intensité. Un avantage du procédé selon l'invention est qu'il est récursif, donc il peut être mis en oeuvre avec une grande économie de calculs parce que, à chaque passage du filtre en un point, des éléments de calcul déjà précédemment effectués sont utilisés. Ainsi, il peut être mis en oeuvre en un temps plus court qu'un procédé itératif. Par exemple, avec les techniques actuelles, le présent procédé nécessite moins d'une minute pour traiter un volume de 1283 voxels. De ce fait, il est beaucoup plus facile à mettre en oeuvre et environ 1000 fois plus rapide que le procédé connu. Un autre avantage du présent procédé est qu'il est adaptatif, c'est-àdire que le filtrage s'adapte automatiquement en chaque point au type de structure rencontrée. Ce procédé applique donc un filtrage différent aux points du fond et aux points des bords d'objets. Donc ce procédé permet de lisser l'image en conservant des bords d'objets bien nets, non déplacés, non déformés, et débarrassés du bruit. Ce procédé est précis et performant. Un autre avantage du présent procédé est qu'il concerne un filtre spatial. Il ne nécessite donc pas d'autres données que celle de la seule image traitée, et en outre il est robuste. Comme il est très rapide, il peut néanmoins être utilisé dans le traitement d'une séquence d'images. Un autre avantage du présent procédé est qu'il est complètement automatique. L'invention est décrite ci-après en détail, en référence aux figures schématiques annexées, dont la liste est la suivante: les FIG. 1A et lB qui représentent respectivement une image tridimensionnelle bruitée et filtrée; la FIG.2 qui représente sous forme de blocs les étapes d'un procédé de filtrage du bruit dans une image tridimensionnelle ou bidimensionnelle; la FIG.3A qui représente sous forme de blocs les étapes générales de détermination du facteur de lissage a, et les FIGs.3B et 3C qui illustrent en détail les déterminations du facteur de lissage a pour traiter respectivement une image tridimensionnelle et une image bidimensionnelle; la FIG. 4A qui représente les orientations discrétisées ei du vecteur gradient pour une image tridimensionnelle dans un repère et la FIG.4B qui représente les droites de balayage,'j 4,j discrétisées utilisées pour filtrer l'image tridimensionnelle; les FIGs.5A, 5B, 5C qui montrent dans trois cas les rapports entre les directions discrétisées des gradients et les droites de balayage perpendiculaires; la FIG.6A qui représente les orientations discrétisées Ei du vecteur gradient dans le plan de l'image pour une image bidimensionnelle et la FIG.6B qui représente les droites de balayage discrétisées,'j $4'j utilisées pour filtrer l'image bidimensionnelle; les FIGS.7A et 7B qui illustrent les balayages causal et anticausal pour filtrer une image bidimensionnelle; les FIGS.SA et 8B qui illustrent les noyaux des filtres Fa. à dimensions variables, utilisés pour filtrer une image tridimensionnelle; la FIG.9 qui montre un appareil à rayons X avec des moyens pour mettre en oeuvre une méthode de traitement d'image. PROCEDE Les figures sont décrites ci-après en détail en relation avec les différentes étapes d'un procédé de traitement spatial récursif complètement automatique d'une image numérique bruitée multidimensionnelle, pour diminuer le bruit tout en conservant avec une grande précision les détails fins et les bords des objets représentés dans cette image. Le procédé vise particulièrement à effectuer un fort lissage des bords pour mieux distinguer les bords du fond de l'image. Le procédé vise donc à appliquer un fort degré de filtrage dans les directions des bords d'objets et un plus faible degré dans les directions perpendiculaires à ces bords. Par détails fins, on entend des détails d'objet ayant une dimension de l'ordre de quelques pixels, par exemple inférieure ou égale à 10 pixels. Par multidimensionnelle, on entend une image à 2 ou 3 dimensions spatiales, et à l'occasion 1 dimension temporelle supplémentaire si l'image est animée. Le procédé trouve une application particulière dans le traitement d'images vasculaires obtenues par le moyen d'un appareil aux rayons X travaillant en mode fluoroscopique. Dans ce mode, les doses de rayons X pour former l'image sont très faibles d'o il résulte que l'image est très bruitée. Les images vasculaires représentent des objets tels que des artères qui sont distingués du fond surtout par leur lignes de bord, des très fins vaisseaux tels que ceux du cerveau, et des outils fins tel qu'un cathéter qui peuvent être assimilés à un bord. Le bruit dans ces image fait que ces bords sont entachés de rugosités, recouverts de petites structures avec des courbures atypiques qui rendent difficile l'exploitation de ces images par un praticien qui veut visualiser en temps réel une opération pratiquée sur un patient. Le présent procédé vise à éliminer ces structures et rugosités pour révéler les bords des objets ou les objets fins dans toute leur netteté. La FIG.1A montre une image numérique JO dans une application représentant une artère en trois dimensions (3-D) envisagée à titre d'exemple. Cette image J0. appelée image d'origine, est très bruitée. Notamment, le bruit forme des taches empêchant un praticien de bien distinguer les limites entre l'artère et le fond. Cette artère, dans l'exemple traité, a une excroissance latérale appelée anévrisme dont le praticien veut diagnostiquer l'importance, ou que le praticien veut éliminer par une opération dans le but d'éviter un problème appelé rupture d'anévrisme. Dans une telle image médicale, il est très important de supprimer le bruit, sans supprimer des détails de l'objet, afin que le praticien puisse finalement opérer l'artère avec une grande précision, c'est- à-dire obturer l'anévrisme sans endommager l'artère elle-même. La FIG.lB représente l'image J débruitée par la méthode de l'invention, correspondant à l'image bruitée J0 de la FIG.1A. Dans le cas précité, il est important de pouvoir traiter l'image multidimensionnelle en un temps très court, afin que le praticien puisse suivre en direct l'évolution de la région de l'artère sous l'action des gestes opératoires. Cette technique de suivi d'une opération au moyen d'un affichage simultané d'une image, ou d'une séquence d'images, permet de réaliser des gestes opératoires de manière non invasive, c'est- à-dire sans ouvrir le patient, par le moyen d'un cathéter introduit dans le système artériel du patient et dont la progression est suivie sur une image de rayons X en mode fluoroscopique. Cette technique permet de traiter des anévrismes du cerveau sans avoir recours à une trépanation qui serait nécessaire pour avoir une vision directe de la zone à opérer. Ceci évite de nombreux problèmes post-opératoires. Il est également important de pouvoir traiter l'image par un procédé automatique. Par procédé automatique de traitement d'image, on entend un procédé de construction d'une image débruitée à partir d'une image bruitée, sans intervention humaine, par l'action autonome d'un système de mise en oeuvre du procédé. Le procédé selon l'invention est décrit ci-après pour le traitement d'une image tridimensionnelle 3-D, ainsi que pour le traitement d'une image bidimensionnelle 2-D. On appelle voxels les éléments volumiques des images tridimensionnelles et pixels les éléments des images bidimensionnelles. Ce procédé comprend en référence avec la FIG.2 les étapes principales de: 1) Acquisition d'une image numérique. En 3-D: Une image numérique J0 est acquise par des moyens 100 de formation d'image numérique sous la forme d'une matrice tridimensionnelle de voxels P(x,y,z), caractérisés par leurs coordonnées x,y,z selon trois axes de coordonnées X,Y,Z relatifs à l'image, et par leur intensité numérisée I(x,y,z), par exemple codée en niveaux de gris. En 2-D: Une image numérique est acquise par les moyens de formation d'image numérique sous la forme d'une matrice bidimensionnelle 2-D de pixels P(x,y) caractérisés par leurs coordonnées x,y selon deux axes de coordonnées X,Y relatifs à l'image et par leur intensité I(x,y) numérisée par exemple codée en niveaux de gris. 2) Evaluation des gradients d'intensité. En 3-D: Une détermination des gradients d'intensité est d'abord effectuée par des moyens 300 qui réalisent des calculs des dérivées d'ordre 1 de la fonction d'intensité I(x,y,z) dans l'image 3-D bruitée J0, selon par exemple les axes orthogonaux X,Y,Z, en chaque voxel courant P(x,y,z) à traiter. Pour réaliser cette étape de manière automatique, de préférence, l'image J0 à traiter est soumise à un balayage systématique et automatique 200, voxel après voxel, selon des droites de balayage parallèles aux axes X,Y,Z. Au cours de ce balayage systématique et automatique 200, l'intensité I(x,y,z) de chaque pixel P(x,y,z) est soumise à un filtrage 300 au moyen de filtres orientés connus sous le nom de "Sovel" et notés ci-après Sovelx selon l'axe X, Sovely selon l'axe Y et Sovelz selon l'axe Z. La mise en oeuvre des filtres de Sovel est décrite dans un ouvrage d'enseignement général intitulé "Computer Vision" par Dana H. Balard et Christopher M. Brown, publié chez PRENTICE-HALL, Inc., Englewood, Cliffs, New Jersey. US, 07632 p.77, 2.6-17 en relation avec la FIG.3.10C. Le résultat des filtres de Sovel permet de déterminer dans une étape 400 en chaque voxel courant P(x,y,z) de l'image J0 tridimensionnelle un vecteur gradient d'intensité Gp caractérisé par son module noté IIGpll et par la direction Cp de son support repéré dans le système orthonormé X,Y,Z. En référence à la FIG.4A, une discrétisation 500 de la direction (p du support du gradient est ensuite effectuée selon des directions O p. Dans le traitement de l'image tridimensionnelle, des moyens 500 associent chaque direction Pp du vecteur gradient à une direction discrète notée p fixée dans un repère tridimensionnel noté CFC dont le centre O est placé en coïncidence avec le voxel courant P(x,y,z). Le repère CFC est un système défini A PRIORI et par exemple stocké dans une mémoire 640. Le repère CFC forme un système cubique à faces et arêtes centrées dont les axes principaux sont parallèles aux axes X,Y,Z de la matrice de l'image et présente une symétrie autour de son centre 0. On définit dans ce repère CFC un nombre de 13 axes Op directs passant par le centre O du système cubique et respectivement par quatre sommets 1, 3, 7, 9, soit 01, (3, 0, 09; par le centre 0 et par trois centres 5, 11, 13 des faces centrées soit 05 11' O13; et par le centre O et par six milieux d'arêtes 2, 4, 6, 8, 10, 12 soit 02' Q4' 06' 08' 10 et 12. On définit également dans ce repère CFC un nombre de 13 axes 0' P indirects: 0' à 0'13 qui ont des indices correspondant aux axes directs quand ils sont situés dans le prolongement des axes directs et qui sont symétriques des axes directs par rapport au centre 0. On dispose donc dans ce repère CFC de 26 directions discrètes fixes. Toute direction (p de vecteur gradient Gp lié à un voxel courant P(x,y, z) de l'image tridimensionnelle est remplacée au cours de l'opération de discrétisation 500 par l'une des 26 directions portées par les axes p ou 0' en sélectionnant la plus proche dans l'espace. A l'issue de l'opération 500, tout vecteur gradient Gp est donc discrétisé, c'est-à-dire caractérisé par les coordonnées x,y,z du pixel courant P(x, y,z) auquel il s'applique, par son module IGp et par une direction p ou 0' de 1 à 13 ou de 0'1 à 0'13. p p 1 3 En 2-D: En référence à la FIG.2, dans l'image bidimensionnelle, la détermination des gradients d'intensité est effectuée par les moyens 300 qui réalisent des calculs des dérivées d'ordre 1 de la fonction d'intensité I(x,y) dans l'image bruitée selon les deux axes de coordonnée X,Y de la matrice de l'image. Cette étape peut être réalisée de manière automatique en balayant l'image systématiquement pixel après pixel selon des droites de balayage parallèles aux axes X,Y par une opération 200. Au cours de ce balayage systématique et automatique l'intensité I(x,y) au pixel courant est soumise au filtrage 300 au moyen de filtres de Sovel déjà décrits et notés Sovelx selon l'axe X et Sovel selon l'axe Y. Ce filtrage détermine en chaque pixel un Y vecteur gradient Gp caractérisé dans une étape suivante 400 par son module IIGpII et par la direction cp de son support repérée dans le système d'axes X,Y. Comme illustré par la FIG.6A, par l'application des filtres Sovelx et Sovely, on obtient selon les axes X et Y, des composantes GX et Gy, du vecteur gradient d'intensité Gp au pixel courant P(x,y). Le calcul du rapport de ces composantes GX et Gy fournit la valeur de la tangente d'un angle (p que fait ce vecteur gradient Gp par exemple avec l'axe X pris comme référence. Par le calcul de la somme des projections des composantes Gx et Gy sur la direction Wp, on déduit également le module 1JGpl| de ce vecteur gradient. On définit ainsi, en chaque pixel courant P(x,y), le vecteur gradient par la direction (p et par son module Gdpll. Une discrétisation 500 de l'orientation yp du gradient d'intensité est ensuite effectuée selon des directions O. p En référence à la FIG.6A, dans le cas d'une image à deux dimensions (2-D), la discrétisation des directions cp des supports des vecteurs gradients est effectuée par les moyens 500 qui associent la direction Cp du vecteur gradient à une direction discrète notée Op fixée dans le repère bidimensionnel orthogonal X,Y dont le centre O est placé en coïncidence avec le pixel courant P(x,y). Dans ce repère bidimensionnel X,Y, déterminé A PRIORI et stocké dans une mémoire 640, les directions p passant par O sont espacées régulièrement avec un pas r/M o M est une constante. Par exemple, la constante M est choisie égale à 8 et le pas d'angle est noté AO = n/8. Comme montré sur la FIG.6A, la direction associée au module du gradient |1Gp|| est alors non plus pp, mais une direction Op ou O' qui est la direction discrète choisie la plus proche de cp, p par exemple 02. Dans le cas d'un pas d'angle AO =,/8, l'angle de la direction O ou &' associée au gradient prend les valeurs: 01 =; 02 =. /8; 03 = n/4; 04 = 3n/8; 05 = r/2; 6 = 5/8; 7 = 3./4; 08 = 7r/8. 3) Lissage de l'image En référence à la FIG.2, le lissage de l'image pour diminuer le bruit est effectué par filtrage spatial récursif. 3.a) Opération de balayage de l'image à traiter On définit A PRIORI un système de balayage automatique de l'image. Ce système est stocké par exemple dans la mémoire 630. Ce système de balayage a l'avantage de résoudre plusieurs problèmes. En 3-D et en 2-D: En référence respectivement aux FIGs.4B et 6B, ce système de balayage est défini par les mêmes moyens que le système de discrétisation des angles des gradients. Ainsi, respectivement en 3-D et en 2-D, on définit dans le repère CFC et dans le plan de l'image, des ensembles de directions de droites de balayages parallèles entre elles, notées I'j i4,j. La direction de chaque ensemble de droites de balayage est alors respectivement parallèle à une direction discrète Oj, O' des gradients. Du fait que les directions de discrétisation des gradients O. et O'j sont décalées de 180 , ces deux directions EOj et O'. sont en coïncidence avec une même droite de balayage 14' 41j. Les directions des droites de balayage sont en outre choisies de même nombre que les directions de discrétisation des gradients et sont réparties dans l'espace avec le même pas. Ces ensembles de droites de balayage sont parcourues pour effectuer le balayage des images 3-D et 2-D de manière telle que chaque voxel ou pixel est visité une fois et une seule pour subir les étapes du procédé. L'homme du métier connait des méthodes de balayage automatiques d'une image tridimensionnelle et d'une image bidimensionnelle. L'essentiel est que le balayage soit systématique, et effectué selon toutes les directions prédéterminées. Il a été trouvé par l'expérimentation que l'ordre dans lequel les directions prédéterminées sont balayées est indifférent à l'obtention du bon résultat du présent procédé. Ces méthodes de balayage ne sont pas illustrées ici pour des raisons de complexité du dessin. 3.b) Définition du pas d'échantillonnage sur les droites de balayage On définit le pas S d'échantillonnage des voxels et pixels sur les droites de balayages dans une opération 631 en fonction de l'orientation de ces droites respectivement dans le repère CFC et dans le plan de l'image. En 3-D: Lorsque la droite de balayage est parallèle à un axe principal X, Y, Z du repère CFC, c'est-à-dire parallèle à un axe du repère CFC qui joint le centre à un sommet du cube, telle que 111' 1' et.'13, alors le pas S = 1. Lorsque la droite de balayage, telle que i42, est parallèle à un axe du repère CFC qui joint le centre au milieu d'une arête du cube, donc bissectrice d'un angle de 90 fait par deux axes principaux du repère CFC, par exemple X,Y, alors le pas S = /2. Lorsque la droite de balayage, telle que '43, est parallèle à un axe du repère CFC qui joint le centre à un sommet du cube, donc fait avec un plan de référence du repère CFC, tel que le plan (X, Y), un angle Arctg(1/42), alors le pas S = /3. Ainsi, le pas S a toujours des valeurs simples: S = 1, ou S = 12, ou S = f3. En référence à la FIG.5A, on note qu'une direction discrétisée de gradient qui se trouve selon un axe principal X, Y ou Z admet un plan perpendiculaire joignant les milieux des arêtes parallèles à cette direction discrétisée de gradient. Ainsi, les directions discrétisées de gradient O5 '5 admettent un plan ' 5adettunpn perpendiculaire 10, 12, 10', 12'. Il en résulte que quatre droites de balayage perpendiculaires à & ' sont possible: tI '1 ' 5 10 p:, I1 il, Ii' 12112 et ' 13 113. Deux de ces droites de balayage IV'll ll et I '13'q, sont parallèles aux grands axes et ont pour pas S = 1, et deux de ces droites de balayage %l'I 0 et 1l 12 sont parallèles à des bissectrices d'axes principaux et ont pour pas S = 1. En référence à la FIG.5B, on note que les directions discrétisées de gradient qui joignent le centre O à un milieu d'arête admettent un plan perpendiculaire joignant deux sommets et deux sommets diamétralement opposés. Ainsi, les directions discrétisées de gradient 06, 0'6 admettent un plan perpendiculaire 7, 1, 7', 1'. Il en résulte que quatre droites de balayage perpendiculaires à 06 O'6 sont possible: 1'4 Y. 'I 1V11, '1f1 et P'7qIF. Une de ces droites de balayage 'lI" l est parallèle aux grands axes et a pour pas S = 1, une de ces droites de balayage Y' '4'4 est parallèle à une bissectrice d'axes principaux et a pour pas S = f2, et deux de ces droites de balayage WP'1'[1 et IV' I". font un angle de 60 avec un axe principal et ont pour pas S = f3. En référence à la FIG.5C, on note que les directions discrétisées de gradient qui joignent le centre 0 à un sommet admettent seulement une droite perpendiculaire joignant deux milieux d'arêtes opposées. Ainsi, les directions discrétisées de gradient telle que 01 0' admettent une seule droite de balayage normale 1' 1 possible: IJ'12lIl2. Cette droite de balayage W'12Wp12 est une bissectrice de grands axes et a pour pas S = f2. En 2-D: Lorsque la droite de balayage est parallèle à un axe principal X, Y du repère du plan de l'image, c'est-à-dire telle que i'1, et X15' alors le pas S = 1. Lorsque la droite de balayage, telle que '13 et 1, est parallèle à une bissectrice d'axes principaux, tels que X,Y, donc bissectrice d'un angle de 90 , alors le pas S = /2. Lorsque la droite de balayage fait avec l'axe de référence X un angle de /8 telle que 1'2 et 414, alors on réalise une approximation en remplaçant r/8 par Arctg(1/2), pour que cette droite de balayage passe par un maximum de pixels. Dans ce cas, le pas S = /5. Ainsi, le pas S a toujours des valeurs simples: S = 1, ou S = f2, ou S = f5. De plus, chaque direction discrétisée de gradient n'admet qu'une seule droite de balayage perpendiculaire. Toutes les autres droites de balayage du système sont autres que perpendiculaires. En 3-D et en 2-D: Un avantage apporté par ce mode de définition du système de balayage est qu'il est simple à construire. Sa complexité n'excède pas celle de la discrétisation des gradients. Un autre avantage de ce système de balayage est que toute direction discrétisée de gradient admet au moins une droite de balayage qui lui est perpendiculaire. Un autre avantage est que l'on peut facilement déterminer les droites de balayage pour qu'elles passent non pas par des points quelconques difficiles à calculer mais par un maximum de voxels ou pixels de manière à ce que le pas de balayage noté S sur chaque droite de balayage ait des valeurs simples à utiliser. 3.c) Opération de balayage de l'image En 3-D et en 2-D: Un avantage du système de balayage proposé est qu'il permet une mise en oeuvre simple du présent procédé de filtrage récursif appliqué sur l'image tridimensionnelle ou bidimensionnelle. Une fois ce système de balayage défini, l'image est parcourue, dans une opération 700, de manière systématique et automatique, dans les deux directions portées par une même droite de balayage: la première direction par exemple Iij appelé sens causal, et la deuxième direction qui est alors i,'j appelé sens anticausal. A cet effet, les voxels P(x,y,z) ou pixels P(x,y) sont numérotés par rapport à un voxel ou pixel origine sur la droite de balayage. La droite est balayée du voxel ou pixel origine u0 en passant par les voxels ou pixels u1, u2,.. uj jusqu'au dernier voxel ou pixel un dans le sens causal, et réciproquement du voxel ou pixel numéroté un jusqu'à u0 dans le sens anticausal. En 2-D: A titre d'exemple, les FIG.7A et 7B illustrent respectivement le balayage causal et anticausal d'une image bidimensionnelle selon la droite de balayage,45'1'5 dans le plan (X,Y) de l'image à traiter. Ces balayages sont effectués de manière standard et automatique comme connu de l'homme du métier. 3.d) Définition de la fonction de filtrage récursif En 3-D et en 2-D: Le débruitage est obtenu par un lissage 800 de l'image durant le balayage 700. A cet effet, l'intensité I(x,y,z) ou I(x,y) de chaque voxel ou pixel courant P(x,y,z) ou P(x,y) respectivement est soumise à une fonction récursive de filtrage spatial calculée en chaque voxel ou pixel, pour fournir une nouvelle intensité lissée Q(x,y,z) ou Q(x,y) en chaque voxel ou pixel. Cette fonction effectue un lissage de l'intensité d'un plus grand degré dans la direction perpendiculaire à la direction discrétisée du gradient d'intensité relatif à ce voxel ou pixel que dans la direction discrétisée du gradient d'intensité. Autrement dit, en un voxel ou pixel donné, la fonction de filtrage produit une réponse impulsionnelle plus étalée dans la direction perpendiculaire au vecteur gradient que dans la direction du vecteur gradient. Pour lisser et débruiter l'image, on applique en chaque pixel courant uj, lors du passage causal dans une direction de balayage donnée 41j, une fonction récursive de traitement spatial ayant la formulation suivante: QC(uj) = I(uj) + a(uj) [QC(uj_1) I(u1)] (1) L'exposant "C" indique que l'intensité Q est obtenue en phase de balayage causal. Q(uj) est l'intensité nouvelle calculée au pixel uj de la droite de balayage dans le sens causal; Q(uj_1) est l'intensité nouvelle calculée au pixel uj_1 déjà balayé sur la droite dans le senscausal; o I(uj) est l'intensité d'origine en ce pixel courant u. et o a(uj) est un facteur de filtrage compris entre O et 1. Ensuite on applique au même pixel courant uj lors du balayage anticausal sur la même droite de balayage, mais dans la direction de balayage opposée de ll,'j, une fonction récursive de traitement spatial ayant la formulation suivante, en gardant la numérotation des pixels de u0 à un déterminée avant l'étape de balayage proprement dite: QA( (uj)QC(uj) + a(uj) [QA(uj+1) - QC(u.)] (2) L'exposant "A" indique que la nouvelle intensité Q est obtenue en phase anticausale. On rappelle que le pixel uo0 reste fixe sur la droite qui porte 4j. et 1'j quelle que soit la direction de balayage. Ainsi le passage causal, dans une première direction sur une droite de balayage d'orientation donnée, transforme l'intensité d'origine I(uj) d'un pixel courant uj en une nouvelle intensité QC(Uj) qui prend en compte l'intensité du pixel adjacent uj_1 qui vient juste d'être balayé et de se voir attribuer une nouvelle intensité causale QC(uj_1). Puis le passage anticausal sur la même droite de balayage, dans la direction de balayage opposée, transforme l'intensité QC(uj) nouvellement obtenue au passage causal A pour ce pixel numéroté uj en une nouvelle intensité QA(uj) qui prend en compte l'intensité du pixel adjacent numéroté uj+1 qui vient juste d'être balayé et de se voir attribuer une nouvelle intensité anticausale QA(uj+1). Le résultat du filtrage est stocké en 900. On analyse l'opération de cette fonction de lissage selon les positions relatives de la droite de balayage et de la direction discrétisée du gradient au point balayé. Dans une première situation. o la droite de balayage n'est pas dans un plan normal à la direction discrétisée du zradient, on applique systématiquement un filtrage appelé isotrope, c'est-à-dire que, conceptuellement, on définit un noyau de filtre isotrope, de même dimension dans toutes les directions. Dans une seconde situation, o la droite de balayaze se trouve dans un plan normal à la direction discrétisée du gradient, on applique svstématiquement un filtrage appelé anisotrope, c'est-à- dire que, conceptuellement, on définit un noyau de filtre anisotrope, ce noyau présentant une plus faible dimension dans la direction discrétisée du gradient que dans le plan normal à la direction discrétisée du gradient o se trouve la droite de balayage. Dans ce cas- là on se trouve en présence d'un gradient fort si le voxel ou pixel balayé est situé sur un BORD. Ce concept revient à effectuer en un voxel ou pixel balayé donné des passages sur toutes les droites de balayage avec un filtre à noyau dont la forme et les dimensions varient selon les droites de balayage et à combiner les résultats. Dans ces passages: Ce noyau est isotrope pour toutes les droites de balayage qui ne sont pas normales à la direction discrétisée du gradient au voxel ou pixel balayé. Ce noyau est anisotrope pour une droite de balayage qui est normale à la direction discrétisée du gradient au voxel ou pixel balayé. La combinaison de ces passages équivaut au passage d'un filtre ayant un noyau en forme de sphère applatie illustrée par la FIG.8B dont le rayon efficace dans le plan perpendiculaire à la direction discrétisée du gradient est appelé rayon équatorial noté oT, et dont le rayon dans la direction discrétisée du gradient est appelé rayon polaire noté oG. Le paramètre qui caractérise les noyaux isotropes dans la première situation o les directions de balayage ne sont pas normales à la direction discrétisée du gradient est appelé rayon des noyaux isotropes et noté a00. Un noyau isotrope de forme sphérique est illustré par la FIG.8A. Lorsque le gradient d'intensité se trouve être substantiellement faible, pratiquement nul, les dimensions ou paramètres du filtre anisotrope diminuent et deviennent identiques dans toutes les directions à un paramètre basique noté aO. 3.e) Détermination des paramètres des noyaux des filtres Dans cette opération, on détermine à partir du paramètre basique a0O les valeurs des paramètres des noyaux des filtres, notée aT pour le rayon efficace du noyau anisotrope dans le plan normal à la direction discrétisée du gradient et notée aGG pour les noyaux isotropes sur les droites de balayage non normales à la direction discrétisée du gradient. En référence à la FIG.3A, on définit par une opération 401 une valeur appelée gradient de référence IIGREFII qui est la valeur de gradient au-dessous de laquelle on déclare A PRIORI que le pixel correspondant n'appartient pas à un bord, c'est-à-dire au- dessous de laquelle le gradient est considéré suffisamment faible pour que l'on utilise le filtrage isotrope dont le noyau a pour rayon le paramètre basique a,. Au-dessus de cette valeur, on déclare que le pixel correspondant appartient à un bord et on applique le filtrage à noyau anisotrope. Par bord, on entend une zone de l'image qui présente une transition d'intensité, telle que la transition entre un objet et le fond. Pour déterminer IIGREFII, on mesure parmi les gradients par exemple stockés dans la mémoire 610, la valeur trouvée la plus élevée notée IIGMAXI, et on fixe la valeur de IIGREFII petite devant celle de IIGAxlI. A titre d'exemple: IIGREFII = IIG AXI/10 ou IIcREFII = IIGMII/20 En référence aux FIGS.8A, 8B et 3A, on définit alors par une opération 404 le paramètre basique c0o et par les opérations 402, 403 et 406 des relations entre le paramètre basique ao et les rayons aT et aG du noyau anisotrope: 2 O0 Oc= -. -(3) + IIGlIl)2 IIdREF 112 2 0 GT = (4) 1 + lIcGIJ IGREF 112 Le degré d'anisotropie AoT est lié aux rayons polaire et équatorial 2 _2 des noyaux anisotropes par la relation: Ao = T G (5) On définit aussi, par l'opération 405, à partir du rayon polaire et par son carré, le paramètre aGG qui caractérise le noyau isotrope de filtrage quand la droite de balayage n'est pas normale à la direction discrétisée du gradient: O2 G (6a) n2 o n1 est le nombre de dimensions de l'image traitée, et o n2 est le nombre des droites de balayage. En 3-i. tGG - (6b) l'image 3-D ayant 3 dimensions et 13 droites de balayage '11 à En 2-D: 2 2 G (6c) 1u 8 l'image ayant 2 dimensions et 8 droites de balayage l'1 à 8' 3.f) Détermination du facteur de lissage a En 3-D et 2-D: En référence à la FIG.3A, le facteur a(uj) appelé ci-après FACTEUR DE LISSAGE a est défini par l'opération 407 pour fournir un lissage d'un degré plus grand dans la direction perpendiculaire à la direction discrétisée du vecteur gradient d'intensité que dans la direction discrétisée de ce vecteur gradient d'intensité. Le FACTEUR DE LISSAGE a est choisi d'une part comme une fonction de la direction discrétisée Op du vecteur gradient GP et colnme une fonction du pas noté S des voxels ou pixels sur la droite de balayage 14j,,'j. Le FACTEUR DE LISSAGE a est choisi d'autre part comme une fonction du module du gradient d'intensité IIGpll au voxel ou pixel courant P(xyz) ou P(x,y). En 3-D et en 2-D: Les caractéristiques d'intensité I(x,y,z) ou I(x, y) et des gradients en chaque voxel P(x,y,z) ou pixel P(x,y) sont d'abord stockées par exemple dans une mémoire 610 et associées aux coordonnées de manière que: En 3-D: à x,y,z, correspondent I(x,y,z), O parmi 01 à 013, ou 1 à e'13 et Idpll, pour l'image tridimensionnelle. En 2-D: Pour l'image bidimensionnelle, à x,y correspondent I(x,y), Q parmi Q1 à 08' ou O' 1 à e 8 et Il0p' En 3-D et 2-D: Finalement l'opération 407 donne le facteur o par la relation (7) que l'homme du métier peut déduire d'un calcul des moments relatifs à une réponse impulsionnelle infinie, en se basant sur l'hypothèse que le gradient varie lentement. s2 -(±-2)2- ( R2 S_q(1 + SR)2 _ (7) comprise entre 0 et 1, o S est le pas d'échantillonnage des voxels ou pixels le long des droites de balayage. R est un facteur calculé en fonction des rayons polaire et équatorial du filtre anisotrope aG et aT, en fonction de la direction p du vecteur gradient Gp vis-à- vis de la droite de balayage,ji et en outre en fonction du pas d'échantillonnage S défini précédemment. En référence aux FIGs.3B et 3C, la détermination du facteur de lissage a revient donc à déterminer le facteur R qui définit le noyau anisotrope utilisé à chaque passage, lors du balayage de l'image. A cet effet, on détermine le rapport R2/S2. Cette détermination est faite dans les deux situations du balayage de l'image rencontrées. En 3-D: (FIG.3B) Première situation, en un voxel donné, P(x,y,z), la droite de balayage 14j n'est Pas normale à la direction discrétisée Op du vecteur gradient Gp. Dans la formule de a (7), on 2 2= 2 s2 exprime: R2/S2 GG/S2 Deuxième situation, en un voxel donné, P(x,y, z), la droite de balayage 4,j est normale à la direction discrétisée du vecteur gradient Gp dont l'orientation p discrétisée est encore à prendre en compte selon trois cas: Dans le premier cas o Op est parallèle à un axe principal X,Y,Z du repère CFC (et Gp perpendiculaire à 4j), on exprime R2/S2 = (02 + AOT)/S2 Dans le second cas o Op est parallèle à une bissectrice de deux axes principaux du repère CFC, donc parallèle à un axe qui joint le centre O au milieu d'une arête en faisant un angle de 45 avec un axe principal, (avec Gp perpendiculaire à 4,;), et en outre: si S = 1, on exprime: R2/S2 ( 2 T ou bien si S = 2, on exprime: R2/S2 = (a2 + 3 2 2 ou bien si S = 13, on exprime: R2/S = (oG + 2 /3 Les coefficients appliqués au paramètre AaT sont prévus pour que sa contribution soit affectée finalement d'un coefficient 1 après la réalisation des 13 combinaisons. Dans le troisième cas o {% est parallèle à un axe du repère CFC joignant le centre à un sommet du cube, cas qui se produit seulement si S = 2, on exprime: R2/S2 = (2 + 2) /2 En 2-D: (FIG.3C) Dans le traitement d'une image bidimensionnelle le noyau de filtrage est simplement un cercle aplati ayant un grand rayon aT dans la direction perpendiculaire à la direction discrétisée du vecteur gradient et un petit rayon aG dans la direction discrétisée du vecteur gradient. Une seule droite de balayage est perpendiculaire à la direction discrétisée du gradient. Première situation, en un pixel donné, P(x,y), la droite de balayage -i'j,'j n'est pas normale à la direction discrétisée Op du vecteur gradient Gp. Dans la formule de a (7), on exprime: R2/S2 =2 /S2 Deuxième situation, en un pixel donné P(x,y) la droite de balayage 14'j 4,j est normale à la direction discrétisée du vecteur gradient Gp. Le pas S sur la droite de balayage est encore à prendre en compte. Dans le premier cas o la direction discrétisée du gradient O est perpendiculaire à une droite de balayage 41'j '1j parallèle à un axe principal X,Y de l'image, o S = 1, on exprime: R2/S2 = ( sG +T)/S2 Dans le second cas o la direction discrétisée du gradient Op est perpendiculaire à une droite de balayage 1' j 4 parallèle à une bissectrice de deux axes principaux de l'image, en faisant un angle de 45 avec un axe principal, cas o: S = 12, on exprime: R2/S2 = (02 + AT)/2 Dans le troisième cas o la direction discrétisée du gradient O est perpendiculaire à une droite de balayage ayant un pas S = r5, on exprime: R2/S = (o2 + T)/5 Les valeurs de a sont alors calculées avec la relation (7) et tabulées en 620 pour le calcul des fonctions (1) et (2). APPAREIL A titre d'exemple, en référence à la FIG.9, un appareil d'imagerie médicale comprend un système de radiographie numérisée ayant des moyens pour mettre en oeuvre la méthode de réduction du bruit décrite plus haut. Cet appareil comprend une source 1 de rayons X, une table 2 pour recevoir le patient, un dispositif récepteur des rayons X ayant traversé le patient, et en particulier un dispositif intensification d'image 3 couplé à un tube vidéo 4 qui fournit des données à un système 5 de traitement d'image comprenant un microprocesseur. Ce dernier comprend plusieurs sorties dont une sortie 6 couplée à un moniteur 7 pour la visualisation des images traitées ou à traiter, ou bien de séquences d'images. L'image radiographique bidimensionnelle numérisée peut contenir 5122 ou 10242 pixels codés sur 8 bits ou 10 bits. Chaque pixel peut être ainsi affecté d'un parmi 256 ou 1024 niveaux d'intensité. Par exemple, les régions sombres ont un niveau d'intensité faible, et les régions claires de l'image ont un niveau d'intensité élevé. L'image radiographique tridimensionnelle numérisée peut contenir 1283 voxels. L'image numérisée peut être obtenue en mode fluoroscopique. L'invention peut trouver son application particulièrement dans le traitement d'images artériographiques. Le traitement de l'image numérique selon le procédé décrit plus haut est mis en oeuvre dans ses différentes étapes et opérations dans le système 5. Les données peuvent être stockées dans une zone mémoire non représentée. Des moyens d'enregistrement non représentés peuvent être en outre utilisés. REVENDICATIONS: 1. Procédé de traitement d'image pour conserver les bords d'objets représentés dans une image d'intensité et pour réduire le bruit, comprenant des étapes de: acquisition de l'image d'intensité sous forme d'une matrice multidimensionnelle de points (2-D, 3-D) détermination du gradient d'intensité en chaque point, défini par sa direction (Op) et son module (IlGpll), passage en chaque point d'un filtre spatial récursif à noyau anisotrope provoquant un degré de lissage plus important perpendiculairement à la direction du gradient que dans la direction du gradient d'intensité. 2. Procédé selon la revendication 1, comprenant des sous-étapes de: discrétisation des directions (Op) des gradients d'intensité dans un premier système d'un nombre N d'axes centrés, balayage automatique de l'image selon un second système de même nombre N de droites centrées ('1" i,lI) respectivement parallèles J] aux directions du premier système d'axes, et en outre telles que chaque direction discrétisée de gradient (Op) du premier système admette un nombre inférieur à N de droites de balayage non perpendiculaires et un nombre complémentaire égal ou supérieur à 1 de droites de balayage perpendiculaires, définition d'un filtrage à noyau variable, incluant un premier ensemble de sous-filtres à noyaux isotropes identiques passés au point courant sur les droites de balayage du second système qui sont non perpendiculaires à la direction discrétisée du gradient, et un second ensemble de sous-filtres à noyaux anisotropes passés au point courant sur les droites de balayage du second système qui sont perpendiculaires à la direction discrétisée du gradient, et combinaison des N résultats des N sous-filtres pour fournir le résultat du filtre spatial récursif. 3. Procédé selon la revendication 2, comprenant pour la mise en oeuvre du filtrage spatial récursif l'évaluation en chaque point courant sur chaque droite de balayage d'une fonction de filtrage spatial récursif comprenant un premier passage dans un sens causal et un second passage dans un sens anticausal, avec une évaluation de l'intensité filtrée (Qj) au pixel courant par un facteur (a) appelé facteur de lissage (0 4. Procédé selon la revendication 3, o l'image est tridimensionnelle et o les systèmes de discrétisation des orientations des gradients (O) et des droites de balayage (Y'J.3j) p 33J comportent 13 axes centrés qui sont les axes d'un système cubique centré à faces et arêtes centrées (CFC). 5. Procédé selon la revendication 3, o l'image est bidimensionnelle et o les systèmes de discrétisation des orientations des gradients (O,) et des droites de balayage (qI'j /.j) comportent 8 axes centrés qui sont régulièrement répartis dans le plan de l'image. 6. Procédé selon l'une des revendications 3 à 5, comprenant des sous-étapes de: évaluation d'un paramètre basique (aO) de dimension d'un noyau de filtrage isotrope pour filtrer une zone de l'image à gradient substantiellement nul ( IGpll0O), évaluation des paramètres de dimensions (aG, aT) du noyau anisotrope des sous-filtres passés perpendiculairement une direction discrétisée du gradient, et du paramètre de dimension (aGG) du noyau isotrope des sous-filtres passés selon des droites de balayage non perpendiculaires à la direction discrétisée du gradient, et de la contribution à l'anisotropie (AuT) en fonction du paramètre basique (a0) et du module du gradient ( IIGpllO), évaluation d'un paramètre de dimension (R) qui définit le noyau anisotrope utilisé à chaque passage en fonction de la contribution à l'anisotropie (AoT) et du paramètre de dimension des noyaux isotropes (aGG), et évaluation du facteur de lissage (a) comme une fonction du paramètre de dimension (R) dudit noyau de filtre anisotrope utilisé à chaque passage. 7. Procédé selon la revendication 6, o les noyaux isotropes sont des sphères et o les noyaux anisotropes sont des sphères aplaties ayant leur plus grandes dimensions perpendiculaires à la direction discrétisée du gradient. 8. Procédé selon l'une des revendications 3 à 7, comprenant des sous-étapes de: évaluation d'un pas d'échantillonnage (S) sur toutes les droites de balayage, et évaluation du facteur de lissage en fonction desdits pas d'échantillonnage. 9. Procédé selon l'une des revendications 3 à 8, comprenant une évaluation de l'intensité filtrée en un pixel donné comme la combinaison des intensités filtrées selon toutes les droites de balayage. 10. Appareil d'imagerie médicale, ayant des moyens pour réduire le bruit dans une image médicale numérique bruitée, comprenant un système d'acquisition d'une image d'intensité sous forme d'une matrice multidimensionnelle de points, et un système de traitement d'image qui a un accès aux données de l'image et à un système d'affichage d'images et qui comprend un microprocesseur pour mettre en oeuvre un procédé selon l'une des revendications 1 à 9.