Autrefois, on passait des radiographies. Maintenant, on va aussi passer un examen par scanner: la technique s’appelle la tomodensitométrie axiale. Dans les deux cas, ce sont des rayons X qui traversent le corps du patient. Quelle est la différence entre les deux? Et pourquoi le scanner a-t-il révolutionné l’imagerie médicale?
Lorsque les rayons X traversent le corps humain, ils perdent de leur intensité. Pour une radiographie classique, les rayons X sortant du corps du patient impressionnent une plaque photographique: le ton de gris de la photo dépend de l’énergie du rayon qui est absorbée par le corps.
Wilhelm Rüntgen
1845-1923
Les rayons X
Les rayons X font partie de la grande famille des rayonnements électromagnétiques, au même titre que la lumière ou les ondes radio. Ils sont constitués, comme la lumière, de photons. Ce qui les en distingue c’est leur très courte longueur d’onde, laquelle varie entre 0,01 et 10 nanomètres pour les rayons X, alors que les longueurs d’onde de la lumière visible varient entre 390 nanomètres pour le violet, et 780 nanomètres pour le rouge, et que les longueurs d’onde des ondes radio sont beaucoup plus longues. La découverte des rayons X par Wilhelm Rüntgen lui valut le prix Nobel de physique en 1901. Il tire quatre grandes conclusions, dont deux sont essentielles en imagerie médicale:
– Les rayons X sont absorbés par la matière selon la masse atomique des atomes.
– Les rayons X impressionnent une plaque photographique.
Voici maintenant un scanner.
La région du corps à inspecter est enfilée dans l’anneau. L’appareil produit des images de coupes transversales du patient.
Comment? Visiblement il n’y a pas de plaque photographique parallèle à la coupe… Considérons une coupe transversale dans un plan: les rayons X ont traversé le corps dans toutes les directions de ce plan. Une partie de l’énergie de chaque rayon a été absorbée lors de sa traversée du corps. C’est cette perte d’énergie que mesure l’appareil: il obtient toute une série de nombres: un nombre pour chaque rayon.
Regardons l’exemple A1. On voit une tache plus dense et un rayon. De l’énergie a été absorbée lorsque ce rayon a traversé la forme associée à cette tache.
Dans l’exemple B1, la densité est uniforme et deux fois moins grande que dans la tache de l’exemple précédent. La même quantité d’énergie a été absorbée le long de ce même rayon.
C’est aussi la même quantité d’énergie qui a été absorbée le long du rayon dans l’exemple C1 à gauche!
Et même en rajoutant des rayons passant par le centre comme en C2, on ne réussirait pas à distinguer ces trois exemples: les mesures enregistrées seraient identiques…
Mais permettons maintenant dans nos trois exemples des rayons horizontaux ne passant pas par le centre.
Maintenant, on voit que la densité n’est pas distribuée de la même façon dans les trois exemples.
En effet, on voit qu’aucune énergie n’est absorbée pour les rayons extrêmes du premier et du troisième exemple (soit A2 et C3). Aussi, dans C3, moins d’énergie est absorbée le long des rayons passant par le centre que pour certains rayons décentrés, comme le rayon bleu. Mais ces rayons ne permettent pas de distinguer l’exemple de la densité uniforme, soit B2, de notre nouvel exemple D1 ici à droite.
Par contre des rayons verticaux du côté droit les distingueraient.
Vous pourriez continuer à étudier d’autres exemples plus compliqués et vous convaincre que:
Plus on utilise de rayons, plus on obtient d’information sur la densité interne.
Johann Radon
1887-1956
C’est le mathématicien autrichien, Johann Radon qui a montré en 1917 que, si on utilise tous les rayons possibles, alors il n’y a pas de perte d’information et on peut reconstruire la distribution de la densité en chaque point. À l’époque il ne s’intéressait pas à l’imagerie médicale et traitait un problème purement mathématique, soit celui de reconstruire une fonction à l’aide de l’ensemble des projections sur toutes les droites du plan.
Son travail est un exemple du fait qu’on ne sait pas d’avance, en général, quels problèmes mathématiques théoriques permettront une percée technologique importante.
Bien sûr, si l’on voulait considérer tous les rayons possibles, il y en aurait un nombre infini. On se contente donc en pratique de regarder un nombre fini de rayons, tous situés dans un même plan et formant un réseau assez dense, comme sur la figure ci-dessous. Pour chacun de ces rayons, on mesure l’énergie absorbée lorsque le rayon traverse le corps.
Le défi est de reconstruire une image à partir de tous ces nombres…
Mais une image de quoi ? L’énergie absorbée en un point du corps est proportionnelle à l’intensité du rayon et au coefficient d’atténuation en ce point (voir l’encadré « Les données recueillies par le scanner » ci-contre). Ce qu’on cartographie est ce coefficient d’atténuation. Sur nos dessins, un ton de gris plus élevé correspondait à une densité plus grande. Dans les images médicales obtenues par scanner, c’est le contraire : plus le coefficient d’atténuation est élevé, plus le ton de gris associé est pâle.
Les données recueillies par le scanner
Considérons un rayon D paramétré par s. En chaque point représenté par s la quantité d’énergie \(\Delta I\) absorbée est proportionnelle au coefficient d’atténuation A(s),
\[\Delta I = – A(s) I(s) \Delta s.\]
Faisons tendre \Delta s\) vers 0. On obtient une équation différentielle à variables séparables:
\[- \frac{dI}{I(s)} =A(s)ds.\]
I(s)
Appelons \(I_e\) l’intensité du rayon à son entrée dans le corps, et \(I_s\) son intensité à la sortie. Intégrons des deux côtés
\[\int_{I_e}^{I_s}-\frac{dI}{I(s)} = \int_D A(s)ds.\]
On obtient
\[-\ln I_s + \ln I_e = \int_D A(s)ds.\]
On voit donc que si on connaît \(I_e\) et que l’on mesure \(I_s,\) alors l’appareil peut calculer
\[\int_D A(s)ds.\]
Le coefficient d’atténuation augmente avec la densité des tissus. On peut l’augmenter artificiellement dans certains organes en faisant absorber au patient des liquides radioactifs : on améliore alors la qualité de l’image obtenue.
Les mathématiques du processus de récupération du coefficient d’atténuation sont avancées. Nous les décrivons brièvement dans l’encadré « De la transformée de Radon à la transformée de Fourier ». Par contre, la modélisation elle-même est plus accessible, et nous prenons le temps de la décrire en détails.
Mettre l’information sous forme de fonction.
Nous allons regarder le cas d’une image correspondant à une coupe dans le plan.
Nous voulons récupérer, à partir des données fournies par le scanner, le coefficient d’atténuation en chaque point : ceci est une fonction, f(x, y), inconnue. Quelle information nous donne le scanner? Il nous donne, pour chaque droite D du plan, la quantité d’énergie absorbée par le corps lors du passage d’un rayon X le long de cette droite. À chaque droite D nous associons donc un nombre. Nous reconnaissons le concept de fonction…
Par contre le domaine de la fonction est insolite…
C’est l’ensemble des droites du plan. En fait, pas toutes les droites, seulement celles, qui coupent le disque de l’image, mais il suffit de dire que la fonction vaut 0 pour les autres droites.
On sait manipuler des fonctions de plusieurs variables. On va donc décider de coder l’ensemble des droites d’un plan par des variables.
Regardons la droite rouge sur la figure.
Elle est uniquement déterminée par le vecteur bleu, \(\vec{v},\) joignant l’origine à sa projection sur la droite. Donc, on a une bijection entre l’ensemble des droites et l’ensemble des vecteurs issus de l’origine! Et on sait comment paramétrer de tels vecteurs: il suffit d’utiliser les coordonnées de l’extrémité du vecteur. Dans le contexte de la tomodensitométrie axiale, on choisit d’utiliser les coordonnées polaires r et \(\theta\) pour paramétrer le vecteur \(\vec{v}.\) Ce que l’appareil mesure, ce sont les valeurs d’une fonction \(F (r, \theta).\)
A-t-on avancé?
On cherche une fonction inconnue, la fonction d’atténuation, \(f (x, y).\) On connaît la fonction \(F (r, \theta),\) où \(F (r, \theta)\) représente la «somme» de \(f(x,y)\) le long du rayon paramétré par \(r\) et \(\theta.\)
Pour mieux comprendre la structure de toutes ces variables et informations, on va pousser un peu plus loin l’abstraction. Ce qu’on a construit, c’est une fonction \(\mathcal{R}.\)
\[f \mapsto F = \mathcal{R}(f),\]
dont le domaine et l’image sont des ensembles de fonctions! Cette fonction \(R\) porte un nom: c’est la transformée de Radon, du nom de Johann Radon qui l’a introduite en 1917.
La transformée de Radon
Appelons \(D_r,\theta\) la droite paramétrée par \((r, \theta).\) La droite passe par le point \((r \cos \theta,r \sin \theta)\) et un vecteur directeur est donné par \((-\sin \theta, \cos \theta).\) La droite est donc l’ensemble des points
\[\{(r \cos \theta – s \sin \theta, r \sin \theta + s \cos \theta), s \in \mathbb{R} \}.\]
Le scanneur mesure la transformée de Radon, soit
\[F(r,\theta)= \int_{D_{r, \theta}} f(x,y)ds\]
c’est-à-dire
\[\int_{-\infty}^{\infty} f(r \cos \theta-s \sin \theta,r \sin \theta+s \cos \theta)ds.\]
Supposons que la fonction \(\mathcal{R}\) soit inversible, et appelons \(\mathcal{R}^{-1}\) son inverse.
Pour récupérer \(f\) de \(F,\) il suffit de calculer \(\mathcal{R}^{-1}(F).\) Effectivement, un inverse de la transformée de Radon existera.
Essayons une idée. On doit récupérer le coefficient d’atténuation, \(f (x, y)\) au point \((x, y).\) Ce point a contribué aux valeurs de \(F (r, \theta)\) pour toutes les droites paramétrées par \(r\) et \(\theta\) et passant par le point \((x,y).\) Quelles sont ces droites?
Regardons la figure. Si une droite passant par \((x, y)\) est paramétrée par le vecteur \(\vec{v},\) alors le produit scalaire du vecteur\( (x,y)\) avec le vecteur unitaire \(\vec{e}. =(\cos \theta, \sin \theta) \)dans la direction du vecteur \(\vec{v}.\) vaut précisément \(r!\) (Voir Section problèmes).
Donc, l’ensemble de ces droites correspond à l’ensemble des solutions \((r, \theta)\) de
\[x \cos \theta + y \sin \theta= r.\]
On connaît maintenant toutes les valeurs de \(F(r, \theta)\) auxquelles \(f(x,y)\) a contribué.
On pourrait essayer de faire la moyenne de ces valeurs pour récupérer \(f(x,y).\) L’idée n’est pas mauvaise mais un peu trop naïve.
Revenons sur deux de nos exemples et intéressons-nous à \(f(0,0).\) Toutes les droites passant par (0,0) ont pour vecteur associé un vecteur \(\vec{v}\) de longueur \(r = 0.\) Dans le premier, A1, on a un disque uniforme de rayon 1/2.
Donc, tout rayon qui passe par l’origine a traversé une zone dense sur un segment de longueur 1. On a donc \(F (0, \theta) =1.\)
Dans l’illustration C1, la forme de l’image est un anneau de rayon intérieur 1/4 et de rayon extérieur 3/4. Ici encore, tout rayon qui passe par l’origine a traversé une zone dense sur un segment de longueur 1. On a donc encore \(F (0, \theta) =1.\)
On voit donc qu’on ne peut espérer récupérer \(f (0, 0)\) à partir de la seule connaissance des \(F (0, \theta) =1.\)
Pourquoi?
On n’a pas assez mélangé les différents \(f (x, y)\) lorsqu’on a calculé \(F (r, \theta)\)…
Il va falloir les mélanger encore plus pour obtenir une fonction \(G(X,Y),\) définie sur le plan, dont la valeur en \((X, Y)\) est une moyenne pondérée de tous les \(f(x,y).\) La pondération dépend de \((X, Y)\) (voir encadré « De la transformée de Radon à la transformée de Fourier »).
La transformation
\[f \mapsto F = \mathcal{R}(f) \mapsto G = \mathcal{F}(f),\]
sera inversible. Elle est connue des mathématiciens: c’est la transformée de Fourier en 2 variables. En appliquant \(\mathcal{F}^{-1}\) à la fonction \(G,\) on récupérera la fonction de départ \(f,\) et on pourra générer notre image. Cette partie devient trop technique, et nous ne rentrerons pas dans les détails.
Même si une formule existe pour inverser la transformée de Radon, il reste des défis mathématiques. En voici deux. Le scanner utilise seulement un échantillon fini de droites: quel est le meilleur choix, compte tenu des contraintes technologiques de fonctionnement de l’appareil? Du point de vue mathématique, comment faire pour que certaines petites erreurs d’approximation ne produisent pas de grosses erreurs lors de la reconstruction mathématique de l’image? On décide souvent d’appliquer un filtre avant de reconstruire l’image pour obtenir une image plus nette et plus fidèle (voir fin de l’encadré « De la transformée de Radon à la transformée de Fourier »).
De la transformée de Radon à la transformée de Fourier
Au coefficient d’atténuation, \(f,\) qui est une fonction de \((x,y),\) on a fait correspondre sa transformée de Radon, \(F,\) qui est une fonction de \((r, \theta).\) On fait correspondre à \(F\) deux fonctions \(G_1\) et \(G_2\) définies par
\[G_1(\rho,\theta)= \int_{-\infty}^{\infty} F(r,\theta) \sin r \rho dr\]
et
\[G_2(\rho,\theta)= \int_{-\infty}^{\infty} F(r,\theta) \cos r \rho dr. \]
Pour une direction \(\theta\) donnée, on peut avoir des oscillations en \(r.\) Ces fonctions viennent mesurer l’amplitude de l’oscillation de période \(2 \pi \rho.\) Elles nous disent en effet que \(F (r, \theta)\) contient une composante
\[\frac{1}{2 \pi} (G_1(\rho, \theta) \sin \rho r+ G_2(\rho, \theta) \cos \rho r). \]
La « somme » de toutes ces composantes pour toutes les valeurs de \(\rho\) (donnée par une intégrale) nous donne \(F (r, \theta).\)
Les erreurs d’approximation dans les composantes pour \(\rho\) grand conduisent souvent à beaucoup de flou dans l’image.
On peut choisir d’appliquer un filtre qui annule ces composantes pour obtenir une image plus nette, par exemple en posant \(G_1(\rho,\theta)=0\) et \(G_2(\rho,\theta)=0\) lorsque \(\rho\) dépasse un certain seuil.
