FR2953964A1 - Procede de traitement d'images obtenues par tomographie ou tomosynthese a faible nombre de projections - Google Patents

Description

DOMAINE TECHNIQUE GENERAL L'invention concerne la reconstruction d'une image bidimensionnelle (2D) ou tridimensionnelle (3D) d'un objet, par exemple une partie d'une région d'intérêt d'un patient à partir d'un ensemble de projections respectivement monodimensionnelles ou bidimensionnelles de la région d'intérêt, obtenues pour différentes positions d'un système d'imagerie autour de la région d'intérêt. Et l'invention trouve notamment application en imagerie médicale par reconstruction tomographique ou de tomosynthèse à faible nombre de projections (en anglais « Few-View Tomography »).

ETAT DE LA TECHNIQUE La tomographie permet de réaliser des images de coupes (en anglais, « slice ») d'une région d'intérêt d'un objet en acquérant des projections. On a illustré sur la figure 1 schématiquement l'acquisition d'images 2D d'un organe et la reconstruction en une image 3D de cet organe par tomosynthèse. Des rayons X, R, issus d'une source S sont émis selon différentes angulations (1, ..., i, ..., n) vers l'organe O. Après avoir traversé l'organe, ils sont détectés par le détecteur Det formant un ensemble de projections 2D s,, ..., s;,..., sn. Il est à noter qu'il y a autant de projections 2D acquises que d'angulations considérées. Si l'ensemble des angulations couvre au moins un demi-cercle, on parle de tomographie. Quand cette condition ne peut être vérifiée (par exemple : arc de cercle inférieur à 180 degrés, ligne droite), on parle de tomosynthèse.

L'acquisition est mise en oeuvre par un détecteur Det situé en regard d'une source de rayons X, par exemple une caméra numérique. Une application de la tomographie est la détection et la caractérisation d'une lésion dans un organe, par exemple une lésion cancéreuse. Dans le cas du sein, on utilise la tomosynthèse.

Les projections 2D acquises sont utilisées pour reconstruire une image 3D de l'objet. Cette image 3D est plus précisément une cartographie 3D des coefficients d'atténuation aux rayons X du milieu traversé. C'est au moyen de cette cartographie que le praticien radiologue interprète cette image en fonction des différences de contraste observées. La reconstruction de l'image 3D est une opération très couteuse et complexe, en particulier si elle met en oeuvre un algorithme itératif. On connait un procédé de reconstruction 3D itératif. Ce procédé est basé sur une expression discrète et matricielle du problème de reconstruction tomographique. De manière plus précise le problème se modélise par l'équation suivante : Rf = s où s est un vecteur des projections acquises, R est un opérateur de projection qui modélise le système d'imagerie tomographique et f est l'image 3D de l'objet à reconstruire. Le problème de la reconstruction tomographique est de déterminer f en connaissant s et R . La reconstruction exacte de l'image 3D f requiert en théorie un ensemble de mesures au moins aussi grand que le nombre d'inconnues formant l'image 3D acquises suivant les contraintes de la tomographie, c'est-à-dire des mesures avec un pas d'échantillonnage très faible et une rotation de 180 degrés autour de l'objet. Ce cas peut en pratique se révéler impossible pour les raisons suivantes. D'une part dans des examens de type tomosynthèse, il n'est pas possible d'obtenir des mesures obtenues à 180° autour de l'objet à imager et d'autre part, plus le pas d'échantillonnage est petit plus le nombre de mesures est important, ce qui augmente la durée des examens et donc la dose de rayons X à appliquer.

Ainsi, on rencontre très souvent en tomographie le cas où la taille de l'image 3D est très largement supérieure à la taille des données, ce qui implique la résolution d'un système ayant un très grand nombre de solutions possibles.

Le choix du procédé de reconstruction permet de trouver une image 3D comme solution acceptable du point de vue de l'application visée. En outre, les procédés connus de reconstruction dans l'espace des images 3D ont la forme d'un algorithme itératif défini par son itération A qui génère une reconstruction f n' de la solution à partir d'une reconstruction f'" telle que : = A[f (n-" ] . A partir d'une reconstruction initiale notée f de la série des reconstructions f , on obtient la reconstruction à l'itération N comme : f (N) = A[ f (N-" ] A[A[f (N 2) ]] = A2[f~N 2)]ù... = A`[f(N a)] ... = AN[f `0']• Cet algorithme requiert le stockage des étapes d'initialisation, de reconstructions intermédiaires et d'archivage de la reconstruction, soit une ou plusieurs images 3D selon que le stockage est fait simultanément ou successivement, et représente une quantité d'informations générées très supérieure à la quantité d'informations acquises (un faible nombre de vues).

Tel que connu, des procédés mettent en jeu des opérations discrètes et matricielles (donc linéaires). Certains sont associés à la résolution du système linéaire Rf = s, d'autres sont associés à la minimisation de l'erreur quadratique Rf ûsr et du système associé R`Rf = R`s où R` est la matrice transposée de R , d'autres encore généralisent l'approche précédente à la 25 résolution d'un système RRf = Rs où R est une matrice telle que RRf est une image 3D. On pourra se référer au document « CIARLET, P.G., Introduction à l'Analyse Numérique Matricielle et à l'Optimisation, Masson, Paris, 1982 ». Plus récemment, il a été proposé de calculer de manière itérative, une fois pour toutes, un filtre obtenu par reconstruction itérative d'un objet test et d'appliquer à tout objet à imager ce filtre aux projections acquises afin d'obtenir le volume reconstruit de l'objet sans calcul itératif. On pourra se référer au document US 7 447 295. Par ce biais, la puissance de calcul requise est diminuée par rapport aux procédés itératifs. Toutefois, la substitution de la reconstruction itérative par un filtrage engendre une perte d'informations. Il s'agit donc d'un procédé approché.

PRESENTATION DE L'INVENTION L'invention propose quant à elle un procédé itératif qui ne nécessite pas le stockage d'images 3D intermédiaires.

Selon un premier aspect, l'invention concerne un procédé de traitement d'images obtenues par tomographie ou tomosynthèse comprenant : une étape d'acquisition d'une pluralité d'images de projection 2D d'un objet, l'acquisition étant définie par Rf = s où s est un vecteur des projections acquises, R est un opérateur de projection qui modélise le système d'imagerie et f est l'image 3D de l'objet à reconstruire connaissant R et s ; une étape de traitement des images de projection 2D acquises. Le procédé selon le premier aspect de l'invention est caractérisé en ce que le traitement des images de projection consiste en une application d'un algorithme itératif défini par son itération V qui génère à chaque itération n =1,...,N au moins un ensemble d'images de projection traitées p L'itération V est définie de manière à ce qu'à chaque itération l'image 3D de l'objet soir une fonction linéaire des images de projection 2D traitées selon la propriété f 0n' = Rp(n' , où R est une matrice telle que RRf est une image 3D.

Compte tenu de la contrainte associée à l'algorithme itératif V , seul un ensemble d'images 2D de projection est manipulé. Ceci évite d'avoir à stocker à quelque moment du procédé que ce soit (initialisation n = 0, calcul intermédiaire 0 < n < N , archivage n = N) l'intégralité d'une image 3D dont la taille est très supérieure à l'ensemble des mesures acquises. En outre, la relation r' = Rp(' est établie de façon suffisamment simple pour n'être calculée que lors de l'étape de visualisation. En particulier, cette relation requiert l'ensemble des images 2D traitées p , mais seulement le calcul et le stockage de la partie de r' qui est visualisée.

Comme il n'est ni possible ni nécessaire de visualiser tous les plans de coupes d'un volume simultanément, mais seulement successivement, il n'est jamais nécessaire de mettre en oeuvre des moyens de stockage couvrant l'intégralité de la taille de l'image 3D r'. En résumé, le stockage de l'intégralité d'une image 3D n'est nécessaire ni lors de l'étape de reconstruction, ni lors de l'archivage, ni lors de l'étape de visualisation. Pour des algorithmes connus définis par leur itération A associée à la matrice RR dans le domaine des images 3D qui requiert donc la manipulation, le stockage et l'archivage d'au moins une image 3D stockée dans son intégralité, on montre qu'il existe un algorithme défini par l'itération V qui génère une suite d'images traitées sous forme d'images de projection 2D et une contrainte linéaire R permettant de générer une reconstruction de l'objet à partir des projections traitées telles que les reconstructions obtenues par V et k sont mathématiquement équivalentes à celle générées par A, sans qu'il ne soit jamais nécessaire de stocker l'intégralité d'une image 3D.

De plus, le temps de calcul est proportionnel à l'échantillonnage de chaque image 3D à reconstruire alors que la taille occupée par le stockage des images de projections 2D est constante. Il est donc possible de régler l'échantillonnage de l'image 3D reconstruite au-delà des capacités de stockage d'une unité de traitement considérée et donc d'implémenter sur cette unité une itération V équivalente à une itération A non-implémentable sur cette unité, mais seulement sur une autre unité au stockage considérablement augmenté.

D'autres aspects du procédé selon le premier aspect de l'invention sont les suivants : le procédé comprend une étape de reconstruction 3D (S2) de l'objet à déterminer f (N) = 17p(N) où N est le nombre d'itérations fixées ; à chaque itération n =1,...,N, l'algorithme V génère un ensemble d'images de projection traitées p(o = V [p(n-»] , l'itération V étant définie pour tout vecteur p dans l'espace des images de projection 2D de la manière suivante V [p] = (I ûpRR )p + ps où p(°) = ps ou bien p(°) = 0 , s est l'ensemble des images de projection 2D acquises, p est un paramètre fixé pour que 11I ûpRR < 1 et ainsi garantisse la convergence du procédé et I est la matrice carrée identité.

à chaque itération, l'algorithme V génère les ensembles d'images de projection traitées (p (n)q(n) t(n))=V[(p(nù1) q(nù1) t(n-l))] suivant la méthode du gradient conjugué, typiquement donné par :

a'" = ù(q(n 1) Rtt(n 1))/(RRtt(n 1) RRtt p(n) ù p(n l) + a(n 1)t(n 1) q(~) = q(~ 1) + a(n 1) Rtt(n 1) /3(n 1) =(q() , RRtq~(nn)))/(q(n 1)~ ,RR t(n) = q(n) + /~(n 1)t(n 1) i~L~

où p , q , t sont tels que p(°) est arbitraire, typiquement nul ou égal à s le vecteur des images de projection 2D acquises, et q(0) = t(0) = RRt p(0) û s . à chaque itération, l'algorithme V®(n) génère un ensemble d'images de projection t' (n) = V®(n) ' '1) ] = V0(n)VO(n-1) [p''] = V0(n)VO(n-1) ...V0(1) [p'') ] avec V®(,,)=(Ion ù p®()I®( )RR)p+po(0Io(os+I,(,,)p où O(n) et '(n) sont respectivement un premier et un second sous-ensemble des index des mesures projectives acquises, définis pour chaque itération n, les deux sous-ensembles étant disjoints, leur réunion indexant la totalité du vecteur s des images 2D de projection acquises, et Ion) et I,(n) sont des matrices dont les coefficients sont tous nuls sauf pour les indices de la diagonale appartenant à O(n) et '(n) respectivement où ils valent 1, de sorte que p = I o(n) p + I,(ä) p , où p(0) = p®(Q)I®(Q)s ou bien p(°' = 0, et po(n) est un paramètre fixé pour que 11I®O ù p®( )RR <1 et ainsi garantisse la convergence du procédé . le procédé comprend une extraction d'une coupe issue de l'image 3D reconstruite ; le procédé comprend un réglage arbitrairement fin de l'échantillonnage de l'image 3D. Selon un second aspect, l'invention concerne un système d'imagerie médicale comprenant des moyens pour la mise en oeuvre d'un procédé selon le premier aspect de l'invention. Selon un troisième aspect l'invention concerne un produit programme 20 d'ordinateur comprenant des instructions de code de programme pour l'exécution des étapes du procédé selon le premier aspect de l'invention.

PRESENTATION DES FIGURES D'autres caractéristiques et avantages de l'invention ressortiront 25 encore de la description qui suit laquelle est purement illustrative et non limitative et doit être lue en regard des dessins annexés sur lesquels outre la figure 1 déjà discutée : - la figure 2 illustre schématiquement un système d'imagerie médicale pour la mise en oeuvre du procédé de l'invention ; la figure 3 illustre schématiquement les étapes du procédé de l'invention ; la figure 4 illustre schématiquement, en développé, une étape du procédé de l'invention. DESCRIPTION DETAILLEE DE L'INVENTION

Système d'imagerie médicale Sur la figure 2, on a illustré, schématiquement, un système d'imagerie 10 médicale 1 pour l'acquisition d'images 2D de projection pour la reconstruction d'une image 3D d'un organe. Un tel système peut être un appareil mammographique pour la détection et la caractérisation de lésions dans le cas du dépistage, du diagnostic et du traitement du cancer du sein. 15 Le système d'imagerie médical 1 comprend un système d'acquisition d'images 3, un système de traitement d'images 5 et un système d'affichage 4. Le système d'acquisition 3 permet d'acquérir une pluralité de projection 2D d'une région d'intérêt ù d'un organe ù d'un patient. Le système 20 d'acquisition 3 est notamment constitué par un détecteur De situé en regard d'une source de rayons X. Le détecteur est par exemple une caméra numérique. Le système d'acquisition est par exemple un système d'acquisition par rayons X, ce dernier comprenant tout moyen connu permettant l'émission de rayons X sur l'objet 2 et l'acquisition d'images 25 résultantes. Le système d'affichage 4 peut être intégré dans le système d'acquisition d'images 3 ou le système de traitement d'images 5, ou être séparé du système d'acquisition 3 et du système de traitement 5. Le système d'affichage 4 est par exemple un écran d'ordinateur, un moniteur, 30 un écran plat, un écran plasma ou tout type de dispositif d'affichage connu5 du commerce. Le système d'affichage 4 permet à un praticien de contrôler la reconstruction et/ou l'affichage des images 2D acquises. Le système de traitement 5 est adapté à la mise en oeuvre de procédés de traitement (par exemple de reconstruction d'une image 3D à partir d'images 2D). Le système de traitement 5 peut être intégré dans le système d'acquisition 3 d'images ou être séparé du système d'acquisition 3 d'images. Le système de traitement 5 est par exemple un/des ordinateur(s), un/des processeur(s), un/des microcontrôleur(s), un/des micro-ordinateur(s), un/des automate(s) programmable(s), un/des circuit(s) intégré(s) spécifique(s) d'application, d'autres circuits programmables, ou d'autres dispositifs qui incluent un ordinateur tel qu'une station de travail. Le système de traitement 5 est couplé à des moyens mémoires 6 qui peuvent être intégrés ou séparés du système de traitement 5. Ces moyens peuvent être formés par un disque dur ou SSD, ou tout autre moyen de stockage amovible et ré-inscriptible (clés USB, cartes mémoires etc.). Ces moyens mémoires peuvent servir à stocker une image 3D de la zone de l'organe visualisée comme une image 2D acquise ou traitée. Il peut s'agir d'une mémoire ROM/RAM du système de traitement 5, une clé USB, une carte mémoire, une mémoire d'un serveur central. Le système de traitement 5 peut comprendre un dispositif de lecture (non représenté) par exemple un lecteur de disquettes ou un lecteur de CD-ROM, pour lire les instructions du procédé de traitement (qui va être décrit dans la suite) d'un support d'instructions (non montré), comme une disquette ou un CD-ROM ou de manière plus générale par tout support de mémoire amovible ou encore via une connexion réseau. En variante, le système de traitement 5 peut comprendre un dispositif de connexion réseau (non représenté) filaire ou sans-fil. En variante, le système de traitement 5 exécute les instructions du procédé de traitement (qui va être décrit dans la suite) stockées dans des micrologiciels (non représentés).30 Description générale du procédé de traitement d'image Le procédé de traitement d'images comprend une étape SO d'acquisition d'une pluralité d'images de projection 2D si i=1,•••,M >1, M étant le nombre d'images de projection 2D acquises, d'un objet et une 5 étape de traitement S1 des images de projection 2D acquises. C'est l'optimisation de l'étape S1 de traitement qui permet de limiter les besoins de mémoire de stockage de la reconstruction 3D. De manière plus précise on applique un algorithme itératif sur les images de projection acquises pour obtenir, à chaque itération, des images 10 de projection 2D traitées p de sorte qu'à chaque itération n, l'image 3D de l'objet est une fonction linéaire des images de projection 2D traitées. En particulier l'algorithme est défini par son itération V qui génère un ensemble d'images de projection traitées p(n) = v[p,n-»]= v n [p,0' ] de manière à ce qu'à chaque itération n=1,•••,N, l'image 3D de l'objet soit une 15 fonction linéaire des images de projection 2D traitées, suivant la propriété : fin' = xp'n' . C'est cette relation qui montre qu'à chaque itération n il est possible de passer simplement d'un ensemble de projection 2D à une image 3D reconstruction de l'objet. On note que l'étape de traitement S1 comprend une étape 20 d'initialisation S10 et une étape S11 au cours de laquelle on applique l'itération V (voir la figure 4). L'étape de traitement ne requiert le stockage que des ensembles d'images (acquises et/ou traitées) au lieu de volumes. Ci-dessous on présente de manière détaillée trois algorithmes 25 implémentés dans le domaine des images 2D. On montre pour chacun de ces algorithmes, son équivalence avec un algorithme de reconstruction itérative connu implémenté dans le domaine des images 3D.

Premier algorithme Un premier algorithme classique dit des « approximations successives », consiste en l'itération A définie de la manière suivante A[f]=(IùpRR)f+pRs. où I est la matrice carrée identité dont les coefficients sont tous nuls sauf sur la diagonale où ils sont égaux à 1, k est une matrice telle que RRf est une image 3D, et p un paramètre fixé a priori pour que I û pkR < 1 pour garantir la convergence du processus. On connaît que l'algorithme converge vers une solution du système RRf = Rs . Si R est inversible, cette solution satisfait aussi Rf = s . Si R = R` , la suite des estimées converge vers une solution telle que limAn[f(0)]= f* =minRf û s2 . En pratique, on fixe un entier N suffisamment grand pour lequel on calcule f (") =A[f(n-"] pour n =1,•••,N avec f=Rp`°' avec p(0' =0 ou p`°'=ps. On constate que, dans cet algorithme standard, à chaque itération, des images 3D sont manipulées ( f,pRRf,pRs,A[f]) ce qui nécessite un stockage d'un nombre important de données (l'intégralité d'au moins une image 3D).

Afin de n'avoir à stocker que des images 2D, on applique l'algorithme défini par l'itération V en calculant la suite des estimées p(") =V[p(n-»] pour n =1,•••,N avec p(0' =0 ou p(0' = ps, c'est-à-dire tel que f=Rp`°'. L'itération V est définie pour tout vecteur p dans l'espace des images de projection 2D de la manière suivante : V[p]=(IùpRR)p+ps.

On note que si f =17p on a R(V[p])=R(I ù pRR)p+ pRs =(IùpRR) f+pRs . =A[f] On a ainsi montré par simple récurrence que si f o (°' l= Rp, es deux algorithmes, l'un défini dans le domaine des images de projection 2D par l'itération V et l'autre dans le domaine des images 3D par l'itération A, donnent la même succession d'images 3D f''') puisqu'on a montré que : =A"[f'°']=Rp'n' =RVn[p(°)]. Il est toutefois à remarquer que ceci n'est pas le cas si on a la contrainte p = Rf puisqu'on a alors : R(A[f 1) = R(I ù pRR) f + pRRs =(IùpR`R)p+pRRs . ≠V[p] Second algorithme Dans le cas particulier où R =R' , le premier algorithme ù ci-dessus décrit ù nous donne une solution à un problème d'optimisation quadratique.

Toutefois, on connaît que l'algorithme du gradient conjugué nous donne une solution plus efficace à ce problème d'optimisation, au sens où un nombre très réduit d'itérations doit être calculé pour obtenir une image visuellement acceptable. En revanche, du point de vue du stockage, le gradient conjugué requiert beaucoup plus de capacité de stockage que le premier algorithme présenté. On démontre ci-dessous que son implémentation peut être effectuée dans le domaine des images de projection 2D aussi bien que dans le domaine des images 3D.

De façon connue, dans le domaine des images 3D, l'algorithme du gradient conjugué peut s'implémenter en définissant l'itération A de la manière suivante r (n) , d (n) / = AL\f (n-1) r (n-1) d (n-1) /J avec a(n-1) _ ù \r (n-1) , ..7 (n-1) )/ \RtRd (n-1)' ..7 (n-1) = f (n-1) + a(n-1)d (n-1) = r(n-1) +a(n-1)R`Rd(n-1) (n-1) ù(r(n),r(n))/\r(n-1) d (n) r (n) + (n-1) r7 (n-1) où R est la matrice modélisant le système d'acquisition et g la matrice transposée, f('') est une image 3D estimée de la solution du problème à l'itération n, r(n) et d(n) sont des images 3D auxiliaires, et( , symbolise le produit scalaire de deux vecteurs de sorte que a et /3 sont deux réels. On initialise les trois vecteurs de l'espace des images 3D de la manière suivante : f (0) = R`p(0) r(0) = d(' = R` (Rf (0) ù s) où s est l'ensemble des images de projection 2D acquises, p(0) est arbitraire mais en pratique choisi égal à 0 ou s . On constate très simplement qu'une telle implémentation nécessite le stockage de quatre volumes à la fois : f ,,d(n),R`Rd(n). Comme précédemment, il est possible de s'affranchir du stockage de 20 ces images 3D par l'intermédiaire des changements de variables suivants : p(n) tel que f (n) = R`p(n) on fixe à cet effet q(n) tel que r (n) = R`q(n) tt(n) t tel que d(n)=R

où (p(n),q(n),t(n)) est un triplet de vecteurs de l'espace des images de projection 2D. On initialise le vecteur p(0) comme pour l'itération A, et q(0) =t( = RRtp(0) ù s L'itération V est définie de la manière suivante : )=VI(p(nù1>t(n-' )] avec a(n 1) = _(q(n 1) RRtt(n 1) )/( tt(n 1) RRtt(n 1) p(n) = p(n 1) + a(n 1)t(n 1) q(n) = q(n 1) + a(n 1)RRtt(n 1) /3(n 1) = (q(n), RRtq(n))/(q(n 1) RRtq(n 1)) t(n) = q(n) + /1(n 1)t(n 1) puis, on détermine f( =Rtp(N), où f est l'image 3D de l'objet à reconstruire. Avec cet algorithme on constate bien qu'à chaque itération on ne stocke que des images 2D ce qui permet d'obtenir une reconstruction moins coûteuse en besoin mémoire. Ainsi, il est possible de remplacer les étapes intermédiaires f (n) ,d(n),R`Rd(n) dans le domaine des images 3D par des étapes intermédiaires dans le domaine des images 2D : s(n) r(n) q(n) RRtq(n) On démontre ci-dessous l'équivalence entre l'itération A et l'itération V de l'algorithme du gradient conjugué. On montre d'abord que pour tout n =1,•••,N, les produits scalaires a et /3 sont les mêmes pour l'itération A dans le domaines des images 3D que pour l'itération V dans le domaine des images 2D puisque :

ù(r(" d(n))/(RtRd(n' d(n') ù(Rtq(n) Rtt(n))/(R`RRtt(n) Rttl ù(q(n),RRtt(n))/(RRtt(n),RRtt(n)) et r(n+1) r(n+1Mr(n) , = (Rtq(n+1) Rtq(n+1))/(Rtq(n), Rtq(n)) = (q(n+1) RRtq(n+1))/(q(n) RRtq(n)) Ceci conduit à obtenir, étant donné : f (n) = Rtp(n),r(n) = Rtq(n),d(n) = Rttn RRttn RRtq(n) les équivalences suivantes entre l'itération V du domaine des images 2D et 5 l'itération A du domaine des images 3D : f(n) = Rtp(n) = Rtq(n),d(n) = Rtt(n),RRtt(n), RRtq(n) q(n) RRtt(n) a(n) p(n+l) = p (n) + a(n)t(n) f (n+l) = f (n) +a(n)d(n) q(n+l) = q(n) +a n)RRtt(n) r(n+l) = r(n) + a(n)RtRd (n) RRtq(n+l) 3(n) t(n+l) = q(n+l) + p(n)t(n) .7(n+l) = r(n+l) + p(n).7(n) RRtt(n+l) = RRtq(n+l) + /1(n)RRtt(n) a(n+l)

on obtient donc p(0) .f (0) = Rtp(0) q(0) =t(0) =RRtp(0) ùs r(0) =d(0) =Rt(Rf(0) û s), détermination de RRtt(0) a(°) RRtq(0) = RRtt(0) p(1) = p(0) +a(0)t(0) f (1) = f (0) +a(0)d(0) q(1) = q(0) +a(0)RRtt(0) r(1) = r(0) +a(0)R`Rd(0) détermination de RRtq(1) 13(0) t(1) = q(l) + f3(0)t(0) d(1) = r(1) + /l(°)d(°) RRtt(1) = RRtq(1) +,(3(1)RRtt(1) Fa'(" On a démontré par récurrence l'équivalence entre l'itération V dans le 10 domaine des images 2D et l'itération A du domaine des images 3D. Enfin comme pour l'algorithme précédent, on obtient l'image 3D de l'objet en appliquant la contrainte linéaire, ici f(') = Rt p(N) . La formule du gradient conjugué donnée ci-dessus est celle de Fletcher et Reeves, mais un raisonnement similaire est possible pour d'autres formulations du gradient conjugué (Polak-Ribière, Graham-Schmidt,...). De même, le résultat s'applique pour une matrice R quelconque où dès que le produit RR est une matrice symétrique définie positive. 5 Troisième algorithme Ce troisième algorithme dit « itératif par blocs » généralise le premier dans le sens que l'on va scinder l'ensemble des indices des mesures en une partition de deux sous-ensembles d'indices. 10 On note O le premier sous-ensemble et cl) le second. Les deux sous-ensembles sont disjoints et leur réunion indexe l'ensemble de toutes les images 2D de projection. Dans une première variante, le sous-ensemble O contient un unique indice (soit un unique point de mesure du détecteur pour un seul angle 15 d'acquisition donné) alors que cl) indexe toutes les autres mesures (c'est-à-dire toutes sauf une). Cette partition est utilisée dans l'algorithme « Algebraic Reconstruction Technique (ART) ». Dans une seconde variante, le sous-ensemble O indexe toutes les images 2D pour un sous-ensemble d'angulations pour lesquelles les 20 projections sont acquises, le sous-ensemble cl) indexant les angulations acquises complémentaires. Cette partition est utilisée dans l'algorithme « Simultaneous Algebraic Reconstruction Technique (SART) ». On note que le cas où O contient tous les indices et cl) est vide correspond au cas du premier algorithme. 25 On note Io (respectivement Id,) la matrice dont les coefficients sont tous nuls sauf sur pour les indices de la diagonale appartenant à e (respectivement cl)) où ils valent 1. (O,c) étant une partition, la matrice identité dans l'espace des images 2D est 4 + I, . On note que le produit I®R (respectivement I~R) est la restriction de 30 R aux indices de O (respectivement cl)). On note de façon similaire que le produit RI, (respectivement RI,) est la restriction de k aux indices de O (respectivement cl)). A partir de ces définitions, l'algorithme dans le domaine des images 3D peut s'implémenter par l'itération Ao définie de la manière suivante : Ao [f] _ (I û p®RI®R) f + p®RIos où / est la matrice identité dans le domaine des images 3D et p® un paramètre fixé a priori pour que I û p®kI®R < 1. Bien évidemment, toute reconstruction fait intervenir toutes les mesures pour déterminer une solution, par conséquent la partition (0,c) d'indices n'est pas constante mais change à chaque itération. Il s'agit donc d'une pluralité d'itérations A®. L'algorithme génère ainsi une suite d'estimée f("' telles que: f (n) = --U(n) [ '"-"] U(n)--U(n-1) [f (n 2) ] = -l^I(n) O(n-1) ... A_ 1) [f(°)] On note qu'à chaque itération n, l'estimée est obtenue à partir d'un ensemble restreint de mesures 0(n), de sorte que le temps de calcul de l'itération est réduit. Il est connu que ces techniques fournissent donc des solutions au problème posé avec un temps de calcul réduit. En revanche, le problème associé au stockage des images 3D est inchangé et indépendant de la partition (0,c) choisie.

On démontre ici de façon nouvelle que l'on peut comme précédemment établir des itérations Vo dans le domaine des images de projections 2D qui conduiront aux mêmes estimées que celle obtenues avec les itérations Ao . Afin de n'avoir à stocker que des images 2D, on applique l'algorithme défini par l'itération Vo en calculant la suite des estimées p(n> =Vo n [p(n-1)] pour n =1,...,N avec p(0' tel que f (0' = R`p(0 et suivant les mêmes partitions O(n) que celles utilisées dans l'itération A, équivalente. L'itération Vo est définie pour tout vecteur p dans l'espace des images de projection 2D et pour toute partition (0,4)) des mesures telles que p=Iop+I,p, de la manière suivante : V0[pl = (Io ùp®I®RR)p+polos+I,p .

Enfin comme pour les deux autres algorithmes, on obtient l'image 3D de l'objet en déterminant f 'N' = Rp' . Ce troisième algorithme û comme on peut le constater û ne modifie à l'itération n que les variables indexées par O . On montre que, quelque soit la partition (0,4)), si f =17p , on peut passer de l'algorithme défini par les itérations k dans le domaine des images 3D à l'algorithme défini par les itérations Vo dans le domaine des images 2D de la manière suivante : R(17,[pl)=R((I® ùp®I®RR)p+polos)+Rl,p = R (Iop + I , p) ù p®RI®RRp + poRlos = f ù p®RI®Rf + poRlos =Ao[f Comme précédemment, l'itération Vo évite d'avoir à stocker à quelque moment que ce soit (initialisation n = 0, calcul intermédiaire 0 < n < N , archivage n= N) l'intégralité d'une image 3D dont la taille est très supérieure à l'ensemble des mesures acquises.

Application à la tomographie à faible nombre de projections Comme déjà mentionné, un problème en tomographie est le faible nombre de projections quand l'on ne peut pas tourner complètement autour de l'objet que l'on veut imager ou alors que l'on doit sous-échantillonner l'acquisition pour limiter le temps de l'examen ou la dose de rayons X au patient. En conséquence, la taille de l'image 3D que l'on veut reconstruire est beaucoup plus importante que la taille de l'ensemble des images de projection 2D acquises.

En utilisant les algorithmes de reconstruction tels que présentés, aucun volume n'est stocké afin que l'échantillonnage du volume puisse être réglé finement et de façon arbitraire sans aucun impact sur les besoins de stockage mémoire.

En effet, puisque l'étape finale f (N' =1p(N' , qui, pour qu'une seule coupe soit visualisée, est simplement un opérateur linéaire s'appliquant à d'un faible nombre d'images de projections 2D acquises traitées, il peut être effectué en temps réel lors de la visualisation de la coupe sans jamais avoir besoin de stocker l'ensemble des coupes du volume.

On note que selon la résolution qui est souhaité pour l'image 3D à reconstruire, le nombre d'itérations peut augmenter sans augmenter les besoins de stockage mémoire. Les procédés présentés sont particulièrement bien adaptés aux GPU (en anglais « Graphic Processor Units ») qui présentent de bonnes performances en ce qui concerne le temps de calcul mais qui sont limités en mémoire.

Claims (9)

1. REVENDICATIONS1. Procédé de traitement d'images obtenues par tomographie ou tomosynthèse comprenant : une étape d'acquisition (S0), au moyen d'un système d'imagerie, d'une pluralité d'images de projection 2D (si) d'un objet, l'acquisition étant définie par Rf = s où s est un vecteur des projections acquises, R est un opérateur de projection qui modélise le système d'imagerie et f est l'image 3D de l'objet à reconstruire connaissant R et s ; une étape de traitement (Si) des images de projection 2D acquises ; caractérisé en ce que le traitement des images de projection consiste en une application d'un algorithme itératif défini par son itération V qui génère à chaque itération n =1,...,N au moins un ensemble d'images de projection traitées p , l'itération V étant définie de manière à ce qu'à chaque itération l'image 3D de l'objet soit une fonction linéaire des images de projection 2D traitées selon la propriété fin' =Rpoù R est une matrice telle que RRf est une image 3D.
2. 2. Procédé selon la revendication 1 comprenant une étape de reconstruction 3D (S2) de l'objet à déterminer f(» = Rp(N) où N est le nombre d'itérations fixées. 25
3. 3. Procédé de traitement selon l'une des revendications 1 à 2, dans lequel à chaque itération n=1,...,N, l'algorithme V génère un ensemble d'images de projection traitées p(n' =V[p~n , l'itération V étant définie pour tout vecteur p dans l'espace des images de projection 2D de la20manière suivante V [p] _ (I ù pRR)p + ps où p(0) = ps ou bien p(0) = 0, s est l'ensemble des images de projection 2D acquises, p est un paramètre fixé pour que 11I ù p17111 <1 1 et ainsi garantisse la convergence du procédé et I est la matrice carrée identité.
4. 4. Procédé de traitement selon l'une des revendications 1 à 2 dans lequel à chaque itération, l'algorithme V génère les ensembles d'images de projection traitées (p n V [(p(n-" q(n-" t(n_1) )] suivant la méthode du gradient conjugué, typiquement donné par : a(n 1) = ù(q (n 1) RRtt(n v~~~RRtt(n v RRtt(n v~ p(n) = p(n 1) + a(n 1)t(n 1) q(n) = q(n 1) + a(n 1) Rtt(n 1) = (q RR q n)/(q(n 1),RRt ( t(n) = q(n) +!' (n-1)t(n-1) où p , q , t sont tels que p(°) est arbitraire, typiquement nul ou égal à s le vecteur des images de projection 2D acquises, et q(0 =t(0) =RRtp(0) ù s. 15
5. 5. Procédé de traitement selon l'une des revendications 1 à 2 dans lequel à chaque itération, l'algorithme V0(n) génère un ensemble d'images de projection p (n) = VO(n) (û1) ] = VO(n)VO(n-1) ' (n-2) ] = VO(n)VO(n-1) ... `%Oc1) [ (°) J avec V®(n) [P] = (I0(n) ù Po(n) I0(n)R11P + P®(n) I D(n) s + I,D(n) p où ®(n) et c(n) sont respectivement un premier et un second sous-ensemble des index des 20 mesures projectives acquises, définis pour chaque itération n, les deux sous-ensembles étant disjoints, leur réunion indexant la totalité du vecteur s des images 2D de projection acquises, et IO(n) et I,(n) sont des matrices dont les coefficients sont tous nuls sauf pour les indices de la diagonale appartenant 10à ®(n) et cl )(n) respectivement où ils valent 1, de sorte que Pù10(n)p+' (n)p, OÙ (°) P = P0(0)10(0) S ou bien p(°' =0, et p0(n) est un paramètre fixé pour que 1110(n) ù p0(n)RR <1 et ainsi garantisse la convergence du procédé .
6. 6. Procédé de traitement selon l'une des revendications précédentes comprenant une extraction d'une coupe issue de l'image 3D reconstruite.
7. 7. Procédé de traitement selon l'une des revendications précédentes 10 comprenant un réglage arbitrairement fin de l'échantillonnage de l'image 3D.
8. 8. Système d'imagerie médicale comprenant : ù un système d'acquisition comprenant une source de rayonnement, un capteur, pour l'acquisition d'une pluralité d'images de 15 projection 2D d'un objet ; ù des moyens de stockage des images de projection 2D acquises ; ù un système de traitement des images de projection 2D acquises destiné à mettre en oeuvre un procédé selon l'une des 20 revendications précédentes.
9. 9. Produit programme d'ordinateur comprenant des instructions de code de programme pour l'exécution des étapes du procédé selon l'une des revendications 1 à 7. 25