Keywords: , teeth, etc

optimization by using the VBA algorithm . . . . . .

Table des matières

The original object f j the j-th voxel in object f f Estimation of ojbect

f ∈ R N ×1 Gradient of f f -r
Vector f except the element f r F (u, v)

2D Fourier transformation of f (x, y)

F 1D 1D Fourier Transformation F 2D 2D Fourier Transformation F -1 1D 1D inverse Fourier Transformation F -1 2D 2D inverse Fourier Transformation g(r, φ) the observed projection data g ∈ R M ×1
The projection sinogram g i

The data observed by the i-th detector

g 0 ∈ R M ×1 Projection data without noise g ∈ R M ×1
Gradient of projection g G(ρ, φ)

1D Fourier Transformation of g(r, φ) G(•)

Gamma distribution g 0 (n), g 1 (n)

Low-pass/high-pass filter pair for lower filter bank of DT-CWT

H ∈ R M ×N projection operator H ij
The length of the i-th projection line in pixel j h 0 (n), h 1 (n)

Low-pass/high-pass filter pair for upper filter bank of DT-CWT I 0 incident beam intensity I attenuated radiation intensity

I n ∈ R n×n Identity matrix J(•) Optimization criterion J(f ) Optimization criterion on f J(•) Optimization criterion ∇J(•)
Derivative of criterion J(•) l trajectory of a radiation Chap 4 : Level of Haar transformation

L(•) Laplacian distribution m ∈ R N ×1
Mean of distribution of Markovian model for f N r (•)

Set of the neighbour voxels N (•)

Normal or Gaussian distribution p(•)

Probability distribution function q(•)

Probability distribution function Q(g, f )

Discrepancy of observed data q ∈ R N ×1

Normalized gradient of object f r distance of l from the origin r Current voxel in the Markovian model r Neighbour voxels in teh Markovian model r ∈ [1, l + 1] Rank of coefficent in l-level

Haar transformation R(f ) Regularization term on f R Radon transform R # Dual Radon transform St(•) Student-t distribution St g (•) Generalized Student-t distribution v ∈ R M ×1
Variance of the additive noise V ∈ R M ×M Covariance matrix of the additive noise

v f ∈ R N ×1 Variance of f V f ∈ R N ×N Covariance matrix of f v z ∈ R N ×1 Variance of z V z ∈ R N ×N Covariance matrix of z v ξ ∈ R N ×1 Variance of ξ V ξ ∈ R N ×N Covariance matrix of ξ v ρ ∈ R M ×1 Variance of ρ V ρ ∈ R N ×N Covariance matrix of ρ v τ ∈ R M ×1 Variance of τ V τ ∈ R N ×N Covariance matrix of τ v c ∈ R M ×1
Variance of f

V c ∈ R N ×N Covariance matrix of f z ∈ R M ×1
Hidden variable Noise of projection of the gradient f 3

1

Résumé en franc ¸ais

Motivation et contexte

Cette thèse est une contribution au problème de reconstruction de tomographie en 3D. Le rayonnement X après sa découverte par Wilhelm Rontgen a été rapidement utilisé dans le domaine de la tomographie. Au début, on utilise uniquement les images radio permettant de voir la projection de l'objet traversé par les rayons X sur la plaque radio. Mais la projection des rayons X dans une seule direction perdait des informations, notamment l'information de profondeur.En 1917, Radon propose un algorithme afin de reconstruire les images, coupes 2D du patient mais qui ne pouvait etre appliqué en pratique avant les années 1970 et l'arrivée des scanners CT (Computed Tomography). Cette technique d'acquisition avec les projections appliquées dans différentes directions est développée avec de nombreuses publications sur les méthodes analytiques et itératives. Dans ces techniques, les données mesurées à différents angles sont utilisées pour reconstuire le volume 3D. À la fin des années 1970, les méthodes de reconstruction à faisceau en éventail (fan beam) ont été développées. La technique de tomographie en rayons X a d'abord été utilisée dans les applications médicales rapidement après son invention puis a été utilisée dans des nombreux autres domaines : la science marine, les géosciences, les applications industrielles, etc.

La base de la technique CT de rayons X est basée sur l'atténuation des rayons X en passant par un objet. Une vue de dessus d'un système de projection de rayons X est illustrée dans la figure 2.1. Différents matériaux avec une densité différente dans l'objet ont un taux d'atténuation différent pour le rayon X. Cette atténuation est définie par la loi de Beer-Lambert :

I = I 0 e -f (l) dl , (1.1)
où l est le rayon X passant par un objet, I 0 est l'intensité du faisceau incident, f (l) est le coefficient d'atténuation linéaire local sur le rayonnement l et I est l'intensité de rayonnement atténuée.

CHAPITRE 1. R ÉSUM É EN FRANC ¸AIS

L'un des défis de la CT est de diminuer le nombre de projections. L'absence de données complètes est justifiée par des raisons différentes selon les applications. Le rayon X a été utilisé pour la première fois dans les applications médicales depuis son introduction dans les années 1970. Il a été utilisé dans le diagnostic de la tête, des poumons, de l'angiographie pulmonaire, du cardiaque, de l'abdomen et du pelvis, des extrémités, etc. En utilisant la tomographie en rayon X, le diagnostic des organes peut être effectué avec une haute résolution afin de distinguer plus de détails. Cependant, cette technique présente de nombreux inconvénients. Le plus important est le risque pour les patients. Le rayon X utilisé peut endommager les cellules du corps, y compris les molécules d'ADN, qui peuvent conduire à cancer, ou à un autre effet secondaire indésirable. Par conséquent, les efforts de recherche récents se concentrent sur les techniques de reconstruction CT avec une dose plus faible de rayonnement, également appelée ALARA (As Low As Reasonably Achievable) et moins de temps d'exposition pour chaque projection. Dans les applications industrielles, le CT à rayons X est souvent utilisé dans le Contrôle Non Destructif (CND). Il est utilisé pour produire des représentations 3D de composants à l'externe et à l'interne, pour la détection des défauts, l'analyse des pannes, etc. Dans les applications industrielles, la dose de rayonnement n'est pas nécessairement faible. Le défi principal de l'application industrielle est la détection dynamique. Dans certains cas, l'objet examiné n'est pas statique, par exemple l'objet est imagé avec un liquide qui s'écoule. Dans certains autres examens industriels, les objets sont placés sur un convoyeur. Dans ces cas, la projection devrait se faire en peu de temps, donc le nombre de projections doit être réduit. Par conséquent, les reconstructions devraient être réalisées avec moins de projections. Dans l'imagerie médicale et le CND industriel, un autre problème auquel nous pouvons faire face est l'angle limité d'acquisition. Dans les cas général, le scanner CT teste un volume en obtenant des rayonnements en projection pour des angles uniformément répartis entre [0 • , 180 • ]. Mais dans certaines applications, la couverture angulaire totale n'a pas pu être obtenue. Un exemple dans la demande médicale est la tomographie dentaire où toutes les projections ne sont pas faites en raison d'encombrement du tomographe. Dans les applications industrielles, ce problème apparaît également lorsque le composant testé a une grande taille et est supérieure à la portée du scanner, on obtient alors un angle de projections limité. Dans l'angle angulaire limité, les projections sont réparties entre [0 • , θ], où θ < 180 • . Dans les travaux précédents, un grand nombre d'articles sont publiés sur la reconstruction de CT de rayons X, en s'appuyant sur différentes méthodes (analytique, algébrique, régularisation, Bayesian, etc.) pour différents types de systèmes de projections (faisceau parallèle, faisceau en éventail, faisceau conique, etc.), compte tenu de différents modèles de bruits (bruit de Poisson, bruit gaussien) et en utilisant différentes méthodes d'estimation (Maximum A Posterior, Posterior Mean, etc). À partir de ces travaux, nous pouvons voir que les méthodes sont plus robustes au caractère incomplet des données lorsque des informations complémentaires sont prises en compte.

Dans la tomodensitométrie 3D, en raison de la grande taille des données, les coûts de calcul sont un frein aux développements de nouvelles méthodes. Au début de la tomographie en rayon X, les méthodes statistiques n'étaient généralement pas utilisées en raison du calcul itératif trop coûteux. Avec le développement du matériel informatique et l'utilisation des aspects de l'accélération du calcul, les méthodes itératives deviennent abordables. Les processeurs many cores de type GPU (Graphic Processing Unit) qui permettent une parallélisation massive des calculs sont à présent largement utilisés dans la reconstruction d'objets 3D.

Un autre défi dans le rayon X est constitué par les artefacts de reconstruction. Il existe de Contributions nombreux types d'artefacts causés par des raisons différentes : les artefacts métalliques, les artefacts en anneaux, les artefacts dû au mouvement, etc. Le bruit d'acquisition est également un défi communément rencontré dans les problèmes de reconstruction de tomographie. Les méthodes itératives sont adaptées pour éliminer ces différents types d'artefacts. Au cours de ma thèse, je ne considérerai que le bruit lié à l'acquisition.

Malgré le fait que les méthodes itératives obtiennent des résultats plus précis et peuvent supprimer les artefacts dans des conditions mal posées avec un ensemble de données insuffisantes, limitées ou tronquées, il y a un autre problème crucial à considérer : la décision d'initialisation des hyper-paramètres. Il existe de nombreux chercheurs qui travaillent sur ce problème. Avec ces techniques, on peut obtenir une valeur optimale pour les paramètres qui conduisent à des bons résultats de reconstruction, mais les coûts de calcul de ces méthodes sont coûteux. Dans la plupart des cas, les personnes choisissent la valeur des paramètres de manière manuelle ou empirique lors de simulations préalables. Dans cette thèse, la méthode bayésienne est utilisée et les paramètres sont estimés simultanément pendant la reconstruction.

L'inférence bayésienne est une fac ¸on couramment utilisée de prendre en compte l'information préalable et de la modélisée par une distribution a priori choisie afin de complexifier le monis possible les calculs. Dans les méthodes bayésiennes, le modèle markovien est fréquemment utilisé pour définir la variable continue par morceaux. La valeur de chaque voxel n'est pas indépendante des autres, mais dépendantes des voxels voisins.

Objet examiné dans la tomographie à rayons X étant composé de plusieurs matériaux homogènes, la propriété continue par morceaux de l'objet est l'une des informations les plus fréquemment considérées dans la reconstruction tomographie de rayons X. Par ailleurs, la parcimonie d'une variable cachée qui est liée à la variable continue par morceaux est également fréquemment considérée. Cette variable cachée qui a une propriété de parcimonie peut être, par exemple, le gradient du volume, le coefficient de transformation d'ondelettes du volume ou un coefficient de la transformée du volume à l'aide d'un dictionnaire. Le mot parcimonie implique que le vecteur ne contient que quelques éléments non nuls. De nombreuses transformations parcimonieuses sont réversibles. On peut reconstruire l'objet avec les coefficients parcimonieux qui contiennent beaucoup moins de valeurs non nulles que l'objet lui-même. Le problème inverse mal posé lié à la reconstruction de l'objet revient alors à l'estimation des coefficients parcimonieux.

Contributions

Un aperc ¸u de la contribution de cette thèse figure à l'annexe B.

Nous proposons dans cette thèse de doctorat principalement deux méthodes dans le contexte bayésien. Ces deux méthodes ont un modèle de système hiérarchique, une variable cachée étant considérée.

En reconstruction tomographique à rayons X, l'objet est typiquement continu par morceaux et, dans certains cas, constant par morceaux. Ainsi, deux aspects du choix de la variable cachée sont considérés dans notre travail : la transformation des ondelettes qui est utilisée dans la méthode HHBM au chapitre 4, et l'operateur Laplacien qui est utilisé dans la méthode ROCC au chapitre 5. Pour ces deux variables cachées, une distribution parcimonieuse leur est associée. Une distribution Student-t généralisée est choisie en raison de sa propriété à queue lourde, CHAPITRE 1. R ÉSUM É EN FRANC ¸AIS modélisant ainsi la propriété de parcimonie, et elle s'exprime comme la loi marginale d'une distribution Gamma Inverse. La structure marginale de deux distributions conjuguées facilite l'estimation des variables dans les méthodes bayésiennes.

Parmi les transformations d'ondelettes, nous choisissons d'utiliser la transformation de Haar discrète. Pour un objet continu par morceaux ou constante par morceaux, la transformation Haar donne un coefficient faible pour la variable. Etant une transformée en ondelettes basique, la transformation de Haar n'est pas cher en termes de calcul et peut être accélérée efficacement à l'aide de processeurs GPU. Le Laplacian de l'objet est considéré comme une variable cachée dans une autre méthode proposée, présentée au chapitre 5. Un modèle markovien non homogène est utilisé pour définir la variable.

Dans cette thèse, le mot "variables" représente les inconnues qui doivent être estimées, par exemple les valeurs des pixels ou des voxels d'objet. Le mot "variables cachées" représente les variables intermédiaires qui sont souvent nécessaires pour modéliser l'objet comme les contours, les étiquettes ou les classifications des régions ou les coefficients de transformation dans les modèles hiérarchiques. Le mot "paramètres" représente les paramètres qui apparaissent dans le modèle a priori des variables telles que la variance du bruit, et ils sont également estimés pendant les itérations. Le mot "hyper paramètres" représente les paramètres des distributions a priori des paramètres, et ils doivent être fixés pour l'initialisation. Dans cette thèse, l'initialisation des hyper paramètres est discutée. Au lieu de tester des différents hyperparamètres et de choisir le meilleur, nous discutons de la fac ¸on de fixer les hyper paramètres directement en condiserant les modèles a priori et les informations préalables sur les variables.

Organisation de th èse

Au chapitre 2, nous présentons le modèle continu et discret du système CT, et le modèle direct de l'opération de projection. Les méthodes conventionnelles, y compris les méthodes analytiques, les méthodes algébriques, les méthodes de régularisation, les méthodes statistiques et l'inférence bayésienne sont présentées. Les avantages et les défauts de chaque type de méthode sont discutés. Ensuite, nous montrons quelques résultats de simulation pour la reconstruction CT en utilisant deux des méthodes les plus représentatives : la technique de rétro-projection filtrée (FBP) et la méthode de Variation Totale (TV). Les résultats de la simulation en utilisant ces deux méthodes avec différents paramètres de projection seront donnés dans ce chapitre. Et je vais montrer l'influence du paramètre dans la méthode de regularisation et démontrer la difficulté de choisir cette valeur.

Dans le chapitre 3, nous présentons le contexte de l'inférence bayésienne. Les propriétés intrinsèques de l'objet (informations a priori) dans les applications industrielles sont discutées. Dans ce chapitre, la propriété de la distribution Student-t généralisée est présentée et l'avantage de l'utiliser est discuté.

Dans la deuxième partie de cette thèse, nous présentons les principales contributions de mon travail.

Au chapitre 4, une méthode HHBM est présentée. Dans cette méthode, tout d'abord, un modèle de synthèse est considéré avec un modèle de bruit non stationnaire. Dans ce modèle, seulement les coefficients de transformation de Haar sont estimés et l'objet est obtenu par un traitement a posteriori à partir des coefficients estimés. Ensuite, dans un modèle plus général, Organisation de th èse l'objet et les coefficients de la transformation de Haar sont simultanément estimés avec un modèle de bruit non stationnaire. Comme le volume se compose de plusieurs blocs homogènes, le coefficient de transformation est faible. Nous avons proposé d'utiliser la distribution Studentt généralisée pour modéliser la répartition préalable du coefficient. Une fois que nous obtenons la distribution a posteriori, un estimateur du Maximum A Posteriori (JMAP) est utilisé pour estimer les variables inconnues. Certaines extensions de cette méthode sont également discutées. Une technique de scission des variables est appliquée sur cette méthode et il est démontré qu'en l'utilisant, on peut considérer le bruit comme le mélange d'un bruit gaussien et d'un bruit parcimonieux. Grâce à cela, les résultats de la reconstruction dépassent la méthode originale. Une autre extension consiste à comparer la transformation de Haar avec la transformation d'ondelettes complexes à double arbre (Dual Tree Complexe Wavelet Transformation DT-CWT). En comparant avec la transformation de Haar, la transformation DT-CWT prend en compte plus de directions d'ondelettes et est plus robuste lorsque les coefficients sont contaminés. Lorsqu'ils sont utilisés dans la méthode bayésienne proposée, le DT-CWT offre une meilleure qualité de reconstruction que la transformation de Haar.

Au chapitre 5, une méthode ROCC est proposée. Dans cette méthode, on considère les relations entre le Laplacien de projection et le Laplacien de l'objet. En utilisant le Laplacien de l'objet comme une variable cachée, nous bénéficions de sa propriété de parcimonie pour préserver les contours de l'objet. Dans cette méthode, les voxels sont modélisés par un modèle markovien non homogène. Dans ce modèle markovien non homogène, chaque voxel est modélisé par une distribution qui dépend des valeurs laplaciennes et des valeurs des voxels voisins. Lorsque son voxel voisin a une valeur Laplacienne plus grande, la valeur voxel voisine aura moins d'influence sur le modèle du voxel actuel.

Les résultats de la reconstruction, correspondant à un fantôme Shepp Logan simulé en 2D ou 3D et un objet Head en 3D, sont présentés au chapitre 4 et 5. Nous avons comparé les résultats obtenues avec nos méthodes et les méthodes de l'état de l'art (FBP, QR et TV). Différents paramètres pour les simulations sont également pris en compte. Le jeu de données avec un nombre insuffisant de projections, avec un angle de projection limité, ou avec un SNR élevé sur l'ensemble de données, sont testés.

Les conclusions et certaines perspectives sont données dans le chapitre 6. Dans des travaux à venir, tout d'abord, d'autres estimateurs peuvent être utilisés, par exemple l'estimateur de la Moyenne a Posteriori (PM). Les estimateurs PM sont généralement réalisés par les approches bayésiennes variationnelles (VBA), avec les détails donnés dans l'annexe A. En utilisant la méthode VBA, les coûts de calcul sont très coûteux pour les simulations de tomographie en 3D, mais notre groupe a récemment trouvé un moyen de simplifier le calcul. Cet algorithme pourrait être mis en oeuvre dans nos travaux à venir. L'analyse de l'influence des hyper-paramètres et des explications de la robustesse des hyper-paramètres dans la méthode HHBM sera poursuivie. D'autres distributions de mélange de variance gaussienne seront également utilisées et comparées avec la distribution de St g . Pour la méthode ROCC, seul le cas de projection parallèle est discuté. Nous chercherons les solutions pour étendre cette méthode à d'autres types de projections.

CHAPITRE 1. R ÉSUM É EN FRANC ¸AIS

La m éthode HHBM

Dans cette méthode, le modèle direct du système est considéré comme suivant :

g = Hf + , (1.2) où g ∈ R M ×1 représente la projection estimé, f ∈ R N ×1 représente l'objet, H ∈ R M ×N
représente l'operateur linear de projection, et ∈ R M ×1 représente le bruit additive. Dans notre travail, on considère un bruit Gaussian dans cette modèle, donc on a :

p( |v ) = N ( |0, V ), where V = diag [v ] , (1.3) où v ∈ R M ×1
représente la variance de bruit, qui est considéré comme un inconnues dans cette méthode semi-supervisé. Considerant la propriété de variance de bruit, il est modélisé par la loi Gamma Inverse :

p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ). (1.4)
Dans cette méthode, une autre relation linear est consideré :

f = Dz + ξ, (1.5) où z ∈ R N ×1 représente le coefficient parcimonieux de f via la transformation représenté par l'opérateur D ∈ R N ×N . ξ ∈ R N ×1 représente les incertitude de cette transformation.
En choisissant le modèle précédent pour ξ, l'objectif est que f et Dz sont le plus approximatifs possible. Par conséquent, nous définissons ξ comme un bruit parcimonieux, avec la plupart des valeurs proches de zéro. La distribution Student-t généralisée est utilisée pour modéliser ξ. La distribution de St g pour ξ est définie par le modèle hiérarchique :

p(ξ|0, v ξ ) = N (ξ|0, V ξ ), where V ξ = diag [v ξ ] , (1.6) p(v ξ |α ξ 0 , β ξ 0 ) = IG(v ξ |α ξ 0 , β ξ 0 ).
(1.7)

À partir de la définition de la loi St g :

St g (x|α, β) = N (x|v)IG(v|α, β) dv, (1.8)
on peut trouver que tous les deux bruit et ξ sont modélisé par la loi St g . Un autre variable modélisé par la loi St g est la variable cachée z :

p(z|v z ) = N (z|0, V z ), where V z = diag [v z ] , (1.9) p(v z |α z 0 , β z 0 ) = IG(v z |α z 0 , β z 0 ).
(1.10) Cette loi est utilisé considerant que la variable z est parcimonieuse.

Avec tous les loi définies, la distribution a postériori peut être obtenu par le théorème de Bayes :

p(f , z, v , v ξ , v z |g) ∝ p(g|f , v )p(f |z, v ξ )p(z|v z )p(v |α 0 , β 0 )p(v ξ |α ξ 0 , β ξ 0 )p(v z |α z 0 , β z 0 ).
(1.11) Puis, l'estimation de Maximum a postériori joindre (JMAP) est utilisé pour estimer tous les inconnus.

La m éthod ROCC

La m éthod ROCC

Dans cette méthode, les Laplacian de projection g et Laplacian d'object f sont consideré, avec les définitions :

f =   0 -1 0 -1 4 -1 0 -1 0   * f = d 2 * f , (1.12) g = [-1 2 -1] * g = d 1 * g, (1.13)
où d 1 et d 2 sont des noyau de convolution, et les convolutions avec eux peuvent représenté par multiplication avec matrices correspondants de la même taille de l'objet, qui sont présenté comme suivant :

g = D 1 g, (1.14) f = D 2 f .
(1.15)

Dans cette method, on considère les relations suivant :

f → RT → g, (1.16) f → RT → g.
(1.17)

En considerant les bruit du système, on a :

g = Hf + , (1.18) g = Hf + τ .
(1.19)

En supposant que tous les deux sont bruit Gaussian, on a la loi de vraisemblance suivant :

p(g|f , v ) = N (g|Hf , V ), V = diag [v ] , (1.20) p(g| f , v τ ) = N (g|H f , V τ ), V τ = diag [v τ ] ,
(1.21) où les variances des bruit sont considéré comme inconnus et il sont modélisés par la loi de Gamma Inverse :

p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ), (1.22) p(v τ |α τ 0 , β τ 0 ) = IG(v τ |α τ 0 , β τ 0 ).
(1.23)

Dans cette méthode, le parcimonie de la variable cachée f est considerée, et il est modélisé par la loi Student-t généralisé via le structure de Gaussian-Gamma Inverse :

p( f |v c ) = N ( f |0, V c ), V c = diag [v c ] , (1.24) p(v c |α c 0 , β c 0 ) = IG(v c |α c 0 , β c 0 ).
(1.25)

Dans cette modèle, la variable f est modélisé par une modèle Markovian, qui depend de la valeur de f estimés. Dans cette modèle, la loi de distribution pour chaque voxel f r est :

p(f r |f r , q, v fr ) = N (f r |m r , v fr ) ∝ v -1 fr exp    - 1 2 v -1 fr f r - r ∈N r(r) (1 -q r ) f r r ∈N r(r) (1 -q r ) 2    , (1.26) CHAPITRE 1. R ÉSUM É EN FRANC ¸AIS
où f -r représente tous les voxels sauf f r , v fr est la variance de voxel r, et :

m r = r ∈N r(r) (1 -q r ) f r r ∈N r(r) (1 -q r )
.

(1.27) La variance v f est considérée comme inconue et être modélisé par la loi de Gamma Inverse :

p(v f |α f 0 , β f 0 ) = IG(v f |α f 0 , β f 0 ).
(1.28)

Dans ce méthode, nous présentons une méthode bayésienne qui utilise un modèle markovien pour la variable continue par morceaux. De la propriété de la transformée de Radon, le Laplacien de l'objet f est lié au Laplacien de la projection sinogram g par un modèle direct. En utilisant simultanément les deux fonctions directes g = Hf + et g = H f + τ , ces deux variables sont estimées. Pendant ce temps, la valeur localisée Laplacienne est utilisée comme coefficient de poids dans le modèle antérieur de f dans chaque itération. Enfin, toutes les variables et les paramètres sont estimés simultanément via la distribution postérieure en utilisant l'algorithme d'estimation JMAP. L'initialisation des hyper paramètres est proposée, et leur relative insensibilité à la reconstruction est prouvée par les simulations. À partir des simulations, nous concluons que la méthode ROCC proposée s'améliore pour préserver les contours en comparaison avec les méthodes conventionnelles de QR et de TV. Lorsque l'ensemble de données est biaisé avec un bruit de SNR=40dB, l'erreur carrée moyenne relative de l'objet reconstruit est plus petite que les méthodes QR et TV, et l'avantage est plus évident lorsqu'il y a moins de projections. La valeur RMSE de la méthode ROCC en utilisant 36 projections est au même niveau que la méthode TV en utilisant 90 projections. Cette méthode est adaptée à la reconstruction CT de projection de faisceau parallèle, où la transformée de radon est utilisée. Le travail futur se concentre sur l'extension du faisceau parallèle aux projections du faisceau conique.

Conclusions et persipectives

Dans cette thèse, nous avons proposé principalement deux méthodes bayésiennes pour la reconstruction 3D CT. Ces méthodes sont basées sur le contexte du Contrôle Non Destructif (CND) dans les applications industrielles. La grande taille des données dans les applications 3D et le coût de calcul sont pris en compte dans les algorithmes.

Dans le chapitre 4, nous avons proposé une méthode bayésienne basée sur le modèle direct du système g = Hf + . Dans cette méthode, la nouveauté est qu'on considère une structure hiérarchique dans le modèle du système, avec f = Dz + ξ. Avec le niveau supplémentaire du système, la variable f est représentée par une variable cachée z. Pour le choix du modèle précédent, la propriété considérée passe de la modélisation de la variable continue par morceaux à la modélisation de la variable cachée qui est parcimonieux. Une distribution Student-t généralisée (St g ) est présentée dans ce chapitre. Grâce au fait que la distribution de St g est marginale d'une distribution Gaussian-Gamma Inverse, les algorithmes d'estimation sont faciles à appliquer.

Dans les algorithmes itératifs, y compris l'inférence bayésienne, l'initialisation des paramètres est cruciale pour les résultats finaux de la reconstruction. Dans la section 4.4.2 et la section 5.6.3, nous avons présenté la stratégie d'initialisation des hyper-paramètres pour l'algorithme d'optimisation dans l'inférence bayésienne avec le modèle proposé. Cette stratégie

Conclusions et persipectives

d'initialisation est démontrée expérimentalement. Dans l'analyse des hyper paramètres, nous avons observé une influence des hyper-paramètres relativement faible sur le comportement de l'algorithme itératif correspondant. L'intérêt de cette dépendance faible relative est qu'il offre un moyen pratique d'initialisation qui n'est généralement pas trivial.

Dans la section 4.5, nous avons présenté une autre transformation orthonormée, la transformation des ondes ondulées à double arbre (DT-CWT). Cette transformation est comparée à la transformation de Haar discrète. Dans la comparaison, nous constatons que la transformation DT-CWT ne provoquera pas les artefacts du bloc dans les figures reconstruites et elle est plus robuste à la contamination des coefficients. Lors de l'utilisation dans la méthode bayésienne proposée, le DT-CWT a une meilleure performance de reconstruction que HT. L'inconvénient est que le calcul du DT-CWT est relativement plus coûteux que l'utilisation de la transformation Haar.

Dans la section 4.8, nous avons proposé d'utiliser un autre modèle avancé dans cette approche bayésienne. Dans ce modèle, le bruit additif est séparé en deux parties, l'un est modélisé par une distribution Gaussien et un autre modélisé par une distribution St g . Le modèle direct considéré pour le système de projection est g = Hf + ρ + .

(1.29)

Dans ce système direct, on considère un modèle de bruit plus compliqué. Les résultats de la simulation montrent que ce modèle, en tenant compte d'un modèle de bruit plus complexe, obtient de meilleurs résultats de reconstruction.

Au chapitre 5, une autre stratégie pour la reconstruction d'objets CT à rayons X est présentée. Dans cette méthode, on considère la relation entre l'objet f , la projection de l'objet g, le Laplacien de l'objet f et le Laplacien de la projection sinogram g. Selon le système de transformation du radon, les relations du système de projection directe que nous prenons en compte sont :

g = Hf + , (1.30) g = H f + τ ,
(1.31) avec quelques relations entre g et g et entre f et f . En considérant ces relations, le Laplacien de l'objet f est considéré comme une variable cachée dans le modèle du système. Un modèle Markovien non homogène est utilisé pour l'objet continu par morceaux f . Dans ce modèle Markovien, la variable cachée f est utilisée comme paramètre. En utilisant ce modèle markovien en considérant la valeur de f , les contours et les zones homogènes sont modélisés avec le même modèle avec des paramètres différents, et donc la continuité des zones homogènes et la discontinuité des contours sont à la fois considérées. Cette stratégie préserve les bords lors de la reconstruction de l'objet.

Pour l'avenir, plusieurs aspects sont pris en considération. Tout d'abord, le temps de calcul devrait prendre en compte. Le processeur GPU n'est utilisé que pour l'opérateur de projection et l'opérateur de rétroprojecteur dans cette thèse. Afin d'accélérer le calcul, plus de calculs dans le programme peuvent être mis en oeuvre en utilisant un processeur GPU. Par exemple, la transformation DT-CWT présentée au chapitre 4 présente de très bonnes propriétés en comparaison avec la transformation de Haar, mais il a un calcul plus compliqué. Une fois que nous implémentons la transformation en utilisant le processeur GPU, l'utilisation du DT-CWT ou des nombreuses autres transformations sera moins coûteuse.

CHAPITRE 1. R ÉSUM É EN FRANC ¸AIS

Dans cette thèse, nous avons utilisé l'estimateur Maximum A Posterior (MAP) pour l'estimation lorsque la distribution postériori est obtenue. Il existe d'autres algorithmes d'estimation, et certains donnent plus d'informations pour les variables. Par exemple, l'estimation de la moyenne postériori peut être utilisée pour estimer les variables. En utilisant la méthode Variation Bayesian Approach (VBA), la distribution a postériori est abordée par une autre distribution séparable :

p(f , θ|g) = q(f , θ) = q 1 (f )q 2 (θ),
(1.32) où θ représente toutes les variables cachées et les paramètres du modèle. La moyenne et la variance de la répartition approximative sont estimées.

Dans la méthode VBA, les distributions approches sont calculées en minimisant la divergence de Kullback Leibler. Au cours du calcul de la méthode VBA, un calcul des éléments diagonaux de H T H de la grande taille est nécessaire, montré à l'Annexe A. Certains résultats de simulation de la méthode VBA sont présentés dans [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF] pour la reconstruction du fantôme 2D Shepp Logan en considérant un bruit stationnaire. En raison de la taille énorme de la matrice H dans le problème 3D, il est très coûteux de calculer directement par projection et rétroprojection pour obtenir chaque élément diagonal. Maintenant, notre groupe est en train de étudier ce sujet et rechercher des algorithmes qui calculent approximativement les éléments diagonaux en utilisant les propriétés de projection et de rétroprojection.

Dans cette thèse, la distribution St g est utilisée pour définir les variables parcimonieuses. Comme nous l'avons mentionné, la distribution de St g a une structure à queue lourde. Dans le même temps, il peut avoir une variance très faible grâce au fait qu'il existe deux paramètres pour contrôler la forme de cette distribution. Une autre propriété très utile dans le contexte bayésien est que, elle peut être exprimée comme un marginal de la distribution de Gaussian-Gamma Inverse bivariée. Comme la distribution normale et la distribution Gamma Inversée sont conjuguées, l'estimation des variables dans les méthodes bayésiennes est plus simple.

Dans le cadre de travaux futurs, au lieu d'utiliser les distributions du mélange Gaussian-Gamma Inverse, on cherche d'autres distributions similaires, par exemple les distributions du mélange de variance normale. On va comparer les différentes distributions de mélange et comparer leurs propriétés, concernant la propriété de la distribution marginale et l'influence des hyper paramètres des distributions marginales. L'un des aspects les plus importants de notre travail futur est d'étudier l'influence de l'initialisation des hyper-paramètres. Dans cette thèse, nous avons observé une influence par semaine relative des hyper paramètres du modèle bayésien selon les modèles antérieurs. Dans le futur travail, nous étudierons la raison de cette faible influence sur les hyper paramètres de ce modèle et voir si la même propriété apparaît lorsqu'on utilise les autres modèles de mélange de variance normale.

Dans un futur travail à long terme, nous chercherons de nouveaux modèles a priori pour la reconstruction de CDT. Nous nous concentrerons sur les cas mal posé où le nombre de projections est extrêmement limité. Pour cela, nous devrions optimiser le modèle a priori et éviter les artefacts apparaissant dans les résultats de la reconstruction. Nous allons également étudier la relation entre les méthodes de régularisation et les méthodes bayésiennes. Il existe des nombreuses méthodes de régularisation qui fonctionnent très bien pour le problème de reconstruction mal posé si un paramètre de régularisation approprié est choisi. Nous pouvons étudier, du point de vue bayésien, le principe de la régularisation. Dans le même temps, nous pouvons bénéficier des avantages de la méthode bayésienne selon laquelle les paramètres sont estimés 2 Introduction

Motivation and Context

This thesis is a contribution to the 3D X-ray Computed Tomography (CT) reconstruction problem. After the discovery of the X radiation by Wilhelm Röntgen, it was quickly used in the domain of tomography. At first, people uses the attenuation projection images to see the object passed by the X radiation. In 1917, the proposition of the Radon transformation [START_REF] Stanley R Deans | The Radon transform and some of its applications[END_REF] solves the problem of the reconstruction. However, X-ray projection in one direction would lose information, for example the depth information. Then the technique of acquisition with the projections applied in different directions is developed in the 1970s, with many publications on both analytic and iterative methods [AGU76, OMA + 76]. In these techniques, the data measured at different angles are used to reconstruct the 3D volume. At the end of 1970s, the fan-beam reconstruction methods were developped. The X-ray CT technique was firstly used in the medical applications soon after its invention [START_REF] Gawler | Computer-assisted tomography (EMI scanner) : Its place in investigation of suspected intracranial tumours[END_REF][START_REF] Rs Ledley | The ACTA-scanner : the whole body computerized transaxial tomograph[END_REF][START_REF] Paxton | The EMI scanner. A brief review of the first 650 patients[END_REF]. Then it began to be used in many other domains, for example the marine science [START_REF] Boespflug | CAT-scan in marine stratigraphy : a quantitative approach[END_REF], the geosciences [WV + 87], the industrial applications [HME + 81], etc. With the development of the reconstruction techniques with full projection data, the interest changed to the study of incomplete data problems.

The foundation of the X-ray CT technique is based on the attenuation of the X rays when passing through an object. A top view of a X-ray projection system is shown in Figure 2.1. Different materials with different density in the object have different attenuation rate for the radiation. This attenuation is defined by the Beer's law :

I = I 0 e -f (l) dl , (2.1)
where l is the monochromatic X-ray passing an object, I 0 is the incident beam intensity, f (l) is the local linear attenuation coefficient along the radiation l and I is the attenuated radiation intensity.

One of the challenges of the X-ray CT is to decrease the number of projections. X-ray CT FIGURE 2.1 -Top view of the X-ray CT projections system.

was firstly used in medical applications since its introduction in the 1970s. It was used in the diagnosis of head, lungs, pulmonary angiogram, cardiac, abdominal and pelvic, extremities, etc [START_REF] Simon R Arridge | Optical tomography in medical imaging[END_REF]. By using the X-ray CT, the diagnosis of the organs can be done with high resolution, in order to distinguish more details. However, this technique has many drawbacks. The most important one is the risks for the patients. The radiation used in CT scans can damage body cells, including DNA molecules, which can lead to cancer, or some other adverse side effect. Therefore recent research efforts focus on techniques for the CT reconstruction with lower dose of radiation, also known as ALARA(As Low As Reasonably Achievable), and less expose time for each projection. In the industrial applications, the X-ray CT is often used in the Non Destructive Testing (NDT) [START_REF] Hanke | X-ray based methods for non-destructive testing and material characterization[END_REF]. It is used to produce 3D representations of components both externally and internally, for the flaw detection, failure analysis, etc. In the industrial applications, the reduction of the radiation dose is not a crucial problem. The main challenge of the industrial application is in the dynamic detection. In some cases the object under examination is not static, for example the object with flowing liquid. In some other industrial examinations, the objects are placed on a conveyor. In these cases, the projection should be done in a short time, hence the number of projections need to be reduced. Hence the reconstructions are expected to be done with less projections.

In both medical imaging and industrial NDT, another problem that we may face to is the limited angle of projections. In typical works, the CT scanner tests a volume by obtaining projecting radiations for angles evenly distributed between [0 • , 180 • ], with full angle coverage for the volume. But in some applications, the full angle coverage could not be obtained. An example in the medical application is the dentistry CT. The projection is sometimes not done for all the angles because of the special structure of the mouth. In industrial applications, this problem normally appears when the component under testing has a large size. When the size of the component is larger than the scope of the scanner, one obtains a limited angle of projections. In the limited angle CT, the projections are distributed between [0 • , θ], where θ < 180 • .

In the previous works, a great number of articles are published on the X-ray CT reconstruction, basing on different methods (analytical [START_REF] La Feldkamp | Practical cone-beam algorithm[END_REF], algebraic [GBH70, AK84, BYL + 00], regularization [CMRS00, GO09, SP08], Bayesian [GCSR14, BMD08, MD96, KVS + 06], etc), for different types of projections systems (parallel-beam [START_REF] Coric | Parallel-beam backprojection : an FPGA implementation optimized for medical imaging[END_REF], fan-beam [START_REF] Reza | Development and validation of MCNP4C-based Monte Carlo simulator for fan-and cone-beam x-ray CT[END_REF], cone-beam [START_REF] La Feldkamp | Practical cone-beam algorithm[END_REF], etc), considering different noise models (Poisson noise [START_REF] Sauer | A local update strategy for iterative reconstruction from projections[END_REF][START_REF] La Feldkamp | Practical cone-beam algorithm[END_REF], Gaussian noise [RLR + 12]), and by using different estimation methods (Maximum A Posterior [QL00], Posterior Mean [START_REF] William A Ericson | A note on the posterior mean of a population mean[END_REF], etc). From these works we can see that the methods are more robust in the ill-posed cases when supplementary prior information are considered.

In 3D X-ray CT, because of the large data size, the computational costs is crucial. At the beginning of the development of X-ray CT, the statistical methods were not typically used because of the very costly iterative computation. With the development of the computational hard-

Contributions

ware and the utilization of the computation acceleration aspects, the iterative methods become feasible in the X-ray CT reconstruction. The GPU processor is now commonly used in the reconstruction of 3D CT objects.

Another challenge in the X-ray CT is the artifacts. There are many types of artifacts caused by different reasons. For example the streak artifacts, ring artifact, motion artifact, etc. The noise of acquisition is also a commonly faced challenge in the X-ray CT reconstruction problems. The iterative methods are adapted to removing these different types of artifacts.

Despite the fact that the iterative methods get more precise results and can reduce the artifacts in ill-posed conditions with insufficient or limited or truncated dataset, there is another crucial problem to consider : the decision for the initialization of hyper-parameters. There are many researchers working on this problem [FSH15, MDA15, SC91, ZLQ97, SBS98]. With these techniques, one can obtain an optimal value for the parameters which lead to good reconstruction results, but the computational costs of these methods are expensive. In most cases, people choose the value of the parameters manually, or empirically with some beforehand simulations. In this thesis, the Bayesian inference is considered and the parameters and unknown variables are estimated simultaneously during the reconstruction.

Bayesian inference is a commonly used way to take into consideration the prior information. In the Bayesian context the prior information is modeled by a prior distribution. The convenience by using the Bayesian inference is that the prior properties could be taken into consideration with a variaty of choice of prior distributions, and among the feasible prior distributions we could choose those who simplifies the computation. In Bayesian inference, the Markovian model is frequently used to define the piecewise continuous variable. The value of each voxel is not independent from the others, but depending on the neighbor voxels.

Thanks to the fact that the object under examination in X-ray CT is composed of several homogeneous materials, the piecewise continuous property of the object is one of the most commonly considered prior information in the X-ray CT reconstruction. Besides, the sparsity of a hidden variable which is related with the piecewise continuous variable is also considered frequently. This hidden variable which has a sparsity property can be, for example, the gradient of the volume, the wavelet transformation coefficient of the volume or a dictionary basis coefficient of the volume. The word "sparsity" implies that the vector contains only few non zero entries. Many sparse transformations are reversible. One can reconstruct the object with the sparse coefficient which contains much less non zero values than the object itself. The ill-posed inverse problem of reconstructing the object can then be transferred to the reconstruction of the sparse coefficients.

Contributions

An outline of the contribution of this thesis is given in Appendix B.

We propose in this PhD thesis mainly two methods in the Bayesian context. Both of these two methods have a hierarchical system model, with a sparse hidden variable being considered. A sparsity enforcing prior distribution is used to model the hidden variable.

In the X-ray CT reconstruction, the object is typically piecewise continuous, and in some cases piecewise constant. So, two aspects for choosing the sparse hidden variable are considered in our work : the wavelet transformation which is used in the HHBM method, Chapter 5,

CHAPITRE 2. INTRODUCTION

and the Laplacian which is used in the ROCC method, Chapter 6. Both of them are sparse, and considered as a hidden variable, associated with a prior distribution which enforces the sparse structure of the hidden variable. A generalized Student-t distribution [START_REF] Dumitru | A Bayesian approach for periodic components estimation for chronobiological signals[END_REF][START_REF] Dumitru | Performance comparison of Bayesian iterative algorithms for three classes of sparsity enforcing priors with application in computed tomography[END_REF] is chosen because of its heavy-tailed property, modeling the sparsity property, and it is expressed as the marginal of a Normal-Inverse Gamma distribution. The marginal structure of two conjugate distributions facilitate the estimation of the variables in the Bayesian methods.

Among the wavelet transformations, we choose to use the discrete Haar transformation. For a piecewise constant or piecewise continuous object, the Haar transformation gives a sparse coefficient for the variable. As one of the basic wavelet transformation, the Haar transformation is computationally cheap, and can be efficiently accelerated using the GPU processor. The gradient of the object is considered as a hidden variable in another proposed method, presented in Chapter 6, and a non-homogeneous Markovian model is used to define the variable.

In this thesis, the word "variables" represents the unknowns which need to be estimated, for example the object's pixels or voxels values. The word "hidden variables" represents the intermediate variables which are often needed to model the object such as contours, region or classification labels or the transformation coefficients in the hierarchical models. The word "parameters" represents the parameters which appear in the prior model of the variables such as the variance of the noise, and they are also estimated during the iterations. The word "hyperparameters" are the parameters of the prior distributions of the parameters, and they need to be fixed for the initialization. In this thesis, the initialization of the hyper-parameters is discussed. Instead of testing for different hyper-parameters and choosing the best one, we discuss the way to fix the hyper-parameters directly by taking into consideration the prior models and the prior information on the variables.

Organization of thesis

In Chapter 3, we present the continuous and discrete model of the CT system, and the direct model of the projection operation. The conventional methods, including the analytical methods, the algebraic methods, the regularization methods, the statistical methods, and the Bayesian inference are presented. The advantages and shortcomings for each different type of methods are discussed. Then we show some simulation results for the CT reconstruction by using two of the most representative methods : the Filtered Back Projection (FBP) technique and the Total Variation (TV) method. The simulation results by using these two methods under different settings of the projection system will be given in this chapter. And I am going to show the influence of the parameter in TV and demonstrate the difficulty to choose this parameter value.

In Chapter 4 we present the context of the Bayesian inference. The prior properties of the object in the industrial applications are discussed. In this chapter, the property of the generalized Student-t distribution is presented and the advantage of using it is discussed.

In the second part of this thesis, we present our main contributions.

In Chapter 5, a HHBM method is presented. In this method, first of all a synthesis model is considered with a non-stationary noise model. In this model, only the Haar transformation coefficients are estimated, and the object is obtained by a post-processing from the estimated coefficients. Then, in a more general model, the object and the Haar transformation coefficients Organization of thesis are simultaneously estimated with a non-stationary noise model. As the volume consists of several homogeneous blocks, the transform coefficient is sparse. We proposed to use the generalized Student-t distribution to model the prior distribution for the coefficients. Once we get the posterior distribution, a Joint Maximum A Posterior estimator (JMAP) is used to estimate the unknown variables. Some extensions of this method are also discussed. A variable splitting technique is applied on this method and it is shown that by using it, we can consider the noise as the mixture of a Gaussian noise and a sparse noise. Thanks to that, reconstruction results surpasses the original method. Another extension is to compare the Haar transformation with the Dual-Tree Complex Wavelet Transform (DT-CWT). Comparing with the Haar transformation, the DT-CWT transformation takes into account more wavelet directions, and is more robust when the coefficient is contaminated. When they are used in the proposed Bayesian method, the DT-CWT offers a better quality of reconstruction than the Haar transformation.

In Chapter 6, a ROCC method is proposed. In this method, the relations between the Laplacian of projection sinogram and the Laplacian of the object is considered. By using the Laplacian of the object as a hidden variable, we benefit from its sparse property to preserve the contours of the object. In this method the voxels are modeled by a non-homogeneous Markovian model. In this non-homogeneous Markovian model, each voxel is modeled by a distribution which depends on the laplacian values and the voxel values of its neighbor voxels. When its neighbor voxel has a bigger laplacian value, the neighbor voxel value will have less influence on the model of the current voxel.

The reconstruction results, corresponding to a simulated 2D or 3D Shepp Logan phantom, and a 3D Head object, are shown in Chapter 5 and Chapter 6. We compared the proposed methods with the FBP, QR and TV methods. Different settings for the simulations are also considered. Observations with insufficient number of projections, or limited angle of projections, with both high and low SNR, are tested.

The conclusions and some perspectives are given in Chapter 7. In the future work, firstly, other estimators can be use, for example the Posterior Mean (PM) estimator. PM are usually realized by the Variational Bayesian Approach (VBA) method, with details given in Appendix A. By using the VBA method, the computational costs is very expensive for the 3D CT simulations, but our group have recently find a way to simplify the computation. This algorithm could be implemented in our future work. The analysis of the influence of hyper-parameters and the explanations of the robustness of the hyper-parameters in the HHBM method will be continued. Some other Gaussian variance mixture distributions [START_REF] Dumitru | Unsupervised sparsity enforcing iterative algorithms for 3D image reconstruction in X-ray computed tomography[END_REF] will also be used and compare with the St g distribution. For the ROCC method, only the parallel projection case is discussed. We will search for the solutions to extend this method to other projection geometries.

3

X-ray CT and the state-of-the-art methods

Introduction

The 3D CT reconstruction is treated as a typical inverse problem but with a huge data size. With the increase of the computing speed of computer hardware, the research interest changed from the analytical methods to the iterative methods. The iterative methods are now more and more frequently used in real applications. This chapter is devoted to the presentation of the reconstruction methods for X-ray CT reconstruction. The Filtered Back-Projection (FBP) based on the Radon transform is still the most commonly used reconstruction algorithm in real applications. Research has been focused on the reconstruction with fewer projections, which makes it an ill-posed inverse problem. In this case, the iterative methods are used to optimize the reconstruction results. Algebraic Reconstruction Technique (ART) [BKK12, BKK12, CM15, GBH70], Simultaneous algebraic reconstruction technique (SART) [START_REF] Anders | Simultaneous algebraic reconstruction technique (SART) : a superior implementation of the ART algorithm[END_REF][START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF] and Simultaneous Iterative Reconstruction Technique (SIRT) [BYL + 00, WLL14] are a few of those firstly proposed iterative methods. Generally we get an under-determined system of equations for the reconstruction of X-ray CT. In this context, the regularization methods are used such that complementary prior information is combined in the criterion. The statistical method are frequently used as the statistical models can be used to define the unknowns and be estimated by using the statistical estimator. The Maximum Likelihood (ML) methods are one type of the statistical methods. Together with their different algorithms such as the Expectation Maximization (EM) algorithm [START_REF] Todd K Moon | The expectation-maximization algorithm[END_REF], the stochastic EM (SEM) [TET + 04] and the Ordered subsets-EM (OS-EM) [START_REF] Malcolm | Accelerated image reconstruction using ordered subsets of projection data[END_REF], the ML methods are commonly used in PET-CT reconstruction problems. Another widely used type of statistical methods for PET or X-ray CT reconstruction are the Bayesian methods [AMD10, BMD08, FR12, KVS + 06, ZHL + 16], and they are attracting more and more attentions in the field of tomography as the computation time is significant reduced by using the GPU processor.

In this chapter, we are going to present the basic analytical methods, the regularization CHAPITRE 3. X-RAY CT AND THE STATE-OF-THE-ART METHODS methods and the Bayesian methods for the X-ray CT reconstruction.

The X-ray Computed Tomography

The technique of the X-ray CT has been mainly used for medical imaging and for Non Destructive Testing (NDT).

The Radon Transform and Analytical reconstruction techniques

The Radon transform is based on the work of the mathematician Johann Radon in 1917 [START_REF] Radon | Uber die bestimmug von funktionen durch ihre integralwerte laengs geweisser mannigfaltigkeiten[END_REF]. The problem of reconstructing a function f (x, y) (in 2D case) or f (x, y, z) (in 3D case) from the integrals along different directions is becoming widely used in many domains, for example the seismology, the medical imaging, the astronomy, the industrial Non Destructive testing, etc. The Radon transformation is expressed by operator R, and g is the corresponding function In tomography, the function f (x, y, z) represents an unknown density, and the integrals of the function is called the projection, or sinogram.

From the Beer's law, Eq.(2.1), by applying the logarithm in both size, we obtain : (3.4)

log I = log I 0 -f (l) dl, ( 3
And we obtain :

g(r, φ) = +∞ -∞ +∞ -∞ f (x, y)δ (r -x cos φ -y sin φ) dx dy. (3.5)
Mathematically, the adjoint Radon transform, also known as the dual Radon transform, is expressed as : where * * is the 2D convolution and h(x, y) = 1

f (x, y) = R [g(r, φ)] = 1 2π π 0 +∞ -∞ g(
x 2 + y 2 . (3.8)
Geometrically, it is simply a propagation of the measured sinogram back into the image space along the projection paths. In the following sections, I'm going to describe in a short and synthetic way the analytical relations between f (x, y) and g(r, φ) which are the basis of many analytical reconstruction methods.

The central slice theorem

The central slice theorem indicates that : the 1D Fourier Transformation (FT) of a projection line

G(ρ, φ) = F 1D(r) {g(r, φ)} is equal to the 2D FT of F (u, v) = F 2D {f (x, y)} evaluated at angle φ : G(ρ, φ) = F (u, v)| u=ρ cos φ,v=ρ sin φ = F (ρ, φ).
(3.9)

The relations between the figure 2D f (x, y), the projection sinogram g(r, φ) and their Fourier Transformations F (u, v), G(ρ, φ) are illustrated in Figure 3.5.

This property leads to a direct reconstruction strategy : the direct Fourier interpolation method. By using this method the object is reconstructed by the following steps :

The Radon Transform and Analytical reconstruction techniques f(x,y) g(r, )

F(u,v) G( , ) Radon R FT FT u=x cos£ , v=y sin£ 2D 1D F 2D F 1D FIGURE 3
.5 -The central slice theorem.

1. Apply 1D FT for each projection. 2. Convert from polar coordinates F (ρ, φ) to Cartesian coordinates F (u, v) using u = ρ cos φ and v = ρ sin φ. 3. Apply inverse 2D FT to F (u, v) and obtain f (x, y).

The Back-projection filtering (BPF) algorithm

As presented previously, the back-projection is the adjoint operator of the Radon transform :

b(x, y) = 1 2π π 0 +∞ -∞ g(r, φ)δ(x cos φ + y sin φ -r) dr dφ. (3.10)
Consider that :

g(r, φ) = F -1 1D(ρ) {F (ρ, φ)} = +∞ -∞ F (ρ, φ) exp {i2πρr} dρ, (3.11)
we have :

b(x, y) = 1 2π +∞ -∞ π 0 +∞ -∞ F (ρ, φ)δ(x cos φ + y sin φ -r) exp {i2πρr} dr dφ dρ = 1 2π π 0 +∞ -∞ F (ρ, φ) exp {i2πρ(x cos φ + y sin φ)} dφ dρ.
(3.12)

From the definition of FT we have F (-ρ, φ) = F (ρ, φ + π), so we obtain :

b(x, y) = 1 2π 2π 0 +∞ 0 F (ρ, φ) ρ exp {i2πρ(x cos φ + y sin φ)} ρ dφ dρ = 1 2π F -1 2D F (ρ, φ) ρ .
(3.13)

From the obtained relation we get :

f BP F (x, y) = F -1 2D ρF 2D f (x, y)
, and the Backprojection filtering method implies the following steps :

1. Back-projection from g(r, φ) to get b(x, y).

2. Apply 2D FT for b(x, y) and we get F (ρ,φ)

ρ = F 2D {b(x, y)}.
3. Filter with ρ and obtain F (ρ, φ).

Apply inverse 2D FT to get f BP F (x, y).

The BPF algorithm is illustrated in Figure 3.6.

-

The Filtered Back-projection

The Filtered Back-Projection (FBP), illustrated in Figure 3.7, is a conventional analytical technique for the reconstruction, and it is the most commonly used technique in real medical and industrial applications. In FBP, the projection data is convoluted by a filter function in the Fourier space, and then a back-projection is applied. Mathematically it is presented by

f F BP = BF -1
1D |ρ| F 1D g where B represents the back-projection. Comparing with back-projection, FBP compensates the high frequency component of the phantom, and derives results with much clearer contours than the back-projection.

f(x,y) g(r, ) G( , ) G F ( , ) g F (r, ) f FBP (x,y) Radon FT Filter with |¢| FT -1
Backprojection FIGURE 3.7 -The Filtered Back-projection reconstruction algorithm.

The algorithm of the FBP method implies the following two steps :

1. Filter the projection with |ρ| in the Fourier domain or by a convolution method.

2. Apply the back-projection.

With these two steps we have :

f F BP (x, y) = π 0 +∞ -∞ F -1 1D {|ρ| F 1D {g(r, φ)}} δ (x cos φ + y sin φ -r) dr dφ = π 0 +∞ -∞ F -1 1D {|ρ| F (ρ, φ)} δ (x cos φ + y sin φ -r) dr dφ = +∞ -∞ π 0 +∞ -∞ |ρ| F (ρ, φ) exp {i2πρr} δ (x cos φ + y sin φ -r) dρ dr dφ = π 0 +∞ -∞ F (ρ, φ) exp {i2πρ (x cos φ + y sin φ)} |ρ| dρ dφ = 2π 0 +∞ 0 F (ρ, φ) exp {i2πρ (x cos φ + y sin φ)} ρ dρ dφ = F -1 2D {F (ρ, φ)} = f (x, y). (3.14)
In the real applications, a low-pass filter, for example the Hanning filter, is used during the filter for the projection in order to reduce the noise.

The Radon Transform and Analytical reconstruction techniques

However, the presented analytical reconstruction methods, for example the most commonly used FBP method, have a good performance only when there is sufficient number of projections and that the projections are distributed from 0 to π. We can see from the results that, when there are 64 projections distributed from 0 to π, the reconstructed figure is clear and we can distinguish the details. When there are only 32 projections, distributed from 0 to π, some obvious artifacts appears in the reconstructed figure. In the case where the angle of projection is limited, the reconstruction has a poor quality. In Figure 3.9, the reconstruction of the Shepp Logan figure by using the FBP algorithm is presented. As a conventional analytical method, we can see that when the number of projections is sufficient and the projections are obtained from 0 to π, the reconstruction result is clear. However when the number of projections or the angle of projections are insufficient, the reconstruction results by using the FBP method are blurred and some obvious artifacts appear in the results. Another shortcoming of the FBP method is its sensitiveness to the noise, because of that the high frequency parameters corresponding to the noise are not filtered in the projection domain.

The algebraic methods

In imaging problems, the function f (x, y) is typically considered in a discrete form, and is presented by a vector f . The projection is presented by vector g. In a discrete system, the integral operation becomes the sum. Figure 3.10 illustrates the Radon transform of a discretized 2D image. After discretization, the image is represented by pixels. The value of each pixel is the average of the continuous function f (x, y) in the pixel area.

x y r ¡ f j H ij projection i f 1 f N g i projection i
Detector FIGURE 3.10 -The Radon transform of discrete system.

From the figure, we can see that the projection line passes the pixels, with different length of segment in different pixels. We use H ij to represent the length of the ith projection line in

The algebraic methods the pixel j. According to the definition, we have :

g i = N j=1 H ij f j , (3.15)
where j ∈ [1, N ] and N is the size of the image f . It can also be expressed as :

g i = [H i1 H i2 • • • H iN ]      f 1 f 2 . . . f N      . (3.16)
We use the vector g ∈ R M ×1 to represent the dataset of the projection, g = [g 1 , g 2 , • • • , g M ] T , and we can obtain :

     g 1 g 2 . . . g M      =      H 11 H 12 • • • H 1N H 21 H 22 • • • H 2N . . . . . . . . . . . . H M 1 H M 2 • • • H M N           f 1 f 2 . . . f N      , (3.17)
and it can be expressed as :

g = Hf . (3.18)
This expression is used as the forward model of the projection system. Each row of matrix H corresponds to a projection line, considering every pixel of the image. Each column of matrix H corresponds to all the projection lines in all angles for one pixel f j .

Consequently, the transpose of the matrix H is used as the back-projection operator :

f BP = H T g. (3.19)
Based on the discrete forward model of the projection system in Eq.(3.18), the algebraic methods are used, such as the algebraic reconstruction technique (ART) [START_REF] Gordon | A tutorial on ART (algebraic reconstruction techniques)[END_REF][START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[END_REF], the related simultaneous ART (SART) [START_REF] Anders | Simultaneous algebraic reconstruction technique (SART) : a superior implementation of the ART algorithm[END_REF] and the Simultaneous Iterative Reconstruction Technique (SIRT) [TL90, BYL + 00]. The update step for a voxel f j in ART, SART and SIRT are respectively expressed as :

ART : f (k+1) j = f (k) j + λ g i -N n=1 H in f (k) n N n=1 H 2 in H ij , (3.20) SART : f (k+1) j = f (k) j + λ g i ∈G φ g i -N n=1 H in f (k) n N n=1 H in H ij g i ∈G φ H ij , (3.21) SIRT : f (k+1) j = f (k) j + λ g i ∈G g i -N n=1 H in f (k) n N n=1 H in H ij g i ∈G H ij , (3.22)
where λ is a relaxation parameter and H ij is the weight with which the voxel j contributes its value to the projection i. In ART, only one ray in one projection direction is considered for one iteration of voxel j. In SART, all rays in one projection G φ in direction φ is considered for one

-

The Filtered Back-projection iteration. In SIRT, all rays in all projections are considered for one iteration. A Discrete Algebraic Reconstruction Technique (DART) is presented in [START_REF] Kees | DART : a practical reconstruction algorithm for discrete tomography[END_REF]. The DART technique considers the reconstruction of the strictly piecewise constant objects, which consists of several different homogeneous materials, and the gray level of each different material is known. For this special types of problems, the DART method has a very good performance with insufficient number of projections and is robust to the noise. However, in our work, we consider a more general case, where we have no information on the number of materials and the gray level of each material.

In Figure 3.11, we show the reconstructed Shepp Logan figure of size 256 2 by using the SIRT method. From the results we can see that, when there are sufficient number of projections, the shape of the reconstructed phantom is clear, but the contours in the phantom is a little blurring. When there is insufficient number of projections, there will be many artifacts in the reconstructed figure. From the definition and the simulation results, we can see that the limitation of the algebraic reconstruction techniques is the very costly computation for the calculation of the weighted sum. What's more, the size of detected data g will have a big influence on the reconstruction result. In the ill-posed reconstruction problems, where the number of observed data is much less than the unknown data, the prior information need to be considered during the reconstructions.

In the following section, I am going to present the regularization methods, which are fundamental and widely used methods concerning the prior information.

The Regularization methods

In X-ray CT, generally the object under testing needs to be reconstructed with precise details and the dimension of voxel is very small compared with the size of object. With a large size of unknown variable, the acquisition of sufficient dataset is too expensive. In most of the cases we get an under-determined system of equations (N > M ). The following direct model is used to account for the modeling error and measurement noise :

g = Hf + , (3.23)
The Regularization methods

where ∈ R M ×1 corresponds to all the additive noise and uncertainties of the projection system.

The regularization methods are derived from this model. In regularization methods, we optimize a criterion which consists of generally two terms : the data consistency Q(g, f ) and the regularization term R(f ) :

J(f ) = Q(g, f ) + λR(f ), (3.24)
where λ is called the regularization parameter.

The data consistency term Q(g, f ) describes the discrepancy between the observation and estimation. Based on direct model defined in Eq.(3.24), it is commonly under form of Q(|g -Hf |). The l 2 norm is mostly used in the regularization methods, with :

Q(g, f ) = g -Hf 2 2 . (3.25)
The L q norm is a general form, with Q(g, f ) = g -Hf q q , for example q = 1. The diversification of the criterion of regularization methods is normally on the regularization term R(f ). This term is fixed by considering the prior information that we known about the unknown variable f . We can see that when R(f ) = 0, we obtain a Least Square (LS) method [YAT + 97, SC00].

Different regularization criterion

The regularization term R(f ) is typically a penalty on the complexity criterion of f . By choosing different regularization function R(f ), different regularization methods are referred to.

When R(f ) = Φ(f ) 2 2 , the ridge penalty, or the Quadratic Regularization (QR) [RSD09] is derived, where the linear operator Φ(•) can be the multiplication by an identical matrix which lead to R(f ) = f 2 2 . Φ(•) can also be the gradient operator. The L 2 form of the regularization term will result in the smoothness in the final results.

The Total Variation (TV) method [CL97, GO09, PTS + 15] is defined by

R(f ) = D T V f T V where D T V is the gradient operator. In the anisotropic form, R(f ) = D x f 1 + D y f 1 for 2D image and R(f ) = D x f 1 + D y f 1 + D z f 1 for 3D object,
where D x , D y and D z are respectively the gradient in the x, y and z direction of the Cartesian coordinate system. L 1 norm is used for sparse estimations, which enforces the sparsity of the D T V f term.

Sometimes there are more than one regularization terms in the criterion. For example in the Elastic-net regularization method [START_REF] Winters | A sparsity regularization approach to the electromagnetic inverse scattering problem[END_REF], the regularization term is

R(f ) = λ 1 f 1 + λ 2 f 2 2 .
The elastic net penalty is a convex combination of the lasso and ridge penalty. When λ 1 = 0, it becomes a simple ridge regression.

More general regularization methods are developed basing on constrained and dual-variable regularization method :

J(f , z) = Q 1 (g, f ) + ηQ 2 (f , z) + λR(z).
(3.26)

In such a model, the penalty regularization term is set on z, which is associated with f via a linear transformation. The loss function

Q 1 (g, f ) and Q 2 (f , z) are generally quadratic, i.e. Q 1 (g, f ) = g -Hf 2 2 and Q(f , z) = f -Dz 2 2 .
They describe respectively the forward model consistency and the prior model consistency. R(z) is the regularization term on z.

-Optimization algorithms

Optimization algorithms

In the problem of optimizing large data size variable, the optimization result sometimes contains the inversion of big size matrix, which is impossible to compute. In these cases, the optimization algorithms are often used. The optimization of the variables are derived by minimizing the loss function J(f ) with respect to the unknowns.

f = arg min f {J(f )} . (3.27)
The optimization algorithms fall in two major categories : the first order optimization algorithms and the second order optimization algorithms.

The first order optimization algorithms minimize or maximize the loss functions by considering their gradient. The Gradient Descent (GD) algorithm [START_REF] Battiti | First-and second-order methods for learning : between steepest descent and Newton's method[END_REF][START_REF] Ma Vorontsov | Adaptive phase-distortion correction based on parallel gradient-descent optimization[END_REF], which is the most important technique and the foundation of the optimization algorithms, belongs to the first order optimization category. The GD algorithm takes steps proportional to the negative of the gradient of the loss function at the current point and iteratively reaches the local minimum or maximum of the loss function :

f (k+1) = f (k) -γ∇J(f ), (3.28)
where ∇J(f ) is the derivative of loss function and γ is called the descent step length.

Though the most widely used algorithm, the GD algorithm has some drawbacks, for example the convergence is slow and the final result might be the local minimum. Many other first order algorithms are developed to overcome these problems, for example the conjugate gradient algorithm [Bat92, Sca87, FB99], the stochastic gradient descent algorithm [Bot10, NWS14, KPSV09], etc.

The second order optimization algorithms are less frequently used, as they use the second order derivative, also known as the Hessian, to minimize or maximize the loss function, and is more costly than the first order derivative. However, in some cases the form of the second order derivative is known, they can be used for the big data size optimization problems. The most famous second order optimization algorithm is the Newton's method [START_REF] Carl | Iterative methods for optimization[END_REF][START_REF] Fischer | A special Newton-type optimization method[END_REF].

In gradient based optimization algorithms, the derivative of the criterion function always need to be calculated. Note that the L 1 norm is not derivable everywhere, the conventional optimization algorithms can not be used. Many methods have been proposed to solve this L 1 norm optimization problem, for example the Primal Dual method [CGM99], the Split Bregman method [GO09] and so on.

Other methods are also used for the large scale regularization methods, for example the augmented Lagragian method [TW09, ABDF11], the Alternating Direction Method of Multipliers (ADMM) [START_REF] Wahlberg | An ADMM algorithm for a class of total variation regularized estimation problems[END_REF], etc. These methods are often applied to solve the constrained optimization problems. The ADMM method considers to minimize Φ(f ) + Ψ(z) subject to Af + Bz = C, and it covers a large amount of estimation forms. An example is when

Φ(f ) = g -Hf 2
2 , Ψ(z) = R(z), A = I, B = -D and C = 0, referring to the above mentioned bi-variable regularization method corresponding to Eq.(3.26).

The Regularization methods

The regularization parameter

In the above mentioned regularization methods, there is always a regularization parameter λ to be fixed. This parameter compromise between not having enough regularization to reduce the noise and too much regularization that the solution becomes too smooth. Sometimes there are more than one regularization terms in the minimization criterion, for example the Elastic-net method. Hence two or even more regularization parameters are introduced. These parameters tune the terms of regularization. The value of these parameters will effect on the optimization result.

Figure 3.12 and Figure 3.13 demonstrate respectively the influence of the regularization parameter λ in the QR and TV method for the reconstruction of the Shepp Logan phantom. From the figures we can see that the optimal regularization parameter value is different for the cases with different number of projections or different SNR. In order to fix the regularization parameter, many method have been proposed. A most intuitive strategy is to fix them manually. In [KVS + 06], the authors proposed to fix the parameters with training dataset with similar distribution structure as the unknown variable. They demonstrated that a fixed parameter usually performs well for similar imaging data. The reconstruction program should be done many times for the training object in order to find the suitable parameter. What's more, the parameter is not ensured to be suitable for the same structured data with different noise level.

The Generalized Cross Validation (GCV) [GHW79, RLR + 12] and L-curve [START_REF] Calvetti | Tikhonov regularization and the L-curve for large discrete ill-posed problems[END_REF][START_REF] Per | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF] methods are also proposed to determine suitable values for this parameter. The prin- cipal strategy of the L-curve method is to set the data consistency term and the regularization term of the criterion as the abscissa and the ordinate respectively, and by setting different value of the parameter λ, we can get a curve which has a 'L' shape. The point which has the biggest curvature corresponds to the suitable parameter.

-The regularization parameter

In the GCV method, a GCV function is defined depending on parameter λ. For example for the Tikhonov regularization J(f ) = g -Hf 2 2 + λ Df 2 2 , the GCV function is :

G(λ) = H H T H + λD T D -1 -I g 2 2 T r I -H H T H + λD T D -1 H T 2 . (3.29)
Then the optimal parameter value is estimated by minimizing the GCV function. However, in many cases the GCV function does not have an explicit minimum point and the parameter obtained by the minimization is not accurate. For the reconstruction problem of large data size object, this method is too costly.

The statistical methods

When data are corrupted by noise, statistical methods are often used. With the development of Computer Science, the statistical inference [START_REF] Wang | Statistical image reconstruction for muon tomography using a Gaussian scale mixture model[END_REF] are used in X ray tomography reconstruction. For the inverse problems which consider a large amount of data, the statistical The statistical methods methods set a probability model. Instead of estimating all the parameters, the parametric model is examined.

In the forward model presented in Eq.(3.23), the noise is modeled by a probabilistic density function (pdf). In tomography, the acquisition of photon count is obtained by detecting the photons arrived at the detectors, which is modeled by a Poisson distribution :

p(g|f ) = M i=1 [Hf ] g i i g i ! e -[H f ] i . (3.30)
where g i is the detected photon rate, as shown in Figure 3.14.

Emitter Detector g T g i FIGURE 3.14 -One X ray projection from emitter to detector.

In X-ray CT, the amount of photon is very big. With the large amount of photon, the Poisson distribution can be approximate by a Gaussian distribution [START_REF] La Feldkamp | Practical cone-beam algorithm[END_REF].

The likelihood plays an important role in the statistical inference. It describes the probability density function of parameters of the model knowing data. For a forward model shown in Eq.(3.23), the likelihood is p(g|f ). In X-ray CT, it is modeled by a Gaussian distribution :

p(g|f ) = N (g|Hf , V ), (3.31)
where

V = diag [v ], v = [v 1 , • • • , v i , • • • , v M ]
and v i is the variance of noise of projection i.

When v i = v , ∀i, it refers to a stationary noise. However for the X-ray CT applications, the noise depends on the projection system and also the object. The shape and the material structure of the object will have an influence of the detection noise. Hence a non-stationary noise is considered. By using the Maximum Likelihood (ML) methods [START_REF] Richard | Mixture densities, maximum likelihood and the EM algorithm[END_REF], the statistical parameters are estimated by maximizing the likelihood of making the observations given the parameters :

f = arg max f {p(g|f )} . (3.32)
In the Emission Computed Tomography, the detection of the photons is modeled by the Poisson distribution, and so the logarithm of the likelihood is :

ln p(g|f ) = M i=1 - N n=1 H in f n + g i ln N n=1 H in f n -ln g i ! .
(3.33)

-The regularization parameter

As demonstrated in [START_REF] Vardi | A statistical model for positron emission tomography[END_REF], the EM algorithm of this ML method can be written as :

f (k+1) j = f (k) j    1 + i∈I j g i -N n=1 H in f (k) n N n=1 H in f (k) n H ij i∈I j H ij    = f (k) j + f (k) j i∈I j g i -N n=1 H in f (k) n N n=1 H in f (k) n H ij i∈I j H ij , (3.34)
where I j is the set of projections to which the pixel j contributes. We can see that this result has a very similar form as the SIRT method shown in Eq.(3.22).

As mentioned before, the Poisson noise is approximated by a Gaussian noise in X-ray CT. By using the Gaussian model, the maximization of the likelihood is equivalent to the weighted Least Square regression method :

f = arg max f {N (g|f , V )} = arg max f |V | -1 2 exp - 1 2 (g -Hf ) V -1 (g -Hf ) = arg min f (g -Hf ) V -1 (g -Hf ) = arg min f g -Hf V -1 , (3.35)
where V -1 is the weight of the Least Square method.

The Expectation Maximization (EM) [START_REF] Todd K Moon | The expectation-maximization algorithm[END_REF] algorithms are used to solve the Maximum Likelihood estimations when the equations are difficult to solve directly. Typically the models depends on some other unobserved hidden variables. The EM iteration alternates between performing an expectation (E) step, and a maximization (M) step. Some of its extensions, for example the stochastic EM (SEM) [TET + 04] or the Ordered subsets-EM (OS-EM) [START_REF] Malcolm | Accelerated image reconstruction using ordered subsets of projection data[END_REF] are also used for the tomography reconstruction problems.

The shortcomings of the Maximum Likelihood methods are that, the numerical estimation is usually non-trivial, except for a few cases where the maximum likelihood formulas are simple, and that they do not consider any prior information. Consequently, they are not robust in the illposed reconstruction systems where the observations are much less the unknowns. Therefore, we present the Bayesian methods, which combine the model and the prior information together to estimate the unknowns.

The Bayesian inference

The Bayesian inference [SB92, GRL + 93, MD96, BMD08, GCSR14, ZFR15] is a type of statistical methods, which is widely used in the inverse problems. In Bayesian inference, Bayes' theorem is used to update the probability for a hypothesis as more evidence or information becomes available. The Bayes' rule derives the posterior distribution as a consequence of the The Bayesian inference prior law and the likelihood, with the basic expression :

p(f |g) = p(g|f )p(f ) p(g) , (3.36) with p(g) = p(g|f )p(f ) df , (3.37)
where p(g|f ) is the likelihood, p(f ) is the prior law and p(f |g) is the posterior law. p(g) is called the marginal likelihood or the evidence. As with the marginalization the variable f doesn't appear in this term, so it is normally not considered in the determination of the variables. Figure 3.15 shows the generative graph of the basic Bayesian model. In this model, the parameters θ 1 and θ 2 are fixed and are considered as known constant in the model.

θ 2 θ 1 f g p(g|f , θ 1 ) p(f |θ 2 ) FIGURE 3.15 -Generative graph of the basic Bayesian model.
In the supervised or semi-supervised systems, some parameters are also estimated. The unknown variables contains not only f but also the parameters θ. In this case, Bayes' rule derives the joint posterior probabilistic law :

p(f , θ|g) = p(g|f , θ 1 )p(f |θ 2 )p(θ) p(g) , (3.38)
where vector θ = [θ 1 , θ 2 ] contains all the parameters to be estimated. This joint posterior law can then be used with different Bayesian point estimation methods to infer on f and θ. For a hierarchical structured model where a hidden variable z appears, a posterior law with all the unknown variables is obtained via the Bayes' rule :

p(f , z, θ|g) = p(g|f , θ 1 )p(f |z, θ 2 )p(z|θ 3 )p(θ) p(g) , (3.39)
where θ = [θ 1 , θ 2 , θ 3 ] corresponds to the parameters which are estimated during reconstruction, and κ = [κ 1 , κ 2 , κ 3 ] are the fixed hyper-parameters. This model is illustrated by the generative graph shown in Figure 3.17.

Bayesian point estimators

In Bayesian inference, the posterior distribution is typically considered for the estimation of parameters and variables. The Bayesian point estimators aims at choosing the "best value"

3.7.1 -Bayesian point estimators κ 2 κ 1 θ 2 f g θ 1 p(g|f , θ 1 ) p(f |θ 2 ) p(θ 2 |κ 2 ) p(θ 1 |κ 1 )
FIGURE 3.16 -Generative graph of the semi-supervised Bayesian model.

κ 3 κ 1 κ 2 θ 3 z f g θ 1 θ 2 p(g|f , θ 1 ) p(f |z, θ 2 ) p(z|θ 3 ) p(θ 2 |κ 2 ) p(θ 1 |κ 1 ) p(θ 3 |κ 3 ) FIGURE 3
.17 -Generative graph of the model with hidden variable z.

The Bayesian inference

of the unknown variable from its posterior distribution. They imply the posterior distribution's statistics of central tendency, for example the mean, the median, or the mode.

Consequently, three types of point estimators are derived : Posterior Mean (PM) [START_REF] William A Ericson | A note on the posterior mean of a population mean[END_REF], Posterior Median [START_REF] Bochkina | On the posterior median estimators of possibly sparse sequences[END_REF] and the Maximum a Posterior (MAP) [START_REF] Qi | Resolution and noise properties of MAP reconstruction for fully 3-D PET[END_REF] estimation. For some distributions, there is no explicit form for the median value. Hence typically only PM and MAP estimation are used for the determination of variables is Bayesian inference.

With the posterior distribution obtained from an unsupervised Bayesian inference [START_REF] Mohammad-Djafari | Joint estimation of parameters and hyperparameters in a Bayesian approach of solving inverse problems[END_REF] as in Eq.(3.38), we distinguish three estimation methods. The first is to integrate out θ from p(f , θ|g) to obtain p(f |g) and then use it to infer on f . Unfortunately, this approach does not often give explicit expression for p(f |g). The second approach is that one can first marginalize it with respect to f to obtain p(θ|g) = p(f , θ|g) df and estimate θ = arg max θ {p(θ|g)}, then use it as it was known. The third and easiest algorithm to implement is the joint optimization, which estimates the variable f and parameter θ iteratively and alternately.

In MAP estimation, the unknown variables are estimated by maximizing the posterior distribution. When there are several variables, they are optimized jointly :

(f , θ) = arg max f ,θ {p(f , θ|g)} . (3.40)
The posterior distribution has normally an negative exponential form, and hence the maximization of the posterior distribution are normally considered as a minimization of the negative logarithm of posterior distribution :

(f , θ) = arg min f ,θ {-ln p(f , θ|g)} . (3.41)
In Posterior Mean estimation, when there are several related variables and the posterior distribution is not separable for the variables, the Variational Bayesian Approximation (VBA) [FR12, TLG08, AMD10] is always used. The VBA method approximates the posterior distribution p(f , θ|g) by a separable distribution q(f , θ) by minimizing the Kullback-Leibler divergence [BMD08, CB01], p(f , θ|g) ≈ q(f , θ) = q(f )q(θ). The conjugate prior distributions are normally used to define the prior model in order to simplify the computation of the posterior model and the parameters.

Link between Bayesian and Regularization methods

By using the MAP estimation for the posterior distribution, the links between the Bayesian inference and the regularization methods are derived. For example we consider the following MAP problem :

p(f ) = exp {-R(f )} , (3.42) p(g|f ) = exp {-Q(g, f )} .
(3.43)

The posterior distribution obtained from the likelihood and prior model is :

p(f |g) ∝ p(g|f )p(f ) = exp {-Q(g, f ) -R(f )} .
(3.44)

-Link between Bayesian and Regularization methods

By minimizing the negative logarithm of the posterior distribution, the variable is estimated via :

f = arg max f {p(f |g)} = arg min f {-ln p(f |g)} = arg min f {Q(g, f ) + R(f )} , (3.45)
which is exactly the general model of the regularization method presented in Eq.(3.24).

As in X-ray CT we the noise is modeled by a Gaussian distribution, we have

p(g|f ) = N (g|Hf , V ) ∝ |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) , (3.46)
and we get :

p(f |g) ∝ exp -g -Hf 2 V -R(f ) , (3.47) and f = arg max f g -Hf 2 V + R(f ) . (3.48)
When the prior distribution of f is a Gaussian model, with zero mean and variance v f :

p(f ) = N (f |0, v f ) ∝ exp - 1 2 f 2 V f , (3.49) we get f = arg min f g -Hf 2 V + f 2 V f . We can see that when V = v I M and V f = v f I N
where I is an identity matrix, we obtain the Quadratic Regularization method.

When the prior distribution of f is a Laplacian model :

p(f ) = L(f |0, V f ) ∝ exp - 1 2 f 1V f , (3.50) we get f = arg min f g -Hf 2 V + f 1V f . When V = v -1 I M and V f = v -1 f I N , we obtain the LASSO regularization method [Tib96]. If p(f ) = L(f |0, D T V f D -1 ), we get f = arg min f g -Hf 2
V + Df 1V f and it corresponds to the TV method.

Conclusions

In this chapter, a framework of the X-ray CT reconstruction methods is presented. Among these conventional methods, the analytical and algebraic ones considering the forward model without noise : g = Hf , are not robust to the noise, and no additional information is considered except for the biased data.

The regularization methods and the statistical methods consider an additive noise in the forward model g = Hf + . These methods, especially the Bayesian ones, take into consideration the prior information as a supplement criterion for the ill-posed reconstruction problems. In the Bayesian methods, the information on the object is modeled by some appropriate probability density function (pdf). An important objective in the Bayesian method is to choose suitable pdfs for the unknown variables.

In the next chapter, I am going to give a brief review of the Bayesian methods, and present some conventional used prior models for the piecewise continuous variables in the X-ray CT problems.

4

Bayesian method and sparse enforcing prior distributions

Introduction

In this chapter, we present the Bayesian inference and its utilization in the CT reconstruction problem. As the most important part of the Bayesian methods, the definition of the prior models are discussed in this chapter. The prior model for defining the piecewise continuous images or objects are typically realized by using a Markovian model or by introducing a sparse hidden variable. The definition of the sparse prior model can be realized by using plenty of distributions. In this chapter we present a generalization of the Student-t distribution and discuss its property.

Bayesian inference has provided a useful tool for the inverse problems in different domains, especially in the ill-posed problems where the number of detected data is insufficient. This advantage of the Bayesian inference comes from the fact that the prior information is considered in the method in terms of the prior model of the variables. The prior information are combined with the likelihood and the detected data in the reconstruction.

We recall the Bayes rule for the unsupervised system model :

p(f , θ|g) = p(g|f , θ 1 )p(f |θ 2 )p(θ 2 ) p(g) ∝p(g|f , θ 1 )p(f |θ 2 )p(θ 2 ). (4.1)
In this formula, p(f , θ|g) is the posterior distribution, p(g|f , θ 1 ) is the likelihood, p(f |θ 2 ) and p(θ 2 ) are the prior distributions.

We can see that the posterior distribution depends on the likelihood and the prior distributions. In the Bayesian methods, an essential work is to define the prior models. The prior information of the variable f and parameter θ = [θ 1 , θ 2 ] are considered and described by some suitable distribution. Some examples of the prior information are the "positivity", "continuity",

CHAPITRE 4. BAYESIAN METHOD AND SPARSE ENFORCING PRIOR DISTRIBUTIONS

"sparsity", etc.

In X-ray CT, the objects under consideration are for example the brain or organs in medical applications, or components of industrial objects in Non Destructive Testing (NDT).

The industrial components normally consist of several different homogeneous materials, for example the steel, iron, air, etc. In medical CT, the brains, bones, teeth or other organs also consist of some homogeneous substances. We can also see from the examples shown in Figure 3.1 and Figure 3.2 that the objects under examination are piecewise continuous.

In this thesis, the piecewise continuous property of the object is particularly considered in the Bayesian methods when choosing the prior models. In the following sections I'm going to present the choice of the prior models adapting to the piecewise continuous prior information.

Gradient and Markovian models

The "homogeneous" property implies that the voxels in the object are not independent. Each voxel depends on other voxels. But the coherence between two voxels decreases rapidly when the distance between them increases. Consequently, only the nearest neighbor voxels are related to the current voxel.

Noted by r the current voxel position, and f (r) the value of the voxel, we use r to represent the neighbor voxels of r, r ∈ N (r).

Four or eight neighbors are the most commonly used cases. In order to consider the coherent relations between voxels, the Markovian method is generally used. r r' In many methods, the gradient of the piecewise continuous variable is used as a sparse hidden variable. In the gradient-based methods, the difference between the neighbour pixels are considered. As shown in Figure 4.2, for a piece-wise constant object, the gradient is zero in most of parts except for the contours areas.

In the Markovian models, the strategy is also to consider the neighbour pixels in the prior models of the piecewise continuous variables. There are mainly two types of Markovian models, known as the homogeneous Markovian models and the inhomogeneous Markovian models. In the homogeneous Markovian models, the big variation of the pixel values is penalized for all the pixels, thus the noise is limited. However, the contours in the image are also obscured. The inhomogeneous Markovian models solves this problem by considering separately the homogeneous areas and the contours.

Homogeneous Markovian model

The homogeneous Markovian model defines the same Markovian model for every voxel, with the following prior model expression :

p(f ) = Ψ (Φ (f (r), f (r ))) . (4.2)
The function Φ (f (r), f (r )) describes the influence of the neighbor pixels on the current pixel. In the Figure 4.1, we show the example of considering 4 and 8 neighbor pixels, and the corresponding Φ(f (r), f (r )) of these two cases are :

Φ (f (r)) = f (r) - 1 4 r ∈N 4 (r) f (r ), (4.3) Φ (f (r)) = f (r) - 1 8 r ∈N 8 (r) f (r ), (4.4)
where N 4 (•) and N 8 (•) represent respectively the set of 4 and 8 neighbor pixels.

The prior distribution Ψ with exponential form are commonly used :

p(f ) = Ψ (Φ (f (r), f (r ))) ∝ exp -λ r Φ (f (r), f (r )) 2 2 , (4.5)
where λ is a parameter of this prior distribution.

-Inhomogeneous Markovian models

The Laplacian Markovian model is also used in some methods :

p(f ) = Ψ (Φ (f (r), f (r ))) ∝ exp -λ r Φ (f (r), f (r )) 1 . (4.6)
By using these two prior models in Bayesian method and optimize via the MAP algorithm, we obtain respectively the conventional Quadratic Regularization method (QR) and the Total Variation method (TV) :

QR : J 1 (f ) = g -f 2 2 + λ Df 2 2 , TV : J 2 (f ) = g -f 2 2 + λ Df 1 , (4.7)
where D is the gradient operator, corresponding to Φ (f (r), f (r )).

Inhomogeneous Markovian models

The homogeneous Markovian models are frequently used, but a shortcoming is that these models do not distinguish the homogeneous areas and the contours in the object, thus the contours in the reconstructed image will not be preserved.

The inhomogeneous Markovian model solves this problem. In these models, one or several hidden variables are defined to identify different areas. For example, the label of different homogeneous zones, or the label of contours.

The word "inhomogeneous" implies that the probability model for the pixels are not the same for all the pixels. As shown in Figure 4.3, the pixels f 1 and f 4 have 2 neighbour pixels on the contours, f 2 has only 1 neighbour pixel on the contours, and f 3 has no neighbour pixel on the contours. In the inhomogeneous Markovian model, these four pixels have different prior distributions, depending on the neighbour pixels in the homogeneous areas (circles) while ignoring the neighbour pixels on the contours (triangles). The inhomogeneous Markovian method depending on the label of different blocks is presented in [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF], and the prior model for the object is :

p(f |z) = R N (µ z (R), v z (R)), (4.8)

Sparse transformation base coefficients

where R ∈ R is the label of different blocks, µ z (R) and v z (R) are respectively the mean and variance of the pixel z in label R. From the model we see that in different blocks, f follows a Gaussian distribution with different parameters,

p(f |z) = R N (µ z (R), v z (R)).
The gradient can also be used to control the Markovian model, separating the model for voxels into homogeneous areas and contours areas. A method using the gradient value in the prior model in proposed and presented in Chapter 6.

Sparse transformation base coefficients

The piecewise continuous variables can be represented by a optimally sparse base. A simple example is the wavelet transformation. The wavelet transformation has been used in many works of the X-ray CT reconstruction [BKW96, CKL01, RVJ + 06] where the variable of interest is piecewise continuous.

In Chapter 5, a Bayesian method basing on a prior model concerning the Haar transform coefficients is presented.

With the wavelet transformation, or dictionary base transform, the piecewise continuous prior property of the object become the sparse structured property of the coefficient of the transformation. A prior model for enforcing the sparse structure is then used.

As it was presented in the previous section, the piecewise continuous variable can be represented by another sparse variable, via some linear or non linear transformation operations. We present here some commonly used sparse representations.

The linear wavelet transformations include for example the Haar transformation [START_REF] Radomir | The Haar wavelet transform : its status and achievements[END_REF], the contourlet transformation [START_REF] Minh | The contourlet transform : an efficient directional multiresolution image representation[END_REF], the curvelet transformation [START_REF] Starck | The curvelet transform for image denoising[END_REF], etc. We present here the definition of an orthogonal transformation.

A linear transformation can be expressed as a multiplication by a matrix. Here we define a linear transformation :

z = Df , (4.9)
where D ∈ R M ×N represents the transformation matrix, f ∈ R N ×1 is a vector and z ∈ R M ×1 is the transform coefficient.

An orthogonal matrix is square, thus D ∈ R N ×N , and z ∈ R N ×1 . In linear algebra, an orthogonal matrix is a square matrix with real entries whose columns and rows are orthogonal unit vectors (i.e., orthonormal vectors), and thus

D T D = I, (4.10)
where I is an identity matrix.

The equivalent characterization of the orthogonal matrix is that : its transpose is equal to its inverse :

D T = D -1 . (4.11)
One of the frequently used orthogonal transformation is the wavelet transformation. For example the Discrete Haar transformation. The N × N Haar transformation matrix is :

D = d T 0 , d T 1 , • • • , d T N -1
T , (4.12)

-Inhomogeneous Markovian models

and

D T = [d 0 , d 1 , • • • , d N -1 ] . (4.13)
d k is called the Haar function, and is defined as :

d k (t) = 1 √ N      2 p/2
(q -1)/2 p ≤ t < (q -0.5)/2 p -2 p/2 (q -0.5)/2 p ≤ t < q/2 p 0 otherwise , (4.14) and d 0 (t) = 1/ √ N . p and q are defined as :

k = 2 p + q -1, (4.15)
and 2 p is the largest power of 2 that is smaller than k.

when N = 2 :

D 2 = 1 √ 2 1 1 1 -1 . (4.16)
When N = 4, we have :

D 4 = 1 2     1 1 1 1 1 1 -1 -1 √ 2 - √ 2 0 0 0 0 √ 2 - √ 2     .
(4.17)

In Figure 4.4, the wavelet shapes show the corresponding Haar functions. We can see that the function has a shape as the step function and the period of the wavelet reduces by half when level of p augments by 1. We can see that according to the definition, we have :

d T i * d j = 1 if i = j 0 if i = j (4.18)
Therefore we obtain : There exist many orthogonal transformations, for example the wavelet transformation [RVJ + 06, BKW96], the curvelet transformation [SCD02, DS07], etc. By using the orthogonal transformations in the Bayesian method, the transformation of the operator can be replaced by the inversion of the operator, which is quick and convenient.

DD T =      d T 0 d T 1 . . . d T N -1      [d 0 d 1 • • • d N -1 ] =        1 0 0 • • • 0 0 1 0 • • • 0 0 0 1 • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • 1        , ( 4

Sparsity enforcing prior distributions

As it was mentioned previously, the sparse transformation coefficients of a piecewise continuous variable is used as a hidden variable in many Bayesian methods. In a sparse vector, most of the elements are zero. The sparse transformations are widely used in compressed sensing problems [START_REF] David | Compressed sensing[END_REF][START_REF] Richard | Compressive sensing [lecture notes[END_REF].

Distributions for sparse variables

In Bayesian inference, a sparse structured variable can be defined by some sparsity enforcing prior distributions. Generally, three types of distributions are used to define the sparse variables :

1. The Generalized Gaussian distributions :

p(z) ∝ exp - 1 2 α z β β . (4.21)
For example when β = 1, it becomes Lapacian distribution. 2. The Mixture distributions :

p(z) ∝ k exp - 1 2 α k z 2 2 . ( 4
lim x→∞ e λx p(X > x) = ∞ ∀ λ > 0. (4.23)
The generalized Gaussian distribution or the Mixture distributions which can be used to define the sparse variables are also heavy tailed.

Generalized Student-t distribution

Among all the sparsity enforcing distributions, the Laplacian distribution is frequently considered. However, the shortcoming is that the computation is complex and slow because of the L1 norm optimization. The Student-t distribution is another well known sparsity enforcing distribution. In my work, I use a two-parameter version of the Student-t distribution, benefiting the convenience that it can be derived via a hierarchical model and this model leads to a semisupervised system model for the Bayesian method.

As a heavy tailed distribution, the standard Student-t distribution has a probability density function as :

St(x|ν) = Γ ν+1 2 √ νπΓ ν 2 1 + x 2 ν -ν+1 2 , ν > 0, (4.24)
where ν is the number of degrees of freedom and Γ represents the Gamma function.

The standard Student-t distribution can also be expressed as an Infinite Gaussian Scale Mixture (IGSM) distribution [START_REF] Volker | Mixtures of Gaussian processes[END_REF][START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF] :

St(x|ν) = ∞ 0 N x|0, 1 z G z| ν 2 , ν 2 dz, (4.25)
where N (•) is the Gaussian distribution and G(•) is the Gamma distribution. A Normal-Inverse Gamma marginalization can also be used to define the Student-t distribution :

St(x|ν) = ∞ 0 N (x|0, z) IG z| ν 2 , ν 2 . (4.26)
Comparing with many other heavy tailed distributions, for example the Laplacian distribution and the Generalized Gaussian distribution, the Student-t distribution is particularly interesting thanks to the IGSM expression. What's more, the IGSM expression is composed by a Normal and a Gamma or Inverse Gamma distribution, which are conjugate likelihood and prior pairs. This property facilitates the computation in Bayesian methods.

In our work, we replace the one-parameter Inverse Gamma distribution IG(z| ν 2 , ν 2 ) by the two-parameter Inverse Gamma distribution IG(z|α, β), and we call it the generalized Student-t distribution. It is expressed as :

St g (x|α, β) = ∞ 0 N (x|0, z) IG (z|α, β) dz. (4.27)
And the explicit expression of the generalized Student-t distribution is :

St g (x|α, β) = Γ α + 1 2 √ 2βπΓ(α) 1 + x 2 2β -(α+ 1 2 ) , α > 0, β > 0.
(4.28)

Sparsity enforcing prior distributions

When α = β = ν 2 , we obtain the standard Student-t distribution. In Figure 4.6, the form of the pdf of Gaussian distribution, Student-t distribution and the generalized Student-t distribution are compared. From the comparison of the standard Student-t distribution and the generalized Student-t distribution, we can see that the standard Student-t distribution has a limitation : the shape of the pdf is either heavy tailed or with narrow peak. These two properties can not be realized simultaneously. The generalized Student-t distribution, however, can have a pdf with both narrow peak and a heavy tail, by choosing suitable values for the parameters α and β. In Bayesian inference, the sparse variables are defined belonging to a generalized Studentt distribution. With the Normal-Inverse Gamma marginalization property, a sparse variable x is defined as a Gaussian distributed variable, with mean equal to zero and variance z, and the variance z is defined belonging to an Inverse Gamma distribution with parameters α and β :

p(x|z) = N (x|0, z), p(z|α, β) = IG(z|α, β).
(4.29)

By choosing different values for the (α, β) pair, we get different shape of the probability distributed function. The standard Student-t distribution can be approached to a Gaussian distribution when ν → ∞.

-Generalized Student-t distribution

Demonstration.

lim ν→∞ p(x|ν) = lim ν→∞ Γ ν+1 2 √ νπΓ ν 2 1 + x 2 ν -ν+1 2 .
(4.30)

By using the theorem :

lim ν→∞ Γ(ν + α) Γ(ν)ν α = 1, (4.31)
we have :

lim ν→∞ Γ ν+1 2 √ νπΓ ν 2 = lim ν→∞ Γ ν 2 + 1 2 Γ ν 2 ν 2 1 2 √ 2π = 1 √ 2π . (4.32)
And with the Taylor series we have :

lim ν→∞ 1 + x 2 ν -ν+1 2 = lim ν→∞ exp - ν + 1 2 ln 1 + x 2 ν = lim ν→∞ exp - ν + 1 2 x 2 ν - 1 2
x 4 ν 2 + 1 3

x 6 ν 3 + • • • = lim ν→∞ exp - 1 2 ν + 1 ν x 2 = exp - 1 2 x 2 .
(4.33)

And we get a Gaussian distribution. For the St g distribution, we can also use the same strategy :

lim α,β→∞ Γ α + 1 2 √ 2πβΓ(α) 1 + x 2 2β -(α+ 1 2 ) = lim α,β→∞ 1 √ 2π α β exp -(α + 1 2 ) x 2 2β - 1 2 x 4 4β 2 + • • • = lim α,β→∞ 1 √ 2π α β exp - 1 2 α + 1 2 β x 2 .
(4.34)

When α and β are both big enough and they are defined on the same scale, the approximation is also a Gaussian distribution :

lim α,β→∞ 1 √ 2π α β exp -(α + 1 2 ) x 2 2β - 1 2 x 4 4β 2 + • • • = 1 √ 2π exp - 1 2 x 2 . (4.35)

Conclusion

We have introduced in this chapter the context of the Bayesian inference and discussed the choice of the prior model for the piecewise continuous variable. The generalized Student-t distribution is presented in this chapter. From the discussion, we can see that the St g distribution could be used as the prior law for the sparse hidden variable in the Bayesian method. The orthogonal transformation is convenient as the inversion of operator and the transpose of operator are identical.

Conclusion

In the next two chapters, we present the proposed Bayesian methods basing on two principal strategies.

The first strategy is to define a hidden variable, z, which is a sparse transformation coefficient of the object, f , in the prior model. A hierarchical model is proposed. We call it the Hierarchical Haar based Bayesian Method (HHBM). In this method, both the coefficient z and the object f are considered as unknown variables in the model, and the generalized Student-t distribution is used as the prior distribution for the sparse coefficient z. This proposed method and some extensions are presented in Chapter 5.

The second strategy (ROCC method), presented in Chapter 6, uses an inhomogeneous Markovian model. In this proposed method, the Laplacian of the object is considered as a hidden variable, and is modeled by the generalized Student-t distribution because of its sparse structure. The optimized Laplacian is used as a hyperparameter in the Markovian model for the object, and the optimized object then is used in the optimization of Laplacian. Iteratively and alternately, the object and the Laplacian are optimized.

5

Sparsity Enforcing Unsupervised Hierarchical Model

Introduction

As we have presented previously, the objects that we take into consideration in Non Destructive Testing (NDT) are commonly piecewise continuous. As presented in the previous chapter, the prior models for the piecewise continuous variables are normally a Markovian model or a hierarchical model with a hidden variable. This hidden variable is for example a sparse transformation coefficient of the piecewise continuous variable.

In this chapter, a Bayesian method with a hierarchical structured prior model is proposed. In this model, the coefficient of a sparse orthogonal transformation is considered as a hidden variable. Meanwhile, the semi-supervised system is realized such that the parameters are estimated simultaneously with the unknown variables. The initialization for the hyperparameters is discussed and demonstrated in this chapter. The initialization of hyperparameters is insensible to the final results, which is convenient comparing with the state-of-the-art methods.

The choice of the sparse transformation

In the big data size applications, the transformation is normally realized by a convolution operator, as the multiplication matrix will have a huge size. In this case, the inverse of a transformation is normally not expensive to compute. We use an orthonormal transformation in our work. Thanks to the fact that the inversion and transpose of an orthonomal transformation are identical, the computation with big data size can be simplified, with the transpose being replaced by the inversion.

Wavelets and multiscale analysis have proved to be useful in a wide range of applications Among these transforms, the discrete Haar transformation is practical in the discrete phantoms, and especially for the piecewise constant or continous images. The Haar transform is one of the earliest transformation. As one of the simplest and oldest orthonormal wavelet, the Haar transformation is used in many applications as for example image coding [START_REF] Antonini | Image coding using wavelet transform[END_REF][START_REF] Roger L Claypoole | Nonlinear wavelet transforms for image coding via lifting[END_REF], edge extraction [START_REF] Gerónimo | Haar wavelets and edge orientation histograms for on-board pedestrian detection[END_REF], and recently image restoration [START_REF] Petrou | Image processing : the fundamentals[END_REF].

In the Bayesian method proposed in this chapter, the Haar transformation coefficient is used as a sparse hidden variable. In this section, the details of the choice of the transformation level are presented. In order to analyze the properties of HT of different levels, we present in Figure 5.2 and Figure 5.3 the histograms of the transformation coefficients for the Shepp Logan object and the Head object respectively. From the histogram we can see that, when l is small, the coefficients are less compressed with bigger values. When l increases, the sparsity rate increases. We can see from the figures that, when the transformation level is equal to or bigger than 4, the histograms

The choice of the sparse transformation are identical and hence the sparsity rate are on the same scale. Thus in our work, we choose to use a 5-level HT in order to get a stable sparse structure. As we have presented previously, the sparse variable can be modeled principally by three types of distributions : the generalized Gaussian distributions for example the Laplacian distribution, the Gaussian mixture distributions for example a mixture of two Gaussian distributions and the heavy tailed distributions for example the Student-t distribution. As presented in A synthesis hierarchical Bayesian model Section 4.4, the generalized Student-t distribution is used to model the sparse transformation coefficients.

l = 1 l = 2 l = 3 l = 4 l = 5 l = 6

A synthesis hierarchical Bayesian model

In the published works, the prior model types can be classified into two types : the analysis type and the synthesis type [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. The analysis based approach derives the likelihood of the target signal from a forward transform applied to it. An example is the gradient in TV method. This analysis type commonly appears as a regularization term in optimization :

f = arg min f {J(f )} with J(f ) = g -Hf 2 2 + λR (Df ) .
(5.1)

The synthesis type of prior, on the other hand, changes the signal of consideration into the coefficient of a forward transform which will be estimated. The coefficient then will be used in a post processing to obtain the final result of the target signal. An example is as following : with the forward model g = Hf + and the forward transform f = Dz, z is estimated by optimizing :

f = D z with z = arg min z {J(z)} with J(z) = g -HDz 2 2 + λR(z), (5.2)
and f = D z is derived as the final result of the reconstruction.

First of all, we consider a synthesis model with which we obtain directly the coefficient of the transformation.

The Hierarchical Model.

Again, we mention here the forward system model :

g = Hf + , (5.3)
where is an additive noise representing the uncertainties of the projection system. In the X-ray CT models, it is modeled by a Gaussian distribution, with mean equal to zero :

p( |v ) = N ( |0, v I).
(5.4)

The variance of the distribution, v , is a fixed scalar in a basic model, with the hypothesis of a stationary noise.

In most of the cases, the variance of the additive noise in the forward model, Eq.(5.3), is unknown. In that case, we consider it as an unknown variable and a prior probability distribution is assigned to it.

To choose a prior model for the variance of a Gaussian likelihood, a conjugate distribution is always preferred. Taking into consideration of the positivity of the variance, the Gamma or Inverse Gamma distributions are commonly used. In our work, we use the Inverse Gamma distribution to take advantage of the Normal-Inverse Gamma bivariate distribution property, which will be presented later.

-The Hierarchical Model.

Définition. The Inverse Gamma (IG) distribution's probability density function over the support x > 0 is :

p(x|α, β) = β α Γ(α) x -α-1 exp - β α , (5.5)
with shape parameter α and rate parameter β. Γ(•) denotes the gamma function :

Γ(x) = ∞ 0
t x e -t dt.

(5.6)

As we can see in Figure 5.4 which presents the pdf of the inverse gamma distribution, the variable is always positive, and most of the values are centralized at small values near zero, which exactly adapt to the property of the variance. Indeed, the inverse gamma distribution is a conjugate distribution for a Gaussian likelihood with unknown variance, the Bayesian inference methods are therefore simplified with a simple and regular posterior distribution. The inverse gamma model for v is denoted by :

p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ). (5.7)
As is presented in Chapter 4, a Gaussian distribution, Eq.(5.4), with the variance belonging to an inverse gamma distribution, Eq.(5.7), leads to a generalized Student-t distribution, which can either approach a Gaussian distribution or tending to a heavy tailed sparse enforcing distribution depending on the choice of the two hyperparameters α 0 and β 0 .

For a more general case, we consider a non-stationary noise, thus for each measurement the variance of the noise is different. This appears for example in X ray CT for a slab material.

A synthesis hierarchical Bayesian model

The noise will be more important in the projection direction along the longer side of the object, as there will be more absorption and scattering of photons when the radiations pass through the object. In these cases, the variance of noise is defined as a vector, each element of which corresponds to one projection direction. Therefore the prior system model and the noise model are :

p(g|f , v ) = N (g|Hf , V ), where V = diag [v ] ,
(5.8) p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ),

(5.9)

where v = [v 1 , • • • , v M ].
By supposing the independence of v i , ∀i, we obtain :

p(v |α 0 , β 0 ) = M i=1 IG(v i |α 0 , β 0 ).
(5.10)

Each element v i belongs to an inverse gamma distribution.

Now, a prior model is defined for the additive noise . For the reconstruction of f , the most important model to be defined is the prior model of the object f . As it is presented previously, in the X ray CT the phantoms under consideration are generally piecewise continuous, and the discrete multilevel Haar transform coefficient is a suitable representation of the phantom. We take into consideration the sparsity of the Haar transform coefficient with the relation :

f = Dz, (5.11)
where z is the multilevel Haar transform coefficient of f and D -1 is the transformation operator, with z = D -1 f . As this transformation operator is orthogonal with D -1 = D T , and we have z = D T f .

With this relation, we could replace the variable f in Eq.(5.3) by the transformation coefficient z : g = HDz + .

(5.12) By doing so, using the MAP estimation will lead to exactly the same criterion as the synthesis method, shown in Eq. (5.2).

As the Haar transformation coefficient, the variable z is sparse. Consequently, we use the generalized Student-t distribution, presented in Chapter 4, to define each element z j :

p(z j |α z 0 , β z 0 ) = Γ α z 0 + 1 2 2πβ z 0 Γ (α z 0 ) 1 + z 2 j 2β z 0 -(αz 0 + 1 2 ) .
(5.13) Thanks to the Gaussian-Inverse Gamma marginal structure of the St g distribution, we could use a hierarchical model to define the hidden variable :

p(z|v z ) = N (z|0, V z ), where V z = diag [v z ] , p(v z |α z 0 , β z 0 ) = N j=1 IG(v z j |α z 0 , β z 0 ), (5.14)
where

v z = [v z 1 , • • • , v z N ],
and its elements are supposed to be a priori independent and identically distributed.

-The Hierarchical Model.

Now we have all the prior models with the hyper-parameters (α 0 , β 0 ) and (α z 0 , β z 0 ) :

p(g|z, v ) = N (g|HDz, V ) ∝ |V | -1 2 exp - 1 2 (g -HDz) T V -1 (g -HDz) , p(z|v z ) = N (z|0, V z ) ∝ |V z | -1 2 exp - 1 2 z T V -1 z z , p(v |α 0 , β 0 ) = M i=1 IG(v i |α 0 , β 0 ) ∝ M i=1 v -(α 0 +1) i exp -β 0 v -1 i , p(v z |α z 0 , β z 0 ) = N i=1 IG(v z j |α z 0 , β z 0 ) ∝ N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j .
(5.15)

Figure 5.5 present the generative graph of the hierarchical system. In this figure, the circles stand for the unknown variables and hidden variables, and the squares stand for the fixed parameters and hyperparameters.

α z 0 , β z 0 α 0 , β 0 v z z g v HD p(g|z, v ) p(z|v z ) p(v z |α z 0 , β z 0 ) p(v |α 0 , β 0 ) FIGURE 5
.5 -Generative graph of the synthesis model.

Then we use the Bayes rule to obtain the posterior distribution :

p(f , v , v z |g) ∝ p(g|z, v )p(z|v z )p(v z )p(v )
(5.16)

A synthesis hierarchical Bayesian model

The posterior distribution is obtained by taking into consideration Equations (5.15) :

p(z, v z , v |g) = |V | -1 2 exp - 1 2 (g -HDz) T V -1 (g -HDz) • |V z | -1 2 exp - 1 2 z T V -1 z z M i=1 v -(α 0 +1) i exp -β 0 v -1 i • N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j = M i=1 v -1 2 i exp - 1 2 M i=1 v -1 i [g -HDz] 2 i • N j=1 v -1 2 z j exp - 1 2 N j=1 v -1 z j z 2 j M i=1 v -(α 0 +1) i M i=1 exp -β 0 v -1 i • N j=1 v -(αz 0 +1) z j N j=1 exp -β z 0 v -1 z j .
(5.17)

From the posterior distribution we can apply a Bayesian point estimation to solve the reconstruction problem.

Bayesian point estimation

The commonly used Bayesian point estimations are the Joint Maximum A Posterior (JMAP) estimation and the Posterior Mean (PM) estimation via Variational Bayesian Approach (VBA). JMAP estimation maximizes the posterior distribution with respect to each variable alternately. This algorithm is convenient for the big data size problems. Meanwhile, the PM algorithm is more complicated. By using a VBA method, the computation could be simplified, but the computational costs are still very expensive for the simulations of 3D phantom. In my work, the JMAP method is used.

In 3D phantom reconstruction, the size of data is so big that the elements of the matrix H are not accessible. The only things that we have access to are the geometric operator of projection Hf and the geometric operator of back-projection H T g. MCMC method is theoretically feasible but is too costly because of the huge data size and the huge amount of variables.

The Joint MAP estimation.

In the JMAP estimation, the maximization of the posterior distribution is changed to the minimization of the negative logarithm of the posterior distribution. This comes from the fact that the distributions are always in the form of an exponential of a negative polynomial, and that the logarithm operation is monotonous :

( z, v z , v ) = arg max z,vz,v {p (z, v z , v |g)} = arg min z,vz,v {-ln p (z, v z , v |g)} .
(5.18)

-The Joint MAP estimation.

From Eq.(5.17), we calculate the criterion of the optimization :

J(z, v z , v ) = -ln p (z, v z , v |g) = 1 2 M i=1 ln v i + 1 2 M i=1 v -1 i [g -HDz] 2 i + 1 2 N j=1 ln v z j + 1 2 N j=1 v -1 z j z 2 j + (α 0 + 1) M i=1 ln v i + β 0 M i=1 ln v i + (α z 0 + 1) N j=1 ln v z j + β z 0 N j=1 ln v z j .
(5.19)

In order to estimate three unknown variables from one joint criterion, we use an alternate optimization method, which estimates each of the unknown variable alternately. When optimizing one unknown variable, the others are fixed. It is called the Joint Maximum A Posteriori (JMAP) estimation in this thesis.

The optimization criterion for variables are respectively :

J(z) = 1 2 M i=1 v -1 i [g -HDz] 2 i + 1 2 N j=1 v -1 z j z 2 j ,
(5.20)

J(v z ) = 1 2 N j=1 ln v z j + 1 2 N j=1 v -1 z j z 2 j + (α z 0 + 1) N j=1 ln v z j + β z 0 N j=1 ln v -1 z j , (5.21)
J(v ) = 1 2 M i=1 ln v i + 1 2 M i=1 v -1 i [g -HDz] 2 i + (α 0 + 1) M i=1 ln v i + β 0 M i=1 ln v -1 i .
(5.22)

By considering v z j , ∀j ∈ [1, N ] and v i , ∀i ∈ [1, M ] separately accounting for the hypothesis that each of them are independent, we obtain :

J(z) = 1 2 M i=1 v -1 i [g -HDz] 2 i + 1 2 N j=1 v -1 z j z 2 j , (5.23)
J(v z j ) = α z 0 + 3 2 ln v z j + β z 0 + 1 2 z 2 j v -1 z j , (5.24) J(v i ) = α 0 + 3 2 ln v i + β 0 + 1 2 [g -HDz] 2 i v -1 i .
(5.25)

The gradient of ∇J(z) is :

∇J(z) = d 1 2 (g -HDz) T V -1 (g -HDz) + 1 2 z T V -1 z z dz = 1 2 -D T H T V -1 g -D T H T V -1 g + 2D T H T V -1 HDz + 1 2 • 2V -1 z z = -D T H T V -1 g + D T H T V -1 HDz + V -1 z z = -D T H T V -1 g + D T H T V -1 HD + V -1 z z.
(5.26) A local minimum value z appears when ∇J(z)| z = 0 :

-D T H T V -1 g + D T H T V -1 HD + V -1 z z = 0.
(5.27)

A synthesis hierarchical Bayesian model

And so :

z = D T H T V -1 HD + V -1 z -1 D T H T V -1 g.
(5.28) However, this doesn't work in the 3D X ray CT problems. As in the huge data size case, it is impossible to calculate the inversion of the matrix D T H T V -1 HD + V -1 z , needless to say to do this computation in each iteration. Some algorithms are proposed and used to solve this problem, for example the gradient descent algorithm [Bat92, VCR97], the conjugate gradient algorithm [FB99, Wu01], etc. In our work we used the gradient descent method to optimize the variable, by iteratively computing :

z (k+1) = z (k) -γ (k) z ∇J( z (k) ), (5.29)
where ∇J( z (k) ) = dJ(z) dz is the gradient of the criterion and γ

(k) z is the descent step length. γ (k) z
can either be fixed as a constant or change in every iteration adapting the gradient. The last case is called the optimized step length strategy [START_REF] Boyd | Convex optimization[END_REF], in which γ

(k)
z is calculated in each iteration by minimizing the criterion :

γ (k) z = arg min γz J( z (k) -γ z ∇J( z (k) )) , (5.30)
where J(z -γ z ∇J z ) is simplified as :

J(z -γ z ∇J(z)) = 1 2 (g -HD (z -γ z ∇J(z))) T V -1 (g -HD (z -γ z ∇J(z))) + 1 2 (z -γ z ∇J(z)) T V -1 z (z -γ z ∇J(z)) = 1 2 g T -z T D T H T + γ z ∇J(z) T D T H T V -1 (g -HDz + γ z HD∇J(z)) + 1 2 z T -γ z ∇J(z) T V -1 z (z -γ z ∇J(z)) = 1 2 γ z g T V -1 HD∇J(z) -γ z z T D T H T V -1 HD∇J(z) + γ z ∇J(z) T D T H T V -1 g -γ z ∇J(z) T D T H T V -1 HDz + γ 2 z ∇J(z) T D T H T V -1 HD∇J(z) + 1 2 -γ z z T V -1 z ∇J(z) -γ z ∇J(z) T V -1 z z + γ 2 z ∇J(z) T V -1 z ∇J(z) + Cst = 1 2 γ z 2∇J(z) T D T H T V -1 g -2∇J(z) T D T H T V -1 HDz -2∇J(z) T V -1 z z + 1 2 γ 2 z ∇J(z) T D T H T V -1 HD∇J(z) + ∇J(z) T V -1 z ∇J(z) + Cst =γ z ∇J(z) T D T H T V -1 g -D T H T V -1 HDz -V -1 z z + 1 2 γ 2 z Y HD∇J(z) 2 2 + Y z ∇J(z) 2 2 + Cst, where Y = V -1 2 and Y z = V -1 2
z . From Eq.(5.26) we have :

J(z) = -D T H T V -1 (g -HDz) + V -1 z z = -D T H T V -1 g + D T H T V -1 HDz + V -1 z z.

-The Joint MAP estimation.

So that :

J(z -γ z ∇J(z)) = -γ z ∇J(z) T ∇J(z) + 1 2 γ 2 z Y HD∇J(z) 2 2 + Y z ∇J(z) 2 2 + Cst = -γ z ∇J(z) 2 2 + 1 2 γ 2 z Y HD∇J(z) 2 2 + Y z ∇J(z) 2 2 .
The gradient of the criterion is :

dJ(z -γ z ∇J(z)) dγ z = -∇J(z) 2 2 + γ z Y HD∇J(z) 2 2 + Y z ∇J(z) 2 2 .
By setting dJ(z -γ z ∇J(z))

dγ z γ (k) z = 0,
we get :

γ (k) z = ∇J(z) 2 2 Y HD∇J(z) 2 2 + Y z ∇J(z) 2 2 .
Therefore, the gradient descent optimization for variable z is summarized as following :

z (k+1) = z (k) -γ (k) z ∇J( z (k) ), (5.31)
γ (k) z = ∇J(z) 2 2 Y HD∇J(z) 2 2 + Y z ∇J(z) 2 2 , (5.32)
with

J(z) = 1 2 M i=1 v -1 i [g -HDz] 2 i + 1 2 N j=1 v -1 z j z 2 j ,
(5.33)

∇J(z) = -D T H T V -1 g + D T H T V -1 HDz + V -1 z z.
(5.34) Back to the Eq.(5.21) and Eq.(5.22), we calculate the gradient of these two criterions :

∇J(v z j ) = α z 0 + 3 2 v -1 z j -β z 0 + 1 2 z 2 j v -2 z j , (5.35) ∇J(v i ) = α 0 + 3 2 v -1 i -β 0 + 1 2 [g -HDz] 2 i v -2 i .
(5.36)

Now we discuss if the local minimum is the global minimum. In order to consider both of these two variables, we consider a general case, for the criterion of the form

J(x) = a ln x + bx -1 ,
where a > 0, b > 0 and x ≥ 0. First of all, we look for the point where the gradient of the criterion equals to zero :

∇J(x) = a x - b x 2 = 1 x a - b x ,
where x = 0 is not defined in this case, so the only zero point is when x 0 = b a .

A synthesis hierarchical Bayesian model

For x > x 0 = b/a, we have :

1 x < a b =⇒ b x < a =⇒ - b x > -a =⇒ a - b x > 0 =⇒ 1 x a - b x > 0 =⇒ ∇J(x) > 0. For 0 < x < x 0 = b/a 1 x > a b =⇒ b x > a =⇒ - b x > -a =⇒ a - b x < 0 =⇒ 1 x a - b x < 0 =⇒ ∇J(x) < 0.
And when x → 0 + =⇒ ∇J(x) → -∞, where 0 + is a positive value infinitely approaching zero.

We have demonstrated that when x < x 0 the gradient is always negative and when x > x 0 the gradient is always positive. So that the critical value x 0 = b/a is exactly the global minimum of the criterion. We get :

v z j = β z 0 + 1 2 z 2 j α z 0 + 3 2 , (5.37)
v i = β 0 + 1 2 [g -HDz] 2 i α 0 + 3 2 .
(5.38)

Algorithm and implementation

The variables z, v z j and v i are optimized alternately and iteratively by using Eq.(5.31)-Eq.(5.34), Eq.(5.37) and Eq.(5.38). The computation algorithm for updating the variable z is shown in Algorithm 1.

Algorithm 1 Gradient descent algorithm for z 1: Input : H, D, g 2: Output : z (k+1) 3: Initialization : z (0) 4: k = 0 5: repeat 6:

k = k + 1 7:
Calculate ∇J( z (k) ) according to Eq.(5.34) After the estimation of the variable z, a post processing is applied in order to obtain the reconstruction of phantom : f = D z.

(5.39)

-Simulation results

Algorithm 2 Summary of the JMAP algorithm for the synthesis Bayesian hierarchical method 1: Fix parameters α z 0 , β z 0 , α 0 , β 0 2: Input : H, D, g 3: Output : z, v z , v 4: Initialization : z (0) = D -1 f F BP , where f F BP is the result of the FBP method. 5: repeat 6:

Update z according to Algorithm 1 7:

Compute v z according to Eq.(5.37) 8:

Compute v according to Eq.(5.38) 9: until convergence or max iteration

Simulation results

We use the projections dataset of the Shepp Logan phantom of size 256 3 to reconstruct the phantom. The dataset is obtained respectively by using : 128 projections with a high SNR=40dB, 128 projections with a low SNR=20dB, 64 projections with a high SNR=40dB and 64 projections with a low SNR=20dB. From the results, we can see that from the reconstructed coefficient, we obtain a reconstruc-

HHBM method

ted phantom with block artifacts. This phenomena comes from the fact that, the inverse wavelet transformation can only reconstruct to a perfect object when the coefficient is not contaminated. We can see that the artifacts has square shapes as the Haar transformation considers only two wavelet direction.

Obviously, this result with the square artifacts is unsatisfactory. We propose another Bayesian method basing on this method. The new Bayesian method estimates the object f and the transform coefficient z simultaneously. The estimated coefficient is used during the reconstruction of the object in each iteration. This proposed method is named the HHBM (Haar transformation based Hierarchical Bayesian Method) in our work.

HHBM method

From the simulation results of the previously presented method, we find that the final result of the reconstruction contains a lot of artifacts. This phenomenon comes from the last step of the process : the transform from z to f . A small error on some of the coefficients, z j , affects the image f greatly, in particular when z j is in the upper level. So we propose another Bayesian method, basing on the previous method, but reconstruct the phantom f simultaneously while estimating the coefficient z and other variables.

The system forward model is the same : g = Hf + .

(5.40) However, instead of using an exact discrete multilevel Haar transform from f to z, we consider a relaxed relation :

f = Dz + ξ, (5.41)
where ξ is an additive noise representing uncertainties of the transformation. While choosing the prior model for ξ, the objective is that f and Dz are as approximate as possible. Therefore we define ξ as a sparse noise with most of the values approaching to zero. The generalized Student-t distribution is used to model ξ, enforcing its sparse structure. The St g distribution for ξ is defined via the hierarchical model :

p(ξ|v ξ ) = N (ξ|0, V ξ ) ∝ |V ξ | -1 2 exp - 1 2 ξ T V -1 ξ ξ , (5.42) p(v ξ |α ξ 0 , β ξ 0 ) = N j=1 IG v ξ j |α ξ 0 , β ξ 0 ∝ N j=1 v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j , (5.43)
where With this forward transform model and the prior model of ξ, we have :

V ξ = diag [v ξ ], v ξ = v ξ 1 , • • • , v ξ j , • • • , v ξ N ,
p(f |z, v ξ ) =N (f |Dz, V ξ ) ∝ |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) .
(5.44)

-Simulation results

Vector z = [z 1 , z 2 , • • • , z N ] represents the discrete multi-level Haar transform coefficient of piece-wise continuous phantom f . We use the generalized Student-t distribution (St g ) [START_REF] Dumitru | A Bayesian approach for periodic components estimation for chronobiological signals[END_REF] to model z, enforcing the sparsity property :

p(z|v z ) = N (z|0, V z ),
(5.45)

p(v z |α z 0 , β z 0 ) = N j=1
IG(v z j |α z 0 , β z 0 ).

(5.46)

Thanks to the fact that the Inverse Gamma distribution is conjugate for a Gaussian likelihood with known experience and unknown variance, using this St g distribution simplifies the calculations when using the Bayesian point optimization methods, for example the Posterior Mean via Variational Bayesian Approximation (VBA) method [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF].

We list out all the proposed likelihood and prior models depending on the prior property of each variable :

p(g|f , v ) ∝ |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) (5.47) p(f |z, v ξ ) ∝ |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) , (5.48)
p(z|v z ) ∝ |V z | -1 2 exp - 1 2 z T V -1 z z , (5.49)
p(v z |α z 0 , β z 0 ) ∝ N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j , (5.50) p(v |α 0 , β 0 ) ∝ M i=1 v -(α 0 +1) i exp -β 0 v -1 i , (5.51) p(v ξ |α ξ 0 , β ξ 0 ) ∝ N j=1 v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j . (5.52) In this model, v z j , v i and v ξ j , ∀i ∈ [1, M ] , ∀j ∈ [1, N ]
, are supposed to be independent and identically distributed (iid). The generative graph of this hierarchical structured model is given in Figure 5.7. By using the Bayes rule, Eq.(3.39), the posterior distribution is obtained :

p(f , z, v , v ξ , v z |g) = p(g, f , z, v , v ξ , v z ) p(g) = p(g|f , v )p(f |z, v ξ )p(z|v z )p(v z )p(v )p(v ξ ) p(g) ∝ p(g|f , v )p(f |z, v ξ )p(z|v z )p(v z )p(v )p(v ξ ).
(5.53) Substituting the prior distribution models in Eq.(5.47)-Eq.(5.52), we get the posterior distri-HHBM method

α z 0 , β z 0 α 0 , β 0 α ξ 0 , β ξ 0 v z z f g v v ξ D H p(g|f , v ) p(f |z, v ξ ) p(z|v z ) p(v z |α z 0 , β z 0 ) p(v ξ |α ξ 0 , β ξ 0 ) p(v |α 0 , β 0 ) FIGURE 5.7 -Generative graph of the HHBM method. bution : p(f , z, v z , v , v ξ |g) ∝ |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) |V z | -1 2 exp - 1 2 z T V -1 z z N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j M i=1 v -(α 0 +1) i exp -β 0 v -1 i N j=1 v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j = M i=1 v -1 2 i exp - 1 2 (g -Hf ) T V -1 (g -Hf ) N j=1 v -1 2 ξ j exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) N j=1 v -1 2 z j exp - 1 2 z T V -1 z z N j=1 v -(αz 0 +1) z j N j=1 exp -β z 0 v -1 z j M i=1 v -(α 0 +1) i M i=1 exp -β 0 v -1 i N j=1 v -(α ξ 0 +1) ξ j N j=1 exp -β ξ 0 v -1 ξ j .
(5.54)

As the reasons presented in the previous section, the JMAP estimation algorithm is used to optimize all the unknown variables.

-Joint Maximum A Posterior estimation

Joint Maximum A Posterior estimation

The Joint Maximum A Posterior (JMAP) method estimate all the unknown variables by maximizing the posterior distribution obtained below :

f , z, v z , v , v ξ = arg max f ,z,vz,v ,v ξ {p(f , z, v z , v , v ξ |g)} .
(5.55)

As the posterior distribution has the form of exponential distribution, it becomes a minimization of the negative logarithm of the posterior distribution :

f , z, v z , v , v ξ = arg min f ,z,vz,v ,v ξ {-ln p(f , z, v z , v , v ξ |g)} = arg min f ,z,vz,v ,v ξ {J(f , z, v z , v , v ξ )} .
(5.56) So, the criterion of optimization is :

J(f , z, v z , v , v ξ ) = -ln p(f , z, v z , v , v ξ |g) = 1 2 M i=1 ln v i + 1 2 (g -Hf ) T V -1 (g -Hf ) + 1 2 N j=1 ln v ξ j + 1 2 (f -Dz) T V -1 ξ (f -Dz) + 1 2 N j=1 ln v z j + 1 2 z T V -1 z z + (α z 0 + 1) N j=1 ln v z j + β z 0 N j=1 v -1 z j + (α 0 + 1) M i=1 ln v i + β 0 M i=1 v -1 i + (α ξ 0 + 1) N j=1 ln v ξ j + β ξ 0 N j=1 v -1 ξ j .
(5.57)

The variables are estimated alternately and iteratively. The optimization criterion for each variable consist of the terms containing the corresponding variable :

J(f ) = 1 2 (g -Hf ) T V -1 (g -Hf ) + 1 2 (f -Dz) T V -1 ξ (f -Dz) , J(z) = 1 2 (f -Dz) T V -1 ξ (f -Dz) + 1 2 z T V -1 z z, J(v z ) = 1 2 N j=1 ln v z j + 1 2 z T V -1 z z + (α z 0 + 1) N j=1 ln v z j + β z 0 N j=1 v -1 z j , J(v ) = 1 2 M i=1 ln v i + 1 2 (g -Hf ) T V -1 (g -Hf ) + (α 0 + 1) M i=1 ln v i + β 0 M i=1 v -1 i , J(v ξ ) = 1 2 N j=1 ln v ξ j + 1 2 (f -Dz) T V -1 ξ (f -Dz) + (α ξ 0 + 1) N j=1 ln v ξ j + β ξ 0 N j=1 v -1 ξ j .
(5.58)

HHBM method

And the optimization results are obtained by minimizing the criterions :

f = arg min f {J(f )} , z = arg min z {J(z)} , v z = arg min vz {J(v z )} , v = arg min v {J(v )} , v ξ = arg min v ξ {J(v ξ )} .
(5.59)

The minimization locate at the critical points, where the gradient of the criterion equal to zero. In order to obtain these critical values, the gradient of the criterion are calculated :

∇J(f ) = -H T V -1 (g -Hf ) + V -1 ξ (f -Dz) ,
(5.60)

∇J(z) = -D T V -1 ξ (f -Dz) + V -1 z z.
(5.61)

And for vectors v z , v and v ξ , each element of them are considered separately because of the independence hypothesis :

∇J(v z j ) = 1 2 v -1 z j - 1 2 z 2 j v -2 z j + (α z 0 + 1) v -1 z j -β z 0 v -2 z j , (5.62)
∇J(v i ) = 1 2 v -1 i - 1 2 [g -Hf ] 2 i v -2 i + (α 0 + 1) v -1 i -β 0 v -2 i , (5.63)
∇J(v ξ j ) = 1 2 v -1 ξ j - 1 2 [f -Dz] 2 j v -2 ξ j + (α ξ 0 + 1) v -1 ξ j -β ξ 0 v -2 ξ j .
(5.64) From Eq.(5.60) and Eq.(5.61) :

∇J(f ) f = 0 =⇒ f = H T V -1 H + V -1 ξ -1 H T V -1 g + V -1 ξ Dz , (5.65) ∇J(z) | z = 0 =⇒ z = D T V -1 ξ D + V -1 z -1 D T V -1 ξ f .
(5.66)

As we have mentioned previously, the inversion of the matrix

H T V -1 H + V -1 ξ and matrix D T V -1 ξ D + V -1 z
are impossible in the reconstruction of 3D objects because of the huge data size. The gradient descent algorithm is used to realize the optimization for these two variables by computing :

f (k+1) = f (k) -γ (k) f ∇J( f (k) ), z (k+1) = z (k) -γ (k) z ∇J( z (k) ), (5.67)
where γ

(k) f and γ (k)
z are the gradient descent step length and are obtained by minimizing J( f

(k) - γ f ∇J( f (k)
)) and J( z (k) -γ z ∇J( z (k) )) respectively :

J(f -γ f ∇J(f )) =γ f ∇J(f ) T H T V -1 (g -Hf ) -V -1 ξ (f -Dz) + 1 2 γ 2 f Y H∇J(f ) 2 2 + Y ξ ∇J(f ) 2 2 + Constant, where Y = V -1 2 and Y ξ = V -1 2 ξ .
(5.68)

-Joint Maximum A Posterior estimation

Substituting Eq.(5.60) we obtain :

J(f -γ f ∇J(f )) = arg min γ f -γ f ∇J(f ) T ∇J(f ) + 1 2 γ 2 f Y H∇J(f ) 2 2 + Y ξ ∇J(f ) 2 2 + Cst .
(5.69)

The gradient of the criterion is :

∇J(f -γ f ∇J(f )) = -∇J(f ) 2 2 + γ f Y H∇J(f ) 2 2 + Y ξ ∇J(f ) 2 2 .
By setting

∇J(f -γ f ∇J(f )) γ (k) f = 0,
we get :

γ (k) f = ∇J(f ) 2 2 Y H∇J(f ) 2 2 + Y ξ ∇J(f ) 2 2 , where Y = V -1 2 and Y ξ = V -1 2 ξ .
(5.70)

For γ

(k)
z the same algorithm is used :

J(z -γ z ∇J(z)) =γ z ∇J(z) T D T V -1 ξ (f -Dz) -V -1 z z + 1 2 γ 2 z Y ξ D∇J(z) 2 2 + Y z ∇J(z) 2 2 + Cst, where Y z = V -1 2 z .
(5.71) Substituting Eq.(5.61) we have :

J(z -γ z ∇J(z)) = arg min γz -γ z ∇J(z) 2 2 + 1 2 γ 2 z Y ξ D∇J(z) 2 2 + Y z ∇J(z) 2 2 + Cst .
(5.72)

The gradient of the criterion is :

∇J(z -γ z ∇J(z)) = -∇J(z) 2 2 + γ z Y ξ D∇J(z) 2 2 + Y z ∇J(z) 2 2 .
(5.73)

By setting ∇J(z -γ z ∇J(z)) γ (k) z = 0,
we get :

γ (k) z = ∇J(z) 2 2 Y ξ D∇J(z) 2 2 + Y z ∇J(z) 2 2
, where

Y z = V -1 2 z .
(5.74)

The optimizations for f and z by using gradient descent algorithm are :

for k = 1 → I 2 : f (k+1) = f (k) -γ (k) f ∇J( f (k)
), (5.75)

for k = 1 → I 2 : z (k+1) = z (k) -γ (k) z ∇J( z (k) ), (5.76)
where I 2 is the number of iterations for the gradient descent algorithm.

The optimization for the other variables : v z j , v i and v ξ j are obtained by minimizing the corresponding criterion. From Eq.(5.62)-Eq.(5.64), and considering the discussion of the global HHBM method minimum point for the criterion of form a ln x + bx -1 in Section 5.3.3, the critical point of each criterion is the minimum :

∇J(v z j ) vz j = 0 =⇒ v z j = β z 0 + 1 2 v -1 z j α z 0 + 3 2 , (5.77) ∇J(v i ) v i = 0 =⇒ v i = β 0 + 1 2 [g -Hf ] 2 i α 0 + 3 2 , (5.78) ∇J(v ξ j ) v ξ j = 0 =⇒ v ξ j = β ξ 0 + 1 2 [f -Dz] 2 j α ξ 0 + 3 2 , (5.79) ∀i ∈ [1, M ] and ∀j ∈ [1, N ].
The algorithms of the implementation of estimation for f and z are shown in Algorithm 3 and Algorithm 4.

Algorithm 3 Gradient descent algorithm for f 1: Input : H, D, g, z, f 0 2: Output : f 3: Initialization :

f ini = f 0 4: k = 0 5: repeat 6: k = k + 1 7: Calculate ∇J( f (k-1)
) according to Eq.(5.65) Update f Update

(k) = f (k-1) -γ (k) f ∇J( f (k-1)
z (k) = z (k-1) -γ (k)
z ∇J( z (k-1) ) 10: until convergence or maximum iteration An integrated algorithm concerning the estimation of all the unknown variables is given in Algorithm 5.

Analysis of hyper-parameters and initialization

In the regularization methods, there are always one or several regularization parameters. One of the shortcoming of the regularization method is that the choice of these parameters is

-Analysis of hyper-parameters and initialization

Algorithm 5 Summary of the JMAP algorithm for HHBM method 1: Fix parameters α z 0 , β z 0 , α 0 , β 0 , α ξ 0 , β ξ 0 , l 2: Input : H, D, g 3: Output : f , z, v z , v , v ξ 4: Initialization : 5: f Compute v ξ j , ∀j ∈ [1 : N ] according to Eq.(5.79) 14: until convergence or maximum iteration crucial in the reconstruction, meanwhile, for different size and different quality of the measured dataset, the optimal value of this parameter normally changes. Thus, for each different case, one should always consider to initialize with suitable values for the parameters.

In the proposed Bayesian method, however, the case is different. Even though we should always initialize the hyperparameters, we get a theoretical basis to fix the initialization values of hyperparameters. The method is therefore semi-supervised, which is a priority property in the ill-posed reconstruction problems.

The initialization for the hyperparameters is considered based on the hierarchical prior model chosen for the variables. For and ξ, both of them are modeled by a Gaussian distribution with an unknown variance belonging to an Inverse Gamma distribution. is supposed to be a Gaussian noise and ξ is supposed to be sparse such that f and z are satisfying the transformation relationship for most of elements. Hyperparameters α 0 , β 0 , α ξ 0 , β ξ 0 are initialized according to these properties. The hidden variable z is sparse, and is defined by the St g distribution. The hyperparameters α z 0 and β z 0 are initialized based on this sparse enforcing prior model.

Fixing α z 0 and β z 0

As we can see in Figure 5.8, there are several different ranks of transformation coefficient in the multilevel discrete Haar transformation coefficient. We use r ∈ [1, l + 1] to represent the rank of the coefficient of the l-level Haar transformation (shown in Figure 5.8). The low frequency part of the coefficient z, in the top-left of the coefficient Figure 5.8 with r = 1, is less sparse than the other coefficients. The sparsity rate augments when when r is bigger. Consequently, we set a different (α z 0 , β z 0 ) pairs for coefficients in different coefficient ranks. For (α z 0 (j), β z 0 (j)) in rank r, we set (α r , β r ). The strategy of the initialization is to fix one hyperparameter, α z 0 , and use different values for hyperparameter β z 0 , to realize the definition of different sparsity rate for different ranks. We set α z 0 = 2.01, so that the variance of the St g distribution will be approximate to β z 0 according to the definition Var[z|α z 0 , β z 0 ] = β z 0 /(α z 0 -1). When α z 0 is fixed, we get a higher sparsity rate with a smaller value for β z 0 = β r . We set

β r = [β 1 , β 2 , • • • , β l+1 ]
where β r = 10 -r+1 . When r = 1 we have β = 1. In the simulations, we set l = 5, hence we have r = [1, 2, 3, 4, 5, 6], and the hyperparameters β z 0 (j) are initialized with 6 ranks : β r = [10 0 , 10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 ]. The noise depends on the SNR of the dataset. The biased dataset is expressed as the sum of uncontaminated dataset and the additive noise :

g = g 0 + .
(5.80)

As the noise and the uncontaminated data g 0 are supposed to be independent, we have :

g 2 2 = g 0 2 2 + 2 2 .
(5.81)

The SNR of the dataset is :

SNR = 10 log g 0 2 2 2 2 = 10 log g 2 2 - 2 2 2 2 .
(5.82)

As the expectation of is

E[ ] = 0, v = E 2 ≈ 2 2 M = g 2 2 M × 1 1 + 10 SN R/10 .
(5.83)

Belonging to an Inverse Gamma distribution, the expectation of the variance E[v i |α 0 , β 0 ] = β 0 α 0 -1 , therefore :

β 0 = g 2 2 M × 1 1 + 10 SN R/10 × (α 0 -1) .
(5.84)

The two hyperparameters α 0 and β 0 are combined according to Eq.(5.84) and hence initialization for one of them is sufficient. Figure 5.9 shows the influence of the value of α 0 on the reconstruction. According to the results, a bigger value for α 0 results to a smaller value on RMSE, for different number of projections and SNR of dataset. For the case with a low SNR and less projections numbers, there is a slight augment of RMSE when the hyperparameter increase, however, it soon arrives at a platform and stay robust. This monotonous property solves the problem of initialization of this hyperparameter, as a big value for α 0 satisfies all cases, since when α 0 is greater than a certain value, the RMSE does not change by using different initialization value α 0 .

-Analysis of hyper-parameters and initialization

In the reconstruction of a real data, we do not know the SNR of the data set. But we can estimate the noise by some experiments. One could estimate the detected data of in two cases : the first is when there is no source of X-ray, and the second is to emit X-rays without any object. By using these two data sets, we could obtain a rough additive noise in order to get the information of v .

Fixing α ξ 0 and β ξ 0 For ξ, we do not have evident information on it, therefore the initialization for both α ξ 0 and β ξ 0 is necessary. Figure 5.10 and Figure 5.11 show respectively the influence of the initialization of α ξ 0 and β ξ 0 on the reconstruction. For α ξ 0 , Figure 5.10, the influence is monotone and the reconstructed phantom has a lower RMSE when α ξ 0 is smaller.

Figure 5.11 shows the influence of the hyperparameter β ξ 0 on RMSE of reconstructed phantom. From the curves we find that the minimum RMSE value appears when β ξ 0 is approaching 0.01. Consequently we use β ξ 0 = 0.01. There is an augment of the RMSE value when β ξ 0 increases, but it soon arrives at a platform which is slightly bigger than the minimum RMSE value. So we conclude that the value of β ξ 0 is relatively insensitive to the reconstruction quality.

-Simulation results with limited number of projections

In [START_REF] Dumitru | A Bayesian approach for periodic components estimation for chronobiological signals[END_REF], it is pointed out that when α and β of the St g distribution are both big, it approaches a Gaussian distribution, which is the case for the additive noise ; if α and β are small (approaching to 0), it is an non-informative distribution (Jeffreys distribution) ; while α and β are both small but not zero, then the St g has the sparsity enforcing property, which is the case for the additive noise ξ.

Simulation results with limited number of projections

We apply respectively 180, 90, 60, 45, 36 and 18 projections evenly distributed in [0 • , 180 • ] for reconstruction of the 3D Shepp Logan phantom of size 256 3 , each projection contains 256 × 256 detectors. The number of projections are chosen such that there are respectively one projection every 1 • , 2 • , 3 • , 4 • , 5 • and 10 • . The iterative methods have all been initialized with the FBP reconstruction which has a high frequency artifacts when limited number of projections are used. In Table 5.1, different metrics of evaluation of the reconstructed 256 3 Shepp Logan phantom are compared. It is shown that the HHBM method performs not always better than the TV method, especially when there are sufficient number of projections. But when there are insufficient detected data, for example in the case of 18 projections, the HHBM method stay more robust than TV. On the other hand, as we all know that the choice of regularization parameter plays an important role in the regularization methods like QR or TV, and the value for the regularization parameter should be selected for each different case, the HHBM method solves this problem. As we can see from Figure 5.9, Figure 5.10 and Figure 5.11, once we have chosen the hyperparameters in a certain interval which is not difficult to fix according to the theoretical base, we can get the appropriate reconstruction results. What is more important, in the Bayesian approach, the prior model can be chosen by many other suitable distributions, which gives more HHBM method possibilities for the models than the conventional regularization methods. We may also choose different point estimators from the posterior distribution, for example the Posterior Mean, etc. Figure 5.12 shows the reconstructed middle slice of "Shepp Logan" phantom and "Head" object by using TV and HHBM methods with 36 projections and SNR=40dB. The red curve illustrate the profile of the blue line position. In the reconstructed Shepp Logan phantom by using TV method, the three small circles on the top of the slice are not evident. By using the HHBM method, we can distinguish these three small circles. By comparing the profiles of the slice of reconstructed Shepp Logan phantom, we can see that by using HHBM method, the contour positions on the profile is closer to the original profile than the TV method. In the reconstructed Head object, there are more details than the simulated Shepp Logan phantom, especially in the zoom area in the second line in Figure 5.12. By comparing the results, we can see that for the type of object which contains some small details, the TV method derives a results with smoother homogeneous areas, but less details at the contour areas than the HHBM method. For the Head object, we can also see that HHBM gives more details than TV, and is closer to the original one. Some of the white material which is dispersed into discontinuous small blocks in Head object are connected in result of TV method. From these images we conclude that in the case of insufficient number of projections the proposed method gives results with clearer contours and details.

Figure 5.13 shows the reconstructed Shepp Logan phantom from 18 projections and SNR=40dB. In this very underdetermined case, the HHBM method can still get a result which is clear enough to distinguish the primary zones and contours of the object.

Figure 5.14 and Figure 5.15 show the comparison between the QR, the TV and the HHBM method with a high SNR=40dB and a low SNR=20dB dataset respectively. The abscissa corresponds to the number of projections evenly distributed from 0 • to 180 • , and the ordinate is the RMSE after 50 iterations. When the number of projections becomes less and less sufficient, the HHBM method outperforms the other two conventional methods.

Simulation results with limited angle of projections

In medical X-ray CT, a commonly faced problem is the limit of projection angles. In this part of the simulation, we treat the dataset of projection evenly distributed in a limited range of angles for the simulated 3D "Shepp Logan" phantom and the 3D "Head" object, both of which have a size of 256 3 . Figure 5.16 shows the middle slice of the reconstructed Shepp Logan phantom and the Head object, from 90 projections distributed between 0 • and 90 • . By using TV method, the reconstructed object is blurry along the the diagonal direction where there is no projection data estimated, and there is an obvious square corner while the object should have a round shape. By using the HHBM method, we get results with closer shape and clearer contours.

• , 45 • ], [0 • , 90 • ], [0 • , 135 • ] and [0 • , 180 • ].
Thus, there is one projection every 1 • . From these two figures, we conclude that the proposed HHBM method stay more robust than the other two conventional methods when there is limited angle of projections.

-Simulation results with limited angle of projections

Original

TV HHBM FIGURE 5.12 -Slice of reconstructed "Shepp Logan" and "Head" phantom of size 256 3 , with dataset of 36 projections and SNR=40 dB, by using TV and HHBM methods respectively. Bottom figures are part of the corresponding top figures.

HHBM method

Original TV HHBM FIGURE 5.13 -Slice of reconstructed Shepp Logan phantom of size 256 3 , with dataset of 18 projections and SNR=40 dB, by using TV (left) and HHBM (right) methods respectively. The red curves are the profile at the position of the corresponding blue lines.

Dual-tree complex wavelet transform (DT-CWT) Based Hierarchical Bayesian Method

In the proposed method, the multilevel Haar transformation coefficient of the phantom is considered as a hidden variable and is estimated during iterations. The reason of choosing Haar transformation is that the phantom under consideration is normally discrete and strictly piecewise continuous. Haar transformation is one of the simpliest and basic transformation among all the orthonormal transformations. Its computation can be efficiently accelerated thanks to a parallelization on GPUs. What's more, it is orthonormal and therefore it simplifies the computation in the Bayesian methods. Recently, there are more and more new transforms proposed with very good properties, for example the curvelet [CDDY06, SCD02], the contourlet [START_REF] Minh | The contourlet transform : an efficient directional multiresolution image representation[END_REF], etc.

When the wavelet transformations are used, some disadvantages should be considered [START_REF] Ivan W Selesnick | The dual-tree complex wavelet transform[END_REF] :

• Oscillations. In the wavelet transformation, the coefficients tend to oscillate positive and negative around the singularities. These oscillations around the singularities complicate the wavelet-based processing, and makes the singularity extraction and signal modeling more complicated. • Shift variance. In the wavelet transformation, a small shift of the signal will perturb the wavelet coefficient oscillation pattern around the singularities because of the positivenegative oscillation property of the coefficients. This property also complicates the wavelet domain processing. • Lack of directionality. While Fourier sinusoids in higher dimensions correspond to highly directional plane waves, the standard tensor product construction of M-D wavelets produces a checkerboard pattern that is simultaneously oriented along only several directions. This lack of directional selectivity complicates the modeling and processing

Dual-tree complex wavelet transform (DT-CWT) Based Hierarchical Bayesian Method

of geometric image features, for example the ridges and edges. To overcome these problems, the Dual Tree Complex Wavelet Transform (DT-CWT), first introduced by Kingsbury in 1998 [START_REF] Nick | The dual-tree complex wavelet transform : a new technique for shift invariance and directional filters[END_REF], is proposed. The dual-tree structure of wavelet filters are used to obtain the real and imaginary parts of complex wavelet coefficients. Figure 5.19 shows the filter bank analysis of the DT-CWT. The DT-CWT employs two real DWTs. The first DWT gives the real part of the transform while the second real DWT gives the imaginary part. The two real DWTs use two different sets of filters. h 0 (n), h 1 (n) denote the low-pass/high-pass filter pair for the upper filter bank, and g 0 (n), g 1 (n) denote the low-pass/high-pass filter pair for the lower filter bank. These two filter banks are realized by a typical DWT with two pairs of adapted filters. The inverse of DT-CWT is also a synthesis of two real DWTs. The real part and the imaginary part are each inverted to obtain two real signals. These two real signals are then averaged to obtain the final signal. If the two real DWTs are orthonormal transforms, then so does the DT-CWT. It means that the inverse of the DT-CWT can be performed using its transpose, so that the operator D is orthogonal.

- - g 1 (n) -↓ 2 - g 0 (n) -↓ 2 - - - h 1 (n) -↓ 2 - h 0 (n) -↓ 2 - - - g 1 (n) -↓ 2 - g 0 (n) -↓ 2 - - - h 1 (n) -↓ 2 - h 0 (n) -↓ 2 - - - g 1 (n) -↓ 2 - g 0 (n) -↓ 2 - - - h 1 (n) -↓ 2 - h 0 (n) -↓ 2 - FIGURE 5
.19 -Analysis Filter Bank for the dual-tree CWT with a different set of filters at each stage.

Figure 5.20 illustrates the corresponding synthesis filter bank of the DT-CWT transformation. Also, it was pointed out in [XYM + 12] that the dictionary redundancy improves the sparsity of representation. In the discrete Haar transform, the size of the coefficient is the same as the size of image. In DT-CWT, the size of coefficients is two times the size of image in 2D cases (see Figure 5.21) and 4 times the size of phantom in 3D cases. According to the definition, dual-tree of wavelet filters is used to obtain the real and imaginary parts of complex wavelet coefficients. So the coefficients can be divided into 2 groups, which correspond to the real and imaginary parts, respectively.

In this section, I am going to show the comparison of the DT-CWT with the HT. Then DT-CWT will be used in the proposed hierarchical Bayesian method and compared with the HHBM method.

-Simulation results with limited angle of projections

g1 (n) - ↑ 2 - g0 (n) - ↑ 2 - - h1 (n) - ↑ 2 - h0 (n) - ↑ 2 - - g1 (n) - ↑ 2 g0 (n) - ↑ 2 - - h1 (n) - ↑ 2 h0 (n) - ↑ 2 - - g1 (n) - ↑ 2 g0 (n) - ↑ 2 - h1 (n) - ↑ 2 h0 (n) - ↑ 2 - - 0.5 FIGURE 5
.20 -Synthesis Filter Bank for the dual-tree CWT with a different set of filters at each stage. Dual-tree complex wavelet transform (DT-CWT) Based Hierarchical Bayesian Method

Adding noise to coefficients z

When the coefficients are treated as a variable in the algorithms, there will be appearence of noise in the coefficient. So in the second simulation, I am going to add noise to the coefficient with different SNR values (30dB, 20dB, 10dB and 5dB) and compare the reconstructed image f and the RMSE of them. The comparison is shown in Figure 5.

-Test 1 : Truncation in the transform coefficients z

From the comparison, we can see that when there is a noise on the coefficient, the inverse transformed image will be biased. By using the DT-CWT, the reconstructed image is better than the HT with a small RMSE value. We can conclude from the comparison that DT-CWT is more robust than HT against the noise added to the transformation coefficients.

Test on different images

The Shepp Logan phantom is artificially created and is strictly piecewise constant. In this case, the Haar transformation is appropriate. The same comparisons have been done by using some real images of different types, for example the image with textures, the piecewisecontinuous image and the image with blocs and textures. The images used are the "Period", "Stones", "Ground" and "Stairs", shown in Figure 5.24.

Test on different images

DT-CWT transformation. The 80%, 90%, 95% and 99% smallest part of the coefficients of Haar transform and the DT-CWT transform are removed respectively. The results demonstrate that when the coefficient is truncated in the algorithms, the DT-CWT transformation has a better performance for the reconstruction from the transform coefficients. 0.0008 4.1652e -5 0.0029 3.0084e -5 0.0017 4.5542e -5 0.0109 1.3460e -4 90% 0.0054 6.5086e -5 0.0094 4.8588e -5 0.0046 6.9940e -5 0.0258 2.0881e -4 95% 0.0242 8.0747e -5 0.0205 6.1184e -5 0.0089 8.5687e -5 0.0458 2.5838e -4 99% 0.1127 9.5165e -5 0.0658 7.2906e -5 0.0231 1.0029e -4 0.1011 3.0478e -4 5.6.2 Test 2 : Adding noise to the transform coefficients z Table 5.3 shows the RMSE of the images obtained by inverse transformation from the biased coefficient of the Haar and the DT-CWT transformation. A Gaussian noise of 5dB, 10dB, 15dB, 20dB and 30dB are added respectively. The same conclusion is drawn : by using the DT-CWT transformation, the reconstructed images are more robust against the noise added to the transform coefficient.

The use of DT-CWT in Bayesian method

In the previously proposed Bayesian method, the Haar transformation is used and its transformation coefficient is considered as a hidden variable. By comparing the DT-CWT with HT, we demonstrated some advantages of the DT-CWT. In this section, we replace the HT coefficient in the HHBM method by the DT-CWT coefficient, and compare the new Bayesian method with the HHBM method.

The system model depends always on the following two forward models :

g = Hf + , (5.85) f = Dz + ξ, (5.86)
where D is the DT-CWT and z the DT-CWT coefficient in this case.

-Conclusions and perspectives

Simulation results

In the HHBM method, the discrete multilevel Haar transformation coefficient is considered as the hidden variable. In this section, we replace the Haar transformation with the DT-CWT transformation. The transformation coefficient, z, is the DT-CWT coefficient of the object f . The generalized Student-t distribution is always used to model the coefficient z to enforce its sparse structure. The method of reconstruction by considering the DT-CWT transformation coefficient is named the BH-DTCWT method.

In the simulation, the 2D Shepp Logan phantom of size 256 2 is considered as the original object. The projections are evenly distributed from 0 • to 180 • . In Figure 5 Figure 5.26 shows the comparison of the middle slice of the reconstructed Shepp Logan phantom by using different reconstruction methods. From the comparison, we can see that, by using the hierarchical structured Bayesian method, the reconstructed phantom is clearer than the TV method when the number of projections is insufficient. Visually, there is no big difference between the results of the HHBM method and the BH-DTCWT method.

Figure 5.27 shows the comparison of a segment of profile around the contour in the Shepp Logan phantom. From the comparison, we see that the BH-DTCWT outperforms the TV method at the contour areas, as it is closer to the original profile.

In order to compare the HHBM and the BH-DTCWT methods, we show in Figure 5.28, Figure 5.29 and Figure 5.30 the evolution of the RMSE during reconstruction of different methods. The comparison shows that the convergence of BH-DTCWT is faster than the HHBM method, and this advantage reduces when there is less number of projections. While comparing with the TV method, the advantage of the HHBM and BH-DTCWT methods is obvious when the projection number is limited.

Conclusions and perspectives

By comparing the Haar transformation and the DT-CWT transformation, we see that the later is more robust when the coefficient is truncated or contaminated. By using DT-CWT trans- formation in the Bayesian method, the reconstructed phantom has clearer contours than the TV method. Visually, the reconstructed phantom by using the Bayesian method with Haar transformation or DT-CWT transformation does not have much difference, but the BH-DTCWT converges faster than HHBM.

In future work we will compare some other wavelet transformations which consider the property of the geometric figures and the direction of the contours, for example the curvelet transformation [CDDY06, SCD02], the contourlet transformation [START_REF] Minh | The contourlet transform : an efficient directional multiresolution image representation[END_REF], etc. Meanwhile, the Bayesian method using DT-CWT transformation is being implemented in 3D simulations. The computation of the DT-CWT transformation is expensive, particularly in 3D implementation. So we will implement the transformation by using GPU processor in order to accelerate the computation.

Variable Splitting Algorithm for HHBM

From the regularization point of view, variable splitting refers to the process of introducing auxiliary constraint variables to separate coupled components in the cost function. This procedure transforms the original minimization problem into an equivalent constrained optimization problem that can be effectively solved using classical constrained optimization schemes. The variable splitting approach has become popular recently for solving reconstruction problems in image processing [ABDF10, RF12, NVL + 13].

-VS in Bayesian inference

Variable Splitting (VS)

Consider an unconstrained optimization problem in which the objective function is the sum of two functions, one of which is written as the composition of two functions :

min f R 1 (f ) + R 2 (Φ(f )), (5.87)
where Φ is a linear function. Variable splitting is a very simple procedure that consists in creating a new variable, z, to serve as the argument of R 2 , under the constraint that Φ(f ) = z. This leads to the constrained problem :

min f ,z R 1 (f ) + R 2 (z),
subject to z = Φ(f ).

(5.88)

The convenience of the variable splitting methods is that it may be easier to solve the constrained problem in Eq.(5.88) than that of Eq.(5.87). The splitting variable has been recently used in many inverse problem methods, for example in [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF].

The constrained problem in Eq.(5.88) can be derived into a quadratic penalty approach, by solving :

min f ,z R 1 (f ) + R 2 (z) + λ 2 Φ(f ) -z 2 2
(5.89) by alternately minimizing f and z respectively.

While splitting the variable f , we get the criterion similar to the form of the MAP optimization criterion of the HHBM method proposed previously, with one penalty term on f , one on z and one penalize the disparity between them.

VS in Bayesian inference

In the regularization methods, the splitting of variables normally leads to the splitting of regularization terms. In the Bayesian inference, on the other hand, it is considered as a splitting of the variables. By using MAP optimization, we can see the links between the regularization methods and a Bayesian inference. However, the Bayesian approach provide the possibilities to estimate the hyperparameters of the inverse problem which are often crucial in real applications.

In X-ray CT models, typically only one variable is used to represent all the additive noise. Consequently, the system noise and the uncertainties of the forward model are all modeled by the same prior distribution, for example the Gaussian distribution or the Student-t distribution that we used previously. In our work, we propose a noise splitting model. This model considers the system noise in different aspects and use different distributions for the modeling of prior variables.

In a real projection system of X-ray CT, the acquisition of the sinogram dataset is obtained mainly by two steps : the projection and the detection. The projection step is the attenuation of the intensity of radiations when passing through the object, and the detection step is the acquisition of the sinogram by the detectors.

Supposing that the true projection data before detection is g 0 , and the detected projection data is g, and theoretically we have g = g 0 if the detection noise is not considered. The noise of Variable Splitting Algorithm for HHBM detection in tomography is typically modeled by a Poison distribution. In X-ray CT, the number of the photon is big and the Poison distribution is approximated to a Gaussian distribution [START_REF] Sauer | A local update strategy for iterative reconstruction from projections[END_REF]. So the detection of the attenuated radiation leads to an additive noise :

g = g 0 + , (5.90)
where is a Gaussian additive noise, with mean equal to zero and variance v .

Before the detection of the radiation, these rays are attenuated by the object. The exact attenuated data g 0 is obtained according to the Beer's law, Eq.(3.3), and the linear system is model by the direct model g 0 = Hf . Nevertheless, the physical phenomenons during the attenuation of X radiation exist, for example the scattering of some photons, and sometimes the photons with low energies would be absorbed by some materials, for example the metals. These kinds of physical phenomenons will cause uncertainties for the linear model g 0 = Hf , and hence an additive noise should be added for the direct model :

g 0 = Hf + ρ.
(5.91)

The uncertainties caused by the scattering or the absorption appears only for few of the photons among a large amount of photons, and hence ρ is sparse. It is modeled by the generalized Student-t distribution in our work. Notice that in our work we consider only the additive noise model.

By splitting the noise, the following model is considered :

g = Hf + ρ + , (5.92)
where represents the measurement noise and ρ can represent the modelling errors. It can be expressed as :

g = g 0 + , g 0 = Hf + ρ, (5.93)
where g 0 is the hidden variable of the system. Therefore, the additive noise stands for the noise of the detectors, and ρ for the sparse noise.

With this model, by using the classical regularization methods, the solution is defined as the optimizer of :

J(f , g 0 ) = g -g 0 2 2 + λ 1 g 0 -Hf 2 2 + λ 2 R(f ), (5.94)
where parameter λ 1 and λ 2 are chosen beforehand, manually or using an ad-hoc algorithm. An alternate optimization algorithm is :

f (k+1) = f (k) + α 1 2λ 1 H T g 0 -Hf (k) -λ 2 ∇R(f (k) ) , g (k+1) 0 = g (k)
0 + α 2 λ 1 g -Hf (k) .

(5.95)

We can compare it with the ADMM method presented in Chapter 3 where g 0 plays the role of µ.

5.8.3 -Gaussian-St g prior model for the variable splitting method (HHBM-GS)

Gaussian-St g prior model for the variable splitting method (HHBM-GS)

In Bayesian inference, prior models are defined according to the prior information. The noise of detector, , is modeled by a Gaussian distribution, with a fixed variance v , wile the uncertainties ρ is modeled by a heavy tailed distribution, enforcing its sparse structure.

α z 0 , β z 0 α ξ 0 , β ξ 0 α ρ0 , β ρ0 v z z f g 0 g v v ρ v ξ D H p(g|g 0 , v ) p(g 0 |f , v ρ ) p(f |z, v ξ ) p(z|v z ) p(v z |α z 0 , β z 0 ) p(v ξ |α ξ 0 , β ξ 0 ) p(v ρ |α ρ0 , β ρ0 ) FIGURE 5.31 -Generative graph of the Gaussian-St g Variable Splitting method.
The statistical model is expressed via the forward model, Eq.(5.93), and the generative graph of this system model is illustrated in Figure 5.31.

Variable Splitting Algorithm for HHBM

p(g|g 0 , v ) = N (g|g 0 , V ) ∝ exp - 1 2 (g -g 0 ) T V -1 (g -g 0 ) , p(g 0 |f , v ρ ) = N (g 0 |Hf , V ρ ) ∝ |V ρ | -1 2 exp - 1 2 (g 0 -Hf ) T V -1 ρ (g 0 -Hf ) , p(f |z, v ξ ) = N (f |Dz, V ξ ) ∝ |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) , p(z|v z ) = N (z|0, V z ) ∝ |V z | -1 2 exp - 1 2 z T V -1 z z where V z = diag [v z ] , p(v ρ |α ρ0 , β ρ0 ) = IG(v ρ |α ρ0 , β ρ0 ) = M i=1 IG(v ρ i |α ρ0 , β ρ0 ) ∝ M i=1 v -(α ρ0 +1) ρ i exp -β ρ0 v -1 ρ i , p(v ξ |α ξ 0 , β ξ 0 ) = IG(v ξ |α ξ 0 , β ξ 0 ) = N j=1 IG(v ξ j |α ξ 0 , β ξ 0 ) ∝ N j=1 v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j , p(v z |α z 0 , β z 0 ) = IG(v z |α z 0 , β z 0 ) = N j=1 IG(v z j |α z 0 , β z 0 ) ∝ N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j .
(5.96)

As previous, the posterior distribution is obtained from these models via the Bayes rule :

p(g 0 ,f , z, v ρ , v ξ , v z |g) ∝p(g|g 0 , v )p(g 0 |f , v ρ )p(f |z, v ξ )p(z|v z )p(v ρ |α ρ0 , β ρ0 )p(v ξ |α ξ 0 , β ξ 0 )p(v z |α z 0 , β z 0 ) = exp - 1 2 (g -g 0 ) T V -1 (g -g 0 ) |V ρ | -1 2 exp - 1 2 (g 0 -Hf ) T V -1 ρ (g 0 -Hf ) |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) |V z | -1 2 exp - 1 2 z T V -1 z z M i=0 v -(α ρ0 +1) ρ i M i=0 exp -β ρ0 v -1 ρ i N j=0 v -(α ξ 0 +1) ξ j N j=0 exp -β ξ 0 v -1 ξ j N j=0 v -(αz 0 +1) z j N j=0 exp -β z 0 v -1 z j .
(5.97)

The negative logarithm of the posterior distribution is calculated to simplify the optimiza-5.8.3 -Gaussian-St g prior model for the variable splitting method (HHBM-GS) tion :

J(g 0 ,f , z, v ρ , v ξ , v z ) = -ln p(g 0 , f , z, v ρ , v ξ , v z |g) = 1 2 (g -g 0 ) T V -1 (g -g 0 ) + 1 2 M i=1 ln v ρ i + 1 2 (g 0 -Hf ) T V -1 ρ (g 0 -Hf ) + 1 2 N j=1 ln v ξ j + 1 2 (f -Dz) T V -1 ξ (f -Dz) + 1 2 N j=1 ln v z j + 1 2 z T V -1 z z + (α ρ0 + 1) M i=1 ln v ρ i + β ρ0 M i=1 v -1 ρ i + (α ξ 0 + 1) N j=1 ln v ξ j + β ξ 0 N j=1 v -1 ξ j + (α z 0 + 1) N j=1 ln v z j + β z 0 N j=1 v -1 z j .
(5.98)

Then the unknowns are optimized alternately by minimizing the criterions :

g 0 = arg min g 0 1 2 (g -g 0 ) T V -1 (g -g 0 ) + 1 2 (g 0 -Hf ) T V -1 ρ (g 0 -Hf ) , (5.99)
f = arg min f 1 2 (g 0 -Hf ) T V -1 ρ (g 0 -Hf ) + 1 2 (f -Dz) T V -1 ξ (f -Dz) , (5.100)
z = arg min z 1 2 (f -Dz) T V -1 ξ (f -Dz) + 1 2 z T V -1 z z , (5.101)
v ρ i = arg min vρ i 1 2 ln v ρ i + 1 2 [g 0 -Hf ] 2 i v -1 ρ i + (α ρ0 + 1) ln v ρ i + β ρ0 v -1 ρ i , (5.102) v ξ j = arg min v ξ j 1 2 ln v ξ j + 1 2 [f -Dz] 2 j v -1 ξ j + (α ξ 0 + 1) v -1 ξ j + β ξ 0 v -1 ξ j , (5.103) v z j = arg min vz j 1 2 ln v z j + 1 2 [z] 2 j v -1 z j + (α z 0 + 1) ln v z j + β z 0 v -1 z j .
(5.104)

For the optimization of g 0 , f and z, the analytical expressions are :

g 0 = V -1 + V -1 ρ -1 V -1 g + V -1 ρ Hf , (5.105) f = H T V -1 ρ H + V -1 ξ -1 H T V -1 ρ g 0 + V -1 ξ Dz , (5.106) z = D T V -1 z D + V -1 z -1 D T V -1 z f .
(5.107)

In these three expressions, the analytical expression for g 0 is realizable because the matrix V -1 + V -1 ρ is a diagonal matrix. The inversion of a diagonal matrix is another diagonal matrix with each element equal to inversion of corresponding element of the original matrix. Hence :

g 0 = diag 1 v -1 i + v -1 ρ i V -1 g + V -1 ρ Hf (5.108)
But for f and z, the terms of matrix inversion is not realizable in the background of big data size images reconstruction. So as we have presented previously, the gradient descent optimization algorithm is used.

Variable Splitting Algorithm for HHBM

The analytical expression for optimization of v ρ i , v ξ j and v z j are also obtained by calculating the critical point of the corresponding criterion.

Consequently, the iterative optimization expressions of all the unknown variables are given bellow :

g 0 = V -1 + V -1 ρ -1 V -1 g + V -1 ρ Hf , for i = 1 : I 2 f (k+1) = f (k) -γ (k) f ∇J( f (k) ), for i = 1 : I 2 z (k+1) = z (k) -γ (k) z ∇J( z (k) ), v ρ i = β ρ0 + 1 2 [g 0 -Hf ] 2 i α ρ0 + 3 2 , v ξ j = β ξ 0 + 1 2 [f -Dz] 2 j α z 0 + 3 2 , v z j = β z 0 + 1 2 z 2 j α z 0 + 3 2 , ( 5.109)
where

∇J( f (k) ) = -H T V -1 ρ (g 0 -Hf ) + V -1 ξ (f -Dz) , (5.110)
∇J( z (k) ) = -D T V -1 z (f -Dz) + V -1 z z, (5.111)
γ (k) f = ∇J(f ) 2 2 Y ρ H∇J(f ) 2 2 + Y ξ ∇J(f ) 2 2
, where

Y ρ = V -1 2 ρ and Y ξ = V -1 2 ξ ,
(5.112)

γ (k) z = ∇J(z) 2 2 Y ξ D∇J(z) 2 2 + Y z ∇J(z) 2 2
, where

Y z = V -1 2 z .
(5.113)

As we have demonstrated, the Normal distribution can be approximated by a St g distribution. So in the variable splitting method, we also consider to use the St g distribution to model the noise , presented in the next section.

St g -St g prior model for the variable splitting method (HHBM-SS)

Another prior model for the variable splitting method is to model both the noise ρ and by a generalized Student-t distribution. As we have mentioned previously, the St g distribution can be either heavy-tailed for modeling a sparse variable or approach to a Gaussian distribution for non-sparse variable by choosing different pair of hyper-parameters α and β.

In this model, both two additive noises, and ρ are model by a generalized Student-t distribution. Comparing with the previous model, is modeled with a Gaussian distribution with an unknown variance v belonging to an Inverse Gamma distribution :

p( |v ) = N ( |0, V ), where V = diag [v ] , p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ).
α z 0 , β z 0 α 0 , β 0 α ξ 0 , β ξ 0 α ρ0 , β ρ0 v z z f g 0 g v v ρ v ξ D H p(g|g 0 , v ) p(g 0 |f , v ρ ) p(f |z, v ξ ) p(z|v z ) p(v z |α z 0 , β z 0 ) p(v ξ |α ξ 0 , β ξ 0 ) p(v ρ |α ρ0 , β ρ0 ) p(v |α 0 , β 0 ) FIGURE 5.32 -Generative graph of the St g -St g Variable Splitting method.
Consequently the posterior distribution is :

p(g 0 ,f , z, v , v ρ , v ξ , v z |g) ∝p(g|g 0 , v )p(g 0 |f , v ρ )p(f |z, v ξ )p(z|v z )p(v |α 0 , β 0 )p(v ρ |α ρ0 , β ρ0 )p(v ξ |α ξ 0 , β ξ 0 )p(v z |α z 0 , β z 0 ) = |V | -1 2 exp - 1 2 (g -g 0 ) T V -1 (g -g 0 ) |V ρ | -1 2 exp - 1 2 (g 0 -Hf ) T V -1 ρ (g 0 -Hf ) |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) |V z | -1 2 exp - 1 2 z T V -1 z z M i=0 v -(α 0 +1) i M i=0 exp -β 0 v -1 i M i=0 v -(α ρ0 +1) ρ i M i=0 exp -β ρ0 v -1 ρ i N j=0 v -(α ξ 0 +1) ξ j N j=0 exp -β ξ 0 v -1 ξ j N j=0 v -(αz 0 +1) z j N j=0 exp -β z 0 v -1 z j ,
(5.114) 5.8.5 -Initialization of hyperparameters for HHBM-GS and HHBM-SS methods gradient algorithm. The iterative optimization expressions of all the unknown variables are :

g 0 = V -1 + V -1 ρ -1 V -1 g + V -1 ρ Hf , for i = 1 : I 2 f (k+1) = f (k) -γ (k) f ∇J( f (k) ), for i = 1 : I 2 z (k+1) = z (k) -γ (k) z ∇J( z (k) ), v i = β 0 + 1 2 [g -g 0 ] 2 i α 0 + 3 2 , v ρ i = β ρ0 + 1 2 [g 0 -Hf ] 2 i α ρ0 + 3 2 , v ξ j = β ξ 0 + 1 2 [f -Dz] 2 j α z 0 + 3 2 , v z j = β z 0 + 1 2 z 2 j α z 0 + 3 2 , ( 5.123)
where

∇J( f (k) ) = -H T V -1 ρ (g 0 -Hf ) + V -1 ξ (f -Dz) , (5.124)
∇J( z (k) ) = -D T V -1 z (f -Dz) + V -1 z z, (5.125) γ (k)
f = ∇J(f ) 2 2 Y ρ H∇J(f ) 2 2 + Y ξ ∇J(f ) 2 2 , (5.126)
γ (k) z = ∇J(z) 2 2 Y ξ D∇J(z) 2 2 + Y z ∇J(z) 2 2 .
(5.127)

Initialization of hyperparameters for HHBM-GS and HHBM-SS methods

The hyperparameters, α 0 , β 0 , α ρ0 and β ρ0 , of the St g distribution are initialized by considering the prior property of the two noises and ρ. According to the system model of X-ray CT, the hyperparameters are initialized such that is sparse and ρ approach to a Gaussian noise.

As presented in Section 5.4.2, the initialization of the hyperparameters considers the SNR of the dataset. The variance of noise v n can be estimated beforehand with the projection data without any object. According to Eq.(5.83), we have :

v n ≈ g 2 2 M × 1 1 + 10 SN R/10 .
(5.128)

In our work, without any further information about the two different noise of system and ρ, we suppose that v = v ρ , and each of them equal to half of the total noise variance :

v = v ρ = 1 2 v n .
(5.129)

• For HHBM-GS method, we fix v = v n , and the relation β ρ0 = v ρ (α ρ0 -1) is considered. We fix α ρ0 = 2.01 in order to have a sparse structure for ρ.

6

Simultaneous Object-Contour Reconstruction Model

Introduction

In the previous chapter, a hierarchical Bayesian method which uses a heavy tailed generalized Student-t distribution to enforce the sparse structure of the discrete wavelet transformation coefficient is used. As it was pointed out, the advantage of the wavelet transformation comparing with the gradient operator, or the Laplacian operator, is that, while the inversion of gradient operator is inconvenient, the wavelet transformation operator is reversible and orthogonal.

However, the most accurate transformation to obtain a sparse structured coefficient of a piecewise continuous image is the gradient operator or the Laplacian operator, because in the wavelet transformation coefficient, the part corresponding to the low-pass filtered coefficients is not sparse.

In this chapter, we present a Bayesian method which combines the reconstruction of the object and the reconstruction of the Laplacian operator of object. For a piecewise continuous object, the Laplacian is sparse. The St g distribution is used as the prior model for the sparse variable. The prior distribution for the variable f is a non-homogeneous Markovian model. In this non-homogeneous Markovian model, the distribution of each voxel f i depends not only on its six neighbor voxels, but also the Laplacian value of the neighbor voxels. By taking into consideration the Laplacian values, we can distinguish whether the element is in a homogeneous area or on a contour, and the weight of impact of this neighbor voxel on the current voxel will depend on this Laplacian value.

Forward model : Radon Transform and notations

As presented in Chapter 3, the Radon transform is analytically expressed as :

g(r, φ) = Rf (x, y) = L r,φ f (x, y)dl, (6.1)
where R represents the operator which transforms f to Radon space, f (x, y) represents the image, g the detected projection data, r the perpendicular length from center point of coordinate and φ the considered X ray angle. L r,φ is the length of ray (r, φ) passing through the image.

We can also use a vector to replace φ to define the direction of the radiation. As shown in Figure 6.1, the direction of the radiation can be defined by the unit vector η, and r is always the vertical distance from the origin to the radiation. The RT can be represented by the following expression :

g(r, η) = +∞ -∞ f (x, y)δ (r -pη) dt (6.2)
where pη = xη 1 + yη 2 = x cos φ + y sin φ, so that :

g(r, η) +∞ -∞ +∞ -∞ f (x, y)δ(r -xη 1 -yη 2 )
dx dy (6.3)

Basic properties of Radon transformation

Lemma 6.3.1. The shifting property :

Rf (x -a, y -b) = g(r -aη 1 -bη 2 , η) (6.4) Demonstration. By using the definition given in Eq.(6.3), we have :

Rf (x -a, y -b) = +∞ -∞ +∞ -∞ f (x -a, y -b)δ (r -xη 1 -yη 2 ) dx dy (6.5)

Basic properties of Radon transformation

By changing the variable : x = x -a and y = y -b, we get :

Rf (x , y ) = +∞ -∞ +∞ -∞ f (x , y )δ ((r -aη 1 -bη 2 ) -x η 1 -y η 2 ) dx dy =g (r -aη 1 -bη 2 , η) (6.6)
Hence, the following properties are derived :

Rf (x -a, y) = g(r -aη 1 , η) Rf (x, y -b) = g(r -bη 2 , η) (6.7) Lemma 6.3.2. Derivation :

R ∂ 2 f ∂x 2 + ∂ 2 f ∂y 2 = ∂ 2 g ∂r 2 (6.8)
Demonstration. The definition of the second order derivative is :

∂ 2 f (x, y) ∂x 2 = lim ε→0 f (x + ε, y) + f (x -ε, y) -2f (x, y) ε 2 ∂ 2 f (x, y) ∂y 2 = lim ε→0 f (x, y + ε) + f (x, y -ε) -2f (x, y) ε 2 (6.9)
The Radon transform of the derivatives are respectively :

R ∂ 2 f (x, y) ∂x 2 = lim ε→0 Rf (x + ε, y) + Rf (x -ε, y) -2Rf (x, y) ε 2 R ∂ 2 f (x, y) ∂y 2 = lim ε→0 Rf (x, y + ε) + Rf (x, y -ε) -2Rf (x, y) ε 2 (6.10)
According to the shifting property, Lemma 6.3.1, we get :

R ∂ 2 f (x, y) ∂x 2 = lim ε→0 g (r + εη 1 , η) + g (r -εη 1 , η) -2g(r, η) ε 2 =η 2 1 lim (η 1 ε)→0 g (r + εη 1 , η) + g (r -εη 1 , η) -2g(r, η) (εη 1 ) 2 =η 2 1 ∂ 2 g(r, η) ∂r 2 (6.11) R ∂ 2 f (x, y) ∂y 2 = lim ε→0 g (r + εη 2 , η) + g (r -εη 2 , η) -2g(r, η) ε 2 =η 2 1 lim (η 2 ε)→0 g (r + εη 2 , η) + g (r -εη 2 , η) -2g(r, η) (εη 2 ) 2 =η 2 2 ∂ 2 g(r, η) ∂r 2
(6.12)

As η 2 1 + η 2 2 = 1, we obtain :

R ∂ 2 f (x, y) ∂x 2 + ∂ 2 f (x, y) ∂y 2 = η 2 1 + η 2 2 ∂ 2 g(r, η) ∂r 2
= ∂ 2 g(r, η) ∂r 2

(6.13)

CHAPITRE 6. SIMULTANEOUS OBJECT-CONTOUR RECONSTRUCTION MODEL

According to the previous demonstrated property, we obtain the following relationships between the objection f and projection g and between the Laplacian of object f and the Laplacian of projection g : f (x, y) →RT → g(r, φ), (6.14) f (x, y) →RT → g(r, φ), (6.15)

where

f (x, y) = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 f (x, y), g(r, φ) = ∂ 2 ∂r 2
g(r, φ). (6.16)

The same property can also be demonstrated in 3D projection.

Figure 6.2 shows, from left to right and top to bottom, the Shepp Logan phantom of size 256 2 , the projection of Shepp Logan phantom with 256 projections, the Laplacian of the Shepp Logan phantom and the Laplacian of the projection.

In the statistical methods, the discretized data are considered. In the discretized Radon transform system, the transformation operation is represented via the multiplier of a matrix H. So that : g = Hf + .

(6.17)

where an additive noise is presented by vector , representing the error of detection and uncertainties of the projection system.

Meanwhile, the Radon transform from Laplacian of object to the derivation of projection, shown in Eq.(6.15), is expressed as :

g = H f + τ (6.18)
where the additive error represents the uncertainties of the transform system and other errors of the Laplacian operation.

Noted by f the discretized f (x, y) and by g discretized g(r, φ), they are approximated by convolution of the corresponding variable with kernels : In the conventional image reconstruction methods, for example the TV method, the gradient of phantom f is used as a criterion to preserve the edges of image. In the proposed Reconstruction of Object and Contour Combined method (ROCC), the Laplacian of the object f is considered as an unknown variable and is estimated during the reconstruction of the object.

g = [-1 2 -1] * g = d 1 * g (6.19) f =   0 -1 0 -1 4 -1 0 -1 0   * f = d 2 * f (6.

Proposed Method of Reconstruction of Object while Considering Contours (ROCC)

In this proposed ROCC method, we propose to use the forward relations in Eq.(6.21) to simultaneously reconstruct the image f and detect the contours f by using a Bayesian approach. A hierarchical inhomogeneous Markovian model is defined for the object f . A hidden variables q, which is the normalized Laplacian, f -min( f ) / max( f ) -min( f ) , is used as a parameter of this prior model. A sparsity enforcing Generalized Student-t prior distribution is used to enforce the sparse structure of f .

Noise model

Starting from the forward model in Eq.(6.17) and Eq.(6.18), by assuming that the additive noise and τ are iid Gaussian distributed with variance v and v τ , we obtain the likelihood :

p(g|f , v ) = N (g|Hf , V ), p(g| f , v τ ) = N (g|H f , V τ ), (6.22)
where V = diag [v ] and V τ = diag [v τ ] are diagonal matrix. Considering that the variances v i and v τ i are positive and are not very big, they are defined by Inverse Gamma distributions :

p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ), p(v τ |α τ 0 , β τ 0 ) = IG(v τ |α τ 0 , β τ 0 ), (6.23)
where α 0 , β 0 , α τ 0 and β τ 0 are hyper parameters.

Another convenience of choosing the Inverse Gamma distribution for the variance is that, when the likelihood is Gaussian distributed, they are therefore conjugate distributions, and the form of the posterior distribution will still be Gaussian and Inverse Gamma distribution.

Sparsity enforcing Student-t prior model for f

As we have shown in Chapter 4, in X-ray CT, the object under consideration is commonly piecewise continuous. The Laplacian of this type of object is sparse, as shown in Figure 6.2, with most of the part approaching to zero except the contours.

As f is sparse, we use a sparsity enforcing prior to define it. As presented in the previous chapter, the generalized Student-t distribution is a suitable distribution and it has many appro-

Proposed Method of Reconstruction of Object while Considering Contours (ROCC)

priate properties in the Bayesian method. By considering its Normal-Inverse Gamma marginalization form, we model the hidden variable f by :

p( f |v c ) = N ( f |0, V c ), p(v c |α c 0 , β c 0 ) = IG(v c |α c 0 , β c 0 ), (6.24)
where

V c = diag [v c ],
and α c 0 and β c 0 are two hyper-parameters of generalized Student-t distribution.

The prior model for object

In the ROCC method, f is defined by an inhomogeneous Markovian distribution basing on the Laplacian f . In this Markovian model, the neighbor pixels are considered. For 2D image, four neighbor pixels are considered, and for 3D object, six neighbor voxels are considered, as shown in Figure 6.4.

f i,j f i-1,j f i,j+1 f i,j-1 f i+1,j f i,j,k f i+1,j,k f i,j+1,k f i,j,k-1 f i-1,j,k f i,j-1,k f i,j,k+1
FIGURE 6.4 -The neighbour pixels of current pixel in 2D image (left) and the neighbour voxels of current voxel in 3D object (right).

In the conventional homogeneous Markovian model, every neighbor pixel will have the same weight on the current pixel. It works well for the pixels in the homogeneous areas. As shown in Figure 6.5, when the pixel is around the contours between different homogeneous areas, the neighbor pixels should have different weight of influence on the current pixel. Here we use f r to represent the current pixel value, and f r to represent its neighbor pixels values. When f r is not on the contour, it means that f r is in the same homogeneous zone as f r , and that they are approach to each other. When f r is on the contour, on the other hand, f r should not have much influence on f r .

By considering these relations, we propose to use a normalized Laplacian value q :

q r = f r -min f r max f r -min f r . (6.25)
So q r is between 0 and 1, ∀r ∈ N r(r) and ∀r ∈ [1, N ], where N r(r) means the neighbor pixels of r. When q r = 0, it corresponds to a minimum Laplacian value, and it means that the neighbor pixels f r is in the same homogeneous areas as f r and f r should depend on the value 6.4.3 -The prior model for object f r . When q r = 1, it corresponds to a maximum Laplacian value, representing that the neighbor pixel f r is on the contour and should not have influence on f r . Consequently, the influence of the neighbors pixels in the Markovian model is inversely proportional to the value of q r .

We define the variable f by the following Gaussian distribution :

p(f |m f , v f ) = N (f |m f , V f ) (6.26)
where the variance v f are supposed belonging to an Inverse Gamma distribution, in order to realize a semi-supervised model. The mean of this Gaussian distribution, m f , depends on the neighbors of each pixel. The definition of the mean of the pixel r is :

m r = r ∈N r(r) (1 -q r ) f r r ∈N r(r) (1 -q r )
. (6.27)

For each neighbor pixel, r , the pixel value f r is considered with a weight (1 -q r ). q r is the normalized Laplacian value of the neighbor pixel r , defined in Eq.(6.25). Term

1 r ∈N r(r) (1-q r )
is added for the normalization.

Therefore, the expression of the prior distribution for each pixel f r is :

p(f r |f -r , q, v fr ) =N (f r |m r , v fr ) ∝v -1 fr exp    - 1 2 v -1 fr f r - r ∈N r(r) (1 -q r ) f r r ∈N r(r) (1 -q r ) 2    , (6.28)
where f -r represents all the elements in vector f except f r .

We show an example in Figure 6.5 as an explanation of this prior model.

In this example, f r is the value of the current pixel, q r is the normalized Laplacian of pixel r. r 1 , r 2 , r 3 and r 4 are the four neighbor pixel of r. Among them, r 1 and r 4 are on the contour with normalized Laplacian values q r equal to 1, while r 2 and r 3 are in the homogeneous areas with q r equal to 0. Thus, q r 1 = q r 4 = 1 and q r 2 = q r 3 = 0. As defined in Eq.(6.27), the mean of prior distribution for f r is :

m r = 1 -q r 1 f r 1 + 1 -q r 2 f r 2 + 1 -q r 3 f r 3 + 1 -q r 4 f r 4 1 -q r 1 + 1 -q r 2 + 1 -q r 3 + 1 -q r 4 = f r 2 + f r 3 2 .
(6.29)

So that f r belong to a Gaussian distribution :

p(f r ) = N f r | f r 2 + f r 3 2 , v fr . (6.30)
In a more complicated case, the value q r is distributed in [0, 1]. Consequently, a neighbor pixel with a bigger value of q r will have less influence on the current pixel.

Semi-supervised ROCC method f r ,q r f r' 4 ,q r' 4 f r'1 ,q r' 1 f r' 3 ,q r' 3 f r' 2 ,q r' 2 : q(r')=1 : q(r')=0 q=1 q=0 FIGURE 6.5 -Example of the inhomogeneous Markovian prior model of f

Semi-supervised ROCC method

In most of the cases, initializing the variance of phantom v f manually demand some prior information, which is typically not accessible in real applications. A semi-supervised model is considered in our work, in which v f is an unknown vector and the Inverse Gamma distribution is defined for it :

p(v f |α f 0 , β f 0 ) = IG(v f |α f 0 , β f 0 ) ∝ N j=1 v -(α f 0 +1) f j exp -β f 0 v -1 f j .
(6.31)

The system model

Basing on the previously presented prior model for the object f , and by considering v f and an unknown variable, belonging to an Inverse Gamma distribution, we get the following system 6.5.1 -The system model model :

p(g|f , v ) = N (g|Hf , V ) ∝ |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) , p(g| f , v τ ) = N (g|H f , V τ ) ∝ |V τ | -1 2 exp - 1 2 g -H f T V -1 τ g -H f , p(f |m f , v f ) = N (f |m f , V f ) ∝ |V f | -1 2 exp - 1 2 (f -m f ) T V -1 f (f -m f ) , p( f |v c ) = N ( f |0, V c ) ∝ |V c | -1 2 exp - 1 2 f T V -1 c f , p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ) ∝ M i=1 v -(α 0 +1) i M i=1 exp -β 0 v -1 i , p(v τ |α τ 0 , β τ 0 ) = IG(v τ |α τ 0 , β τ 0 ) ∝ M i=1 v -(ατ 0 +1) τ i M i=1 exp -β τ 0 v -1 τ i , p(v f |α f 0 , β f 0 ) = IG(v f |α f 0 , β f 0 ) ∝ N j=1 v -(α f 0 +1) f j N j=1 exp -β f 0 v -1 f j , p(v c |α c 0 , β c 0 ) = IG(v c |α c 0 , β c 0 ) ∝ N j=1 v -(αc 0 +1) c j N j=1 exp -β c 0 v -1 c j , (6.32)
where m f = r ∈Nr (1 -q r )f r r ∈Nr (1 -q r ) , (6.33) and

q = f -min f max f -min f . (6.34)
According to the Bayes rule, the posterior distribution is :

p(f , f , v , v τ , v f , v c |g, g) ∝p(g|f , v )p(g| f , v τ )p(f |m f , v f )p( f |v c )p(v |α 0 , β 0 ) p(v τ |α τ 0 , β τ 0 )p(v f |α f 0 , β f 0 )p(v c |α c 0 , β c 0 ) = |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) |V τ | -1 2 exp - 1 2 g -H f T V -1 τ g -H f |V f | -1 2 exp - 1 2 (f -m f ) T V -1 f (f -m f ) |V c | -1 2 exp - 1 2 f T V -1 c f M i=1 v -(α 0 +1) i M i=1 exp -β 0 v -1 i M i=1 v -(ατ 0 +1) τ i M i=1 exp -β τ 0 v -1 τ i N j=1 v -(α f 0 +1) f j N j=1 exp -β f 0 v -1 f j N j=1 v -(αc 0 +1) c j N j=1 exp -β c 0 v -1 c j .
(6.35)

-JMAP Estimation

for all the variables :

f = H T V -1 H + V -1 f -1 H T V -1 g + V -1 f m f , f = H T V -1 τ H + V -1 c -1 H T V -1 τ g, v i = β 0 + 1 2 [g -Hf ] 2 i α 0 + 3 2 , v τ i = β τ 0 + 1 2 g -H f 2 i α τ 0 + 3 2 , v f j = β f 0 + 1 2 [f -m f ] 2 j α f 0 + 3 2 , v cj = β c 0 + 1 2 f 2 j α c 0 + 3 2 . (6.39)
The updating of the variables v i , v f j and v cj are explicit, but as we have pointed out in Chapter 5 that, the inversion of the non-diagonal big size matrix

H T V -1 H + V -1 f and H T V -1 H + V -1 c
are not realizable, we use the gradient descent algorithm to optimize the unknown variables f and f :

f (k+1) = f (k) -γ f ∇J( f (k) ), f (k+1) = f (k) -γ f ∇J( f (k) ), (6.40)
where ∇J(•) represents the derivation of the corresponding criterion J(•), and γ • represents the descent step length. The descent step length is obtained by minimizing the following criterion :

γ f = arg min γ f J( f -γ f ∇J( f )) , γ f = arg min γ f J( f -γ f ∇J( f )) , (6.41)
and the expressions of γ f and γ f are respectively :

γ f = ∇J(f ) 2 Y H∇J(f ) 2 + Y f ∇J(f ) 2 , where Y = V -1 2 and Y f = V -1 2 f , (6.42) γ f = ∇J( f ) 2 Y H∇J( f ) 2 + Y c ∇J( f ) 2 , where Y c = V -1 2 c .
(6.43)

The optimization of all the variables by using ROCC method follows the following iterative

Experiment results

updating rules :

iter : f (k+1) = f (k) -γ f ∇J( f (k)
), (6.44)

iter : f (k+1) = f (k) -γ f ∇J( f (k)
), (6.45)

v i = β 0 + 1 2 [g -Hf ] 2 i α 0 + 3 2 , v τ i = β τ 0 + 1 2 g -H f 2 i α τ 0 + 3 2 , (6.46) v f j = β f 0 + 1 2 [f -m f ] 2 j α f 0 + 3 2 , (6.47) v cj = β c 0 + 1 2 f 2 j α c 0 + 3 2 . (6.48)
In this proposed method, the estimated Laplacian in each iteration is used in the optimization of the phantom, via the non-homogeneous Markovian model. The normalized Laplacian is considered for the weighted mean in the Gaussian distribution. By doing so, the neighbor pixels are distinguished via the value of the Laplacian, hence the homogeneous areas and the contours areas are separately considered. On the other hand, the optimized object f is used to construct the initialization variable for the optimization of Laplacian f in each iteration. Figure 6.6 shows the Directed Acyclic Graph (DAG) which presents the principal relationship of the updating of all the variables in the iterative updating rule in Eq.(6.44)-Eq.(6.48).

g

-Eq.(6.44)

-f ? D 2 f (0) ? g -Eq.(6.45) - f 6 q 6 - f -f ROCC D 1 ? ? FIGURE 6
.6 -DAG of proposed ROCC model.

Experiment results

Here we show the experimental results of 3D X-ray CT reconstruction of both simulated phantom and real object. Comparison with conventional methods are also derived.

6.6.2 -Simulation results of reconstructing 3D phantom with ROCC method

Implementation

In the simulation, the 3D "Shepp Logan" phantom of size 256 3 and the 3D "Head" object of size 256 3 are considered. Parallel projection and back-projection are implemented via the AS-TRA toolbox [vAPC + 16]. For each projection, 256 2 detectors are used to receive the attenuated radiation, and hence an image of size 256 2 is obtained for each projection. In our notation, N and M represent respectively the size of object f and projection dataset g. n p represents the number of projections. According to the definition, we have M = 256 2 n p .

In the proposed ROCC method, two types of gradient are used : the one of the object f , and the one of the data g. First of all, we consider the 2D case. The 3D case then can be derived with the same token.

6 -1 -1 -1 -1 -1 -1 FIGURE 6.7 -Laplacian kernel in 3D
In 3D cases, the Laplacian of f is the convolution of the object with a 3D Laplacian operator, shown in Figure 6.7. In this operator, there is a value 6 in the middle, and a -1 on the left, right, top, bottom, front and behind. It is represented by the operator :

∂ 2 ∂x 2 + ∂ 2 ∂y 2 + ∂ 2 ∂z 2 .
(6.49)

The derivation of g in 3D cases is the 2D derivation of each projection direction.

In each iteration, after the optimization of f and f , the corresponding m f , Eq.(6.33), and q, Eq.(6.34), are calculated. They are used in the optimization of the other variables.

The implementation algorithm of the ROCC method is illustrated in Algorithm 6.

Simulation results of reconstructing 3D phantom with ROCC method

The simulated Shepp Logan phantom and the real Head object are used in the simulations. Both of these two objects have 256 3 voxels. The datasets are obtained by projecting the object using parallel beam radiations. For each projection direction, 256 2 detectors are used to receive the attenuated photons and hence an image of size 256 × 256 is obtained for each projection direction. The projections are distributed evenly between 0 • and 180 • .

Experiment results

Algorithm 6 Summary of the JMAP algorithm for ROCC method 1: Input : H, D, g 2: Output : f , f 3: Initialization : 4: g = Laplacian(g)

5: f (0)
← normalized FDK 6: repeat 7:

Compute m f according to Eq.(6.33) 8:

Compute q according to Eq.(6.34) 9: Compute v c j , ∀j ∈ [1 : N ] according to Eq.(6.39) 24: until convergence or maximum iteration 6.6.3 -Initializations

Initializations

The initialization of the variables and hyperparameters of the ROCC method is discussed in this section. The reconstruction result of the object by using the FBP technique is used to initialize the variable f . Another hidden variable, f , is initialized as the Laplacian of the object f :

f = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 + ∂ 2 ∂z 2 f . (6.50)
The initialization of hyperparameters α 0 , β 0 , α τ 0 , β τ 0 , α f 0 , β f 0 , α c 0 and β c 0 depends on the chosen prior model and the prior property of the parameters.

Initialization of α 0 and β 0 α 0 and β 0 are the hyperparameters of the prior model of the noise . In the simulation, a Gaussian noise is added to datasets and consequently belong to a Gaussian distribution. However, in the real applications, the model of the noise is not known. In positron tomography the noise is normally modeled by a Poisson distribution, but in X-ray CT the Gaussian distribution is used to model noise. So that α and β are chosen such that St g is approximate to a Gaussian distribution.

As presented in Section 4.4.2, the St g distribution approaches to a Gaussian distribution when the value of α 0 and β 0 are both big. Also, a relationship between α 0 and β 0 is obtained via the variance of the noise v . In the simulations, the value of v is obtained via the value of the SNR and the variance of the projection dataset v g . As presented in Eq.(5.84), we have :

β 0 = g 2 2 M × 1 1 + 10 SN R/10 × (α 0 -1) . (6.51)
Therefore it is sufficient to initialize only the value of α 0 .

Figure 6.8 shows the influence of the value of hyperparameter α 0 on the performance of reconstruction in terms of RMSE of the reconstructed object. We can see from the figures that, for different cases of number of projections or different noise in the dataset, we always have better results when the value of α 0 has a bigger value. Furthermore, when the value of α 0 is bigger than a threshold value, its influence on the RMSE is relatively insensitive. This relative insensitiveness facilitates the initialization of the hyperparameter α 0 .

Initialization of α τ 0 and β τ 0

The noise τ is also modeled by a Gaussian distribution. For the same reason as for , the hyperparameters α τ 0 and β τ 0 are related by :

β τ 0 = g 2 2 M × 1 1 + 10 SN R/10 × (α τ 0 -1) , (6.52)
where g is the Laplacian of dataset g. Figure 6.9 shows the influence of the value of α τ 0 on the performance of reconstruction in terms of RMSE of reconstructed object. In these figures, we can see that with different value of α τ 0 , the RMSE of the reconstruction is the same. In other words, the value α τ 0 has almost no influence on the final results of reconstruction.

Experiment results

Initialization of α f 0 and β f 0

The hyperparameters α f 0 and β f 0 appear in the Inverse Gamma distribution of the parameter v f . Consequently, the variable f belongs to a St g model. A rough value of the variance of the variable f can be obtained via the initialization of the variable f :

v f = Var[ f ].
According to the property of the St g distribution, we have v f |α f 0 , β f 0 = β f 0 /(α f 0 -1), therefore we obtain :

β f 0 = v f × (α f 0 -1) . (6.53)
As there exist a relationship between these two hyperparameters, it is sufficient to initialize α f 0 . In our work, the object f under consideration, which is modeled by the St g distribution, is not sparse. Theoretically the value of α f 0 and β f 0 should be initialized with a bigger value in order to have better results.

Figure 6.10 shows the RMSE of the reconstruction results with different values of α f 0 with different number of projections and different noise level in dataset. The simulation results proved our proposition. When the dataset has a high SNR=40dB, the initilization of α f 0 is not sensitive. While there is a low SNR=20dB, the reconstruction has a better result when α f 0 is bigger.

-Initializations

Experiment results

Initialization of α c 0 and β c 0

In ROCC method, the Laplacian of the object, f , is modeled by a St g distribution depending on two hyper-parameters α c 0 and β c 0 . Here we analysis the initialization of these two hyperparameters simultaneously. The variable f is supposed to have a sparse structure, and could be realized with small values for both α c 0 and β c 0 . Consequently, the ROCC method is expected to have a better result when α c 0 and β c 0 are smaller. Figure 6.11 and Figure 6.12 show the influence of the initialization of α c 0 and β c 0 . In Figure 6.11, different colors correspond to different values of β c 0 , and in Figure 6.12 different colors correspond to different values of α c 0 . We can see from the results that, for different number of projections and different noise levels in dataset, the value of α c 0 and β c 0 have almost no influence on the RMSE of reconstruction. This property is very interesting. As we have mentioned above, when the value of α c 0 and β c 0 are smaller, the variable f will have a sparser structure. In the ROCC method, the variable f is modeled by this St g distribution. However, f also depends on the observed data g. The simulation results show us that the values of α c 0 and β c 0 don't have much influence of the quality of reconstruction because the relationship with g plays a more important role in this method.

Simulation results

The initialization of the variables and hyperparameters are defined as following :

• f (0)
is the reconstructed object from the dataset g by using the FBP method.

• f (0) is the Laplacian of the initialized f

. • α 0 = 100, and

β 0 = v × (α 0 -1) = g 2 2 M × 1 1+10 SN R/10 × (α 0 -1).

-Simulation results

10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 β c 0 • α τ 0 = 100, and

β τ 0 = v τ × (α τ 0 -1) = g 2 2 M × 1 1+10 SN R/10 × (α τ 0 -1). • α f 0 = 100, and β f 0 = v f × (α f 0 -1).
• α c 0 = 0.01 and β c 0 = 0.01. Table 6.1 shows the comparison of different metrics of evaluation of reconstruction by using different methods. Reconstruction with 180, 90, 60, 45, 36 and 18 projections are tested, and they correspond to respectively 1 projection in every 1

• , 2 • , 3 • , 4 • , 5 • and 10 • .
The ROCC method is compared with the QR and TV methods. To evaluate the simulation results, the Relative Means Square Error (RMSE), the Improvement of SNR (ISNR), the Peak SNR (PSNR) and the Structure Similarity of Image (SSIM) of the reconstruction are measured and compared. The results show that TV performs better when a low SNR dataset is used in the reconstruction. However, in the high SNR cases, and when the number of projections is insufficient, the ROCC method outperforms the TV method. The ROCC method is more robust than the two state-of-the-art methods when number of projections is limited. Figure 6.13 represents the middle slice of the reconstructed Shepp Logan phantom by using different methods, using projection sinogram of 60 projections and SNR=40dB. From the figures we can see that in the reconstructed phantom by using QR method, the edge of the phantom and the contours of different blocks in the phantom are blurred, and in the background the artifacts are obvious. The reconstructed phantom using TV and ROCC methods have clearer contours than the QR method. When comparing with TV method, we see that TV method leads to results with more blurry contours and smoother homogeneous areas. The ROCC method gives clearer contours than the two state-of-the-art methods.

In Figure 6.14, the middle slices of the reconstructed Shepp Logan phantom using different methods are presented. In this simulation, the projection sinogram with 36 projections Experiment results In this table, we can conclude that the ROCC method has a much faster convergence rate than the other two methods. For the cases of 128 projections and 64 projections, the HHBM method is slower than the TV method. However, in order to fix the regularization parameter for TV method, a test has to be done for each different parameter value, and normally more than ten times tests are needed. Consequently, the computation time for TV method is much more than the presented time.

Conclusion

In this chapter we present a Bayesian method which uses a Markovian prior model for the piecewise continuous variable. From the property of the Radon transform, the Laplacian of the object f is related with the Laplacian of the projection sinogram g by a direct model. By using simultaneously the two direct functions g = Hf + and g = H f + τ , these two variables are estimated. Meanwhile, the updated normalized Laplacian value is used as a weight coefficient in the prior model of f in each iteration. At last, all the variables and parameters are estimated simultaneously via the posterior distribution by using the JMAP estimation algorithm.

The initialization for the hyperparameters is proposed, and the insensitiveness of them to the reconstruction is proved by the simulations.

From the simulations, we conclude that the proposed ROCC method performs better for preserving the contours comparing with the conventional QR and TV methods. When the dataset is biased with a noise of SNR=40dB, the relative mean square error of the reconstructed object is smaller than QR and TV methods, and the advantage is more obvious when there is less number of projections. The RMSE value of the ROCC method by using 36 projections is the same level as the TV method by using 90 projections. This method is suitable for the parallel beam projection CT reconstruction, where the Radon Conclusion transform is used. The future work focus on the extension from the parallel beam to the cone beam projections.

7

Conclusions and Perspectives

The 3D X-ray CT reconstruction has been a hot topic in the computational imaging and medical imaging during several decades, and has been developed rapidly thanks to the development of computation science. In both medical and industrial applications, researchers' attention are changing to the aspect of the data limitation problems. The limited angle or limited number of detectors involves to an ill-posed problem. Motivated by these challenges, in this thesis, we developed several Bayesian methods for X-ray CT reconstruction.

Conclusions

In this thesis, we proposed mainly two Bayesian methods for 3D CT reconstruction. These methods are based on the context of the Non Destructive Testing (NDT) in industrial applications. The large data size in 3D applications and the expensive computational costs are considered in the algorithms.

In Chapter 5, we proposed a Bayesian method based on the forward system model g = Hf + . In this method the novelty is that we considered a hierarchical structure in the system model, with f = Dz + ξ. With the additional level of the system, the variable f is represented by a hidden variable z. For the choice of the prior model, the considered property then changes from the modelization of the piecewise continuous variable to the modelization of the sparse hidden variable. A generalized Student-t distribution is introduced in this Chapter. Thanks to the fact that the St g distribution is a marginal of a Normal-Inverse Gamma distribution, the estimation algorithms are easy to apply.

In the iterative algorithms, including the Bayesian inference, the initialization for the parameters is crucial for the final results of the reconstruction. In Section 5.4.2 and Section 6.6.3, we presented the strategy of initializing the hyper-parameters for the optimization algorithm in the Bayesian inference with the proposed model. This initialization strategy is experimentally demonstrated. In the analysis of the hyper-parameters, we observed a relative weak influence of the hyper-parameters on the behavior of the corresponding iterative algorithm. The interest of this relative weak dependency is that it offers a practical way to insure the initialization of the algorithm which is typically not trivial.

In Section 5.5, we presented another orthonormal transformation, the Dual Tree Complex Wavelet Transformation (DT-CWT). This transformation is compared with the discrete Haar transformation. In the comparison, we find that the DT-CWT transformation will not cause the bloc artifacts in the reconstructed figures and is more robust to the contamination in the coefficients. When using in the proposed Bayesian method, the DT-CWT has a better performance of reconstruction than HT. The disadvantage is that the computation of the DT-CWT is relatively more costly than using the Haar transformation.

In Section 5.8, we proposed to use another forward model in this Bayesian approach. In this model, the additive noise is separated into two parts, one is modeled by a Gaussian distribution and another modeled by a generalized Student-t distribution. The forward model considered for the projection system is g = Hf + ρ + . (7.1) In this forward system, a more complicated noise model is considered. The simulation results show that this model with the consideration of a more complex noise model derives better reconstruction results.

In Chapter 6, another strategy for the reconstruction of X-ray CT objects is presented. In this method, the relationship between the object f , the projection sinogram of the object g, the Laplacian of the object f and the Laplacian of the projection sinogram g are considered.

According to the Radon transform system, the relations of the direct projection system that we take into consideration are : g = Hf + , (7.2) g = H f + τ , (7.3) with some relations between g and g and between f and f . By considering these relationships, the Laplacian of the object f is considered as a hidden variable in the system model. A non-homogeneous Markovian model is used for the piecewise continuous object f . In this Markovian model, the hidden variable f is used as a parameter. By using this Markovian model considering the value of f , the contours and homogeneous areas are modeled with the same model with different parameters, and therefore both the continuity of the homogeneous areas and the discontinuity of the contours are both considered. This strategy preserves the edges when reconstructing the object.

Perspectives

For the future work several aspects are considered. First of all, the computational time should take into consideration. The GPU processor is only used for the projection operator and back-projection operator in this thesis. In order to accelerate the computation, more computations in the program can be implemented by using GPU processor. For example, the DT-CWT transformation presented in Chapter 5 has very good properties comparing with the Haar transformation, but has a more complicated computation. Once we implement the transformation by Perspectives using the GPU processor, the utilization of the DT-CWT or many other transformations will be less costly.

In this thesis, we used the Maximum A Posterior estimator for the optimization when the posterior distribution is obtained. There exists some other estimation algorithms, and some gives more information for the variables. For example, the Posterior Mean estimation can be used to estimate the variables. By using the Variational Bayesian Approach (VBA) method, the posterior distribution is approached by another separable distribution : p(f , θ|g) = q(f , θ) = q 1 (f )q 2 (θ), (7.4)

where θ represents all the hidden variables and parameters of the model. The mean and variance of the approximating distribution are estimated.

In the VBA method, the approaching distributions are computed by minimizing the Kullback Leibler divergence. During the computation of VBA method, a calculation of the diagonal elements of the large size matrix H T H is needed, shown in Appendix A. Some simulation results of the VBA method are presented in [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF] for the reconstruction of 2D Shepp Logan phantom by considering a stationary noise. Because of the huge size of the matrix H in the 3D problem, it is very costly if we compute directly by projection and back-projection to obtain each diagonal element. Now our group is studying on this subject and searching for algorithms which calculate approximately the diagonal elements by using the projection and back-projection properties.

In this thesis, the generalized Student-t distribution is used to define the sparse variables. As we have mentioned, the St g distribution has a heavy tailed structure. At the same time, it can have a very small variance thanks to the fact that there are two parameters for the control of the shape of this distribution. Another very useful property in the Bayesian context is that, it can be expressed as a marginal of the bivariate Normal-Inverse Gamma distribution. As the Normal distribution and the Inverse Gamma distribution are conjugate, the estimation of variables in the Bayesian methods are more simple.

In future work, in place of using the Normal-Inverse Gamma mixture distributions, we look for other similar distributions, for example the Normal variance mixture distributions. We will compare different mixture distributions and compare their properties, concerning the sparsity enforcing property of the marginal distribution and also the influence of the hyper-parameters of the marginal distributions.

One of the most important aspects in our future work is to study the influence of the initialization of hyper-parameters. In this thesis, we have observed a relative weak influence of the hyper-parameters of the Bayesian model under the prior models. In the future work we will study on the reason of this weak influence of the hyper-parameters of this model, and see if the same property appears when using the other Normal variance mixture prior models.

In a long term future work, we will search for new prior models for the X-ray CT reconstruction. We will concentrate on the ill-posed cases where the number of projection is extremely limited. For that, we should optimize the prior model and avoid the artifacts appearing in the reconstruction results. We are also going to investigate the relation of the regularization methods and the Bayesian methods. There are many regularization methods which works very well for the ill-posed reconstruction problem if a suitable regularization parameter is chosen. We can study, from the Bayesian point of view, the principle of the regularization. At the same time, we can benefit the advantages of the Bayesian method that the parameters are estimated during A Posterior Mean estimation via Variational Bayesian Approach (VBA) algorithm for HHBM method

The Posterior Mean estimation for the hierarchical Haar transformation based Bayesian method is demonstrated in this Chapter. First of all, the posterior distribution is calculated via the likelihood and the prior distributions :

p (f , z, v , v ξ , v z |g) ∝ p (g|f , v ) p (f |z, v ξ ) p (z|v z ) p (v z |α z 0 , β z 0 ) p (v |α 0 , β 0 ) p (v ξ |α ξ 0 , β ξ 0 ) ∝ N (g|Hf , V ) N (f |Dz, V ξ ) N (z|0, V z ) IG (z|α z 0 , β z 0 ) IG (v |α 0 , β 0 ) IG (v ξ |α ξ 0 , β ξ 0 ) = |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) |V z | -1 2 exp - 1 2 z T V -1 z z N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j M i=1 v -(α 0 +1) i exp -β 0 v -1 i N j=1 v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j = M i=1 v -1 2 i exp - 1 2 (g -Hf ) T V -1 (g -Hf ) N j=1 v -1 2 ξ j exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) N j=1 v -1 2 z j exp - 1 2 z T V -1 z z N j=1 v -(αz 0 +1) z j N j=1 exp -β z 0 v -1 z j M i=1 v -(α 0 +1) i M i=1 exp -β 0 v -1 i N j=1 v -(α ξ 0 +1) ξ j N j=1
exp -β ξ 0 v -1 ξ j .

(A.1) q 1j (f j )

A.1 q 1j (f j )

First step is to estimate the distribution on f j :

q 1j (f j ) ∝ exp    - 1 2 v -1 g -Hf 2 2 - 1 2 v -1 ξ f -Dz 2 2 q 1,m m =j
(fm)q 2 (z)q 3 (vz)q 4 (v )q 5( v ξ )

   ∝ exp - 1 2 v -1 q 4 (v )
g -Hf 2

2 q 1,m m =j (fm) -

1 2 v -1 ξ q 5( v ξ ) f -Dz q 1,m m =j
(fm)q 2 (z)

.

(A.8)

A.1.1 Computation of g -Hf 2

2 q 1,m m =j (f m )

We have :

g -Hf = |g 1 -H 11 f 1 -H 12 f 2 -• • • -H 1N f N | 2 + |g 2 -H 21 f 1 -H 22 f 2 -• • • -H 2N f N | 2 + • • • + |g M -H M 1 f 1 -H M 2 f 2 -• • • -H M N f N | 2 =H 2 1j f 2 j + 2H 1j f j (H 11 f 1 + H 12 f 2 + • • • + H 1N f N -H 1j f j ) -2H 1j f j g 1 + H 2 2j f 2 j + 2H 2j f j (H 21 f 1 + H 22 f 2 + • • • + H 2N f N -H 2j f j ) -2H 2j f j g 2 + H 2 M j f 2 j + 2H M j f j (H M 1 f 1 + H M 2 f 2 + • • • + H M N f N -H M j f j ) -2H M j f j g M + C = H 2 1j + H 2 2j + • • • + H 2 M j f 2 j -2 (H 1j g 1 + H 2j g 2 + • • • + H M j g M ) f j + 2f j {(H 1j H 12 + H 2j H 22 + • • • + H M j H M 2 ) f 2 + • • • + (H 1j H 1j-1 + H 2j H 2j-1 + • • • + H M j H M j-1 ) f j-1 + (H 1j H 1j+1 + H 2j H 2j+1 + • • • + H M j H M j+1 ) f j+1 + • • • + (H 1j H 1N + H 2j H 2N + • • • + H M j H M N ) f N } + C (A.9)
Computation of q 2j (z j )

where

J f = v -1 q 4 (v ) A f f 2 j + 2B f f j + v -1 ξ q 5( v ξ ) f 2 j -2 (Dµ z ) j f j = v -1 q 4 (v ) A f + v -1 ξ q 5( v ξ ) f 2 j + 2 B f v -1 q 4 (v ) -v -1 ξ q 5( v ξ ) (Dµ z ) j f j = v -1 q 4 (v ) A f + v -1 ξ q 5( v ξ )    f 2 j + 2 v -1 q 4 (v ) B f -v -1 ξ q 5( v ξ ) [Dµ z ] j v -1 q 4 (v ) A f + v -1 ξ q 5( v ξ ) f j    .
So we find that f j belongs to a Gaussian distribution, with mean and variance equal to :

µ f j = - v -1 q 4 (v ) B f -v ξ -1 q 5( v ξ ) [Dµ z ] j v -1 q 4 (v ) A f + v ξ -1
q 5( v ξ ) (A.15)

σ f j = 1 v -1 q 4 (v ) A f + v ξ -1 q 5( v ξ ) (A.16)
A.2 Computation of q 2j (z j )

From the criterion we have :

q 2j (z j ) ∝ exp    - 1 2 v -1 ξ f -Dz 2 2 - 1 2 z T V -1 z z q 1( f )q 2,m m =j
(zm)q 3 (vz)q 4 (v )q 5( v ξ )

   ∝ exp - 1 2 v -1 ξ q 5( v ξ ) f -Dz 2 2 - 1 2 v -1 z j q 3j (vz j ) z 2 j .
(A.17)

With the same calculation as we've used for g -Hf 2

2 q 1,m m =j (fm) , we obtain : So that we have :

f -Dz 2 2 q 1( f )q 2,m m
q 2j (z j ) ∝ exp - 1 2 v -1 ξ q 5( v ξ ) A z z 2 j + 2B z z j - 1 2 v -1 z j q 3j (vz j ) z 2 j ∝ exp - 1 2 J z (A.21)
The iterative optimization by using the VBA algorithm where

A f = H T H jj , B f = -[H T (g -H µ f )] j -H T H jj µ f j , A z = D T D jj , B z = -D T µ f j + D T D µ z j -D T D jj µ z j .
The computation of v -1 z j q 3j (vz j )

, v -1 q 4 (v ) and v ξ -1

q 5( v ξ ) depends on the estimated values of the hyper-parameters α z j , β z j , α , β and α ξ , β ξ . It is shown in the results that the three variances v z j , v and v ξ are all belonging to an Inverse Gamma distribution. Now let x belonging to an Inverse Gamma distribution, so :

p(x|α, β) = β α Γ(α)
x -α-1 exp -β x dx.

(A.36) So we have :

x -1 = 1 x β α Γ(α) x -α-1 exp - β x dx = β α Γ(α) x -α-2 exp - β x dx = β α Γ(α)
x -(α+1)-1 exp -β x dx So we have :

v -1 z j q 3j (vz j ) = α z j / β z j , v -1 q 4 (v ) = α / β , v ξ -1 q 5( v ξ ) = α ξ / β ξ .
(A.38)

ANNEXE B. MY CONTRIBUTIONS

other transformations may optimize the property of the proposed Bayesian method. In 7 [START_REF] Wang | Bayesian method with sparsity enforcing prior of dual-tree complex wavelet transform coefficients for X-ray CT image reconstruction[END_REF], the DT-CWT transformation coefficient is used instead of Haar transformation. A comparison of the performance of the two transformations is presented.

In 8 [START_REF] Wang | X-ray Computed Tomography simultaneous image reconstruction and contour detection using a hierarchical Markovian model[END_REF], the ROCC method is proposed. In this method, We proposed a nonhomogeneous Markovian prior model to define the piecewise homogeneous object. In this model, the prior distribution of the variable f depends on the gradient of the neighbours. Meanwhile, the gradient of the object is considered as a hidden variable in this model. In this article, a stationary noise model is considered and the 2D CT reconstruction is simulated. The 3D reconstruction with a non-stationary noise model is extended in the thesis.

In the submitted 9 IEEE journal paper, the HHBM method is considered in a 3D framework, for both simulated and real datasets. The semi-supervised extension for the initialization is demonstrated. The influence of the initialization for hyper-parameters and the theoretical basis for initializing the hyper-parameters are presented.

The following table shows the guideline.

5. 28

 28 NMSE of reconstructed image along iterations from sinogram with 128 projections and a high SNR=40dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.29 NMSE of reconstructed image along iterations from sinogram with 64 projections and a high SNR=40dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.30 NMSE of reconstructed image along iterations from sinogram with 32 projections and a high SNR=40dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.31 Generative graph of the Gaussian-St g Variable Splitting method. . . . . . . . . 5.32 Generative graph of the St g -St g Variable Splitting method. . . . . . . . . . . . 5.33 Comparison of the HHBM method and the variational splitting methods in terms of the evolution of the RMSE during reconstruction from 128 projections and SNR=40dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.34 Comparison of the HHBM method and the variational splitting methods in terms of the evolution of the RMSE during reconstruction from 64 projections and SNR=40dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 The Radon Transform with the direction of radiation defined by vector η. . . . 6.2 The original Shepp Logan figure of size 256 2 (top-left) f (x, y), its Laplacian (bottom-left) f (x, y), the projection sinogram (top-right) g(r, φ) and the Laplacian of sinogram (bottom-right) g(r, φ). . . . . . . . . . . . . . . . . . . . . . 6.3 The relations between the original figure, the projection, the Laplacian of the object and the Laplacian of the projection. . . . . . . . . . . . . . . . . . . . . 6.4 The neighbour pixels of current pixel in 2D image (left) and the neighbour voxels of current voxel in 3D object (right). . . . . . . . . . . . . . . . . . . . 6.5 Example of the inhomogeneous Markovian prior model of f . . . . .

  (α, β) Pair of parameters in St g (α 0 , β 0 ) Pair of parameters for St g distribution of v (α ξ0 , β ξ0 ) Pair of parameters for St g distribution of v ξ (α z0 , β z0 ) Pair of parameters for St g distribution of v z (α ρ0 , β ρ0 ) Pair of parameters for St g distribution of v ρ (α τ0 , β τ0 ) Pair of parameters for St g distribution of v τ (α c0 , β c0 ) Pair of parameters for St g distribution of v c λ Regularization parameter ∈ R M ×1 Additive noise of projection system ξ ∈ R N ×1 Uncertainties of transformation ρ ∈ R M ×1 Error of projection i The noise of the i-th detection φ the angle the normal vector to l makes with the x axis ν Degrees of freedom in St θ Parameters in the Bayesian model κ Hyper-parameters in the Bayesian model γ • , γ f , γ z Gradient descent step length Γ(•) Gamma function η Unit vector of r η 1 Projection of η onto axis x η 2 Projection of η onto axis y τ ∈ R M ×1

  Figure 3.1 and Figure 3.2 illustrate some examples of the X-ray CT reconstruction images in industrial and medical applications.

FIGURE 3

 3 FIGURE 3.1 -The industrial components tested by X-ray CT reconstruction. Top : the outward appearance of the components, bottom : the corresponding reconstructed objects. (Image on the left : http ://www.materialstoday.com/hardmetals-and-ceramics/features/2c-ceramics-moves-into-the-industrialreality-zone/ ; image on the right : http ://www.npl.co.uk/science-technology/dimensional/x-ray-computed-tomography.)

The

  FIGURE 3.2 -The medical images detected by X-ray CT reconstruction. Brain, bone, teeth, etc. (Image on the top : http ://www.southlakeregional.org/Default.aspx ?cid=784 ; image on the bottom left : http ://emedicine.medscape.com/article/355892-overview ; image on the bottom right : http ://texas-dental-implants.com/wp-content/uploads/2015/02/PAN-X-rays-Small.jpeg.)

Figure 3

 3 Figure 3.3 illustrate the Radon transform of a continuous function f (x, y). For any point (x, y) on the integral line we have : r = x cos φ + y sin φ.(3.4)

  Figure 3.4 illustrate the Radon transform and its dual transform pair. The dual Radon transform is also called the backprojection (BP) in CT, which is considered as the adjoint operator of the projection. The relation of the backprojection and the real object is : b(x, y) = f (x, y) * * h(x, y), (3.7)

FIGURE 3

 3 FIGURE 3.4 -The Radon transform and Backprojection operator.

FIGURE 3

 3 FIGURE 3.6 -The Back-projection Filtering reconstruction algorithm.

Figure 3 .

 3 Figure 3.8 shows the original Shepp Logan figure of size 256 2 and the projection sinogram with 64 projections evenly distributed in [0, π]. The reconstruction results by using the conventional Back-projection (BP) and Filtered Back-projection (FBP) methods are also shown.

3.3. 3 -

 3 FIGURE 3.9 -The reconstructed figure by using 64 projections distributed in [0, π], 32 projections distributed in [0, π] and 64 projections distributed in [0, π 2 ].

FIGURE 3 .

 3 FIGURE 3.11 -Reconstruction of 2D Shepp Logan phantom of size 256 2 by using the SIRT method, with respectively 180 projections (left) and 60 projections (right).

  FIGURE 3.12 -The influence of the regularization parameter on the relative mean square error of the reconstruction results in the Quadratic Regularization method with 128 projections (top) or 64 projections (bottom) with a noise of SNR=40dB (left) or SNR=20dB (right).

  FIGURE 3.13 -The influence of the regularization parameter on the relative mean square error of the reconstruction results in the Total Variation method with 128 projections (top) or 64 projections (bottom) with a noise of SNR=40dB (left) or SNR=20dB (right).

Figure 3 .

 3 16 illustrate the basic Bayesian model corresponding to Eq.(3.38). In this model, the parameters θ are unknown and depends on the fixed hyper-parameters κ = [κ 1 , κ 2 ].

FIGURE 4

 4 FIGURE 4.1 -The positions of current pixel r and neighbor pixels r by considering left : 4 neighbor pixels or right : 8 neighbor pixels.

FIGURE 4

 4 FIGURE 4.2 -The piecewise constant objects and their gradient.

FIGURE 4

 4 FIGURE 4.3 -Example of inhomogeneous Markovian model.

FIGURE 4

 4 FIGURE 4.4 -The wavelet shape of the Haar functions.

. 19 )

 19 Sparsity enforcing prior distributions and we obtain : DD T = I. (4.20) Consequently D T = D -1 . In Figure 4.5, the Shepp Logan figure of size 256 2 and its multilevel Haar transformation coefficients are shown. The scale of the transformation coefficients are adapted such that they are visually clearly.

  FIGURE 4.5 -The Shepp Logan figure of size 256 2 (left), its 2-level Haar transformation coefficient (middle) and its 4-level Haar transformation coefficient (right).

FIGURE 4 . 6 -

 46 FIGURE 4.6 -The comparison of the pdfs of the Normal distribution, the Student-t distribution and the generalized Student-t distribution.

CHAPITRE 5 .

 5 SPARSITY ENFORCING UNSUPERVISED HIERARCHICAL MODEL [FBS02, DB95, RVJ + 06, LNDD07]. Some commonly used transforms are for example the discrete wavelet transform (DWT) [AU96, TSP + 08], the discrete cosinus transformation (DCT) [RY14, SK13], the discrete Haar transformation [SF03], the discrete curvlet transformation [CDDY06], the dual-tree Complex Wavelet Transform (DT-CWT) [SBK05, Kin98].

Figure 5 .

 5 Figure 5.1 shows the 2 -level, 4 -level and 5 -level Haar transformation coefficient of Shepp Logan phantom and the Head object, both of size 256 × 256. The scale of these figures have been adapted to be visually more clearly. In general the degree of sparsity increases when the transformation level augments.

FIGURE 5 . 2 -

 52 FIGURE 5.2 -The histogram of the multilevel Haar transformation of the Shepp Logan phantom. Comparison of the histogram of different transformation levels.

FIGURE 5 . 4 -

 54 FIGURE 5.4 -The pdfs of Inverse Gamma distribution with different parameter values.

  z (k+1) = z (k) -γ (k) z ∇J( z (k) ) 10: until convergence or max iterationThe entire optimization algorithm which estimate all the variables are shown in Algorithm 2.

  FIGURE 5.6 -Middle slice of reconstructed Shepp Logan phantom of size 256 3 from projection sinogram with 128 projections (top) or 64 projections (bottom) and a noise of SNR=40dB (left) or SNR=20dB (right).

  and v ξ j are iid. In this model, both z and ξ are sparse and belong to a St g distribution.

) 10 :

 10 until convergence or maximum iteration Algorithm 4 Gradient descent algorithm for z 1: Input : H, D, g, f , z 0 2: Output : z 3: Initialization : z ini = z 0 4:

  z (0) ← l-level Haar transformation of f (0) 7: V z , V and V ξ according to Eq.(5.77)-Eq.(5.79) v z j , ∀j ∈ [1 : N ] according to Eq.(5.77) 12:Compute v i , ∀i ∈ [1 : M ] according to Eq.(5.78) 13:

FIGURE 5

 5 FIGURE 5.8 -The illustration of rank r of a 5-level discrete Haar transformation coefficient.

FIGURE 5

 5 FIGURE 5.9 -Influence of hyperparameter α 0 on RMSE of final reconstruction results for different number of projections and noise.

FIGURE 5 FIGURE 5

 55 FIGURE 5.10 -Influence of hyperparameter α ξ 0 , for different values of β ξ 0 , on RMSE of reconstruction results for different number of projections and noise. Each different color corresponds to a different value of β ξ 0 .

Figure 5 .

 5 Figure 5.17 and Figure 5.18 show the comparison of the performance in terms of RMSE of different method with a high SNR=40dB and a low SNR=20dB. In this comparison, four cases of limited angle of projections are considered. They are respectively : 45, 90, 135 and 180 projections evenly distributed in respectively [0• , 45 • ], [0 • , 90 • ], [0 • , 135 • ] and [0 • , 180 • ].Thus, there is one projection every 1 • . From these two figures, we conclude that the proposed HHBM method stay more robust than the other two conventional methods when there is limited angle of projections.

Dual-tree complexFIGURE 5

 5 FIGURE 5.16 -Slice of reconstructed 3D Shepp Logan phantom and 3D Head object, with 90 projections evenly distributed in [0 • , 90 • ].

FIGURE 5

 5 FIGURE 5.18 -The performance of different methods for reconstructing Shepp Logan phantom in terms of RMSE with different limited angles of projections and a low SNR=20dB.

FIGURE 5

 5 FIGURE 5.21 -The 5-level Haar transformation coefficient and 5-level DT-CWT transformation coefficient of the 2D Shepp Logan phantom of size 256 × 256. The scales are adapted in order to be visually clearer.

FIGURE 5

 5 FIGURE 5.23 -Inverse transformed images with different noise added to coefficients by using different transformations, with the Relative Mean Square Error (RMSE) of image shown below.

  FIGURE 5.24 -The images with different shape of textures.

  .25, the Original Shepp Logan figure and its projection sinogram are shown. For the illustration of reconstruction phantom, we show the zone of the phantom which is in the red block areas shown in Figure 5.25.

  FIGURE 5.25 -Original Shepp Logan image of size 256 2 and its (b) sinogram with 128 projections and (c) the same data with additive noise with SNR=20dB.

FIGURE 5 FIGURE 5 FIGURE 5

 555 FIGURE 5.26 -The middle slice of reconstructed Shepp Logan phantom by using different reconstruction methods. The reconstructed form projection dataset with top : 64 projections and bottom : 32 projections respectively.

Figure 5 .

 5 Figure5.32 illustrates the Generative graph of this Variable Splitting model.

FIGURE 5 FIGURE 5

 55 FIGURE 5.33 -Comparison of the HHBM method and the variational splitting methods in terms of the evolution of the RMSE during reconstruction from 128 projections and SNR=40dB.

  FIGURE 6.1 -The Radon Transform with the direction of radiation defined by vector η.

FIGURE 6

 6 FIGURE 6.2 -The original Shepp Logan figure of size 256 2 (top-left) f (x, y), its Laplacian (bottom-left) f (x, y), the projection sinogram (top-right) g(r, φ) and the Laplacian of sinogram (bottom-right) g(r, φ).

  γ f according to Eq.(6.42) 11: Compute ∇J( f ) from Eq.(6.38) 12: f = f -γ f ∇J( f ) γ f according to Eq.(6.43) 17: Compute ∇J( f ) from Eq.(6.38) 18: f = f -γ f ∇J( f ) 19: until convergence or maximum iteration 20: Compute v i , ∀i ∈ [1 : M ] according to Eq.(6.39) 21: Compute v τ i , ∀i ∈ [1 : M ] according to Eq.(6.39) 22:Compute v f j , ∀j ∈ [1 : N ] according to Eq.(6.39) 23:

FIGURE 6

 6 FIGURE 6.8 -Influence of hyperparameter α 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise.

FIGURE 6

 6 FIGURE 6.9 -Influence of hyperparameter α τ 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise.

FIGURE 6

 6 FIGURE 6.10 -Influence of hyperparameter α f 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise.

FIGURE 6

 6 FIGURE 6.11 -Influence of hyperparameter α c 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise.

FIGURE 6

 6 FIGURE 6.12 -Influence of hyperparameter β c 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise.

  =j (zm) = A z z 2 j + 2B z z j + C, (A.18)whereA z = D T D jj (A.19) B z = -D T µ f j + D T Dµ z j -D T D jj µ z j . (A.20)

  ANNEXE C. AUTHOR'S PUBLICATIONS DURING 3 YEARS OF PHD Conference on (pp. 883-887), IEEE. • [WMDG17b] L.Wang, A.Mohammad-Djafari and N.Gac, Bayesian X-ray computed tomography using a three-level hierarchical prior model, in AIP Conference Proceedings (Vol. 1853, No. 1, p. 060003), AIP Publishing. • [WMDG17c] L.Wang, A.Mohammad-Djafari and N.Gac, X-ray Computed Tomography simultaneous image reconstruction and contour detection using a hierarchical Markovian model, in Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on (pp. 6070-6074), IEEE. • [WMDGD17] L.Wang, A.Mohammad-Djafari, N.Gac and M.Dumitru, 3D X-ray Computed Tomography reconstruction using sparsity enforcing hierarchical model based on Haar transformation, in The 2017 International Conference on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine. • [WMDG17a] L.Wang, A.Mohammad-Djafari and N.Gac, Bayesian method with sparsity enforcing prior of dual-tree complex wavelet transform coefficients for X-ray CT image reconstruction, in The 25th European Signal Processing Conference (EUSIPCO). • [DWMDG17] M.Dumitru, L.Wang, A.Mohammad-Djafari and N.Gac, Model selection in the sparsity context for inverse problems in Bayesian framework, in 37th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering. • [DLGMD17] M.Dumitru, L.Wang, N.Gac and A.Mohammad-Djafari, Performance comparison of Bayesian iterative algorithms for three classes of sparsity enforcing priors with application in computed tomography, in Image Processing (ICIP), 2017 IEEE International Conference on, IEEE. • [DGWMD17] M.Dumitru, N.Gac, L.Wang and A.Mohammad-Djafari, Unsupervised sparsity enforcing iterative algorithms for 3D image reconstruction in X-ray computed tomography, in The 2017 International Conference on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine. Papers published in national conferences : • [WGMD17] L.Wang, N.Gac and A.Mohammad-Djafari, Reconstruction 3D en tomographie à rayons X à l'aide d'un modèle a priori hiérarchique utilisant la transformation de Haar, in Colloque GRETSI 2017. • [DWGMD17] M.Dumitru, L.Wang, N.Gac and A.Mohammad-Djafari, Comparaison des performances d'algorithmes itératifs bayésiens basés sur trois classes de modeles a priori parcimonieux appliqués à la reconstruction tomographique, in Colloque GRETSI 2017. Titre : Reconstruction d'image en tomographie 3D pour des applications en Contrôle Non Destructif (CND) Mots clés : Tomographie, Problème inverse, Méthode Bayesian, Structure hiérarchique Résumé : La tomographie 3D est largement utilisée dans l'imagerie médicale ainsi que dans le Contrôle Non Destructif (CND) pour l'industrie. Pour ces deux applications, il est nécessaire de réduire le nombre de projections. Dans certains cas, la reconstruction doit être faite avec un nombre d'angle de projections limité. Nous sommes donc presque toujours dans la situation de problèmes inversés mal posés. En plus, les données mesurées sont toujours avec des erreurs. Le rôle des méthodes probabilistes et de la modélisation a priori devient crucial.

  The middle slice of reconstructed Shepp Logan phantom by using different reconstruction methods. The reconstructed form projection dataset with top : 64 projections and bottom : 32 projections respectively. . . . . . . . . . . . . . . .5.27Comparison of a zone of profile with TV method, from dataset of 64 projections on the top, and 32 projections on the bottom. . . . . . . . . . . . . . . . . . . .

	5.26
	Table des figures	VII
	Liste des tableaux
	1 Résumé en franc ¸ais

  . . . . . 6.6 DAG of proposed ROCC model. . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Laplacian kernel in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Influence of hyperparameter α 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise. . . . . . 6.9 Influence of hyperparameter α τ 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise. . . . . . .12 Influence of hyperparameter β c 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise. . . . . . 6.13 The reconstructed Shepp Logan phantom of size 256 3 with 60 projections and SNR=40dB, by using different reconstruction methods. Red curves represent the profile at the position of the blue line. . . . . . . . . . . . . . . . . . . . . 6.14 The reconstructed Shepp Logan phantom of size 256 3 with 36 projections and SNR=40dB, by using different reconstruction methods. Red curves represent the profile at the position of the blue line. . . . . . . . . . . . . . . . . . . . . 6.15 The profiles of the middle slice of reconstructed Shepp Logan phantom of size 256 3 with dataset of 60 projections and SNR=40dB, by using different reconstruction methods. The bottom figure is a zone of contour of the top figure. . . . 6.16 The profiles of the middle slice of reconstructed Shepp Logan phantom of size 256 3 with dataset of 18 projections and SNR=40dB, by using different reconstruction methods. The bottom figure is a zone of contour of the top figure. . . . 6.17 The Relative Mean Square Error (RMSE) of the reconstructed Shepp Logan phantom with the dataset of different numbers of projections and SNR=40dB, by using different reconstruction methods. . . . . . . . . . . . . . . . . . . . . Comparison of relative error of reconstructed phantom with 50 global iterations (10 iterations of gradient descent in each global iteration) and the computation time of each iteration of different methods. The value of regularization parameter are respectively λ QR = 10 for SNR=40dB, λ QR = 600 for SNR=20dB, λ

	Notations :
	Liste des tableaux Alternating Direction Method of Multipliers Algebraic Reconstruction Technique Back-projection Conjugate Gradient Back-projection result Back-projection operator Linear transformation operator Chap4 : Multilevel Haar transformation operator Convolution kernel for 1D derivation Convolution kernel for 2D derivation D 1 ∈ R M ×M Convolution operator for 1D derivation ADMM ART BP CG b(x, y) B D d 1 d 2 5.1 Table of abbreviations CT D 2 ∈ R N ×N Convolution operator for 2D derivation Computed Tomography f (x, y) attenuation coefficient in a 2D image DT-CWT Dual-Tree Complex Wavelet Transformation DWT f (x, y, z) attenuation coefficient in a 3D object Discrete Wavelet Transformation EM f Estimation of object Expectation-Maximization f ∈ R N ×1
	FBP	Filtered Back-projection
	FT	Fourier Transformation
	GCV	Generalized Cross Validation
	GD	Gradient Descent
	GG	Generalized Gaussian
	GPU	Graphics Processing Unit
	HHBM	Hierarchical Haar based Bayesian Method
	HT	Haar transformation
	ISNR	Improvement of Signal to Noise Ratio
	JMAP	Joint Maximum A Posterior
	KL	Kullback-Leibler divergence
	LASSO	Least Absolute Shrinkage and Selection Operator
	LS	Least Square
	ML	Maximum Likelihood
	NDT	Non Destructive Testing
	OS-EM	Ordered Subset Expectation Maximization
	PET-CT Positron Emission Tomography-Computed Tomography
	pdf	probability density function
	PM	Posterior Mean
	PSNR	Peak Signal to Noise Ratio
	QR	Quadratic Regularization
	RMSE	Relative Mean Square Error
	ROCC	Reconstruction of Object by Considering Contours
	RT	Radon Transform
	SART	Simultaneous Algebraic Reconstruction Technique
	SEM	Stochastic Expectation Maximization
	SIRT	Simultaneous Iterative Reconstruction Technique
	SNR	Signal to Noise Ratio
	SSIM	Structural SIMilarity
	St-g	Generalized Student-t distribution
	St-t	Student-t distribution
	TV	Total Variation
	VBA	Variational Bayesian Approach
	VS	Variable Splitting
	2

6.10 Influence of hyperparameter α f 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise. . . . . . 6.11 Influence of hyperparameter α c 0 in ROCC method on RMSE of final reconstruction results for different number of projections and different noise. . . . . . 66.18 The Relative Mean Square Error (RMSE) of the reconstructed Shepp Logan phantom with the dataset of different numbers of projections and SNR=20dB, by using different reconstruction methods. . . . . . . . . . . . . . . . . . . . . T V = 50 for SNR=40dB, λ T V = 100 for SNR=20dB. . . . . . . . . . . . . . 80 5.2 Truncating the Transform Coefficients z . . . . . . . . . . . . . . . . . . . . . 93 5.3 Adding noise to the Transform Coefficients z . . . . . . . . . . . . . . . . . . 93 6.1 Comparison of RMSE, ISNR, PSNR and SSIM of reconstructed phantom with 50 globle iterations (10 iterations for gradient descent in each of globle iterations) by using QR, TV or ROCC methods. The values of regularization parameters are respectively λ QR = 10 and λ T V = 50 for SNR=40dB, λ QR = 600 and λ T V = 100 for SNR=20dB. . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.2 Comparison of computation time by using different methods. * The threshold of RMSE implies the critical condition of the end of iteration. When RMSE of the reconstruction is smaller than the threshold, the iteration stops. ** For the TV method, the regularization parameter is an optimized parameter, which is fixed by comparing the RMSE of reconstruction with different parameter values. In the table, the number in the parentheses represents the number of iterations till convergence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 B.1 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

  The heavy tailed distributions[START_REF] Borisovič | Heavy tailed distributions[END_REF], for example the Student-t distribution. The tails of this type of distributions are not exponentially bounded, with the definition :

	4.4.2 -Generalized Student-t distribution
	3.
	.22)

TABLE 5 .

 5 1 -Comparison of relative error of reconstructed phantom with 50 global iterations (10 iterations of gradient descent in each global iteration) and the computation time of each iteration of different methods. The value of regularization parameter are respectively λ QR = 10 for SNR=40dB, λ QR = 600 for SNR=20dB, λ T V = 50 for SNR=40dB, λ T V = 100 for SNR=20dB.

							256*256*256
				180 projections					90 projections
			40dB			20dB			40dB			20dB
		QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM
	RMSE 0.0236 0.0114	0.0069	0.1309	0.0209	0.0755 0.0401 0.0212	0.0092	0.1558	0.0491	0.1117
	ISNR 5.5584 8.7217 10.9346 7.2024 15.1775 10.2162 6.6136 9.3832 12.9973 8.4583 13.4765 9.9056
	PSNR 30.0675 33.2308 35.4437 22.6318 30.6069 25.0209 27.7743 30.5439 34.1579 21.8754 26.8937 23.3227
	SSIM	0.9999 0.9999	1.0000	0.9992	0.9999	0.9995 0.9997 0.9999	0.9999	0.9990	0.9997	0.9993
				60 projections					45 projections
			40dB			20dB			40dB			20dB
		QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM
	RMSE 0.0636 0.0321	0.0107	0.1656	0.0753	0.1293 0.0904 0.0474	0.0132	0.1854	0.0901	0.1414
	ISNR 9.3826 12.3480 17.1346 9.1492 12.5701 10.2226 10.3301 13.1308 18.6839 10.0137 13.1476 11.1916
	PSNR 25.7693 28.7347 33.5214 21.6116 25.0325 22.6849 24.2404 27.0412 32.5942 21.1195 24.2535 22.2974
	SSIM	0.9996 0.9995	0.9999	0.9990	0.9995	0.9992 0.9994 0.9997	0.9999	0.9988	0.9994	0.9991
				36 projections					18 projections
			40dB			20dB			40dB			20dB
		QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM
	RMSE 0.1177 0.0680	0.0169	0.1957	0.1116	0.1500 0.2581 0.2104	0.0574	0.2907	0.2313	0.2014
	ISNR 10.6591 13.0424 19.0933 10.8633 13.3032 12.0187 10.7122 11.5992 17.2373 10.8088 11.8022 12.4036
	PSNR 23.0949 25.4783 31.5292 20.8865 23.3264 22.0420 19.6263 20.5133 26.1514 19.1085 20.1020 20.7033
	SSIM	0.9993 0.9996	0.9999	0.9988	0.9993	0.9990 0.9983 0.9987	0.9996	0.9981	0.9985	0.9987

4 -Simulation results with limited angle of projections

		0.3
		0.25				QR TV
	RMSE	0.15 0.2				HHBM
		0.1
		0.05
		0
		18	36	45	60	90	180
			number of projections
	FIGURE 5.15 -The performance of different methods for reconstructing Shepp Logan
	phantom in terms of RMSE with different number of projections evenly distributed in
			[0
		0.3
		0.25				QR TV
	RMSE	0.15 0.2				HHBM
		0.1
		0.05
		0
		18	36	45	60	90	180
			number of projections

FIGURE 5.14 -The performance of different methods for reconstructing Shepp Logan phantom in terms of RMSE with different number of projections evenly distributed in [0 • , 180 • ] and a high SNR=40dB.

5.4.• , 180

• ] and a low SNR=20dB.

TABLE 5 .

 5 2 -Truncating the Transform Coefficients z

	Period	Stones	Ground	Stairs
	truncate Haar	DT-CWT Haar	DT-CWT Haar	DT-CWT Haar	DT-CWT
	80%

TABLE 5 .

 5 3 -Adding noise to the Transform Coefficients z

		Period		Stones		Ground	Stairs
	noise Haar	DT-CWT Haar	DT-CWT Haar	DT-CWT Haar	DT-CWT
	5dB	0.3166 0.1596	0.3180 0.1586	0.3142 0.1592	0.3156 0.1563
	10dB 0.0990 0.0495	0.1008 0.0501	0.1004 0.0497	0.0993 0.0498
	15dB 0.0316 0.0158	0.0313 0.0159	0.0315 0.0160	0.0315 0.0159
	20dB 0.0100 0.0050	0.0100 0.0050	0.0100 0.0050	0.0100 0.0050
	30dB 0.0010 0.0005	0.0010 0.0005	0.0010 0.0005	0.0110 0.0005

TABLE 6 .

 6 1 -Comparison of RMSE, ISNR, PSNR and SSIM of reconstructed phantom with 50 globle iterations (10 iterations for gradient descent in each of globle iterations) by using QR, TV or ROCC methods. The values of regularization parameters are respectively λ QR = 10 and λ T V = 50 for SNR=40dB, λ QR = 600 and λ T V = 100 for SNR=20dB. .1775 9.9937 6.6136 9.3832 10.1605 8.4583 13.4765 10.7756 PSNR 30.0675 33.2308 33.2443 22.6318 30.6069 25.4231 27.7743 30.5439 31.3211 21.8754 26.8937 24.1928 ISNR 9.3826 12.3480 14.0209 9.1492 12.5701 10.9191 10.3301 13.1308 15.7625 10.0137 13.1476 11.5572 PSNR 25.7693 28.7347 30.4076 21.6116 25.0325 23.3814 24.2404 27.0412 29.6728 21.1195 24.2535 22.6631 ISNR 10.6591 13.0424 16.5786 10.8633 13.3032 12.1617 10.7122 11.5992 16.8020 10.8088 11.8022 11.8279 PSNR 23.0949 25.4783 29.0145 20.8865 23.3264 22.1849 19.6263 20.5133 26.2758 19.1085 20.1020 20.

							256*256*256
				180 projections					90 projections
			40dB			20dB			40dB			20dB
		QR	TV	ROCC	QR	TV	ROCC	QR	TV	ROCC	QR	TV	ROCC
	RMSE 0.0236 0.0114	0.0114	0.1309	0.0209	0.0688 0.0401 0.0212	0.0177	0.1558	0.0491	0.0914
	ISNR 5.5584 8.7217 7.2024 15SSIM 8.6386 0.9999 0.9999 0.9999 0.9992 0.9999	0.9996 0.9997 0.9999	0.9999	0.9990	0.9997	0.9994
				60 projections					45 projections
			40dB			20dB			40dB			20dB
		QR	TV	ROCC	QR	TV	ROCC	QR	TV	ROCC	QR	TV	ROCC
	RMSE 0.0636 0.0321	0.0218	0.1656	0.0753	0.1102 0.0904 0.0474	0.0259	0.1854	0.0901	0.1300
	SSIM	0.9996 0.9995	0.9999	0.9990	0.9995	0.9993 0.9994 0.9997	0.9998	0.9988	0.9994	0.9992
				36 projections					18 projections
			40dB			20dB			40dB			20dB
		QR	TV	ROCC	QR	TV	ROCC	QR	TV	ROCC	QR	TV	ROCC
	RMSE 0.1177 0.0680	0.0301	0.1957	0.1116	0.1451 0.2581 0.2104	0.0566	0.2907	0.2313	0.2252
													2758
	SSIM	0.9993 0.9996	0.9998	0.9988	0.9993	0.9991 0.9983 0.9987	0.9996	0.9981	0.9985	0.9986

FIGURE 6.13 -The reconstructed Shepp Logan phantom of size 256 3 with 60 projections and SNR=40dB, by using different reconstruction methods. Red curves represent the profile at the position of the blue line. FIGURE 6.15 -The profiles of the middle slice of reconstructed Shepp Logan phantom of size 256 3 with dataset of 60 projections and SNR=40dB, by using different reconstruction methods. The bottom figure is a zone of contour of the top figure.

FIGURE 6.16 -The profiles of the middle slice of reconstructed Shepp Logan phantom of size 256 3 with dataset of 18 projections and SNR=40dB, by using different reconstruction methods. The bottom figure is a zone of contour of the top figure.

TABLE 6 .

 6 2 -Comparison of computation time by using different methods. * The threshold of RMSE implies the critical condition of the end of iteration. When RMSE of the reconstruction is smaller than the threshold, the iteration stops. ** For the TV method, the regularization parameter is an optimized parameter, which is fixed by comparing the RMSE of reconstruction with different parameter values. In the table, the number in the parentheses represents the number of iterations till convergence.Logan phantom of size 256 3 , by applying respectively 128, 64 and 32 projections. The computation time represents the total time needed such that the RMSE of estimation f is smaller than the threshold of RMSE. For TV method, the optimal regularization parameter is fixed by comparing the RMSE of reconstruction results with different parameter values. The threshold of RMSE of the critical condition is chosen according to the RMSE of the final reconstruction results of all the reconstruction methods.

	Projection number	128	64	32
	Threshold of RMSE*	0.02	0.04	0.1
	TV**	722.8s (40) † 754.65s (43) 682.08s (42)
	HHBM	1532.3s (13) 822.56s (8)	496.3s (5)
	ROCC	673.2s (20)	440.7s (15)	82.65s (3)

  Pour la modélisation a priori, en particulier dans les applications NDT, l'objet à l'examen est composé de plusieurs matériaux homogènes, avec plusieurs blocs continus séparés par des discontinuités et des contours. Ce type d'objet est dit continu par morceaux. L'objet de cette thèse est sur la reconstruction des objets continus ou constants par morceaux, ou plus généralement homogènes par morceaux. En résumé, deux méthodes principales sont proposées dans le contexte de l'inférence bayésienne.La première méthode consiste à reconstruire l'objet en imposant que sa transformée de Haar soit parcimonieuse. Un modèle bayésien hiérarchique est proposé. Dans cette méthode, les variables et les paramètres sont estimés et les hyperparamètres sont initialisés selon la définition des modèles antérieurs.La deuxième méthode reconstruit les objets en estimant simultanément les contours. L'objet continu par morceaux est modélisé par un modèle markovien non-homogène, qui dépend du gradient de l'objet, et le gradient dépend aussi de l'estimation de l'objet. Cette méthode est également semisupervisé, avec les paramètres estimés automatiquement.Ces méthodes sont adaptées aux reconstructions de grande taille de données 3D, dans lesquelles le processeur GPU est utilisé pour accélérer les calculs. Les méthodes sont validées avec des données simulées et des données réelles, et sont comparées avec plusieurs méthodes classiques.

Dual-tree complex wavelet transform (DT-CWT) Based Hierarchical Bayesian Method

Comparison with Haar transform

The DT-CWT transform and the Haar transform is compared by analysing several different cases of biased transform coefficient. Error of the coefficient will have an influence on the inverse transformation results. We search for a transformation which is more robust against the noise of coefficient, thus the corresponding inverse-transformed data depends less on the noise of the coefficient.

In Figure 5.21, the HT coefficient and the DT-CWT coefficient are presented. It shows that the DT-CWT transformation import a redundancy to the coefficient. Although the size of the coefficient is two times the HT coefficient, thanks to the development of the computational capacity, the computation can be fulfilled.

In order to compare the HT and DT-CWT, we consider two types of biased transformation coefficients. First of all we truncate the coefficient and compare the recontructed image. Then, a Gaussian noise is added to the coefficient and we compare the reconstructed results of these two transformations.

Truncating parts of coefficients z

In image processing, the wavelet transformations are frequently used in many applications, for example the denoising problem, the restoration, etc. In some of these applications, the coefficient of transformations are truncated. So here we analyse the influence of the truncation of coefficient to the image.

In this simulation, the transformation coefficients are truncated. Noting that for Haar transformation the coefficients matrix has size N = N x × N y , and for DT-CWT is 2N = N x × 2N y . For both of the cases, the 80%, 90%, 95% and 99% smallest coefficients are removed. The image of the inverse transformation is obtained and compared.

The first comparison is given by doing the following steps :

1. computing the transform of image z = D -1 f , 2. truncating parts of the minimum coefficients in z (80%, 90%, 95% and 99%), 3. using the truncated coefficients to reconstruct the image f . The Relative Mean Square Error (RMSE) of f are given below the corresponding images which are given in Figure 5.22. From the results, we can see that when 80% of the smallest coefficients are truncated, the reconstructed image is still clear by using both transformation. When 90% coefficients are truncated, the reconstructed figures is still clear by using the Haar transformation, but the figure becomes blurry when the DT-CWT is used. When more coefficients are truncated, for example 95% and 99%, both of the two transformations will reconstruct a blurry image. The discrete Haar transformation generates square textures in the reconstructed image, and the DT-CWT transformation gives a smoother image at the contours. From this simulation we can conclude that, the coefficients of HT are more compressed than DT-CWT. The results demonstrate that the proportion of the non-null coefficient is less than 10% while in DT-CWT it is more than 10%. However, when there are more truncated coefficients, the DT-CWT outperforms the HT. Visually, there is also less artifacts in the results of the inverse DT-CWT transformation.

Variable Splitting Algorithm for HHBM and the negative logarithm of the posterior distribution is :

(5.115)

The unknowns are optimized alternately :

(5.119)

(5.122)

As indicated previously, the estimation of variables f and z are realized by using the descent

Variable Splitting Algorithm for HHBM

• For HHBM-SS method, the following relations are considered :

(5.131)

In our work, we initialized different values for α 0 and α ρ0 in order to realize different prior distributions. We set α ρ0 = 2.01 in order to have a sparse structure for ρ, and we set α 0 = 100 in order to have a Gaussian structured noise . Next we use these two different noise splitting methods for the reconstruction of Shepp Logan figure and compare the simulation results.

Simulation results of the noise splitting method

In the simulations, we compared the proposed HHBM method with the variable splitting method. In our work, we considered to split the noise into two parts, and two strategies of the prior modeling are proposed. In the first strategy, called HHBM-GS, the projection noise ρ is modeled by a St g distribution while the detection noise is modeled by a Gaussian distribution. In the second strategy, called HHBM-SS, two St g distributions are used to model these two noises, and ρ, with different values for the initialization of the corresponding hyperparameters.

In the simulations, the Shepp Logan phantom of size 256 2 is used as the original figure. The projections are evenly distributed from 0 • to 180 • . The initialization for the figure f

is the reconstructed figure of the FBP method.

Figure 5.33 shows the comparison of the convergence of the RMSE of the reconstruction by using the HHBM method, the variable splitting method with two St g models (HHBM-SS) and the variable splitting method with a Gaussian and a St g model (HHBM-GS), with 128 projections and 64 projections, respectively.

From the comparisons, we can see that, by considering two St g noises, the variable splitting method (HHBM-SS) has a slightly better reconstruction result than the HHBM method in terms of RMSE value. By modeling one of the noise by a Gaussian distribution and the other by the St g distribution, the reconstruction results surpasses that of the HHBM method and the HHBM-SS method, with a faster convergence than the other methods.

Conclusions and perspectives

In this chapter, the HHBM method is presented. In this method, a hierarchical Bayesian model is defined, in which a hidden variable is introduced. The posterior distribution in the Bayesian method is obtained from the hierarchical structured prior model. From the posterior distribution, the unknown variables and parameters are estimated via MAP algorithm. Compared with the regularization approaching methods, the proposed Bayesian method estimates the parameters and variables simultaneously. The hyperparameters which defines the prior distribution of the parameters are fixed during initialization. A proposition of the initialization for the hyperparameters is presented in this chapter, and the simulation results proved a relative sensitiveness of the hyperparameters' values. Simulations with both limited number of projections and limited angle of projections are done, and it is proved that the proposed Bayesian method

Conclusions and perspectives

is more robust than the conventional quadratic regularization method or total variation method when the number or angle of projection is limited.

We also compared the Dual-Tree complex wavelet transformation (DT-CWT) with the Haar transformation. From the comparison, we observed that the DT-CWT transformation is more robust to the contamination than the Haar transformaion. Thanks to that, the Bayesian method BH-DTCWT based on this wavelet transformation outperforms the HHBM method using Haar wavelet transformation. The short term future work devote to find some other transformation which is computationally faster and considers more wavelet directions instead of only two directions in Haar transformation. The DT-CWT could be implemented in GPU for 3D simulations in the future in order to accelerate the computation.

Based on the proposed HHBM method, a noise splitting method is proposed and presented in this chapter. The splitting of noise considers the noise in two parts : a Gaussian noise and a St g noise. By considering these two noise, the reconstruction surpasses the HHBM in terms of convergence speed and the quality of reconstruction.

Semi-supervised ROCC method

JMAP Estimation

By using the JMAP estimation, the unknown variables are estimated iteratively. The maximization of the posterior distribution is transferred into the minimization of the negative logarithm of the posterior distribution, considering the exponential form of the distribution :

where

(6.37)

By estimating the unknown variables alternately, we estimate the unknowns by minimizing their corresponding criterions :

By calculating the zero points of the derivation of each criterion, we get the updating rules

-Simulation results

Original QR TV ROCC FIGURE 6.14 -The reconstructed Shepp Logan phantom of size 256 3 with 36 projections and SNR=40dB, by using different reconstruction methods. Red curves represent the profile at the position of the blue line.

and SNR=40dB is used. As shown in the results, by using the QR and TV, the contours of the reconstructed phantom are blurry, and in the object there exist some circle artifacts. In the reconstructed phantom by using ROCC method, the contours are clearer than the other state-ofthe-art methods.

Figure 6.15 and Figure 6.16 present a profile of the middle slice of the reconstructed Shepp Logan phantom. By comparing the ROCC method with the state-of-the-art methods, we can see that the shape of the reconstructed phantom is closer to the original one at the contour areas by using the ROCC method. When the number of projections is insufficient, ROCC method stays more robust than the other two state-of-the-art methods.

In Figure 6.17 and Figure 6.18, we compared the ROCC method with QR and TV method in terms of the RMSE of reconstructed object. The evaluation are compared with projection sinograms with different number of projections. As we can see in the figures, the RMSE is smaller when the number of projections is larger. When there are sufficient number of projections, both these three methods perform very well. When there are insufficient number of projections, the ROCC method stays more robust than QR and TV methods. In the cases where there is a high SNR, the ROCC method stays robust when there are more than 36 projections.

Comparison of computation time

The computation time is a crucial problem in the real applications, especially for the iterative methods for the three dimensional objects. In this section, we show the comparison of the computation time by using the TV, HHBM and ROCC methods.

In Table 6 First of all we consider a simplified model, in which the variance of noise and ξ are scalars, so we have :

With the VBA method, the posterior law is approximated by a separable law

The approximate distribution q (f , z, v z , v , v ξ ) is obtained via the minimization of the Kullback Leibler divergence.

df dz dv z dv dv ξ .

Now we want to minimize the term KL (q : p), for a hierarchical problem, we will optimize the parameters alternatively.

(zm)q 3 (vz)q 4 (v )q 5( v ξ ) (A.4)

where the defenition of the notation u

In the following expressions, the notation C represents a constant.

Where C is a polynome which has no relation with f j , so that we have :

where

We have :

From the previous computation we obtain :

where

From the formulation we can see that z j belong to a Gaussian distribution with mean and variance equal to :

A.3 Computation of q 3j v z j From Eq.(A.5), we have :

From the definition of variance we have z 2 j q 2j (z j ) = Var(z j ) + z j 2 q 2j (z j ) = µ 2 z j + σ z j , so that :

Computation of q 4 (v )

A.4 Computation of q 4 (v )

From Eq.(A.6) we have :

The computation begin with :

By using the property that : T r (AB) = T r (BA), we obtain :

Finally we find that v belong to an Inverse Gamma distribution with parameters equal to :

where

A.5 Computation of q 5 (v ξ )

From Eq.(A.7), we obtain :

By using the same method as for the calculation of g -Hf 2 2 q 1( f ) , we have :

And finally we have that v ξ belongs to an Inverse Gamma distribution with the parameters :

where

A.6 The iterative optimization by using the VBA algorithm

From the previous calculation, we obtain the updating formulations for all the parameters of the approximated distribution of each variable :

B

My contributions

This Appendix shows my main contributions during the thesis. The publications and the details are listed in Table B.1.

My work begins from the use of the basic Bayesian methods in CT reconstruction. 1 [WGMD15] used the Bayesian method in 2D case with a stationary noise. Several conventional prior models are used and compared. In this paper, the work mainly concentrate on the use of the Bayesian method on the X-ray CT reconstruction problem. 2D simulations are realized by using Bayesian method.

In 2 [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF], a hierarchical prior model is proposed to account for the sparse structure of the Haar Transformation coefficient of the image. We call this method the Hierarchical structured Haar based Bayesian Method (HHBM). In this method, the transformation coefficient is considered as a hidden variable, and is estimated simultaneously during reconstruction. The noise is modeled by a stationary Gaussian model with unknown variance, which is modeled by an Inverse Gamma distribution. The JMAP and the VBA methods are used for the optimization of the variables. The VBA method is realized in the 2D simulation with stationary noise model, however it is still computationally expensive.

In 3 [START_REF] Wang | Bayesian X-ray computed tomography using a three-level hierarchical prior model[END_REF], the HHBM method is developed, by adding an additive uncertainty variable in the Haar Transformation prior model : f = Hz + ξ. The simulation for 2D CT reconstruction is done. In this paper, a non-stationary noise model is considered. The variance of the noise of the forward projection model and the variance of the uncertainties of the forward transformation model are considered as unfixed parameters, and are modeled by Inverse Gamma distributions. The details and extensions of this conference paper is presented in the published journal article 4 [START_REF] Wang | X-ray Computed Tomography using a sparsity enforcing prior model based on haar transformation in a Bayesian framework[END_REF]. Then, in 5 [START_REF] Wang | 3D Xray Computed Tomography reconstruction using sparsity enforcing hierarchical model based on Haar transformation[END_REF], the HHBM method is developed for 3D CT reconstruction, in which the ASTRA toolbox is used in order to accelerate the computation. A French version of this article is presented in 6 [WGMD17].

In the above published papers, the discrete Haar transformation coefficient is considered because it is computationally cheap and is suitable for the strict piecewise-constant or piecewisecontinuous objects. However, sometimes the object is not strictly piecewise-constant, and some

C

Author's publications during 3 years of PhD

The list of my publications can be found online : https ://cv.archives-ouvertes.fr/li-wang-l2s.

Articles published in journals :

• In both domains, there is a need to reduce the number of projections. In some cases we may also be limited in angles. Consequently, we are always in the situation of ill-posed inverse problems. What's more, the measured data is always contaminated by noise. Facing these problems, the role of the probabilistic methods and the prior modelling becomes crucial.

For prior modelling, in particular in NDT applications, objects under examination are piecewise continuous, composing of several homogeneous materials with several continuous blocs separated by some discontinuities and contours. The focus of this thesis on the reconstruction of the piecewise continuous or constant, or more generally piecewise homogeneous objects. In summary two main methods are proposed in the context of the Bayesian inference.

The first method consists in reconstructing the object while enforcing the sparsity of the discrete Haar transformation coefficients of the object. A hierarchical Bayesian model is proposed. In this method, the unknown variables and parameters are estimated and the hyperparameters are initialized according to the definition of prior models.

The second method reconstructs objects while the contours are estimated simultaneously. The piecewise continuous object is modeled by a non-homogeneous Markovian model, which depends on the Laplacian of the object, while the Laplacian also depends on the estimation of the object. In this method, a semisupervised system model is also achieved, with parameters estimated automatically.

Both methods are adapted to the 3D big data size reconstructions, in which the GPU processor is used to accelerate the computation. The methods are validated with both simulated and real data, and are compared with several conventional stateof-the-art methods.