MASTER 1 Reconstruction tomographique Patrick Bourguet Rennes

1. Traitement et analyse d’images MASTER 1 Reconstruction tomographique Patrick Bourguet Rennes Support de cours adapté d’Irène Buvat buvat@imed.jussieu.fr octobre 2003 http://www.guillemet.org/irene

2. Plan du cours • A- Introduction • - Problème de reconstruction tomographique • - Tomographie en transmission • - Tomographie en émission • - Spécificité du problème de reconstruction tomographique • B- Notions de base • - Projection • - Transformée de Radon • - Sinogramme • - Rétroprojection • C- Méthodes de reconstruction analytique • - Principe • - Théorème de la coupe centrale • - Rétroprojection filtrée • - Filtres • D- Méthodes de reconstruction itérative • - Principe • - Opérateur de projection R • - Méthodes ART • - Méthodes SIRT • - Méthodes de descente • - Méthodes statistiques • - Caractéristiques des méthodes itératives • - Régularisation • - Interprétation probabiliste • E- Discussion • - Reconstruction analytique ou itérative ? • - Historique • - Quelle méthode pour quelle application ?

3. Découvrir l’intérieur… RX

4. X Tomographie axiale

5. A - Introduction : la reconstruction tomographique projections 2D sous différentes incidences angulaires Reconstruction tomographique sagittale transverse coronale coupes d’orientation quelconque : imagerie 3D Reconstruction tomographique = problème inverse : estimer la distribution 3D à partir des projections 2D mesurées

6. X d d   0 0 Introduction : tomographie de transmission • Source (X ou g) externe au patient • Mesures : projection du rayonnement ayant traversé le patient  intégrale des coefficients d’atténuation • Objet à reconstruire : cartographie 3D des coefficients d’atténuation m du milieu traversé N0 N0 N0 g b+ N N N scanner X (tomodensitométrie) SPECT PET N0 ln (l) dl N = N0 exp(- (l) dl) N

7. d  0 Introduction : tomographie d’émission • Source g ou b+ interne au patient • Mesures idéales (sans atténuation) : intégrale de l’activité le long des raies de projections • Objet à reconstruire : cartographie 3D de la distribution d’activité n dans l’organisme détecteur * détecteur détecteur * * * * détecteur PET SPECT N = n(l) dl y z x

8. Factorisation du problème de reconstruction • Un ensemble de projections 2D • reconstruction d’un objet 3D • Un ensemble de projections 1D • reconstruction d’un objet 2D (coupe zi) • ensemble de coupes zi = volume d’intérêt détecteur en position q une projection q z x N(x,z) volume 3D étudié y z x une ligne de projection détecteur en position q z N(x) x x coupe axiale zi

9. Reconstruction : non unicité de la solution • Pas de solution unique : toujours plusieurs objets compatibles avec un ensemble fini de projections • Unicité de la solution pour une infinité de projections seulement 1 projection : plusieurs solutions possibles direction de projection 2 projections : plusieurs solutions possibles directions de projection

10. Reconstruction : problème inverse mal posé • Pas de solution du fait du bruit entachant les données • Non unicité de la solution du fait du manque d’informations (nombre de projections fini) • typiquement, 64 ou 128 projections <<  • Instabilité de la solution : petite différence sur les projections peut provoquer des écarts importants sur les coupes reconstruites projection idéale projection bruitée

11. B - Notions de base !

12. Opération de projection y u = x cos + y sin v = -x sin + y cos f(x,y) v u  x   p(u,) = f(x,y) dv - u projection 

13.   p(u,) = f(x,y) dv - Transformée de Radon • ensemble des projections pour q = [0, p ] • = transformée de Radon de f(x,y) • R[f(x,y)] = p(u,q) dq p  0 f(x,y) p(u,q) domaine spatial espace de Radon Problème de reconstruction tomographique : inverser la transformée de Radon, i.e., estimer f(x,y) à partir de p(u,q)

14. p  0 C - Méthodes de reconstruction analytique • Inversion analytique de la transformée de Radon • Expression continue du problème de reconstruction • tomographique • Méthode la plus courante : rétroprojection filtrée • FBP : Filtered Back Projection • Méthodes rapides • Méthodes disponibles sur tous les dispositifs d’acquisition commercialisés (scanner X, SPECT, PET) • R[f(x,y)] = p(u,q) dq

15. Opération de rétro-projection « simple » y u = x cos + y sin v = -x sin + y cos f(x,y) v u  x p  f*(x,y) = p(u,) d 0 u rétroprojection Attention : la rétroprojection ‘’simple’’ n’est pas l’inverse de la transformée de Radon

16. p  0 Limites de la rétro-projection simple f*(x,y) = p(u,) d u rétroprojection simple artefacts d’épandage en étoile nombre de projections 3 1 4 image originale 16 32 64 La rétroprojection n’inverse pas la transformée de Radon

17. p  0 p  0 Rétro-projection filtrée : principe f*(x,y) = p(u,) d u rétroprojection simple f*(x,y) = p(u,) d projection filtrée u rétroprojection filtrée réduction des artefacts inversion de la transformée de Radon

18.           - - - - - changement de variable : (u,v) (x,y) Théorème de la coupe centrale (TCC) transformée de Fourier (TF)   P(r,) = p(u, )e-i2pru du p(u,) = f(x,y) dv - y u = x cos + y sin v = -x sin + y cos v u x = cos  y =sin  du.dv = dx.dy  x P(r,) = f(x,y)e-i2pru du dv = f(x,y)e-i2p(xrx + yry) dx dy TF monodimensionnelle d’une projection par rapport à u = TF bidimensionnelle de la distribution à reconstruire

19.             - - - - - - p  0 p  0 Rétro-projection filtrée P(,) = F(x,y) TF-1 f(x,y) = F(rx,ry)ei2p (xrx + yry) drx dry = P(r,q)ei2p(xrx + yry) drx dry = P(r,q) |r|ei2pru dr dq = p’(u,q) dq avec p’(u,q) = P(r,q) |r|ei2pru dr TCC x = cos  y =sin   = (x2 + y2)1/2 dx.dy = .d.d u =x cos + y sin changement de variable : (rx, ry) (r, q) projections filtrées filtre rampe objet f(x,y) à reconstruire = rétroprojection des projections filtrées

20.   - f(x,y) = p’(u,q) dq avec p’(u,q) = P(r,q) |r|ei2pru dr p  0 Algorithme de rétro-projection filtrée projections images reconstruites p(u,) f(x,y) transformée de Fourier rétroprojection Transformée de Fourier inverse filtrage P(r,) |r| P(r,) p’(u,)

21.   - f(x,y) = p’(u,q) dq avec p’(u,q) = P(r,q) |r|ei2pru dr 1 1 1 p  0 0 0 0 0 0,8 0 0,8 0 0,8 Insuffisance du filtre rampe filtre rampe |r| |r|w(r) fenêtre d’apodisation Exemple : |r| w(r) |r|w(r) fenêtre d’apodisation de Hann filtre de Hann résultant filtre rampe w() = 0,5.(1+cos/c) si< c = 0si ≥ c domaine fréquentiel 1 2 1 1/16 2 4 2 si c= N 1 2 1 domaine spatial

22. Filtres classiques : filtre de Hann • Filtre rampe • meilleure résolution spatiale mais forte amplification du bruit haute fréquence • Filtre de Hann • modifie les moyennes fréquences w() = 0,5.(1+cos/c) si< c = 0si ≥ c 0,5 0,4 0,3 0,2 0,1 fréquence de coupure c plus faible est la fréquence de coupure, moins on préserve les détails “haute fréquence”, i.e., plus fort est le lissage

23. Filtres classiques : filtre de Hamming • Filtre rampe • Filtre de Hamming w() = 0,54+0,46(cos/c) si< c = 0si ≥ c 0,5 0,4 0,3 0,2 0,1 fréquence de coupure c plus faible est la fréquence de coupure, moins on préserve les détails “haute fréquence”, i.e., plus fort est le lissage

24. Filtres classiques : filtre gaussien • Filtre rampe • Filtre gaussien (domaine spatial) c(x) = (1/s 2p).exp[-(x-x0)2/2s2] si< c = 0si ≥ c 0 1 2 3 4 FWHM = 2 2ln2 s (pixel) plus grande est la dispersion du filtre gaussien (FWHM ou s), moins on préserve les détails “haute fréquence”, i.e., plus fort est le lissage

25. Filtres classiques : filtre de Butterworth • Filtre rampe • Filtre de Butterworth w() = 1/[1+(/c)2n] si< c = 0si ≥ c 2 paramètres : fréquence de coupure c et ordre n 2 4 6 8 10 ordre n, c=0,25 plus faible est l’ordre moins on préserve les détails “haute fréquence”, i.e., plus fort est le lissage

26. |r| |r|w(r) Implantation du filtrage • Multiplication du filtre rampe par une fenêtre d’apodisation • filtrage 1D (direction x) • Filtrage des projections puis reconstruction avec un filtre rampe • filtrage 2D (directions x,z) • Reconstruction avec un filtre rampe puis filtrage 2D des coupes reconstruites • filtrage 2D (directions x,y) • Reconstruction avec un filtre rampe puis filtrage 3D du volume reconstruit • filtrage 3D (directions x,y,z) z x coupe zi y

27. p1 p2 p3 p4 r11 r14 r41 r44 f1 f2 f3 f4 = D - Méthodes de reconstruction itératives • Expression discrète et matricielle du problème de reconstruction tomographique • Problème de reconstruction : système d’équations de grande taille • 128 projections 128 x 128 • 2 097 152 équations à autant d’inconnues • Inversion itérative du système d’équations

28. p1 p2 p1 p2 p3 p4 r11 r14 r41 r44 f1 f2 f3 f4 = p3 f1 f2 p4 fk f3 f4 pi projection Expression discrète du problème de reconstruction p1 = r11 f1 + r12f2 + r13f3 + r14 f4 p2 = r21 f1 + r22f2 + r23f3 + r24 f4 p3 = r31 f1 + r32f2 + r33f3 + r34 f4 p4 = r41 f1 + r42f2 + r43f3 + r44 f 4 p = R f projections acquises opérateur de projection objet à reconstruire Problème : déterminer f connaissant p et R

29. Expression de l’opérateur de projection R p = R f fk pi projection • Modélisation géométrique • - choix d’un modèle de distribution de l’intensité des pixels • - modélisation de la géométrie de détection • Modélisation de la physique de détection • * atténuation • * diffusion • * résolution limitée du détecteur

30. Modélisation géométrique de l’opérateur R • Modèle de distribution de l’intensité des pixels : détermination de la contribution de chaque pixel i à une raie de projection k • Modèle de la géométrie de détection bins de projection 0 1 0 0 modèle exact = modèle uniforme modèle de Dirac l modèle de longueur de raies modèle du disque concave bins de projection parallèle en éventail

31. p1 p2 f1 f2 p1 p2 f3 f4 p3 p4 f1 f2 f3 f4 Modélisation physique de l’opérateur R • Atténuation • Diffusion • Réponse du détecteur p1 = r11 f1 exp(-1d1) + r13f3 exp(-3d3) carte des m sans modélisation de la diffusion : p1 = r11 f1 + r13 f3 avec modélisation de la diffusion : p1 = r11 f1 + r12 f2 + r13 f3 + r14 f4 p1 p2 sans modélisation de la fonction de réponse de la caméra : p1 = r11 f1 + r13 f3 avec modélisation : p1 = r11 f1 + r12 f2 + r13 f3 + r14 f4 f1 f2 f3 f4

32. fk pi projection fk pi rétroprojection Opérateurs de projection R et de rétroprojection p1 p2 f1 f2 p3 f3 f4 p4 p1= f1 + f3 p2= f2 + f4 p3= f1 + f2 p4= f3 + f4 f1= p1 + p3 f2= p2 + p3 f3= p1 + p4 f4= p2 + p4 1 0 1 0 0 1 0 1 1 1 0 0 0 0 1 1 1 0 1 0 0 1 1 0 1 0 0 1 0 1 0 1 R Rt = =

33. Résolution du problème inverse p = R f Recherche d’une solution f minimisant une distance d(p,Rf), p et R étant connus estimée initiale de l’objet à reconstruire fn=0 p = R f facteur de correction cn comparaison ^ ^ pn pn vs p fn+1 projection correspondant à l ’estimée fn définit la méthode itérative : additive : fn+1 = fn + cn multiplicative : fn+1 = fn . cn

34. 7/8 1/8 1/8 7/8 R = Exemple • Deux pixels et deux raies de projection • système de deux équations à deux inconnues p1 p2 f1 f2 f1 = 10 f2 = 15 Opérateur de projection (connu) Valeurs de projection mesurées p1 = 7/8f1 + 1/8f2 = 85/8 p2 = 1/8f1 + 7/8 f2 = 115/8 Inconnues à déterminer: f1 et f2

35. Approche ART (Algebric Reconstruction Technique) p1 = 7/8f1 + 1/8f2 = 85/8 (équation 1) p2 = 1/8f1 + 7/8 f2 = 115/8 (équation 2) • Utilisation d’une seule équation de raie par itération, et mise à jour de chaque inconnue à partir de cette équation : • itération 1 : estimation de f1 à partir de l’équation 1 slmt • estimation de f2 à partir de l’équation 1 slmt • itération 2 : estimation de f1 à partir de l’équation 2 slmt • estimation de f2 à partir de l’équation 2 slmt • itération 3 : estimation de f1 à partir de l’équation 1 slmt • estimation de f2 à partir de l’équation 1 slmt • etc… • Modification à chaque itération proportionnelle à l’erreur par rapport à la projection vraie • ART additive (erreur = pk - pkn) • ou ART multiplicative (erreur = pk / pkn)

36. ART additive p1 = 7/8f1 + 1/8f2 = 85/8 = 10,6 p2 = 1/8f1 + 7/8 f2 = 115/8 = 14,4 • Repose sur la méthode de Kaczmarz • fin+1 = fin + (pk - pkn ) rki/Sjrkj2 • Exemple (initialisation f10 = f20 = 0  p10 = p20 = 0 ) : • itération n=1, raie k=1, estimation de f1 et f2 : • f11 = 0 + (85/8-0).(7/8)/(49/64+1/64) = 119/10 = 11,9 • f21 = 0 + (85/8-0).(1/8)/(49/64+1/64) = 17/10 = 1,7 •  p21 = 119/40 = 3,0 • itération n=2, raie k=2, estimation de f1 et f2 : • f12 = 11,9 + (115/8-119/40).(1/8)/(49/64+1/64) = 13,7 • f22 = 1,7 + (115/8-119/40).(7/8)/(49/64+1/64) = 14,5 •  p11 = 13,8 • itération 1 2 3 4 5 • f1 11,9 13,7 10,1 10,3 10,0 • f2 1,7 14,5 13,9 14,9 14,9 • A chaque itération, mise à jour successive de toutes les inconnues (ici f1 et f2) à partir de l’équation correspondant à une seule raie de projection contribution du pixel i à la raie de projection k écart entre valeurs du bin de projection estimée et observée

37. Approche SIRT (Simult. Iterative Reconst° Tech.) p1 = 7/8f1 + 1/8f2 = 85/8 (équation 1) p2 = 1/8f1 + 7/8 f2 = 115/8 (équation 2) • SIRT : Simultaneous Iterative Reconstruction Technique • Utilisation de toutes les équations à chaque itération et mise à jour de chaque inconnue : • itération 1 : estimation de f1 en résolvant l’équation 1 estimation de f2 en résolvant l’équation 2 • itération 2 : estimation de f1 en résolvant l’équation 1 estimation de f2 en résolvant l’équation 2 • itération 3 : estimation de f1 en résolvant l’équation 1 estimation de f2 en résolvant l’équation 2 • etc… • Modification à chaque itération proportionnelle à l’erreur par rapport à la projection vraie • Méthode de Jacobi • Méthode de Gauss-Seidel • Nombre d’itérations réduit par rapport aux méthodes ART

38. SIRT : méthode de Jacobi p1 = 7/8f1 + 1/8f2 = 85/8 = 10,6 p2 = 1/8f1 + 7/8 f2 = 115/8 = 14,4 contribution du pixel i à la raie de projection i • fin+1 = fin + (pk - pkn ) / rii • Exemple (initialisation f10 = f20 = 0  p10 = p20 = 0 ) : • itération n=1 : • estimation de f1 en résolvant l’équation 1 • f11 = (85/8)/(7/8) = 85/7 = 12,1 • estimation de f2 en résolvant l’équation 2 • f21 = (115/8)/(7/8) = 115/7 = 16,4 • itération n=2 : • estimation de f1 en résolvant l’équation 1 • f12 = [85/8-(1/8)*(115/7)]/(7/8) = 9,8 • estimation de f2 en résolvant l’équation 2 • f22 = [115/8-(1/8)*(85/7)]/(7/8) = 14,7 • itération 1 2 3 • f1 12,1 9,8 10,0 • f2 16,4 14,7 15,0 écart entre valeurs du bin de projection estimée et observée

39. SIRT : méthode de Gauss-Seidel p1 = 7/8f1 + 1/8f2 = 85/8 = 10,6 p2 = 1/8f1 + 7/8 f2 = 115/8 = 14,4 • Identique à la méthode de Jacobi mais en tirant immédiatement parti des estimations obtenues aux itérations précédentes • fin+1 = fin + (pk - pkn ou (n+1) ) / rii • Exemple (initialisation f10 = f20 = 0  p10 = p20 = 0 ) : • itération n=1 : • estimation de f1 en résolvant l’équation 1 • f11 = (85/8-0)/(7/8) = 85/7 = 12,1 • estimation de f2 en résolvant l’équation 2 avec f11=85/7 • f21 = [115/8-(1/8)*(85/7)]/(7/8) = 14,7 • itération n=2 : • estimation de f1 en résolvant l’équation 1 avec f21=14,7 • f12 = (85/8-14,7/8)/(7/8) = 10,0 • estimation de f2 en résolvant l’équation 2 avec f12=10,0 • f22 = (115/8-10,0/8)/(7/8) = 15,0 • itération 1 2 • f1 12,1 10,0 • f2 14,7 15,0 estimé à partir de toutes les valeurs courantes des fi  convergence plus rapide que Jacobi

40. Méthodes de descente • Méthodes ART et SIRT : • fin+1 = fin + a (erreurn) • Méthodes de descente : • fin+1 = fin + an(erreurn) • modification du coefficient de pondération à chaque itération • Méthode du gradient : Iterative Least Squared Technique (ILST) • Méthode du gradient conjugué : meilleure méthode de descente coefficient de pondération constant e.g., pk - pkn coefficient de pondération optimisé

41. Méthodes de descente : gradient conjugué • Adaptée à la résolution d’un système d’équations dont la matrice est symétrique : • p = R f  Rtp = RtR f • Formule de mise à jour : • fin+1 = fin + andn • Vitesse de descente optimisée à chaque itération : • an= ||Rtp - RtR fn||2 / (Rtp - RtR fn)t R (Rtp - RtR fn) • Direction de descente optimisée à chaque itération : • itération 1 : d1 = Rtp - RtR f0 • itération n : dn = (Rtp - RtR fn) + bn dn-1 • optimisation de la convergence • convergence rapide • méthode additive • utilisée en SPECT matrice symétrique coefficient de pondération optimisé à chaque itération (vitesse de descente) direction de descente optimisée à chaque itération

42. Méthode statistique : MLEM • MLEM = Max. Likelihood Expectation Maximization • Utilise une formulation probabiliste du problème de reconstruction : • - modèle probabiliste : • Les pk sont des variables aléatoires de Poisson de paramètres pk, d’où l’expression de la vraisemblance de f : • prob(p|f) = P exp(- pk) . pk / pk ! • - détermine la solution f qui maximise la vraisemblance (ou log-vraisemblance), i.e., prob(p|f) par rapport au modèle probabiliste choisi. pk k

43. Algorithme MLEM • Deux étapes : • - calcul de l’espérance de la log-vraisemblance compte tenu des projections pk mesurées et de l’estimation courante des fi. • - maximisation de l’espérance en annulant les dérivées partielles par rapport à fi. • Formule de mise à jour : • fin+1 = fin . [S rki (pk/ pkn)] / S rki • fn+1 = fn . Rt[ p/ pn ] • * méthode multiplicative • * solution toujours positive ou nulle • * nombre d’événements conservé au fil des itérations • * convergence lente • * méthode itérative la plus utilisée en SPECT • (dans sa version accélérée OSEM) k k opérateur de rétroprojection erreur

44. Version accélérée de MLEM : OSEM • OSEM = Ordered Subset Expectation Maximisation • Tri des P projections en sous-ensembles ordonnés • Exemple : • Application de MLEM sur les sous-ensembles : • - itération 1 : • estimation de f1à partir de l’initialisation f0et • des projections p1 correspondant au sous-ensemble 1 • f1 = f0 . Rt[ p/ p1 ] • estimation de f’1à partir de f1et des projections p’1 • correspondant au sous-ensemble 2 • f’1 = f1 . Rt[ p/ p’1 ] • - itération 2 : • estimation de f2à partir de f’1et des projections p2 • correspondant au sous-ensemble 1 • f2 = f’1 . Rt[ p/ p2 ] • estimation de f’2à partir de f2et des projections p’2 • correspondant au sous-ensemble 2 • f’2 = f2 . Rt[ p/ p’2 ] • etc. 8 projections 2 sous-ensembles de 4 projections

45. Caractéristiques de OSEM • OSEM avec S sous-ensembles et I iterations •  SI itérations de MLEM • Facteur d’accélération ~ nombre de sous-ensembles • méthode itérative la plus utilisée en SPECT MLEM 1 16 24 32 40 itér. OSEM 1 4 6 8 10 itér. 4 ss-ens. OSEM 1 2 3 4 5 8 ss-ens.

46. Caractéristiques des méthodes itératives • Plus élevé est le nombre d’itérations, meilleure est la restitution des hautes fréquences • Problème du choix du nombre d’itérations • - convergence vers la solution puis divergence de • la procédure lors de la reconstruction des très • hautes fréquences du fait de la présence de bruit • (haute fréquence) • nécessité de régulariser OSEM 1 2 3 4 5 8 ss-ens. gradient 1 2 3 4 5 conjugué

47. Importance de la régularisation • Régularisation implicite : arrêt précoce des itérations • Régularisation explicite : introduction d’un a priori sur la solution : • - solution non régularisée : • minimisation de d(p,Rf) • - solution régularisée : • minimisation de d1(p,Rf) + ld2(f,fa) • solution compromis entre la fidélité aux mesures et l’a priori • Exemples d’a priori régularisants : • - distribution de f connue (Poisson, Gauss) • - image lisse • - image présentant des discontinuités a priori régularisant

48. Effet de la régularisation sans régularisation image idéale avec régularisation itérations : 0 30 80 sans régularisation avec régularisation

49. Interprétation probabiliste • Problème de reconstruction : • Chercher f la plus probable compte tenu des projections p observées • Interprétation probabiliste : • Maximiser prob(f|p) : probabilité avoir l’image f quand les projections valent p • Théorème de Bayes : • prob(f|p) = prob(p|f) . prob(f) / prob(p) • projections connues  prob(p) = 1 • pas d’hypothèse a priori sur f  prob(f) = 1 • alors : • prob(f|p) = prob(p|f) • maximiser prob(f|p) = maximiser la vraisemblance prob(p|f) • = minimiser l’écart entre projections calculées et observées probabilité a priori de l’image f probabilité a priori des projections p probabilité de mesurer les projections p pour une image f = vraisemblance des projections p

50. Algorithmes associés à l’interprétation probabiliste • maximiser prob(f|p) = maximiser la vraisemblance prob(p|f) • = minimiser l ’écart entre projections calculées et observées • pk ~ loi de Poisson • algorithme MLEM • pk ~ loi de Gauss • algorithme gradient conjugué • Régularisation • prob(f|p) = prob(p|f) . prob(f) • méthodes MAP (maximum a posteriori) • versions régularisées : • MLEM  MAP-EM • Gradient Conjugué (GC)  MAP-GC probabilité a priori de l’image f ≠ 1 probabilité a posteriori de l’image f