Même si l'ensemble des procédures ont été écrites dans un souci de rapidité, nous allons insister plus particulièrement sur 2 procédures: la rétroprojection et le calcul du produit matriciel impliqué dans la reconstruction algébrique. Ces deux procédures en effet, montrent comment on peut tirer parti de la puissance de calcul du CRAY à l'aide des simples notions que nous avons introduites précédemment.
La rétroprojection est un mécanisme essentiel dans la reconstruction d'images. Si on raisonne en 2D et si l'image est décrite sur une base de pixels, sa taille est de dimensions . Un algorithme conventionnel comme l'algorithme de Shepp et Logan requiert multiplications et additions. En passant en 3D, la complexité passe à . Ce surcoût de calcul reste encore à l'heure actuelle un handicap.9.4 C'est pourquoi, il s'avère impératif de réduire le temps de calcul nécessaire à la rétroprojection. Les tentatives faites en ce sens peuvent être regroupées en trois catégories:
Concernant la deuxième catégorie, des solutions spécifiques ont été mises en place. Par exemple, un circuit spécialement conçu pour l'algorithme STRETCH a été utilisé par [96]. Plus récemment, un circuit permettant d'effectuer une interpolation bilinéaire nécessaire à la reconstruction 3D a été proposé [47]. Les temps de rétroprojection sont alors inférieurs à la minute.
Finalement le recours au parallélisme est une idée séduisante. D'une part par les fermes de PCs qui voient le jour actuellement. Et plus particulièrement au centre CYCERON puisqu'elle peut être facilement adaptée de manière performante à la machine dont nous disposons.
Rappelons que la projection bidimensionnelle analytique nous avait amené à écrire:
Inversement, la rétroprojection consiste à affecter en chaque point d'une image la valeur de toutes les projections pour tous les couples tels que Ceci correspond à une étape de localisation qui consiste à chercher quels sont les couples ou ``faisceaux'' qui voient ce pixel précis.
L'interpolation représente la deuxième étape majeure de la rétroprojection. Il s'agit de la pondération des contributions des différents faisceaux à un pixel particulier.
La localisation consiste à récupérer, pour un angle de vue particulier, l'indice de la projection la plus proche et donc du faisceau le plus proche contenant le point considéré (Fig.10.6).
Ensuite comme l'abscisse de la projection de suivant cette direction ne conduit pas au centre du détecteur référencé par , on définit la distance et on procède ensuite à une interpolation (classiquement une interpolation linéaire). La valeur affectée en est alors:La rétroprojection complète s'obtient ensuite par une sommation sur tous les angles de vue possibles:
Cet algorithme s'appuie sur le fait que notre image est discrétisée dans une
base de pixels. Dans cette base, les déplacements pour passer du centre d'un
pixel au centre d'un autre pixel suivant la direction (ou suivant la
direction ) sont quantifiés par la la taille (ou )
du pixel suivant cette direction. Tout déplacement correspond à un multiple
de (ou ). De ce fait, pour un angle de vue
donné, ce déplacement suivant (resp. ) une fois projeté conduit
à un incrément suivant la direction tel que:
(resp.
). Il en découle alors l'algorithme
Alg.7
Même si les auteurs ne précisent pas l'extension au cas 3D, il se déduit assez
facilement de l'algorithme bidimensionnel. Dans l'étape de localisation, il
sera alors nécessaire de repérer les valeurs et les
plus proches du point
correspondant à la projection
d'un voxel ainsi que les distances
et
.
L'interpolation deviendra une interpolation bilinéaire:
Pour implémenter l'algorithme de Shepp et Logan en 3D, il faut également définir,
suivant les directions et , les incréments associés à la projection
des déplacements unitaires sur le volume suivant les directions ,
et .
0 | |||
Le calcul des cosinus et sinus étant une opération coûteuse en temps de calcul, on calcule au préalable des tables de cosinus et sinus que l'on mémorise. On remplace ainsi le calcul par un accès mémoire.
On peut donc maintenant donner une version de l'algorithme de rétroprojection
3D (Alg.8).
L'interpolation intervient au coeur de cinq boucles imbriquées. On comprend
alors l'importance d'avoir diminué le nombre d'opérations requises pour son
calcul. D'autre part, la condition concernant la validité du capteur provient
du fait que certains voxels ne sont pas vus par tous les détecteurs (géométrie
cylindrique d'acquisition et caractère fini du sinogramme).
On constate que la valeur de correspondant à la coordonnée suivant dans le plan de projection est uniquement dépendante des coordonnées du voxel projeté (). L'idée est donc d'inverser l'ordre des boucles par , de sorte qu'avant de rentrer dans la boucle sur les , nous puissions d'ores et déjà calculer . Si cette valeur conduit à un indice n'appartenant pas au sinogramme discret, nous pouvons ainsi éviter une boucle sur les . Ceci est important puisque c'est la boucle la plus externe qui est parallélisée. Après inversion de l'ordre, la boucle sur les correspond à la boucle la plus externe.
Les accès mémoires sont importants et l'inversion de l'ordre de parcours du
volume conduit à un adressage qui n'est plus optimal avec la façon dont l'ordinateur
lit en mémoire. C'est pourquoi dans la rétroprojection, nous construisons un
volume
que nous permutons à la fin de la procédure
pour obtenir . L'algorithme de rétroprojection devient Alg.9.
Nous pouvons également, dans certains cas, faire intervenir des symétries permettant
de minimiser la nombre d'itérations des différentes boucles [20].
D'une manière générale, dans l'ensemble des procédures, nous essayons de limiter
le nombre de divisions à l'intérieur des boucles, car cette opération coûteuse
requiert trois multiplications et une réciproque.
Rappelons que la reconstruction algébrique implique le système linéaire:
Cette opération revient à faire la projection d'une image. En effet, partant d'un vecteur de dimensions (égales au nombre de voxels dans un volume émetteur), on obtient par l'utilisation d'une matrice de dimensions un vecteur de dimensions (égales au nombre de détecteurs dans le sinogramme).
Pour une matrice connue, la réponse impulsionnelle sur les détecteurs pour une source ponctuelle s'obtient en considérant un vecteur (une image) non nul pour une seule de ses composantes (sorte de dirac discret) et par l'application de Eq.10.3. Inversement, dans la colonne de est stockée la réponse impulsionnelle sur les détecteurs pour une source ponctuelle unitaire placée en (Fig.10.7).
|
Cette opération consiste à faire la rétroprojection (backprojection) d'une série de mesure sur les détecteurs. En effet, partant d'un vecteur de dimensions (égales au nombre de détecteurs dans le sinogramme), on obtient par l'utilisation d'une matrice de dimensions un vecteur de dimensions (égales au nombre de voxels dans un volume émetteur)
|
Notons que dans le cas tridimensionnel, la matrice est de dimension conséquente ( ). C'est pourquoi nous envisageons la reconstruction algébrique après rebinning. Le sinogramme sera donc de taille fois plus petite. Mais surtout le rebinning décompose le problème de reconstruction 3D en une série de reconstructions 2D par tranche. Autrement dit, on peut considérer une matrice par tranche transaxiale. Chaque matrice n'est plus alors que de dimension . Si du plus, nous considérons une invariance de la matrice suivant la direction , nous avons alors une seule matrice (notée par la suite) pour les plans à reconstruire. De ce fait, nous travaillons maintenant en 2D. La reconstruction consiste alors partant d'un sinogramme ``rebinné'' à estimer le volume en calculant images 2D. Ceci nous permet de paralléliser avec une forte granularité, car nous obtenons une indépendance entre les reconstructions. Nous pouvons donc maintenant poursuivre le raisonnement en 2D en envisageant la reconstruction d'un plan transaxial précis.
Par la suite, nous faisons les hypothèses supplémentaires:
Une fois l'invariance de la réponse impulsionnelle posée, on peut chercher les redondances présentes dans la matrice produit et donc réduire calculs et espace mémoire en construisant directement cette matrice. On s'affranchit ainsi de la construction de la matrice et du produit matriciel .
Considérons deux composantes et du vecteur image. Exprimer le produit , c'est chercher à exprimer la rétroprojection en de la projection d'une source ponctuelle placée en . La valeur de est la somme des contributions de chaque détecteur pondérées par la réponse impulsionnelle d'une source ponctuelle placée en (Fig.10.8 et Fig. 10.9 ).
Fixons une direction de projection et notons par la distance qui sépare les droites de projection suivant relatives à et . représente toujours l'abscisse de la projection suivant de (Fig. 10.9 ). comme sont dépendants de ( ). Nous avons donc :
on intègre sur :
Puis on intègre le résultat suivant :
apparaît donc comme la fonction d'autocorrélation
de la réponse impulsionnelle à la distance . 9.5
Pour l'intégration de (10.5) suivant , il est nécessaire de définir explicitement en fonction des écarts et séparant les deux pixels considérés.
Avec les notations explicitées Fig.10.10, on obtient simplement:
On pose alors qui nous permet de réécrire la distance comme:
En rappelant que la fonction modifiée de Bessel de premier ordre est telle que:
Un calcul analogue à ce qui vient d'être fait montre que
La fonction peut être prolongée par continuité puisque :
Pour construire la matrice partant de , il n'est pas nécessaire de conserver toutes les valeurs de pour les différents couples possibles dans l'image à reconstruire. En effet, pour construire la valeur de pour un couple particulier, on peut utiliser un tableau à deux dimensions. caractérise les écarts possibles entre pixels suivant la direction pour l'image envisagée et ceux suivant . En effet, les déplacements sur l'image sont quantifiés par la taille du pixel et conduisant à et . On remplit le tableau en utilisant Eq.10.12 et Eq.10.13. Autrement dit . Pour les indices et dans l'image tel que et . On obtient la valeur de par . On peut donc ainsi de manière très simple et très rapide recalculer n'importe quelle ligne où colonne de la matrice produit (Voir Fig.10.12). De plus, on diminue la taille de la matrice de manière importante. En effet, pour une image ECAT HR+ classique, . La taille globale de serait normalement , ici elle n'est plus que de soit fois plus petite. De plus, on s'affranchit du produit matriciel . Pour des raisons de symétrie de révolution de la fonction , on peut encore diminuer la taille du tableau par 4 en travaillant sur les valeurs absolues des écart et .
On obtient alors les algorithmes Alg.10
pour la construction du
et Alg.11
pour le produit matriciel partant de cette matrice.