Maintenant que nous avons vu les principales étapes nécessaires à la reconstruction, nous allons pouvoir les implémenter.

1. L’implémentation
1. Matériel et langage
1. Matériel

Si l’utilisation d’une méthode directe pour la résolution a été envisagée, cela vient certainement du fait que le centre CYCERON possède un supercalculateur qui autorise le calcul sur des volumes de données importants.

Le Centre CYCERON dispose d’un supercalculateur CRAY J916 appartenant à la série des J90. Ce système possède l’architecture vectorielle des systèmes CRAY. Il comporte aujourd’hui 4 modules de 4 processeurs. Le J916 est typiquement une machine à mémoire partagée.

Le CRAY est intégré au réseau VIKMAN et peut ainsi être utilisé en tant que serveur de calcul. Son interface FDDI garantit une grande vitesse d’accès au réseau.

La période d’horloge est de 10 ns (100mhz). Grâce à l’architecture vectorielle, la puissance en crête atteint 200 MegaFlops (nombres d’opérations sur des réels flottants par seconde), soit une puissance en crête de 3.2 GigaFlops pour le système complet. (Une station Silicon Graphics peut atteindre 40 MegaFlops).

La plupart des problèmes nécessitant l’emploi d’un supercalculateur traitent d’un volume de données important et font de nombreux accès à la mémoire. L’efficacité ne dépend pas alors que de la rapidité de calcul proprement dite mais aussi de la bande passante de la mémoire. Pour le Cray J916, tous les processeurs peuvent accéder à n’importe quel mot de la mémoire. Même s’il existe des limitations au accès mémoire, le débit maximal, la bande passante, est de 32 mots mémoires par période d’horloge.

2. Langage

1. Pour les données

A ce niveau, au bout du compte, on souhaite disposer d’un Sinogramme. Le sinogramme est un tableau bidimensionnel dont l’ordonnée donne la direction de projection et l’abscisse un capteur de la barrette placée perpendiculairement à cette direction. A chaque couple (Direction, Capteur) est associée une valeur représentant le nombre de coups pour ce capteur. Dans la pratique, c’est sous cette forme que sont obtenus les résultats de mesure à l’issue d’un examen TEP.

Les Sinogrammes que nous utilisons sont principalement de deux sortes :

0. Sinogrammes Calculés

1. Sinogrammes Mesurés

1. Pour la Reconstruction

1. Pour les Tests
0. Rms Distance[RB8] [RB13]

Dans le cas où les images reconstruites ont été obtenues à partir de sinogrammes calculés, on connaît l’image exacte. Ainsi pour réaliser une évaluation quantitative, ou plus précisément pour mesurer l’écart de l’image reconstruite avec l’image originale, nous calculerons un écart quadratique moyen relatif ( Root Mean Square distance ) :



où :

f_i représente la valeur du pixel i sur l’image reconstruite

t_i représente la valeur du pixel i sur l’image originale

N le nombre de pixels considérés

t_a la valeur moyenne sur cet ensemble de pixels

Cette distance peut être calculée sur la totalité ou seulement sur une partie de l’image suivant les besoins.

1. Moments

1. Résultats et discussion
0. Fantômes utilisés pour les sinogrammes calculés
0. disque

Sur une image de dimension 128*128 pixels, on considère un disque uniforme d’un rayon de 50 pixels centré sur l’image. Ce disque correspond à une zone d’activité constante dont l’intensité variera et sera précisée au cours des expériences.

1. fant1

Le fantôme fant1.128 correspond à une image 128*128 sur laquelle une région elliptique décentrée représente une activité de fond. L’intensité sur cette région est de 5. Deux ouvertures rectangulaires ont été mises à zéro à l’intérieur de cette zone. Les deux régions circulaires haut à gauche et la région elliptique en bas représentent des zones chaudes d’intensités respectives 5, 7.5 et 11.

2. irm

Ce fantôme irm.128 est réalisé à partir d’une image irm de dimension 128*128 segmentée substance blanche/substance grise, sur laquelle on a ajouté 2 points chauds. L’intensité de la substance grise va être de 1, celle de la substance blanche de 4. Les points chauds sont, un disque de rayon 5 pixels d’intensité 14, et un de rayon 2 pixels d’intensité 7.

Les deux derniers fantômes sont représentés dans la table des illustrations Figure 5 et 6.

1. Premières reconstructions et premières constatations
0. Illustration

A titre d’exemple, pour la reconstruction directe, nous allons projeter le disque de différentes manières. La valeur moyenne sur ce disque est alors de 4.82 dans le cas d’une projection par bande et de 2.16 dans le deuxième cas. Pour notre sinogramme, nous projetons suivant 222 directions avec des barrettes de 64 capteurs d’une largeur de 2 pixels. Après reconstruction, les images auront une résolution de 128*128 pixels.

Les résultats sont représentés figure 7 de la table des illustrations.

Le fait de pouvoir choisir une base de reconstruction, un mode de projection, une méthode de projection fournit un grand nombre de combinaisons possibles. Dans l’exemple précédent, on cherchait juste à donner une idée de la reconstruction pour différentes bases de reconstruction. C’est pourquoi seule la reconstruction utilisant une factorisation QR a été envisagée et le mode de remplissage de H était limité au cas le plus simple d’une projection par bande en parcourant l’image pixel par pixel. Nous allons par la suite, chercher à montrer parmi les nombreuses combinaisons quelles sont celles ayant un sens et un intérêt pratique.

Nous avons, dans l’illustration précédente, uniquement utilisé une factorisation QR... Pourquoi ?

1. Factorisation SVD ou QR

Nous avons vu dans le chapitre précédent, qu’à l’aide d’une factorisation QR Householder avec permutation des colonnes, on pouvait trouver une solution au sens des moindres carrés de notre système. Cette solution n’était pas une solution de norme minimale mais s’en approchait. La factorisation SVD permettait également de trouver une solution, cette fois-ci de norme minimale. Dans les deux cas, une certaine régularisation du problème était possible. Alors entre SVD et QR, laquelle choisir ?

Une première remarque, qui ne nous aidera pas dans notre choix, est que les reconstructions obtenues par l’une ou l’autre de ces deux factorisations, sont d’aussi bonne qualité.

On va donc privilégier la méthode présentant le temps de calcul et l’encombrement mémoire les plus faibles.

0. Temps de calcul

A partir d’un sinogramme de 222 directions avec des barrettes de 64 capteurs d’une largeur de 2 pixels, reconstruisons, pour chacune des deux méthodes, une image sur une base de pavés(= base de pixels de décomposition de résolution n*n pixels de représentation). Pour celle-ci, H est calculé en accord avec une projection suivant des bandes.

Dans le cas d’une reconstruction par SVD, le temps global comprenant les saisies et le chargement des images est de 208 secondes, le temps net de calcul est alors de 206 secondes. Dans le cas d’une reconstruction QR, ces temps sont respectivement de 90 secondes et 52 secondes. Ce gain de 4 entre les deux temps de reconstructions est loin d’être un facteur négligeable. Ce gain provient principalement du fait que la factorisation SVD n’est pas parallélisée et donc n’est effectuée que sur un processeur. La factorisation QR, elle, est parallélisée s’effectue donc sur les 16 processeurs du CRAY. Dans notre problème des moindres carrés, la matrice H a plus de lignes que de colonnes. Pour gagner du temps, lors de la factorisation SVD, on calcule donc le produit H^TH puisque cette matrice est plus petite que H. Cependant ce produit, bien que parallélisé, prend du temps.

Une reconstruction sur une base de 64*64 paves prend 608 secondes à l’aide d’une factorisation QR. On est loin des 3 secondes que prend la rétroprojection filtrée.

1. Espace mémoire

La matrice H étant en général de grande taille, cela pose rapidement des problèmes : 222 directions de projections, 64 barrettes de capteurs, une base de 32*32 pixels, cela représente une matrice H de 14208 lignes et de 1024 colonnes, soit pas moins de 14.5 millions d’éléments. On atteint vite les capacités limites de la machine quand on augmente la dimension de la base de reconstruction.

Pour la SVD, Le calcul de H^TH, du moins comme il a été implémenté, suppose la présence simultanée de H^T et H. De plus, lors de la décomposition SVD de H^TH, on doit mémoriser les matrices U, L , et V. Heureusement le calcul de H^TH produit une matrice carrée de taille utile fortement inférieure à celle du nombre d’éléments de la base. En effet, H a des colonnes qui sont nulles. Pour une résolution par SVD l’encombrement mémoire est donc d’un peu plus de 2 fois H.

Dans le cas de la reconstruction QR, on s’est affranchi du produit H^TH et lors de la factorisation, la matrice Q n’est pas calculée explicitement. Seuls les vecteurs v sont conservés (voir chapitre précédent) et du fait des zéros qui apparaissent sous la diagonale lors de la factorisation, ces vecteurs peuvent être directement mémorisés dans la matrice H. La reconstruction par une factorisation QR nécessite donc un peu plus de H en place mémoire. Ce gain de 2 n’est pas non plus négligeable.

Pour les deux raisons précédentes, dans le cas d’une solution au sens des moindres carrés , on préférera donc la factorisation QR Householder avec permutation des colonnes..

Dans les reconstructions qui suivront, on n’envisagera, sauf mention particulière, que ce type de résolution. De plus pour limiter le nombre de possibilités, on n’envisage le remplissage de H qu’en accord avec un protocole de projection par bande et nous nous limitons à une base de pavés.

2. Reconstruction et géométrie d’acquisition
0. Influence du nombre de capteurs

Pour l'image fant1.128, on va calculer les sinogrammes en conservant un nombre de projections constant (212 directions), mais en faisant varier le nombre de capteurs sur une barrette. On s'arrange pour que la totalité de l'image soit vue en permanence. A l ’aide de ces sinogrammes et pour un pavage constant, on va reconstruire successivement les images. Le fantôme original est injecté comme information à priori et le calcul de la Rms distance est effectué sur toute l'image. On obtient alors les résultats suivants :

NbCapteurs	16.00	32.00	64.00	128.00
Nb Pixels	8.00	4.00	2.00	1.00
Rms dist. 16 Pavés	4.07	1.23	0.93	0.89
Rms dist. 32 Pavés	6277	9.14	2.89	2.68
Rms dist. 48 Pavés	66.51	23.24	6.28	5.07

Si maintenant, on reconstruit sans fournir d'information a priori et que pour le calcul de la Rms distance on ajuste la résolution de l'image source au pavage de la reconstruction, on obtient les résultats suivants :

NbCapteurs	16.00	32.00	64.00	128.00
Nb Pixels	8.00	4.00	2.00	1.00
Rms dist. 16 Pavés	52.11	17.09	14.58	12.99
Rms dist. 32 Pavés	337.21	37.44	12.69	9.65
Rms dist. 48 Pavés	165.70	61.52	15.58	10.17

Même si les distances mentionnées dans ce tableau ne peuvent être appréhendées dans l'absolu, on peut en dégager une tendance générale. On voit que pour une base de dimension fixée (NbPavés), la distance après reconstruction diminue à mesure que le nombre de capteurs augmente jusqu'à ce que ce nombre soit supérieur à la résolution. La distance devient alors constante.

Cette constatation est importante et il est nécessaire de s’y arrêter plus longuement. Pour cela, on va considérer le fantômes irm.128. On va réaliser une série de sinogrammes utilisant le modèle de projection par bande, en faisant varier le nombre de directions, et le nombre de capteurs. Dans chaque cas, les sinogrammes seront bruités suivant un modèle de bruit poissonien.

Le tableau suivant regroupe les caractéristiques du sinogramme, la rms distance calculée vis à vis de l’image originale et les valeurs extrêmes obtenues après reconstruction sur une base de 32*32 (ou de 64*64 pavés). La rms distance est calculée seulement sur une masque elliptique englobant la zone active et en référence à l’image originale ajustée à la résolution de 32*32 (ou de 64*64 pavés)

 

 

 

Sinogrammes	NbCapteurs	32(4pixels)			64(2pixels)			128
	Nb Directions	100	500	1000	40	250		150
	NbPavés	32(4pixels)			32	32	64	64
Reconstruction	Rms distance	353.7	188.6	130.9	50.3	24.6	140.3	29.9
	Minimum	-23.6	-18.15	-7.58	-3.31	-0.67	-9.37	-1.58
	Maximum	28.13	24.14	13.47	12.24	11.56	20.01	16.38
Figure		Figure 9					Figure 10
		1	2	3	4	5	1	2

 

Dans le cas de la reconstruction de l’image à partir du sinogrammes à 128 capteurs, la rms distance vis à vis de l’image originale (si cette dernière n’est pas ajustée au pavage de 64*64) est de 76.5. A titre de comparaison, la rms distance obtenue par une rétroprojection filtrée est de 64.37.

Le tableau ci dessus est assez explicite. Lorsque le nombre de pixels par pavés de reconstruction est inférieur, à la limite égal, au nombre de pixels que voit un capteur (image 1,2,3), le système malgré un nombre de données toujours croissant (jusqu’à 1000 directions) aura du mal à retrouver l’image initiale.

Ce problème n’est donc pas lié à un nombre de redondances insuffisant et se manifeste, quelle que soit la base de reconstruction, dès que l’on cherche à reconstruire une image avec une résolution (au sens de la représentation) inférieure à celle que voit le capteur.

C’est une limitation forte de notre système, qui vient certainement d’une violation du théorème d’échantillonnage dans ce cas là. Il s’avère alors nécessaire de penser différemment la géométrie d’acquisition. [RB11]

Par la suite, on suppose que les différentes reconstructions sont en accord avec cette règle :

Résolution de l’image reconstruite < Résolution des Capteurs

Pour une base de pavés, cette condition équivaut à :

Nb de Pavés < Nb de Capteurs

1. Nombres moyens de données par éléments de la base

Cette règle étant posée, on peut maintenant chercher à savoir comment la reconstruction de sinogrammes bruités va évoluer en fonction du nombre de données disponibles.

Pour un fantôme, on va réaliser une série de sinogrammes, pour lesquels varie le nombre de directions. Le nombre de capteurs d’une largeur de 2 pixels est fixé à 64 . Chacun de ces sinogrammes, calculés par une projection par bande, va ensuite être bruité par un bruit additif gaussien dont la variance sera successivement 2,5,9,14,20% du nombre de coups maximum sur le sinogramme. On va reconstruire les images sur une base de pavés dont le nombre d’éléments sera successivement 16², 24², 32², 40², 48², 56² . Lors de la reconstruction, on injecte comme information à priori le fantôme lui même. On calcule ensuite la rms distance, la moyenne et la variance sur la zone uniforme de ce fantôme par rapport à l’image originale. On va s’intéresser plus particulièrement aux critères de rms distance et de variance en fonction du nombre moyen de données. Par nombre moyen de données on comprend le rapport du nombres de données de projections (Nbdirection*NbCapteurs) sur le nombre d’éléments de la base (NbPavés²).

On obtient pour le fantôme fant1 :





et pour le fantôme disque uniforme où l’intensité est de 4.82.



Ce critère nombre moyen de données permet donc de traduire les redondances du système. Si ce nombre est assez élevé, la reconstruction même bruitée sera assez proche de l’image originale.

Pour une distance à l’image originale fixée, on peut ainsi connaître, à partir d’une telle courbe, le nombre de données, ie le nombre de projections et de capteurs requis.

3. Reconstruction et information a priori
0. Ajuster une surface et non un plan.

Dans le cas d’une résolution sur une base de pavés, on a déjà vu précédemment que le pavage ne pouvait être choisi de manière quelconque et qu'il était étroitement lié au nombre de capteurs dont on dispose lors de la mesure. (paragraphe 1.2.3.1).

Dans le cas où la base de pavés n’est pas contrainte par une matrice d’information a priori, la résolution de l’image reconstruite est celle du nombre de pavés. Il semble normal alors qu’en augmentant le pavage, l'image reconstruite va se rapprocher de l'image originale qui est de résolution 128*128.

Mais une autre constatation s'impose...

Pour cela, à l'aide du sinogramme issu du fantôme fant1.128 projeté suivant 200 directions sur des barrettes de 60 capteurs de 4 pixels, on va reconstruire des images en augmentant régulièrement le nombre de pavés. La matrice H est conditionnée par l'information a priori exacte fant1.128. On va s'intéresser alors à la distance quadratique moyenne ainsi qu'aux propriétés statistiques sur une partie uniforme de l'image fant1.128.masque.

les résultats obtenus sont les suivants :

NB Pavés	16	20	24	28	32	48
Rms Dist.	0.53	0.95	1.43	1.71	2.42	4.88
Moyenne	4.95	4.95	4.95	4.95	4.95	4.95
Variance	0.172	0.175	0.175	0.180	0.184	0.232

 

On s'aperçoit que lorsque le nombre de pavés augmente, la Rms distance mais également la variance sur une zone homogène de l'image augmente. En extrapolant on serait tenté de déduire que la reconstruction est meilleure dans le cas où le pavage est unitaire!!! C'est effectivement le cas...

Lors de la reconstruction conditionnée par de l’information a priori, on cherche, sur chaque pavé, le coefficient par lequel il est nécessaire de multiplier cette portion de surface injectée en information pour la faire coller le mieux qui soit en moyenne à la portion de surface correspondante de l'image originale

De fait, si l'information a priori est exacte (si la surface injectée en information a priori est identique en forme et en intensité à l'image à reconstruire), le fractionnement de cette image pour ajuster chaque partie induit, de par ce découpage, plus d'erreurs que si l'on cherchait à ajuster l'image a priori dans son ensemble. Ajuster l’image dans son ensemble consiste simplement à affecter le coefficient 1 à cette image.

Le découpage de l'image en pavés se justifie dans la pratique par le fait qu’on ne connaît évidemment pas exactement l'image à reconstruire. Autrement dit, on ne peut pas limiter la base de décomposition (les portions de surface) à un seul vecteur (une seule image ou un seul pavé), et on est obligé d'augmenter la dimension de l'espace de reconstruction (Nb de Pavés) pour que la décomposition sur cette base fournisse une bonne estimation de l'image originale.

1. Taille de l'information a priori

1. Avec ou sans Information a Priori ?

1. Information a priori Exacte ou Approchée