B.2 Modèle poissonien et l'algorithme OSEM-OSL

L'émission étant un phénomène statistique, trouver un algorithme itératif qui intègre cette information dans sa résolution ne pouvait que séduire. C'est le cas de l'algorithme EM (pour Expectation Maximisation) [31,55]. Cet algorithme est de loin le plus connu et le plus utilisé des algorithmes itératifs de reconstruction.

B.2.1 l'algorithme EM-OSL[41]

Nous avons vu qu'en reconstruction et dans le cas d'une modélisation poissonienne et dans le cadre d'une régularisation markovienne, le critère du MAP que nous devions minimiser s'écrivait:

$$J(\mathbf{f}) = J_{1}(\mathbf{f})+\beta ^{2}J_{2}(\mathbf{f})$$

$$\Downarrow$$

$$J(\mathbf{f}) = \sum _{m}p_{m}\log (\hat{p}_{m})-\sum _{m}\hat{p}_{m}+\beta ^{2}J_{2}(\mathbf{f})$$

La constante $\beta$, que nous rajoutons, nous permet de controler le poids de la régularisation ($J_{2}$) vis à vis de l'importance accordée aux données ($J_{1}$). La plupart du temps, on introduit à ce niveau une variable auxiliaire non-observée $z_{mn}$, les données complètes représentant le nombre de photons émis par le voxel $n$ et captés par le détecteur référencé par $m$. Nous avons alors le nombre de coups détectés par un capteur qui s'écrit comme:

$$p_{m}=\sum _{n}z_{mn}$$

L'algorithme consiste alors en deux phases: l'estimation des $z_{mn}$ (Expectation) et la maximisation de $I\! \! P(\mathbf{z}\vert\mathbf{f})$ (Maximisation). En combinant les deux étapes, on arrive à une expression [41]:

$$\{f_{n}\}_{\alpha +1}=\frac{\{f_{n}\}_{\alpha }}{\sum _{m}H_{mn}+... ...)}{\delta f_{n}}\right) }\sum _{m}H_{mn}\frac{p_{m}}{\{\hat{p}_{m}\}_{\alpha }}$$

où $\hat{p}_{m}=\sum _{n}H_{mn}\{f_{n}\}_{\alpha }$. Evidemment dans cette expression, le terme $\frac{\delta J_{2}(\mathbf{f})}{\delta f_{n}}$ qui représente la dérivée locale de notre information a priori, nous pose problème (cf Par.9.2.4.5). Green propose, à l'itération $\alpha +1$, de calculer ce terme en partant de l'estimée à l'itération $\alpha$ (d'où le nom de l'algorithme OSL pour One Step Late). Nous avons alors, en utilisant l'expression de la dérivée locale (Eq.9.5):

$$\{f_{n}\}_{\alpha +1}=\frac{\{f_{n}\}_{\alpha }}{\sum _{m}H_{mn}-... ...a })\mathbf{f}_{\alpha }}\sum _{m}H_{mn}\frac{p_{m}}{\{\hat{p}_{m}\}_{\alpha }}$$ (B.1)

En l'état, l'implémentation de cette expression conduit à l'algorithme EM-OSL. Un inconvénient quant à l'utilisation de cet algorithme est la lenteur de la convergence. C'est pourquoi, dans un souci d'accélérer cette convergence nous utiliserons la technique Ordered Subset (OS), qui conduit à l'algorithme que nous avons implémenté OSEM-OSL.

B.2.2 l'algortithme OSEM-OSL

La technique OS comprend deux concepts qui permettent l'accélération: le réarrangement des projections (Ordered) et leur partitionnement (Subset)

B.2.2.0.1 Réarrangement des projections.

Ce réarrangement peut se comprendre de manière intuitive. Notre sinogramme correspond à un ensemble de plan $\Pi$ tournant autour de l'objet à imager. Le nombre de vue $N_{\theta }$ est grand (144), autrement dit la différence angulaire entre deux plans succesifs (2 $\theta$ successif) est faible (1.25$°$). Il va sans dire que les informations apportées par chacun de ces deux plans sur l'image $\mathbf{f}$ sont proches. Il en va tout autrement, si entre deux itérations successives, les angles de vue sont éloignés (grande différence angulaire). Si donc, on veut accroitre la vitesse de convergence de l'algorithme, il faut fournir d'une itération à l'autre une information qui soit la plus différente possible de celle qu'il vient de recevoir. Intuitivement, on se dit qu'un écart angulaire $\Delta \theta =\theta _{k+1}-\theta _{k}=90$$°$ serait optimal. Nous envisageons le réarrangement tel qu'il est décrit par Herman [43]. Supposons que nous ayons 10 directions de projections, plutôt que de parcourir les projections dans l'ordre 1,2,3,... nous les parcourons de sorte d'avoir toujours le plus grand écart angulaire (Tab.B.1).

Table: Réarrangement des projections

Parcours usuel
projection	1	2	3	4	5	6	7	8	9	10
$\theta$	0	18	36	54	72	90	108	126	144	162
$\Delta \theta =18$$°$
Parcours après réarrangement
projection	1	6	2	7	3	8	4	9	5	10
$\theta$	0	90	18	108	36	126	54	144	72	162
$\Delta \theta$		90	72	90	72	90	72	90	72	90

B.2.2.0.2 Création de subset

Toujours en travaillant sur les redondances liées au différentes vues, on peut se poser la question de savoir ce que donnerait la reconstruction en sous échantillonant le sinogramme (Ne prendre qu'un angle de vue tous les $K$ ). On crée un sous ensemble du sinogramme (subset). Avec un sinogramme on peut évidemment construire $K$ subsets.

$$S_{1}=\left\{ p_{i,j,1},p_{i,j,K+1},\ldots \right\} ,S_{2}=\left\... ...2},\ldots \right\} ,\ldots ,S_{K=}\left\{ p_{i,j,K},p_{i,j,2K},\ldots \right\}$$

Partant de cette subdivision, on va reconstruire le volume pour un subset. L'image $\mathbf{f}_{S_{1}}$ obtenue sert de point de départ pour une reconstruction sur le deuxième subset, et ainsi de suite. Le passage sur tous les subsets va constituer une itération (un passage de toutes les projections), que l'on répète jusqu'à convergence. En parcourant les subsets dans l'ordre, on obtient juste une technique de réarrangement qui consiste à prendre les projections dans l'ordre $p_{i,j,1},p_{i,jK+1},\ldots ,p_{i,j,2},p_{i,j,K+2},\ldots$. En revanche, on peut combiner les deux approches en définissant un ordre de parcours des subsets. De sorte que l'information apportée par deux subsets successifs soit la plus différente possible. Cette approche définit l'algorithme Ordered Subset (Alg.13)