4.2 La Projection Analytique.

Dans ce type d'approche, on va chercher à estimer de manière globale, le nombre d'émissions vues par un couple de détecteurs et donc par un quadruplet particulier $(r,s,\theta ,\phi )$ du sinogramme.

4.2.1 Sinogramme et projection.

Toute portion radioactive de l'espace émet des photons $\gamma$ en opposition de manière isotrope. Pour cette portion de l'espace, une petite partie seulement de ces émissions va se produire suivant la direction $\vec{n}$. Pour que, de plus, ces photons $\gamma$ soient comptabilisés par le couple de détecteurs $(A,B)$ (on suppose pour l'instant l'absence de phénomènes perturbateurs), il faut qu'ils soient émis sur la droite $(AB)$ joignant ces deux capteurs: la portion de l'espace envisagée doit appartenir à cette droite. Ainsi, tout couple de détecteurs $(A,B)$ va enregistrer toutes les paires de photons $\gamma$ crées dans une portion de l'espace appartenant à la droite $(AB)$ et émises suivant la direction $\vec{n}$. Calculer un sinogramme revient à affecter à tout quadruplet $(r,s,\theta ,\phi )$ caractérisant un couple de détecteurs, le nombre $p$ de photons $\gamma$ comptés suivant cette ligne de coïncidence durant toute la durée d'un examen:

$$p(r,s,\theta ,\phi )=\int _{k_{1}\in \mathbb{R}}f(\overrightarrow{O'P}+k_{1}.\vec{n})dk_{1}$$ (4.1)

Une valeur dans le sinogramme $p(r,s,\theta ,\phi )$ correspond à la somme des contributions de l'objet radioactif sur la ligne de coïncidence.

4.2.2 Choix de l'origine $O'$.

Lorsque nous avons défini le choix de l'origine $O'$ (origine dans le plan $\Pi$ de projection cf. Ch.4), nous avons dit qu'elle pouvait se situer n'importe où sur la ligne de coïncidence. C'est à dire que si nous décalons l'origine en $O_{2}'$ tel que $\overrightarrow{O'O_{2}'}=k_{2}.\vec{n}$ alors:

$$\int _{k_{1}\in \mathbb{R}}f(\overrightarrow{O_{2}'P}+k_{1}.\vec{n})dk_{1} = \int _{k_{1}\in \mathbb{R}}f(\overrightarrow{O'P}+(k_{1}-k2).\vec{n})dk_{1}$$

$$= \int _{k_{3}\in \mathbb{R}}f(\overrightarrow{O'P}+k_{3}.\vec{n})dk_{3}$$

$$= p(r,s,\theta ,\phi )$$

avec $k_{3}=k_{2}-k_{1}$. On définit bien là, au sens mathématique du terme, une projection suivant la direction $\vec{n}$, du fait de cette invariance de $p(r,s,\theta ,\phi )$ lorsqu'on déplace l'origine $O'$ suivant $\vec{n}$. Sa position est donc sans importance sur les valeurs du sinogramme.

4.2.3 Projection et fonction caractéristique.

Jusqu'à présent, nous avons toujours considéré que seules les portions de l'espace situées sur la droite $(AB)$ et dans la direction $\vec{n}$ contribuaient à la valeur $p(r,s,\theta ,\phi )$ du sinogramme. Considérons maintenant l'expression plus générale:

$$p(r,s,\theta ,\phi )=\int \int \int _{\mathbb{R}^{3}}f(\vec{x}).\Phi _{r,s,\theta ,\phi ,\vec{x}}(\vec{x})d\vec{x}$$ (4.2)

où $\Phi _{r,s,\theta ,\phi ,\vec{x}}$ est une fonction caractéristique marquant les portions de l'espace vues par le couples de détecteurs $(r,s,\theta ,\phi )$. Dans l'expression Eq.5.2, si la fonction $\Phi _{r,s,\theta ,\phi ,\vec{x}}$ vaut 1 uniquement pour les $\vec{x}$ appartenant à la ligne de coïncidence définie par $(r,s,\theta ,\phi )$ et 0 partout ailleurs, nous retrouvons l'expression Eq.5.1 de la projection précédemment citée. Cette fonction $\Phi$ peut devenir très générale, car elle spécifie la portion de l'espace vue par les détecteur $(A,B)$. Elle dépend du couple de détecteur $(r,s,\theta ,\phi )$ considéré mais aussi du point de la position $\vec{x}$ dans l'espace. Par la suite, on considère que cette fonction caractéristique ne dépend plus de la position $\vec{x}$ et se réduit donc à $\Phi _{r,s,\theta ,\phi }$ (Invariance spatiale de la fonction caractéristique).

4.2.4 Choix d'un algorithme de projection.

Le calcul d'un sinogramme passe par le choix d'une méthode pour estimer la projection d'un volume sur un plan suivant une direction $\vec{n}$. Pour expliciter les algorithmes liés à la projection et afin d'alléger notations et figures, nous allons aborder le problème en deux dimensions.^4.1 Le cas tridimensionnel s'en déduit naturellement, les principes restant les mêmes. Notre objet devient donc bidimensionnel $f(x,y)$. Sa version discrétisée reste $\textbf {f}$. Le sinogramme correspondant se définit simplement par la projection de cet objet sur une droite $\Delta$. Cette droite est référencée dans l'espace par la donnée du couple $(r,\theta )$ (Fig.5.2).

Figure: Représentation schématique de la projection en 2D.

La projection à calculer devient:

$$p(r,\theta )=\int \int _{\mathbb{R}^{2}}f(x,y)\Phi _{r,\theta }(x,y)dxdy$$

Ce qui, dans un monde discret, se traduit par:

$$p_{m}=\sum _{n\in [1,N_{x}]\times [1,N_{y}]}H_{mn}f_{n}$$ (4.3)

avec $m=r_{i}+N_{r}\times \theta _{k}$. L'équation Eq.5.3 constitue une équation particulièrement importante que nous retrouverons longuement lors de la reconstruction algébrique (Ch.9). Dans cette équation, $H_{mn}$ traduit la fraction de $f_{n}$ qui est vue par le capteur $m$.

4.2.4.0.1 Remarque:

A ce niveau, Il semble indispensable de faire une remarque sur le parti pris par la suite pour la représentation de l'objet. Nous avons exposé brièvement qu'il existait de nombreuses représentations de l'image. Parmi celles-ci, la représentation spatiale $(x,y,z)$ dans un contexte continu (celle des voxels en discret) n'est qu'une alternative. Sa popularité vient de sa simplicité d'emploi qui donne sans intermédiaire une représentation visuelle de l'objet. En effet, Il y a équivalence entre les coefficients décrivant l'image et la perception que nous en avons. Dans l'équation Eq.5.3, les coefficients $f_{n}$ sont étroitements liés à la base de représentation de notre objet. D'ailleurs $H_{mn}$ représente la contribution d'une fonction élémentaire de la base de représentation à la projection. Le calcul de $H_{mn}$ pour différentes bases de représentations peut être trouvé dans [56].

Pour construire des projections partant d'un volume représenté par des voxels, deux approches ont été envisagées.

4.2.4.1 Approche par le parcours de l'image.

Cette approche est la plus simple à mettre en oeuvre (Fig.5.3)

Figure 5.3: Projection par le parcours de l'image.

. Pour une direction de projection possible (une direction $\vec{n}$ fixée par un $\theta _{k}$), on parcourt tous les pixels de l'image. Un pixel d'indice $n$ correspond à une position bien précise dans l'espace $(x,y)$. Il s'agit en général du centre du pixel. On détermine où ce point se projette sur la droite $\Delta$ par une projection suivant $\theta _{k}$.

$$r=-x.\sin \theta _{k}+y.\cos \theta _{k}$$

La valeur entière de $r$ nous conduit à un indice $i$ du sinogramme-tableau $r_{i}$ tel que $r\in [r_{i},r_{i+1}]$. On cherche les indices $m$ du sinogramme $\textbf {p}$ correspondant au couples de détecteurs vus de ce pixel (réponse impulsionnelle) et on incrémente les valeurs $p_{m}$ des contributions $H_{mn}$ apportées par $f_{n}$.

$$p_{m}=p_{m+}H_{mn}f_{n}$$

Cette méthode élémentaire se décline sous de nombreuses variantes en fonction de la forme de la réponse impulsionnelle choisie. On distribue l'information contenue dans le pixel $n$ en fonction de la distance $d_{i}$ séparant le centre d'un détecteur d'indice $i$ du point réel de projection $d_{i}=r-r_{i}$. A titre d'illustration, on donne deux exemples correspondant à deux réponses impulsionnelles différentes (Fig.5.4)

Figure: Distribution de l'information d'un pixel pour 2 réponses impulsionnelles différentes.

[Interpolation linéaire.] [Interpolation Gaussienne.]

.

4.2.4.1.0.1 Interpolation linéaire [15].

Dans ce cas, on ne répartit la contribution $f_{n}$ sur le détecteur $r_{i}$ le plus proche par valeur inférieure ($r_{i}<r$) et sur son successeur $r_{i+1}$, on effectue alors une simple interpolation linéaire

• la contribution $H_{mn}=1-\frac{d_{i}}{r_{i+1}-r_{i}}$ est affectée à $p(r_{i},\theta _{k})$.
• la contribution $H_{mn}=\frac{d_{i}}{r_{i+1}-r_{i}}$ est affectée à $p(r_{i+1},\theta _{k})$.

4.2.4.1.0.2 Interpolation Gaussienne.

Il est admis qu'une gaussienne est une bonne approximation permettant de modéliser la contribution des photons directs [35,8]. Sa forme générale dépend principalement de l'angle solide sous lequel est vu le pixel par les détecteurs. Toutefois, si on considère une réponse impulsionnelle gaussienne d'écart type $\sigma$ fixe (invariance spatiale), on affecte une contribution au détecteur $(r_{i},\theta _{k})$ en fonction de

$$H_{mn}=\frac{1}{\sqrt{2\pi }\sigma }\exp [-\frac{(r-r_{i})^{2}}{2\sigma ^{2}}]$$

4.2.4.1.1 Inconvénient de la méthode.

L'inconvénient majeur de cette méthode vient du fait que l'image comme les projections sont échantillonnées. Quand on parcourt l'image, le déplacement se fait par des incréments régulier de $dx$ où $dy$ pixels. Ces déplacements, une fois projetés, se traduisent par un pas d'échantillonnage suivant la direction $\vec{r}$ tel que $dr=-dx\sin \theta _{k}$ ou $dr=dy\cos \theta _{k}$. Le sinogramme est lui même déjà échantillonné avec un pas $r_{i+1}-r_{i}$. Le ratio $\frac{dr}{r_{i+1}-r_{i}}$ n'est pas constant et dépend de la direction de projection $\theta _{k}$. De plus, il existe toujours une direction de projection qui conduit à une situation critique produisant un effet d'aliasing (Fig. 5.5).

Figure: Effet d'aliasing dû au double échantillonnage. Il existe toujours une direction de projection conduisant à un effet d'aliasing. Sur l'exemple, pour une simple ligne d'un pixel d'épaisseur, soit un pixel soit deux se projettent sur un détecteur.

Dans ce cas, une ligne continue et uniforme sur l'image conduit à un sinogramme dentelé. L'interpolation linéaire, et a fortiori l'interpolation gaussienne atténue ce phénomène par l'introduction d'un lissage, mais ne le corrige pas. Pour remédier à ce problème, lié à l'échantillonnage régulier de notre image, nous allons bruiter légèrement les positions de pixels $(x+\delta x,y+\delta y)$. $(\delta x,\delta y)$ est un vecteur de translation tiré au hasard dont la norme est de dimension inférieure à la taille du pixel. A titre d'illustration, nous donnons Fig.5.6.a (resp. Fig.5.6.b) un sinogramme calculé sans bruitage (resp avec bruitage) de la position des voxels.

Figure: Un plan $\Pi$ d'un sinogramme calculé par le parcours du volume.

On voit très nettement Fig.5.6.a, l'effet lié à l'aliasing qui se traduit par l'apparition de bandes verticales sur le plan $\Pi$ du sinogramme considéré.

4.2.4.2 Approche par le parcours du sinogramme.

Dans l'approche précédente, nous parcourions le volume pour savoir à quel endroit chaque voxel se projetait dans le sinogramme. A l'inverse, on veut savoir ici quels sont les voxels qui contribuent à une projection donnée. On va donc parcourir tous les dexels. Chaque dexel du sinogramme, i.e. tout indice $m$, référence un point $P$ particulier dans le plan $\Pi$ de projection. D'autre part, cet indice $m$ fixe également une direction par la donnée des deux angles $(\theta ,\phi )$. Partant du point $P$, on parcourt l'espace suivant cette direction (sorte de lancer de rayon suivant la ligne de coïncidence) et on intègre en $p_{m}$ les contributions, en terme de nombre de coups, de tous les éléments de volume rencontrés sur notre parcours de la ligne de coïncidence. Lors de ce parcours , il faut être certain d'avoir traverser tout le volume émetteur. Or, nous avons vu que les plans de projections $\Pi$ sont invariants pour un déplacement suivant la direction de projection $\vec{n}$. On va donc placer le plan $\Pi$ tangent à la sphère $\mathcal{C}$ englobant tout le volume servant à la projection (Fig.5.7.a).

Figure 5.7: Projection d'un volume par le parcours du Sinogramme.

[Choix de l'origine du plan $\Pi$] [Parcours de la droite dans l'espace discret du volume.]

Ceci fixe une position $(x_{P},y_{P},z_{P})$ précise du point $P$. Ce point $P$ est ensuite référencé dans une grille analogue à celle servant à échantillonner le volume (grille régulière constituée par les centres des voxels). Le parcours du volume suivant la direction de projection se résume alors à suivre la ligne de coïncidence dans un espace discret (Fig.5.7.b). La question est alors: ``Quel sont les voxels du volume qu'il faut compter comme appartenant à la droite de projection?'' Ce problème est classique, il s'agit de représenter une ligne non pas dans un monde continu, mais dans un monde discret présentant une trame plus ou moins grossière. Ce problème a été résolu par Jack Bresenham, employé à IBM, qui travaillait à l'époque sur les imprimantes matricielles. L'algorithme de Bresenham (voir Annexe) tout d'abord écrit en 2D fait référence en la matière. C'est cet algorithme que nous utiliserons.