Imagerie et transformée de Radon

Reconstruction d’images et problèmes inverses

Pour cette partie nous allons nous placer dans un cadre classique de l’étude des problèmes inverses linéaires, c’est à dire dans le cadre des espaces de Hilbert où nous définissons deux ensembles F et G et un opérateur linéaire R : F → G tel que : Rf = g, (I.1) avec f ∈ F la fonction à reconstruire et g ∈ G le jeux de données que nous possédons. Imagerie et transformée de Radon 8 I. Imagerie et transformée de Radon I.1.1 Définition d’un problème mal-posé Un problème est dit bien-posé si le problème vérifie trois conditions : – l’existence, le problème admet au moins une solution pour tout ensemble de données. ∀g ∈ G ∃f : Rf = g (I.2) – l’unicité, la solution est unique, i.e. le noyau de R est : Ker(R) = 0 (I.3) – la stabilité, la solution dépend continûment des données, i.e. un faible écart sur les données ne doit pas entrainer une solution radicalement différente. On dit qu’un problème est mal-posé si au moins l’une de ces conditions n’est pas vérifiée. En général, les problèmes directs sont bien posés et les problèmes inverses sont souvent mal posés. En reconstruction d’image, le problème est quasiment systématiquement malposé ne serait ce qu’à cause de la limitation du jeu de données (fenêtres d’exploration limitées, durée de l’acquisition, nombre de prise de vue limitée, etc). L’étude de tels problèmes mal posés a conduit au développement de nombreuses méthodes de régularisation. 

Techniques de reconstruction d’images et de régularisation

De façon générale en problème inverse, nous pouvons distinguer deux approches pour la reconstruction d’image et donc pour la résolution du problème inverse lié : dans la première nous travaillons directement dans le cadre fonctionnel et en continu et nous étudions avant tout les fonctions et les opérateurs, ce sont les méthodes analytiques sur lesquels nous reviendrons plus en détail par la suite dans le cadre de la transformée de Radon ; dans la deuxième, nous nous plaçons dès le départ dans un cadre discret pour les équations dont je ferai ici un bref aperçu. De manière générale, toutes les méthodes de reconstruction ont besoin d’une régularisation. La régularisation d’un problème mal-posé consiste à le transformer en un problème bien posé, c’est-à dire à définir une solution unique pour toutes les mesures possibles dans l’espace des observations et à en assurer la stabilité vis-à-vis des erreurs sur ces mesures[8, 68]. C’est pourquoi les méthodes de reconstruction se basent ou incorporent une méthode de régularisation. 

Méthodes algébriques

Bien que ce manuscrit de thèse s’orientera avant tout sur les résolutions analytiques des différents problèmes, il existe d’autres façons de résoudre les problèmes inverses liés aux transformées de Radon généralisées et de façon plus générale aux problèmes inverses. Les méthodes algébriques sont basées sur l’introduction d’un système linéaire caractérisant le problème et de le résoudre. Cela permet de prendre en compte initialement la discrétisation des formules et de prendre en compte de nombreux phénomènes en amont assez facilement (atténuation, caractéristiques physiques, erreurs de mesure, géométrie du phénomène). Nous retrouvons donc le système linéaire g = Rf en posant une matrice de projection R qui est un ensemble de fonctionnelles linéaires traduisant le problème et transformant une carte de l’image f (ce qui est à reconstruire) en des mesures g. I.1. Reconstruction d’images et problèmes inverses 9 Les méthodes déterministes correspondent à la résolution du système par l’inversion directe de la matrice de projection R. Pour cela nous pourrons par exemple travailler le modèle d’inversion généralisée. Le système à résoudre sera alors R t g = R tRf. Si le rang de la matrice (R tR) de taille n × n est bien de n, alors nous aurons bien l’existence de la solution inverse généralisé : ˆf = (R tR) −1R t g (I.4) qui est unique et qui minimise l’erreur quadratique (g − Rf) t (g − Rf). R ∗ = (R tR) −1R t est l’inverse généralisé de R. Le problème restant étant la stabilité. Pour étudier cette stabilité, nous étudierons la décomposition de la matrice R en valeurs singulières où la source de l’instabilité se retrouvera pour les faibles valeurs singulières qui amplifieront le bruit. Il faudra alors régulariser pour obtenir un meilleur conditionnement de la matrice (R tR). Mais l’une des méthodes les plus efficaces pour ce faire reste la technique de la décomposition en valeurs singulières (SVD). Les valeurs singulières σi d’une matrice R sont les racines carrées des valeurs propres de la matrice carrée R tR. Les valeurs propres sont liées aux directions invariantes par la transformation R, tandis que les valeurs singulières contiennent la métrique sur cette transformation. Cette méthode permet de décomposer une matrice rectangulaire Rm×n (si le rang de la matrice R est rang(R) = min(n, m)) dont les coefficients appartiennent au corps K, où K = R ou K = C, de dimension m par n en un produit de matrice : Rm×n = Um×mSm×nV t n×n , (I.5) où Um×m est une matrice unitaire m × m sur K qui contient un ensemble de vecteurs de base orthonormés de Km, dits « de sortie », V t n×n est la matrice transposée de Vn×n, matrice unitaire n × n sur K qui contient un ensemble de vecteurs de base orthonormés de Kn , dits « d’entrée » ou « d’analyse » et Sm×n est la matrice diagonale qui contient les valeurs singulières de la matrice Rm×n. Néanmoins si le rang de la matrice R est rang(R) = k < min(n, m), la décomposition s’écrit alors : Rm×n = Um×kSk×kV t k×n , (I.6) A partir de la décomposition SDV, on peut obtenir la matrice pseudo-inverse de R : R −1 = V S−1U t . (I.7) Cette décomposition permet aussi de calculer le conditionnement du système (rapport entre le minimum et maximum des valeurs singulières) et si le système est trop mal conditionné, on pourra régulariser le problème par la troncature des valeurs singulières trop petites, ce qui reste l’un des principaux avantages de cette méthode. Dans le cadre de ces méthodes, pour effectuer une régularisation du système, nous pouvons aussi ajouter des contraintes à ce dernier (l’exemple classique en imagerie est la positivité) ou effectuer une régularisation avec a priori de douceur. On cherche ici la solution approchée bf qui minimise la fonctionnelle J du type : minJ(f) = ||Rf − g||2 + λΓ(f) (I.8) dans laquelle Γ(f) représente une contrainte de douceur sur f et λ est un multiplicateur de Lagrange. Par exemple, nous pouvons prendre le cas simple Γ(f) = ||f||2 qui est un 10 I. Imagerie et transformée de Radon cas particulier de la méthode de Tikhonov-Phillips avec Γ(f) = Z X k pk(t)| ∂ kf ∂tk | 2dt où k et les fonctions pk positives, continues et continûment dérivables restent à choisir. Les méthodes itératives quant à elles, permettent de trouver une solution par une succession d’estimations. Pour chaque estimation de l’objet, les projections correspondantes sont calculées et comparées aux projections mesurées, et le résultat de la comparaison est utilisé pour modifier l’estimation. La plus connue de ces méthodes est ART (algebraic reconstruction method) proposée par Gordon [27], mais il en existe d’autres comme l’algorithme SIRT par exemple. L’arrêt d’une méthode itérative en choisissant judicieusement le nombre d’itérations et la variation du poids de chaque correction permettent de régulariser la méthode. 

Méthodes statistiques

Les approches statistiques se basent sur un modèle direct qui considère les mesures de g du flux de photons arrivant sur le détecteur dépendant de phénomènes stochastiques (bruit poissonnien, bruit gaussien. . . ). Le problème de reconstruction devient l’estimation probabiliste des valeurs de f. Nous résolvons généralement ce problème d’estimation par des méthodes de maximum de vraisemblance [49] basées par exemple sur l’algorithme Expectation Maximisation (EM) [69]. La majeure partie de ces méthodes incorpore une régularisation bayésienne [24,29], servant à introduire une modélisation a priori de l’objet. Le principe de l’approche stochastique est de construire des distributions de probabilité pour modéliser le type de solution attendue et le processus de formation des projections. Nous notons p(f) le modèle a priori de l’objet, et p(g|f) le modèle du processus de mesure qui tient compte du bruit. L’application du théorème de Bayes fournit alors la loi a posteriori p(f|g) (qui nous permettra ’dobtenir notre solution) donnée par : p(f|g) = p(g|f)p(f) p(g) . (I.9) Parmi les différents modèles de l’objet introduit par le terme de probabilité p(f), nous trouvons notamment – le modèle gaussien – le modèle de maximum d’entropie – le modèle markovien.