Analyse et comparaison de procédures d’optimisation en tomodensitométrie

Maxime Toussaint, Christian Thibaudeau, Roger Lecomte et Jean-Pierre Dussault

Université de Sherbrooke Courriel: maxime.toussaint@usherbrooke.ca

8 décembre 2011

Résumé

La tomodensitométrie est une modalité d’imagerie médicale utilisée pour représenter la structure interne d’un être vivant. La reconstruction de cette image doit être rapide et stable. En outre, le résultat doit représenter le plus fidèlement que possible la réalité. Ce document propose de comparer les performances de quelques algorithmes itératifs en termes de qualité d’image reconstruite et de durée de calcul. Parmi les algorithmes testés, certains ont été développés spécialement pour le problème de tomodensitométrie alors que d’autres sont des solveurs généraux. De plus, chacun des algorithmes sera testé avec une variante incluant un terme de régularisation. Ce document propose de tester quelques variantes de gradient conjugué non-linéaire afin d’observer leur performance en tomodensitométrie.

Table des matières

1 Introduction 3

2 Algorithmes utilisés 4 2.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Concept commun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 APT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 SG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Description des images tests utilisés 13 3.1 Shepp-logan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Shepp-logan facile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 H.G.D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4 Catphan CTP401 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

4 Analyses accomplies 14 4.1 Observer l’évolution des solutions au fil des itérations . . . . . . . . . . . . . . . . . . . . . . 15 4.2 Comparer les solutions obtenues avec un critère d’arrêt . . . . . . . . . . . . . . . . . . . . . . 16

5 Résultats obtenus 17 5.1 Comportement des algorithmes sur 80 itérations . . . . . . . . . . . . . . . . . . . . . . . . . 17 5.2 Étalon des solutions obtenues à l’arrêt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

6 Idée à explorer 19 6.1 Les paramètres de la fonction de pénalisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6.2 La majoration de la courbure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6.3 Optimisation séquentielle de sous espaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6.4 Approfondir l’expérience . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

7 Remerciements 20

8 Annexe 20 8.1 Majoration utilisée dans pasAdmissible2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 8.2 Zone d’importance pour l’évaluation du LMH et de la variance . . . . . . . . . . . . . . . . . 21 8.3 Configuration d’acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 8.4 Mires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 8.5 Code des mires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 8.6 Comparaison d’images obtenues avec quelques algorithmes . . . . . . . . . . . . . . . . . . . . 25 8.7 Figures utilisées dans l’analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

1 Introduction

La tomodensitométrie est une procédure qui consiste à exploiter le taux d’absorption de rayons-x des tissus pour représenter la structure interne d’un être vivant. Ce procédé se divise en trois étapes, soit l’acquisition, le traitement des données et la reconstruction. Ce document se concentra sur cette dernière étape dans le cas particulier d’une géométrie d’acquisition en seulement 2 dimensions.

Dans la littérature, deux familles de méthode sont souvent reliées à la tomodensitométrie. Les méthodes analytiques qui reposent sur l’inverse de la transformée Radon et les méthodes itératives qui optimisent des modèles correspondant au problème. Ces deux procédés sont possibles grâce à un théorème démontré par M. Johann Radon en 1917 [RP86].

Théorème de projection de Radon : Si la totalité des projections d’un objet est disponible et que ces dernières sont acquises selon des droites parallèles, il est alors possible de reconstruire la fonction réelle à deux variables (position spatiale en x et en y) représentant cet objet.

En pratique, il est impossible d’obtenir toutes les projections. Toutefois, un sous échantillon judicieusement choisi permet d’acquérir suffisamment de données pour conserver l’idée du théorème. De plus, les limitations des ressources informatiques ne permettent pas de représenter l’objet dans sa forme continue. La discrétisation en voxels permet de pallier à cette lacune. Pour plus d’informations concernant l’acquisition et le traitement des données, le lecteur est invité à consulter [KS88] et [H80].

Deux caractéristiques importantes en tomodensitométrie différencient les algorithmes analytiques de ceux dits itératifs : la durée de reconstruction et la flexibilité de la méthode. Pour se représenter la différence, il faut comprendre leur fonctionnent. La reconstruction par calcul analytique consiste à appliquer la fonction inverse sur les résultats obtenus. L’idée de base veut que la méthode d’acquisition des données peut être représentée par une fonction qui, en utilisant le théorème de projection de Radon, admet un inverse. Le principe clé des algorithmes itératifs est de représenter la problématique sous forme d’un problème d’optimisation et de chercher, à chaque étape, une solution qui se rapproche de la vérité. Il s’agit ensuite d’arrêter lorsque cette solution est jugée satisfaisante.

Pour comparer la vitesse de reconstruction des familles d’algorithmes, il faut introduire l’élément le plus coûteux à manipuler, soit la matrice système. Cette matrice permet de relier les données acquises à la version discrète de l’objet. Elle devient rapidement imposante puisque sa forme explicite est composée de N_{p} × N_{d} × N_{v} éléments (où N_{p} représente le nombre deprojections, N_{d} le nombre de détecteurs et N_{v} le nombre de voxels). Ainsi, la vitesse de reconstruction des méthodes dépendra largement de la taille de cette matrice et du nombre de fois qu’elle doit être utilisée. La méthode analytique demande l’application d’une formule, soit l’inverse de la fonction de Radon, ce qui consiste à n’utiliser qu’une seule fois la matrice système pour approximer notre objet. La méthode itérative demande l’application de cette matrice au minimum une fois par itération. Donc, la méthode itérative demande beaucoup plus de calcul que la méthode analytique.

Toutefois, les algorithmes itératifs sont connus pour être plus flexibles que les algorithmes analytiques. D’une part, il est possible d’inclure dans ces premiers la mathématique liée à la physique du processus d’acquisition (telle que la nature Poissonnienne du comptage de photons). D’un autre côté, alors que les algorithmes analytiques doivent travailler avec des données spatialement équidistribuées (et souvent uniquement en géométrie parallèle), les méthodes itératives ne souffrent d’aucunes limitations de ce type. En modélisant correctement la géométrie du système d’acquisition, par le biais de la matrice système, toutes les configurations possibles de scanners peuvent être utilisées pour produire des images sans ré-interpolation des données. Les algorithmes itératifs sont principalement composés de deux éléments. D’abord, il faut un modèle qui se comporte bien, c’est-à-dire continu et dérivable, et dont l’image recherchée correspond à la solution optimale. Puis, il faut une méthode d’optimisation qui converge globalement vers la solution désirée. Les détails seront expliqués dans la section suivante.

Le but de ce document est de comparer les performances de différentes méthodes itératives. Afin d’offrir une vue d’ensemble sur les méthodes disponibles en tomodensitométrie, des algorithmes reconnus dans le domaine ont été employées ainsi que des variantes de deux solveurs généraux qui semblent prometteuses. Ainsi, la liste est composée de méthodes développées uniquement pour résoudre des problèmes de tomographie et de solveur général. Il est à noter que dans ce travail, ces deux groupes sont respectivement nommés APT^{1} et SG^{2}. De plus, chacun de ces algorithmes sera testé avec et sans un terme de régularisation, ce qui permettra d’observer leurs performances comparatives. Pour finir, les expériences incluent un APT qui se sert de sous-groupes ordonnés, une méthode reconnue dans la littérature pour sa rapidité de convergence.

Avant de présenter les résultats, ce document propose d’introduire les éléments testés. Pour commencer, la motivation et la forme explicite des modèles et des algorithmes inclus dans les expériences seront présentées. Ensuite, les mires utilisées dans les simulations seront décrites. Suivra la description de la méthodologie appliquée pour effectuer les tests et les comparaisons. Les résultats, appuyés par des graphiques, viendront illustrer les performances des algorithmes. Finalement, une évaluation des performances et une ouverture sur ce qui pourrait être fait pour approfondir le sujet se trouvent en conclusion de ce document. Notez que la plupart des figures se trouvent en annexe. Celle présentées dans le corps du texte ont simplement pour but de faciliter la compréhension lors de la lecture.

2 Algorithmes utilisés

Cette section décrit les algorithmes manipulés pendant ce projet. Tout d’abord, la notation utilisée pour décrire les algorithmes sera introduite. Ensuite, les concepts communs aux APT et aux SG seront présentés. Puis, chacun des algorithmes seront expliqués. Notez que les trois premiers ont été introduits dans [LF95] et que les preuves de leur convergence sont expliquées en détail dans celui-ci.

1. APT : Algorithmes Produits pour la Tomodensitométrie

2. SG : Solveurs généraux

2.1 Notations

Posons : I  Nombre de projections N_{p} multiplié par le nombre de détecteur N_{d}
J  Nombre de voxels
μ  Vecteur colonne de J éléments tel que μ_{j} est la moyenne d’absorption dans le voxel j
μ^{∗}  Valeur courante de l’estimé de μ à une itération donnée
l  Matrice I × J tel que l_{ij} est la longueur d’intersection entre la projection i et le pixel j
i  Indice utilisé pour représenter les différentes projections
j  Indice utilisé pour représenter les différents voxels
k  Indice utilisé pour représenter les voxels qui partagent au moins un vertex avec la

projection i

b  Indice utilisé pour représenter les différents détecteurs
Y  Vecteur colonne de I éléments tel que Y_{i} est le nombre de photons détecté à la

projection i lorsque l’objet est présent

Ȳ  Vecteur colonne de I éléments tel que Ȳ_{i} = Y_{i}/d_{b} où b = i mod N_{d}
Ỹ  Vecteur colonne de I éléments tel que Ỹ_{i} = ln(Ȳ_{i})
d  Vecteur colonne de N_{d} éléments tel que d_{b} est le nombre de photons

qui est compté par le détecteur b lorsque l’objet est absent

11  Vecteur ligne de I éléments tel que I_{i} = 1 ∀i
γ  Constante utilisée dans la fonction pénalisation
w_{j_{1}j_{2}}  Constantes utilisées dans la fonction pénalisation
κ  Ensemble de J valeurs tel que κ_{j} représente le poids relatif du voxel j
ϵ  Constante utilisée dans la fonction potentiel
x_critique Valeur utilisée pour évaluer ϵ
η  Valeur utilisée pour évaluer ϵ

2.2 Concept commun

Modèle utilisant la log-vraisemblance

Les APT présentés dans ce document reposent tous sur le même modèle. Seulement la méthode permettant de trouver la solution optimale varie. Il faut donc commencer par présenter ce modèle.

Le modèle en question est la vraisemblance^{3} logarithmique des données obtenues. Cette méthode possède des propriétés théoriques intéressantes et est reconnue pour son efficacité dans divers domaines. Elle permet en autre de modéliser correctement le fait que le nombre de photons compté suit une loi de Poisson. Cette approche permet d’obtenir une fonction objectif qui répond à la problématique. Il est évident, par la définition de la vraisemblance, que la structure recherchée est une solution optimale de ce modèle. Il est aussi possible de démontrer, avec l’hypothèse que les données ne sont pas trop bruitées, que le maximum du logarithme de la vraisemblance produit une fonction objectif qui a pour unique solution optimale celle recherchée. Pour finir, les propriétés du logarithme naturel vont permettre de simplifier la fonction objectif. Notez que le modèle inclut a posteriori le comportement Poissonnien du taux d’absorption. Ainsi, l’algorithme devrait être en mesure de produire des images satisfaisantes même lorsque le nombre de rayon-x utilisé est statistiquement faible.

Origine de la fonction objectif : Posons g(Y, μ) la vraisemblance de μ sachant Y et f_{i}(Y_{i}, μ) la vraisemblance de μ sachant Y_{i}.

3. La vraisemblance est une fonction qui correspond à se demander : "Quel devrait être la valeur des paramètres considérant les résultats obtenus et le modèle statistique qu’ils suivent ?"

Notez que la probabilité qu’un photon réussit à traverser un objet obéit à une loi de Poisson.

ln(g(Y, μ)) = ln(∏_{i=0}^{I} f(Y_{i}, μ))  Y_{i} sont indépendants
= ∑_{i=0}^{I} ln(f(Y_{i}, μ))  Prop. de ln
= ∑_{i=0}^{I} ln(\frac{e^{−d_{i}e^{−⟨l_{i},μ⟩}} (d_{i}e^{−⟨l_{i},μ⟩})^{Y_{i}}}{Y_{i}!} )  Déf. loi de poisson
= ∑_{i=0}^{I} ln(e^{−d_{i}e^{−⟨l_{i},μ⟩}}) + ln((d_{i}e^{−⟨l_{i},μ⟩})^{Y_{i}}) − ln(Y_{i}!)  Prop. de ln
= ∑_{i=0}^{I} −d_{i}e^{−⟨l_{i},μ⟩} + Y_{i}(ln(d_{i}) + ln(e^{−⟨l_{i},μ⟩})) − ln(Y_{i}!) Prop. de ln
= ∑_{i=0}^{I} −d_{i}e^{−⟨l_{i},μ⟩} + Y_{i}ln(d_{i}) − Y_{i}⟨l_{i}, μ⟩ − ln(Y_{i}!)  Prop. de ln
= ∑_{i=0}^{I} −d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}⟨l_{i}, μ⟩ + constante  μ est la seule variable

Ainsi, l’image recherché est la solution du problème d’optimisation suivant :

arg max_{μ} ln(g(Y, μ)) = arg max_{μ} ∑_{i=0}^{I} −d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}⟨l_{i}, μ⟩

Notez : L(μ) = ∑_{i=1}{−d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}⟨l_{i}, μ⟩}
\frac{∂L(μ)}{∂μ_{j}} = ∑_{i=1} l_{ij}[d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}]
\frac{∂^{2}L(μ)}{∂μ_{j}∂μ_{k}} = − ∑_{i} l_{ij}l_{ik}d_{i}e^{−⟨l_{i},μ⟩}

Fonction de pénalisation

Dans un monde parfait, les données obtenues lors de l’acquisition ne seraient pas affectées par des limitations physiques telles que la diffusion des rayons-X et le bruit de mesure des détecteurs. De plus, il est admis que les structures reconstruites en tomodensitométrie sont généralement composées de zones assez uniformes. Il est donc parfois avantageux de forcer les algorithmes à obtenir des solutions qui respectent ces propriétés. Une démarche reconnue en optimisation est d’ajouter une fonction de pénalisation à la fonction objectif de départ. Par exemple, le but des APT est de maximiser la log vraisemblance. Donc, il faut une fonction qui nuit à la maximisation lorsque la solution est peu uniforme. Son effet doit augmenter lorsque des voxels voisins possèdent des valeurs différentes tout en évitant de favoriser le rapprochement de voxels voisins très différents. Pour ce projet, le choix est tombé sur une fonction suivant l’idée de la variation totale. Voici la forme générale de la fonction de pénalisation implantée :

Ψ(μ) = γ (∑_{j}^{J} ∑_{k}^{K_{j}} κ_{j}κ_{k}w_{jk}ϕ(μ_{j} − μ_{k}))

tel que

ϕ(μ) = \sqrt{ϵ + x^{2}}

est la fonction potentiel.

Dans les SG de ce projet, il faut additionner cette fonction à la fonction objectif pour satisfaire le problème de minimisation. De son côté, la constante gamma établit la force de la régularisation et le compromis entre la résolution spatiale et le bruit dans l’image reconstruite. L’epsilon, produit une barrière qui dicte la valeur de x_critique, où la fonction potentiel passe d’une croissance quadratique à une croissance linéaire. Ainsi, lorsque la différence entre deux voxels dépasse cette valeur, la fonction considère que les voxels ne sont pas de la même valeur et aura moins tendance à pénaliser leur différence. L’idée est de s’assurer que les délimitations entre deux zones uniformes ne soient pas trop lissées. Voici la formule pour calculer la valeur d’epsilon pour les tests présentés dans ce document^{4} :

ϵ = (\frac{x_critique}{η})^{2} (1 − η)^{2}

4. Cette formule correspond à trouver ϵ tel que ϕ′(x_{critique}) = η

Où η ∈ (0, 1) tel que plus η est proche de 1, plus la fonction potentiel va tenter de forcer le rapprochement de voxels voisin. Les omégas servent à pondérer la valeur de la fonction potentiel selon la proximité des voxels j et k. Dans ce travail, les poids associés aux voxels partageant deux vertex, un vertex et zéro vertex sont respectivement 1, 1/\sqrt{2} et 0. Enfin, l’essence des kappas est de pondérer chacune des sommes afin d’obtenir une résolution uniforme dans l’image. La preuve qu’une pénalisation standard produit des images de résolution non uniforme, la motivation des kappas et plus encore sont disponible dans [FR96]. De plus, les auteurs de ce rapport technique présentent une méthode pour produire une table de gamma, propre à chaque configuration d’acquisition, qui permet de choisir la résolution désirée.

Voici la formule pour calculer les kappas pour les tests présentés dans ce document :

κ_{j}(μ) = \sqrt{∑_{j} l_{ij}^{2}Y_{i}/ ∑_{j} l_{ij}^{2}}

2.3 APT

L’aticle suivant [LF95] explicite les motivations et les preuves de convergences des algorithmes APT présentés dans le présent document.

Algorithme EM

Pour commencer, constatez que le problème est plus facile à résoudre si l’ensemble des données est disponible. C’est-à-dire que non seulement le nombre de photons sortant de la source, d_{i}, et sortant de l’objet, Y_{i}, sont connus, mais aussi celui qui entre et sort de chaque voxel, pour chaque projection. Notons U_{ij} et V_{ij} étant respectivement le nombre de photon rentrant et sortant au voxel j sur la projection i. L’algorithme EM consiste à répéter deux étapes. Pour commencer, il calcule l’espérance de la fonction objectif en utilisant le dernier estimé de μ. Puis, il trouve le μ qui maximise la fonction obtenue à la première étape. Pour formuler la fonction obtenue dans la première étape, notez Q(μ|μ^{n}) où μ^{n} est le dernier estimé. Il devient nécessaire de connaître M_{ij} = E(U_{ij}|Y_{i}, μ^{n}) et N_{ij} = E(V_{ij}|Y_{i}, μ^{n}). Or, [LC84] démontre que M_{ij} = Y_{i} + d_{i}e^{− ∑_{k∈S_{ij}} l_{ik}μ_{k}^{∗}} − d_{i}e^{−⟨l_{i},μ^{∗}⟩}∀i, j et que N_{ij} = Y_{i} + d_{i}e^{− ∑_{k∈S_{ij} ∪j} l_{ik}μ_{k}^{∗}} − d_{i}e^{−⟨l_{i},μ^{∗}⟩} ∀i, j. Notez que S_{ij} est l’ensemble des voxels se trouvant entre la source et le pixel j, sur la projection i. Ainsi, Q(μ|μ^{∗}) = ∑_{ij}[N_{ij}l_{ij}μ_{j} + (M_{ij} − N_{ij})ln(1 − e^{−l_{ij}μ_{j}})].

Toutefois, il est difficile d’expliciter le μ qui maximise cette fonction. Malgré tout, il a été démontré que la fonction Q(μ|μ^{∗}) admet seulement un extremum et qu’il s’agît d’un maximum global. Il existe plusieurs méthodes pour approximer la solution cherchée. Celle choisie pour l’algorithme EM est d’approximer l’élément problématique, soit (e^{−l_{ik}μ_{k}} − 1)^{−1} de la fonction par une approximation de sa série de Taylor. Après avoir pris en considération la précision voulue et le nombre de calculs, les auteurs de [LC84] ont décidé d’utiliser une approximation d’ordre 3. C’est-à-dire que \frac{1}{s} − \frac{1}{2} + \frac{s}{12} +O(s^{3}) est utilisé pour approximer \frac{1}{(e^{s}−1)} puisque \frac{1}{(e^{s}−1)} = \frac{1}{s} − \frac{1}{2} + \frac{s}{12} + (s^{3}) où s = l_{ik}μ_{k}. Ainsi, les valeurs approximatives qui maximisent Q(μ|μ^{∗})

sont égales à μ_{j} = \frac{[11(M_{.j}−N_{.j})]}{\frac{1}{2} [M_{.j}+N_{.j}]^{t}l_{.j}} ∀j.

Espérance-Maximisation (μ, Ȳ, l, S_{ij})

{ Donnés : μ_{j} > 0 ∀j; }
{ Indice des pixels entre la source i et le pixel j : S_{ij} }
{ Q(μ|μ^{∗}) = ∑_{ij}[N_{ij}l_{ij}μ_{j} + (M_{ij} − N_{ij})ln(1 − e^{−l_{ij}μ_{j}})] }

répéter

M_{ij} ← Ȳ_{i} + e^{− ∑_{k∈S_{ij}} l_{ik}μ_{k}^{∗}} − e^{−⟨l_{i},μ^{∗}⟩} ∀i, j ; N_{ij} ← Ȳ_{i} + e^{− ∑_{k∈S_{ij} ∪j} l_{ik}μ_{k}^{∗}} − e^{−⟨l_{i},μ^{∗}⟩} ∀i, j ; μ_{j}^{∗} ← \frac{[11(M_{.j}−N_{.j})]}{\frac{1}{2} [M_{.j}+N_{.j}]^{t}l_{.j}} ≈ max_{μ>0} Q(μ|μ^{∗}) ∀j

jusqu’à ( Critère de convergence satisfait)

Résultats ← μ^{∗}

Algorithme 1 – trml

L’algorithme 2 est une modification de la méthode de Newton appliquée à la fonction objectif pénalisée. Notez que la fonction de pénalité exploitée est une majoration de Ψ qui respecte l’idée des algorithmes EM et que seulement les éléments diagonales du Hessien sont utilisés.

Majoration de Ψ :  V (μ|μ^{∗}) = \frac{γ}{2} (∑_{j} ∑_{k}(w_{jk}ϕ(μ_{j} − μ_{k}^{∗} − μ_{j}^{∗})))
+\frac{γ}{2} (∑_{j} ∑_{k}(w_{jk}ϕ(μ_{k} − μ_{k}^{∗} − μ_{j}^{∗})))
Fonction optimisée : −Υ(μ|μ^{∗}) = −Q(μ|μ^{∗}) + V (μ|μ^{∗})

Espérance-Maximisation pénalisé (μ, Ȳ, l, S_{ij})

{ Donnés : μ_{j} > 0 ∀j; }
{ Q(μ|μ^{∗}) = ∑_{ij}[N_{ij}l_{ij}μ_{j} + (M_{ij} − N_{ij})ln(1 − e^{−l_{ij}μ_{j}})] }

répéter

M_{ij} ← Ȳ_{i} + e^{− ∑_{k∈S_{ij}} l_{ik}μ_{k}^{∗}} − e^{−⟨l_{i},μ^{∗}⟩} ∀i, j ; N_{ij} ← Ȳ_{i} + e^{− ∑_{k∈S_{ij} ∪j} l_{ik}μ_{k}^{∗}} − e^{−⟨l_{i},μ^{∗}⟩} ∀i, j ; μ_{j}^{∗} ← μ_{j}^{∗} − \frac{(∇Υ(μ|μ^{∗})|_{μ=μ^{∗}})_{j}}{(∇^{2}Υ(μ|μ^{∗})|_{μ=μ^{∗}})_{j}} ∀j

jusqu’à ( Critère de convergence satisfait)

Résultats ← μ^{∗}

Algorithme 2 – trml_P

Algorithme du gradient

L’architecture de cette procédure provient d’une observation faite sur l’algorithme EM dans le cas de l’émission. En effet, elle peut être représenté sous la forme suivante : λ^{n+1} = λ^{n} + D(λ^{n})dL(λ^{n})′. Où λ est l’équivalent du coefficient d’atténuation dans le cas de l’émission et D(λ^{n}) est une matrice diagonale telle que D(λ^{n})_{jj} = \frac{λ_{j}^{n}}{∑_{i∈J_{j}} C_{ij}}. Pour la transmission, cette reformulation est complexe. Toutefois, la forme générale

de la formulation dans le cas de l’émission permet de développer un algorithme pour la transmission. Pour commencer, il faut changer λ et la dérivée de sa fonction objectif par leur équivalent dans le cas de la transmission. Trouver l’équivalent de D(λ^{n}) pour la situation étudiée exige une observation de son comportement dans sa formulation de départ. Afin que le comportement de D(λ^{n}) soit semblable dans les deux cas, il faut que D(μ^{n})_{jj} = \frac{μ_{j}^{n}}{∑_{i∈J_{j}} Y_{i}l_{ij}}. Donc, la formule obtenue est de la forme suivante : μ_{j} = μ_{j}^{n} + D(μ_{j}^{n})_{j}∇L(μ_{j})_{j}

∀j. Cet algorithme demande moins de calcul que la version EM puisqu’il nécessite le calcul de moins d’exponentielles. Cependant, il n’assure pas que L(μ) augmente pour chaque itération. Cette propriété est obtenue en calculant un pas admissible à chaque itération.

Descente par le gradient (μ, Ȳ, l, UseArmijo)

{ Donnés : μ_{j} > 0 ∀j; }

répéter

dir_{j} ← \frac{μ_{j}^{∗}}{Y^{t}l_{i.}} ∇L(μ^{∗})_{j} ∀j ; f ← ∑_{i=1}\{−d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}⟨l_{i}, μ⟩\}|_{μ=μ^{∗}} ; g ← ∑_{i=1} l_{ij}[d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}]|_{μ=μ^{∗}} ; θ ← trouveP asAdmissible(f, g, dir) ; μ_{j}^{∗} ← μ_{j}^{∗} + θ × dir_{j} ∀j

jusqu’à ( Critère de convergence satisfait) Résultats ← μ^{∗}

Algorithme 3 – Gradient

Où trouvePasAdmissible renvoie θ tel que θ respecte le critère d’Armijo^{5}. La méthode commence par tester θ = 1 et divise par 2 le θ jusqu’à ce qu’il trouve un θ satisfaisant.

trouvePasAdmissible (f, g, dir)

{ Donnés : f valeur de la fonction objectif, g sa dérivée et dir la direction de la mise à jour de l’algo ; }
{  F est la fonction objectif }

τ_{0} ← 0.25 prodGradDir ← ⟨g, dir⟩ tantque ([F (μ + θ ∗ dir) − f] < τ_{0} ∗ θ ∗ prodGradDir )

θ ← θ/2

Résultats ← θ

Algorithme 4 – trouvePasAdmissible

L’algorithme qui suit est une modification de la méthode de Newton appliquée sur la fonction objectif pénalisée. L’équation maximiser est la somme de −Ψ(μ) et de Q(μ|μ^{∗}) tel que

Q(μ|μ^{∗}) = ∑_{j} (μ_{j}(l_{ij}e^{−⟨l_{i},μ⟩}) − \frac{μ_{j}^{2}}{2} \frac{∑_{i} Ȳ_{i}l_{ij}}{μ_{j}^{∗}} )

5. Voir condition d’arrêt de trouvePasAdmissible

Descente par le gradient pénalisé (μ, Ȳ, l, UseArmijo)

{ Donnés : μ_{j} > 0 ∀j; }
{ Δ(μ) = Q(μ|μ^{∗}) − Ψ(μ) }

répéter

dir_{j} ← \frac{μ_{j}^{∗}}{∑_{i} Y_{i}l_{ij}+μ_{j}^{∗}(∇^{2}Ψ(μ))_{jj}|_{μ=μ^{∗}}} (∇Δ(μ))_{j}|_{μ=μ^{∗}} ∀j ; f ← ∑_{i=1}\{−d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}⟨l_{i}, μ⟩\}|_{μ=μ^{∗}} ; g ← ∑_{i=1} l_{ij}[d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}]|_{μ=μ^{∗}} ; θ ← trouveP asAdmissible(f, g, dir) ; μ_{j}^{∗} ← μ_{j}^{∗} + θ × dir_{j} ∀j

jusqu’à ( Critère de convergence satisfait)

Résultats ← μ^{∗}

Algorithme 5 – Gradient_P

Algorithme convex

L’algorithme convexe a été adapter à partir de son homonyme en tomographie par émission. Il consiste à formuler, pour chaque itération, une équation équivalente à la fonction objectif et de trouver le μ correspondant à son maximum. Soit f_{i}(t) = d_{i}e^{−t} + Y_{i}t. Cette équation est strictement convexe. Ainsi, la fonction objectif peut être reformulée comme suit : L(μ) = − ∑_{i} f_{i}(⟨l_{i}, μ⟩). Insérez les idées de De Pierro dans le cas de la tomographie par émission et la procédure est complète :

L(μ) = − ∑_{i} f_{i}(⟨l_{i}, μ⟩)
= − ∑_{i} f_{i}(∑_{j} \frac{l_{ij}μ_{j}^{n}}{⟨l_{i},μ^{n}⟩} \frac{μ_{j}}{μ_{j}^{n}} ⟨l_{i}, μ^{n}⟩)
≥ − ∑_{i} ∑_{j} \frac{l_{ij}μ_{j}^{n}}{⟨l_{i},μ^{n}⟩}f_{i}(\frac{μ_{j}}{μ_{j}^{n}} ⟨l_{i}, μ^{n}⟩)
= Q(μ|μ^{n})

La forme de Q(μ|μ^{n}) garantit que le μ maximisant ce dernier correspond à celui maximisant aussi la fonction objectif. Donc, il ne reste plus qu’à trouver μ qui maximise Q(μ|μ^{n}). Puisque f_{i}(t) est une fonction convexe, le seule point stationnaire de Q(μ|μ^{n}) est un maximum global. Or, il est difficile de résoudre analytiquement l’équation 0 = \frac{∂}{∂μ_{j}} Q(μ|μ^{n}). Les auteurs de [LF95] ont décidé de se servir de la méthode de Newton pour trouver la solution optimale. Notez que cet algorithme ne garantit pas que μ_{j}^{n} ≥ 0. Toutefois, l’argument respecte habituellement cette restriction dans la pratique.

Algorithme convexe (μ, Ȳ, d, l)

{ Donnés : μ_{j} > 0 ∀j; }
\{ Q(μ|μ^{∗}) = − ∑_{ij} \frac{l_{ij}μ_{j}^{∗}}{⟨l_{i},μ^{∗}⟩} \frac{μ_{j}}{μ_{j}^{∗}} ⟨l_{i}, μ^{∗}⟩ \}
{ f_{i}(t) = e^{−t} + Ȳ_{i}t }

répéter

μ_{j}^{∗} ← μ_{j}^{∗} + \frac{∇Q(μ|μ^{∗})_{j}|_{μ=μ^{∗}}}{∇^{2}Q(μ|μ^{∗})_{jj}|_{μ=μ^{∗}}} ∀j jusqu’à ( Critère de convergence satisfait) Résultats ← μ^{∗}

Algorithme 6 – Convex

La modification de l’algorithme convexe pour tenter de diminuer le bruit est la même que celle de l’algorithme EM.

Algorithme convexe pénalisé (μ, Ȳ, l)

{ Donnés : μ_{j} > 0 ∀j; }
{ f_{i}(t) = e^{−t} + Ȳ_{i}t }
\{ Q(μ|μ^{∗}) = − ∑_{ij} \frac{l_{ij}μ_{j}^{∗}}{⟨l_{i},μ^{∗}⟩}f_{i}(\frac{μ_{j}}{μ_{j}^{∗}} ⟨l_{i}, μ^{∗}⟩) \}
\{ V (μ|μ^{∗}) = \frac{γ}{2} (∑_{j} ∑_{k}(w_{jk}ϕ(μ_{j} − μ_{k}^{∗} − μ_{j}^{∗}))) + \frac{γ}{2} (∑_{j} ∑_{k}(w_{jk}ϕ(μ_{k} − μ_{k}^{∗} − μ_{j}^{∗}))) \}
{ Υ(μ|μ^{∗}) = Q(μ|μ^{∗}) − V (μ|μ^{∗}) }

répéter

μ_{j}^{∗} ← μ_{j}^{∗} − \frac{(∇Υ(μ|μ^{∗})|_{μ=μ^{∗}})_{j}}{(∇^{2}Υ(μ|μ^{∗})|_{μ=μ^{∗}})_{jj}} ∀j jusqu’à ( Critère de convergence satisfait) Résultats ← μ^{∗}

Algorithme 7 – Convex_P

En reconstruction itérative, il existe une classe d’algorithme qui est reconnue pour sa rapidité de reconstruction. Il s’agit des algorithmes utilisant les sous-ensembles ordonnés, notés OS. Le principe est de séparer les projections en quelques groupes et de mettre à jour, à chaque itération, la solution avec chacun des sous-groupes. De manière à diversifier l’information contenue dans un sous-groupe, la séparation angulaire entre ses éléments est maximisée. La mise à jour peut être faite avec n’importe quelle méthode APT. La référence [KB05] rapporte les performances pour différentes grandeurs d’OS^{6}. Toutefois, cette méthode tournera autour de la solution optimale du problème sans l’atteindre puisqu’elle tente de résoudre le sous-problème courant. Malgré tout, les performances des OS méritent d’être analysés plus en détails. Ainsi, l’algorithme convexe et convexe pénalisé tous deux avec 10 OS seront inclus dans les tests. Dans ce document, convexos va être utilisé pour désigné l’utilisation de l’algorithme convex avec 10 OS.

2.4 SG

Les SG vont être appliqués sur deux fonctions objectif, soit L(μ) et G(μ). L(μ) est la fonction utilisée avec les APT. G(μ) est une approximation quadratique de L(μ). La démarche derrière le développement de cette formule se retrouve dans l’appendice A de [SB93].

G(μ) = \frac{1}{2} ∥ Ỹ − lμ ∥_{Σ}^{2} où Σ = diag\{Ȳ_{1}, ..., Ȳ_{J}\}

6. OS vient de Ordered subsets

Gradient conjugué non linéaire

La méthode du gradient conjugué non linéaire est souvent utilisée pour résoudre des problèmes d’optimisation avec des systèmes éparses. Cette technique fait partie de la famille des algorithmes à direction de descente. En particulier, sa direction de descente est composée de l’inverse additif du gradient à laquelle est ajouté une fraction, le Bêta, de la direction appliqué dans l’itération précédente. Pour plus d’information sur cet algorithme, voir [HZ05].

Posons F (μ) la fonction objectif. Alors F (μ) est composé de −L(μ) ou de G(μ) additionné à une fonction barrière. Il suffit d’ajouter la fonction Ψ pour avoir un modèle qui converge vers une image moins bruitée.

Cette procédure n’a pas été construite pour prendre en compte des contraintes. Une fonction barrière est ajoutée à l’objectif afin d’assurer que la solution finale soit dans le domaine réalisable. La seule restriction en tomodensitométrie est qu’un coefficient d’atténuation ne peut pas être négatif. La fonction choisie dans ce projet est la barrière logarithme. Donc, elle a la forme suivante : ρ ∑_{j} ln(μ_{j}). Puisque l’objectif est de minimiser et que logarithme naturel d’un nombre près de 0 est négatif, il faut soustraire cette fonction au modèle. Le rho régularise la force de la barrière. Cette dernière doit être diminuée lorsque la solution optimum du problème courant est obtenu. La solution du problème respectant les contraintes est atteinte lorsque rho est suffisamment petit. Voir [D95] pour une explication détaillée de la méthode. Notez que Bêta doit être remis à 0 lorsque rho est diminué puisque le problème n’est plus le même.

Algorithme gradient conjugué Formule de Polak et Ribière(μ, Ȳ, l)

{ Donnés : μ_{j} > 0 ∀j; }

β ← 0

g ← derfonctObj

répéter dir ← −derfonctObj +β × dir f ← fonctObj θ ← trouveP asAdmissible2(f, g, dir) g_prec ← g μ_{j} ← μ_{j} + θ × dir si (μ est la solution optimale du problème courant ) alors

{ρ ← 0.5ρ \atop β ← 0}

sinon g ← derFonctObj y ← g − g_prec β ← Formule Beta de Pilak et Ribière ou de Hager et Zhang β ← max(0, β)

jusqu’à ( Critère de convergence satisfait)

Résultats ← μ

Algorithme 8 – Grad_Conj

Formule Bêta de Polak et Ribière :

β = \frac{⟨g, y^{t}⟩}{⟨g, g^{t}⟩}

Formule de Hager et Zhang :

β = \frac{1}{⟨y, dir⟩} (y − 2 × dir^{t} \frac{⟨y, y^{t}⟩}{⟨y, dir⟩}) g

L’implémentation de trouvePasAdmissible2 est décrite dans [CMI09]. Elle a été construite spécifiquement dans le but de trouver un pas optimal dans le cas où la fonction objectif inclut une fonction barrière. Les performances de cette méthode peuvent être modifiées par l’utilisateur à l’aide de deux paramètres : le nombre d’itérations et la majoration de la fonction objectif sans la barrière. Un plus grand nombre d’itérations permet de trouver un meilleur estimé de pas optimal. Tandis qu’un choix judicieux de majoration^{7} permet de faire les approximations sur un modèle proche de celui qui doit être maximisé. L’implémentation de cette méthode dans le cadre des tests de ce projet limite le nombre d’itérations de trouvePasAdmissible2 à 6. Ce choix a été motivé en observant l’évolution du pas à chaque itération pour une batterie de test et par la volonté d’utiliser le moins possible la matrice système. Les majorations utilisées et leur motivation peuvent être trouvés à la section 8.1.

L-BFGS-B

Le L-BFGS-B correspond à Limited-memory Broyden Fletcher Goldfarb Shanno Bounded. Elle fait partie de la famille des méthodes quasi-Newton d’optimisation. Il s’agît d’une variante de L-BFGS qui prend en compte les restrictions de la forme a_{i} ≤ x_{i} ≤ b_{i}. Cette méthode est reconnue pour ses bonnes performances dans des problèmes d’optimisation avec un grand nombre de variables. De plus, les résultats présentés dans [HGD11] indique que cet algorithme est efficace dans des problèmes de tomodensitométrie. Les tests ont été fait avec le code présent dans EMAN2^{8}. Pour ce projet, L-BFGS-B garde en mémoire les 13 dernières itérations effectuées.

Posez F (μ) la fonction objectif. Alors F (μ) est composé de −L(μ) ou de G(μ). Il suffit d’ajouter la fonction Ψ(μ) pour avoir un modèle qui converge vers une image moins bruitée.

3 Description des images tests utilisés

Afin de comparer les performances des algorithmes dans différentes situations, les données vont être simulées à partir de 4 objets synthétiques. Les mires de test ont été produites par CTSim^{9} et les sinogrammes ont été créés par simulation à l’aide du programme SpectroSin développé par M. Christian Thibaudeau. La simulation a été produite dans quatre conditions différentes. Soit 100, 1 000 et 10 000 photons envoyés par détecteur par projection en ajoutant un bruit de nature Poissonienne^{10} et une version sans bruit. Ainsi, il est possible d’effectuer des analyses plus intéressantes.

Pour faciliter et uniformiser l’évaluation des critères de comparaison, les objets ont été construits avec quelques propriétés. Tout d’abord, les objets peuvent tous être insérés dans un plan de 20 par 20 mm. Le calcul de certains des critères de comparaison utilisés dans ce projet exige une zone uniforme, notée A, et une mince ligne, notée B, avec un coefficient d’atténuation élevé. Voir la section 8.2 pour observer leur position dans les mires. Les valeurs des coefficients d’atténuation des mires ont été ajustées afin que la zone qui inclut A ait la même amplitude dans les quatre mires. Les mires sont affichées à la sections 8.4.

Le code qui a permis de construire les mires et des images de ceux-ci se retrouvent respectivement aux sections 8.5 et 8.4 respectivement.

7. C’est-à-dire une valeur proche du maximum en valeur absolue de la courbure de la fonction objectif.

8. http :blake.bcm.edueman2doxygen_htmllbfgsb_8cpp_source.html

9. Programme open source de simulation de tomographie, http ://ctsim.org/

10. Le rapport de signal-sur-bruit est donné par N/sigma = N/sqrt(N)

3.1 Shepp-logan

La mire du Shepp-logan, dont la définition peut être trouvée dans [KS88], représente une coupe du crâne humain. À l’intérieur de celui-ci, se trouvent des zones d’atténuation différentes correspondant à des formes typiques de ce que l’algorithme doit être capable de reproduire.

3.2 Shepp-logan facile

La structure de cette mire est identique à la précédente. Cependant, les valeurs des coefficients d’atténuation ont été modifiées afin de faciliter l’observation des constrastes. Ils consistent à augmenter de 25 ou 50 fois la différence entre les coefficients de la structure de base et ceux des autres zones.

3.3 H.G.D.

L’image a été inspirée de la mire utilisée dans [HGD11] afin de comparer les performances de certains algorithmes de reconstruction en tomographie. Elle est composée de plusieurs disques qui varient en taille et en atténuation. Cette mire permet d’évaluer les limites des algorithmes. Elle permet d’évaluer la limite d’amplitude et de surface qu’une zone uniforme doit avoir pour être visible dans la solution.

3.4 Catphan CTP401

Les mires Catphan sont utilisées pour tester et ajuster la qualité des reconstructions. Habituellement, elles servent à tester l’équipement. Après une lecture d’un de leurs catalogues^{11}, il a été décidé d’inclure dans les tests le module CTP401 afin d’obtenir un aperçu de la performance des algorithmes sur une structure reconnue dans le domaine. La mire utilisée est une approximation du module CTP401.

4 Analyses accomplies

L’objectif principal de ce document est de comparer l’efficacité des algorithmes. Or, il existe plus d’une technique pour comparer leurs performances. L’approche choisie dans ce projet est divisée en deux parties. Pour commencer, elle propose d’observer le comportement des méthodes au fil des itérations. Ensuite, elle confronte les résultats des algorithmes afin d’en tirer un classement selon différents critères. Cependant, cette dernière approche inclut quelques difficultés. Tout d’abord, la conduite et les solutions obtenues par les méthodes sont potentiellement sensibles à la configuration de l’acquisition, à la solution initiale, à l’objet analysé et à la qualité des données obtenues. Les analyses de ce projet se limitent à une seule configuration^{12}. Tous les tests sont initialisés avec une image dont toutes les valeurs d’atténuation sont égales à 0.02. La section 3 introduit les 16 problèmes que doivent résoudre les méthodes. Soit 4 mires avec chacun 4 niveaux de bruit. Cet éventail d’expériences sera utilisé pour tenter de présenter des données indépendantes des mires et/ou de la qualité des acquisitions. D’autres obstacles se présentent lorsque les solutions des méthodes sont comparées. Il faut choisir un critère d’arrêt pour les algorithmes et définir une méthode pour évaluer le meilleur résultat. Notez que dans ce travail, la valeur du gamma est de 0.1 pour tous les algorithmes sauf lorsque les OS sont utilisés, où il vaut 0.01. Le gamma a été choisi après des tests effectués sur les APT. La modification apportée pour les OS est motivée par un désir de garder la même force de pénalisation pour tous les algorithmes et le fait que l’utilisation de 10 sous-groupes divise le problème de départ en 10 sous-problèmes. La force du epsilon sera établie à 0.8 et le x_critique correspondra à la plus petite différence d’atténuation non-nulle. Soit 0.05, 0.001, 0.005 et 0.089, respectivement pour Shepp-Logan facile, Shepp-Logan, H.G.D et CTP401.

11. Il est accessible à la page web www.phantomlab.com/documents.php

12. Voir section 8.3.

4.1 Observer l’évolution des solutions au fil des itérations

L’objectif derrière l’observation du comportement des méthodes est de fournir suffisamment d’information afin de pouvoir les classer. En vue d’accomplir cette tâche, cet article propose d’observer trois variables : le MSE^{13} entre la solution à une itération donnée et la mire pixelisée, la LMH^{14} de la ligne B et la variance de la zone A^{15}. Les valeurs affichées sont les moyennes des 16 tests. Les variables ont été enregistrées à toutes les 10 itérations, et ce jusqu’à la 80ième itération. C’est en observant que les APT produisaient des images convenables en 80 itérations et le désir de ne pas faire trop d’itérations qui a motivé le choix cette borne. Notez que dans ce document, une itération est définie comme une mise à jour de la solution courante. Donc, une application de la méthode convexe OS 10 est comptée pour 10 itérations. Afin que les résultats soient semblables entre les algorithmes, ils sont tous compilés aux 10 itérations. Les méthodes vont être séparées en deux graphiques selon s’ils ont un terme de régularisation ou non, et ce, pour chaque variable. Ils vont être respectivement nommés algorithmes pénalisés et algorithmes standards. L’ajout de _P à un nom d’algorithme indique qu’il fait partie des algorithmes pénalisés. Le but est d’éviter de surcharger les graphiques et de permettre une comparaison rapide entre ces deux types d’algorithmes. Les graphiques sont tous d’une forme semblable à la figure 6a. Chaque algorithme est représenté par une ligne brisée d’une couleur indiquée dans la légende. Cette courbe représente la moyenne de l’évolution de la variable en question sur les 16 tests selon les itérations.

Figure 1 – MSE entre la solution courante et la mire selon le nombre d’itérations

13. Mean square error ou erreur quadratique moyenne

14. La largeur à mi-hauteur est évaluée sur une ligne verticale de 18 voxels. Afin de diminuer l’effet du bruit sur l’évaluation du LMH, cette ligne est la moyenne de 22 lignes de même longueur et dont chacune intersecte perpendiculairement la ligne B.

L’intersection divise les 22 lignes en deux parties égales.

15. Voir la section 8.2

4.2 Comparer les solutions obtenues avec un critère d’arrêt

Pour être en mesure de comparer les résultats produits par chacune des méthodes, il faut déterminer un critère d’arrêt. Habituellement, l’objectif principal en tomodensitométrie est de reconstruire une image qui a une résolution satisfaisante. Une approche reconnue pour évaluer la résolution est de mesurer la LMH. De plus, des analyses préliminaires montrent que les algorithmes itératifs commencent par produire la forme générale de la mire^{16}. Donc pour cette batterie de tests, les algorithmes itèrent jusqu’à ce que la LMH de la ligne B soit suffisamment petite. Toutefois, certains algorithmes convergent trop lentement vers une LMH acceptable ou même qu’ils commencent à croître au fil des itérations. Ainsi, les algorithmes ont été fixés à une limite de 200 itérations. Ceci correspond à 2.5 fois plus que ce qui est nécessaire pour les algorithmes APT avant d’obtenir une solution très proche, à vue d’oeil, de leur mire respective. De plus, des expériences préliminaires indiquent que les algorithmes avec un terme de régularisation ont des valeurs de LMH plus élevées. Cette différence est induite par l’effet de lissage du terme de pénalisation. Suite à d’autres tests, les choix de LMH sont fixés à 2 voxels et à 5 voxels pour les algorithmes standard et les algorithmes_P respectivement. Le critère d’arrêt étant choisi, il ne reste plus qu’à choisir les critères de comparaison. L’une des caractéristiques les plus importantes lorsque vient le temps de comparer des méthodes reconstruction est la durée de calcul. Afin que les conclusions soient les plus générales possibles, les données affichées correspondront aux moyennes des 16 cas disponibles. Ensuite, le MSE et la variance dans la zone A seront calculés pour permettre d’évaluer qualitativement les solutions obtenues. Des diagrammes à bandes, tel que présenté en figure 13a, vont être utilisés pour faciliter les comparaisons. La hauteur d’une bande représente la moyenne de la variable pour les tests présentés par cette dernière. Dans la figure 13a, la bande à l’extrémité gauche représente la moyenne des MSE pour les 4 mires lorsque reconstruites en utilisant l’algorithme convexe avec 10 OS, et ce, pour des données non bruitées. Les bandes sont triées en ordre croissant de bruit. Il s’agit des mêmes sous-groupes pour les graphiques sur la variance. Dans le graphique qui représente la durée de reconstruction des algorithmes, la bande à gauche correspond à la version standard de l’algorithme.

Figure 2 – Valeurs moyennes des MSE sur les mires selon les algorithmes standard et le bruit

16. Ce constat n’est pas toujours vrai, et en conséquence, il a fallu demander un minimum d’itérations, 15 dans cette expérience, avant l’arrêt de l’algorithme.

5 Résultats obtenus

Avant d’évaluer la capacité des algorithmes à produire une solution exacte et rapide, il faut observer comment ils se comportent généralement. Préalablement à la lecture des résultats, il faut spécifier certains détails. Les figures nommées dans cette section se retrouvent toutes dans la section 8.7. Les algorithmes vont être appelés par les noms utilisés dans les légendes des figures de cette section. Les notations "FObj1 " ou "*" font référence à la fonction objectif présentée dans la section 2.4.

5.1 Comportement des algorithmes sur 80 itérations

MSE sur 80 itérations

Les figures 6a et 7a montrent que les solutions courantes des variantes du GC commencent loin de la solution de départ et peinent à rattraper les autres méthodes. Cette observation n’est pas sans fondement. Ces algorithmes sont les seuls qui doivent inclure un terme de barrière. Or, lors des tests, il n’y avait pas suffisamment d’information pour gérer la force de la barrière. Donc, il est possible qu’il existe un meilleur choix de ρ_{0} qui améliorerait leur performance. De plus, il est possible de constater que les performances sont pires lorsqu’un terme de régularisation est ajouté. Cette constatation peut être expliquée par le choix de majoration qui n’est pas optimale. Puisqu’il est difficile de voir l’évolution des autres algorithmes dans ces figures, une figure correspondante, avec une échelle logarithmique a été jointe. La figure 6b montre une augmentation du MSE des algorithmes standards ce qui indique qu’ils s’éloignent des mires^{17}. Ce résultat est attendu puisque la solution optimale d’un sinogramme bruité est une image bruitée. L’ajout d’un terme de régularisation permet de corriger cette divergence (Figure 7a). Remarquez que la courbe de convex n’est pas visible dans ces figures. Cela est dû au fait qu’elle est trop près de la courbe de convexos10. Ainsi le fait d’utiliser environ 10 fois moins d’information par itération a permis d’obtenir des résultats quasi identiques à la méthode utilisant la totalité des données. On note aussi dans la figure 6b que les variantes de L-BFGS-B se comportent correctement au début, mais divergent de la solution voulue dès la 20ième itération. Cette procédure est reconnue pour sa grande vitesse de convergence, ce qui n’est pas nécessairement voulue lorsque les données sont bruitées. La figure 7a montre que l’ajout d’un terme de régularisation rend ses performances équivalentes aux APT. En résumé, les algorithmes qui se comportent le mieux pour le MSE, avec ou sans fonction pénalisation, sont convexos10 et convex.

LMH sur 80 itérations

Avec les figures 8a et 9a, il est possible de voir que l’évolution de la LMH selon les itérations. Bien que la LMH semble être quasi indépendant du terme de régularisation, il converge vers des valeurs différentes, soit 2 et 5 voxels respectivement pour tous les algorithmes sauf les GC. D’ailleurs, le comportement des variantes du GC dans les 30 premières itérations est complètement différent de ceux des autres algorithmes. Les hypothèses retenues, après analyse des images à chaque 10 itérations, sont que la fonction barrière est trop forte au départ ou que sa force n’est pas suffisamment diminuée au fil des itérations. Ces hypothèses sont motivées par le fait que les premières estimées de ces algorithmes ont des coefficients d’atténuation élevés dans des zones où ils devraient être nuls. Pour être en mesure de départager les procédures selon la LMH, il faut regarder les figures 8b et 9b qui représentent un zoom des figures précédentes sur les meilleurs algorithmes. Tout d’abord, notez que gradient et gradient_P semblent être meilleur que leurs confrères. Toutefois, cette avance semble être favorisée par la limite de 80 itérations puisque L-BFGS-B et L-BFGS-BFObj1 semblent être en mesure de le rattraper. De plus, les méthodes convex et convexos10 n’ont pas eu le même genre de départ et convexeos10 a toujours un peu d’avance sur convex. Pour finir, les APT semblent tous avoir un comportement semblable, et ce, dans le cas avec et sans le terme de régularisation. Cette constatation est logique avec le fait qu’ils ont tous plus ou moins été construits sur les mêmes bases. En outre, les variantes de L-BFGS-B gèrent plus difficilement l’ajout d’un terme régulateur que les APT. Considérant la vitesse de convergence de ce dernier, il faudrait peut-être songer à diminuer la force de la fonction de pénalisation lorsqu’il a atteint une certaine valeur de la fonction objectif.

17. Ce constat n’est pas vrai pour les variantes du gradient conjugué. Cette contradiction provient surement du fait qu’elles ne sont pas encore suffisamment près de leur solution.

Variance sur 80 itérations

L’analyse de l’évolution de la variance de la zone B^{18} selon les itérations confirme certaines des conclusions obtenues en observant le MSE. La figure 10a montre que lorsque les algorithmes itèrent trop longtemps, ils commencent à inclure du bruit dans la solution. De plus, elle appuie le fait que L-BFGS-B et L-BFGS-BFObj1 débutent tôt à reproduire le bruit. La figure 11a, pour sa part, apporte un peu plus d’information. Ce graphique semble indiquer que les APT ont plus tendance que les SG à reconstruire le bruit localement, ce qui ne se voyait pas dans les graphiques de MSE. Pour finir, la figure 10b, qui est un zoom sur les lignes brisées les plus basses de la figure 10a, semble indiquer que les variantes du GC ont un contrôle acceptable du bruit localement lorsqu’il n’inclut pas un terme de régularisation. Cette dernière remarque supporte l’idée que ces algorithmes ont du potentiel et qu’une meilleure gestion de l’interaction entre la fonction barrière et de pénalisation, suivie d’une meilleure majoration, pourraient améliorer leurs performances.

5.2 Étalon des solutions obtenues à l’arrêt

Variance avec critère d’arrêt sur FWHM

L’idée de cette sous-section est de présenter le meilleur algorithme lorsque le but est d’obtenir des solutions dont les zones uniformes ont le moins de bruit possible. Si les figures 12a et 12b sont utilisées pour comparer les méthodes, il en ressort que les algorithmes trml_P et GC_P sont les plus performants. En outre, les algorithmes avec un terme régulateur se comportent mieux que leur confrère. Toutefois, il faut savoir que GC_P a de la difficulté à atteindre le critère d’arrêt. Donc, il est possible qu’il n’ait simplement pas convergé suffisamment pour reconstruire le bruit. Dans une situation où utiliser un algorithme avec un terme de pénalisation est déconseillé, notez que l’algorithme trml_P est le chef de file dans cette sous-catégorie.

MSE avec critère d’arrêt sur LMH

Maintenant, les algorithmes seront classés selon le MSE entre la solution à l’arrêt de l’algorithme et la mire originale. En observant les figures 13a et 13b, il est possible de constater que les résultats sont beaucoup plus stables avec les méthodes pénalisées, mais que les résultats dans les cas peu bruités sont meilleurs avec les procédures standard. Donc, si le niveau de bruit est relativement faible, il est avantageux de prendre une méthode standard. Plus particulièrement, la figure 13a montre qu’utiliser L-BFGS-B ou L-BFGS-BFOb1, lorsque les données sont prises dans des bonnes conditions, est un choix judicieux. Cependant, si des solutions stables sont recherchées, mieux vaut prendre des algorithmes pénalisés. En outre, convexos_P a obtenue des résultats supérieurs et plus stables que ses confrères pénalisés. Toutefois, il ne faut pas écarter L-BFGS-B_P et L-BFGS-BFOb1_P qui sont très proche de celui-ci. Afin d’appuyer les résultats précédents, la section 8.6 montre des images obtenues avec certains algorithmes pour la mire Shepp-Logan-easy lorsque la source envoie théoriquement 1000 photons. Les reconstructions obtenues par les algorithmes pénalisés sont flous ce qui indique que le gamma choisi était un peu trop fort.

Temps de calcul avec critère d’arrêt sur LMH

Le dernier critère de comparaison utilisé dans ce projet est la durée de reconstruction. Les données obtenues pour chacun des algorithmes sont illustrées dans la figure 14a. La figure montre que convexos se distingue des autres méthodes. Ensuite, il y a les variantes de L-BFGS-B suivi de près par convex et Gradient. Notez que la durée affichée n’a pas été pénalisée lorsque la méthode est arrêtée avant d’atteindre la LMH visé. Toutefois, les variantes de GC sont celles qui ont le moins réussi à atteindre cette cible, ce qui ne change pas le classement actuel. L-BFGS-B_P et L-BFGS-BFOb1_P ont eue, eux aussi, quelques problèmes à arriver à l’objectif. Donc, leurs valeurs sont sous-évaluées dans cette figure. La figure 14b n’inclut que les tests qui ont atteint la LMH critique. Cette figure montre seulement que les variantes du GC utilisées dans ce travail ont de la difficulté à converger et que les algorithmes pénalisés n’atteignent pas toujours le critère d’arrêt. Cette dernière remarque pourrait provenir de la force du gamma.

18. Voir section 8.2

Conclusion

Pour terminer, un rappel sur certains des résultats obtenus. Pour commencer, l’utilisation de méthodes standards devrait être limitée aux données qui sont peu contaminées par du bruit puisque la qualité de la reconstruction est très affectée par la présence de bruit. Sinon il vaut mieux utiliser les méthodes robustes que sont les algorithmes pénalisés. De plus, malgré les réserves faites sur la convergence des APT utilisant des OS, les figures de ce projet semblent indiquer que ces algorithmes ont un comportement égal ou même supérieur aux autres méthodes. Ensuite, les résultats de ce projet démontrent encore que les variantes de L-BFGS-B sont très efficaces en tomodensitométrie. Toutefois, les résultats de ce document déconseillent de l’utiliser sur des données relativement bruit ées sans terme de pénalisation. Enfin, les résultats de l’implantation de variantes de gradient conjugué n’ont pas été fructueux. Or, ils indiquent que cette méthode possède du potentiel. Il suffirait de la retravailler sur les points notés précédemment.

6 Idée à explorer

Voici quelques idées qui seraient intéressantes d’explorer.

6.1 Les paramètres de la fonction de pénalisation

Le gamma et l’epsilon présents dans la pénalisation déterminent l’efficacité de cette dernière. Dans ce travail, le choix du gamma a été fixé à une valeur choisie a posteriori. Or, il serait intéressant de voir si ce paramètre a une influence différente sur les performances de reconstruction selon les algorithmes. De plus, la formule de l’epsilon utilisé par l’auteur est construite dans l’idée de garder le plus petit détail d’intérêt dans un objet. Cependant, il faudrait démontrer l’efficacité de cette dernière ou trouver une meilleure méthode pour atteindre cet objectif.

6.2 La majoration de la courbure

Les majorations utilisées pour trouvePasAdmissible2 ne sont pas optimales. C’est-à-dire qu’il existe des plus petites majorations. Le fait que ∇^{2}h(θ) soit toujours positif dans son domaine et que cette dernière soit une fonction convexe permet de trouver une meilleure majoration, plus particulièrement dans le cas de G(μ). Améliorer ces majorations, surtout dans les cas avec pénalisation, permettrait surement d’obtenir des résultats satisfaisants en utilisant les variantes de gradient conjugué non linéaire.

6.3 Optimisation séquentielle de sous espaces

Cet algorithme est reconnu pour sa grande vitesse de convergence dans des problèmes tel la tomographie. De plus, cette méthode a été testée dans ce domaine et il semblerait que les résultats obtenus soient très satisfaisants pour une durée de reconstruction courte. Bref, il s’agit d’une procédure qui mérite d’être prise en compte.

6.4 Approfondir l’expérience

Il y a beaucoup d’avenues qui méritent d’être étudiées. Par exemple, voir si les conclusions sont semblables dans la reconstruction d’objet en 3 dimensions. Ensuite, les figures montrent que certains des algorithmes débutent lentement. Cependant, une fois suffisamment proche de la solution, elle converge rapidement. Or, l’inverse de Radon produit un bon estimé de la solution et à moindre coût. Donc, il serait intéressant de refaire les analyses avec cette solution initiale et de voir si le classement des algorithmes est modifié.

7 Remerciements

Cet article provient d’un travail qui a été financé par une bourse offerte par la Faculté de médecine et des sciences de la santé de l’Université de Sherbrooke (FMSS), les Instituts de recherche en santé du Canada (IRSC) et le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG). L’auteur tient à remercier M. Roger Lecomte pour lui avoir permis de travailler sur ce projet, M. Jean-pierre Dussault pour les discussions qui ont motivé l’ajout des SG et M. Christian Thibaudeau pour son aide dans le projet. Les critiques de ceux-ci ont permis d’améliorer la rédaction de ce rapport technique.

8 Annexe

8.1 Majoration utilisée dans pasAdmissible2

La courbure d’un graphe, h(θ) est évalué avec la formule suivante :

|\frac{∇^{2}h(θ)}{(1 + (∇h(θ))^{2})^{1.5}} |

L’algorithme trouvePasAdmissible2 demande une valeur majorant cette expression pour θ ∈ [0, ᾱ[. Ces bornes proviennent du fait que la valeur ᾱ représente le plus grand theta tel que μ_{j} + θ ∗ dir_{j} ≥ ϵ ∀_{j} et qu’un pas dans une direction ne peut pas être négatif. Notez que h(θ) = F (μ+θ∗dir) où F (x) est la fonction à optimisé sans la barrière logarithme. Dans ce projet, ces fonctions sont L(μ) et G(μ).

Notez que −∇^{2}L(μ) > 0∀μ et que ∇^{2}G(μ) > 0∀μ Aisni, (1 + (∇h(θ))^{2})^{1.5} ≥ 1 ⇒ \frac{∇^{2}h(θ)}{(1+(∇h(θ))^{2})^{1.5}} ≤ ∇^{2}h(θ)

Rappel : L(μ) = ∑_{i=1}\{−d_{i}e^{−⟨l_{i},μ⟩} − Y_{i}⟨l_{i}, μ⟩\} et G(μ) = \frac{1}{2} ∥ Ỹ − lμ ∥_{Σ}^{2}.

Majoration avec L(μ)

Soit h(θ) = L(μ + θ ∗ dir) alors ∇^{2}h(θ) = ∑_{i}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−θ⟨l_{i},dir_{i}⟩}]

∇^{2}h(θ) = ∑_{i}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−θ⟨l_{i},dir_{i}⟩}]
= ∑_{i|dir_{i}≥0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−θ⟨l_{i},dir_{i}⟩}] + ∑_{i|dir_{i}<0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−θ⟨l_{i},dir_{i}⟩}]
≤ ∑_{i|dir_{i}≥0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}] + ∑_{i|dir_{i}<0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−θ⟨l_{i},dir_{i}⟩}]
≤ ∑_{i|dir_{i}≥0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}] + ∑_{i|dir_{i}<0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−ᾱ⟨l_{i},dir_{i}⟩}]

Donc, K = ∑_{i|dir_{i}≥0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}] + ∑_{i|dir_{i}<0}[⟨l_{i}, dir_{i}⟩^{2}e^{−⟨l_{i},μ⟩}e^{−ᾱ⟨l_{i},dir_{i}⟩}] pour h(θ) = L(μ + θ ∗ dir)

Majoration avec G(μ)

Soit h(θ) = G(μ + θ ∗ dir) alors ∇^{2}h(θ) = ∑_{i}[⟨l_{i}, dir_{i}⟩^{2}Y_{i}]

Donc, K = ∑_{i}[⟨l_{i}, dir_{i}⟩^{2}Y_{i}] pour h(θ) = G(μ + θ ∗ dir)

Majoration avec l’ajout de la fonction de pénalisation

Soit Ψ(μ) = γ (∑_{j}^{J} ∑_{k}^{K_{j}} κ_{j}κ_{k}w_{jk}ϕ(μ_{j} − μ_{k})). Notez que

∇_{θ}^{2}Ψ(μ + θ ∗ dir) = γ (∑_{j}^{J} ∑_{k}^{K_{j}} κ_{j}κ_{k}w_{jk} \frac{ϵ(dir_{k} − dir_{k_{j}})^{2}}{(ϵ + [(μ_{k} + θ ∗ dir_{k}) − (μ_{k_{j}} + θ ∗ dir_{k_{j}})])^{1.5}} )

On remarque que le domaine de cette fonction est stricement positif. Donc, il suffit de trouver une valeur, K_{p} qui majorise cette expression pour obtenir une majoration des courbures des fonctions pénalisés La majoration des fonctions pénalisés utilisés dans ce projet est K + K_{p} avec K correspondant à la fonction objectif courante.. La majoration choisie est K_{p} = γ (∑_{j}^{J} ∑_{k}^{K_{j}} κ_{j}κ_{k}w_{jk} \frac{(dir_{k}−dir_{k_{j}})^{2}}{(ϵ)^{0.5}} )

8.2 Zone d’importance pour l’évaluation du LMH et de la variance

L’image qui suit est la superposition des 4 mires à laquelle a été ajouté la position de la zone A et celle de la ligne B.

Figure 3 – Zone utilisée pour les tests

8.3 Configuration d’acquisition

Dimension de l’image  20cm x 20cm
Discrétisation  256 x 256
Mode d’acquisition  Géométrie à angles égaux
Nombre de détecteurs  520
Nombre de projections  360
Largeur angulaire du panneau des détecteurs 60^{∘}
Distance entre la source et le point de rotation 28.2843 cm
Rayon de courbure des détecteurs  56.5686 cm

8.4 Mires

(a) Shepp-Logan  (b) Shepp-Logan-easy
(c) H.G.D.  (d) CTP401

Figure 4 – Mires : Noir = 0 et blanc = 2.5

8.5 Code des mires

Une ligne de code représente un objet en CTsim. Signification du code : OBJET PositionX PositionY HauteurX HauteurY AngleRotation CoeffD’aténuation

H.G.D. :

Structure de base :

ELLIPSE 0 0 10 10 0 0.1

Cercle extérieur ELLIPSE 6.6 0 0.55 0.55 0 0.005 ELLIPSE 3.3 5.5 0.55 0.55 0 0.020 ELLIPSE −3.3 5.5 0.55 0.55 0 0.035 ELLIPSE −6.6 0 0.55 0.55 0 0.050 ELLIPSE −3.3 −5.5 0.55 0.55 0 0.065 ELLIPSE 3.3 −5.5 0.55 0.55 0 0.080

Cercle du milieu ELLIPSE 2.2 1.1 0.033 0.033 0 0.150 ELLIPSE 1.1 2.2 0.044 0.044 0 0.150 ELLIPSE −1.1 2.2 0.055 0.055 0 0.150 ELLIPSE −2.2 1.1 0.066 0.066 0 0.150 ELLIPSE −2.2 −1.1 0.077 0.077 0 0.150 ELLIPSE −1.1 −2.2 0.088 0.088 0 0.150 ELLIPSE 1.1 −2.2 0.099 0.099 0 0.150 ELLIPSE 2.2 −1.1 0.11 0.11 0 0.150

Cercle intérieur ELLIPSE 1.1 0 0.0275 0.0275 0 0.3 ELLIPSE 0 1.1 0.00198 0.00198 0 0.3 ELLIPSE −1.1 0 0.0022 0.0022 0 0.3 ELLIPSE 0 −1.1 0.0242 0.0242 0 0.3

Ligne A ELLIPSE 0 4 1.5 0.02 0 0.4

CTP401 :

Structure de base :

ELLIPSE 0 0 10 10 0 0.12

Segment droite :

RECTANGLE 5.5 1.5 0.025 1.4 0 0.2 RECTANGLE 5.5 −1.5 0.025 1.4 0 0.2

Segment gauche :

RECTANGLE −5.5 1.5 0.025 1.4 0 0.2 RECTANGLE −5.5 −1.5 0.025 1.4 0 0.2

Ligne A :

ELLIPSE 0 4 0.02 3 90 0.4

Segment bast :

RECTANGLE 1.5 −5.5 0.025 1.4 90 0.2 RECTANGLE −1.5 −5.5 0.025 1.4 90 0.2

Cercle extérieur ds segments : ELLIPSE 8 0 1 1 0 0.089 ELLIPSE 0 8 1 1 0 0.157 ELLIPSE −8 0 1 1 0 0.436 ELLIPSE 0 −8 1 1 0 −0.12 Cercle entre deux segments :

ELLIPSE 4 4 0.25 0.25 0 0.436 ELLIPSE 4 −4 0.25 0.25 0 −0.12 ELLIPSE −4 −4 0.25 0.25 0 0.436 ELLIPSE −4 4 0.25 0.25 0 −0.12

Cercle qui forment un pentagone : ELLIPSE 0 −2 0.8 0.8 0 0.157 ELLIPSE 1.90 −0.62 0.6 0.6 0 0.157 ELLIPSE 1.18 1.62 0.4 0.4 0 0.157 ELLIPSE −1.18 1.62 0.2 0.2 0 0.157 ELLIPSE −1.90 −0.62 0.1 0.1 0 0.157

Shepp-easy :

Structure de base : ELLIPSE 0 0 10 7 90 0.2 ELLIPSE 0 −0.0184 9.46 6.624 90 −0.1

Grande zone ELLIPSE 2.2 0 3.1 1.1 72 −0.05 ELLIPSE −2.2 0 4.1 1.6 108 −0.05 ELLIPSE 0 3.5 2.5 2.1 90 0.005 ELLIPSE 5.538 −3.858 0.33 2.06 −18 0.075

Petite zone ELLIPSE −0.8 −6.05 0.46 0.23 0 0.05 ELLIPSE 0 1 0.46 0.46 0 0.05 ELLIPSE 0 −1 0.46 0.46 0 0.05 ELLIPSE 0 −6.05 0.23 0.23 0 0.05 ELLIPSE 0.6 −6.05 0.46 0.23 90 0.05

Ligne A ELLIPSE 0 4 1.5 0.02 0 0.4

Shepp-logan :

Structure de base : ELLIPSE 0 0 10 7 90 0.2 ELLIPSE 0 −0.0184 9.46 6.624 90 −0.1

Grande zone ELLIPSE 2.2 0 3.1 1.1 72 −0.002 ELLIPSE −2.2 0 4.1 1.6 108 −0.002 ELLIPSE 0 3.5 2.5 2.1 90 0.001 ELLIPSE 5.538 −3.858 0.33 2.06 −18 0.0015

Petite zone ELLIPSE −0.8 −6.05 0.46 0.23 0 0.001 ELLIPSE 0 1 0.46 0.46 0 0.002 ELLIPSE 0 −1 0.46 0.46 0 0.001 ELLIPSE 0 −6.05 0.23 0.23 0 0.001 ELLIPSE 0.6 −6.05 0.46 0.23 90 0.001

Ligne A ELLIPSE 0 4 1.5 0.02 0 0.4

8.6 Comparaison d’images obtenues avec quelques algorithmes

(a) ConvexOs_10  (b) Convex  (c) L-BFGS-B  (d) trlm

(e) ConvexPOs_10  (f) Convex_P  (g) L-BFGS-B_P  (h) trml_P

Figure 5 – Images obtenues pour différents algorithmes avec la mire Shepp-Logan-easy et en utilisant 1000 photons. Les algorithmes standards et pénalisés ont été arrêtés lorsque la ligne B de la solution courante avait un LMH de 2 et de 5 respectivement. Un pixel noir et un pixel blanc représentent respectivement un coefficient d’atténuation de 0 et de 2.5

8.7 Figures utilisées dans l’analyse

Soit MSE_{k}^{(i,j)} le MSE entre la mire et la solution à l’itération k d’un algorithme appliqué sur un sinogramme de la mire i avec un niveau de bruit j^{19}. Soit Var_{k}^{(i,j)} et LMH_{k}^{(i,j)} respectivement la variance de la zone A et la largeur à mi-hauteur de la ligne B^{20} de la solution à l’itération k d’un algorithme appliqué sur un sinogramme de la mire i avec un niveau de bruit j. Alors, MSE_{k}^{∗}, Var_{k}^{∗} et LMH_{k}^{∗} sont respectivement les moyennes des MSE_{k}^{(i,j)}, des Var_{k}^{(i,j)} et des LMH_{k}^{(i,j)} pour les 4 mires avec chacun des 4 niveaux de bruit.

Soit MSE^{(i,j)} et Var^{(i,j)} respectivement MSE_{k^{∗}}^{(i,j)} et Var_{k^{∗}}^{(i,j)} où k^{∗} représente l’itération à laquelle la solution courante a atteint le critère d’arrêt. Soit Recon_Time^{(i,j)} le temps nécessaire pour atteindre le critère d’arrêt d’un algorithme appliqué sur un sinogramme de la mire i avec un niveau de bruit j.

Alors, MSE^{∗,j} et Var^{∗,j} sont respectivement la moyenne des MSE(i, j) et la moyenne des Var^{(i,j)} sur les 4 mires à un niveau de bruit j. L’ordre, de gauche à droite, des niveaux de bruit est noNoise, 10000c, 1000c et 100c^{21}. Recon_Time^{∗} est la moyenne des Recon_Time^{(i,j)} pour les 4 mires avec chacun des 4 niveaux de bruit.

19. Voir section 3

20. Voir section 8.2

21. Voir la section 3

(a) échelle standard  (b) échelle log

Figure 6 – MSE_{k}^{∗} à chaque 10 itérations pour chacun des algorithmes standards

(a) échelle standard  (b) échelle log

Figure 7 – MSE_{k}^{∗} à chaque 10 itérations pour chacun des algorithmes pénalisés.

(a) échelle standard  (b) Zoom in

Figure 8 – LMH_{k}^{∗} à chaq ue 10 itérations pour chacun des algorithmes standards

(a) échelle standard  (b) Zoom in

Figure 9 – LMH_{k}^{∗} à chaque 10 itérations pour chacun des algorithmes pénalisés

(a) échelle standard  (b) Zoom in

Figure 10 – Var_{k}^{∗} à chaque 10 itérations pour chacun des algorithmes standards

(a) échelle standard

Figure 11 – Var_{k}^{∗} à chaque 10 itérations pour chacun des algorithmes pénalisés

(a) LMH critique : 2  (b) LMH critique : 5

Figure 12 – Var^{∗,j} selon les algorithmes pénalisés et le bruit représenté avec une échelle log. Les algorithmes sont arrêtés lorsque la solution courante a atteint la LMH critique ou lorsque le nombre d’itérations dépasse 200. Les couleurs de gauche à droite représentent respectivement noNoise, 10000c, 1000c et 100c.

(a) LMH critique : 2  (b) LMH critique : 5

Figure 13 – MSE^{∗,j} selon les algorithmes pénalisés et le bruit représenté avec une échelle log. Les algorithmes sont arrêtés lorsque la solution courante a atteint la LMH critique ou lorsque le nombre d’itérations dépasse 200. Les couleurs de gauche à droite représentent respectivement noNoise, 10000c, 1000c et 100c.

(a) Recon_Time^{∗} sur tous les tests  (b) Recon_Time^{∗} sur les tests qui ont atteint la LMH critique

Figure 14 – Recon_Time^{∗} selon les algorithmes et l’inclusion d’un terme de régularisation. Les algorithmes ont été arrêtés lorsque la LMH de la solution courante atteint la valeur critique ou lorsque le nombre d’itérations dépasse 200. La valeur critique pour les algorithmes standards(Bleu) et les algorithmes pénalisés(Cyan) sont respectivement de 2 et de 5.

Références

[CMI09] É. Chouzenoux, S. Moussaoui & J. Idier, "A majorize-minimize line search algorithm for barrier function", version 4, 31 p., 2009.

[D95] J-P. Dussault, "Numerical stability and efficiency of penalty algorithms", SIAM Journal on Numerical Analysis, Vol 32, No. 1, pp. 296-317, 1995.

[FR96] J.A. Fessler & W.L. Rogers, "Resolution properties of regularized image reconstruction methods", Technical Report No. 297, 25 p., 1996.

[H80] G.T. Herman, "Image reconstruction from projections : The fundamentals of computerized tomography", IEEE Transactions on Signal Processing, vol. 41, no 2, pp. 534-548, 1980.

[HGD11] B. Hamelin, Y. Goussard & J-P. Dussault, "Comparison of optimization techniques for regularized statistical reconstruction in x-ray tomography", IEEE IPTA, pp. 87-91, 2011.

[HZ05] W.H. Hager & H. Zhang, "A survey of nonlinear conjugate gradient methods", 21 p. , 2005.

[KB05] J.S. Kole & F.J. Beekman, "Evaluation of the ordered subset convex algorithm for cone-beam CT", Physics in Medicine and Biology, pp. 613-623, 2005.

[KS88] A.C. Kak & M. Slaney, "Principles of Computerized Tomographic Imaging", IEEE, 327 p., 1988.

[LC84] K. Lange & R. Carson, "EM reconstruction algorithms for emission and transmission tomography", Journal of Computer Assisted Tomography pp. 306-316 1984.

[LF95] K. Lange & J.A. Fessler, "Globally convergent algorithms for maximum a posteriori transmission tomography", IEEE Trans on Image Processing pp. 1430-1438, 1995.

[RP86] J. Radon & P.C. Parks, (translator), "On the determination of functions from their integral values along certain manifolds", IEEE Transactions on Medical Imaging 5, pp. 170–176, 1986.

[SB93] K. Sauer, & C. Bouman , "A local update strategy for iterative reconstruction from projections," IEEE Transactions on Signal Processing, Vol. 41, no. 2, pp. 534-548, Feb. 1993.