8.1 Introduction.

8.1.1 Monde continu versus monde discret.

Attardons nous quelques instants sur la problématique de la reconstruction! Nous avons d'un coté un jeu de données $p(r,s,\theta ,\phi )$ fournies par notre imageur, et de l'autre coté la volonté de connaître la valeur de la distribution $f(x,y,z)$ introduite dans l'imageur. Dans notre description du sinogramme, nous avons montré comment partant de $f$ nous pouvions construire un jeu de données $p$. Il s'agissait du problème direct. Lors de la reconstruction, nous cherchons à inverser ce processus.

Dans une approche analytique, le lien qui unit la projection au volume est établi dans un monde continu. Nous envisageons simplement la projection comme une intégration suivant une ligne de coïncidence (Fig.9.1.a). La reconstruction représente alors une inversion analytique de cette intégration. Or, les données de projection et le volume doivent être discrétisés pour pouvoir être mémorisés. Il faut échantillonner $f$ et $p$ (discrétisation). La nature discrète des données n'est envisagée qu'au moment de l'archivage/relecture des données.

Figure 9.1: Monde Continu versus Monde discret.

[Monde Continu.] [Monde discret.]

Dans une approche algébrique, le problème direct est appréhendé dans un contexte discret (Fig.9.1.b). Distribution et projections sont alors deux vecteurs monodimensionnels $\mathbf{f}$ et $\mathbf{p}$. Sur une base de voxels, chaque élément $\{f\}_{n}$ traduit la valeur moyenne sur un élément de volume (voxel) référencé par $n$. De même $\{p\}_{m}$ représente le nombre de photons $\gamma$ détectés par un élément de surface d'un plan de projection référencé par $m$ (dexel). De surcroît, chaque composante de l'un quelconque de ces deux vecteurs ne peut prendre qu'un nombre fini de valeurs. Cette valeur correspond à un état particulier dans lequel se trouve le voxel ou le dexel (quantification). Il y a $\Lambda _{f}$ états possibles pour chaque voxel $f_{m}$ et $\Lambda _{p}$ pour chaque dexel $p_{m}$. $\Lambda _{f}$ (ou $\Lambda _{p}$) représente l'espace des états pour les voxels (les dexels). Ce nombre est fonction du type (short, float, ...) dans lequel sont stockées distribution et projections. De par cette nature discrète, il existe un nombre fini $\Lambda ^{N}_{f}$ de volumes possibles, où $N$ représente le nombre de voxels.

Reconstruire, c'est donc choisir parmi ces $\Lambda ^{N}_{f}$ volumes possibles celui qui semble le plus approprié au vu des projections mesurées.

On définit une fonction de coût $J(\mathbf{f})$ qui, pour un jeu observé de projections, donne une idée de la pertinence de chaque image $\mathbf{f}$ considérée. La reconstruction revient à minimiser cette fonction de coût. Elle s'inscrit donc dans le domaine plus général de la minimisation de fonction. Il est évident que le choix de la fonction $J$ conditionne la qualité de la reconstruction. La richesse de l'approche algébrique provient justement de ces multiples fonctions $J$ envisageables. Le lien entre les projections et le volume n'est plus rigide, comme il pouvait l'être dans le cas d'une reconstruction analytique.

Nous limitons notre étude à la modélisation de la projection par un lien linéaire. Les projections $\{\hat{p}\}_{m}$ sont alors une somme pondérée des différentes composantes $\{f\}_{n}$ de l'image, ou en envisageant une écriture matricielle:

$$\hat{\mathbf{p}}=\mathbf{H}\mathbf{f}$$ (8.1)

8.1.2 Approche bayésienne.

Pour bien comprendre le sens de la reconstruction algébrique, nous l'introduisons dans un contexte statistique. Dans une approche statistique, l'image $\mathbf{f}$ est considérée comme un champ aléatoire. La reconstruction se résume alors à ceci:

Comment, pour le jeu de données $\mathbf{p}$ observées et parmi les $\Lambda ^{N}_{f}$réalisations possibles de $\mathbf{f}$, choisir la distribution la plus probable. Autrement dit, on veut maximiser la probabilité de $\mathbf{f}$ sachant $\mathbf{p}$ ou probabilité a posteriori: $I\! \! P(\mathbf{f}\vert\mathbf{p})$. D'après Bayes, cette probabilité peut s'exprimer différemment:

$$I\! \! P(\mathbf{f}\vert\mathbf{p})=\frac{I\! \! P(\mathbf{p}\vert\mathbf{f}).I\! \! P(\mathbf{f})}{I\! \! P(\mathbf{p})}$$ (8.2)

où $I\! \! P(\mathbf{p}\vert\mathbf{f})$ représente la vraisemblance des projections pour une image $\mathbf{f}$ particulière, $I\! \! P(\mathbf{f})$ la probabilité d'observer la dite image $\mathbf{f}$, et $I\! \! P(\mathbf{p})$ la probabilité d'obtenir ce jeu de projection. Dans une optique de reconstruction, le jeu de données de projection est fixé puisqu'il correspond à une acquisition. Autrement dit, la probabilité $I\! \! P(\mathbf{p})$ est constante. Elle intervient seulement comme une constante de normalisation dans l'Eq.9.2. Par conséquent, la maximisation de $I\! \! P(\mathbf{f}\vert\mathbf{p})$ ne fait intervenir que la vraisemblance des données (probabilité d'observer le jeu de projection $\mathbf{p}$ pour la distribution $\mathbf{f}$ connue) et la probabilité de $\mathbf{f}$ ou probabilité a priori (ce que je sais de $\mathbf{f}$). Pour maximiser la probabilité a posteriori, il faut construire un estimateur statistique $\mathbf{e}$ de $\mathbf{f}$. Pour choisir cet estimateur, il est nécessaire d'introduire une fonction de coût qui mesure la qualité de l'estimation $\mathbf{e}$ par rapport à la solution exacte. Dans nos travaux, nous n'envisageons que l'estimateur au sens du Maximum A Posteriori (MAP) qui nous conduit à choisir $\mathbf{f}$ tel que:

$$\mathbf{f}=arg\max _{\mathbf{e}}I\! \! P(\mathbf{e}\vert\mathbf{p... ...f}=arg\max _{\mathbf{e}}I\! \! P(\mathbf{p}\vert\mathbf{e})I\! \! P(\mathbf{e})$$

Maximiser ce critère revient à minimiser $-\log I\! \! P(\mathbf{e}\vert\mathbf{p})$ car la fonction logarithme est croissante et monotone. Nous obtenons donc, au sens de MAP un critère de choix de l'image $\mathbf{f}$ tel que:

$$\underbrace{-\log I\! \! P(\mathbf{f}\vert\mathbf{p})}_{J(\mathbf... ...{J_{1}(\mathbf{f})}+\underbrace{-\log I\! \! P(\mathbf{f})}_{J_{2}(\mathbf{f})}$$ (8.3)

8.1.3 Reconstruction algébrique.

Dans l'expression Eq.9.3, pour trouver l'image $\mathbf{f}$ qui optimise le critère global $J(\mathbf{f})$, il est nécessaire de minimiser de façon conjointe deux termes $J_{1}(\mathbf{f})$ et $J_{2}(\mathbf{f})$. A toute image $\mathbf{f}$, le premier terme $J_{1}$ attribue la vraisemblance que la projection $\mathbf{p}$ soit issue de cette image. Ce critère s'appuie donc sur le lien pour passer de l'image à ses projections. C'est au niveau de ce terme que la modélisation du problème direct va intervenir. C'est un terme d'attache aux données qui donne un coût sur la façon dont on projette une image. Ce critère dépend donc:

• de la modélisation du problème direct (comment, partant de $\mathbf{f}$, je peux construire un jeu de projection $\hat{\mathbf{p}}$). En considérant un lien linéaire, nous avons $\hat{\mathbf{p}}=\mathbf{H}\mathbf{f}$.
• de la fonction de coût caractérisant la similarité de cette projection $\hat{\mathbf{p}}$ et du jeu de données réelles $\mathbf{p}$.

Le deuxième terme $J_{2}$ correspond à la probabilité d'observer une image $\mathbf{f}$ particulière. Tout image $\mathbf{f}$ prise au hasard parmi les $\Lambda ^{N}_{f}$ possibles est-elle susceptible d'être une image d'émission. Ce critère caractérise donc la connaissance a priori sur les images. C'est à ce niveau que tout information annexe sur le volume (indépendamment des projections) peut être introduite. Quelles sont, a priori, sans avoir connaissance des projections, les images les plus probables dans notre reconstruction? Si je ne sais absolument rien de l'image à obtenir, alors toutes les configurations possibles ont la même chance de survenir. Le terme $J_{2}$ qui correspond à cette équiprobabilité ne guide alors en rien dans le choix d'une image particulière. Le terme $J_{2}$ devenant constant, il ne nous reste plus que les données pour trancher dans le choix d'une estimée. Une telle approche correspond à choisir $\mathbf{f}$ au sens du maximum de vraisemblance (ML pour Maximum Likelihood).

Nous avons donc pour établir une reconstruction algébrique de multiples degrés de liberté. Le choix de l'algorithme de reconstruction est dépendant de la réponse à ces 3 questions:

1. Quel modèle vais-je choisir pour construire une projection à partir d'un volume émetteur quelconque ?
2. Comment vais-je traduire la similarité ($\approx$ distance) entre mes projections calculées et les données réelles ? Ceci pose la question, de manière déguisée, du modèle de bruit que nous adopterons !
3. Quel type d'information annexe vais-je introduire ?

La dernière question est une question piège, car il serait utopique de croire que nous pouvons nous en passer. En effet, quelle que soit la réponse aux deux premières questions (réponse qui sera toujours imparfaite), le système basé sur $J_{1}$ uniquement est mal posé. Cela veut dire qu'il n'existe pas d'image (=solution) $\mathbf{f}$ unique. La seule connaissance des données observées est insuffisante pour assurer l'existence, l'unicité et la stabilité d'une solution. Autrement dit, une reconstruction au sens de maximum de vraisemblance est illusoire. Il nous faut une contrainte supplémentaire qui nous guide dans le choix d'une estimée. Cette régularisation sera d'autant plus nécessaire que la réponse aux deux premières questions sera imparfaite.

L'avantage de la fonctionnelle (Eq.9.3) sur laquelle repose la reconstruction algébrique est qu'elle dissocie la reconstruction à proprement parler ($J_{1}$) de la régularisation ($J_{2}$). On peut donc s'attacher à décrire la fonctionnelle de régularisation $J_{2}$ d'abord, puis celle d'attache aux données et les algorithmes qui en découlent et qui furent envisagés. C'est pourquoi, nous allons d'abord décrire le choix d'une fonctionnelle $J_{2}$ qui s'appuie sur des principes qui furent d'abord énoncés en segmentation. Nous pourrons ensuite lui ajouter le terme d'attache aux données qui nous mènera directement aux algorithmes que nous avons envisagés dans nos travaux.