FR2978273A1

FR2978273A1 - Procede de reduction du bruit dans une sequence d'images fluoroscopiques par filtrage temporel et spatial

Info

Publication number: FR2978273A1
Application number: FR1102306A
Authority: FR
Inventors: Jeremie Pescatore; Carole Amiot
Original assignee: Thales SA
Current assignee: Thales SA
Priority date: 2011-07-22
Filing date: 2011-07-22
Publication date: 2013-01-25
Anticipated expiration: 2031-07-22
Also published as: WO2013014110A1; FR2978273B1; EP2734975B1; US9058654B2; US20140153842A1; JP6120843B2; JP2014520641A; EP2734975A1

Abstract

L'invention se situe dans le domaine de l'imagerie médicale par fluoroscopie. Elle concerne un procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X. Le procédé selon l'invention comprend les étapes successives suivantes pour chaque image (x ) de la séquence : ▪ appliquer un filtrage temporel (22) sur l'image (x ) acquise à un instant n, ledit filtrage temporel comportant un traitement de réduction du bruit fluoroscopique (FNR), ▪ appliquer un filtrage spatial (23) sur l'image (y ) acquise à l'instant n, ledit filtrage spatial comportant les sous-étapes suivantes : - transformer (231) l'image (y ) acquise à l'instant n du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image transformée étant représentée par un ensemble de coefficients, - seuiller (232) les coefficients de l'image (y ) par une fonction de seuillage, la fonction de seuillage annulant les coefficients inférieurs à un troisième seuil prédéterminé, et conservant ou ajustant les coefficients supérieurs au troisième seuil prédéterminé, - transformer (233) l'image (z ) dont les coefficients ont été seuillés du domaine des curvelets au domaine spatial par une transformation en curvelets inverse.

Description

PROCÉDÉ DE REDUCTION DU BRUIT DANS UNE SÉQUENCE D'IMAGES FLUOROSCOPIQUES PAR FILTRAGE TEMPOREL ET SPATIAL L'invention se situe dans le domaine de l'imagerie par rayons X et, plus spécifiquement, celui de l'imagerie médicale par fluoroscopie. Elle concerne un procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X.

Dans le domaine de l'imagerie médicale, l'imagerie par fluoroscopie à rayons X permet de fournir un flux d'images d'un patient à un médecin lors d'interventions chirurgicales dites minimalement invasives, par exemple la chimio-embolisation du foie, la vertébroplastie, la cathérisation des anévrismes ou le traitement des rétrécissements vasculaires. Les images peuvent notamment aider le guidage d'instruments chirurgicaux. Une intervention guidée par fluoroscopie est typiquement réalisée en faisant passer un cathéter à l'intérieur du réseau vasculaire du patient. Un agent de contraste peut ou non avoir été injecté préalablement dans le réseau vasculaire afin de l'opacifier et d'améliorer la visibilité des vaisseaux. Une telle intervention est généralement relativement longue et la dose de rayons X auxquels le patient est soumis doit être limitée afin d'éviter d'occasionner des lésions ou des brûlures des tissus. Du fait de cette limitation de la dose de rayons X utilisée, les images fluoroscopiques obtenues comportent un niveau de bruit relativement élevé, et donc un rapport contraste sur bruit relativement faible, rendant les images difficilement lisibles.

Des traitements de filtrage sont mis en oeuvre afin de diminuer le bruit quantique présent dans ces images et augmenter leur rapport contraste 25 sur bruit. Un traitement connu sous le nom de traitement de réduction du bruit fluoroscopique, ou traitement FNR selon le sigle de l'expression anglo- saxonne "Fluoroscopic Noise Reduction", est ainsi généralement réalisé sur les images acquises par le détecteur à rayons X. L'objectif de ce traitement 30 FNR est de filtrer le bruit présent dans l'image tout en préservant le contraste de l'information présente dans cette image. Le traitement FNR est réalisé en appliquant un filtre temporel aux régions des images où il n'y a pas de déplacement. L'existence ou l'absence de déplacements dans des régions de l'image est détectée à partir de la variation individuelle de l'intensité de chacun des pixels considérés séparément. Un pixel est considéré comme étant en déplacement lorsque sa variation d'intensité entre deux images dépasse un seuil lié à l'écart type du bruit. Ces pixels déterminés comme étant en déplacement ne sont pas, ou très peu, filtrés. En revanche, un traitement par filtrage temporel récursif est appliqué aux pixels déterminés comme étant fixes. Le traitement FNR présente l'avantage de ne pas colorer le bruit spatialement. Autrement dit, il ne fait pas apparaître d'objets fictifs résultant de groupements de pixels d'intensités voisines entre elles. Cependant, le traitement FNR présente une capacité de débruitage relativement limitée et tend à faire disparaître les objets d'intérêts dans l'image. A côté du traitement par filtrage temporel FNR, les traitements par filtrage spatial présente des propriétés intéressantes. Parmi les différentes variantes de traitement par filtrage spatial, le filtrage par transformation de l'image dans le domaine des ondelettes est l'une des techniques les plus utilisées. La théorie des ondelettes, introduite dans les années 1980, a été développée à partir de la transformée de Fourier. La transformée de Fourier est utilisée pour décrire de nombreux phénomènes physiques. Elle constitue en effet un outil mathématique performant pour décrire de nombreux phénomènes physiques, notamment les phénomènes physiques modélisés par des signaux stationnaires. En revanche, les représentations des phénomènes physiques modélisés par des signaux non stationnaires ne sont généralement pas satisfaisantes. En particulier, la transformée de Fourier est incapable de localiser les portions du signal où la fréquence change subitement. Pour y remédier, une transformée de Fourier à fenêtres a été développée. Elle consiste à multiplier le signal par une fenêtre de dimension donnée. L'analyse devient donc locale, la dimension de la fenêtre déterminant la résolution temporelle obtenue, ou la résolution spatiale dans le cas de signaux à deux dimensions. La théorie des ondelettes s'appuie sur la transformée de Fourier à fenêtres. Elle s'en différencie de par le fait qu'elle permet, au cours de la transformation, de changer la dimension de la fenêtre. Pour cette raison, on parle d'analyse multirésolution ou analyse multi-échelle.

Une autre différence est que le signal n'est plus décomposé en une somme de sinus et de cosinus, mais selon des fonctions appelées "ondelettes". Les ondelettes présentent l'avantage de transformer l'image en une représentation dite creuse. Par représentation creuse, on entend une représentation nécessitant un faible nombre de paramètres pour représenter fidèlement l'image. Cependant, les ondelettes sont mal adaptées à la représentation de discontinuités de type ligne, telles que les contours d'objets, car leur forme ne leur permet pas d'exploiter la forme des objets. En conséquence, l'information sur ces discontinuités est contenue dans les coefficients de nombreuses échelles. Pour obtenir une faible erreur de reconstruction, il faut donc calculer un grand nombre de coefficients, entraînant un coût de calcul important. Pour une application de fluoroscopie où les images doivent être traitées "en temps réel", c'est-à-dire dans des durées inférieures à la période de rafraichissement des images, ce coût de calcul implique des ressources matérielles importantes.

D'autres théories ont été développées pour réduire la redondance des informations. Il s'agit notamment des ondelettes isotropes et anisotropes, ou encore des arbres duaux orthogonaux ou non. Au début des années 2000, la théorie des curvelets est apparue. Cette théorie est également fondée sur le principe de l'analyse mufti- échelles. La géométrie des curvelets leur permet de tirer profit de la régularité des contours des objets. Une image peut ainsi être représentée par un nombre relativement faible de coefficients. Cependant, la transformation en curvelets présente l'inconvénient d'introduire des artéfacts dans l'image, notamment une coloration spatiale du bruit. L'introduction d'artéfacts est d'autant plus forte que le rapport contraste sur bruit est faible. Or, dans le domaine de l'imagerie par fluoroscopie, le rapport contraste sur bruit est souvent inférieur à 1. La présence d'artéfacts apparaît donc comme étant inévitable et a limité jusque là son utilisation à des images dont le rapport contraste sur bruit est bien plus élevé.

Un but de l'invention est notamment de proposer un traitement de débruitage adapté aux images fluoroscopiques dont le rapport contraste sur bruit est relativement faible. En particulier, le traitement de débruitage doit présenter une forte capacité de débruitage, engendrer un faible coût de calcul, limiter la perte des objets d'intérêt dans l'image, et ne pas introduire d'artéfacts. A cet effet, l'invention a pour objet un procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X, chaque image étant formée dans le domaine spatial d'une matrice de pixels, chaque pixel ayant une valeur représentative d'un niveau de signal, le procédé comprenant les étapes successives suivantes pour chaque image : ^ appliquer un filtrage temporel sur l'image acquise à un instant n, l'application dudit filtrage temporel comportant les sous-étapes suivantes pour chaque pixel de l'image : - déterminer une différence entre la valeur du pixel considéré dans l'image acquise à l'instant n et la valeur du pixel correspondant dans l'image acquise à un instant n-1 sur laquelle le filtrage temporel a été appliqué, - si la différence est inférieure à un premier seuil prédéterminé, 15 corriger la valeur du pixel correspondant avec un premier ensemble de paramètres, - si la différence est supérieure à un deuxième seuil prédéterminé, corriger la valeur du pixel correspondant avec un deuxième ensemble de paramètres, 20 ^ appliquer un filtrage spatial sur l'image acquise à l'instant n, l'application dudit filtrage spatial comportant les sous-étapes suivantes : - transformer l'image acquise à l'instant n du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de 25 coefficients, - seuiller les coefficients de l'image par une fonction de seuillage, la fonction de seuillage annulant les coefficients inférieurs à un troisième seuil prédéterminé, et conservant ou ajustant les coefficients supérieurs au troisième seuil prédéterminé, 30 - transformer l'image dont les coefficients ont été seuillés du domaine des curvelets au domaine spatial par une transformation en curvelets inverse.

L'étape de transformation de l'image acquise à l'instant n dans le domaine des curvelets utilise avantageusement une transformée en curvelets à 6 ou 9 échelles, et/ou 16 orientations.

Selon une première forme particulière de réalisation, l'étape de transformation de l'image acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par USFFT. Selon une deuxième forme particulière de réalisation, l'étape de 10 transformation de l'image acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par wrapping. La fonction de seuillage du filtrage spatial peut être une fonction de seuillage dur, c'est-à-dire une fonction de seuillage selon laquelle les 15 coefficients sont annulés s'ils sont inférieurs au troisième seuil prédéterminé, et conservés sinon. Le troisième seuil prédéterminé TS du filtrage spatial peut être déterminé par les étapes successives suivantes : 20 ^ créer une image Y de mêmes dimensions que l'image acquise à l'instant n et de valeur uniforme 1 sur tous les pixels, ^ transformer l'image Y par une transformée de Fourier, ^ normaliser l'amplitude F du pic dans la transformée de Fourier de l'image Y par la relation : 6 F IFI _ 5 M x P 25 où M et P désignent les dimensions de l'image Y, ^ transformer la transformée de Fourier de l'image Y du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients ci,k, 30 ^ déterminer une norme c i des coefficients c l,k par la relation : Y z y _ El cl,l,k M x P où M et P désignent les dimensions de l'image Y et l désigne les orientations des curvelets, ^ déterminer le seuil prédéterminé TS par la relation : TS=r.c j.Q où o- est l'écart type du bruit dans l'image acquise à l'instant n et sur laquelle le filtrage temporel a été appliqué, et où le facteur r prend la valeur 0 pour l'échelle la plus grossière des curvelets, la valeur 4 pour l'échelle la plus fine, et la valeur 1 pour les autres échelles.

Selon une forme particulière de réalisation, l'étape d'application du filtrage temporel sur l'image acquise à l'instant n comprend, à la suite de la sous-étape de détermination de la différence entre la valeur du pixel considéré dans l'image acquise à l'instant n et la valeur du pixel correspondant dans l'image acquise à l'instant n-1 sur laquelle le filtrage temporel a été appliqué, la sous-étape supplémentaire suivante : ^ si la différence est comprise entre le premier seuil prédéterminé et le deuxième seuil prédéterminé, corriger la valeur du pixel correspondant avec un troisième ensemble de paramètres.

En résumé, le procédé de réduction de bruit selon l'invention comporte successivement un filtrage temporel et un filtrage spatial par transformation en curvelets. L'ordre dans lequel sont effectués les filtrages présente une importance en ce qui concerne la coloration du bruit. En effet, le filtrage spatial par transformation en curvelets introduit des artéfacts s'il n'est pas précédé du filtrage temporel. Par ailleurs, l'invention a notamment pour avantage de ne pas nécessiter d'a priori de formes pour le filtrage spatial.

L'invention sera mieux comprise et d'autres avantages apparaîtront à la lecture de la description qui va suivre, faite en regard de 30 dessins annexés sur lesquels : - la figure 1 représente un exemple de points pris en compte pour déterminer le contraste d'un objet dans une image fluoroscopique ; - la figure 2 représente des étapes possibles du procédé de réduction de bruit selon l'invention ; - les figures 3A et 3B représentent par des courbes une fonction de seuillage dur et une fonction de seuillage doux, respectivement ;

- la figure 4 illustre, par un graphique, la performance de ce procédé de réduction de bruit pour une séquence d'images fluoroscopiques représentant un objet ayant un mouvement de translation et de rotation. Pour la suite de la description, on considère une séquence d'images fluoroscopiques xn acquises par un détecteur à rayons X, où n désigne les différents instants d'acquisition dans la séquence. Chaque image est composée d'une matrice de pixels ayant chacun une valeur xn(m,p) représentative d'un niveau de signal, avec m et p désignant respectivement la ligne et la colonne de la matrice. La valeur est par exemple une valeur numérique comprise entre 0 et 255. La valeur numérique est ainsi codée sur 8 bits. Elle peut également être codée sur 14 bits. Le rapport contraste sur bruit (CNR) est défini comme suit : 'contraste - bruit' CNR = a où le bruit est de préférence déterminé sur une zone de l'image où il n'y a pas d'objets d'intérêt, où le contraste est de préférence déterminé à partir de points correspondant à un objet d'intérêt dans l'image, et où Q désigne l'écart type du bruit, c'est-à-dire l'écart type des valeurs de pixels correspondant à une zone de l'image ne comportant pas d'objets d'intérêt. La figure 1 représente les points considérés (en blanc) pour déterminer le contraste d'une endoprothèse vasculaire, mieux connue sous la dénomination anglo-saxonne de "stent".

La figure 2 représente des étapes possibles pour le procédé de réduction de bruit dans la séquence d'images fluoroscopiques. Le procédé est décrit pour l'image xn acquise à l'instant n. Il peut cependant être appliqué successivement sur chacune des images de la séquence.

Dans une première étape 21, l'image xn est traitée par une chaîne de prétraitement. Le prétraitement peut notamment comporter des regroupements de pixels adjacents dans l'image. Les pixels peuvent être regroupés par exemple en fonction de la résolution souhaitée ou du niveau de bruit souhaité. Le prétraitement peut également comporter une phase d'étalonnage du détecteur à rayons X. En particulier, le gain et l'offset à appliquer à chaque pixel du détecteur peuvent être déterminés. En outre, l'étalonnage peut également comprendre une identification des pixels prenant une valeur aléatoire ou une valeur constante quel que soit le rayonnement reçu. Ces pixels défectueux peuvent être corrigés en interpolant leur valeur par les valeurs des pixels voisins. Le prétraitement peut aussi comporter une remise à l'échelle des images, c'est-à-dire un ajustement global de la valeur de tous les pixels de chaque image, de manière à ce que la valeur moyenne des pixels soit constante pour toutes les images de la séquence. Enfin, le prétraitement peut comporter une sous-étape d'adaptation des images. Étant donnée la physique des rayons X, le bruit présent dans chaque image xn est un bruit poissonien. Or un tel bruit est inadapté aux traitements de réduction de bruit. De plus, chaque image xn est une image représentative des atténuations des rayons X, alors qu'une image des épaisseurs radiologiques est plus naturelle à analyser. La formule 15 d'Ascombe permet d'obtenir une image des épaisseurs radiologiques avec un bruit pseudo-gaussien : xn (m, p) _ où xo,n(m, p) est la valeur du pixel (m, p) avant la sous-étape d'adaptation de l'image xn. Dans une deuxième étape 22, l'image xn est filtrée par un filtre 20 temporel. Le filtrage temporel est par exemple un filtrage appelé "Fluoroscopic Noise Reduction" (FNR). Ce filtrage est décrit notamment dans le brevet US 6,314,160 B1. Soit yn l'image xn obtenue après le filtrage temporel. On note 0(m, p), la différence entre la valeur xn(m, p) dans l'image dans l'image filtrée xn acquise à l'instant n et la valeur yn_1(m, p) 25 yn_1 correspondant à l'acquisition à l'instant n-1 : 0(m, p) = xn (m, p) - yn-1(m, p) Pour alléger les notations, on notera simplement xn la valeur d'un pixel de l'image xn, yn la valeur d'un pixel de l'image yn et 0 la différence entre les deux valeurs xn et yn_1 du pixel considéré. Le filtrage FNR se décline selon deux modes, appelés mode adaptatif et mode non adaptatif. Le mode 30 adaptatif est appelé ainsi du fait de sa prise en compte du mouvement des objets. Le filtrage FNR consiste à comparer pixel à pixel, la différence 0 à des seuils prédéterminés. Le premier seuil SI correspond à une différence 0 relativement faible. Le deuxième seuil S2 correspond à une différence A relativement forte. Le filtrage FNR est alors paramétré par des coefficients a et [3 en fonction du résultat de cette différence A et du mode considéré. On note a la force du filtre : ~2 a= y a + Qx où a-x est l'écart type dans l'image xn et où Qy est l'écart type dans l'image Yn-1 On définit aussi les constantes AA, BB, CC et DD, comme suit : ^ AA = 1=a 2.ax ^ BB=2a-1 Q(1-a) ^ CC = 2 ^ DD = i Q(1- a)

Pour le mode adaptatif, les seuils SI et S2 sont définis comme suit : ^ Ss=a" ^ S2=2Q et les coefficients a et R sont déterminés comme suit : ^ si IAI < S1, lIa = a 1(3=0 AA.AZ+BB.I oI +CC ^ si SI < 1AI < S2, a = lol /3=0

a=1-II ^ sinon, Pour le mode non adaptatif, les seuils SI et S2 sont définis comme suit : 20 ^ Ss = 26 ^ S2 Ja et les coefficients a et [3 sont déterminés comme suit : ^ si 14I < S1, fa= a = 0 a _ a.°Z ^ si SI < IA1 < S2, 32.a.)2 4 fi' = 3a.a a=1 25 ^sinon, _ a 4.v ~ 3 va Le filtrage FNR proprement dit consiste à déterminer une nouvelle valeur yn pour chaque pixel de l'image à partir de sa valeur xn et des coefficients a et R: yn=a.xn+(1-a)yn-1+fi' Autrement dit, la valeur xn de chaque pixel est corrigée à partir des coefficients a et R. Sur la figure 2, l'étape 22 de filtrage FNR est représentée pour le mode adaptatif. A titre d'exemple, dans une première sous-étape 221, la valeur absolue de la différence 0 est comparée pour chaque pixel xn au seuil S1, à savoir l'écart type a. Si la valeur absolue de la différence à est inférieure au seuil S1, le pixel considéré représente probablement un niveau de bruit et un filtrage dit fort est appliqué dans une sous-étape 222 avec les paramètres correspondants. Sinon, dans une sous-étape 223, la valeur absolue de la différence 0 est comparée pour chaque pixel xn au seuil S2, à savoir 2a. Si la valeur absolue de la différence 0 est supérieure au seuil S2, le pixel considéré correspondant probablement à un point de l'objet ayant bougé et un filtrage dit faible est appliqué dans une sous-étape 224 avec les paramètres correspondants. Enfin, si la valeur absolue de la différence 0 est comprise entre les seuils S1 et S2, une incertitude existe. Afin d'obtenir une fonction continue entre le filtrage fort et le filtrage faible, un filtrage moyen est appliqué dans une sous-étape 225 avec les paramètres correspondants. Dans une troisième étape 23, l'image yn obtenue après filtrage temporel est filtrée par un filtre spatial basé sur la transformée en curvelets. La théorie des curvelets utilise le principe des espaces multirésolution ou multi-échelle, une échelle correspondant à un niveau de partition de l'espace. L'espace est partitionné en curvelets dont les enveloppes sont de dimensions 2j x 2j/2, où j est un entier positif désignant l'échelle. Chaque curvelet dépend du nombre d'échelles considérées, de sa position dans l'image et de son orientation. Pour chaque échelle, la meilleure approximation possible de l'image est calculée. D'une échelle à l'autre, seules les améliorations obtenues sont représentées. Le nombre d'échelles détermine la qualité de la reconstruction de l'image. Avec les curvelets, un nombre d'échelles compris entre 6 et 10 permet de fournir des résultats très satisfaisants en termes de rapport contraste sur bruit.

Dans une première sous-étape 231 de l'étape de filtrage spatial 23, l'image yn est transformée du domaine spatial au domaine des curvelets par une transformation en curvelets. De préférence, une transformation en curvelets discrète est utilisée. Une telle transformation est notamment décrite dans E. Candès, Demanet L., Donoho D. and Ying L., "Fast Discrete Curvelet Transforms", Multiscale Model. Simul., vol. 5, no. 3, 2006, pp. 861-899. Sous leur forme discrète, les curvelets ne peuvent pas être échantillonnées selon les grilles de l'espace de Fourier. Le document de Donoho et al propose ainsi deux solutions pour translater les curvelets à échelle et orientation données sur la grille spatiale. Une première solution pour transformer l'image dans le domaine des curvelets discrètes est appelée transformée en curvelets discrète via USFFT, où USFFT est le sigle anglo-saxon pour "unequally-spaced fast Fourier transform" (USFFT), que l'on peut traduire en français par "transformée de Fourier rapide pour données à répartition irrégulière". Dans cette solution, la grille de translation est inclinée de manière à être alignée dans la direction des curvelets à échelle et angle donnés. Une deuxième solution est appelée transformée en curvelets discrète via wrapping, où le terme "wrapping" peut être traduit en français par "enroulement". Dans cette solution, la grille de translation est identique pour tous les angles d'un même quadrant à une échelle donnée. Cela implique que le nombre d'angles utilisés pour la transformée soit un multiple de 4. La transformée en curvelets via USFFT est plus fidèle à la théorie continue des curvelets que la transformée en curvelets via wrapping. Cependant, son coût de calcul est plus important. Dans les deux solutions, la transformée en curvelets de l'image yn peut être représentée par un ensemble de coefficients c°(j, 1, k), où l'exposant D fait référence à la transformée discrète, et où j, l et k désignent respectivement l'échelle, l'orientation et la position de la curvelet considérée. Selon une forme préférentielle de réalisation du filtre spatial, le nombre d'échelles utilisées est égal à 9. Selon une autre forme de réalisation, le nombre d'échelles utilisées est égal à 6. Cette dernière forme de réalisation implique un nombre plus faible de coefficients à calculer et donc une durée de traitement inférieure. Par ailleurs, le nombre d'orientations utilisées est avantageusement égal à 16. Ce nombre présente un bon compromis entre la résolution en orientation et le coût de calcul des curvelets, notamment pour des images réelles dont les courbes nécessitent un grand nombre d'orientations pour les approcher. De plus, le nombre 16 étant un multiple de 4, il facilite la transformée en curvelets par wrapping. Dans une deuxième sous-étape 232 de l'étape de filtrage spatial 35 23, les coefficients c°(j, 1, k) de l'image yn sont seuillés par une fonction de seuillage. L'étape de seuillage 232 est basée sur l'hypothèse selon laquelle certains coefficients cD(j, 1, k) représentent le bruit alors que d'autres contiennent l'information utile. L'étape de seuillage 232 consiste donc à identifier et à supprimer les coefficients ne contenant pas d'informations utiles. La fonction de seuillage se caractérise d'une part par sa forme et d'autre part par son seuil T. Il existe principalement deux fonctions de seuillage, la fonction de seuillage dit doux et la fonction de seuillage dit dur. En simplifiant la notation des coefficients par en et en appelant en les coefficients seuillés, la fonction de seuillage dur peut être modélisée comme suit : 1cn sinon e _ n sinon La fonction de seuillage doux peut être modélisée comme suit : _ si Icni T en - IO cn - T. signe(cn) sinon La figure 3A représente par une courbe 31 la fonction de seuillage dur et la figure 3B représente par une courbe 32 la fonction de seuillage doux. La fonction de seuillage doux a une capacité de débruitage plus importante que la fonction de seuillage dur. En revanche, elle modifie la valeur des pixels en fonction du seuil T choisi. Il en résulte une perte de précision sur certains détails de l'image. La détermination du seuil T est une étape cruciale car elle détermine la performance du filtrage spatial, tant en ce qui concerne la capacité de débruitage que l'introduction d'artéfacts.

Selon une première forme de réalisation, le seuil T est déterminé comme étant le seuil universel Tu décrit notamment dans D.L. Donoho and Johnstone, I.M., "Ideal spatial adaptation via wavelet shrinkage", Biometrica, vol. 81, 1994, pp. 425-455. Ce seuil Tu est défini par la relation suivante : logN Tu _ 2 N où N est le nombre d'échantillons, c'est-à-dire ici le nombre de pixels dans l'image yn. Le seuil universel TU n'est pas optimal et peut introduire des artéfacts. Néanmoins, des seuils proportionnels au seuil universel TU peuvent présenter de bons résultats si le facteur de proportionnalité est bien choisi. Selon une deuxième forme de réalisation, le seuil T est déterminé par la méthode décrite dans J.L. Starck, Candès I.J. and Donoho L., "The Curvelet Transform for Image Denoising", IEEE Transactions on Image Processing, vol. 11, no. 6, 2002, pp. 670-684. Cette méthode consiste à créer une image, notée Y, de mêmes dimensions que l'image yn à débruiter et de valeur uniforme 1 sur tous les pixels. L'image Y est transformée par la transformée de Fourier. La fréquence de coordonnées (0,0) dans l'image Y est centrée et son amplitude est normalisée par les dimensions de l'image. Soit F l'amplitude de la fréquence centrale. Son amplitude normalisée, notée IFI, est par exemple réalisée selon la relation suivante : F IFI = 5 VM x p où M et P désignent les dimensions de l'image Y, c'est-à-dire également de l'image xn et de l'image yn.

On obtient ainsi une image avec un pic de Dirac centré sur l'image, dont l'amplitude est fonction des dimensions de l'image. Le seuil de cette méthode, noté TS dépend de l'échelle j et de l'orientation 1. Il est déterminé en fonction de la norme des coefficients ci,k des curvelets de l'image Y pour une échelle j et une orientation l données. Cette norme c est déterminée par la relation suivante : 2 y El cjY,I,k - M x P

La fonction de seuillage est une fonction de seuillage dur, définie comme suit : Cn = 1cn sinon Le facteur multiplicatif r est déterminé en fonction de l'échelle j. Pour l'échelle la plus fine (par exemple j = 9), le facteur r est égal à 4. Pour les autres échelles, le facteur r est égal à 3. Ces facteurs présentent une bonne capacité de débruitage pour des images présentant un fort rapport contraste sur bruit. Cependant, pour des images fluoroscopiques, des artéfacts sont introduits. Ces artéfacts peuvent être évités en adaptant le facteur r. En particulier, pour l'échelle la plus grossière (j = 1), le facteur r peut être égal à 0. Il n'y a donc pas de seuillage. Pour l'échelle la plus fine, le facteur r peut encore être égal à 4. Et pour les autres échelles, le facteur r peut être égal à 1. Dans une troisième sous-étape 233 de l'étape de filtrage spatial 23, l'image yn est transformée du domaine des curvelets au domaine spatial (matriciel) par une transformation en curvelets inverse. A cet effet, les 0si Icnl <r.ci .6 coefficients seuillés sont utilisés. L'image obtenue après filtrage temporel et filtrage spatial est notée zn.

La figure 4 illustre par un graphique la performance du procédé de réduction de bruit selon l'invention pour une séquence d'images représentant un stent ayant un mouvement de translation et de rotation. Sur le graphique, l'axe des abscisses représente le CNR moyen de la séquence d'images xn et l'axe des ordonnées représente le CNR moyen de la séquence d'images obtenue après un traitement de réduction de bruit. Une première courbe 41 représente le CNR moyen de la séquence d'images Zn obtenue par le procédé de réduction de bruit selon l'invention, c'est-à-dire par un filtrage temporel suivi d'un filtrage spatial par transformation en curvelets. Une deuxième courbe 42 représente le CNR moyen de la même séquence d'images après traitement par le même filtrage temporel suivi d'un filtrage spatial par transformation en ondelettes. Le graphique montre que pour une séquence d'images à traiter dont le CNR est relativement faible (0,5 - 0,7 - 1), le procédé de réduction de bruit selon l'invention présente une amélioration notable en termes de CNR. La réduction de bruit est encore bien meilleure pour une séquence d'images dont le CNR est plus élevé (1,5 - 2 - 2,5). De plus, le procédé de réduction de bruit selon l'invention permet de réduire significativement l'introduction d'artéfacts.

Claims

REVENDICATIONS1. Procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X, chaque image (xn) étant formée dans le domaine spatial d'une matrice de pixels, chaque pixel ayant une valeur représentative d'un niveau de signal, le procédé comprenant les étapes successives suivantes pour chaque image (xn) : ^ appliquer un filtrage temporel (22) sur l'image (xn) acquise à un instant n, l'application dudit filtrage temporel comportant les sous-étapes suivantes pour chaque pixel de l'image : - déterminer (221, 223) une différence entre la valeur du pixel considéré dans l'image (xn) acquise à l'instant n et la valeur du pixel correspondant dans l'image (yn_1) acquise à un instant n-1 sur laquelle le filtrage temporel (22) a été appliqué, - si la différence est inférieure à un premier seuil prédéterminé, corriger (222) la valeur du pixel correspondant avec un premier ensemble de 15 paramètres, - si la différence est supérieure à un deuxième seuil prédéterminé, corriger (224) la valeur du pixel correspondant avec un deuxième ensemble de paramètres, ^ appliquer un filtrage spatial (23) sur l'image (yn) acquise à l'instant n, 20 l'application dudit filtrage spatial comportant les sous-étapes suivantes : - transformer (231) l'image (yn) acquise à l'instant n du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients, 25 - seuiller (232) les coefficients de l'image (yn) par une fonction de seuillage, la fonction de seuillage annulant les coefficients inférieurs à un troisième seuil prédéterminé, et conservant ou ajustant les coefficients supérieurs au troisième seuil prédéterminé, - transformer (233) l'image (yn) dont les coefficients ont été 30 seuillés du domaine des curvelets au domaine spatial par une transformation en curvelets inverse.
2. Procédé selon la revendication 1, dans lequel l'étape (231) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets à 6 échelles.
3. Procédé selon la revendication 1, dans lequel l'étape (231) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets à 9 échelles.
4. Procédé selon l'une des revendications précédentes, dans 10 lequel l'étape (231) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets comprenant 16 orientations.
5. Procédé selon l'une des revendications précédentes, dans 15 lequel l'étape (231) de transformation de l'image (A) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par USFFT.
6. Procédé selon l'une des revendications 1 à 4, dans lequel 20 l'étape (231) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par wrapping.
7. Procédé selon l'une des revendications précédentes, dans 25 lequel la fonction de seuillage est une fonction de seuillage dur, c'est-à-dire que les coefficients sont annulés s'ils sont inférieurs au troisième seuil prédéterminé, et conservés sinon.
8. Procédé selon la revendication 7, dans lequel le troisième 30 seuil prédéterminé TS est déterminé par les étapes successives suivantes : ^ créer une image Y de mêmes dimensions que l'image (yn) acquise à l'instant n et de valeur uniforme 1 sur tous les pixels, ^ transformer l'image Y par une transformée de Fourier, ^ normaliser l'amplitude F du pic dans la transformée de Fourier de 35 l'image Y par la relation :6 F 1FI _ 5 M x P où M et P désignent les dimensions de l'image Y, ^ transformer la transformée de Fourier de l'image Y du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients Cjl,k, ^ déterminer une norme ê I des coefficients c jl,k par la relation : Y y zI Cl I,k a sil- M x P où M et P désignent les dimensions de l'image Y et l désigne les orientations des curvelets, ^ déterminer le seuil prédéterminé TS par la relation : TS = r. ci. Q où or est l'écart type du bruit dans l'image (A) acquise à l'instant n et sur laquelle le filtrage temporel (22) a été appliqué, et où le facteur r prend la valeur 0 pour l'échelle la plus grossière des curvelets, la valeur 4 pour l'échelle la plus fine, et la valeur 1 pour les autres échelles.
9. Procédé selon l'une des revendications précédentes, dans lequel l'étape (22) d'application du filtrage temporel sur l'image acquise (xn) à l'instant n comprend, à la suite de la sous-étape (221, 223) de détermination de la différence entre la valeur du pixel considéré dans l'image acquise (xn) à l'instant n et la valeur du pixel correspondant dans l'image (yn_1) acquise à l'instant n-1 sur laquelle le filtrage temporel (22) a été appliqué, la sous-étape supplémentaire suivante : ^ si la différence est comprise entre le premier seuil prédéterminé et le deuxième seuil prédéterminé, corriger (225) la valeur du pixel correspondant avec un troisième ensemble de paramètres.25