6.1 Reconstruction d'images : solution générale.

La solution analytique du problème d'inversion, permettant de retrouver une image 2D $f(x,y)$ à partir de ses projections $p(r,\theta )$, est un problème connu et résolu depuis plus de 80 ans par un mathématicien autrichien (Radon 1917). Si on prend le temps d'en rappeler les rudiments ici, c'est que le cas bidimensionnel contient en lui même les fondements du problème plus complexe que constitue la reconstruction 3D. En effet, les opérations de filtrage, de rétroprojection, le théorème de la section centrale sont des composantes présentes aussi bien en 2D qu'en 3D. La reconstruction 2D représentant une version allégée du cas 3D, nous traiterons ultérieurement les problèmes inhérents au seul cas du 3D.

6.1.1 Reconstruction analytique 2D : Une solution unique.

6.1.1.1 La projection.

Dans la première partie de cette thèse, nous avons, de manière approfondie, cherché à modéliser le lien qui unit l'objet radioactif à ses projections. Lorsqu'on se place dans un contexte continu, nous avons montré que chaque élément du sinogramme $p(r,\theta )$ comptabilisait les nombres de paires de photons $\gamma$ émises par l'objet radioactif $f(x,y)$ suivant la ligne de coïncidence référencée par $(r,\theta )$. Il s'agissait d'une projection. En termes mathématiques, on considère que $p(r,\theta )$ intègre l'information relative à cet objet sur une droite (Fig.7.1):

$$p(r,\theta )=\int \int _{x,y}f(x,y)\delta (r-x\cos \theta -y\sin \theta )dxdy$$

où $\delta$, distribution de Dirac traduit l'intégration sur la droite de projection.

Figure: Modélisation de la projection analytique en 2D.

Nous sommes donc capable de construire un sinogramme partant d'une distribution de radioactivité. Comme en pratique, les imageurs TEP ne nous donnent accès qu'au sinogramme, il faut bien, pour obtenir des images, résoudre le problème inverse. La reconstruction consiste donc partant d'un sinogramme $p(r,\theta )$ à estimer la distribution $f(x,y)$. Il devient donc indispensable de définir la rétroprojection.

6.1.1.2 La rétroprojection.

L'opération de rétroprojection correspond en fait à un épandage. La rétroprojection élémentaire permet de construire à partir d'une projection suivant un angle $\theta$ une fonction $e_{\theta }(x,y)$. Cette fonction est obtenue en affectant à chaque point $(x,y)$ le nombre de paires de photons émises $p(r,\theta )$. $r$ indique l'abscisse où se projette ce point suivant une direction $\theta$ ( $r=x\cos \theta +y\sin \theta$). La rétroprojection élémentaire conduit donc à:

$$e_{\theta }(x,y)=p(r,\theta )\vert _{r=x\cos \theta +y\sin \theta }$$

Ce processus est illustré Fig.7.2.

Figure: Illustration de la rétroprojection en 2D.

Cette opération élémentaire nous permet de définir la rétroprojection, opération qui consiste à sommer l'ensemble de ces rétroprojections élémentaires pour toutes les directions de projections possibles:

$$e(x,y)=\int ^{\pi }_{0}e_{\theta }(x,y)d\theta =\int ^{\pi }_{0}p(r,\theta )\vert _{r=x\cos \theta +y\sin \theta }d\theta$$

Cette sommation nous conduit à une image $e(x,y)$, image d'épandage, qui même si elle n'est pas loin de l'image que l'on cherche (à une convolution près) n'est pas notre distribution finale. Il nous manque un élément: Le théorème de section-centrale (coupe-tranche).

6.1.1.3 Théorème de section centrale.

Ce théorème, dont la preuve peut être trouver dans [66], traduit un lien entre la transformée de Fourier 1D( $P(\nu _{r},\theta )=TF_{1D}[p(r,\theta )]$) des projections $p(r,\theta )$ et la transformée de Fourier 2D ( $F(\nu _{x},\nu _{y})=TF_{2D}[f(x,y)]$) de l'objet $f(x,y)$.

Figure: Théorème de section centrale en 2D.

[Espace de l'objet.]	[Espace des fréquences.]

Ce théorème (Fig.7.3) stipule que la transformée de Fourier 1D des projections échantillonne suivant une direction définie par $\theta$ l'espace de Fourier de l'objet. Si $\tilde{F}(\rho ,\theta )=F(\rho \cos \theta ,\rho \sin \theta )$ traduit la transformée de Fourier 2D de l'objet exprimée en coordonnées polaires, alors le théorème de section centrale se traduit mathématiquement par:

$$P(\nu _{r},\theta )=\tilde{F}(\rho ,\theta )\vert _{\rho =\nu _{r}}$$

Nous disposons maintenant de toutes les définitions nécessaires à la reconstruction de notre image $f(x,y)$.

6.1.1.4 La reconstruction.

Pour reconstruire notre image partant de sa transformée de Fourier exprimée en coordonnées polaires, nous utilisons la définition de la transformée de Fourier:

$$f(x,y)=\int ^{2\pi }_{0}\int ^{+\infty }_{0}\rho \tilde{F}(\rho ,\theta )e^{2i\pi (\rho \cos \theta .x+\rho \sin \theta .y)}d\rho d\theta$$

En divisant cette intégrale en deux:

$$f(x,y)= \int ^{\pi }_{0}\int ^{+\infty }_{0}\rho \tilde{F}(\rho ,\theta )e^{2i\pi (\rho \cos \theta .x+\rho \sin \theta .y)}d\rho d\theta$$

$$+ \int ^{2\pi }_{\pi }\int ^{+\infty }_{0}\rho \tilde{F}(\rho ,\theta ')e^{2i\pi (\rho \cos \theta '.x+\rho \sin \theta '.y)}d\rho d\theta '$$ (6.1)

Si dans le deuxième terme de cette expression, je change $\theta '$ par $\theta +\pi$, j'obtiens pour celle-ci:

$$\int ^{\pi }_{0}\int ^{+\infty }_{0}\rho \tilde{F}(\rho ,\theta +\pi )e^{2i\pi (\rho \cos (\theta +\pi )x+\rho \sin (\theta +\pi )y)}d\rho d\theta$$

Pour des raisons de symétrie évidente (Fig.7.3), nous avons $\tilde{F}(\rho ,\theta )=\tilde{F}(-\rho ,\theta +\pi )$, et comme $\cos (\theta +\pi )=-\cos \theta$ et $\sin (\theta +\pi )=-\sin \theta$, nous obtenons pour ce second terme:

$$\int ^{\pi }_{0}\int ^{+\infty }_{0}\rho \tilde{F}(-\rho ,\theta )e^{2i\pi (-\rho )(\cos \theta .x+\sin \theta .y)}d\rho d\theta$$

$$= \int ^{\pi }_{0}\int ^{0}_{-\infty }\vert\rho \vert\tilde{F}(\rho ,\theta )e^{2i\pi \rho (\cos \theta .x+\sin \theta .y)}d\rho d\theta$$

En rassemblant maintenant les deux termes de 7.1, on obtient pour l'expression de la transformée de Fourier:

$$f(x,y)=\int ^{\pi }_{0}\int ^{+\infty }_{-\infty }\vert\rho \vert... ...(\rho ,\theta )e^{2i\pi (\rho \cos \theta .x+\rho \sin \theta .y)}d\rho d\theta$$ (6.2)

De cette expression, par application du théorème de section centrale et en mettant en évidence certaines portions, on arrive finalement:

$$f(x,y) = \int ^{\pi }_{0}\underbrace{\int ^{+\infty }_{-\infty }\vert\nu _... ..._{r}}_{p_{f}(r,\theta )=TF_{1D}(\vert\nu _{r}\vert P(\nu _{r},\theta ))}d\theta$$

$$= \int ^{\pi }_{0}p_{f}(r,\theta )d\theta$$ (6.3)

L'image $f(x,y)$ s'obtient donc comme la rétroprojection, non pas des projections simples, mais des projections filtrées $p_{f}(r,\theta )$.(Eq.7.3). Le filtre est simple, et nous avons directement son expression dans le domaine fréquentiel: il s'agit du filtre monodimensionnel rampe $h(\nu _{r})=\vert\nu _{r}\vert$. L'algorithme de reconstruction comporte donc deux étapes:

• Filtrage monodimensionnel des projections dans l'espace de Fourier $P_{f}(\nu _{r},\theta )=P(\nu _{r},\theta )\times h(\nu _{r)}$.
• Rétroprojection des projections filtrées

Il porte donc le nom de rétroprojection des données filtrées ou filtered-backprojection (FBP).

D'un point de vue totalement théorique, les opérateurs utilisés étant tous linéaires, il n'y a pas de contradiction à rétroprojeter les projections et à filtrer l'image épandue $e(x,y)$ ensuite par un filtre rampe bidimensionnel $H(\nu _{x},\nu _{y})=\sqrt{\nu ^{2}_{x}+\nu ^{2}_{y}}$, nous avons alors:

$$F(\nu _{x},\nu _{y})=H(\nu _{x},\nu _{y})\times E(\nu _{x},\nu _{y})$$

Cet algorithme porte le nom de rétroprojection filtrée ou backprojection-filtering (BPF). En pratique, on filtre les projections avant de les rétroprojeter car la rétroprojection est ainsi limitée par le volume d'intérêt [84].

N'en déplaise aux puristes de la langue française, il faut noter une faiblesse de la langue française pour traduire la différence entre ces deux algorithmes de manière simple, c'est pourquoi les acronymes utilisés (FBP et BPF) seront construits sur les noms des algorithmes exprimés en langue anglaise.

6.1.1.5 Filtre de reconstruction et filtre d'apodisation.

Nous venons de voir que pour reconstruire notre image, nous étions amenés à utiliser le filtre rampe $h(\nu _{r})=\vert\nu _{r}\vert$. L'utilisation de ce filtre rampe pose deux problèmes majeurs:

• Le support de la fonction représentant le filtre rampe n'est pas borné, ce qui pose évidemment des problèmes dans un milieu discret.
• Cette fonction n'est pas majorée.

Afin de remédier à ces deux problèmes, l'idée consiste à conserver la forme de la rampe fréquentielle pour les petites fréquences et de l'atténuer pour les fréquences plus élevées. Il s'agit de choisir judicieusement la fréquence seuil $\nu _{c}$ à partir de laquelle le filtre atténuera la rampe. Il faut garder à l'esprit qu'en pratique les calculs sont effectués sur des mesures discrètes. Par conséquent un spectre trop étendu provoquera des phénomènes de repliement de spectres. Pour éviter cela nousdevons tronquer les composantes spectrales situé au dessus de la fréquence de Shannon. D'autre part, une troncature brutale dans l'espace de Fourier conduit à des ondulations dans l'espace de l'objet. C'est pourquoi généralement, on utilise le filtre $h(\nu _{r})=\vert\nu _{r}\vert\times A(\nu _{r})$ où $A$ est une fonction d'apodisation. Il existe plusieurs fonctions d'apodisation possibles. Celle utilisé en routine sur la caméra ECAT HR+ à CYCERON est la fenêtre de Hanning:

$$\left\{ \begin{array}{ll} A(\nu _{r})=\frac{1}{2}\left( 1+\cos \f... ...\, \vert\nu _{r}\vert<\nu _{c}\\ A(\nu _{r})=0 & ailleurs \end{array}\right.$$

D'un point de vue physique, la fenêtre d'apodisation $A(\nu _{r})$ permet de ne pas trop amplifier les hautes fréquences et donc le bruit dans les images reconstruites.

6.1.2 Reconstruction analytique 3D.

Dans le cas de la reconstruction d'un objet $f(x,y)$ bidimensionnel, les projections étaient bidimensionnelles $p(r,\theta )$. Le système était donc totalement défini analytiquement et la solution de notre système était unique. En 3D cependant, l'intégration sur les lignes de coïncidence est spécifiée par 4 paramètres $(r,s,\theta ,\phi )$. Or le volume que l'on doit reconstruire ne comporte que trois variables $(x,y,z)$, le problème d'inversion est donc surdéterminé. En effet, en se limitant aux projections situées dans des plans perpendiculaires à l'axe du cylindre ($\phi =0$), nous disposons de suffisamment d'information pour reconstruire l'objet, du moins en faisant une reconstruction 2D tranche par tranche (i.e. en faisant varier $s$). D'un point de vue purement analytique, l'information supplémentaire apportée par les acquisitions suivant les autres inclinaisons ( $\phi \neq 0$) est superflue. Enfin, cela reste vrai tant qu'on ne fait pas intervenir de considérations statistiques ! De toute façon, cette surdétermination nous conduit à une solution qui n'est plus unique. Idéalement, toutes les solutions conduisent à la même estimation du volume. Malheureusement dans le monde réel, c'est-à-dire en présence de bruit, la solution varie sensiblement. Suivant une approche proche de celle de Defrise [30], nous allons voir les éléments qui constituent les fondements de la reconstruction analytique 3D.

6.1.2.1 Ensemble $\Omega$ des projections.

Tout plan de projection $\Pi$ est défini par une direction de projection $\vec{n}$, donc par deux angles $(\theta ,\phi )$. Ce vecteur $\vec{n}$ étant un vecteur unitaire, sa norme est égale à 1. Si on prend pour représentant de ce vecteur, le vecteur $\overrightarrow{ON}=\vec{n}$ dont une extrémité $O$ est située sur l'origine du repère et si on s'intéresse au lieu géométrique que dessine l'autre extrémité $N$ lorsque $(\theta ,\phi )$ varient dans $[0,\pi ]\times [-\frac{\pi }{2},\frac{\pi }{2}]$, on constate que ce point $N$ décrit une sphère de rayon unité. Nous appellerons cette sphère, la sphère des projections.

Or l'ensemble des valeurs que peuvent prendre les deux angles $\theta$ et $\phi$ est lié à la géométrie d'acquisition. En effet, ces angles sont fixés par les couples de détecteurs. Par conséquent, le lieu géométrique que dessine le vecteur $\vec{n}$ est fonction de la géométrie. C'est en général un sous ensemble de la sphère des projections. C'est ce sous ensemble que nous désignons par $\Omega$, qui constitue l'ensemble des projections. La figure Fig.7.4

Figure 7.4: Exemples d'ensemble de projections $\Omega$.

[$\Omega$ quelconque.]	[$\Omega$ relatif à une acquisition 2D.]
[$\Omega$ relatif à une acquisition sphérique (Full 3D).]	[$\Omega$ relatif à une acquisition cylindrique.]

donne quelques exemples de cet ensemble $\Omega$. Il faut noter que l'ensemble des projections correspondant au cas Full 3D conduit à une géométrie d'acquisition qui serait sphérique. Ce qui est impossible techniquement. On note aussi que l'acquisition 2D est un cas limite, où l'ensemble $\Omega$ se réduit à un cercle équatorial (cercle dont le plan support passe par le centre de la sphère de projection).

6.1.2.2 Extension du théorème de section centrale.

Comme précédemment le théorème de section centrale traduit un lien, dans l'espace de Fourier, entre les projections et l'objet radioactif. Un plan de l'espace de Fourier 3D de l'objet ( $F(\nu _{x},\nu _{y},\nu _{z})=F(\vec{\nu })=TF_{3D}[f(x,y,z)]$) et perpendiculaire à la direction de projection $\vec{n}$ correspond à la transformée de Fourier 2D des projections ( $P(\nu _{r},\nu _{s},\theta ,\phi )=P(\vec{\nu }_{\bot },\vec{n})=TF_{2D}[p(r,s,\theta ,\phi )]$) pour cette même direction de projection (Fig.7.5).

Figure 7.5: Passage des plans de projections dans l'espace de Fourier de l'objet

[Dans l'espace de l'objet.]	[Dans l'espace de Fourier.]

$$F(\vec{\nu })=P(\vec{\nu }_{\bot },\vec{n})\vert _{\vec{\nu }.\vec{n}=0}$$ (6.4)

6.1.2.2.1 Remarque:

Nous pouvons toujours décomposer le vecteur $\vec{\nu }$, comme la somme de deux vecteurs $\vec{\nu }_{\bot }$ et $\vec{\nu }_{//}$. Le premier $\vec{\nu }_{\bot }$ représente un vecteur situé dans le plan perpendiculaire au vecteur de projection $\vec{n}$ (plan $\Pi$ de projection) . Le deuxième $\vec{\nu }_{//}$ est colinéaire à la direction de projection $\vec{n}$. Nous avons:

$$\vec{\nu } = \vec{\nu }_{r}+\vec{\nu }_{//}$$

$$= \vec{\nu }_{\bot }+(\vec{\nu }.\vec{n}).\vec{n}$$

Dans ce cas, il est clair que si $\vec{\nu }.\vec{n}=0$, le vecteur $\vec{\nu }$ est un vecteur situé dans le plan perpendiculaire à la direction de projection. De la même manière, nous pouvons définir dans l'espace de l'objet:

$$\vec{x}=\vec{x}_{\bot }+(\vec{x}.\vec{n}).\vec{n}$$

Partant de $\vec{x}$, nous avons accès aux coordonnées de la projection dans le plan perpendiculaire à $\vec{n}$ par $\vec{x}_{\bot }=\vec{x}-(\vec{x}.\vec{n}).\vec{n}$.

6.1.2.3 Ensemble $\Omega ^{\bot }$.

Dans l'Eq. 7.4, il est nécessaire de bien comprendre l'expression:

$$\vec{\nu }.\vec{n}=0$$ (6.5)

et de voir ce qu'elle représente. Dans l'espace de Fourier de l'objet, dont nous définissons le repère $(O,\vec{\nu }_{x},\vec{\nu }_{y},\vec{\nu }_{z})$, si nous fixons une direction $\vec{n}=\overrightarrow{ON}$, le produit scalaire précédent définit un ensemble de fréquences qui doivent être situées dans le plan perpendiculaire à la direction $\vec{n}$ et passant par l'origine. Inversement en fixant une fréquence $\vec{\nu }$ , on devrait construire un plan perpendiculaire à $\vec{\nu }$ qui contient les directions de projection $\vec{n}$ vérifiant Eq.7.5. Or, d'après le paragraphe Par.7.1.2.1, le lieu géométrique suivi par les directions de projection est également situé sur la sphère unité. De fait, le lieu géométrique défini par les directions vérifiant Eq.7.5 correspond à l'intersection d'un plan perpendiculaire à $\vec{\nu }$ passant par $O$ et de la sphère des projections.

Figure: Lieu géométrique correspondant à $\vec{\nu }.\vec{n}=0$.

Ce lieu géométrique $\Omega ^{\bot }(\vec{\nu })$ (Fig.7.6) correspond à un cercle équatorial .

$\Omega ^{\bot }(\vec{\nu })$ correspond donc à l'ensemble des directions de projection $\vec{n}$ pour lesquelles la fréquence $\vec{\nu }$ appartient au plan $\Pi$ de projection construit sur $\vec{n}$.

Nous avons vu également que la géométrie d'acquisition limitait les directions de projection à $\Omega$. Ainsi de manière générale, l'équation Eq.7.5 fournit un jeu de directions qui correspond à l'intersection de $\Omega$ et $\Omega ^{\bot }(\vec{\nu })$. Muni de ces deux ensembles, nous pouvons maintenant comprendre les conditions suffisantes qu'il faut imposer pour pouvoir reconstruire.

6.1.2.4 Les conditions d'Orlov. [70]

Du fait du théorème de section centrale, l'objet à imager $f(\vec{x})$ est déterminé, si dans l'espace de Fourier, nous disposons de l'information $F(\vec{\nu })$ pour toutes les fréquences $\vec{\nu }$. De par l'équation Eq.7.5, il faut pour chaque fréquence $\vec{\nu }$ que l'intersection entre $\Omega$ et $\Omega ^{\bot }(\vec{\nu })$ soit non nulle.

Il doit exister au moins une direction $\vec{n}$ qui construit dans l'espace de Fourier un plan contenant cette fréquence !

Nous venons de le voir, toute direction correspond à une cercle équatorial unique, Le théorème d'Orlov est donc le suivant:

La reconstruction de l'objet est possible si tout cercle équatorial à au moins une intersection non vide avec l'ensemble des projections possibles $\Omega$.

Si on reprend les exemples d'ensembles de projections illustrés sur la Fig.7.7,

Figure 7.7: Exemples d'ensemble de projections $\Omega$.

[$\Omega$ quelconque.]	[$\Omega$ relatif à une acquisition 2D.]
[$\Omega$ relatif à une acquisition sphérique (Full 3D).]	[$\Omega$ relatif à une acquisition cylindrique.]

on constate que le premier exemple Fig.7.7.a ne permet pas de reconstruction car on peut construire des cercles équatoriaux qui n'intersectent pas le domaine $\Omega$. En revanche, les autres cas autorisent tous la reconstruction de l'objet.

Arrêtons nous quelques instants sur le cas limite de l'acquisition 2D. Pour toute fréquence $\vec{\nu }$ ( $\vec{\nu }\neq \vec{\nu }_{z}$), le cercle équatorial $\Omega ^{\bot }$ intercepte l'ensemble des projections $\Omega$ (qui correspond lui aussi à un cercle équatorial) en 2 points. Cela revient à dire qu'il existe 2 directions de projection pour lesquelles le plan $\Pi$ contient la fréquence considérée. Ces deux directions sont en opposition (180$^{\circ }$ l'une de l'autre) et sont donc équivalentes. Elles correspondent à un même couple de détecteurs. Cela veut dire que nous avons, pour chaque fréquence, une seule direction de projection qui contribue à cette fréquence. Nous avons donc un rapport d'une direction pour une fréquence.

Dans le cas Full 3D, l'intersection de $\Omega ^{\bot }$ et $\Omega$ conduit toujours à $\Omega ^{\bot }$. Nous avons donc un rapport de $2\pi$ directions pour une fréquence. Ceci illustre bien la redondance des données lorsqu'on réalise une acquisition 3D.

Dans ces 2 cas (acquisition 2D ou Full 3D), le rapport du nombre de projection(s) par direction est constant quelle que soit la fréquence considérée. En revanche, pour une acquisition cylindrique, ce ratio n'est plus constant et dépend de la fréquence considérée. En effet, l'intersection de $\Omega ^{\bot }$ et $\Omega$ conduit à un arc de cercle, dont la longueur dépend de la fréquence. Or le ratio fréquence/coupure est justement lié à cette longueur. Ainsi, dans ce cas pratique fréquent, il est nécessaire de prendre en compte la variabilité de la redondance au moment de la reconstruction, i.e. l'intégrer à notre filtre de reconstruction.

6.1.2.5 Expression générale du filtre de reconstruction.

Nous considérons le cas de la reconstruction par FBP. Evidemment et comme son nom l'indique, nous avons d'abord une étape de filtrage suivie d'une rétroprojection. Le filtrage des données ne consiste plus uniquement en un simple filtre rampe car il doit intégrer la variabilité des redondances.

6.1.2.5.1 Filtrage.

Chaque projection est filtrée. $H(\vec{\nu }_{\bot },\vec{n})$ représente le noyau de filtrage, que nous spécifierons par la suite, exprimé dans l'espace de Fourier. Le filtrage se traduit juste par une multiplication:

$$P_{f}(\vec{\nu }_{\bot },\vec{n})=H(\vec{\nu }_{\bot },\vec{n})\times P(\vec{\nu }_{\bot },\vec{n})$$ (6.6)

Où la transformée de Fourier s'obtient par l'application de :

$$H(\vec{\nu }_{\bot },\vec{n})=\int \int _{R^{2}}h(\vec{x}_{\bot },\vec{n})e^{-2i\pi \vec{\nu }_{\bot }.\vec{x}_{\bot }}d^{2}\vec{\nu }_{\bot }$$

6.1.2.5.2 Rétroprojection.

Comme précédemment, la rétroprojection consiste à affecter à la position $\vec{x}$ la valeur du sinogramme en $\vec{x}_{\bot }$ (endroit où se projette ce point) pour l'ensemble $\Omega$ des directions.

$$e(\vec{x})=\int \int _{\Omega }p_{f}(\vec{x}_{\bot },\vec{n})d^{2}\vec{n}$$ (6.7)

où $\vec{x}_{\bot }$ représente la projection de $\vec{x}$ sur un plan de projection perpendiculaire à $\vec{n}$.

6.1.2.5.3 Noyau de filtrage $H$.

Le noyau de filtrage n'est pas unique. Pour le trouver, écrivons la fonction de transfert de l'ensemble du processus que représente la mesure des projections, leur filtrage et leur rétroprojection. Pour cela, réécrivons la rétroprojection Eq.7.7 dans l'espace de Fourier:

$$E(\vec{\nu })=\int \int _{\Omega }P_{f}(\vec{\nu },\vec{n})\delta... ...{2}\vec{n}=\int _{\Omega ^{\bot }(\vec{\nu })}P_{f}(\vec{\nu },\vec{n})d\vec{n}$$

où la fonction de Dirac $\delta$ nous permet de ne sélectionner dans l'espace de Fourier, uniquement les fréquences $\vec{\nu }$ perpendiculaire à la direction de projection $\vec{n}$. Pour une fréquence $\vec{\nu }$, cela nous conduit à une intégration curviligne sur $\Omega ^{\bot }(\vec{\nu })$. En utilisant l'expression du filtrage des projections Eq.7.6 et le théorème de section centrale, nous obtenons:

$$E(\vec{\nu }) = \int \int _{\Omega }H(\vec{\nu },\vec{n}).P(\vec{\nu },\vec{n}).\delta (\vec{\nu }.\vec{n})d^{2}\vec{n}$$

$$= \int \int _{\Omega }H(\vec{\nu },\vec{n}).F(\vec{\nu }).\delta (\vec{\nu }.\vec{n})d^{2}\vec{n}$$

$$= F(\vec{\nu }).\underbrace{\int \int _{\Omega }H(\vec{\nu },\vec{n}).\delta (\vec{\nu }.\vec{n})d^{2}\vec{n}}_{T(\vec{\nu })}$$

Dans l'espace de Fourier, pour que la rétroprojection des données filtrées nous permette de remonter à l'objet, il suffit de poser $T(\vec{\nu })=1$ pour toutes les fréquences $\vec{\nu }$. Le filtre de projection prend donc en compte le fait que les arcs $\Omega ^{\bot }(\vec{\nu })$ ne sont pas égaux pour toutes les fréquences. Ainsi, la solution générale conduit à toute fonction $G(\vec{\nu },\vec{n})$ normalisée par l'intégration de ces valeurs sur l'arc $\Omega ^{\bot }(\vec{\nu })$ [30]:

$$\forall \vec{\nu }\in \mathcal{R}^{3}\; \frac{G(\vec{\nu },\vec{n... ...nt \int _{\Omega }G(\vec{\nu },\vec{n})\delta (\vec{\nu }.\vec{n})d^{2}\vec{n}}$$ (6.8)

Choisir un filtre correspond donc à choisir une fonction $G$ qui vérifie la propriété Eq.7.8. Parmis ces filtres, on distingue les filtres $G(\vec{\nu },\vec{n})=_{1}H(\vec{\nu }).w(\vec{n})$ qui sont à variables séparables. Pour ces filtres factorisables, les opérations de filtrage et rétroprojection peuvent être permutées et le choix de la fonction $w(\vec{n})$, fonction positive est indifférent. Tous les filtres conduisent à la même reconstruction en l'absence de données bruitées. En présence de bruit, les filtres se comportent de manières différentes. Nous allons maintenant décrire le filtre que nous utiliserons par la suite dans nos reconstructions analytiques par rétroprojection des données filtrées.

6.1.2.6 Le filtre de Colsher.

Le filtre que nous utiliserons fut développé au départ par Colsher [21]. Ce filtre nécessite au départ un angle d'acceptance suivant la direction $z$ qui soit constant et limité. Nous avons vu Ch.4 ce que représentait l'angle d'acceptance axial. Cet angle n'est pas réellement constant dans notre cas, puisque suivant le segment envisagé, les plans pairs et les plans impairs n'intègrent pas l'information sur le même nombre de plans de détecteurs. D'autre part, du fait de la géométrie d'acquisition cylindrique, certaines projections sont incomplètes. Supposons néanmoins que l'angle d'acceptance est borné par $\phi _{max}$ de la sorte que toutes les directions de projections soient telles que l'angle $\phi$ associé à la direction de projection vérifie $\phi <\phi _{max}$. Dans ce cas, le filtre se met sous la forme:

$$H_{Colsher}(\vec{\nu },\vec{n})=_{1}H_{Colsher}(\vec{\nu }).\unde... ...u }}\vert}} & Si\, \vert\Psi _{\vec{\nu }}\vert>\phi _{max} \end{array}\right.$$

Dans cette expression $\Psi _{\vec{\nu }}$ désigne l'angle que fait le vecteur $\vec{\nu }$ avec l'axe porté par $\overrightarrow{\nu _{z}}$ dans l'espace des fréquences (Fig.7.6). On trouvera de plus amples renseignements sur la façon dont est obtenue ce filtre dans [21] et sur l'implémentation dans [20].

6.1.2.7 Vues manquantes et reconstruction.

Jusqu'à présent, nous avons supposé que les projections suivant une direction $\vec{n}$ sont totalement mesurée. Or, du fait de la géométrie cylindrique d'acquisition, certaines vues sont manquantes. Pour détourner ce problème, Kinahan et Rogers ont proposé un algorithme basé sur une reprojection d'une estimée initiale [51]. On la calcule en utilisant uniquement les projections du sinogramme correspondant à des directions transaxiales ($\phi =0$). L'ensemble de ces projections répond aux conditions d'Orlov. Partant de ces projections et par l'utilisation d'un algorithme de reconstruction par rétroprojection 2D des données filtrées, on reconstruit un premier volume (=estimée initiale). Ce volume est ensuite projeté afin de compléter les vues manquantes du sinogramme. Le sinogramme alors complet sert ensuite à la reconstruction 3D. Cette procédure est illustrée Fig.7.8.

Figure 7.8: Principe de la reconstruction 3D avec prise en compte des vues manquantes[4].

6.1.2.8 Algorithme de reconstruction.

Par la suite, lorsque nous parlerons de reconstruction analytique standard, nous entendrons la rétroprojection filtrée utilisant le filtre de Colsher. Nous nous placerons toujours juste après les différentes corrections, et l'algorithme utilisé pour obtenir une image partant de ce sinogramme corrigé est illustré Alg.1.