REDUCTION DE BRUIT DANS UNE IMAGE NUMERIQUE PAR TRANSFORMEE EN COSINUS DISCRETE

DOMAINE TECHNIQUE DE L' INVENTION

La présente invention concerne les procédés et dispositifs de filtrage d'images numérisées, en particulier les procédés et dispositifs pour éliminer le bruit dans une image numérisée.

Certaines images sont très bruitées, ce qui en réduit la lisibilité. Tel est le cas, par exemple, des images obtenues à partir de signaux à bas niveau, dans lesquels les signaux utiles sont fortement perturbés par des parasites et autres phénomènes physiques qui altèrent la transmission des signaux.

Tel est le cas, également, de certaines images d'imagerie médicale, telles que les images scintigraphiques .

La technique de scintigraphie en médecine nucléaire consiste à injecter à des patients des molécules biologiques marquées par des isotopes radioactifs émetteurs gamma. Les photons gamma sont détectés à l'aide d'une gamma caméra.

Les photons gamma sont émis de façon aléatoire par les particules émettrices injectées dans le corps humain. Les signaux captés dépendent d'une part de la quantité de particules émissives dans une zone donnée, et d'autre part des aléas d'émission de chaque particule. Dans une image scintigraphique, qui est l'image de la distribution des photons gamma provenant d'un organisme, chaque élément de l'image ou pixel est un nombre entier qui est le nombre de particules, c'est-à-dire le nombre de photons détectés en face de lui. On constate que ces valeurs entières sont distribuées autour d'une valeur moyenne, selon la loi de Poisson.

La distribution de la loi de Poisson suit la formule μ _-n Pr (X=k) k !

Pr (X=k) est la probabilité pour que le pixel X prenne la valeur k. μ est la valeur moyenne de la distribution.

La figure 1 donne un exemple d'une distribution selon la loi de Poisson pour une moyenne μ = 5.

Un autre cas de distribution selon la loi de Poisson est la distribution des photons X tels que ceux utilisés en

tomodensitométrie (scanner X) qui consiste à visualiser les structures anatomiques d'un organisme en utilisant les différences d'atténuation des rayons X par les tissus biologiques.

Dans tous les cas, les images sont entachées d'un bruit poissonien, résultant du caractère aléatoire des émissions, et ce bruit est d'importance relative d'autant plus grande que le nombre d'événements détectés est petit. En effet, l'erreur relative de la mesure est donnée par la formule : Erreur relative = [ (variance) (X) ] ^1/2 / moyenne (X) = Vμ / μ = 1/Vμ

L'erreur relative varie donc en sens inverse du nombre d'événements détectés.

Pour améliorer la lisibilité d'une image, on est donc naturellement conduit à réduire le bruit, et pour cela à augmenter le nombre d'événements détectés.

Une première façon est d'acquérir les images pendant un temps plus long. Mais cela peut être difficile lorsque le sujet bouge, par exemple en radiologie à cause des mouvements respiratoires et cardiaques. De plus, l'augmentation des temps d'acquisition immobilise les appareils de détection plus longtemps, ce qui pose des problèmes importants de planning en imagerie médicale .

Une seconde façon d'augmenter le nombre de photons détectés est d'augmenter le flux des particules. En imagerie médicale, cela implique l'augmentation de la dose radioactive injectée ou du flux de photons X. Mais cela va augmenter d'autant la dose de rayonnement reçue par le patient, ce qui est en contradiction avec les principes de radioprotection.

Enfin, on peut imaginer utiliser des appareils de détection plus efficaces, mais cela pose immédiatement les problèmes de coût.

Pour illustrer l'influence de l'augmentation du nombre de photons détectés sur la qualité de l'image, on a reproduit sur la figure 2 une image d'un fantôme numérique à lignes (en haut à gauche) . A partir de celle-ci, trois acquisitions scintigraphiques ont été simulées en ajoutant un bruit de Poisson. L'image en haut à droite correspond à un nombre de coups moyen de 30 (bruit élevé) . L'image en bas à gauche correspond à un nombre de coups moyen de 60

(bruit moyen). L'image en bas à droite correspond à un nombre de coups moyens de 120 (bruit faible) . La qualité de cette dernière image est très supérieure.

On - comprend que le bruit poissonien consiste en une variation aléatoire entre un pixel et les pixels adjacents, variation qui n'est pas liée à la différence réelle de concentration de particules émissives entre la zone face au pixel considéré et les zones face aux pixels adjacents. Il s'agit ainsi d'un bruit à haute fréquence, c'est-à-dire introduisant des variations brusques d'un pixel à l'autre dans la suite des pixels d'une image.

Comme il est souvent très difficile d'augmenter le nombre de photons détectés, l'invention se propose de réduire le bruit à haute fréquence contenu dans une image par l'utilisation d'un dispositif de filtrage particulier.

Il a déjà été proposé de filtrer les images, dans le but de réduire l'influence du bruit. De très nombreux filtres ont été proposés .

Le filtre le plus simple et le plus connu consiste à remplacer la valeur de chaque pixel par la moyenne des valeurs des pixels adjacents. On a également proposé le filtre gaussien, consistant à donner un poids de distribution spatiale gaussienne aux différents pixels voisins, et en remplaçant le pixel central par une combinaison linéaire des pixels voisins.

Un exemple est décrit dans le document WO 01/72032, qui définit des coefficients particuliers permettant à la fois le filtrage, et conservant si possible un certain contraste de 1 ' image .

Ce document exprime d'ailleurs le problème de tous les filtres pondérés ou non : tous ces filtres améliorent le rapport signal/bruit, mais au prix d'une forte dégradation de la résolution spatiale et du contraste de l'image. En effet, on conçoit aisément qu'un filtre effectuant une moyenne pondérée entre plusieurs pixels voisins réduit nécessairement les variations brusques d'intensité entre les pixels voisins, et réduit donc les effets de contraste.

A titre d'exemple, on a illustré sur la figure 3 le traitement d'une image de fantôme numérique bruitée (bruit moyen de

la figure 2 en bas à gauche) par deux filtres : la vue en haut à gauche est l'image brute non filtrée ; la vue en haut à droite est l'image filtrée par un filtre 3 x 3 pondéré ; la vue en bas à gauche est l'image filtrée par un filtre 3 x 3 médian, qui remplace donc chaque pixel par la médiane des pixels adjacents. On voit, sur les images filtrées, que les contours ne sont pas très nets, par suite d'une réduction de la résolution spatiale et du contraste.

On a également proposé de filtrer des images par des méthodes d'analyse statistique multivariée appliquées à un tableau de nombres exprimant chacun le degré de luminosité d'un pixel correspondant dans une image numérisée.

Ainsi, le document XP 002212651 AEBERSOLD, STADELMANN, ROUVIERE décrit un traitement des images de microscopes électroniques. Ce traitement est effectué au moyen d'une analyse factorielle des correspondances, appliquée sur un tableau constitué par les pixels de l'image numérisée. L'analyse est appliquée sur l'image numérisée tout entière, en recherchant sur cette image entière les facteurs de correspondance les plus représentatifs de l'image, c'est-à-dire ceux qui correspondent à des valeurs propres les plus grandes, puis en reconstruisant l'image à partir de ces seuls facteurs représentatifs. Si en théorie la méthode permet d'obtenir un filtrage de bonne qualité d'une image donnée, en pratique elle ne permet pas d'adapter le filtrage à la structure de l'image à filtrer, de sorte que le résultat reste décevant.

Le document XP 001074312 HANNEQUIN, LIEHN, VALEYRE, propose d'appliquer l'analyse factorielle des correspondances à une série d'images entières de scintigraphie, et suggère d'utiliser le test du rapport de vraisemblance pour déterminer les facteurs de correspondance à utiliser pour reconstruire la série d'images. La méthode est appliquée à nouveau à la totalité des images . Le résultat du filtrage est un peu amélioré, mais reste décevant pour certaines structures d'image à filtrer.

Les mêmes difficultés se retrouvent dans les documents XP 008007773 et XP 008007666, qui appliquent des traitements statistiques sur la totalité d'une image.

L'inventeur a proposé plus récemment, dans le document WO 03/050760, de filtrer une image par l'analyse factorielle des

correspondances (AFC) appliquée de manière séparée à des portions adjacentes de l'image. Cela a permis d'améliorer sensiblement la résolution de l'image filtrée (figure 3, image en bas à droite), mais au prix d'une augmentation du volume et de la durée des calculs .

On constate encore une dégradation de l'image dans les parties d'image à haute résolution, c'est-à-dire dans les parties inférieure droite et supérieure gauche de chaque image en comparant les images non bruitées (en haut à gauche de la figure 2) et filtrées AFC (en bas à droite de la figure 3) .

Il y a encore un besoin pour un dispositif de filtrage plus efficace et rapide, permettant d'éliminer le bruit de haute fréquence sans réduire sensiblement la résolution ou le contraste de l ' image .

EXPOSE DE L'INVENTION

La présente invention a ainsi pour but de permettre un filtrage à la fois plus rapide et plus efficace d'images à bruit poissonnien telles que des images d'imagerie médicale, de scintigraphie, de tomodensitometrie, produisant rapidement des images de meilleure qualité, à l'aide de moyens particulièrement peu onéreux.

En outre, selon l'invention, ce filtrage est de préférence adaptatif, permettant alors de réduire au mieux la perte d'informations résultant du filtrage.

Pour atteindre ces buts ainsi que d'autres, l'invention prévoit un procédé de traitement d'une image numérisée constituée d'un tableau T de nombres x _U3) exprimant chacun le degré de luminosité d'un pixel de rang i, j correspondant (par exemple le nombre de particules détectées dans le pixel correspondant) , le procédé comprenant la réduction du bruit à haute fréquence (par exemple un bruit statistique poissonnien) par les étapes suivantes : a) décomposer le tableau T en une suite continue de p tableaux élémentaires de même dimension ayant chacun n pixels, b) ordonner les données de la suite des tableaux élémentaires en un tableau de traitement X de p lignes et n

colonnes, chaque ligne i étant formée de la suite ordonnée des pixels du tableau élémentaire de rang i, d) effectuer sur le tableau de traitement X une transformation orthogonale dans l'espace des fréquences, considérant que les n colonnes sont les variables, pour en extraire les n fonctions orthogonales représentatives, associées à leur coefficient, f) générer un tableau de traitement reconstitué XR de nombres xr _(1/ J ₁ en reconstruisant indépendamment chaque ligne i en prenant en compte seulement les fonctions ayant un poids significatif avec la ligne i, et en rétablissant les degrés de luminosité absolus, puis générer un tableau reconstitué TR constituant l'image numérisée reconstituée dans laquelle le bruit à haute fréquence a ainsi été réduit.

Selon l'invention, lors de l'étape d) , on utilise une transformation orthogonale préétablie à n coefficients orthogonaux préétablis. La transformation orthogonale prend un ensemble de pixels de l'image et les transforme en une représentation équivalente dans l'espace fréquentiel.

Le traitement du tableau de traitement X par transformation orthogonale préétablie utilise des tables de calcul préenregistrées, qui sont identiques quelles que soient les images à filtrer. Il n'est ainsi pas nécessaire de les recalculer à chaque opération de filtrage, de sorte que les calculs sont plus rapides et moins volumineux que dans une méthode d'analyse statistique multivariée telle que l'analyse factorielle des correspondances.

Grâce à la reconstruction ligne par ligne, l'invention adapte le filtrage à chaque sous-ensemble de l'image, ce qui améliore considérablement le résultat obtenu.

Selon un mode de réalisation avantageux, le procédé utilise, à l'étape d) , une transformée en cosinus discrète (DCT), selon laquelle on calcule, pour chaque ligne du tableau de traitement X, les coefficients correspondant aux n fonctions orthogonales .

De préférence, lors d'une étape e) , on calcule les cosinus carrés des lignes sur les n fonctions orthogonales, et on utilise

les cosinus carrés comme test du poids des fonctions représentatives dans la ligne i.

On constate que la transformée en cosinus discrète améliore encore la qualité du filtrage et notamment la résolution des images filtrées obtenues. Le résultat est illustré sur la figure 4, à l'issue d'un filtrage appliqué aux images bruitées de la figure 2.

Selon un mode de réalisation avantageux, l'invention propose un procédé permettant l'auto-adaptation de la reconstruction, afin d'éliminer au mieux le bruit de haute fréquence sans affecter la qualité de l'image. Le moyen utilisé est de ne conserver que la variance du signal en éliminant la variance du bruit. Pour cela, selon l'invention, on arrête la reconstruction d'une ligne du tableau reconstitué XR dès qu'on a utilisé un nombre de fonctions suffisant pour que la variance de la ligne reconstruite soit supérieure à la variance du signal initial de laquelle on a soustrait la variance estimée du bruit.

Le test s'applique particulièrement bien par exemple au traitement d'une image à bruit poissonien. Dans ce cas, une propriété de la loi de Poisson est que sa moyenne est égale à sa variance : variance (X) = moyenne (X) = μ

II en résulte qu'une bonne estimation de la variance du bruit consiste dans le calcul de la moyenne des signaux de la ligne i .

Selon un autre aspect, l'invention propose un dispositif de traitement d'images numérisées, comprenant une mémoire, une unité de calcul, un dispositif d'entrée-sortie pour recevoir les données constituant l'image numérisée à traiter, des moyens de visualisation et/ou d'impression pour visualiser l'image traitée, et un programme enregistré en mémoire et adapté pour mettre en œuvre le procédé tel que défini ci-dessus.

L'invention s'applique notamment à une installation d'imagerie médicale comprenant un tel dispositif.

DESCRIPTION SOMMAIRE DES DESSINS

D'autres objets, caractéristiques et avantages de la présente invention ressortiront de la description suivante de modes

de réalisation particuliers, faite en relation avec les figures jointes, parmi lesquelles :

- la figure 1 illustre _^ une distribution de Poisson de moyenne 5 ; la figure 2 illustre l'amélioration d'une image lors de l'augmentation de la durée d'observation d'un phénomène poissonien ;

- la figure 3 illustre les résultats obtenus par des filtres pondérés, médians ou encore par analyse factorielle des correspondances de l'art antérieur ;

- la figure 4 illustre le résultat que l'on peut obtenir avec un filtre selon un mode de réalisation de l'invention ;

- la figure 5 illustre les principales étapes d'un procédé de filtrage selon la présente invention ;

- la figure 6 illustre un principe de décalage utilisé pour limiter les effets de bord selon le procédé de l'invention ; et

- la figure 7 illustre schématiquement une installation d'imagerie incorporant un dispositif de traitement d'images selon l'invention.

DESCRIPTION DES MODES DE REALISATION PREFERES On considérera tout d'abord la figure 7 où l'on a illustré schématiquement une installation d'imagerie médicale comprenant un capteur de rayons gamma 1 tel qu'une gamma caméra se déplaçant dans deux directions face au sujet 2 à observer, et captant les rayons gamma 3 issus des particules radio-émissives préalablement injectées dans le corps du patient 2. Le capteur de rayons gamma 1 envoie à une unité de calcul 4 la suite des signaux images des photons reçus sur chaque zone élémentaire ou pixel du capteur de rayons gamma 1, l'unité de calcul 4 mémorisant dans une mémoire 5 les nombres de photons de chaque pixel, correspondant à l'intensité du pixel. On retrouve ainsi dans la mémoire 5 une image numérisée constituée d'un tableau T de nombres x (i, j) exprimant chacun le nombre de photons détectés (ou degré de luminosité) d'un pixel de la ligne i et de la colonne j de la zone observée 1. L'installation comprend en outre, selon l'invention, un programme enregistré dans la mémoire 5 et permettant de piloter l'unité de calcul 4 pour filtrer l'image ainsi numérisée et pour produire sur un dispositif de visualisation ou d'impression 6 une image filtrée de bonne qualité et dans laquelle on a extrait le bruit de haute fréquence.

On va maintenant décrire un procédé de filtrage d'image selon un mode de réalisation de la présente invention, en relation avec la figure 5.

On retrouve le tableau T, constituant l'image numérisée. En pratique, les images numérisées contiennent un grand nombre de pixels. Pour la simplicité de l'exposé, on considère une image constituée par 8 x 8 pixels, de forme carrée, chaque pixel étant illustré par un petit carré.

La première opération a) du procédé de réduction de bruit selon l'invention consiste à décomposer le tableau T en une suite continue de p tableaux élémentaires de même dimension ayant chacun n pixels, n étant une puissance de 2. Dans l'exemple illustré sur la figure 5, on considère quatre tableaux élémentaires Tl, T2, T3 et T4, ayant chacun 16 pixels.

Ensuite, selon une étape b) , on ordonne les données de la suite des tableaux élémentaires T1-T4 en un tableau de traitement X de p lignes et n colonnes, chaque ligne i étant formée de la suite ordonnée des pixels du tableau élémentaire de rang i. On retrouve ainsi, sur la première ligne du tableau X, les pixels 1 à 16 du tableau Tl, rangés en ordre. De même, on retrouvera dans la seconde ligne du tableau X les pixels rangés en ordre du tableau T2, et ainsi de suite. Dans l'exemple, le tableau de traitement X a ainsi quatre lignes de 16 colonnes.

Dans l'exemple de la figure 5, on a décomposé le tableau T en quatre tableaux élémentaires carrés T1-T4 ayant chacun 4 lignes et 4 colonnes. On pourra toutefois, sans sortir du cadre de l'invention, décomposer le tableau T en une suite de tableaux carrés ou rectangulaires ayant tous la même dimension mais dont le nombre de lignes et/ou de colonnes diffère de 4.

Ensuite, selon une étape c) on peut éventuellement normaliser le tableau de traitement X pour obtenir une matrice normalisée Xn dont chaque élément xni _j de ligne i et de colonne j est pondéré par une transformée utilisant la moyenne des valeurs des éléments de la ligne i et la moyenne des éléments de la colonne j-

Ensuite, selon l'étape d) , on calcule, par la transformée en cosinus discrète (DCT) , pour chaque ligne, ses n coefficients correspondant aux n fonctions orthogonales .

Pour le calcul par transformée en cosinus discrète (DCT) , on pourra utiliser la formule :

N-I N-I (2x+l) u π (2y+l) v IT T(u,v) = σ σ f (x,y) α(u) α(v) cos [ ] cos [ ] x=0 y=0 2N 2N

N est la taille du tableau élémentaire f (x,y) est la valeur du pixel au point de coordonnées x et y-

T(u,v) est l'équivalent dans l'espace des fréquences. α(u) = V (1/N) pour u=0

= V (2/N) pour u=l, ,N -1 et de même pour α(v)

On calcule ensuite les cosinus carrés des p lignes sur les n fonctions orthogonales .

Au cours d'une étape e) , on classe ensuite les n fonctions en ordre décroissant en fonction de leur poids respectif.

Au cours d'une étape f) , on génère un tableau de traitement reconstitué XR de nombres xr (i, j) en utilisant seulement les q premières fonctions représentatives de chaque ligne.

Enfin, on génère un tableau reconstitué TR, constituant l'image numérisée reconstituée, dans laquelle on a réduit le bruit à haute fréquence, par exemple le bruit statistique poissonnien.

Selon l'étape e) , on peut calculer le cosinus carré par la formule : cos ² _k (i)= c _k (i) c _k (i) où C _k (i) est le coefficient de la ligne i pour la fonction k selon la formule : n c _k (i) = σ x±j • fk(j)

J= ¹ dans laquelle fk(j) est la j ^ieme valeur de la fonction fk.

Au cours de l'étape e) , on peut avantageusement classer les n facteurs en fonction de leur cosinus carré ainsi ^' calculé.

Selon l'invention, la reconstitution du tableau de traitement reconstitué XR s'effectue ligne par ligne, de façon indépendante, en prenant en compte seulement les q facteurs ayant le cosinus carré maximal pour la ligne i.

Supposant que l'on prenne en considération les q premières fonctions, la valeur reconstruite xri _j (q) de l'élément du tableau de traitement reconstitué XR de ligne i et de colonne j est calculée par la formule : q xri _j (q) = σ c _k (i)fk(j) k=l où fk(j) est la j ^ièrae composante de la k ^ième fonction orthogonale.

A partir du tableau de traitement reconstitué XR, sur la figure 5, la reconstruction du tableau reconstitué TR s'effectue ligne par ligne, la première ligne du tableau de traitement reconstitué XR constituant les pixels du premier tableau élémentaire TRl, et ainsi de suite.

Selon l'invention, on peut avantageusement automatiser l'adaptation du dispositif de filtrage au contenu de l'image. Cette automatisation se fait pour chaque zone de l'image correspondant chacune à l'un des tableaux élémentaires Tl à T4. Pour cela on reconstruit le tableau reconstitué XR ligne par ligne. On effectue le calcul des valeurs reconstruites xry d'une ligne i d'éléments du tableau reconstitué XR par un calcul pas à pas :

- on calcule successivement la valeur des éléments xr^ _j pour des valeurs q croissantes du nombre de fonctions pris en compte, on calcule à chaque fois la variance résiduelle var_res (q) de la ligne i,

- on compare la variance résiduelle à la variance estimée du bruit à réduire,

- et on arrête le calcul de la ligne i lorsque la variance résiduelle de la ligne i n'est plus statistiquement supérieure à la variance estimée du bruit de la ligne i dans l'image de départ, obtenant ainsi une image finale (Im finale) estimée sans bruit.

En pratique la variance résiduelle Var_res (q) de la ligne i est la variance de la différence entre la ligne i du tableau de traitement X et la ligne i du tableau de traitement reconstitué XR tel que reconstruit avec q fonctions.

Le test de comparaison de la variance résiduelle et de la variance estimée du bruit peut avantageusement être effectué par : a) calculer la variable t par la formule t = (Var__bruit)xhi (ddl) /ddl dans laquelle xhi (ddl) est la valeur donnée par la table du χ ² pour un risque de 5 % et un nombre ddl de degrés de liberté, ddl est le nombre de degrés de liberté, avec ddl = n-q-1 q étant le nombre de facteurs pris en compte, b) arrêter la reconstruction lorsque la variable résiduelle Var_res (q) est inférieure à t.

Dans le cas où l'on applique le procédé au traitement d'une image ayant un bruit qui suit une loi de Poisson, la variance estimée du bruit de la ligne i est prise égale à la moyenne des éléments x^ de la ligne i du tableau de traitement X.

Afin de diminuer encore l'influence du bruit sur les résultats, et de réduire les effets de bords, la procédure décrite ci-dessus peut être répétée plusieurs fois sur la même image en ^• décalant à chaque fois de un pixel le découpage en tableaux élémentaires . La figure 5 illustre la première procédure pour un décalage de 0 en x et de 0 en y. La figure 6 illustre la seconde procédure pour un décalage de 1 pixel en x et 1 pixel en y : le tableau élémentaire T ¹ I est décalé de 1 pixel vers la droite et de 1 pixel vers le bas dans le tableau T. Par exemple, pour un découpage en carrés 4 x 4, on réalisera 16 fois la procédure, avec des décalages en x allant de 0 à 3 et des décalages en y allant de 0 à 3. L'image finale (Im_finale) , estimée sans bruit, sera la moyenne des seize images ainsi reconstituées.

Il peut être tenu compte dans cette moyenne du nombre de fois où chaque pixel de l'image est réellement inclus dans le traitement, afin de ne pas faire apparaître d'effets de bord. Un

' autre avantage de la répétition est de s'affranchir des artefacts

géométriques qui peuvent apparaître en raison du découpage en rectangles élémentaires .

La figure 4 illustre le résultat d'un filtrage selon l'invention pour l'image de fantôme numérique non bruitée (en haut à gauche) et pour les trois images bruitées présentées sur la figure 2.

La présente invention n'est pas limitée aux modes de réalisation qui ont été explicitement décrits, mais elle en inclut les diverses variantes et généralisations contenues dans le domaine des revendications ci-après.