Rapport de stage

6 mai - 17 janvier 2003

 

Effectué par: Eugène Lakinsky

Lieu du stage: Centre d'Imagerie Métabolique et Fonctionnelle, Département de Médecine Nucléaire et Radiobiologie, Université de Sherbrooke

 

Introduction

J'ai effectué mon stage dans le Centre d'Imagerie Métabolique et Fonctionnelle du Département de Médecine nucléaire et Radiobiologie (Université de Sherbrooke), situé au Centre de Recherche Clinique du Centre Hospitalier Universitaire de Sherbrooke.

L'équipe de recherche en Médecine nucléaire et Radiobiologie, qui comprend un groupe des Instituts de recherche en santé du Canada en sciences des radiations, privilégie une approche multidisciplinaire innovatrice permettant l'étude de problèmes reliés aux radiations sous plusieurs angles. Le programme de recherche vise la description de l'action des rayonnements dans le milieu biologique et son utilisation pour le diagnostic et le traitement de maladies. Le programme compte également un volet préventif permettant l'évaluation des risques et dangers résultant de la radioactivité. Sur le plan des applications, l'équipe développe de nouvelles approches d'imagerie médicale basées sur la tomographie d'émission par positrons (TEP) qui sont utilisées aussi bien pour le diagnostic clinique que pour la recherche biomédicale sur modèle animal. Les interactions du groupe avec l'industrie ont mené à la mise au point d'instruments et de nouvelles technologies d'étude scientifique et clinique.

Mon stage faisait la partie de recherche sur le problème de diffusé dans –la tomographie d’émission par positrons (TEP). Le problème global était l’estimation par calculs analytiques du rayonnement diffusé en TEP.

 

1. Bref description due problème traité

Les images obtenues en TEP sont basées sur la détection simultanée de deux rayonnements gamma émis à 180° l’un par rapport à l’autre. Cependant, un photon peut interagir avec la matière et pourrait diffuser. Le diffusé, non seulement il réduit le signal dans l’image, mais aussi il génère des artefacts et des imprécisions. Différentes méthodes de correction sont actuellement utilisées pour approximer la distribution du diffusé: par des sources ponctuelles; par estimation à l’extérieur de l’objet mesuré; par calcul analytique en supposant une seule interaction d’un seul photon de la paire de photons émis. Et pour toutes ces méthodes, les photons absorbés ou diffusés sont d’abord compensés avec l’aide de mesures de la transmission, ensuite le diffusé est soustrait des images.

Ce projet de recherche consiste à intégrer les corrections de la résolution spatiale et du rayonnement diffusé dans les algorithmes de reconstruction tomographique en TEP. Le système de détection est influencé par la géométrie du tomographe et par la conception des détecteurs. La réponse du système est mesurée avec des mires, puis elle est intégrée dans les algorithmes de reconstruction. La simulation Monte Carlo permet d’observer les différents ordres du diffusé dans l’objet qui devrait être corrigé.

2. Les tâches effectuées

Pendant ce stage, j'ai écris un programme pour calculer la probabilité de diffusion des photons dans le tomographe en ne considérant qu’un seul anneau de détecteurs, et en connaissant au préalable l’origine de la position de la diffusion du photon dans le l’objet.

Le programme que j'ai écrit considère les coordonnées d'un point et retourne une matrice des probabilités sur les angles de diffusion (j'ai aussi écrit la version pour les coordonnées polaires).

Le programme a été écrit en Matlab. Le choix de ce langage était déterminé par les raisons suivantes:

    1. C'est le langage traditionnellement utilisé par mon équipe de recherche. Ainsi son utilisation facilitait notre collaboration.
    2. Matlab donne un très grand nombre d’outils essentiels pour un programmeur: facilité a) de déclaration et de re-déclaration des variables et des tableaux; b) d'affichage des images; c) de sauvegarde des donnés et de leur chargement dans la mémoire pendant l'exécution d'un programme, etc.. En plus les programmes sont faciles à changer et à corriger, car il n'y a pas de compilation comme telle.
    3. Les programmes écrits en Matlab peuvent être facilement convertis dans le code C.

 

Figure 1. La circonférence interneverte délimite le champ utile, espace réservé aux objets mesurés. La circonférence externebleue représente l'anneau de détecteurs.

Pour tout point (y0, z0) du champs utile du tomographe (fig.1), où le diffusé est possible, et pour toute direction g d'un photon avant ce diffusé on doit calculer des probabilités du fait que ce photon sera enregistré par chacun des détecteurs de l'anneau.

 

3. Les approches utilisées

La probabilité de chaque angle de diffusion est calculée avec la formule de Klein-Nishina:

(1)

Cela nous permet de calculer un tableau des probabilités de chaque angle de diffusion q entre 0° et 180° (fig. 2a et 3).

Ainsi, pour effectuer notre tâche pour chaque détecteur il faut trouver tous les angles q qui correspondent à lui.

 

 

 

 

(a) (b)

Figure 2. a) Diffusé dans le champ utile du tomographe. Aprèsr avoir changé sa direction par un angle q , ce photon sera enregistré par un détecteur . b) Pour chaque détecteur il y a un certain nombre d'unités des' angles (D q ) qui lui correspondent à lui. Les valeurs de D q dépendent des coordonnées du point de diffusioné, surtout de distance entre ce point et les détecteurs.

 

 

 

Figure 3. DistributionLe tableau des probabilités de chaque angle de diffusion q entre 0° et 360° , calculé avec la formule de Klein-Nishina. Mon programme sauvegarde ces données dans une variable appelée Klein_Nishina_Array .

Ensuite on peut trouver la somme des probabilités qui correspondent à ces angles. En bref:

(2)

Alors pour chaque detX il faut trouver son D q (fig. 2b).

J'ai essayé les trois approches. La première prévoyait le grand nombre de calcules trigonométriques. Bien que les résultats étaient assez précis, cette méthode prenait trop de temps pour être utilisée. Parfois même 8 heures n'étaient pas suffisantes pour calculer la matrice des probabilités pour un seul point. Deuxième approche : diviséere l'anneau de détecteurs àen quelques sections ( 4, 8 ou 16). Ensuite pour les détecteurs de chaque section on trouvait les angles D q d'une façons proportionnelle. Malheureusement, cette approche, bien que rapide, ne pouvait pas donner des résultats précis.

Finalement, j'ai réussi à écrire un programme qui donne des bons résultats et ne prendne que quelques secondes pour les calculsêtre executé.

 

4. Brève description du programme

Il faut remarquer d'abord que ce programme a été fait pour le tomographe particulier (Scanner TEP Animal de Sherbrooke), mais peut être facilement adapté pour un autre tomographe de même type.

On voit facilement que le tableau des D q (dans notre programme on l'appelle "thetas") sera le même pour chaque point du champs utile. C'est plus facile de travailler avec des coordonnées polaires. Par exemple un point devient un point . Alors on calcule "array_of_thetas36000" (la matrice des D q ) qui contient le tableau des unités des angles par détecteur pour chaque point du champ utile.

Pour calculer le tableau "thetas" pour chaque point , j’utilise l’équation suivante:

(3)

Note: est la longueur d'un seul détecteur et est le rayon du tomographe.

Les calculs de "array_of_thetas36000" ne prennent que 3 minutes. Une fois "array_of_thetas36000" est calculé, on le sauvegarde dans la mémoire. Matlab permet le stockage et le chargement rapide pour les tableaux des données. Ainsi, on garde l'information sur les nombres des unités des angles pour chaque détecteur, et on le charge pendant l'exécution des programmes suivants. Cela réduit économise dule temps de calcul.

La prochaine étape c'est de calculer les probabilités pour chaque détecteur, en utilisant l’équation (2).

5. Les résultats obtenus

(a)(b)

(c)

Figure 4.a) Tableau "thetas" qui correspond au point (50 mm, 80° ) du champ utile.

b) Courbe des probabilités pour g =0° obtenue à partir de "thetas" et Klein_Nishina_Array" pour le même point. c) Matrice des probabilités pour g =0° -..360° pour ce même point (50 mm, 80° ).

 

Dans la figure 4 on peut voir des résultats, obtenus grâce au programme "angleprobability" et les programmes reliés. La figure 5 présente quelques matrices. Comme dans fig. 4.c., les probabilités sont indiquées par les couleurs, tandis que les axes X et Y indiquent les détecteurs et les valeurs d'angles g (voir fig. 1 et fig. 2.a).

(a)(b)

(c)

Figure 5. a) Matrice des probabilités pour g =0° -..360° pour le diffusé produit auarrivé dans le point avec des coordonnées =0 mm, =0 mm. b) Matrice des probabilités pour g =0° -..360° pour le diffusé produit auarrivé dans le point avec des coordonnées =20 mm, =30 mm. c) Matrice des probabilités pour g =0° -..360° pour le diffusé produit auarrivé dans le point avec des coordonnées =50 mm, =-10 mm.

Après avoir ajoutéer de nombreuses améliorations, l'exécution du programme "angleprobability" ne prend que 10 secondes.

 

 

Conclusions

Le travail effectué jusque là n'est qu'un prologue pour les travaux à venir. Il faut créer des outils dont la vitesse d'exécution est encore plus rapide. Bien que très rapide, la version actuelle n'en est pas assez. Le champ utile du tomographe étudié contient quelques 21240 points. Ainsi, les calcules de matrice des probabilités pour chacun d'entre eux avec "angleprobability" prendra plus que 59 heures. Il faut alors chercher des moyens pour augmenter la vitesse d'exécution.

La deuxième tâche d'avenir c'est de rendre le programme plus intégral et plus facile à utiliser. Aujourd'hui, pour obtenir une matrice comme celles de la figure 5, il faut lancer successivement quelques programmes. Dans la version future il faut créer un programme qui va gérer tous ces modules, pour que l'utilisateur puisse les exécuter sans en penser trop et sans perdre le temps en passant d'un module à l'autre. À cette fin je vais mieux profiter des avantages de la programmation orientée objet. Une fois aussi développé,er l'interface complet et facile à utiliser.

Ce stage m'a permit de:

    1. Apprendre Matlab
    2. Se familiariser avec la technologie de la recherche scientifique
    3. Améliorer mes connaissances dans le domaine de la radiobiologie
    4. Se familiariser avec des applications de moyens informatiques
    5. Développer mes côtés créatifs en cherchant des solutions de problèmes et réels

 

Annexe A

Le code écrit:

%Written by E.Lakinsky

%To evaluate probability for each point (dist,betyz) of field of view

%The probability is calculated only for the first scatter.

%The photon's initial energy is 511 keV

%Polar coordinates version

function sinoprob=angleprobability(dist,betyz);

tic;

rfov = 59;

rtomo = 155; % Ring's radius

detring = 256;

det0=todet(asin(dist*sin(betyz*pi/180)/rtomo)*detring/(2*pi));

load my_klein_nishina_prob_array36000;

load array_of_thetas36000;

thetas(1,:)=array_of_thetas36000(round(dist)+1,betyz+1, :);

detring360=detring/360;

sinoprob=zeros(360,detring);

for gamma=0:359;

detn=todet(gamma*detring360);

det0detn=todet(det0+detn);

firstangle=1;

for x=det0detn:detring+det0detn;

if x<=detring;

thisangle=sum(thetas(det0detn:x));

else;

thisangle=sum(thetas(det0detn:detring))+sum(thetas(1:x-detring));%sum(thetax(1:todet(x)))

end;

if thisangle<=36000;

pr(todet(x))=sum(my_klein_nishina_prob_array36000(firstangle:thisangle));

firstangle=thisangle;

else;

pr(todet(x))=sum(my_klein_nishina_prob_array36000(firstangle:36000))+...

sum(my_klein_nishina_prob_array36000(1:thisangle-36000));

firstangle=thisangle-36000;

end;

end;

pr=pr/sum(pr);

sinoprob(gamma+1,:)=pr(1,:);

end;

figure;imagesc(sinoprob);title(['Probabilité pour le point r=' num2str(dist) ' alpha=' num2str(betyz)]);xlabel('Détecteurs');ylabel('Angles gamma, degrés');

toc

% To compute the Array_of_Thetas

function array_of_thetas=calc_thetas;

tic;

rfov = 59;

rtomo = 155; % Ring's radius

detring = 256;

%hnus=511;

detw = rtomo*2*pi/detring; % Detector's width with package.

k1=detw/2;

k2=rtomo^2;

k3=2*rtomo;

k4=pi/180;

k5=2*pi/detring;

k18000pi=18000/pi;

array_of_thetas=zeros(rfov+1,360,detring);

array_of_thetas(1,:,:)=140;

for r=1:rfov;

r1=r+1;

k2r2=k2+r^2;

k3r=k3*r;

for alpha=0:359;

alpha1=alpha+1;

alphak4=alpha*k4;

for x=1:detring;

array_of_thetas(r1,alpha1,x)=2*asin(k1/sqrt(k2r2-k3r*cos(alphak4-x*k5)));

end;

end;

end;

%size(array_of_thetas) = 60 360 256

save array_of_thetas;

toc

function array_of_thetas36000=calc_thetas36000;

tic;

rfov = 59;

rtomo = 155; % Ring's radius

detring = 256;

k18000pi=18000/pi;

array_of_thetas=zeros(rfov+1,360,detring);

load array_of_thetas;

array_of_thetas36000(1,:,:)=array_of_thetas(1,:,:);

for r=1:rfov;

r1=r+1;

k2r2=k2+r^2;

k3r=k3*r;

for alpha=0:359;

alpha1=alpha+1;

a(1,:)=array_of_thetas(r1,alpha1,:)*k18000pi;

array_of_thetas36000(r1,alpha1,:)=floor(a(1,:)*36000/sum(a));

end;

end;

%size(array_of_thetas36000) = 60 360 256

save array_of_thetas36000;

toc

Hosted by www.Geocities.ws

1