11.3.4 Méthode de rétroprojection des projections filtrées

La méthode de rétroprojection des projections filtrées se déduit de l'observation suivante: l'image f_B($\vec{{x}}\,$) dans la méthode du filtrage de la rétroprojection est non-nulle en-dehors du support de f ($\vec{{x}}\,$), et donc diffère de f ($\vec{{x}}\,$), uniquement parce que les mesures g($\theta$, s) sont strictement positives. Avant d'effectuer l'opération de rétroprojection, il faut effectuer un rehaussement des hautes fréquences dans les projections afin d'y introduire des composantes négatives. De cette façon, on peut espérer une certaine compensation de valeurs lors de l'étape de rétroprojection, pour obtenir une image qui est plus proche de f ($\vec{{x}}\,$), avec des valeurs nulles en-dehors du support de f ($\vec{{x}}\,$). La figure 11.14 illustre la situation.

Figure 11.14: Rétroprojection des projections filtrées.

 

On utilise l'expression mathématique

$$\vec{{x}}\, \int_{{0}}^{{\pi}} \theta \theta \vec{{x}}\, \vec{{\theta}}\,$$ (11.13)

pour représenter la reconstruction obtenue en utilisant la méthode de rétroprojection des projections filtrées. Cette expression fait intervenir un filtre de convolution 1D, h(s), qu'il faut appliquer aux données g($\theta$, s) à $\theta$ fixé pour trouver f_h($\vec{{x}}\,$) proche de f ($\vec{{x}}\,$). Pour des raisons de symétrie, on impose à h(s) d'être une fonction paire.

Le filtre h(s) s'obtient à partir de la relation entre f_h($\vec{{x}}\,$) et f ($\vec{{x}}\,$). Suivant les mêmes étapes de calcul qu'à la section précédente, on trouve en effet que

$$\vec{{x}}\, \int_{{-\infty}}^{{+\infty}} \int_{{-\infty}}^{{+\infty}} \pi \vec{{x}}\, \vec{{X}}\, \mathcal {F} \vec{{X}}\, {\frac{{(\mathcal{F}_{1}h)(\vert\vec{X}\vert)}}{{\vert\vec{X}\vert}}}$$ (11.14)

À partir de cette relation, on voit qu'il suffit de choisir h(s) tel que ($\mathcal {F}$₁h)($\sigma$) = |$\sigma$| pour obtenir une reconstruction exacte. Notons que ce filtre de convolution est uniquement défini au sens des distributions.

En pratique, la reconstruction exacte de f ($\vec{{x}}\,$) n'est jamais possible parce que les projections ne sont connues qu'avec un certain pas d'échantillonnage $\Delta$s. Le mieux que l'on puisse obtenir est une reconstruction avec résolution $\Delta$x = $\Delta$y = $\Delta$s et le filtre h(s) doit être adapté pour prendre cet aspect en compte. Cela conduit au choix suivant pour h(s):

$$\mathcal {F} \sigma \left\vert\vphantom{\sigma}\right. \sigma \left.\vphantom{\sigma}\right\vert \left(\vphantom{\frac{\vert\sigma\vert}{\sigma_{max}}}\right. {\frac{{\vert\sigma\vert}}{{\sigma_{max}}}} \left.\vphantom{\frac{\vert\sigma\vert}{\sigma_{max}}}\right)$$ (11.15)

où $\sigma_{{max}}^{}$ = 1/(2$\Delta$s) est défini par la fréquence de coupure de NYQUIST des projections et W(t) est par exemple choisi tel que

$$\left\{\vphantom{ \begin{array}{ccc} 1 & & si\,\left\vert t\right\vert<1\\ 0 & & si\,\left\vert t\right\vert\geq1\end{array}}\right. \begin{array}{ccc} 1 & & si\,\left\vert t\right\vert<1\\ 0 & & si\,\left\vert t\right\vert\geq1\end{array}$$ (11.16)

(On parle de fenêtre de coupure rectangulaire). Avec une telle définition de h(s), on obtient ce que l'on appelle une version à bande limitée de f ($\vec{{x}}\,$), à savoir

$$\vec{{x}}\, \int_{{-\infty}}^{{+\infty}} \int_{{-\infty}}^{{+\infty}} \pi \vec{{x}}\, \vec{{X}}\, \mathcal {F} \vec{{X}}\, \left(\vphantom{\frac{\vert\vec{X}\vert}{\sigma_{max}}}\right. {\frac{{\vert\vec{X}\vert}}{{\sigma_{max}}}} \left.\vphantom{\frac{\vert\vec{X}\vert}{\sigma_{max}}}\right)$$ (11.17)

Il s'agit d'une version de f ($\vec{{x}}\,$) où tous les détails correspondants aux fréquences supérieures à $\sigma_{{max}}^{}$ en module sont supprimés. Un aspect intéressant de ce type de reconstruction est que l'on peut espérer une résolution uniforme au sein de f_h($\vec{{x}}\,$) car f_h($\vec{{x}}\,$) est relié à f ($\vec{{x}}\,$) par une opération de convolution.

Finalement, notons que l'on peut choisir différentes expressions pour la fonction W(t), suivant le niveau de bruit dans les mesures. Lorsque les données sont bruitées, on modifie généralement l'expression de W(t) de façon à réduire la contribution des hautes fréquences et donc la propagation du bruit au sein de l'image. Les expressions typiquement choisies sont la fenêtre sinus cardinal (sinc),

$$\left\{\vphantom{ \begin{array}{ccc} \frac{\sin(\pi t/2)}{\pi t/2} & & si\,\left\vert t\right\vert<1\\ 0 & & sinon\end{array}}\right. \begin{array}{ccc} \frac{\sin(\pi t/2)}{\pi t/2} & & si\,\left\vert t\right\vert<1\\ 0 & & sinon\end{array}$$ (11.18)

et la fenêtre de HANNING

$$\left\{\vphantom{ \begin{array}{ccc} \cos^{2}(\pi t/2) & & si\,\left\vert t\right\vert<1\\ 0 & & sinon\end{array}}\right. \begin{array}{ccc} \cos^{2}(\pi t/2) & & si\,\left\vert t\right\vert<1\\ 0 & & sinon\end{array}$$ (11.19)

Comparée à la fenêtre de HANNING, la fenêtre SINC propage plus le bruit au sein de l'image mais permet aussi d'obtenir des images avec meilleure résolution.

Figure 11.15: Quelques fenêtres de coupures.

 

Figure 11.16: Quelques résultats de reconstruction.

 

La figure 11.15 montre l'allure de ($\mathcal {F}$₁h)($\sigma$)/$\sigma_{{max}}^{}$ pour les différents choix de W(t). La figure 11.16 montre des résultats de reconstruction obtenus pour le fantôme de SHEPP-LOGAN sur une grille de 256×256 pixels avec un sinogramme de taille 400×256, en utilisant successivement la fenêtre rectangulaire, la fenêtre sinc, puis la fenêtre de HANNING.