A simple Tomography program in Javascript

Author : Nicolas Lebbe -

Image originale à scanner Reconstruction après le scan
Nombre rayons parallèles :   Nombre angles : Résolution x :   Résolution y :
\(\renewcommand{\rar}[1]{\overrightarrow{#1}} \newcommand{\phi}{\varphi} \newcommand{\fourier}[1]{\mathcal{F}\left(#1\right)} \newcommand{\ifourier}[1]{\mathcal{F}^{-1}\left(#1\right)} \renewcommand{\Re}{\mathbb{R}} \definecolor{lightblue}{rgb}{.8, .8, 1} \newcommand{\highlight}[1]{\fcolorbox{black}{lightblue}{$\displaystyle #1$}}\)

Explications techniques :


keywords : tomography, Beer-Lambert law, Radon transform, FBP, Fourier transform, javascript.

notation : \(\fourier{f}(\xi)\) la transformée de Fourier de \(f\) en \(\xi\).

Sommaire

1. Introduction
2. La transformée de Radon
3. Récupération de \(\mu\)
a. The projection slice theorem
b. Filtered back projection
4. Intégration numérique

1. Introduction

La loi de Beer-Lambert donne la relation entre l'intensité de lumière en entrée \(I_0\) et en sortie \(I_D\) d'un rayon qui a traversé durant une longueur \(l\) un matériau de coefficient d'absorption \(\alpha\) : $$ \begin{eqnarray*} I_D & = & I_0 e^{-\alpha l}. \end{eqnarray*}$$ Supposons que l'on envoie maintenant un rayon d'intensité \(I_0\) connue en direction d'un détecteur. Le rayon passe à travers différents matériaux de longueurs \((l_i)_{i=1..n}\) et de coefficients d'absorptions \((\mu_i)_{i=1..n}\).
En appliquant alors la loi de Beer-Lambert \(n\) fois on trouve la relation : $$ \begin{eqnarray*} I_D & = & I_0 e^{-\sum_{i=1}^n \mu_i l_i}. \end{eqnarray*}$$ Si le rayon traverse maintenant un matériau de coefficient d'absorption variable dans l'espace \(\mu(\rar{x})\) avec \(\rar{x} \in B(0, 1)\), en faisant tendre le nombre de matériaux \(n\) vers l'infini dans l'expression précédente, on trouve le résultat suivant pour un rayon suivant la ligne \(L\) : $$ \begin{eqnarray} \highlight{ I_D = I_0 e^{-\int_{L} \mu(\rar{x}) d\rar{x}} } \end{eqnarray}$$ Le principe de la tomographie est d'envoyé un très grand nombre de rayons selon différentes lignes \(L_i\) afin d'obtenir différentes intensités \(I_D(L_i)\).
Connaissant ces intensités, on utilise une méthode de « projection inverse » (FBP) qui permet de trouver le coefficient d'absorption \(\mu(\rar{x})\) en tout point \(\rar{x}\).

2. La transformée de Radon

Avec \(\varphi \in [0, 2\pi[\) et \(s \in [-1, 1]\) on a $$ \begin{eqnarray*} \rar{\theta}(\phi) = \left( \begin{array}{c} \cos{\phi}\\ \sin{\phi} \end{array} \right) & \text{et} & \rar{\zeta}(\phi) = \left( \begin{array}{c} -\sin{\phi}\\ \cos{\phi} \end{array} \right). \end{eqnarray*} $$ D'ou la définition suivante d'une ligne (rayon) pour un angle \(\phi\) à une distance \(s\) du centre : $$ \begin{eqnarray*} L_{\phi, s} & = & \{\rar{x} \in \Re^2, \rar{x} \cdot \rar{\theta}(\phi) = s\}\\ & = & \{\rar{x} = s \rar{\theta}(\phi) + l \rar{\zeta}(\phi),\ \forall\ l \in \Re\}\\ & = & s \rar{\theta}(\phi) + \Re \rar{\zeta}(\phi). \end{eqnarray*} $$ La transformée de Radon d'une fonction \(\mu \in L^1(\Re^2)\) noté \(R_\mu\) est alors défini par : $$ \begin{eqnarray*} \highlight{ R_\mu(\phi, s) = \int_{L_{\phi, s}} \mu(\rar{x}) d\rar{x} } \end{eqnarray*} $$ Aussi notée \( R_\mu(\phi, s) = \displaystyle\int_{\Re^2} \mu(\rar{x}) \delta_0(\rar{x} \cdot \rar{\theta}(\phi) - s) d\rar{x}\).
(Avec \(\delta_0\) la fonction dirac en 0)
On voit donc qu'en envoyant un rayon (d'angle \(\phi\) et de distance au centre \(s\)) on récupère l'intensité \(I_D(\phi, s)\) qui, connaisant \(I_0\), permet de récupérer la valeur de la transformée de Radon de \(\mu\) en \((\phi, s)\) : $$ \begin{eqnarray*} \highlight{ -\ln\left(\frac{I_D(\phi, s)}{I_0}\right) = R_\mu(\phi, s) } \end{eqnarray*} $$

3. Récupération de \(\mu\)

On cherche maintenant à obtenir une formule donnant \(\mu \in L^1(\Re^2)\) à partir des données \(R_\mu(\phi,s)\) récupérées.

a. The projection slice theorem

$$ \begin{eqnarray*} \forall \phi \in [0, 2\pi], \sigma \in \Re, \quad \fourier{R_\mu(\phi, \cdot)}(\sigma) & = & \fourier{\mu}(\sigma \rar{\theta}(\phi)) \end{eqnarray*}$$
Démonstration $$ \begin{eqnarray*}\begin{array}{rcll} \fourier{R_\mu(\phi, \cdot)}(\sigma) & = &\displaystyle \int_{\Re} R_\mu(\phi, s)e^{-2i\pi \sigma s} ds & \text{par transformée de Fourier}\\ & = &\displaystyle \int_{\Re} \int_{\Re} \mu(s\rar{\theta}(\phi) + l\rar{\zeta}(\phi)) dl e^{-2i\pi \sigma s} ds & \text{par définition de } L_{\phi,s}\\ & = &\displaystyle \int_{\Re} \int_{\Re} \mu(s\rar{\theta}(\phi) + l\rar{\zeta}(\phi)) e^{-2i\pi \sigma s} ds dl & \text{avec le théorème de Fubini}\\ & = &\displaystyle \int_{\Re^2} \mu(\rar{x}) e^{-2i\pi \sigma \rar{x}\cdot\rar{\theta}(\phi)} ds dl & \text{chdvr. } \rar{x} = s\rar{\theta}(\phi) + l\rar{\zeta}(\phi)\\ & = &\displaystyle \fourier{\mu}(\sigma \rar{\theta}(\phi)) \end{array}\end{eqnarray*}$$

b. Filtered back projection

$$ \begin{eqnarray*} \highlight{ \mu(\rar{x}) = \int_0^\pi \int_{\Re} R_\mu(\phi, t) r(\rar{x}\cdot\rar{\theta}(\phi) - t) dt d\phi } \end{eqnarray*}$$
Démonstration $$ \begin{eqnarray*} \begin{array}{rcll} \mu(\rar{x}) & = &\displaystyle \int_{\Re^2} \fourier{\mu}(\xi)e^{2i\pi \rar{x}\cdot\rar{\xi}}d\rar{\xi} & \text{par transformée de Fourier}\\ & = &\displaystyle \int_0^\pi \int_\Re \fourier{\mu}(\sigma \rar{\theta}(\phi))e^{2i\pi\sigma\rar{x}\cdot\rar{\theta}(\phi)}|\sigma|d\sigma d\phi & \text{chdvr. polaire } \rar{\xi} = \sigma \rar{\theta}(\phi)% = (\sigma \cos(\phi), \sigma \sin(\phi))\\ \\ & = &\displaystyle \int_0^\pi \int_\Re \fourier{R_\mu(\phi, \cdot)}(\sigma)e^{2i\pi\sigma\rar{x}\cdot\rar{\theta}(\phi)}|\sigma|d\sigma d\phi & \text{d'après le théorème précédent.}\\ & \simeq &\displaystyle \int_0^\pi \int_\Re \fourier{R_\mu(\phi, \cdot)}(\sigma) \times |\sigma|\chi_{[-f,f]}(\sigma) e^{2i\pi\sigma\rar{x}\cdot\rar{\theta}(\phi)}d\sigma d\phi & \text{avec $f$ une haute fréquence} \end{array}\end{eqnarray*}$$ Supposons maintenant l'existence d'une fonction \(r\) dont la transformée de Fourier soit telle que \(\fourier{r}(\sigma) = |\sigma|\chi_{[-f,f]}(\sigma)\).
Posons désormais \(g_\phi(\sigma) = \fourier{R_\mu(\phi, \cdot)}(\sigma) \times \fourier{r}(\sigma)\) de telle sorte que : $$\begin{eqnarray*} \mu(\rar{x}) & = & \int_0^\pi \int_\Re g_\phi(\sigma)e^{2i\pi\sigma\rar{x}\cdot\rar{\theta}(\phi)}d\sigma d\phi = \int_0^\pi \ifourier{g_\phi}(\rar{x}\cdot\rar{\theta}(\phi))d\phi. \end{eqnarray*}$$ En utilisant le résultat suivant sur la transformée de Fourier on peut réécrire la fonction \(g_\phi\) : $$\begin{eqnarray*} \forall f,g \in L^1(\Re), \quad \fourier{f}\times\fourier{g} & = & \fourier{f \ast g} \quad \text{ avec }\ast\text{ le produit de convolution.} \end{eqnarray*}$$ $$\begin{eqnarray*} \Rightarrow g_\phi & = & \fourier{R_\mu(\phi, \cdot) \ast r}\\ & = & \fourier{ \int_{\Re} R_\mu(\phi, t) r(\cdot - t) dt} \end{eqnarray*}$$ Donc finalement \(\ifourier{g_\phi}\) se simplifie et on obtient : $$\begin{eqnarray*} \mu(\rar{x}) & = & \int_0^\pi \int_{\Re} R_\mu(\phi, t) r(\rar{x}\cdot\rar{\theta}(\phi) - t) dt d\phi \end{eqnarray*}$$ Qui est le résultat attendu. Il est à noter que l'intégrale sur \(\Re\) se simplifie en une intégrale sur \([-1,1]\) car c'est le support de la fonction \(R_\mu(\phi,\cdot)\).

La fonction \(r\) est obtenue par transformée de Fourier inverse et on obtient $$\begin{eqnarray*} r(s) & = & \int_{-f}^f |t|e^{2 i \pi s t} dt = \frac{2 \pi s f \sin(2 \pi s f) + \cos(2 \pi s f) - 1}{2\pi^2 s^2}. \end{eqnarray*}$$

4. Intégration numérique

En pratique on ne peut envoyer qu'un nombre fini de rayons.
On va donc pouvoir obtenir la transformée de Radon de \(\mu\) pour seulement \(n_\phi\) angles \((\phi_i)_{i=0..(n_\phi-1)}\) et \(n_s\) distances au centre \((s_j)_{j=0..(n_s-1)}\).
On a donc pour données la matrice \((R_\mu(\phi_i,s_j))_{\substack{i=0..(n_\phi-1)\\ j=0..(n_s-1)}}\).
Afin de pouvoir utiliser le thèorème FBP donné dans la section précédente, il va donc falloir réaliser une intégration numérique. Pour ce faire, le plus simple est d'utiliser des angles et distances au centre à pas constant. Autrement dit on aura \(\phi_i = \frac{i}{n_\phi}\pi\) et \(s_j = -1 + 2\frac{j}{n_s-1}\).
On utilise alors une méthode de quadrature du type suivant : $$\begin{eqnarray*} \int_a^b f(x) dx & \simeq & \frac{b-a}{n} \left( \sum_{k=0}^{n} f\left(a+k\frac{b-a}{n}\right) \right) \end{eqnarray*}$$ Ce qui donne alors dans notre cas en \(\rar{x} = (x,y)\) : $$\begin{eqnarray*}\highlight{ \mu(x,y) \simeq \frac{2\pi}{n_s n_\phi} \sum_{i=0}^{n_\phi} \sum_{j=0}^{n_s} R_\mu(\phi_i, s_j) r(x\cos(\phi_i) + y\sin(\phi_i) - s_j) }\end{eqnarray*}$$ Deux choses intéressantes à remarquer :