Méthodes standards de reconstruction tomographique en SPECT/PET

Méthodes standards de reconstruction tomographique en SPECT/PET Philippe Ciuciu(CEA/SHFJ) ciuciu@shfj.cea.fr http://www.madic.org/people/ciuciu

Cours préparé à partir des cours de Master de physique médicale, Univ. Paris Sud (Orsay) d‘Irène Buvat (CNRS, INSERM U678) Et de Claude Comtat (CEA/SHFJ)

Plan • Le problème de la reconstruction tomograhique • Méthodes de reconstruction analytique • Rétroprojection filtrée • Techniques de reconstruction itératives • Méthodes « algébriques » • Méthodes statistiques ou pénalisées

Problématique : images détectées par la gamma caméra Intégrale du rayonnementgémis dans différentes directions

Problématique : signaux détectés par le tomographe à EP plans droits plans droits plans croisés Lignes de réponse dans toutes les directions

Problématique : estimer la distribution 3D du radiotraceur coupe transaxiale (ou transverse) coupe sagittale coupe coronale … à partir de mesures intégrales de cette distribution dans différentes directions

1 projection 2D détecteur en positionq 1 projection 1D détecteur en positionq coupe axiale reconstruction d’un volume 3D Reconstruction tomographique : factorisation du problème Ensemble de projections 2D Ensemble de projections 1D volume 3D étudié reconstruction d’une coupe 2D y z juxtaposition des coupes x volume 3D

3 1 4 16 32 64 Principe de la reconstruction tomographique   r projection rétroprojection

 rétroprojection filtrée r  projection filtrée Principe de la reconstruction tomographique

y v  x Opérateurs impliqués en reconstruction tomographique u = x cos  + y sin  projection rétroprojection     p(u,) = f(x,y) dv f *(x,y) = p(u,) d - 0

+ +   p(u,) = f(x,y) dv P(,) = p(u,) e-i2u du - - Théorème de la tranche centrale transformée de Fourier 1D y v u = x cos  + y sin  x =  cos  y =  sin  du.dv = dx.dy p(u,)  x + + + +     P(,) = f(x,y) e-i2u du.dv = f(x,y) e-i2(xx+yy) dx.dy - - - -

Rétroprojection filtrée P(,) = F(x, y) + +   x =  cos  y =  sin   = x2 + y2 d x.d y = d.d u = x cos + y sin f(x,y) = F(x, y) ei2(xx+yy) d x. d y - - + +   = P(, ) ei2(xx+yy) d x. d y - -  +   = P(, ) |  |ei2u d . d  - 0  +   = p’(u, ) d avec p’ (u, ) = P(, ) |  |ei2u d  - 0

  f(x,y) = p’(u, ) d avec p’ (u, ) = P(, ) || ei2u d  - 0 Algorithme de rétroprojection filtrée  + projections images reconstruites p(u,) f(x,y) TF rétroprojection filtrage TF-1 P(,) || P(,) p’(u,)

Algorithme de rétroprojection filtrée • Filtered Back-Projection FBP : • 1) calculer la transformée de Fourier 1D d’une projection pour un angle  fixé • 2) multiplier par le filtre rampe || • 3) calculer la transformée de Fourier inverse 1D de la projection filtrée • 4) rétroprojeter la projection filtrée • 5) répéter les étapes 1 à 4 pour chaque angle 

Algorithme de rétroprojection filtrée p(xr,f) Convolution par le filtre rampe h(xr) = F-1{|vxr|} p(xr,f)*h(xr) xr xr Opération de projection effectuée par le scanner Rétroprojection

  f(x,y) = p’(u, ) d avec p’ (u, ) = P(, ) || ei2u d - 0 ||w Insuffisance du filtre rampe  + || w ||w 1 1 1    0 0 0 0 0,8 0 0,8 0 0,8 filtre rampe x filtre de Hann filtre résultant

4 4 8 8 1 1 2 2 Sans filtrage,#= 144 36 36 16 16 72 72 144 144 Avec filtrage,# = 144 Effet du Filtrage

SPECT cérébral HMPAO Tc-99m Syndrome de fatigue chronique Avant traitement Après traitement IRM anatomique

Reconstruction analytique vs algébrique rétroprojection filtrée reconstruction algébrique

reconstructions analytiques reconstructions itératives Modèle : Intégrale ligne itérative Améliorer le Compromis Résolution  Bruit • Modèles plus complexes : • bruit statistique • dispositif d’acquisition analytique

Reconstructions Itératives en TEP Modèle « réaliste » Trop complexe pour une inversion analytique directe Approcher la solution pas à pas, en partant d’une image initiale

Reconstructions Itératives • Plutôt que d’estimer la distribution 3D par une implémentation discrète d’une solution analytique telle que FBP (opération directe), on l’estime par une succession d’affinages : • meilleure modélisation du dispositif discret d’acquisition des données que le modèle de l’intégrale ligne • incorporation d’un modèle statistique de bruit dans les données • Intégration durant le processus de reconstruction d’informations a priori connuessur l’image.

Reconstructions Algébriques et Statistiques • On distingue deux classes d’approches itératives: • 1) algébriques [1] • incluent un modèle discret du dispositif d’acquisition • 2) statistiques • incluent un modèle discret du dispositif d’acquisition • incluent un modèle statistique du bruit • peuvent inclure un a priori (bayésiennes)

Reconstructions Itératives : Caractéristiques • Une reconstruction itérative est caractérisée par [Fessler98]: • Un modèle paramétrique de l’image, l = {lj | j = 1,...,n} ; • un modèle des mesures, reliant les données discrètes mesuréesy = {yi | i = 1,...,m}à l’imagel: • un modèle du bruit, c’est-à-dire une distribution de probabilité poury ; • Un critère à minimiser et un algorithme itératif de minimisation

Paramétrisation non unique de l’image Paramétrisation de l’Image • Voxel • Blobs:fonctions à symétrie sphérique (fonctions généralisées de Kaiser-Bessel) qui se chevauchent [Lewitt90,Lewitt92]

Modèle d’acquisition : Matrice système  approche matricielle si(x) probabilité de détecter dans la LOR i une désintégration survenant enx. Aijprobabilité de détecter dans la LOR i une désintégration survenant dans le voxelj Apeut être complexe à souhait…

paire de détecteurs i Pixel j Matrice Système géométrique • i) projection géométrique • Aij volume d’intersection entre le voxel j et la LOR i

paire de détecteurs i paire de détecteurs i Pixel j Pixel j’ aij aij’ Matrice Système géométrique (suite) • ii) angle solide • Aij angle solide sous lequel la paire de détecteurs i est vue depuis le voxel j aij’ > aij

Pixel j Pixel j Modéliser dans la matrice système • profondeur de pénétration dans le cristal • diffusions dans les cristaux • non-colinéarité de la paire de photons • LOR non mesurées (gaps) • …

Factorisation de la matrice système • Matrice système peut être gigantesque (tera-octet) • 1) calculer les éléments de la matrice en ligne (routine clinique, généralement uniquement géométrique) • 2) techniques de factorisation [Qi98] : • A = Asens. dét. Arép. dét.  Aattn  Agéom. • Asens. dét. : sensibilité des détecteurs (diagonale) • Arép. dét. : fonction de réponse des détecteurs • Aattn : atténuation (diagonale) • Agéom. : projection géométrique

Modèle du Bruit (TEP) • Hypothèse de Poisson pour l’acquisition si : • le nombre d’atomes du traceur injecté suit une distribution de Poisson • les localisations spatiales des atomes du traceur, en tout instant, sont des variables aléatoires indépendantes • les instants auxquels surviennent la désintégration des atomes du traceur sont des variables aléatoires indépendantes, qui suivent une loi exponentielle de moyenne= t1/2/ln2 • les processus de détection de chaque désintégration sont des processus aléatoires indépendants (pas de temps mort)

Modèle du Bruit (suite) • Données pré-traitées : • ré-échantillonage, correction de l’efficacité de détection, de l’atténuation, des coïncidences fortuites et diffusées, du temps mort,… •  Données pas Poisson • Données Poisson : corrections dans la matrice système • Pour satisfaire un modèle de Poisson, il faut inclure TOUTES les corrections dans la matrice système et ne reconstruire que des données brutes. Souvent, non réalisable.

Choix d’un modèle statistique de bruit • Aucun :y = A·l • Bruit blanc gaussien : Avec bruit :y = A·l + b • Bruit gaussien corrélé/coloré • Bruit poissonnien

Distribution de Poisson : • variable entière positive • un seul paramètre (moyenne = variance) • Distribution de Gauss : • variable réelle • deux paramètres (moyenne et variance) Distributions de Poisson et de Gauss

Fonction de Coût/Critère (y,Al) terme d’attache de l’image laux données de projection mesurées y, ie vraisemblance des observations U(l) terme d’attache delà un modèle a priori de l’imageie terme de régularisation Contrainte (non négativité, …)

Attache aux données :vraisemblance Mesure la distance entre les mesures et le modèle • Bruit blanc gaussien : terme de moindres-carrés • Bruit gaussien corrélé : moindres-carrés pondérés (WLS) [Fessler94] • Bruit poissonnien : distance de Kullback-Leibler [Titterington87]

Maximum de vraisemblance (ML), sans a priori Maximum de vraisemblance Approche statistique : terme d’attache aux données  maximum de vraisemblance Modèle de bruit poissonnien :

Maximum a posteriori (MAP) [Hebert89] Maximum a Posteriori Approche statistique : modèle a priori règle de Bayes Typiquement :

Exemple d’A Priori : Champ de Markov Champ de Markov : modèle de corrélations spatiales locales [Geman84] N3D

Incorporation Information Anatomique • Modèle simple : • désactiver (wjk = 0) les corrélations locales si deux voxels j et k appartiennent à deux structures anatomiques différentes Image anatomique obtenue par : IRM ou Rayons X (CT)

Maximisation itérative du critère • Une fois définis : • un modèle de l’image (voxels, blobs) • un modèle de l’acquisition (matrice A) • un modèle statistique du bruit (Poisson, Gauss) • un critère (Moindres-Carrés, ML, MAP) • Il nous reste à choisir un algorithme permettant de le maximiser • Algorithme le plus classique : • EM-ML et variations (OSEM, RAMLA, …)

Algorithme Expectation Maximization (EM) [Dempster77] Critère : fonction de vraisemblance (ML) Algorithme EM : recherche du maximum de vraisemblance A) Absence de maxima locaux 2L : semi-définie négative car  n,nT ·2L·n  0 L(l) est concave : maximisation locale converge vers un maximum global

EM-ML (2) • Espace incomplety: • pas de maximisation analytique carp(yi|l) dépend de {lj | j = 1,...,n} • Espace complet x (y+données manquantes) : • contribution du voxel j à la mesure i{xij | j = 1,...,n; i = 1,...,m; yi = jxj} • p(xij|lj) ne dépend que delj

EM-ML (3) Maximisation analytique désormais possible Mais x inconnu !!! • Principe de l’algorithme EM [Dempster77] • Plutôt que de maximiser log p(x |l), on maximise l’espérance sur x de log p(x |l), en utilisant une estimation l(p)del pour estimer x

Deux étapes dans EM : • 1) calcul de l’espérance (E-step) • 2) maximisation de l’espérance (M-step) • On montre queL(l(p+1))  L(l(p)) EM-ML (4)

EM-ML (5) E-step : M-step :

Méthodes standards de reconstruction tomographique en SPECT/PET