next up previous contents index
Next: 7.4 Ce qu'il faut Up: 7. Reconstruction et Rebinning. Previous: 7.2 Dans l'espace de   Contents   Index

Subsections

7.3 Dans l'espace de Fourier.


7.3.1 FOurier REbinning (FORE).

Le rebinning par l'algorithme FORE proposé dès 1995 par Defrise et collaborateurs [28] est une approche attrayante, déjà intellectuellement par la théorie sur laquelle il s'appuie, mais aussi parcequ'il permet de s'affranchir des problèmes de précisions (SSRB) ou de mauvais comportement en présence de bruit (MSRB). En revanche, l'algorithme FORE ne repose pas sur des concepts aussi simples que les autres méthodes de rebinning. Nous l'envisageons seulement dans une approche géométrique, mais une expression exacte du rebinning développée par Defrise [29] permet d'établir le coeur de l'algorithme et ses limitations. L'approche théorique s'offre même le luxe d'englober le SSRB.

7.3.1.1 Approche géométrique.

Nous avons montré que dans le cas idéal, le sinogramme revenait à intégrer les photons $ \gamma $ sur la LOR $ A_{i}B_{j} $. Nous avons montré que le plan de projection $ \Pi $ est invariant par translation suivant le vecteur $ \vec{n} $. Le point $ P $ de projection peut donc être choisi n'importe où sur la LOR. Par la suite ce point est choisi à mi-chemin entre les deux détecteurs $ A_{i} $ et $ B_{j} $. Ceci conduit à des coordonnées faciles à calculer.

Lorsque, partant de la géométrie illustrée Fig.8.2,

Figure: Géométrie pour la description de l'algorithme FORE.
[Une LOR $ A_iB_i $.] \resizebox*{0,45\columnwidth}{!}{\includegraphics{imgps/rre_fig15.ps}}

[Notation dans le plan transaxial passant par $ P $.] \resizebox*{0,45\textwidth}{!}{\psfrag{theta}[][][2]{\textcolor{blue}{\( \theta ...
...ag{n}[][][2]{\textcolor{red}{\( \vec{n} \)}}\includegraphics{imgps/rre_fig3.ps}} [Notation dans le plan contenant axe $ z $ et le point $ P $.] \resizebox*{0,45\textwidth}{!}{\psfrag{theta}[][][2]{\textcolor{blue}{\( \theta ...
...ag{n}[][][2]{\textcolor{red}{\( \vec{n} \)}}\includegraphics{imgps/rre_fig4.ps}}

nous cherchons à lier les coordonnées du point $ P $ avec les données auxquelles nous avons accès dans un sinogramme (le quadruplet $ (r,s,\theta ,\phi ) $), nous trouvons :

$\displaystyle \overrightarrow{OP}\left\{ \begin{array}{c}
r\sin \theta \\
-r\cos \theta \\
\frac{s}{\cos \phi }
\end{array}\right. $

7.3.1.1.1 Expression de la transformée de Fourier des projections.

Exprimons maintenant le vecteur directeur de la LOR:

$\displaystyle \vec{n}\left\{ \begin{array}{c}
\cos \theta \cos \phi \\
\sin \theta \cos \phi \\
\sin \phi
\end{array}\right. $

Comme $ \phi $ est constant et différent de $ \frac{\pi }{2} $, nous réécrivons ces coordonnées différemment:

$\displaystyle \vec{n}\left\{ \begin{array}{c}
\frac{1}{\cos \phi }\cos \theta \...
...{\cos \phi }\sin \theta \\
\frac{1}{\cos \phi }\tan \phi
\end{array}\right. $

Si on traduit l'intégration sur la LOR, nous obtenons:
$\displaystyle p(r,s,\theta ,\phi )$ $\displaystyle =$ $\displaystyle \int _{t_{1}}f(r\sin \theta +\frac{t_{1}}{\cos \phi }\cos \theta ...
...hi }\sin \theta ,\frac{s}{\cos \phi }+\frac{t_{1}}{\cos \phi }\tan \phi )dt_{1}$  
  $\displaystyle =$ $\displaystyle \cos \phi .\int _{t}f(r\sin \theta +t\cos \theta ,-r\cos \theta +t\sin \theta ,s+t\tan \phi )dt$ (7.1)

D'un autre coté, exprimons la transformée de Fourier 2D de $ p(r,s,\theta ,\phi ) $ suivant les variables $ r $ et $ \theta $.

$\displaystyle P(\nu _{r},s,\nu _{\theta },\phi )=\int ^{2\pi }_{0}\int ^{\infty...
...ty }p(r,s,\theta ,\phi )e^{-2\pi i(r.\nu _{r}+\theta .\nu _{\theta })}drd\theta$ (7.2)

En remplaçant dans Eq.8.2, $ p(r,s,\theta ,\phi ) $ par sa valeur dans Eq.8.1, nous obtenons:

$\displaystyle P(\nu _{r},s,\nu _{\theta },\phi )=\cos \phi .\int _{\theta }\int...
...\theta ,s+t\tan \phi )e^{-2\pi i(r.\nu _{r}+\theta .\nu _{\theta })}dtdrd\theta$ (7.3)

On réécrit tout cela en considérant les variables $ x $ et $ y $ et on choisit $ x=\psi _{1}(r,t) $ et $ y=\psi _{2}(r,t) $ de telle sorte que:

$\displaystyle \left\{ \begin{array}{c}
x=r\sin \theta +t\cos \theta \\
y=-r\c...
...sin \theta -y\cos \theta \\
t=x\cos \theta +y\sin \theta
\end{array}\right. $

Comme nous reconnaissons le changement de variable associé à une rotation, nous en déduisons que le déterminant de la matrice Jacobienne $ \vert\mathcal{J}\vert $ de ce changement de variable vaut 1 ( $ \vert\mathcal{J}\vert=1 $). Sachant que:

$\displaystyle \int _{x}\int _{y}f(x,y)dxdy=\int _{r}\int _{t}f(\psi _{1}(r,t),\psi _{2}(r,t))\underbrace{\vert\mathcal{J}\vert}_{=1}drdt$

Nous en déduisons pour Eq.8.3 en permutant l'ordre des intégrations,

$\displaystyle P(\nu _{r},s,\nu _{\theta },\phi )=\cos \phi .\int \int _{x,y}\un...
... -y\cos \theta ).\nu _{r}+\theta .\nu _{\theta })}d\theta \right] }_{I_{1}}dxdy$ (7.4)

Regardons de plus près ce que représente l'intégration sur les $ \theta $.

7.3.1.1.2 Théorème de phase stationnaire

Nous pouvons considérer l'intégrale $ I_{1}$ (Eq.8.4) comme une somme vectorielle de complexes $ \bar{z}(\theta )=A(\theta )e^{-iB(\theta )}=A(\theta )[cos(B(\theta ))+i.sin(B(\theta ))] $ (Fig.8.3.a). Notons que le changement de signe de $ A(\theta )$ fait basculer la phase. Considérons alors les cas suivants:

Figure 8.3: Principe de phase stationnaire.
[représentation du complexe $ \bar{z}(\theta)$.] [Cas où $ A(\theta )$ et $ B(\theta )$ varient rapidement.] [Cas où seul $ B(\theta )$ varie.] [Cas où $ A(\theta )$ et $ B(\theta )$ varient lentement.]
\resizebox*{0,2\textwidth}{!}{\psfrag{A(theta)}[br][][2]{\( A(\theta ) \)} \psfrag{B(theta)}[][][2]{\( B(\theta ) \)}\includegraphics{imgps/rre_fig8.ps}} \resizebox*{0,2\textwidth}{!}{\psfrag{A(theta)}[br][][2]{\( A(\theta ) \)} \psfrag{B(theta)}[][][2]{\( B(\theta ) \)}\includegraphics{imgps/rre_fig9.ps}} \resizebox*{0,2\textwidth}{!}{\psfrag{A(theta)}[br][][2]{\( A(\theta ) \)} \psfrag{B(theta)}[][][2]{\( B(\theta ) \)}\includegraphics{imgps/rre_fig10.ps}} \resizebox*{0,2\textwidth}{!}{\psfrag{A(theta)}[br][][2]{\( A(\theta ) \)} \psfrag{B(theta)}[][][2]{\( B(\theta ) \)}\includegraphics{imgps/rre_fig11.ps}}

Ainsi, en faisant l'hypothèse que la fonction $ f $ est relativement lisse, donc qu'elle varie peu, du moins par rapport à la variation de la phase, nous pouvons dire que les contributions à l'intégrale $ I_{1}$ proviennent principalement des endroits où la phase est stationnaire. Cette propriété correspond justement au théorème de phase stationnaire [26]. Il s'agit donc de trouver les $ \theta $ qui rendent cette phase stationnaire. Pour cela, on cherche angles de vues pour lesquels $ \frac{\delta B(\theta )}{\delta \theta }=0 $.

$\displaystyle \frac{\delta [(x\sin \theta -y\cos \theta )\nu _{r}+\theta \nu _{\theta }}{\delta \theta }=(x\cos \theta +y\sin \theta )\nu _{r}+\nu _{\theta }=0$

Notons que lorsque $ \vert\nu _{\theta }\vert\leq \vert\nu _{r}\vert\sqrt{x^{2}+y^{2}} $, il y a deux solutions dans l'intervalle $ [0,2\pi ] $. Nous avons alors:

$\displaystyle -\frac{\nu _{\theta }}{\nu _{r}}=x\cos \theta +y\sin \theta =t$

Cette propriété, constatée par Edholm, est remarquable. Sur la figure Fig.8.4.a, nous avons représenté la variation de la phase $ B(\theta )$ et de sa dérivée $ \frac{dB(\theta )}{d\theta } $. Nous allons regarder la contribution à l'intégrale $ I_{1}$ des $ \theta $ appartenant soit à $ [3.6,3.8] $, soit à $ [2.1,2.3]. $ Le premier intervalle correspond à des valeurs de $ \theta $ qui rendent la dérivée de la phase maximale. Le deuxième, à des $ \theta $ qui donnent une variation de phase proche de $ 0 $. On constate dans le premier cas illustré Fig.8.4.b. que les vecteurs $ \bar{z}(\theta)$ échantillonnent régulièrement le cercle unité. On retrouve le cas décrit Fig.8.3.c. La contribution de cet intervalle est donc négligeable. En revanche, sur le deuxième intervalle (phase stationnaire), les vecteurs $ \bar{z}(\theta)$ sont confinés sur une petite portion du cercle unité (Fig.8.4.c), on retrouve donc le cas décrit Fig.8.3.d. La contribution de cet intervalle va donc être prépondérante.

Figure 8.4: Phase stationnaire pour l'algorithme Fore.
[Expression de $ B(\theta )$ (vert) et de sa dérivée (bleu) quand $ \theta \in [0,2\pi [ $] [Vecteur $ \bar{z}(\theta)$ lorsque $ \theta \in {[} 3.6,3.8 {]} $] [Vecteur $ \bar{z}(\theta)$ lorsque $ \theta \in {[} 2.1,2.3 {]} $]

Ainsi, les 2 points qui contribuent dans l'espace de Fourier au sinogramme sont situés sur la LOR à la distance $ t=-\frac{\nu _{\theta }}{\nu _{r}} $. En pratique, le sinogramme intègre l'information sur de multiples LOR.

Cette relation de distance fréquentielle établit que toute valeur de $ P $ de fréquence $ (\nu _{r},\nu _{\theta }) $ reçoit principalement des contributions de points situés à une distance $ t=-\frac{\nu _{\theta }}{\nu _{r}} $ sur les LOR.

Cette approximation peut en quelque sorte être vue comme un temps de vol virtuel [29].

En appliquant ce théorème de phase stationnaire à notre équation (Eq.8.4), i.e. en remplaçant le troisième argument de $ f $ et en sortant ce terme de l'intégrale nous obtenons:

$\displaystyle P(\nu _{r},s,\nu _{\theta },\phi )\approx \cos \phi .\int \int _{...
...sin \theta -y\cos \theta )\nu _{r}+\nu _{\theta }.\theta )}d\theta \right] dxdy$ (7.5)

Comme ce troisième argument ne dépend ni de $ x $, ni de $ y $, il n'est donc pas affecté par le changement de variable que nous avons effectué. Autrement dit, si je calcule $ P(\nu _{r},s-\frac{\nu _{\theta }}{\nu _{r}}\tan \phi ,\nu _{\theta },0) $, j'obtiens exactement la même expression que 8.5, c'est à dire la clé de voûte de l'algorithme FORE:

$\displaystyle P(\nu _{r},s,\nu _{\theta },\phi )\approx \cos \phi .P(\nu _{r},s-\frac{\nu _{\theta }}{\nu _{r}}\tan \phi ,\nu _{\theta },0)$ (7.6)

Au vu des angles $ \phi $ utilisés, on peut en première approximation considérer le cosinus de cet angle comme égal à l'unité. En conséquence, nous pouvons légèrement simplifier l'Eq.8.6 qui devient:

$\displaystyle P(\nu _{r},s,\nu _{\theta },\phi )\approx P(\nu _{r},s-\frac{\nu _{\theta }}{\nu _{r}}\tan \phi ,\nu _{\theta },0)$ (7.7)

Cette expression correspond vraiment au coeur de l'algorithme FORE puisqu'elle nous permet de passer d'une acquisition 3D fournissant un sinogramme à 4 dimensions (membre de gauche de l'Eq.8.7) en un sinogramme fictif qui proviendrait d'une acquisition 2D ne comptant donc plus qu'une seule inclinaison (la variable $ \phi $ est constamment nulle dans le membre de droite de l'Eq.8.7). Ce sinogramme fictif intègre donc bien les données issues des différents segments dans des plans transaxiaux. Ces données sont intégrées à une hauteur particulière ( $ s-\frac{\nu _{\theta }}{\nu _{r}}\tan \phi $) dans le sinogramme fictif. Nous pouvons donc directement en déduire un algorithme.

7.3.1.2 Implémentation.

Elle se traduit par l'utilisation de l'Eq.8.7 dans laquelle nous cherchons à construire le sinogramme fictif $ p_{2D} $ correspondant à une acquisition 2D. Lorsqu'il décrit le principe de cet algorithme, Defrise montre qu'il est nécessaire de considérer différemment certaines régions de l'espace de Fourier associées à $ \theta $ et $ r $. Nous renvoyons le lecteur à l'article fondateur pour ne conserver que le cas général [29]. L'algorithme correspondant à l'Eq.8.7 est donné Alg.2.
\begin{algorithm}
% latex2html id marker 11791\begin{enumerate}
\item Initiali...
... ) \).
\end{enumerate}\par\caption{FORE
}\selectlanguage{french}
\end{algorithm}


7.3.2 FOurier Simple Averaging (FOSA).

7.3.2.1 Principe.

En reconstruction, du fait du théorème de coupe-tranche, nous savons que les projections échantillonnent de manière particulière l'espace de Fourier de l'objet à reconstruire. Ainsi, dès 1981, Stark [87] propose en 2D une reconstruction qui s'appuie sur ce théorème. L'idée est de prendre la transformation de Fourier 2D des projections $ (r,s,\theta ,\phi ) $, d'affecter par interpolation la valeur de chaque fréquence $ (\nu _{r},\nu _{s},\theta ,\phi ) $ ainsi obtenue dans la grille de l'espace de Fourier de l'objet $ (\nu _{x},\nu _{y},\nu _{z}) $ et de faire la transformée de Fourier 3D inverse pour retrouver l'objet [89,90]. L'inconvénient majeur de cette méthode vient de l'interpolation nécessaire pour passer des projections à l'espace de Fourier de l'objet. On considère en général, une grille cartésienne pour tirer parti de la vitesse d'exécution des FFT. L'interpolation des valeurs sur cette grille conduit à des distorsions de l'image. On trouve donc différentes méthodes d'interpolation pour passer des projections à la grille cartésienne: plus proches voisins, bilinéaires, sinus cardinal [75,88,87,94]. D'autres proposent des méthodes par transformation du système de coordonnées [19]. C'est pourquoi, en nous arrêtant au rebinning nous évitons ce passage à une grille cartésienne. Autrement dit, on va rebinner les données dans l'espace de Fourier des projections en affectant les fréquences $ (\nu _{r},\nu _{s},\theta ,\phi \neq 0) $ au niveau des fréquences $ (\nu _{r},\nu _{s},\theta ,\phi =0) $ issues d'un sinogramme ne présentant qu'un seul segment. Les angles $ \phi $ restant petits l'interpolation se fait entre 2 grilles voisines.

7.3.2.2 Implémentation.

Sur la Fig.8.5.a, les plans de projection $ \Pi _{1} $, $ \Pi _{2} $,$ \Pi _{3} $ récoltent les photons pour 3 angles de vues différents sous l'inclinaison $ \phi =0 $. Il s'agit de plans de notre premier segment. Le plan $ \Pi _{4} $ intègre les photons sous une inclinaison $ \phi \neq 0 $.

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

[Dans l'espace de l'objet.] \resizebox*{0,45\textwidth}{!}{\psfrag{Cercle}[][][2]{\( \Omega \)} \psfrag{x}[]...
...\)} \psfrag{PiF4}[][][2]{\( \Pi ^{F}_{4} \)}\includegraphics{imgps/rre_fig5.ps}} [Dans l'espace de Fourier.] \resizebox*{0,45\textwidth}{!}{\psfrag{Cercle}[][][2]{\( \Omega \)} \psfrag{x}[]...
...\)} \psfrag{PiF4}[][][2]{\( \Pi ^{F}_{4} \)}\includegraphics{imgps/rre_fig6.ps}}

C'est l'information de ce plan que nous allons intégrer dans les plans du premier segment. Pour cela, prenons la transformée de Fourier 2D suivant les variables $ r $ et $ s$ de ces plans de projections. Nous obtenons les transformées de nos 4 plans de projection $ \Pi ^{F}_{1} $, $ \Pi ^{F}_{2} $, $ \Pi ^{F}_{3} $, $ \Pi ^{F}_{4} $. Du fait du théorème de section centrale, ces derniers échantillonnent l'espace de Fourier comme sur la Fig.8.5.b. Considérons un point particulier $ P(_{P}\nu _{r},_{P}\nu _{s},\theta _{4},\phi _{4}) $ de la transformée de Fourier 2D $ \Pi ^{F}_{4} $ d'un plan de détecteurs incliné. De manière générale, ce point est ``coincé'' entre les transformées de deux plans du premier segment. Sur notre figure, il s'agit des plans $ \Pi ^{F}_{2} $ et $ \Pi ^{F}_{3} $. Ces plans sont échantillonnés. Le point $ P $ est donc situé dans un héxaèdre formé par les 8 échantillons fréquentiels les plus proches de $ \Pi ^{F}_{2} $ et $ \Pi ^{F}_{3} $ (Fig.8.6).

Figure 8.6: Interpolation dans l'espace de Fourier.
\resizebox*{0,3\textwidth}{!}{\psfrag{PiF2}[][][2]{\( \Pi ^{F}_{2} \)} \psfrag{PiF3}[][][2]{\( \Pi ^{F}_{3} \)}\includegraphics{imgps/rre_fig7.ps}}

On envisage donc l'algorithme (Alg.3) suivant pour rebinner les données.
\begin{algorithm}
% latex2html id marker 11980\begin{enumerate}
\item Initiali...
... ) \).
\end{enumerate}\par\caption{FOSA
}\selectlanguage{french}
\end{algorithm}


next up previous contents index
Next: 7.4 Ce qu'il faut Up: 7. Reconstruction et Rebinning. Previous: 7.2 Dans l'espace de   Contents   Index
Lecomte Jean François 2002-09-07