Keywords: problème inverse, BAYES, segmentation, hyper-paramètre, déconvolution Inversion for textured images : unsupervised myopic deconvolution, model selection, deconvolution-segmentation Inverse problems, BAYES, segmentation, Hyper-parameters, deconvolution

This thesis is addressing a series of inverse problems of major importance in the field of image processing (image segmentation, model choice, parameter estimation, deconvolution) in the context of textured images. In all of the aforementioned problems the observations are indirect, i.e., the textured images are affected by a blur and by noise. The contributions of this work belong to three main classes: modeling, methodological and algorithmic.

From the modeling standpoint, the contribution consists in the development of a new non-Gaussian model for textures. The Fourier coefficients of the textured images are modeled by a Scale Mixture of Gaussians Random Field. The Power Spectral Density of the texture has a parametric form, driven by a set of parameters that encode the texture characteristics.

The methodological contribution is threefold and consists in solving three image processing problems that have not been tackled so far in the context of indirect observations of textured images. All the proposed methods are Bayesian and are based on the exploiting the information encoded in the a posteriori law. The first method that is proposed is devoted to the myopic deconvolution of a textured image and the estimation of its parameters. The second method achieves joint model selection and model parameters estimation from an indirect observation of a textured image. Finally, the third method addresses the problem of joint deconvolution and segmentation of an image composed of several textured regions, while estimating at the same time the parameters of each constituent texture.

Last, but not least, the algorithmic contribution is represented by the development of a new efficient version of the Metropolis Hastings algorithm, with a directional component of the proposal function based on the "Newton direction" and the Fisher information matrix. This particular directional component allows for an efficient exploration of the parameter space and, consequently, increases the convergence speed of the algorithm.

To summarize, this work presents a series of methods to solve three image processing problems in the context of blurry and noisy textured images. Moreover, we present two connected contributions, one regarding the texture models and one meant to enhance the performances of the samplers employed for all of the three methods.

Titre : Inversion pour image texturée : déconvolution myope non supervisée, choix de modèles, déconvolution-segmentation

Résumé : Ce travail est dédié à la résolution de plusieurs problèmes de grand intérêt en traitement d'images : segmentation, choix de modèle et estimation de paramètres, pour le cas spécifique d'images texturées indirectement observées (convoluées et bruitées). Dans ce contexte, les contributions de cette thèse portent sur trois plans différents : modéle, méthode et algorithmique.

Du point de vue modélisation de la texture, un nouveau modèle non-gaussien est proposé. Ce modèle est défini dans le domaine de Fourier et consiste en un mélange de Gaussiennes avec une Densité Spectrale de Puissance paramétrique.

Du point de vueméthodologique, la contribution est triple -troisméthodes Bayésiennes pour résoudre de manière :

-optimale -non-supervisée -des problèmes inverses en imagerie dans le contexte d'images texturées ndirectement observées, problèmes pas abordés dans la littérature jusqu'à présent. Plus spécifiquement, 1. la première méthode réalise la déconvolution myope non-supervisée et l'estimation des paramètres de la texture, 2. la deuxième méthode est dédiée à la déconvolution non-supervisée, le choix de modèle et l'estimation des paramètres de la texture et, finalement, 3. la troisième méthode déconvolue et segmente une image composée de plusieurs régions texturées, en estimant au même temps les hyperparamètres (niveau du signal et niveau du bruit) et les paramètres de chaque texture. La contribution sur le plan algorithmique est représentée par une nouvelle version rapide de l'algorithme Metropolis-Hastings. Cet algorithme est basé sur une loi de proposition directionnelle contenant le terme de la "direction de Newton". Ce terme permet une exploration rapide et efficace de l'espace des paramètres et, de ce fait, accélère la convergence.

Pȃrint ¸ilor mei, Paula si Cornel Surorii mele, Laura I warmly thank my colleagues in the Signal and Image Processing Team at IMS for the good moments we have shared in the 5 years I have been around. I have seen people come and go, but I always took pleasure in the company of Valerie, Céline, Jean-Pierre, Christian, Guillaume F., Pascal and Marc. Thank you Audrey for always being so positive and keeping us up to date with who's dating who and thank you Lionel for not making public the secrets you have on all of us. I wish to also thank my office colleagues Ioana, Olivier and Guillaume R. for accepting me among them, even if I wasn't a "teledetectionist". Special thanks to Olivier for helping me clean up the mess I often made by spilling my coffee and for not being afraid of being seated so close to me. And to Julien who made sure our lunch breaks were never too short and who took his coffee break with us without ever having coffee. I also thank Guillaume R. for only throwing my boots across the halway and never off the window. I am also very happy for having survived several troubling events in the team, such as the kidnapping of Julien's cactus by the "El Cartel" (safely returned to its owner in the end), the transformation of scientific posters into nudity revealing material, dangerous chair races and many more... I would also like to mention the support of a few very special friends who have been close to me throughout this time. Among them, Maria, Denisia, Roxana, Irina, Serena, Vincent, thank you for listening and encouraging me. I am particularly grateful to my parents for their unconditional love and support and for all the energy they have invested in me. My sister has also been extremely supportive and motivated to stay in contact with my work, so she read all my scientific papers (without having any image processing background). I am so fortunate to have you ! ! ! ! A very special thought goes to Guillaume, for all the patience and helping me keep my balance in my moments of doubt. Thank you for being there for me, listening and encouraging me.

Résumé

Ce travail est dédié à la résolution de plusieurs problèmes de grand intérêt en traitement d'images : segmentation, choix de modèle et estimation de paramètres, pour le cas spécifique d'images texturées indirectement observées (convoluées et bruitées). Dans ce contexte, les contributions de cette thèse portent sur trois plans différents : modéle, méthode et algorithmique.

Du point de vue modélisation de la texture, un nouveau modèle non-gaussien est proposé. Ce modèle est défini dans le domaine de Fourier et consiste en un mélange de Gaussiennes avec une Densité Spectrale de Puissance paramétrique.

Du point de vue méthodologique, la contribution est triple -trois méthodes Bayésiennes pour résoudre de manière :

-optimale -non-supervisée -des problèmes inverses en imagerie dans le contexte d'images texturées indirectement observées, problèmes pas abordés dans la littérature jusqu'à présent. Plus spécifiquement, 1. la première méthode réalise la déconvolution myope non-supervisée et l'estimation des paramètres de la texture, 2. la deuxième méthode est dédiée à la déconvolution non-supervisée, le choix de modèle et l'estimation des paramètres de la texture et, finalement, 3. la troisième méthode déconvolue et segmente une image composée de plusieurs régions texturées, en estimant au même temps les hyperparamètres (niveau du signal et niveau du bruit) et les paramètres de chaque texture. La contribution sur le plan algorithmique est représentée par une nouvelle version rapide de l'algorithme Metropolis-Hastings. Cet algorithme est basé sur une loi de proposition directionnelle contenant le terme de la "direction de Newton". Ce terme permet une exploration rapide et efficace de l'espace des paramètres et, de ce fait, accélère la convergence.

Résumé en franc ¸ais du contenu du mémoire

Contexte et formulation du problème L'énorme quantité de données numériques produites actuellement peuvent être utilisées dans des domaines très variés, néanmoins, leur très grand volume rend impossible leur traitement par des spécialistes. Par conséquent, il y a un besoin réel de méthodes de traitement automatique de données pour le contrôle de la qualité en industrie, l'aide au diagnostique dans des applications biomédicales, systèmes de vidéo-surveillance et nombreuses autres applications. Dans ce contexte, le traitement d'images est devenu une axe centrale de recherche scientifique.

Par ailleurs, la majorité des images naturelles comportent de la texture. En général, la texture contient une grande quantité d'information concernant la scène. Pour cette raison, dans le domaine du traitement d'images, les problématiques liées à la texture ont une importance particulière. Par conséquent, des nombreuses méthodes de modélisation, synthèse et analyse de textures ont été développées. Pour donner juste quelques exemples d'applications, ces méthodes sont utilisées dans des domaines comme :

• l'industrie : pour l'analyse et l'inspection de la structure de matériaux, la détection des défauts, l'analyse des modèles des tissus, le contrôle de la production ;

• imagerie biomédicale : segmentation d'organes ou de vaisseaux, l'analyse structurale des tissus afin de détecter des anomalies pouvant indiquer la présence de maladies, la détection de tumeurs, l'estimation du mouvement pour étudier le cycle cardiaque ou respiratoire ou la structure de différentes structures anatomiques ;

• traitement d'images satellite : segmentation d'images pour des raisons cartographiques, classification de régions (urbaine, rurale, exploitation agricole, environnement naturel), classification de cultures (vigne, forêt), estimation et prédiction de la productivité d'une culture, évaluation de dommages (tempête, incendie) ;

• investigations géologiques : caractérisation des structures géologiques, à partir de différentes techniques d'imagerie.

Les applications dédiées à la texture consistent en la définition d'un modèle pour la texture et l'apprentissage des caractéristiques de ce modèle à partir d'un certain nombre de textures de référence. Cette information peut être utilisée pour réaliser deux tâches principales :

• analyse de la texture : i.e., l'extraction de l'information pertinente décrivant les principales caractéristiques de la texture. Cette analyse peut être basée sur l'interaction entre les pixels dans différents espaces de représentation (spatial, ondelettes), sur la distribution des fréquences dans le domaine de Fourier, sur l'orientation des

x Contents vecteurs ou tenseurs, sur la saillance, sur l'entropie, etc. . .. Plus spécifiquement, c ¸a revient à trouver des mesures qui décrivent l'image, appelées attributs, qui permettent d'identifier les textures similaires et de distinguer les textures différentes.

• synthèse de la texture : en utilisant modèle déterminé à partir des textures de référence, des nouvelles images ayant des caractéristiques similaires et appartenant à la même famille de textures peuvent être générées.

Un aspect très important est que pas tous les modèles sont génératifs, par exemple, les méthodes d'analyse de textures basées sur des attributs peuvent être utilisées seulement pour classifier ou comparer des textures et pas pour synthétiser une nouvelle texture.

Ce travail représente un étude sur des divers problèmes inverses en traitement d'images, appliqués au cas spécial de la texture. Les problèmes inverses qui sont traités sont la déconvolution, le débruitage, la segmentation, l'estimation des paramètres et le choix de modèle. Plus spécifiquement, trois méthodes ont été développées : I. une méthode de déconvolution myope, qui estime les hyperparamèters (est nonsupervisée), les paramètres de l'instrument et les paramètres de la texture, II. une méthode pour le choix de modèle de la forme de la matrice de covariance pour des images texturées indirectement observées, i.e., floutées et bruitées, III. une méthode pour la déconvolution-segmentation d'images texturées, qui estime au même temps le niveau du bruit et les paramètres des textures.

Tous ces problèmes sont très mal-posés, c'est-à-dire que les données observées ne suffisent pas pour obtenir une solution unique. Pour cela, de l'information supplémentaire est nécessaire. Une manière simple et naturelle est la formulation Bayésienne où l'information supplémentaire est incluse sous la forme des lois a priori. En fait, toutes les méthodes présentées dans ce manuscrit sont basées sur une formulation Bayésienne, où l'information contenue par les données est exploitée via le terme d'adéquation aux données (la vraisemblance), tandis que l'information a priori constitue le terme de régularisation.

Dans le cadre des approches Bayésiennes, toute l'information est encodée par la loi a posteriori. Des estimateurs sont définis à partir de cette loi, leur rôle étant d'exploiter l'information afin de trouver une solution du problème. En fonction de l'estimateur qui est utilisé, la solution peut avoir certaines propriétés d'optimalité. Par exemple, dans ce travail, la Moyenne A Posteriori est utilisée très souvent, sachant que cet estimateur minimise l'erreur quadratique moyenne. Pour calculer les estimés, la Moyenne A Posteriori utilise des échantillons de la loi a posteriori. Néanmoins, comme dans le cas de la plupart des problèmes non-triviales, cette loi a une forme compliquée et ne peut être échantillonnée directement. Par conséquent, des méthodes de Monte Carlo par Chaîne de Markov doivent être employées, plus spécifiquement l'échantillonnage de Gibbs et, pour les lois conditionnelle non-standard, l'algorithme de Metropolis-Hastings. Les problèmes d'inversion mentionnées précédemment sont traitées dans le contexte d'images texturées. De plus, le focus est mis sur une classe particulière de textures, modélisée par des champs aléatoires centrés et stationnaires, avec une Densité Spectrale de Puissance paramétrique. Le modèle de texture basé sur de champs aléatoires gaussiens, Contents xi qui représente le premier modèle exploité dans ce travail, possède des caractéristiques très avantageuses qui permettent son intégration dans le paradigme d'inversion. D'un autre coté, même si ce modèle est relativement versatile, il souffre de quelques limitations :

• il est basé seulement sur les statistiques d'ordre 2,

• les coefficients de Fourier de la texture sont indépendants. L'extension de ce modèle à une loi non-gaussienne n'est pas une tâche facile, car le but est d'augmenter sa capacité de représentation, tout en gardant sa complexité à un niveau qui permet son intégration dans notre formalisme de traitement d'images. Dès lors, il y a un compromis à faire entre la versatilité et la complexité du modèle, i.e., afin de pouvoir gérer facilement le modèle, certains limitations vont être gardés.

Contributions

Les principales contributions présentées dans ce manuscrit reflètent la dualité de ce travail. D'abord, le focus sur les images texturées a ouvert la question sur le besoin d'un modèle de texture qui est à la fois facilement gérable et qui a des capacités de représentation satisfaisantes. Ensuite, en utilisant ce modèle, de problèmes compliquées en traitement d'images ont été abordées dans le contexte des images texturées. Finalement, la résolution de ces problèmes de manière optimale a nécessitée l'utilisation de l'échantillonnage. En général, l'échantillonnage n'est pas un outil très souvent utilisé en traitement d'images, à cause de la grande dimension des problèmes. Néanmoins, grâce aux propriétés du modèle de texture proposé, l'échantillonnage n'est pas trop couteux dans ce cas. De plus, pour améliorer les performances de ces méthodes, nous avons développé un algorithme efficace d'échantillonnage pour des lois cible compliquées.

Ces contributions sont de nature assez différente. La première contribution est liée à la modélisation : Mo1. Proposition d'un modèle non-gaussien de texture, basé sur les moments d'ordre 2, et son intégration dans un problème compliqué de traitement d'images, les prochains trois contributions sont de nature méthodologique -développement et mise en oeuvre de méthodes pour :

Me1. la déconvolution myope, avec estimation des hyperparamètres, paramètres instrument et paramètres de la texture,

Contenu de la thèse

Chapitre 2 : Problèmes inverses, cadre bayésien et échantillonnage stochastique

Ce chapitre est dédié à l'introduction des notions qui seront récurrentes pendant ce travail, par exemple l'idée d'observations indirectes et les divers problèmes inverses qui peuvent être abordés dans ce contexte. Le système d'observation choisi introduit un flou et un bruit blanc gaussien additif. De plus, le flou est considéré comme ayant une forme paramétrique avec des paramètres η et la loi du bruit est pilotée par le paramètre γ n . Sachant que ces paramètres sont inconnus, ce travail va viser à résoudre le problème de la déconvolution dans un cas non-supervisé (estimer γ n ) et myope (estimer η).

Le problème complète qui peut se poser dans ce cas pour des images texturées est décrit dans la suite. Celui-ci combine les aspects non-supervisés et myopes, associés aux observations indirectes, avec la sélection du modèle de texture et la segmentation d'images texturées. Ce problème est fortement mal posé à cause du déficit d'information (le nombre d'inconnues est nettement plus grand que celui des observations). Ce problème inverse sera traité dans un cadre bayésien, où l'information contenue par les données est extraite via la loi de vraisemblance et la régularisation est réalisée en imposant des modèles pour les inconnues et des lois a priori sur les paramètres de ces lois. La loi a posteriori peut être exprimée à partir de ces ingrédients et des estimateurs pour les inconnues seront formulés. Vu que ces estimateurs ne peuvent pas être calculés de manière analytique, ils seront calculés numériquement en utilisant un échantillonneur de Gibbs, qui consiste à échantillonner de manière itérative les lois a posteriori conditionnelles.

Dans ce contexte, ce chapitre introduit la contribution algorithmique, un échantillonneur efficace du type Metropolis-Hastings pour des lois cible compliquées. Cet échantillonneur est basé sur une proposition directionnelle similaire à la direction de Newton. La spécificité de cette proposition est que la matrice Hessienne est remplacée par la matrice d'information de Fisher, afin d'éliminer certaines instabilités numériques. De plus, dans notre cas spécifique de loi pour les paramètres de la texture, ce principe nous permet de profiter des avantages d'une proposition directionnelle de deuxième ordre en calculant seulement les dérivées de premier ordre.

Chapitre 3 : Modélisation de la texture Cette thèse est focalisé sur le traitement d'images texturées. La plupart des travaux qui portent sur ce sujet sont basés sur le principe d'analyse de la texture, sous la forme de l'extraction d'attributs texturaux et la comparaison de ces attributs. L'approche choisie Contents xiii dans ce travail et d'utiliser des modèles génératifs, capables au même temps d'indiquer si une image est bien décrite par ce modèle et de produire de réalisations de textures ayant certaines caractéristiques. Un modèle souple et facile à intégrer dans le cadre de notre problématique est le modèle basé sur des Champs Aléatoires Gaussiens (CAG). Ce modèle repose sur l'hypothèse que les images sont stationnaires, ce qui implique, par approximation de Whittle, une structure circulante-bloc-circulante de la matrice de covariance. Cette structure particulière implique que la matrice est diagonalisable par Transformée de Fourier (TF) et que les coefficients de Fourier de l'image sont indépendants (décorrelés et gaussiens). Grâce à cette propriété, les coefficients de Fourier de l'image peuvent être traités en parallèle, par conséquent, de manière très efficace. De plus, dans ce travail, des formes paramétriques sont employés pour la Densité Spectrale de Puissance (DSP) des images. Ces formes sont pilotées par un nombre réduit de paramètres, ce qui implique une forte compressibilité du modèle.

La gaussianité et l'indépendance des coefficient de Fourier de l'image peuvent être considérés comme des limitations assez fortes du modèle. Néanmoins, ces propriétés jouent un rôle crucial dans la souplesse et la trac ¸abilité du modèle. Pour cette raison, l'effort d'augmenter la capacité de représentation du modèle va s'articuler autour du principe de garder, au moins en partie ces deux propriétés. Le nouveau modèle de texture proposé dans ce chapitre est basé sur l'idée de la gaussianité conditionnelle par rapport à un ensemble de variables auxiliaires de telle manière que marginalement par rapport à ces variables, la distribution ne sera plus gaussienne. Plus spécifiquement, un ensemble de variable auxiliaires sont introduites, une par pixel, et les coefficients de Fourier de la texture on une loi conditionnelle gaussienne et une loi marginale non-gaussienne. Ce modèle permet une extension des capacités de représentation du modèle (qui peut mieux décrire des texture plus stochastiques), grâce à une DSP enrichie, tout en gardant les avantages de l'indépendance et de la gaussianité conditionnelle des coefficients de Fourier.

Chapitre 4 : Déconvolution myope et estimation des paramètres d'images texturées

Le premier chapitre dédié à la résolution de problèmes inverses dans le contexte d'images texturées indirectement observées vise l'estimation des paramètres de l'image texturée. Cette estimation sera réalisée conjointement avec une déconvolution myope (estimation des paramètres instrument) et non-supervisée (estimation du niveau du bruit). Le formalisme bayésien offre un cadre unifié où tous ces paramètres et leur inter-dépendances sont représentés. Dans cette approche, l'information contenue par les données est exploitée via la loi de vraisemblance et la régularisation de ce problème mal-posé est réalisée à travers l'utilisation de modèles pour l'image x et la réponse instrument H. L'image x est modélisée par un CAG avec DSP paramétrique ayant un profil laplacien. Ce profil est piloté par l'ensemble de paramètres θ, contenant les deux fréquences centrales et les dispersions fréquentielles correspondantes dans l'espace des fréquences réduites. La réponse instrument est représentée par un filtre passe-bas, avec profil de Dirichlet, centré et isotrope, piloté par le paramètre de dispersion η.

Afin d'avoir une description complète du problème, des lois a priori pour les pa-xiv Contents ramètres γ n , θ et η seront choisies. Pour γ n une loi Gamma a été employée grâce à sa forme conjuguée par rapport à la loi de vraisemblance. En ce qui concerne les θ et η, à cause de la dépendance très compliquée, il n'existe pas de forme conjuguée. Par conséquent, une loi uniforme sur le domaine de définition de ces paramètres a été choisie.

A partir de tous ces ingrédients, la loi a posteriori peut être exprimée. L'information contenue par cette loi est exploitée en utilisant l'échantillonnage de Gibbs. Pour les lois conditionnelles non-standard des θ et η, des étapes de Metropolis-Hastings ont été intégrées dans l'échantillonneur de Gibbs. En utilisant les échantillons ainsi obtenus, des estimateurs EAP sont calculés pour chaque inconnue. Comme sous-produit, cette méthode fournit aussi une estimation de l'image originale.

Chapitre 5 : Sélection de modèle pour des images texturées indirectement observées

Ce chapitre porte sur le problème de sélection de modèle de texture pour des images floutées et bruitées. Plus exactement, les images texturées sont modélisées par des CAG ou MCAG avec des formes de PSD paramétriques et l'objectif est de sélectionner la forme de la DSP parmi un ensemble de modèles candidats. L'aspect données indirectes est traité de manière non-supervisée (le paramètre γ n est estimé), mais il ne s'agit plus d'une déconvolution myope, car la réponse instrument est considérée connue.

Ce problème est abordé dans un cadre bayésien, à travers une stratégie optimale du point de vue du risque moyen de classification. Cette stratégie est basée sur le calcul de l'évidence (autrement appelées la vraisemblance marginale). Néanmoins, l'évidence a une expression très compliquée et, par conséquent, ne peut pas être calculée de manière analytique. Pour cette raison, des méthodes numériques doivent être utilisées. Parmi les différentes options, nous avons employé la moyenne harmonique et l'approximation de Laplace-Metropolis. Il est démontré expérimentalement que les deux approximations sont équivalentes et fournissent les mêmes résultats de sélection. De plus, cette application intègre l'échantillonneur efficace Fisher Metropolis-Hastings pour la loi conditionnelle compliquée des paramètres des DSPs θ.

Chapitre 6 : Segmentation d'images texturées indirectement observées

Le problème de la segmentation d'images texturées floutées et bruitées est traité de manière non-supervisée dans un cadre bayésien. Les images sont composées par un nombre inconnu R de régions, chacune de ces régions contenant une texture d'une certaine classe k. Le nombre de classes K est connu. Chacune des K textures est modélisée pas un CAG avec DSP paramétrique ayant un profil connu, mais des paramètres inconnus.

La composition de l'image à partir des K textures est pilotée par un champs d'étiquettes cachées. Cet ensemble d'étiquettes est modélisé par un champs de Potts. La méthode est basée sur une approche probabiliste pour estimer conjointement les étiquettes, le niveau du bruit et les paramètres de chaque classe de texture. Les estimés sont obtenus en utilisant les échantillons de la loi a posteriori de la manière suivante : MAP marginalisé pour les

Conclusion et perspectives

Ce travail a présenté une gamme variée de problèmes inverses de grande importance dans le traitement d'images, appliqués dans le cas particulier d'images texturées. Dans ce contexte, nous avons proposé une série de contributions sur un plan algorithmique, modélisation et méthodologique.

La contribution algorithmique est représentée par le développement d'une version efficace de l'algorithme Metropolis-Hastings. La loi de proposition de ce nouvel échantillonneur est basée sur une composante inspirée de la direction de Newton, dans laquelle la matrice hessienne est remplacée par la matrice de Fisher. Par conséquent, les problèmes liés à l'inversion de la matrice hessienne sont évités et la composante directionnelle de la proposition a surement la direction de descente de gradient. En outre, dans le cas particulier de notre loi conditionnelle a posteriori, la loi de proposition s'écrit seulement à base de la première dérivée de la loi cible.

En ce qui concerne les perspectives de ce travail, cet algorithme peut être inclus dans une grande variété d'applications, dans tous les cas où des lois de probabilité compliquées doivent être échantillonnées. ********** La contribution liée à la modélisation consiste dans le développement d'un modèle pour l'analyse et la synthèse de textures, basé sur un Mélange de Champs Aléatoire Gaussiens avec Densité Spectrale de Puissance paramétrique. La forme paramétrique de la Densité Spectrale de Puissance et les valeurs de ses paramètres rassemblent les caractéristiques texturales et permettent la classification de la texture ou la synthèse de nouveaux échantillons de texture ayant les mêmes attributs que l'original. Ce modèle non-gaussien est construit en utilisant un ensemble de variables auxiliaires, de telle manière que la loi de l'image soit gaussienne, conditionnellement à ces variables, et non-gaussienne marginalement. Dans ce travail, afin d'assurer des bonne performances en terme de vitesse, les variables auxiliaires sont considérés indépendantes.

Les perspectives dans ce direction de recherche visent l'utilisation des priors corrélés pour les variables auxiliaires. Cela signifie que les variables auxiliaires ne seront plus indépendantes, résultant probablement dans l'augmentation des capacités de représentation de notre modèle de texture. Cependant, l'effet secondaire est l'inefficacité calculatoire de l'échantillonnage d'une telle texture. Une autre piste est l'utilisation de formes plus complexes pour la Densité Spectrale de Puissance et l'étude de la possibilité d'exploiter la phase des coefficients de Fourier afin de construire des modèles plus puissants. Ce travail peut être facilement enrichi en utilisant une forme plus compliquée pour le filtre convolutif. Une autre extension pourrait être l'utilisation de notre modèle de texture non-gaussien. *** La deuxième contribution méthodologique est une méthode de sélection de modèle basée sur le calcul de l'évidence dans le contexte d'observations indirectes d'images texturées. Comme dans le cas précédent, les images texturées sont affectées par un flou et par du bruit et le but est de sélectionner le modèle et les valeurs des paramètres de la Densité Spectrale de Puissance de la texture. Dans ce cas, le filtre convolutif est connu, i.e., nous ne positionnons pas dans un cadre déconvolution myope. Toutefois, le niveau du bruit n'est pas connu et doit être estimé. Les textures sont modélisées soit par des "Scale Mixture of Gaussian Random Fields", soit par des champs aléatoires gaussiens, et la méthode choisit le modèle qui a plus probablement généré la réalisation de la texture. La méthode est basée sur de l'échantillonnage et sa mise en oeuvre inclut l'échantillonneur efficace Fisher Metropolis-Hastings.

Dans la continuité de ce travail, les dictionnaire de formes pour la Densité Spectrale de Puissance pourrait être étendu. Une autre perspective est l'estimation du filtre convolutif. *** Une dernière contribution méthodologique est une méthode pour la déconvolutionsegmentation d'images texturées, indirectement observées. Comme précédemment, le filtre convolutif est connu et le niveau du bruit est estimé. Cette méthode arrive à fournir de très bons résultats pour un problème de difficulté significative. Dans le contexte où la segmentation de textures est en soi-même une tâche difficile, il n'y a pas de travaux dans la déconvolution-segmentation d'images texturées. Notre méthode utilise un modèle de Potts pour les étiquettes, avec un paramètre de température réglé à la main, et un modèle de champs aléatoire gaussien pour les textures. Pour chaque texture, le modèle de sa Densité Spectrale de Puissance est connu, tandis que ses paramètres doivent être estimés.

Hormis les défis théoriques significatives de ce problème, sa mise en oeuvre nous a confrontés avec la tâche difficile d'échantillonner les étiquettes et les pixels de manière efficace, afin d'assurer un coût calculatoire acceptable. Cette contrainte nous a guidé graduellement vers cette version finale de formulation du problème, présentée dans ce manuscrit. Néanmoins, cette version est le résultat d'une série d'alternatives que nous avons explorées. Parmi ces formulations, celle que nous avons choisie est la seule qui nous a permis au même temps de représenter correctement les dépendances entre les variables et d'obtenir de bonnes performances algorithmiques.

Contents xvii

Ce sujet peut être encore développé en intégrant l'estimation du paramètre de température. Dans ce cas, la méthode sera capable de s'adapter automatiquement à la topologie de l'image traitée. Une autre perspective serait d'estimer le filtre convolutif.

Abstract

This thesis is addressing a series of inverse problems of major importance in the field of image processing (image segmentation, model choice, parameter estimation, deconvolution) in the context of textured images. In all of the aforementioned problems the observations are indirect, i.e., the textured images are affected by a blur and by noise. The contributions of this work belong to three main classes: modeling, methodological and algorithmic.

From the modeling standpoint, the contribution consists in the development of a new non-Gaussian model for textures. The Fourier coefficients of the textured images are modeled by a Scale Mixture of Gaussians Random Field. The Power Spectral Density of the texture has a parametric form, driven by a set of parameters that encode the texture characteristics.

The methodological contribution is threefold and consists in solving three image processing problems that have not been tackled so far in the context of indirect observations of textured images. All the proposed methods are Bayesian and are based on the exploiting the information encoded in the a posteriori law. The first method that is proposed is devoted to the myopic deconvolution of a textured image and the estimation of its parameters. The second method achieves joint model selection and model parameters estimation from an indirect observation of a textured image. Finally, the third method addresses the problem of joint deconvolution and segmentation of an image composed of several textured regions, while estimating at the same time the parameters of each constituent texture.

Introduction 1.1 Context and Problem statement

Humanity currently produces an exponentially increasing amount of information under the form of digital data, a significant part of which consists in images and videos. This huge amount of information could be exploited in various domains, nevertheless, the huge volume makes the processing by human specialists prohibitive. Consequently, while this data has allowed us to document and better understand our world, it has also given rise to a need for automatic processing to perform relevant data selection and provide preliminary results. We are thus witnessing an increasing need to automatize the processing in fields such as industrial quality control, medical diagnosis, image analysis and video surveillance systems. In this context, image processing has become a central axis of scientific research.

Moreover, most natural images contain texture, which generally encodes a great amount of information concerning the scene. For this reason, in the field of image processing, the texture related topics represent a field of significant importance. The considerable attention that has been payed to these topics has resulted in the development of a multitude of modeling, synthesis and analysis methods. To list only a few of the possible applications, let us refer to the fields of:

• industry: material structure analysis and inspection, flaws detection, fabric pattern analysis, production control;

• biomedical imaging: organ or vessel segmentation, tissue analysis for determining whether there are abnormalities that may indicate the presence of a disease, tumor detection, motion estimation to study the cardiac or the breathing cycle;

• satellite image processing: image segmentation for cartographic purposes, region recognition (urban, rural, crop, natural), crop recognition (vineyards, forests), crop production estimation and predictions, damage evaluation (what percentage of a crop / forest has been damaged by a storm / fire).

• geological investigations: characterizing the geological structures based on various scanning techniques.

The texture related applications consist in defining a model for the texture and acquiring information regarding that model based on observing a certain number of reference textures. Then, this information can be used to perform the two main tasks:

• texture analysis: or description, i.e., the extraction of pertinent information describing the main textural characteristics. It can be based on pixel interaction in various representation spaces (spatial, wavelet), on the frequencies distribution in the Fourier domain, on orientation vectors or tensors, on saliency, on entropy, etc. . . . More specifically, it consists in finding measures that best describe the texture, called features, that permit us to identify similar textures and to distinguish the textures that are different.

• texture synthesis: using the model that has been determined using the reference textures, new images can be generated that have similar characteristics and that belong to the same family of textures, as it is defined by the model.

An important aspect is that not all the models are generative, i.e., some of them (the feature based texture analysis methods) can only be used for comparison and classification purposes and cannot be used to synthesize a new texture.

This thesis presents a study on inverse problems in image processing applied to the special case of textures. The inverse problems we are dealing with are deconvolution, denoising, segmentation, parameter estimation and model choice. More specifically we will present: I. a myopic deconvolution method, that also estimates the hyperparameters, the instrument parameters and the texture parameters, II. a model choice method for the form of the covariance matrix for an indirectly observed texture, i.e., a blurry and noisy version of the texture, III. a method for the joint deconvolution-segmentation of a textured region, that also estimates the noise level and the texture parameters.

All of the previously mentioned problems are very ill-posed, meaning that the data alone does not contain sufficient information as to provide a unique solution. For this reason, supplementary information must be provided and this can be done in several manners. A simple and natural way to do this is to use a Bayesian framework and embed all the extra information in the a priori laws. In fact, all the methods presented in the following rely on a Bayesian formulation, where the information contained by the observations is exploited via the data adequacy term (the likelihood), while the a priori information forms the regularization term.

In the previously described Bayesian formulation for our problems, all the available information is embedded in the a posteriori law. Estimators are then defined based on this law, their role being to exploit the information in order to yield a solution for the problem. Depending on the estimator that is being used, the solution can have some optimality properties. For instance, in this work, the Minimum Mean Square Error estimator is predominant, which, as specified by its name, is the estimator that minimizes the mean square error. In order to compute the estimates, the Minimum Mean Square Error relies on samples of the a posteriori law. However, as it is the case in most non-trivial problems, this law has a complicated form and thus cannot be sampled directly. To this end, Monte Carlo Markov Chain methods will be employed, more specifically, Gibbs sampling and, for the parameters with non-standard conditional laws, the Metropolis-Hastings sampler.

The inversion problems previously described are applied to textured images. Moreover, we focus on a particular type of textures, modeled by stationary, zero-mean Random Fields with structured covariance, i.e., a parametric form for the Power Spectral Density. The texture model based on Gaussian Random Fields, which represents the first model we have explored, presented a series of extremely advantageous features, enabling its integration in the inversion paradigm. Although this model is rather versatile and yields a large variety of textures, it suffers from several limitations:

• it is only based on the first and second order statistics,

• the texture Fourier coefficients are independent.

Extending the model to a non-Gaussian law, in a manner that increases its representation capacity, but at the same time keeps its complexity at a level that permits its integration in the larger image processing formalism, is not an easy task. Consequently, there is a compromise to be made between model versatility and complexity and this means that in order to keep the model handleable, some of the aforementioned limitations will be preserved. Figure 1.1 shows several texture realizations for a model with a parametric Power Spectral Density.

Contributions

The main contributions presented in this manuscript reflect the duality of our work. On the one hand, our focus on textured images has confronted us with the need for a texture model that is at the same time tractable and has strong representation capabilities. On the other hand, using on the aforementioned texture model, we were able to tackle image processing problems that are very complicated in the context of textured images. Finally, solving these problems in an optimal manner has required the use of sampling. Generally, sampling in the context of image processing is not a popular tool, due to the large dimension of the unknown. Nevertheless, due to the properties of our texture model, sampling is not prohibitively costly in our case. To further enhance the performances of our methods, we have devised an efficient sampling technique for complicated target laws.

As it can easily be noticed, our contributions are of relatively different nature. To be more exact, the first contribution is model-related:

Mo1. Development of a non-Gaussian texture model, based on the first and second order moments, and its integration in a complex image processing problem, the following three contributions are of methodological nature:

Me1. Devising and implementing a myopic deconvolution method with joint estimation of the hyperparameters, instrument parameters and texture parameters, Me2. Devising and implementing a method for model selection for the Power Spectral Density of a textured image, using an indirect observation, Me3. Devising and implementing a deconvolution-segmentation method for textured images with joint estimation of the hyperparameters and texture parameters, while the last contribution is algorithmic:

A1. Development of and efficient version of the Metropolis-Hastings sampler, based on the Fisher Information Matrix and its integration in our algorithms.

Manuscript Structure

This manuscript is organized as follows: This Chapter also provides a description of the numerical algorithms employed for sampling and the various versions of these samplers, along with our first contribution, the efficient sampler Fisher Metropolis-Hastings.

• Chapter 3 -Texture Modeling presents our new non-Gaussian texture model based on a Scale Mixture of Gaussian Random Fields (SMGRF). The Power Spectral Density of these Random Fields is parametric and this Chapter explains the significance of the parameters and provides priors. This represents in fact the second contribution of our work.

• Chapter 4 -Unsupervised Myopic Deconvolution of a Textured Image focuses on the first methodological contribution. The method presented here relies on a fully parametric formulation of the observation system and of the image model. Consequently, the unsupervised myopic deconvolution problem becomes a threefold parameter estimation problem: observation system parameters, image model parameters and hyperparameters.

• Chapter 5 -Model Choice for the Law and the PSD of a Textured Image is devoted to our second methodological contribution, a method to select the Power Spectral Density model for a texture among a series of candidates, based on blurred and noisy observations.

• Chapter 6 -Deconvolution Segmentation for Textured Images presents the third methodological contribution: a method for textured image segmentation from blurred and noisy data. This Bayesian method assigns a prior to the labels and provides a global approach for estimating the unknowns (the labels, the original image, the noise precision).

• Chapter 7 -Conclusion and Perspectives lists the most important aspects that have been clarified throughout this thesis, the most interesting remarks and the main questions that remain unanswered concerning the addressed problems.

In support of the technical information comprised in the aforementioned Chapters, we have provided detailed information regarding some connected topics. In order to preserve the natural flow of our presentation, these topics have made the object of a series of Appendices:

• Appendix A -Fisher information for indirectly observed SMGRF textures gives a detailed description on how to compute the Fisher information regarding the various parameters of the observation system and of the non-Gaussian texture model.

• Appendix B -Optimal Bayesian Estimation analyzes in parallel the optimality property of the Minimum Mean Squared Error estimator and of the evidence-based classifier, analysis relying on the Bayesian risk minimization.

• Appendix C -Potts Model explicits the form of the Potts prior for modeling the pixel interaction in image processing.

Inverse Problems in Image Processing

Inverse problem is the denomination given to any attempt to infer on the inputs of a physical system, starting from the data observed at its output.

I O Transformation

Some of the most commonly encountered inverse problems in signal and image processing are deconvolution, denoising, parameter estimation, segmentation, super-resolution, optical flow estimation, source detection, particle image velocity, model choice.

In general, the unknowns are obtained as the solution of a system with N equations, where N is the number of observations, O. Most often than not, they are ill-posed problems, since the number of data is smaller than the number of unknowns and thus it is not Chapter 2. Inverse Problems, Bayesian Framework and Stochastic Sampling straightforward to determine a unique solution. In order to overcome this impasse, supplementary constraints must be imposed on the solutions. This is done by introducing information concerning the structure of the expected solution. This information is called regularization information, it is acquired independently on the observation process and it represents the knowledge we possess concerning the physical quantity. This information can introduce for instance a smoothness or a sparsity constraint, a restriction of the range of possible values or favor certain solutions.

The Bayesian framework represents a natural setting for such problems. Its principle consists in probabilizing the unknowns and defining a likelihood law and a priori laws. The likelihood is formulated based on the direct physical model of the system and represents the probability density of the data realizations, conditionally on the unknowns. The a priori laws contain the regularization information under the form of probability distributions for the unknowns. The solution to the inverse problem can be found by maximizing a criterion based on two terms: a data term, corresponding to the likelihood law in the Bayesian setting, and the regularization term.

Among the numerous and varied aforementioned problems, this work focuses on parameters estimation, model choice and segmentation, all this in the context of indirect observations. The indirect observations mean that the data are a blurred and noisy version of the original image on which we should base the further inference to determine the quantities of interest. Consequently, throughout this work, deconvolution is a recurrent topic, bringing an additional layer of complexity to the problems. These problems are exclusively addressed in the advantageous Bayesian framework, where the aspect of indirect data can naturally be merged with the actual problem of estimating the underlying parameters, selecting the most probable model or segmenting the original image.

Indirect Data

Deconvolution is a topic generating vivid interest in the signal processing community. It is encountered in most applications involving indirect observations, since there are no ideal observation systems. Let us represent the convolution based observation system as follows:

X + N Y h
This convolutive system is described by the equation:

Y = h * X + N (2.1)
where X, N and Y are the unobserved original image, the noise and the observations, respectively, and h is the convolution filter. The original images are of size N × N and P = N 2 is the total number of pixels. In the following, the mathematical developments 2.1. Inverse Problems in Image Processing are based on the alternative expression of (2.1):

y = Hx + n (2.2)
where the P × 1 size vectors x, n and y are obtained by lexicographically ordering the pixels of the X, N and Y matrices. Furthermore, H represents the P × P convolution matrix representing the point spread function (PSF) of the blurring filter.

In the context of such an observation system, various problems may be considered:

-for a known PSF (known H), if the eigenvalues of H are not tending to zero, i.e., there are no instabilities, the problem is well-posed and can be solved by applying the inverse filter. When dealing with noise corrupted data, even for known PSF, the problem becomes ill-posed. Nevertheless, it can be solved through Wiener filtering.

-for an unknown PSF (unknown H), the instrument characteristic is estimated along with the original image. This ill-posed problem can be solved by adding regularization information in order to constrain the solution. One possibility is to make no assumptions on the form of H, except for some regularity constraints, meaning that the PSF estimation implies the estimation of all its elements -blind deconvolution. The alternative option is to consider that the PSF has a parametric form (H η ), in which case its estimation comes down to estimating the parameters (η) driving this form -myopic / semi-blind deconvolution. The highlighted aspects describe the observation system considered in this work: data corrupted by noise and an unknown, parametric PSF, driven by the parameter set η.

In this context, we have focused on an additive, Gaussian model for the noise. This model is driven by the covariance R n (γ n ), which depends on the parameters γ n . Although the methods presented in the following are adapted for any type of additive Gaussian noise, for algorithmic simplicity, we are only considering the case of white noise. This is the assumption made in the vast majority of the deconvolution applications [MMK06,BMK08, BMK09, BWMK10, BMK11] and consists in setting R n = γ n I, where γ n is the inverse variance / precision parameter for the noise. In this case, based on the aforementioned assumptions about the noise, the likelihood law can be expressed as:

f (y|γ n , η, x) = (2π) -P/2 γ P/2 n exp - 1 2 γ n y -H η x 2 (2.3)
The operation of convolution has the special property of being written as a product in the Fourier domain. Moreover, in the case where the H η has a circulant form, it is diagonalizable by Discrete Fourier Transform (DFT). This assumption has minimal impact when the size of the image is large with respect to the support of the convolution filter, which is true in the case of most observation systems.

Let us denote by

• y,

• n,

•

x and • h the DFTs of the observations, noise, original textured image and PSF, respectively. Based on the white noise and circulant PSF assumptions, the law for the noise in the Fourier domain is separable, making way for significant mathematical simplifications and, consequently, performance enhancements. From the noise law, Chapter 2. Inverse Problems, Bayesian Framework and Stochastic Sampling the likelihood can be straightforwardly deduced and writes:

f ( • y|γ n , η, • x) = (2π) -P γ P n exp -γ n P p=1 | • y p - • h p (η) • x p | 2 (2.4) where • h p (η)
shows explicitly the parametric form of the Transfer Function (TF), the frequency domain counterpart of the PSF.

In order to exploit this law separability, the mathematical developments are done in the frequency domain. Consequently, we will be dealing with complex quantities:

• y, • n, • x,
• h ∈ C P . Moreover, for model simplicity and computational efficiency, we have also considered that the spatial domain variables are complex: y, n, x, h ∈ C P . In this manner, the supplementary conditions of spectrum symmetry are avoided.

Remark: The expression in Equation (2.4) takes into account the fact that we are dealing with complex quantities and has the form of a complex Gaussian. This likelihood law describes our direct problem and is the basis of the inference to determine the unknown quantities. Nevertheless, it is straightforward that this law is not sufficient for the estimation, since the number of unknowns is higher than the number of observations. For this reason, regularization is necessary in order to determine a unique solution. In the following, we will give more detailed considerations concerning the regularization and the estimation process.

Full Inversion Problem

Let us consider the original image x is composed of one or several regions R. Each of these regions consists in a single stationary texture patch, belonging to one of the K classes of textures. No assumption is made on the relation between R and K, since none of them is probabilized. K is known and R is not represented in the mathematical developments.

In order to fully describe the original image x, a set of hidden variables (the hidden label field) z can be introduced. The variable z p = k indicates for each position p that the pixel x p is extracted from the full texture x k . Consequently, the prior for the image can be written as f (x|z, x 1 , ..., x K ).

Let k = 1...K be the model index. Then, each texture x k can be fully described by its model, indexed by the discrete variables M k , and the corresponding model parameters χ k . Further details regarding the texture models employed in this work are given in Chapter 3, for the moment the prior law for each textured image being written generically as f (x k |M k , χ k ). In fact, x is obtained deterministically from x 1 ...x K based on the variables z. For this reason, x is not represented in the variable hierarchy reresented in Figure 2.1. This global hierarchical model for our complete problem is based on the aforementioned interdependencies and the observation model presented in the previous section.

By analyzing this hierarchical dependency, the elevated complexity of this problem is obvious. Starting from the observations y, the goal is to estimate:

y η γ n x 1 χ 1 M 1 x K χ K M K z Figure 2
.1: Hierarchical model for the complete problem.

-the unknowns of the observation system γ n and η, -the original image x, -the hidden labels z, -the model M k and parameters χ k of the textures x k . This problem can only be solved by imposing a series of constraints, this being a typical procedure in the field of inverse problems. We have chosen to use a Bayesian approach, where the regularization is achieved through the priors on the parameters.

Bayesian Approach -Prior, Joint and Posterior Laws

The Bayesian approach [START_REF] Gelman | Bayesian Data Analysis[END_REF][START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF] consists in writing the a posteriori law for all the unknowns, given the observations, and using this law in order to obtain estimates for the unknowns. The posterior law can be determined from the joint law:

f (γ n , η, z, x 1...K , M 1...K , χ 1...K |y) = f (y, γ n , η, z, x 1...K , M 1...K , χ 1...K ) f (y) ∝ f (y, γ n , η, z, x 1...K , M 1...K , χ 1...K ) (2.5)
This law, based on the variable dependency presented in Figure 2.1, writes:

f (y, γ n , η, z, x 1...K , M 1...K ,χ 1...K ) = f (y|γ n , η, z, x 1...K ) • π(γ n ) • π(η) • π(z) • k f (x k |M k , χ k ) • k π(M k ) • k π(χ k ) (2.6)
In order to fully specify this law, a priori laws for the unknowns must be chosen, these priors having the role of regularizing the problem. These laws must be completely unrelated to the observations and must reflect general knowledge about the unknowns.

Regarding the hidden labels, the only available prior information is that, irrespective of the observations, the pixels exhibit an aggregation tendency. This means that the labels have a tendency of forming regions with the same value. This can be formalized by using a Potts model for the labels (see Appendix C), with a temperature parameter β that drives the strength of this aggregation behavior. The value of the β parameter determines the mean size of the regions of pixels with the same label.

Uninformative Priors. Conjugacy

The previous considerations concerning the labels have illustrated the influence of the priors on the entire estimation process. Consequently, in cases we do not possess prior information regarding the parameters, in order to avoid biasing the estimation, it is preferable to use uninformative priors. This is the case of γ n , η, β, M 1...K and χ 1...K .

Nevertheless, keeping in mind that the posterior law is used to determine the estimates, it is important from a computational point of view for this law to have standard forms. For this reason, when possible, the priors should have a conjugate form with respect to the likelihood or the other models depending on that parameter. For instance, γ n is the precision parameter in the likelihood law. The conjugate form with respect to this dependency is a Gamma law, i.e.,:

π(γ n |α n , β n ) = β αn n Γ(α n ) γ αn-1 n exp -β n γ n = G (α n , β n ) (2.7)
where α n , β n are considered known. This distribution becomes uninformative in the limit case when α n → 0 and β n → 0. In this case, the Gamma distribution becomes an uninformative Jeffreys law π(γ n ) = 1/γ n .

Remark: The problem in this case is that the Jeffreys law is an improper prior (it is not integrable), which may lead to an improper posterior and may give rise to problems in the estimation process. Nevertheless, in some situations the posterior law is proper despite the presence of improper priors [START_REF] Kass | The Selection of Prior Distributions by Formal Rules[END_REF], this being an open topic of research.

The priors for the rest of the unknowns will be defined in the following chapters, once the texture models f (x k |M k , χ k ) and the dependency h(η) are given explicitly, however, they will all have uninformative forms.

Remark:

The uninformative priors might give the wrong impression that our ill-posed problem is not regularized. In fact, this is done through highly structured models for the PSF, image and labels and not through the priors on the parameters. Moreover, without anticipating too much, we can assert that the posterior law has a complicated form and its dependency with respect to some of the parameters does not have a standard form.

Estimation

The information contained by the aforementioned posterior law concerning the parameter values can be exploited using various estimators: Maximum A Posteriori (MAP), marginalized MAP (mMAP), Median A Posteriori (MeAP) or Posterior Mean (PM). Among these estimators, the PM is the optimal one from the Mean Squared Error (MSE) point of view. This is in fact the reason for which we have chosen it among the various options. Similarly, the evidence based classifier is optimal from the mean classification risk viewpoint. The mathematical proof of these properties is given in Appendix B, where we have analyzed in parallel the two estimators.

Addressed Problems

The previous sections have presented a very complex inverse problem of myopic deconvolution and segmentation of an indirectly observed image consisting in several textured regions. Our goal is to devise a method that can achieve this and at the same time estimate the noise and signals levels, the model and the parameters of each of the textures present in the image. This is a very difficult task and thus we have divided it into three sub-problems: i. myopic deconvolution and parameter estimation of a textured image;

ii. model selection and parameter estimation for a blurred and noisy textured image;

iii. deconvolution segmentation and parameter estimation for a blurred and noisy image composed of different textures.

All these methods, which will be detailed in the following chapters, rely on formulating the posterior law and extracting the information encoded by this law. In every case the form of the posterior will be too complicated, thus numerical methods will be employed. Among the available options, we have chosen to use a Monte Carlo Markov Chain (MCMC) method, more specifically, Gibbs sampling. This iterative method sequentially samples every variable, conditionally on the values of the rest of the variables. In this context, in cases where the conditional a posteriori law for a certain variable has a nonstandard form, more advanced samplers must be embedded within the Gibbs algorithm.

The following section is devoted to the analysis of several numerical methods and the presentation of our algorithmic contribution.

Efficient Metropolis-Hastings Samplers

One of the most commonly employed samplers for complicated laws is the Metropolis-Hastings (MH) algorithm. Then, the solution is to integrate an MH step for θ in the Gibbs sampler. This Metropolis within Gibbs strategy is convergent, as proven by [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF], thus provides samples from the posterior distribution, in our case π(Ψ|y). Since the performances of the sampler directly affect the estimation process, the following section is devoted to the study of efficient versions of the MH algorithm.

The recent literature on efficient sampling is abundant and includes methods that rely either on optimal tuning of the standard MH samplers [START_REF] Haario | An adaptive Metropolis algorithm[END_REF][START_REF]Optimal proposal distributions and adaptive MCMC[END_REF], or on formulating the proposal law based on the target in order to achieve an efficient exploration.

An MH algorithm, as described by [MRR + ], [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], relies on a transition kernel Chapter 2. Inverse Problems, Bayesian Framework and Stochastic Sampling consisting of two ingredients: a transition law, q(θ c , θ p ), and an acceptance probability:

α(θ c , θ p ) = min 1, π(θ p ) π(θ c ) • q(θ p , θ c ) q(θ c , θ p ) (2.8)
where θ c is the current and θ p is the proposed value for the parameters and let Ξ = q(θp,θc) q(θc,θp)

be the ratio of the jumps. The functioning of an MH algorithm is the following:

1. Initialize the iteration counter j = 1 and set θ (0) c . 2. Propose a new value θ p for the parameter, generated from the density q(θ (j-1) c , •).

3. Evaluate the acceptance probability α(θ (j-1) c , θ p ) given by Equation (2.8) and according to this value accept the proposal and update the parameter θ (j) c = θ p , or reject it and keep the old value for the parameter θ (j) c = θ (j-1) c .

4. Update the counter and return to step 2 until convergence.

There are numerous options for formulating the proposal and the convergence speed and mixing properties are directly influenced by the adequacy between the proposal law and the target. The simplest version of this algorithm is the Independent MH (IMH), whose proposal does not depend on the current value θ c . Nevertheless, this algorithm's simplicity is reflected in its inability to use any of the previously acquired information to render the sampling more efficient. For this reason, its use is limited to cases where the target has very complicated forms with mass distributed in a large area of the parameter space. Except for the "independent" version of this algorithm, the other versions can be considered as being Random Walks (RWMH), since their proposals contain the current value of the chain θ c and a stochastic term. The simplest version is the isotropic RWMH, which proposes an isotropic displacement around the current value of the chain.

Moreover, choosing a transition kernel q that embeds information about the shape of the target can significantly enhance the algorithm performances, especially when the target is very peaked, as in our problem. Advanced methods include a directional component for the proposal, meaning that the proposal is built by making an isotropic move around the current value plus the directional component. In the class of directional MH methods, a first idea is to build the proposal based on first order derivatives of the target, this being the case of the Metropolis Adjusted Langevin Algorithm (MALA) [RS03, KKS + 10, GC11, MR12] and of the Hamiltonian methods [START_REF] Girolami | Riemannian manifold Hamiltonian Monte Carlo (with discussion)[END_REF][START_REF] Zhang | Quasi-Newton methods for Markov chain Monte-Carlo[END_REF][START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF][START_REF] Beskos | Hybrid Monte Carlo on Hilbert spaces[END_REF]. Another class of methods that improve the sampling performances is the quadratic (Newton-like) approach, based on second order derivatives (Hessian matrices for multidimensional parameter spaces) and exploiting the information regarding the target curvature. Such methods have first been employed in optimization theory, but they have recently been successfully adapted to sampling [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF][START_REF] Bui-Thanh | A Scaled Stochastic Newton Algorithm for Markov Chain Monte Carlo Simulations[END_REF][START_REF] Martin | A Stochastic Newton MCMC Method for Large-Scale Statistical Inverse Problems with Application to Seismic Inversion[END_REF].

In this work, we have focused our attention on the second order derivative-based proposals, which we have improved by building a new transition kernel based on the Fisher information matrix. Our contribution is modifying the Hessian algorithm presented in [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF] to eliminate the concerns regarding the positiveness of the Hessian matrix. This is achieved by its replacement with the Fisher information matrix. The results obtained for the proposed problem are encouraging, confirming the stable behavior and speed perfor-mances of the second order derivative-based samplers.

All of the analyzed samplers are convergent and accurately explore the target, the differences between them being given by the different formulation of the transition kernel and, implicitly, by the expression of the acceptance probability. These differences translate into:

-different exploration of the target, this influencing how fast the sampler reaches the high probability regions, -differences between the time needed by each sampler to produce one sample. In the following, we will present the specific aspects of these algorithms, adapted to our inverse problems. The target law for these algorithms is the a posteriori law for some of the texture parameters θ. In our case, this law has a very peaked exponential form and its values may overpass Matlab's numerical precision. For this reason, the evaluation of the acceptance probability can lead to numerical problems and indeterminations.

By applying the logarithm to Equation (2.8) we avoid handling exponential quantities. Consequently, in order to evaluate this probability, we only have to compute a difference between the Log-Posterior (LP) taken in the proposed value and in the current value, where LP = log π(θ|y). Moreover, instead of computing the second ratio in Equation (2.8), it suffices to compute the value of Ξ:

log α = min {0, LP(θ p ) -LP(θ c ) + log Ξ} (2.9)
Furthermore, let us consider that the prior law for θ is uniform π(θ) = U [θm,θ M ] (θ). We can then write:

LP(θ) = LL(θ) + log π(θ) -log f (y) (2.10)
where LL denotes the Log-Likelihood, and log π(θ) = -log (θ M -θ m ) and log f (y) are constant with respect to θ. Then, (2.9) can be rewritten as:

log α = min {0, LL(θ p ) -LL(θ c ) + log Ξ} (2.11)
In the following, we will present the various versions of the MH algorithm and their specificities in terms of computational complexity and efficient parameter space exploration.

Independent

The proposal for this algorithm writes:

θ p = U [θm,θ M ] (θ) (2.12)
where the values of θ m , θ M are the limits of the parameter space. In this case, the transition kernel is independent of the current value θ c and writes:

q(θ c , θ p ) = U [θm,θ M ] (θ p ) (2.13)
Chapter 2. Inverse Problems, Bayesian Framework and Stochastic Sampling

The acceptance probability becomes:

log α = min 0, LL(θ p ) -LL(θ c ) + log U [θm,θ M ] (θ c ) -log U [θm,θ M ] (θ p ) = min {0, LL(θ p ) -LL(θ c )} (2.14)
since the transition kernel is uniform and has the same value for any θ.

As shown by the previous equations, this algorithm consists in very simple and inexpensive iterations. Nevertheless, in most cases this proposal is not adapted to the form of the target, leading to high rejection rates and, consequently, to very slow convergence.

Standard Random-Walk

Taking into account the information that the current value θ c has been considered adequate enough to be accepted can improve the performances of the MH sampler. This information is embedded in the proposal, which becomes:

θ p = θ c + ε N (0, I) (2.15)
In this case, the value of ε sets the size of the jumps. This is very important since a small value of the step results in a slow exploration of the parameter space and strongly correlated samples, while a large value may lead to high rejection rates, as in the case of the IMH. The optimal value for ε, in the sense that it provides the best compromise between the amplitude of the transitions and the acceptance rate is the one that yields an acceptance rate of approximately 24% [START_REF] Gelman | Bayesian Statistics[END_REF].

In this case, the transition kernel writes:

q(θ c , θ p ) = exp - 1 2ε 2 θ p -θ c 2 (2.16)
The acceptance probability becomes:

log α = min 0, LL(θ p ) -LL(θ c ) - 1 2ε 2 θ c -θ p 2 + 1 2ε 2 θ p -θ c 2 = min {0, LL(θ p ) -LL(θ c )}
(2.17) due to the symmetry of the transition kernel.

The isotropic RWMH algorithm also has a reduced cost per iteration. Nevertheless, depending on the initialization and the value of the tuning parameter, it may be slow to converge. Moreover, it does not exploit the available information concerning the form of the target in order to increase sampling efficiency.

Langevin adapted Random-Walk

Langevin algorithms are derived from diffusion approximations and rely on the principle of using the information concerning the target density, in the format ∇ log π, in order to build a proposal distribution well-adapted to the problem in question [START_REF] Gelman | Bayesian Statistics[END_REF]. The Langevin-based MH algorithm proposes an RW-like transition of the form [START_REF] Besag | Bayesian computation and stochastic systems (with discussion)[END_REF]:

θ p = θ c - ε 2 2 g(θ c ) + ε N (0, I) (2.18)
where g(θ c ) = ∂LL(θ) ∂θ | θ=θc . The acceptance probability can be obtained from Equation (2.11), by using the transition kernel:

q(θ p , θ c ) = exp - 1 2ε 2 θ c -θ p - ε 2 2 g(θ p ) 2 (2.19)
Then, this acceptance probability writes:

log α = min 0, LL(θ p ) -LL(θ c ) - 1 2ε 2 θ c -θ p - ε 2 2 g p 2 -θ p -θ c - ε 2 2 g c 2
and the log ratio of the jumps becomes:

log Ξ = - 1 2 (θ p -θ c ) t (g p -g c ) + ε 2 2 (g t p g p -g t c g c ) (2.20)
As compared to the non-directional MH methods, the complexity of the Langevin MH is increased due to the form of the acceptance probability and to the necessity of evaluating the gradient for every new proposal. However, the increased computation time per iteration is compensated by the smaller number of iterations needed to reach convergence.

In regions far from the maximum of probability, the gradient is large (the directional component is dominant), i.e., the algorithm approaches the high probability regions with high amplitude jumps. Near the maximum of probability, the gradient is small, thus the stochastic component is dominant and it permits for the region to be explored.

Hessian adapted Random-Walk

This section is devoted to a sampling method seldom explored and whose presence in the literature is scarce. The directional component of the proposal is in this case formulated using Newton's direction, which, for a quadratic law, indicates the maximum. In [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF] a version of this sampler has been tested and compared to methods such as Gibbs and optimal marginal data augmentation (DA) samplers on a probit regression problem, proving that the performances of this sampler are superior. The transition for the Hessian sampler is of the form:

θ p = θ c + ε Σ(θ c ) g(θ c ) + N (0, Σ(θ c )) (2.21)
where Σ(θ c ) = -H(θ c ) -1 and H(θ c ) is the Hessian matrix of the target law, computed in the value θ c . The acceptance probability is obtained from (2.11), for:

q(θ p , θ c ) = N (θ c -θ p + ε Σ(θ p ) g(θ p ), Σ(θ p )) (2.22)
Chapter 2. Inverse Problems, Bayesian Framework and Stochastic Sampling

In this case, the acceptance probability has the form:

log α = min 0, LL(θ p ) -LL(θ c ) - 1 2 [log (det (Σ(θ p ))) -log (det (Σ(θ c )))] - 1 2 θ c -θ p + ε Σ(θ p ) g(θ p ) 2 Σ(θp) + 1 2 θ p -θ c + ε Σ(θ c ) g(θ c ) 2 Σ(θc)
and the log ratio of the jumps in this case writes:

log Ξ = - 1 2 [log (det (Σ(θ p ))) -log (det (Σ(θ c )))] - 1 2 θ c -θ p 2 Σp-Σc -2ε(θ c -θ p ) t (g p + g c ) + ε 2 g t c Σ(θ c )g c -g t p Σ(θ p )g p (2.23)
The advantage of the method is that, for quadratic and quasi-quadratic distributions, the regions of high probability are approached in a very small number of iterations (ideally, a single one) and then explored with the contribution of the stochastic component, a quadratic law of variance Σ(θ c ), which is an accurate approximation of the target. However, it is clear that this method is also rather complex, as each iteration translates in the computation of the gradient and the Hessian matrix and the evaluation of the acceptance probability.

In the context of quadratic directional methods, a major concern is the need to perform the inversion of the Hessian, as this may be problematic if the Hessian is not positive definite.

Fisher adapted Random-Walk

In order to avoid the aforementioned matrix inversion problems, in the present work we have explored the idea of replacing the Hessian by the Fisher information matrix (see Appendix A), which is by definition positive definite:

I pq (θ) = E x|θ - ∂ 2 ∂θ p ∂θ q LL(θ) θ (2.24)
The Fisher matrix quantifies the mean amount of information that the observations contain regarding the parameter θ. We have been able to apply such an approximation, as in the case of our problem, we have a great amount of independent observations, thus, a scenario close to the asymptotic case. This efficient sampler has been presented in [START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF] under the name of Fisher adapted MH (FMH).

A very interesting sampler, developed independently from our own is the manifold MALA (mMALA) [START_REF] Girolami | Riemannian manifold Hamiltonian Monte Carlo (with discussion)[END_REF], dedicated to highly-dimensional and highly-correlated targets. Its development is based on the symmetric Kullback Leibler divergence, D S (p||q), a first order Taylor approximation, f (y|θ+δθ) ≈ f (y|θ)+δθ t ∇ θ f (y|θ) and the approximation log(1 + ε) ≈ ε:

D S [f (y|θ + δθ)||f (y|θ)] = = δθ t • E y|θ ∇ θ log f (y|θ) • ∇ θ log f (y|θ) t I (θ) •δθ (2.25)
The mMALA is obtained by defining the Langevin diffusion with invariant measure π(θ|y). This algorithm formulates the proposal:

θ pr = θ c + ε 2 2 I -1 (θ c ) • ∇ θ LP(θ c ) + ε I -1 (θ c ) • z c (2.26)
with LP(θ) = log π(θ|y) the log-posterior and z c ∼ N (0, I) an isotopic displacement.

Our FMH algorithm, that uses the same proposal law as the one in Equation (2.26), is in fact based on the idea of quasi-Newton proposals [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF] for a fast exploration of the parameter space and superior mixing properties and replacing the Hessian with the Fisher information matrix. The use of this proposal proved advantageous from multiple points of view. Firstly, this exploits the target curvature similarly to the Newton step from the optimization theory. Secondly, this made way for a series of algorithmic simplifications and performance enhancement:

• for our GRF modeled textures (described in Chapter 3), the second order derivatives vanish under the expectation. This allows the formulation of the efficient proposal only based on first order derivatives,

• the Fisher matrix is positive definite. Hence, if its eigenvalues are not zero, there are no instabilities when taking the inverse, such as those mentioned in [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF] for the use of the Hessian,

• also due to the positive definite Fisher matrix, the Newton term always has the direction of gradient ascent, thus the algorithm only makes efficient steps. previously presented. As expected, the quadratic methods require the least number of samples to converge. Moreover, due to the simplification achieved using the Fisher matrix, the overall convergence time is the most reduced for the FMH.

Conclusion

This chapter has presented the formulation of our full inverse problem, theoretical aspects regarding the Bayesian framework that is employed, its mathematical formulation and practical issues related to the implementation. We have specified our choice of estimator and the reason behind this choice. In this context, sampling is used to numerically compute the estimator and several versions of the MH sampler are presented.

Our contribution consists in developing the new efficient sampler FMH, based on the Fisher information matrix. This sampler will be employed in one of the inverse problems presented in the following chapters in order to enhance the speed of our sampling based method. As a perspective on this topic, this sampler can be integrated in numerous other sampling based applications in order to accelerate the sampling speed. Texture represents one of the main elements of the framework explored in this work. Inverse problems such as myopic deconvolution, parameter estimation, model choice, joint deconvolution + segmentation, of high interest in signal and image processing, are addressed in the special case of textures. Our work is set apart from the existing literature by the fact that we are not only dealing with texture, a rather challenging aspect in itself, but we are also considering indirect observations, thus introducing an additional complexity level.

Generally, ill-posed inverse problems are regularized through priors on the parameters. In this work, however, the highly structured image models (texture models) provide the needed constraints for the regularization.

Introduction and State of the art

Referring to texture in general, it is an omnipresent aspect in everyday life, the Oxford Dictionaries defining it as "the feel, appearance, or consistency of a surface or a substance". Consequently, texture is one of the main characteristics of any material object. Moreover, the notion of texture has also been extended to more abstract quantities and has become a central aspect in signal processing. For instance, we can find a texture in an audio signal or in a data sequence. Nevertheless, most of the texture-related applications can be found in image processing, [START_REF] Stockman | Computer Vision[END_REF] defining the image texture as the "set of metrics calculated in image processing designed to quantify the perceived texture of an image". This is only one of the numerous definitions and descriptions that have been formulated in the literature, each influenced by the application, the author's background or the approach being employed.

Texture modeling represents an important aspect in image processing applications, the goal being to find comprehensive models that allow a good description, but are relatively easy to handle and integrate in a more complex image processing formalism. Consequently, there is always a compromise to make between model versatility and complexity.

In the context of a growing interest towards texture in image processing applications, in the late 1990s several extensive works have tried to make an inventory of the texture types, approaches and applications. Amongst the most well-known, [START_REF] Tuceryan | Texture Analysis[END_REF] provides a valuable overview of texture definitions, incentives for taking interest in texture, fields of applications, types of modeling and image processing problems that can arise in the context of textured images. The authors classify the approaches to texture analysis into four main classes:

1. statistical methods -co-occurrence matrices, autocorrelation features, 2. geometrical / structural methods -Voronoi tessellation features, 3. model based methods -Random Field (RF) models, fractals, 4. signal processing / transform-based methods -spatial domain filters, Fourier domain filters, Gabor and wavelet methods.

This classification, also used in [START_REF] Materka | Texture Analysis Methods -A Review[END_REF], has been enriched over time by the development of new modeling techniques. However, the main classes remain the same, as proven by its use in the more recent [START_REF] Bharati | Image texture analysis: methods and comparisons[END_REF].

A great amount of attention has been dedicated to feature extraction. For instance, [START_REF] Sharma | Evaluation of texture methods for image analysis[END_REF] provides an evaluation of five feature extraction methods: autocorrelation, edge frequency, primitive-length, Laws' method [START_REF] Laws | Rapid Texture Identification[END_REF], and co-occurrence matrices. [START_REF] Srinivasan | Statistical Texture Analysis[END_REF] offers a more comprehensive overview of statistical texture analysis. The authors include in their comparison transform-based methods such as Fourier, Gabor functions and wavelet transforms. Although in contradiction with the previous classification, which separates the statistical and the transform approaches, this point of view is coherent and generalizes the statistical approach beyond the spatial representation domain.

Texture modeling in a transformed space, as for instance the wavelet domain, has been explored for decades, [START_REF] Chang | Texture analysis and classification with tree-structured wavelet transform[END_REF][START_REF] Arivazhagan | Texture segmentation using wavelet transform[END_REF] pointing out the main advantage of the treestructured wavelet transform that provides non-redundant representations, with wavelet coefficients that take at most the storage space of the original image. Moreover, as opposed to the Laws or Gabor filters, the filter coefficients do not require tuning and do not differ from one texture to another. The wavelet based approach has been used mostly in texture classification and retrieval applications, as is the case of [START_REF]Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF], where a signature is defined for each texture and the classification is based on the Kullback-Leibler distance, or of [START_REF] Do | Rotation invariant texture characterization and retrieval using steerable wavelet-domain hidden Markov models[END_REF], where a rotation invariant method is devised based on wavelet domain Hidden Markov Models (HMM). [FX03] uses a similar wavelet domain HMM principle and exploits the cross correlation across the sub-bands. A texture synthesis algorithm based on steerable pyramids with complex coefficients is given in [START_REF] Portilla | A Parametric Texture Model based on Joint Statistics of Complex Wavelet Coefficients[END_REF].

In the same class of transform based methods, frequency analysis represents a very important approach. [START_REF] Brown | Texture synthesis and prediction error filtering[END_REF] presents the idea behind Fourier transform based texture synthesis, i.e., the fact that all stationary texture realizations should have the same amplitude spectrum. Since the phase spectrum localizes the frequency components in space, in the case of stationary textures, this can differ between the realizations. In the same line of work, the recent [START_REF] Galerne | Random Phase Textures: Theory and Synthesis[END_REF] shows two very interesting aspects: firstly, that "random phase textures and random shift textures generated from the same sample are indistinguishable" and, secondly, that "random phase textures are perceptually invariant under a multiplicative noise on the Fourier modulus". Moreover, texture features based on the Fourier transform can be very efficient in classification tasks, as shown in [START_REF] Zhou | Texture feature based on local Fourier transform[END_REF], where a feature based on histograms of local Fourier coefficient maps is tested for image retrieval and in [START_REF] Ahmad | Texture features based on Fourier transform and Gabor filters: an empirical comparison[END_REF], in a comparison with Gabor based features. In the dynamic textures field, [START_REF] Abraham | Dynamic texture with Fourier descriptors[END_REF] reports a significant speed enhancement of the learning process and the need for a smaller learning set due to the use of a set of Fourier descriptors on the frames, instead of the raw sequence, for the linear dynamic system identification. [START_REF] Li | A new texture generation method based on pseudo-DCT coefficients[END_REF] is also worth mentioning since it provides an interesting method for 2D texture synthesis starting from a 1D signal, by using the Fourier coefficients of the 1D signal and arranging them to form a so called pseudo Discrete Cosine Transform (DCT) and then applying and inverse DCT to obtain the texture.

Among the stochastic representations, the Wold decomposition model has also been used for texture modeling [START_REF] Sriram | Texture coding using a Wold decomposition model[END_REF]. It consists in decomposing the texture field into a sum of two orthogonal components: a deterministic and an indeterministic component. The deterministic component accounts for the structural properties, and is further decomposed into a harmonic and a finite number of evanescent components, while the structureless indeterministic component accounts for the randomness of the texture. [START_REF] Liu | Periodicity, directionality, and randomness: Wold features for image modeling and retrieval[END_REF] defines the main dimensions of human texture perception: "periodicity", "directionality" and "randomness", and provides a more robust method for texture modeling and classification. Moreover, [START_REF] Campisi | Nonlinear Prediction in the 2D Wold Decomposition for Texture Modeling[END_REF] shows that the indeterministic component can be predicted from the deterministic one using suitable nonlinear schemes. The more recent [START_REF] Stitou | 3-D Texture Characterization Based on Wold Decomposition and Higher Order Statistics[END_REF] presents a parametric model for representing 3D textures, which describes both spatially and spectrally each component of the Wold decomposition.

The statistical approaches are not only adapted for texture analysis, by inferring on the underlying model parameters, but also for texture synthesis. This can be achieved by characterizing a set of visually similar textures through a probability distribution on a RF. It is a two step process, firstly, the feature extraction phase uses a filter set to capture the texture features characteristic to that class. The histograms of the filtered images are used to estimate the marginal distributions on the RF. Secondly, the feature fusion phase consists in determining the maximum entropy distribution that matches the previously determined marginals. This approach, based on the minimax entropy principle has been presented in [START_REF] Zhu | Minimax Entropy Principle and Its Application to Texture Modeling[END_REF] and then revisited in [START_REF] Zhu | Filters, Random fields And Maximum Entropy (FRAME) -Towards a Unified Theory for Texture Modeling[END_REF]. The maximum entropy based texture analysis has been employed in applications such as estimating the crystallite orientation distribution function [START_REF] Bohlke | Application of the maximum entropy method in texture analysis[END_REF], texture and object recognition [START_REF] Lazebnik | A Maximum Entropy Framework for Part-Based Texture and Object Recognition[END_REF].

In a more structural oriented approach, [START_REF] Malik | Contour and Texture Analysis for Image Segmentation[END_REF] tackles the image segmentation problem using a method based jointly on a texton based texture analysis and on contours. Another interesting application is the texture synthesis method in [START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF] inspired from the Markov RF (MRF) texture models and relying on a deterministic searching process.

The special properties of textures somewhat resemble those of fractals, which are defined by self-similarity and scale invariance. The theory of fractals has thus inspired the use of new descriptors for texture analysis, namely the fractal dimension and lacunarity [CKC93]. The problem of texture segmentation based on the fractal dimension has been presented in [START_REF] Chaudhuri | Texture Segmentation Using Fractal Dimension[END_REF], while the more recent [START_REF] Lee | A robust algorithm for the fractal dimension of images and its applications to the classification of natural images and ultrasonic liver images[END_REF] provides a robust and efficient method for computing the fractal dimension and its application in image classification. However, the fractal dimension characterizes self-similarity only in ideal cases, this being the reason for using multifractal analysis, this idea being explored in [START_REF] Xia | Morphology-based multifractal estimation for texture segmentation[END_REF] for texture segmentation.

One of the most accessible methods, equally adapted to synthesis and analysis tasks, is RF modeling, with good performances for both stochastic and deterministic-like textures. This model has been explored for decades, [START_REF] Kashyap | Univariate and multivariate random field models for images[END_REF] introducing the idea of using RFs for image modeling and [START_REF] Sklansky | The use of Markov Random Fields as models of texture[END_REF] the use of MRFs for texture modeling. These ideas were further explored in [START_REF] Cross | Markov Random Field Texture Models[END_REF] that investigates the representation capabilities of MRFs for the synthesis of microtextures, regular textures and inhomogeneous textures. Furthermore, [START_REF] Chellappa | Texture synthesis using 2-D noncausal autoregressive models[END_REF] shows that non-causal Auto Regressive (AR) models often have the oscillatory behavior characteristic to textures and that they can be successfully used for texture synthesis, and [START_REF] Chellappa | Classification of textures using Gaussian Markov random fields[END_REF] provides two feature extraction methods based on the Gaussian MRF (GMRF) model for texture classification. [START_REF] Levina | Texture synthesis and nonparametric resampling of Random Fields[END_REF] presents a texture synthesis algorithm based on Markov mesh models, a special case of MRF. Recent works such as [START_REF] Cariou | Unsupervised texture segmentation/classification using 2-D autoregressive modeling and the stochastic expectation-maximization algorithm[END_REF], based on the AR model for image segmentation and texture classification, also illustrate the interest of this class of methods, which are rather simple, with a reduced number of parameters and easy to handle.

Valuable insight into the use of Gibbs RFs for texture modeling is provided by [START_REF] Gimel'farb | Image Textures and Gibbs Random Fields[END_REF], along with the theoretical connections between the MRFs and the Gibbs RFs, in a local versus global expression of image energy. The Gibbs model is defined based on the clique potentials, and these can be tailored in such a manner as to model multiple translation invariant pairwise pixel interactions, by making each pixel appear simultaneously in several cliques. Moreover, this extensive analysis covers the difference between Gibbs-Markov fields and Gibbs-non-Markov fields. From an algorithmic point of view, practical aspects related to texture synthesis or learning are reviewed and put into context, such as simulated annealing, the Maximum Likelihood Estimator (MLE), which are then used in applications of texture synthesis, classification, both in the case of synthetic and natural textures.

In the same spirit of RFs, [START_REF] Bouman | A Multiscale Random Field Model for Bayesian Image Segmentation[END_REF] solves the problem of image segmentation using a multiscale RF model for the image. The very recent [START_REF] Bouman | A Guide to the Tools of Model Based Image Processing[END_REF] provides an overview of model based image processing, centered on stochastic models, either causal Gaussian, such as the AR model, or non-causal, such as the GMRF, and their comparison. Furthermore, it explores the non-Gaussian MRFs and their connection to the Gibbs RF, and also the various types of clique potentials ρ(∆), either convex:

• Gaussian : |∆| 2 /2,
• Total Variation : |∆|,

• Generalized GMRF : |∆| p /p,

• Huber :

∆ 2 /2 for|∆| < T T |∆| -T 2 /2 for|∆| > T, or non-convex • Blake-Zisserman : min {|∆| 2 , T 2 }, • Geman-McClure : ∆ 2 ∆ 2 + T 2 , • Hebert-Leahy : log(∆ 2 + T 2 ),
• Geman-Reynolds :

|∆| |∆| + T ,
Furthermore, the Ising model is also presented, as a particular case of a discrete MRF. Moreover, this work provides valuable information regarding the algorithmic implementation of stochastic sampling, i.e., implementation options for the described models.

A series of works devoted to finding accurate representations for "natural/photographic" images have also been carried out in the wavelet domain and tried to determine adequate decompositions and models for the wavelet coefficients. For instance, [START_REF] Wainwright | Random Cascades on Wavelet Trees and Their Use in Analyzing and Modeling Natural Images[END_REF] defines random cascades on trees of multiresolution coefficients, these cascades reproducing a Gaussian Scale Mixtures (GSM), a special class of random variables. This model is also presented in [START_REF] Eero | Statistical Modeling of Photographic Images[END_REF] in an overview of statistical modeling techniques. The GSM wavelet coefficients image model is further used in an image denoising application [PSWS03] and refined in [START_REF] Hammond | Image Modeling and Denoising With Orientation-Adapted Gaussian Scale Mixtures[END_REF], to deal with orientations, and in [START_REF] Lyu | Modeling Multiscale Subbands of Photographic Images with Fields of Gaussian Scale Mixtures[END_REF], for multiscale subband modeling. Although not specifically devoted to texture representation, this model successfully captures the characteristics of natural images, which in most cases contain texture. The GSM modeling in the wavelet domain has represented one of the inspirations for our non-Gaussian RF based texture model that will be presented in the following section. *** Far from being exhaustive, this section has presented a brief review of the main branches of approaches to texture representation and reference papers in the texture related literature. More or less present in the literature, each of these models has its advantages and its own set of texture classes and applications for which it outperforms from a certain viewpoint the concurrent approaches. Consequently, the question of which is the most appropriate texture model is not trivial and does not have an unique answer.

The intention of embedding the texture model into a complex image processing application has definitely biased our choice. Firstly, we have focused our attention towards models that have both analysis and synthesis capabilities and, secondly, towards highly structured models that manage to capture the texture specificities in a relatively reduced set of parameters.

Texture Modeling by Random Fields

The statistical approach towards the addressed inverse problems has only made it coherent that a statistical model should also be employed for the texture. Moreover, the previous section has presented a part of the work carried in this domain, which confirms the extended representation capabilities of this type of models, both in texture analysis and in texture synthesis tasks.

Our attention has focused on indirect observations, i.e., textures corrupted by a blur and by noise. This convolution based observation system, under the circulant assumptions detailed in Chapter 2, can be expressed as a simple multiplication and addition in the Fourier domain. This has inspired us into defining a texture model directly in the Fourier domain. This model is able to perform both analysis and synthesis tasks, as opposed to the majority of the works based on this representation domain, which are mostly focused on extracting relevant features for texture analysis.

The main idea is to consider that the Fourier coefficients are modeled by an RF and to assign them a probability distribution.

Remark: A very important assumption is that the textures are stationary.

The direct implication of this assumption is that the covariance matrix has a Toeplitz-block-Toeplitz structure. Moreover, by Whittle approximation [P.W], this matrix has a Circulantblock-Circulant form, thus is diagonalizable by Fourier transform. Consequently, the Fourier coefficients,

• x p , p = 1...P , are decorrelated.

Remark: Furthermore, in the case where the RF is Gaussian, decorrelation and independence are equivalent, thus the Fourier coefficients are actually independent.

From an algorithmic point of view, the independence of the Fourier coefficients implies computational efficiency, since all the computations can be done in parallel for all the coefficients. Nevertheless, from a modeling viewpoint, this represents a limitation.

The Gaussian RF (GRF) model is fully described by the second order statistics of the coefficients. The literature on texture modeling records a series of works carried in the field of human perception, which tried to establish the perception sensitivity to the higher order statistics. The first conjecture has been formulated in the early work [Jul62] and stated that humans were unable to preattentively distinguish between textures having the same second order statistics and different higher order statistics, called isodipoles. However, this conjecture has been proven wrong, by Julesz himself and a series of other researchers. Preattentively distinguishable isodipoles and even isotrinomes (textures with identical third order statistics and different higher order statistics) have been synthesized [JGSF73, Jul80, Vic94, MN01], see Figure 3.1 for an example. Nevertheless, all these textures share a common characteristic, the fact that they are very structured, deterministiclike.

The conclusion that can be drawn from this extensive work is that, although not capable of describing all texture classes, the second order statistics encode however a great amount of information regarding textural content. This shows that the GRF model, despite its limitations, is an interesting tool in texture processing.

In this work we are only dealing with stationary textures, consequently, we will not be interested in the phase field [START_REF] Brown | Texture synthesis and prediction error filtering[END_REF]. Thus, the useful information will be encoded in the modulus, and the phase field may differ among realizations of textures with the same modulus field. In order to provide a simple modeling we chose to impose constraints directly on the PSD, without separately controlling the real and the imaginary components.

In the following we will present two texture models based on GRFs: 1. a model based on a Gaussian law for the texture Fourier coefficients, 2. a model based on a non-Gaussian law for the texture Fourier coefficients, with enhanced representation capabilities due to the thicker tails of the law. The models are fully specified by the mean and the PSD, or by the mean, PSD and auxiliary variables law, respectively. In this work we will be treating the zero-mean case. Both of these GRF based models consist in formulating a Gaussian law for the Fourier coefficients of the image, conditionally on the PSD.

Then, the conditional law for the image in the spatial domain can be written as:

f (x|R x ) = (2π) -P |R x | -1 • exp -x 2 Rx (3.1)
where R x is the Circulant-block-Circulant covariance matrix. Due to this particular structure, the Fourier coefficients are independent and thus the previous expression has a separable counterpart in the frequency domain:

f ( • x|χ) = (2π) -P P p=1 χ p • exp - P p=1 χ p | • x p | 2 (3.2)
where χ p represents the element of the inverse PSD at position p. More specifically, χ p = s p λ p is the product of two components: -s p is the scale parameter, -λ p is the shape component. In other words,

• x p |χ p ∼ N 0, χ -1 p (3.3)
i.e., the Fourier coefficients are independent and non-identically distributed (inid).

Gaussian Model

The GRF texture model is very simple, while at the same time capable of representing rather complex textural contents. This is in fact the exact model employed in our work on fast texture parameters sampling [START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF], and on pixel interaction model choice [START_REF] Vacar | Bayesian Texture Model Selection using Harmonic Mean[END_REF]. In the following, we will present this model in detail, along with typical texture realizations.

Chapter 3. Texture Modeling

The particularity of this model is that all the scale parameters are equal: s p = s, for p = 1...P . Then, (3.2) becomes:

f ( • x|s, λ) = (2π) -P s P P p=1 λ p • exp -s P p=1 λ p | • x p | 2 (3.4)
i.e., the conditional law of the Fourier coefficients differs only in the shape component of the PSD, λ p :

• x p |s, λ p ∼ N 0, (sλ p ) -1 . (3.5)
In order to fully specify the model, the prior for s must also be set. For conjugacy reasons, the chosen law has a Gamma form. Then, the GRF texture model can be hierarchically represented as:

x λ sb Figure 3.2: GRF texture model
A texture realization can be obtained using virtually any form for the PSD, from the simplest uniform texture, obtained from a PSD with a single impulsional component at the null-frequency, to a random noise, corresponding to a uniform PSD. The complexity of the texture is directly related to the richness and the structure of its PSD. Figure 3.3 shows a series of texture realizations with some of the simplest possible PSDs. We have chosen combinations of one, two, four and nine pure frequency components, respectively. The goal is to illustrate the type of patterns that can be obtained with such a reduced number of components. Obviously, the number and position of these frequency components determine the characteristics of the corresponding texture. Consequently, a richer spectral content would provide more complex patterns. Moreover, the goal is to find an efficient manner to encode this richer spectral content and one option is the use of parametric models.

The idea of using parametric models for the PSD is very practical, since textures with more complex spectral contents can be obtained using a reduced number of parameters. Figure 3.4 shows a series of parametric PSD shapes, along with the corresponding texture realizations. Notice the effect of changing the position of the component, on the textures in Figures 3.4a and 3.4b, and the effect of changing the parametric shape, between the textures in Figures 3.4a and 3.4c.

In these cases, the models depend on at most four or five parameters, which, along with the model for the PSD shape, fully encode the information regarding the texture. As it can be seen in the presented realizations, these textures are more complex and resemble natural textures. The textures obtained using such a continuous PSD are very stochastic-like, with patterns where it is hard to identify a structural element.

Nevertheless, more structured, deterministic-like textures can also be obtained using this model. This can be achieved starting from the previously mentioned continuous forms, by selecting only a part of the frequency components. Such textures are shown in Figure 3.5, being obtained using a uniform sampling lattice over the entire frequency domain.

To illustrate the effect of this sampling, the texture in Figure 3.5c is obtained from the same PSD as the texture in Figure 3.4e, the difference between the two realizations being significant. Furthermore, the continuous PSD shapes used to obtain the last two textures consist in the sum of two Lorentzian components and two Laplacians, respectively.

These uniformly sampled PSDs are also parametric and, as compared to their continuous counterparts, only contain two supplementary parameters indicating the sampling Figure 3.3: Texture realizations for the zero-mean GRF model. On the first and third rows are illustrated various forms for the PSD and underneath each PSD, the corresponding texture. The spectral content of these textures is very reduced, consisting in a small number of punctual frequency components (impulse like PSD).

Chapter 3. Texture Modeling frequencies on the two axes. Moreover, more complex sampling lattices can be used, lattices that can themselves be parametric.

Remark: These textures have a more obvious geometrical pattern, with a structural element that repeats itself, this giving them a more deterministic-like character.

The textures we have presented cover only a small part of the large variety of texture classes that can be modeled using GRFs and second order statistics. However, they show that rather complex textural content can be encoded in a reduced number of parameters and motivate the interest we have taken in this model. The λ p from Equations (3.2) and (3.4) are the elements of the PSD field and depend on the model, M = k, and on the parameters of that model, θ k . From this point forward, to explicitly show this dependency the notation λ k p (θ k ) is used. Although the use of the aforementioned parametric PSDs means that the entire information regarding the textural content is encoded in a small number of parameters, the problem does not become trivial. The texture parameters will not be easy to estimate, due to the non-linearity of the data with respect to the parameters. To be more exact, λ k p (θ k ) has a highly non-linear, non-standard form. Nevertheless, this model has its limitations. Among these limitations, the Gaussianity of the Fourier coefficients can be overcome almost effortlessly. This will be represented in the following section as our non-Gaussian texture model.

Non-Gaussian Model

The previously presented Gaussian model has the main advantage that the texture Fourier coefficients are independent and follow Gaussian laws. This is extremely advantageous from a practical point of view, since it allows for parallel processing, of standard laws, on the entire texture. The idea behind the new texture model is to obtain non-Gaussian distributed Fourier coefficients, while keeping the aforementioned advantages. This aspect has already been widely explored in statistics, but also in image processing [START_REF] Geman | Nonlinear Image Recovery with Half-Quadratic Regularization[END_REF][START_REF] Geman | Constrained Restoration and the Recovery of Discontinuities[END_REF][START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF] to define the regularization terms in problems of deconvolution and denoising. In such problems, the data adequacy term is Gaussian with respect to the object of interest (the original image) and the goal is to conserve the overall Gaussianity, for computational efficiency. At the same time, the aim is to eliminate the limitations that a Gaussian regularization term imposes, all this being achieved by the introduction of a set of auxiliary variables. This model is called the Gaussian Scale Mixture (GSM).

A similar principle has been presented in [START_REF] Aykroyd | Bayesian Estimation for Homogeneous and Inhomogeneous Gaussian Random Fields[END_REF] with Markov GRFs for image modeling. Furthermore, extensive work [WSW01, PSWS03, Sim05, HS08, LS09] has been carried for finding accurate models for natural images. The approach explored in this case relies on modeling the coefficients of the wavelet decomposition using GSMs. Although the idea of conditional Gaussianity and marginal non-Gaussianity has been previously used, to the best of our knowledge, it has not been exploited so far for texture modeling.

For our new texture model, the set of auxiliary variables is represented by the scale parameters s, so that the conditional law f (x|s, λ) is Gaussian, but the marginal law f (x|λ) is no longer Gaussian. In this case, the set s contains one variable per Fourier coefficient. Let this model be called a Gaussian Scale Mixture Random Field (GSMRF).

As previously stated, our focus is on stationary textures, which, under a Whittle approximation, have decorrelated Fourier coefficients. For the GSMRF model, due to the Gaussianity of the conditionals, these coefficients are independent conditionally on the scale parameters. At this point, there are two modeling alternatives, depending on the type of prior used for the scale parameters:

• independent scale parameters -as in the case of the GRF model, the Fourier coefficients will be both conditionally and marginally independent. Although this reduces the representation capabilities of the model, it renders it very computational efficient.

• interdependent scale parameters -the texture coefficients are no longer marginally independent. Since the number of scale parameters is equal to the number of pixels, the use of an interdependent set could imply a prohibitive processing cost for the synthesis of such a texture. In this case, although the texture coefficients can be sampled in parallel, conditionally on the scale parameters, the scale parameters themselves must be sampled sequentially in a costly process.

GSMRF with independent scale parameters

This is the simplest modeling alternative and ensures the maximum model tractability.

Remark:

The terminology of marginal law as opposed to conditional law, in the context presented here, only refers to the dependency with respect to the scale parameters s p . All the laws are however conditional with respect to the shape components λ p . Thus, the Fourier coefficient

•

x p has the conditional law f ( • x p |s p , λ p ) and the marginal law f (

• x p |λ p ).
In fact, between the two texture models only the form of χ p changes:

• GSMRF -χ p = s p λ p • GRF -χ p = sλ p
This difference can be induced through the pdf of the scale parameters. Consequently, the form of the prior for s will make the difference between the GRF and the GSMRF models and through these parameters we will be able to switch between the two laws, without changing the PSD.

For the same conjugacy considerations evoked in the previous section, the prior distribution of each s p is chosen to have a Gamma form. Then, the s p are iid:

f (s p |α s , β s ) ∝ s αs-1 p exp [-β s s p ] (3.6)
Based on the previously defined prior for the scale parameters, the marginal law f (x|λ) has a Student's t form:

f ( • x p |λ p , α s , β s ) = sp f ( • x p |s p , λ p ) • π(s p |α s , β s ) ds p ∝ 1 + λ p β s | • x p | 2 -αs-1 (3.7)
Moreover, when comparing the GRF and the GSMRF models, we considered it useful to be able to switch between models without inducing other changes except for the probability distribution of the coefficients. In order to evaluate the influence of this change of law on the texture characteristics, we can impose the condition of having identical second order moments for both cases, i.e., having the same marginal variance of the Fourier coefficients. This allows us to make one change at a time: for instance, a change of the coefficients pdf will not affect their variance, and, conversely, a change in the PSD will not change the type of pdf. Thus, for the two types of laws, we want to have the same variance for the Fourier coefficients, conditionally on the λ p :

var GSM RF • x p |λ p = β s 2α s -1 • 1 λ p var GRF • x p |λ p = 1 λ p
Chapter 3. Texture Modeling this imposing a constraint on the parameters of the prior for s p : β s = 2α s -1, with α s > 0.5. In this manner, the two types of textures will have the same PSD, but different laws.

To summarize, in this setting, the texture synthesis will be a relatively effortless process, consisting in two sampling stages: 1. the s p : independent, with Gamma pdf, 2. the • x p |s p : independent, with Gaussian pdf.

Realizations of the GRF and GSMRF models, with the same PSDs, is given in Figure 3.6. On the first column there are realizations of the GRF texture model with a certain PSD, while on the second column are the realizations of the corresponding GSMRF models, with the same PSDs. Due to their common features:

• independent Fourier coefficients,

• second order statistics,

• the same PSD, the GRF and the GSMRF models yield rather similar, primarily stochastic textures. However, the GSMRF is able to generate slightly more complex patterns. This is due to the fact that, although both models have the same PSD, the equality holds only under the expectation. For one sample, the frequency content of the GSMRF textures will be more complex due to the scale parameters.

GSMRF with dependent auxiliary variables

The choice of using independent scale parameters is a limitation for the model and it is meant to ensure that the model is tractable and the texture synthesis does not become prohibitively expensive. There are however options for defining a certain dependence among the scale parameters without dramatically increasing the complexity.

One of these options is to model the scale parameters by a Potts field (see Appendix C for a description of this model). If the cliques are limited on the 4-neighborhood, the s field can be sampled in only two steps, by sampling in parallel all the variables not being neighbors. If the size of the neighborhood increases, so will the number of steps.

However, the use of such a prior for the auxiliary variables may not bring a considerable change for the texture characteristics. This is illustrated by Figure 3.7 where we show two texture realizations of the GSMRF model, the first corresponding to a Potts prior for the scale parameters and the second to a Gamma separable prior. Surprisingly, although the scale parameters fields are significantly different, the corresponding periodograms are very similar. This is reflected in the textures themselves, which are practically indistinguishable.

In the light of this observation, the use of the Gamma prior for the scale parameters is not only justified by its ease of implementation, but also by the fact that more complicated, non-separable models might not bring any gain in terms of representation capability.

Conclusion and Perspectives

This chapter has presented the mathematical formulation and a series of realizations of our SMGRF texture model. This rather simple model is formulated in the Fourier domain and relies on the principle of conditional Gaussianity with respect to a set of auxiliary variables which renders the marginal law non-Gaussian. Due to this property, this model is very tractable, which makes it fit the constraints of the inverse problems in which it will be integrated.

Despite the series of choices that have been done to achieve this easy to handle form, the model has good representation capabilities, especially for stochastic textures. Moreover, some of these limitations may not even have a strong impact on the model's representation capabilities, as shown in the previous section, where the use of a non-separable prior for the scale parameters did not enhance the representation capabilities of the model and the texture realizations were indistinguishable from those obtained using the separable prior.

Nevertheless, more complex priors can be explored and this is the main perspective of this work on texture modeling. Chapter 2 has presented a complex inversion paradigm where the data y is a blurred and noisy version of the original image. This original image is made up of K regions, each region containing a different texture x k , k = 1...K. The hidden label variables z indicate for each pixel the region to which it belongs.

In Chapter 3 we have provided a detailed description of our GRF based texture models and their properties, which make them very suited to be used in the inversion problem. Consequently, we consider that the original images x k are modeled by GRFs.

The problem becomes estimating the original textured images, determining the underlying model M k and its parameters for each image and the hidden label variables z in an unsupervised manner. Due to the large number of the unknowns and the multiple layers of difficulty, we have decided to divide the problem into a series of subproblems of lower complexity. Notwithstanding, even these subproblems have a high level of complexity and, moreover, are of significant interest since they address important topics in image processing.

Problem Statement

One of the subproblems that can be considered is the unsupervised myopic deconvolution of a textured image with PSD parameter estimation. The term unsupervised refers to the fact that the hyperparameters γ n (noise precision) and γ x (global scale parameter for the PSD of the textured image) are also estimated, while the term myopic indicates that the PSF has a parametric form. This topic has made the subject of our publication [START_REF]Bayesian Texture and Instrument Parameter Estimation From Blurred and Noisy Images Using MCMC[END_REF]. The particularities of this problem are the following:

• original image x consisting in a single texture. The texture is modeled by a GRF (no auxiliary variables s), with known PSD model, driven by the unknown parameters θ,

• the PSF has a known parametric model, driven by the unknown parameters η. In this setting, the dependencies between the data and the unknowns can be graphically represented as shown by the graph in Figure 4.1.

y η γ n α n , β n η m , η M θ m , θ M α x , β x x θ γ x η θ x γ x γ n
The goal is to estimate the PSD and PSF parameters, θ and η, respectively. We have devised an unsupervised method that provides optimal estimates for these parameters. Therefore, our method provides, in addition to the original textured image, estimates for the: 1. texture parameters, θ, 2. instrument parameters, η, 3. signal and noise levels, γ x , γ n . This distinguishes our method from the existing works. These methods also perform image deconvolution and denoising by formulating regularization terms and estimating the noise parameter, however, they do not use parametric regularization terms.

The problem formalization, in the context of the textured images presented in Chapter 3, can be represented as:

x (θ,γx) H η + n (γ n )
y where x (θ,γx) is the textured image and H η is the parametric PSF, having a circulant-block circulant structure, due to the Whittle approximation detailed in Chapter 2. Let us denote by Ψ = {θ, η, γ x , γ n } the full set of unknowns.

The observation model is convolutional and can be described by the equation:

y = H η • x (γx,θ) + n (γ n ) (4.1)
x and n are modeled by zero-mean, stationary GRFs, independent, of covariance matrices R x (θ) and R n (γ n ). From the fields stationarity and the Whittle approximation, these covariance matrices are circulant-block circulant and are thus diagonalizable by Discrete Fourier transform (DFT). Hence, the Fourier coefficients of the data are independent:

• y p = • h p • • x p + • n p (4.2)
and can be computed by Fast Fourier Transform (FFT).

In order to fully exploit the laws separability in the mathematical developments, the PSF, the noise and texture PSDs are written in the frequency domain. From this point forward we will be dealing with the frequency domain counterpart of the PSF, the Transfer Function (TF).

Deconvolution

The literature on image deconvolution is vast and covers various regularization forms and various representation domains for the image, either in a deterministic or in a probabilistic approach. Furthermore, from the PSF model point of view, there are two concurrent approaches to the deconvolution problem, each with its advantages and its weaknesses.

Blind deconvolution

Blind deconvolution is the classic approach to image deconvolution, most of the existing works in the literature being devoted to this type of modeling. It consists in considering the PSF as having a non-parametric model, represented by an impulse response matrix H with unknown elements. Since the convolution is a linear operation, the dependency of the data with respect to the elements of this matrix is also linear. Hence, the estimation process in this context consists in determining the PSF elements, from a linear system of equations.

The problem is still ill-posed, since the number of unknowns is higher than the number of observations and consequently regularization is used in order to achieve a unique solution. In this framework, the regularization is done on the image and on the PSF and is based on the idea that in most cases, the images are piecewise smooth, with a reduced number of discontinuities (contours). Thusly, an adapted regularization is to impose a low frequency model with sparse transitions for the image and a low frequency model for the PSF. This is achieved by imposing various constraints through either stochastic models such as Simultaneous Auto-Regressions, Student's t, etc., or deterministic penalties such as L1 norm, L2-L1 norm, L2 norm, Total Variation. Such regularization forms for the image are limited in terms of representation capability and make the methods appropriate mostly for dealing with piecewise-smooth or piecewise-constant images.

The regularization terms can be formulated in a fully deterministic manner, as is the case in [FBD10, [START_REF] Bioucas-Dias | Multiplicative Noise Removal Using Variable Splitting and Constrained Optimization[END_REF][START_REF] Afonso | An Augmented Lagrangian Approach to the Constrained Optimization Formulation of Imaging Inverse Problems[END_REF], where the solution is determined by convex optimization. The alternative approach is the Bayesian formulation, used in works such as [MMK06, BMK09, BWMK10, TLG07, TLG09], where the aforementioned regularization forms are introduced as priors on the image and the PSF. In both of these cases, the image regularization can be formulated in various representation domains, the vast majority of the literature being devoted either to the spatial representation [MMK06, BMK09, BWMK10] or to the wavelet domain [VU08, TLG07, TLG09, FBD10, BDF10, ABDF11].

The aforementioned methods all have a high degree of complexity and do not not allow for an explicit analytical solution. For this reason, numerical methods must be employed in order to solve these problems. Examples of such methods are for instance the algorithms of Alternate Method of Multipliers for the deterministic cases, or the Variational Bayes for the probabilistic approaches.

The advantage of the blind deconvolution method consists in the linearity of the criterion with respect to the elements of the PSF. This meaning that conjugate priors can be employed and, consequently, that the posterior law for these elements has a standard form.

[CE07] offers a detailed review of the blind deconvolution techniques and their applications.

Myopic deconvolution

Myopic, or semi-blind, deconvolution is the alternative to the blind approach. It is based on a parametric model for the PSF, driven by a set of parameters of reduced size. In practice, information regarding the PSF is often available, especially about its form. Moreover, it is not uncommon for this form to be parametric.

The non-linearity of the PSF with respect to its coefficients is for instance typical in astronomy. The problem of myopic deconvolution for astronomical images is tackled in [START_REF] Orieux | Bayesian estimation of regularization and Point Spread Function parameters for Wiener-Hunt deconvolution[END_REF] in a Bayesian framework using an L2-L1 regularization on the image.

In another line of work, [PDH12, PDH14] address the problem of myopic deonvolution from a different perspective, by overcoming the non-linearity of the PDF through its decomposition onto a basis. At a first glance, myopic deconvolution might seem like a simpler problem as compared to the blind deconvolution case. Indeed, the number of unknowns becomes smaller, however, usually, the form of the dependency with respect to the parameters of the PSF is no longer linear and in most cases is rather complicated. Consequently, there will be no conjugate priors to enable a straightforward sampling procedure. In such a situation, adapted sampling algorithms must be employed, the downside being that this can be rather costly.

In other words, the main difference between the blind and myopic approaches is the manner in which they constrain an ill-posed problem. For the blind case, the supplementary information needed to constrain the solution is represented by the regularization term over the PSF coefficients, while for the myopic case the parametric model itself is strongly structuring.

Data, image and TF models

The previous section has presented the two approaches for image deconvolution, blind and myopic. In the following we will present our method based on the latter approach, we will provide the models we have used and the corresponding mathematical developments.

From the noise law and (4.2) follows the relation for the Fourier coefficients of the observations

• y, texture • x, filter's TF • h and hyperparameters γ x , γ n : f ( • y p | • x p , η, γ n ) ∝ µ p (γ n ) • exp -µ p (γ n ) • y p - • h p (η) • x p 2 (4.3) f ( • x p |θ, γ x ) ∝ γ x λ p (θ) • exp -γ x λ p (θ) • x p 2 (4.4)
where µ p are the eigenvalues of R -1 n . The • h p (η) and λ p (θ) from (4.3) and (4.4), respectively, represent the TF values and the inverse PSD, respectively, computed on the discretized frequency domain at position p. λ (ν x , ν y , θ) represents the inverse of the texture's PSD. (ν x , ν y ) belongs to the reduced frequency domain and (ν n , ν m ) = (n∆ν, m∆ν) is pixel's p location.

An interesting aspect is the fact that the method is adapted to any TF and any PSD form for the texture and the noise. However, to illustrate the mathematical developments and numerically evaluate the algorithm, we have chosen:

• TF -low pass filter (Dirichlet kernel), of width η x = η y = η:

• h(ν x , ν y , η) = 1 2η sin 2πη ν 2 x + ν 2 y sin π ν 2 x + ν 2 y (4.5)
• noise -white noise of covariance R n (γ n ) = γ -1 n I, with γ n inverse variance / precision parameter.

• image -exponential model for the PSD:

λ -1 (ν x , ν y , γ x , θ) = exp - |ν x -ν 0 x | u x + |ν y -ν 0 y | u y (4.6)
with θ = ν 0 x , ν 0 y , u x , u y ∈ R 4 . ν 0 x , ν 0 y are the central frequencies and u x , u y are the PSD widths. λ p (θ) from (4.4) is in fact obtained as λ p (θ) = λ(ν n , ν m , θ).

Myopic deconvolution for textured images

As seen in (4.5) and (4.6), the TF and PSD are driven by the unknown parameter sets η and θ. The coefficients of both characteristics have a highly non-linear dependency with respect to the parameters, i.e., the estimation process is rather challenging.

Information and qualitative estimation performance analysis

The parameter estimation performances are directly related to the amount of available information regarding each parameter. Consequently, the Fisher information is extremely useful in evaluating the different situations that may occur, in function of the SNR and the parameter values. Different scenarios are illustrated by Figure 4.2, as 1D cross-sections of the 2D frequency domain:

(a) for narrow TFs (small η) and high frequency PSDs (large (ν 0

x , ν 0 y )), the spectral contents cancel each other, i.e., the data is informative on γ n , but not on η and θ;

(b) for wide TFs and narrow PSDs, the input stimulus is incapable to induce an adequate system perturbation, i.e., the information is insufficient to estimate η;

(c) ideal situation information-wise, i.e., partial overlap. The information available for the estimation depends on this overlap.

Appendix A presents the analytical expressions for the Fisher information for each of the parameters. Based on these analytical forms, we can provide a detailed qualitative and quantitative analysis of the mean available information concerning each parameter of the problem.

The hyperparameters γ n and γ x have an important role in the estimation process since their ratio γ n /γ x represents the Signal to Noise Ratio (SNR). A high SNR means that the signal level is high enough so that the noise corruption cannot mask the useful information. Thus, in such a case, the estimation performances for all the parameters are high. The only exception is the noise level itself, which is better estimated when the SNR is low.

The following subsections provide detailed considerations concerning the available information about the noise parameter, a PSD central frequency and a PSD width. The variation of this quantity is taken with respect to one parameter at a time, in order to show the strong dependencies and to identify the invariant situations.

Noise parameter γ n

The analytical expression of I (γ n ) is:

I (γ n ) = γ -2 n p 1 + γ n /γ x • g p (η) s p λ p (θ) -2 (4.7) where g p (η) = | • h p (η)| 2 Figure 4.
3 illustrates the variation of the amount of mean information concerning the noise parameter. On the x axis we have represented a very wide range of possible values for the noise, since no prior information is available about this parameter. The representation is done in logarithmic scale for both axes.

The SNR sweep, which is equivalent to modifying the γ x , shows that an SNR increase triggers a decrease of I (γ n ), i.e., the estimation of the noise parameter will be more difficult. With respect to the widths of the TF and the PSD, η and u, respectively, the dependency is similar, that is the higher the width, the less the information on γ n .

One of the most interesting considerations is related to the relative positioning of the TF and the PSD. Figure 4.3d illustrates the Fisher information for the three cases previously presented in Figure 4.2. It is interesting to notice that while the case of total overlap coincides with a lower information, there is practically no difference between the cases with no or partial overlap. This means that even in a case with no overlap, γ n should be well estimated.

A more quantitative-oriented analysis can also be made. For this purpose, let us take a closer look at Figure 4.3a, for instance. For γ n = 10 -2 , the corresponding Fisher information is I (10 -2 ) = 4 • 10 7 , while I (10 -1 ) = 4 • 10 5 , meaning that the gain of an order of magnitude for the noise level implies the gain of two orders of magnitude in terms of information amount. Remark: A general trend can be noticed in Figure 4.3 and it consists in a higher sensitivity to the parameter variation for lower noise levels (high γ n ), while for the high noise levels this is almost indistinguishable. This means that in cases where the noise level is high enough, its estimation will be invariant to the rest of the factors involved.

PSD central frequency ν 0 x

The Fisher information for ν 0

x , in the case of a Laplace shaped PSD, writes:

I (ν 0 x ) = 1 u 2 x p g p (η)/λ p (θ) γ x /γ n + g p (η)/λ p (θ) 2 (4.8)
The Fisher information variation for the texture parameters θ is illustrated by Figure 4.4. On the x axis it is represented the range of possible values for the central frequency, i.e., the reduced frequency domain. The representation is done in logarithmic scale on the y axis.

The plot representing the SNR influence clearly shows that the amount of information is strongly influenced by the SNR. Thusly, we may assert that the estimation performances should be significantly superior in the cases with higher SNR and the performance degradation when the SNR decreases should be rather significant.

Figure 4.4b shows the filter width influence and gives valuable insight for the entire information analysis. First of all, no matter the value of η, the Fisher information has very similar values around the null frequency. The difference intervenes at higher central frequencies. The reason for this is the presence of the blurring filter, which has a lowpass character. This means that, irrespectively of the value of η, if ν 0 is close to the null frequency, it will be in the filter bandpass and thus the amount of information is high. However, if ν 0 is high, the filter width becomes a factor and we notice a decrease of the information for small values of η. This is confirmed by the case where the filter TF is very wide, corresponding to a situation where there is almost no blurring. In this case, I (ν 0 ) is almost constant over the entire interval, as it can be seen in Figure 4.4b.

Remark:

The previous considerations apply to the variation with respect to the other pa- rameters as well and explain this generalized behavior with a maximum in the null frequency and local minima in the frequencies where the filter TF is zero.

The corresponding PSD width u x also influences I (ν 0

x ), a larger width meaning larger uncertainty and thus less information about the position. Nevertheless, the width u y has a very reduced impact on I (ν 0

x ). The central frequency ν 0 y intervenes in the expression of I (ν 0 x ) through the g p (η)/λ p (θ) ratio and its influence is due to the presence of the filter. Consequently, the dependency of I (ν 0

x ) will have an oscillating character, in function of the PSD positioning with respect to the local maxima or zeros of the filter TF.

PSD width u x

The Fisher information for u x , for the Laplace shaped PSD, has the following expression:

I (u x ) = 1 u 4 x p ν x p -ν 0 x g p (η)/λ p (θ) γ x /γ n + g p (η)/λ p (θ) 2 (4.9)
Figure 4.5 illustrates the variation of I (u x ) with respect to some of the variables that appear in (4.9). The x axis represents a wide range of positive values, since we do not have any prior information about this parameter. The plot is done in logarithmic scale for both axes.

The shape of I (u x ) itself, when all the other parameters are fixed, is hard to explain, due to the cumbersome dependency with respect to u x , which also intervenes in the terms λ p (θ).

The SNR sweep shows that a high SNR corresponds to a higher amount of information. Figure 4.5b shows that I (u x ) does not vary significantly with the width of the filter, in a case where ν 0

x is close to the null frequency (and implicitly lies in the filter passband). Figure 4.5c gives some insight concerning the information available in the various PSD and TF relative positioning cases. This plot shows that there will be two regimes: for very narrow PSDs the most advantageous situation information-wise is the case of total overlap (which implies a high η). In the second regime (roughly when u x > 10 -3 ) I (u x ) is significantly higher in the case of a partial overlap. As for the other two situations, when u x < η the information is higher in the total overlap case. Last, but not least, when u x > η, the case of total overlap becomes the more challenging information-wise. *** A quantitative analysis of the Fisher information for the texture and TF parameters is harder to make as compared to the signal and noise levels. This is due to the more complex dependency, with multiple levels of non-linearity and interconnections. Despite this difficulty, we have provided an analysis that will aid in anticipating and understanding the method's performances.

Bayesian setting: priors, posterior and conditional posteriors

From the hierarchical conditioning of the variables, shown in Figure 4.1, the joint law for y, x and Ψ writes:

π(y, x, Ψ) = f (y|x, η, γ n ) • f (x|θ, γ x ) • π(Ψ) (4.10)
where π(Ψ) represents the joint prior for the unknowns.

Using Bayes rule and taking into account that f (y) is constant with respect to the unknowns, the posterior is proportional to the joint law (4.10):

π(x, Ψ|y) ∝γ P n exp -γ n P p=1 • y p - • x p • h p (η) 2 • γ P x P p=1 λ p (θ) • exp -γ x P p=1 | • x p | 2 λ p (θ) • π(Ψ) (4.11)
Due to the small amount of available information, there is no indication of a dependency between the priors of the variables, thus we can consider that the priors are independent of each other:

π(Ψ) = π(θ) • π(η) • π(γ x ) • π(γ n ) (4.12)
The forms for the priors can be judiciously chosen by analyzing (4.11). For γ n the uninformative Jeffreys prior will be used, as explained in Section 2.2. It can be noticed that γ x intervenes as a precision parameter in a Gaussian law and that in this case the Gamma law is the conjugate form:

π(γ x |α x , β x ) ∝ γ αx-1 x exp (-β x γ x ) (4.13)
Since there is no prior information concerning the values of the signal level, we should use an uninformative prior. Similarly as for γ n a Jeffreys prior is obtained by setting (α, β) → (0, 0). For θ and η, the complicated dependency means that there is no conjugate form that can be used. Moreover, since there is no prior information regarding the values of these parameters, uninformative priors should also be used for the texture and instrument parameters. Consequently, uniform priors, defined over all the possible range of values, will be used for each parameter:

π(θ) = U [θm,θ M ] (θ) π(η) = U [ηm,η M ] (η) (4.14)
Using the previously mentioned priors results in the following forms for the conditionals a posteriori:

x ∼ exp - P p=1 γ n | • y p - • h p (η) • x p | 2 + γ x | • x p | 2 λ p (θ)
f (x|y, γ n , θ) -separable and quadratic in the Fourier domain;

γ n ∼ γ P +αn-1 n • exp -γ n β n + P p=1 | • y p - • h p (η) • x p | 2 f (γ n |y, x, α n , β n ) -Gamma form; γ x ∼ γ P +αx-1 x • exp -γ x β x + P p=1 | • x p | 2 λ p (θ) f (γ x |x, θ, α x , β x ) -Gamma form; θ ∼ P p=1 λ p (θ) • exp -γ x | • x p | 2 λ p (θ) • U (θ)
f (θ|x, γ x ) -independent on the observations, but has a very complicated dependency;

η ∼ P p=1 exp -γ n | • y p - • h p (η) • x p | 2 • U (η)
f (η|y, x, γ n ) -very complicated dependency.

Remark:

The conditional law of the image with respect to the rest of the unknowns is closely related to Wiener filtering. In fact, the maximum of this law is the signal obtained through Wiener filtering. In the special case with no noise γ n → ∞ the maximum is obtained by applying the inverse filter.

In order to determine the estimates the PM estimator is used, which is MSE optimal, as detailed in Appendix B. The estimates are obtained by calculating the integral ψ|y ψπ(ψ|y)dψ. The parameter dependency is complicated, making the integral intractable. In order to overcome this difficulty, and give the complicated laws for θ and η, the Metropolis within Gibbs sampling technique will be employed to approximate the integral. Among the multiple versions for the MH step, the RWMH is chosen due to its simplicity and ease of implementation.

Results

As stated in the previous section, the iterative algorithm consists in sampling the parameters Ψ and x and averaging the samples to obtain the PM estimates. The algorithm associated to our method is given in Algorithm 1. We consider this algorithm has converged once all the recursive means variations from one iteration to the next are less than 0.1% of the nominal parameter value.

Algorithm 1: Unsupervised Myopic Deconvolution for Textured Images Algorithm input : Data y output: Estimates for the texture parameters θ, instrument parameters η, hyperparameters γ n and γ x . In addition, samples for the texture x (t)

% generate samples of π(θ, η, γ n , γ x , x|y) ; t = 1; initialization θ (t) ,η (t) , x (t) = y;

while |variation| > ε do t = t + 1; γ (t) n ∼ f (γ n |y, x (t-1) , h p (η (t-1) ), α n , β n ) γ (t) x ∼ f (γ x |x (t-1) , θ (t-1) , α x , β x ) θ (t) -RWMH with target f (θ|x, γ x ) η (t) -RWMH with target f (η|y, x, γ n ) x (t) ∼ f (x|y, θ (t) , η (t) , γ (t) x , γ (t) n )
end % Compute the parameter estimates by PM; % let BI be the iteration where we consider the burn in period is over;

θ = PM θ (BI...T ) η = PM η (BI...T ) γn = PM γ (BI...T ) n γx = PM γ (BI...T )
x Table 4.1 lists the results of our estimation method, expressed in percentages (%). These results are a mean relative error on 20 realizations of each scenario for different parameter values. The SNR = 10 log γ n /γ x represents the original signal to noise ratio. However, the Blurred SNR (used in [START_REF] Babacan | Variational Bayesian Blind Deconvolution Using a Total Variation Prior[END_REF] to quantify the problem difficulty) is significantly smaller and depends on the positioning of the TF and PSD.

The data is more informative for some parameters, as shown by our Fisher information analysis in Appendix A and Section 4.2.1. This explains the difference of estimation er- rors that can be observed in Table 4.1. The overall estimation performance for the central frequencies ν 0 x , ν 0 y are superior to those for the widths {u x , u y } and η. This is in direct correlation with the Fisher information shown in Figures 4.4 and4.5. The Fisher information for the central frequencies has a relatively high value over the entire range of possible values, while in the case of the widths, it has a very high value for a range of values and it drops dramatically for the rest. The noise level is also an very important parameter, the accuracy of its estimation having a direct impact on the rest of the parameters. The higher noise levels are more accurately estimated, but the estimation of the rest of the parameters is very difficult in low SNR cases. From the signal level point of view, a high value means high SNR and triggers an accurate estimation for all the parameters, except for the noise level. Overall, a higher sensitivity to the noise than to the convolution can be noticed and this is clearly indicated by the estimation errors in Table 4.1.

We must stress the fact that despite the difficult cases, globally, the method yields a small error for all the parameters, proving the fact that it is well adapted for this problem and that it is able to provide accurate estimates even when the amount of available information is very reduced.

Figure 4.6 presents visual results for our algorithm. More specifically, the cases presented in Figure 4.2 are illustrated by showing on the first column the characteristics in the Fourier domain, on the second column, a texture corresponding to each case, x|θ * , on the third column, the degraded observations, y, and finally, on the last column, the reconstructed image x| θ, from which the blur and noise have been eliminated. The SNR = 25dB for all the cases, at this value the amount of information being sufficient to provide both a very good visual quality of the reconstruction and the restoration of the gray level scale. This shows the method's capacity to restore a textured image affected by a significant visual degradation, by correctly estimating the blur parameters and the noise level.

Figure 4.7 illustrates the method sensitivity to noise, by presenting the deconvolution results for the same texture, blurred by the same TF, but with different noise levels. On the first column we show the characteristics in the Fourier domain and the original texture x|θ * . On the following three columns we present successively the observations y (on top) and the deconvolution results x| θ (below). These cases correspond, from left to right to levels of SNR of 15dB, 20dB and 25dB.

From the information point of view, this is the most advantageous case, meaning that we should expect satisfactory results. Indeed, it can be noticed that the deconvolved images all resemble the original texture in what concerns the orientation and the frequency, meaning that the central frequencies of the PSD have been accurately estimated. However, for the cases with SNR = 15dB and SNR = 20dB, the textural content is very similar, but richer than the original texture, which indicates that the widths of the PSD have been slightly overestimated, resulting in a wider PSD. The reason for this phenomenon is that the algorithm attempts to adjust the thickness of the PSD tails, to account partly for the noise, in cases with high noise corruption.

Another aspect is related to the gray range, i.e., to the estimation of γ x . In the first two cases, the algorithm has a tendency to slightly overestimate the dynamic range of the observation, which means the γ x is not as accurately estimated in these cases. However, in the SNR 25dB case, γ x is correctly estimated.

Conclusion and Perspectives

The main contribution presented in this chapter is the use of parametric models for the image and for the PSF, in order to tackle the deconvolution problem in a context of textured images with non-linear data dependency with respect to the PSF unknowns. To the best of our knowledge, there is no other method to address this problem. Moreover, by using the optimal PM estimator guarantees that at least from the MSE viewpoint, no other method can provide better results.

We have proposed a detailed analysis of the Fisher information for each unknown, which is directly related to the difficulty of the estimation. This has allowed us to anticipate the performances and then confirm them with numerical results.

Further developments to this problem consist in using more complex models for the PSF, which would be able to describe other types of distortions as well. For instance, this could allow us to address the problem of motion blur.

CHAPTER 5

Model Choice for the Law and the PSD of a Textured Image Chapter 4 has presented a first subproblem in our textured images restoration context. This chapter is devoted to another interesting problem: texture model selection from indirect data. Our main contribution is tackling a problem that has not been dealt with so far by developing a method based on an optimal risk decision.

Model Choice -State of the art

Model choice applications span over a wide range of fields, for instance microbiology, proteomics, economics, statistics, signal and image analysis. The large majority of these methods are based on the log-likelihood value computed for the Maximum Likelihood Estimate (MLE), θ and contain different penalization terms for the model dimension.

The Generalized MLE (GMLE) method for model selection consists in determining the model with the highest likelihood for the MLEs of the parameters. However, this

M Model k Expression of λ k (ν x , ν y , θ k ) θ k 1 Lorentz 1 + (νx-νx 0 ) 2 σ 2 x + (νy-νy 0 ) 2 σ 2 y ν x 0 , ν y 0 , σ x , σ y 2 Generalized 1 + (νx-νx 0 ) 2 σ 2 x + (νy-νy 0 ) 2 σ 2 y q ν x 0 , ν y 0 , σ x , σ y , q Lorentz 3 Exponential exp 1 2 |νx-νx 0 | σx + |νy-νy 0 | σy ν x 0 , ν y 0 , σ x , σ y 4 Generalized exp 1 2 |νx-νx 0 | q σ q x +
|νy-νy 0 | q σ q y ν x 0 , ν y 0 , σ x , σ y , q Gauss Table 5.1: PSD models, the expression of the λ coefficients and the corresponding parameters.

One of the most important features of any efficient model selection method is the ability to penalize the model complexity, more specifically to be capable of selecting the less complex model for similar measures of model fit. In order to test the method's performances from the complexity penalization point of view we have chosen to include embedded models in our dictionary of shapes for the PSD.

The texture PSD models are built in the reduced frequency domain. There is virtually no constraint on the form, except for positivity: it can be constant, corresponding to a white noise, or even a single frequency field. Here, the focus is on parametric, unimodal functions, with a relatively reduced number of parameters.

The chosen parametric shapes are: Lorentzian (M = 1), Generalized Lorentzian (M = 2), Exponential (M = 3), and Generalized Gaussian (M = 4). Each model k is driven by the corresponding θ k . Table 5.1 shows explicitly the dependencies λ k (θ k ). The parameters of these parametric forms are: the central frequencies ν x0 , ν y0 and the widths σ x , σ y . q represents the power parameter and is specific to models M = 2 and M = 4.

We have chosen to use embedded models, for instance M = 1 and M = 3 are nested in M = 2 and M = 4, respectively. This enables an analysis on the method's capacity to penalize model complexity when the extra parameters do not trigger a significant model fit increase.

The textured images are spatially discrete, thus the PSD is defined on the reduced frequency domain, i.e., the variables (ν x , ν y ) ∈ [-0.5, 0.5] 2 . Furthermore, let us consider that the Fourier coefficient p has the (ν m , ν n ) position in the discrete reduced frequency domain. The λ p from (5.4) are the elements of the PSD field at these discrete positions and depend on the model, M = k, and on the parameters θ k . From this point forward, to explicitly show this dependency, the notation λ k p (θ k ) is used and, more precisely:

λ k p (θ k ) = λ k (ν m , ν n , θ k ).
The use of parametric models has the advantage of reducing the number of unknowns. On the other hand, this model defines a highly non-linear dependency of λ k p (θ k ) with respect to θ k , as shown in Table 5.1. This complicated dependency means that there is no Moreover, the prior information about the parameters is very reduced, thus, uninformative priors will be used. Consequently, a uniform prior is employed:

π(θ k |M = k) = U [θ m k ,θ M k ] (θ k ).

Probabilistic model choice

This model choice problem is twofold and refers to:

1. the form of the law for the texture coefficients, i.e., whether we are dealing with a GRF or a SMGRF, 2. the parametric form of the PSD.

The difficulty of this problem is mainly due to the fact that the image models (presented in Chapter 3) are highly non-linear and have a very complicated dependency with respect to the parameters. Moreover, the parameters driving these models are unknown. Consequently, in order to infer on the models, we also have to infer on their values. Therefore, this is far from trivial even in the case of direct observations. Furthermore, the fact that we are dealing with blurred and noisy observations adds a supplementary layer of complexity.

The indirect observations aspect is not circumvented by the fact that the PSF is known, since the noise level γ n and the signal level γ are unknown. This is all the more challenging especially in our textured image context, where the blur and noise corruption yield a new textured image. This new image can be mistakenly considered either as being a realization of the same texture model as the original image, but with different parameter values θ, to account for the blur and noise, or even of a different PSD model, with thicker tails, to account for high levels of noise. Consequently, the additional difficulty consists in not being mislead by the distortion and correctly distinguishing the true model M * , the true parameter values θ * and the true hyperparameter values γ * n and γ * . Figure 5.1 shows the variable conditioning for this problem. The red singles out the texture model, which is in fact the quantity of interest. In gray, we have represented the rest of the unknowns, whose estimation is not our primary goal. However, our model choice algorithm will provide samples for these unknowns (conditionally on the model), as an additional, auxiliary result.

The detailed presentation of this problem and the method we propose for solving it has made the subject of our paper [VGB].

This method relies on a probabilistic framework. Using Bayes' rule, the posterior model probabilities are:

Pr(M = k|y) = f (y|M = k) • p k f (y) (5.1)
and require the computation of two quantities.

i. The probability distribution of the data, f (y). Fortunately, it does not depend on the model, thus can be calculated by normalization.

ii. The evidence, e k = f (y|M = k), obtained from the joint law of data and unknowns, given the model, by marginalizing the unobserved texture, the noise precision and the texture parameters:

e k = Ψ f (y, Ψ|M = k)dΨ (5.2)
where Ψ represents all the unknowns.

The optimal classifier is built in a Bayesian framework, where each model's posterior probability is determined from the model evidences. These evidences are intractable and thus are numerically computed by MCMC methods.

From an algorithmic viewpoint, this work embeds the FMH algorithm introduced in [START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF] and presented in Section 2.4.5. Consequently, an important performance increase is achieved, due to the specific nature of our problem. In our case, there is no need for computing second order derivatives and thus a Newton type proposal is built using only first order derivatives.

Joint law and priors for the model, image and parameters

Our model selection method relies on a probabilistic formulation of the problem in a Bayesian framework and determines the a posteriori probability of each model Pr(M = k|y). The a priori distribution for the model is fully described by the p k = Pr(M = k) probabilities. In our numerical study, we have employed an uninformative prior, i.e., a priori equiprobable models:

p k = 1/K, k = 1...K.
The parametric texture models are driven by the parameter set χ. Thus, for each model, k, we have the corresponding parameter set χ k and the law f (x|χ k , M = k).

Let Ψ = {γ n , x, χ k } represent all the unknowns. The joint law is written using the conditioning rule and the hierarchy shown in Figure 5.1:

f (y, Ψ|M = k) = f (y|x, γ n ) • f (x|χ k , M = k) • π(χ k |M = k) • π(γ n ) (5.3)
The f (x|χ k , M = k) represents the law for the texture, driven by the parameters χ k = {γ k , s k , θ k } and the model index M . Then, the expression of this law is:

f (x|γ k , s k , θ k , M = k) = γ P x P p=1 s kp λ k p (θ k ) • exp - P p=1 γ x s kp λ k p (θ k )| • x p | 2 (5.4)
The priors π(γ k ), π(s k ) and π(θ k |M = k) are explicitly given in Section 5.2.1, while π(γ n ) is specified in Section 2.2.

Evidence calculation

The full description of our Bayesian model choice method relies (i) on the data model and (ii) on the priors for the unknowns. The hierarchical direct model and the texture model specificities are shown in Figure 5.1. In the following, the emphasis is on the SMGRF, which encompasses the GRF. By writing (5.2) as:

e k = x,γ k ,γn,θ k ,s k f (y|Ψ, M = k) • π(Ψ|M = k)ds k dθ k dγ k dγ n dx
(5.5) and plugging in (5.11) the intractability of the integral is obvious. For this reason, it must be calculated numerically and the solution chosen here is sampling.

A natural idea is to straightforwardly determine the evidence e k from samples of the prior π(Ψ|M = k):

Ψ (t) = x (t) , γ (t) n , γ (t) k , s (t) k , θ (t) k , t = 1...T as follows: ēk = 1 T T t=1 f y|Ψ (t) , M = k (5.6)
It consists in sampling the priors and computing the arithmetic mean of the corresponding likelihood values -Arithmetic Mean Approximation (AMA). Evidence computation based on prior samples can also be done by nested sampling [START_REF] Skilling | Nested Sampling for General Bayesian Computation[END_REF].

Nevertheless, when the likelihood is very peaked, as in the current case, most of these samples will have weak likelihood, i.e., an insignificant contribution and thus, the algorithm is slow to converge. For this reason, it is more suitable to compute (5.5) by sampling the posterior π(Ψ, M = k|y).

Evidence approximations based on posterior samples

Our model selection method is based on evidence approximation from posterior samples and the proposed method will be referred to as the Classifier based on Evidence

Evidence calculation 59

Approximation from Posterior Samples (CEAPS). This method can be formulated using two different approximations for the evidence, presented in the following: the Harmonic Mean Approximation (HMA) [START_REF] Newton | Approximate Bayesian Inference with the Weighted Likelihood Bootstrap[END_REF] and the Laplace-Metropolis Approximation (LMA) [START_REF] Raftery | Hypothesis Testing and Model Selection via Posterior Simulation, Practical Markov Chain Monte Carlo[END_REF].

Harmonic Mean Approximation

Let us consider that

Ψ (t) = x (t) , γ (t) n , γ (t) k , s (t) k , θ (t) k
with t = 1...T are samples from the a posteriori law. Then, the evidence can be computed as:

ẽk ≈ 1 T T t=1 f (y|Ψ (t) , M = k) -1 -1
(5.7)

i.e., the harmonic mean of the likelihood values for the samples Ψ (t) .

Although ẽk converges almost surely to the true value e k when T → ∞ [Pak99], it does not generally satisfy the central limit theorem [START_REF] Robert | Computational methods for Bayesian model choice[END_REF]. Occasionally, a Ψ (t) with significant a priori probability, but very low likelihood, may occur. Its contribution in the harmonic mean is high and this may trigger infinite variances [START_REF] Newton | Approximate Bayesian Inference with the Weighted Likelihood Bootstrap[END_REF]. Solutions to stabilize this estimator have been provided in [START_REF] Raftery | Estimating the integrated likelihood via posterior simulation using the harmonic mean identity[END_REF]. Nevertheless, we have not encountered this difficulty in neither of our model choice works [START_REF] Vacar | Bayesian Texture Model Selection using Harmonic Mean[END_REF][START_REF]Bayesian Texture Model Choice from Indirect Observations using Fast Sampling[END_REF], where the priors are uninformative on a finite interval and the likelihood is very peaked. In these cases, the posterior samples are distributed in the regions where the likelihood has significant values. Consequently, the situations where the HMA may diverge or converge too slowly are avoided.

Moreover, its computational efficiency and ease of implementation have lead to its use in a series of topic modeling papers such as [START_REF] Griffiths | Finding scientific topics[END_REF][START_REF] Wallach | Topic modeling: beyond bag-of-words[END_REF].

Laplace-Metropolis Approximation

The evidence (5.5) can also be expressed as:

e k = Ψ exp {log [f (y|Ψ, M = k) • π(Ψ|M = k)]} F k (Ψ,y)

dΨ

(5.8) with F k (Ψ, y) the log-posterior computed for observation y.

Remark: F k (Ψ, y) is related to the Fisher information I k and indicates the mean amount of available information, for the given observation y. F k (Ψ, y) is related to the Fisher information I k introduced in Chapter 2 and Appendix A through (A.1) as follows:

I k (Ψ) = -E y [F k (Ψ, y)]
(5.9)

Under the hypothesis that F k (Ψ, y) is twice differentiable with a unique maximum in Ψ * , a Laplace approximation can be applied to evaluate the integral (5.8):

êk ≈ exp [P F k (Ψ * k , y)] • (2π/P ) D k /2 • -F k (Ψ * k , y) -1/2
(5.10)

where Ψ * k represents the MAP value for model k, F k (Ψ * k , y) is the Hessian of the logposterior, evaluated at Ψ * k and D k is the dimension of model k. The last factor in (5.10) is the determinant of the observed information matrix. In fact, computing the evidence in this manner consists in determining the MAP value, Ψ * , i.e., the value for which F (Ψ, y) is maximum, and replacing this value in relation (5.10). [START_REF] Kass | Bayes Factors[END_REF] reviews the Laplace based methods for evidence computation. These approximation methods have relative errors of order O(P -1 ).

This approximation can also be performed based on MCMC samples from the posterior. In this case the method is called Laplace Metropolis Approximation (LMA) [START_REF] Raftery | Hypothesis Testing and Model Selection via Posterior Simulation, Practical Markov Chain Monte Carlo[END_REF]. The LMA can be based in the MAP, the PM, or the MedAP and the Hessian can be approximated by the covariance matrix of the samples. In our case, this approximation is performed using the PM and the value of the Hessian computed for the PM.

Remark: The LMA explicitly penalizes complex models due to the second factor that decreases exponentially with model dimension, while for the HMA the model complexity penalization is implicit, achieved through the likelihood values used to compute the evidence.

Posterior sampling

The posterior law is proportional to the joint law:

f (y, Ψ, M = k) =C • exp -γ n P p=1 | • y p - • h p • x p | 2 • γ P +αn-1 n exp [-β n γ n ] • P p=1 s kp • λ k p (θ k ) • γ P +αx-1 k exp [-β x γ k ] exp -γ k P p=1 | • x p | 2 s kp λ k p (θ k ) • U [θ m k ,θ M k ] (θ k ) • P p=1 s k αs-1 p • exp -β s P p=1 s kp (5.11)
where the normalization constant is

C = K -1 • (2π) -3P • β αn n • Γ -1 (α n ) • β αx x • Γ -1 (α x ) • β P αs s • Γ(α s ) -P • θ M k -θ m k -1 .
However, this law cannot be directly sampled, thus MCMC methods will be employed, more precisely, Gibbs sampling. This can be performed via two types of algorithms.

• Across-model approach -joint sampling of the model index and its parameters. The algorithm jumps from one model to another and explores the joint model index plus parameter space, yielding a joint chain of model indexes and parameter values. The most representative algorithm of this type is Reversible Jump MCMC (RJMCMC) [START_REF] Green | Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination[END_REF].

• Within-model approach -consists in exhaustively visiting the candidate models and parameter sampling conditionally on the model. It provides K chains of parameter values, one for each model. For a detailed description see [START_REF] Gelfand | Bayesian Model Choice: Asymptotics and Exact Calculations[END_REF][START_REF] Newton | Approximate Bayesian Inference with the Weighted Likelihood Bootstrap[END_REF] and the more recent survey [START_REF] Robert | Computational methods for Bayesian model choice[END_REF].

Despite the conceptual differences, for a finite candidate models set, the two approaches yield the same result (provided they have reached convergence) but, under two different forms.

The RJMCMC algorithm is especially interesting for very large numbers of models, when an exhaustive sequential sampling of all the models may be prohibitively expensive. Nevertheless, this algorithm may pose problems when the models are very different. In this case, when skipping from one model to another, a transformation must be applied in order to determine the current values of the parameters for the new model, based on the current parameter values for the old model. This transformation may be difficult, or even impossible to determine analytically. Furthermore, applying an incorrect transformation may trigger high rejection rates and may lead to an inefficient exploration of the modelparameter space.

Since in our problem the number of concurrent models is rather reduced, the withinmodel sampling is the best strategy. This avoids as well the non-trivial RJMCMC problems concerning the parameter transformation when switching models. Moreover, the withinmodel approach guarantees that all models have been thoroughly explored and the model selection is not affected by the sampling algorithm. For this reason, our model choice method is based on within-model posterior sampling.

Gibbs within-model posterior sampling

The samples from the posterior law are obtained using Gibbs sampling. Among the various strategies, we have chosen to sample γ n , γ k , x, θ k and s k . The advantage of this approach is that we obtain rather standard targets and we can perform parallel sampling for x and s. Then, the a posteriori conditional laws for the parameters are:

x ∼ p exp -γ n | • y p - • h p • x p | 2 + γ k | • x p | 2 s kp λ k p (θ k ) = p N (m p , v p ) with m p = γ n • y p • h p v p v p = γ n g p + γ k s kp λ k p (θ k ) -1
where

g p = | • h p | 2 .
-separable in the Fourier domain, i.e., parallel sampling is possible, -computation cost equivalent to sampling the a priori law;

s k ∼ p s k αs p • exp - P p=1 s kp γ x | • x p | 2 λ k p (θ k ) + β s = p G a s , b sp with a s = α s + 1 b sp = β s + γ k | • x p | 2 λ k p (θ k )
-separable, independent on the observations, allowing for parallel sampling;

γ n ∼ γ P +αn-1 n • exp -γ n β n + p | • y p - • h p • x p | 2 = G (a n , b n ) with a n = α n + P b n = β n + p | • y p - • h p • x p | 2 γ k ∼ γ P +αx-1 k • exp -γ k β x + p | • x p | 2 s kp λ k p (θ k ) = G (a x , b x ) with a x = α x + P b x = β x + p | • x p | 2 s kp λ k p (θ k ) θ k ∼ p λ k p (θ k ) • exp -γ k s kp | • x p | 2 λ k p (θ k ) • U (θ k )
-very complicated dependency.

Remark:

The alternative to sampling all the unknowns is to integrate a part of them, in order to avoid certain sampling steps:

• x marginalization and sampling the rest of the unknowns, thus no texture sampling, but even more cumbersome dependency on χ k and γ n . This strategy is similar to the collapsed Gibbs sampler method used in [START_REF] Kail | Blind Deconvolution of Sparse Pulse Sequences under a Minimum Distance Constraint: A Partially Collapsed Gibbs Sampler Method[END_REF]. s remains separable, thus parallel sampling is feasible, but not of standard Gamma laws. Furthermore, the posterior for γ n and γ k no longer have Gamma forms either:

f (s kp , γ n , θ k | * ) ∝ exp 1 r p | • y p | 2 g p γ k s kp λ p (θ k ) • • exp -r p γ n γ k s kp λ p (θ k )   • x p - 1 r p • y * p • h p γ k s kp λ p (θ k )   2
• s k marginalization, resulting in a complicated law for x (loss of the upper hand of SMGRF texture model -the conditional Gaussianity for • x p ). However,

•

x p remain independent, thus can be sampled in parallel, but by more sophisticated samplers.

The new conditional laws are:

f (x p | * ) ∝ 1 + γ k λ k p (θ k ) β s | • x p | 2 -αs-1 exp -γ n | • y p - • h p • x p | 2
• integrate both x and s k and sample only γ n , γ k and θ k , which will be distributed under complicated laws.

Moreover, by variable marginalization, the resulting law is more diluted, thus what is gained by eliminating sampling steps, may be lost in terms speed of convergence. We have chosen not to integrate any of them since, although this implies more sampling steps, the sampled laws are easier to handle.

Since the law for θ has a non-standard, complicated form, it can not be sampled directly and a Metropolis-within-Gibbs strategy, similar to the one in Chapter 4, is employed. However, in this case, the efficient FMH is used.

Hence, it is obvious that sampling the a posteriori conditional laws implies an extra computational effort, both in terms of deployed algorithms, and from the law complexity point of view. Nevertheless, the advantages are significant, the samples no longer having insignificant likelihood. This considerably reduces the number of samples needed to reliably compute the evidences.

Implementation issues

The algorithm implementation has raised a series of numerical problems. Since the likelihood has an exponential form for each Fourier coefficient and consists in a multiplication over all the coefficients, this quantity often exceeds Matlab's representation capabilities and is set to infinity. To overcome this obstacle, the Negative Log-Likelihood (NLL) is computed instead of the likelihood for evaluating the acceptance probabilities in the MH algorithm. However, the likelihood is still needed for the evidence computation (both by the HMA and the LMA) and the indetermination problems caused by the variables set to infinity or zero [START_REF] Vacar | Bayesian Texture Model Selection using Harmonic Mean[END_REF] are again encountered. The employed solution is to determine the minimum value of each NLL chain, subtract it from all the NLL of the chain and compute the evidence as the harmonic mean or the LMA based on the "offset" values. The normalization is reversed in the final stage of posterior probability computation. The resulting CEAPS algorithm is given in Algorithm 2.

Experimental Results

The present section is devoted to describing and discussing the performances of our CEAPS model selection method for blurred and noisy textured images. A series of tests will be presented in the following:

Algorithm 2: Classifier based on Evidence Approximations from Posterior Samples (CEAPS) algorithm input : Data y, models dictionary for M = k, k = 1..K output: Evidences ẽk + samples for texture x (t) , noise parameter γ

(t) n , texture parameters γ (t) k , s (t) k , θ (t)
k (t = 1..T ) % prior model probabilities, e.g.,:

p k = 1/K; % generate samples of π(θ k , x, s, γ n |y, M = k) ; for k = 1 to K do % Gibbs sampler for (γ n , γ k , s k , θ k , x), M = k fixed: t = 1; initialization θ (t) k , x (t) = y, s (t) ; NLL k (t) = NegLogL(x (t) , s (t) , θ (t) k ); m(k) = min(NLL k ); nNLL k = NLL k -m(k); % compute the evidence using HMA or LMA ẽk = ComputeEvidence(nNLL k ); e rec (t) = 0; while |e rec (t) -ẽk | > ε do t = t + 1; γ (t) n ∼ f (γ n |y, x (t-1) , α n , β n ) γ (t) k ∼ f (γ k |x (t-1) , θ (t-1) k , s (t-1) k , α x , β x ) s (t) k ∼ f (s k |x (t-1) , θ (t-1) k , α s , β s ) θ (t) k -FMH with target f (θ k |x (t-1) , s (t) k ) x (t) ∼ f (x|y, γ (t) k , γ (t) n , s (t) k , θ (t) k
) NLL k (t) = NegLogL(x (t) , γ (t) k , s (t) k , θ (t)
k ); m(k) = min(NLL k ); nNLL k = NLL k -m(k); ẽk = ComputeEvidence(nNLL k ); e rec (t) = RecursiveEvidence(e rec (t -1), ẽk ); end end % determine the posterior model probabilities; for k ← 1 to K do Pr(M = k|y) = p k • ẽk K l=1 p l • ẽl • exp [m(k) -m(l)]
% Compute the parameter estimates by PM; θk = PM(θ k ); end 1. The first study compares the two MH samplers: the standard RWMH and the FMH.

2. The second set of tests compares the two evidence approximations based on posterior samples, i.e., the HMA and LMA.

3. The evaluation of our CEAPS classifier represents the third tests set. The classification performances for the CEAPS are first presented for various PSDs and then compared to those of the GMLE classifier, in a simplified scenario.

4. Results concerning the deconvolution are given in a visual form, by presenting the original textures, the observations and the deconvolved images.

RWMH vs FMH

The first set of tests investigates the speed performances of two sampling algorithms, the isotropic RWMH and our efficient version FMH, in the context of the complicated laws for θ.

Our tests indicate that the use of the FMH yields an algorithmic speed increase by a factor of at least 10 as compared to the RWMH. This is due to the FMH directional form of the proposal, which permits the algorithm to attain the high probability region in a very small number of iterations. Once in this region, the directional component will have negligible values and the algorithm will explore this high probability area of the parameter space due to the stochastic component of its proposal. This translates into a very short burn-in period, as opposed to the isotropic RWHM, which has a significantly longer burnin period, depending on the initialization. This efficiency is illustrated in Figure 5.2 where the NLL chains for the two samplers are represented.

HMA and LMA

The posterior samples obtained using the FMH within Gibbs sampler are used to compute the evidences. This is achieved via two approximations: the PM-based LMA and the HMA.

These approximations are computed using the same set of posterior samples in order to evaluate their accuracy in the same conditions. The numerical results show that the difference between the two evidence approximations is less than 0.1%, thus confirming that both approximations are viable for the problem in question.

Moreover, since the sampling itself is the most costly part of the evidence computation, the choice of approximation does not affect the overall speed performance. Consequently, the two computational methods for our CEAPS demand roughly the same amount of time.

In the tests presented in the following section, the CEAPS is based on the HMA.

CEAPS performances

Let us now present the performances of the selection method itself. The experimental setup consists in testing our method on synthetic textures, using 20 sets of parameter values for the PSD. Each set was used for each PSD model to generate both GRF and SMGRF texture realizations. The observations are obtained in a scenario with Gaussian blur of standard deviation w = 0.3 and SNR = 20dB. This corresponds to a partial overlap configuration of the texture PSD and the TF.

CEAPS

The algorithm was run on each texture realization and Tables 5.2 and 5.3 summarize the classification results for GRF and SMGRF textures, respectively. We observe on the main diagonal of both tables the percentages of correct classifications. As expected, the CEAPS chooses the correct PSD model in most cases. There are, however, situations where the method chooses another model. As anticipated by the Fisher information analysis, in the majority of cases, there is enough available information on the central frequencies to ensure their accurate estimation. Nevertheless, the information concerning the widths is more sensitive to the noise level and the PSD model and thus more prone to estimation errors. These errors are important for the method functioning, since they trigger misclassifications. The majority of cases where the CEAPS fails are due to high noise levels and consist in mistakenly considering a PSD with thicker tails as the most likely model. In this situation, the thicker tails account for the noise and the noise level is underestimated.

Nevertheless, in the context of our model choice problem, where the nested models help testing the method's ability to penalize model complexity, choosing another model can be regarded as not necessarily a failure. In this setting, the underlined percentages from Tables 5.2 and 5.3 represent the "good" miss-classifications, for instance, a Generalized Gaussian with q = 1 that is classified as an exponential. This illustrates the method's capacity to penalize model dimension, i.e., eliminate the parameters that do not significantly increase the model fit.

The method is not only able to distinguish between the different PSD forms, but also between the laws for the Fourier coefficients. More specifically, in 82% of the cases the algorithm correctly determined if the texture was from the GRF or the SMGRF class. This means that, on the one hand, the method has the ability to discriminate among a GRF and a SMGRF having the same form for the parametric part of the PSD. On the other hand, the PSD models themselves are structured enough to allow the algorithm to simultaneously identify the PSD model and whether all the PSD coefficients are identically scaled or not.

CEAPS vs GMLE

A crucial point is that the GMLE cannot solve the problem of interest. This is due to the presence of indirect data introduced by the: a) non-Gaussian texture model, b) blurred and noisy observations. Although the comparison cannot be performed on our problem, it is done on a simplified version of the problem consisting in direct observations (no noise and no convolution) of Gaussian textures. (b) Zoom on the NLL chains for the most likely models Figure 5.3: NLL chains to illustrate a typical situation where GMLE fails to select the good model each type of PSD shape. More specifically, this represents an averaging over the PSD models and PSD parameters. The lower classification performance of the GMLE is due to the fact that it does not have any mechanism of model complexity penalization, thus it chooses the most complex among the embedded models. On the contrary, as previously explained, the CEAPS penalizes model dimension and selects the less complex model that fits the data. Table 5.4 shows that the GMLE has a lower success rate and the reason for this is illustrated in Fig. 5.3. In this figure, we plotted a case where GMLE selects the Generalized Gaussian model, since its minimum NLL is the global NLL minimum among all models. However, the minimum neg-log PM is that of the Laplacian model, which is indeed the true model, this being a typical failure situation for GMLE.

As already stated in Section 5.1, the evidence based classifier is optimal from the risk point of view, which can be seen in this table through the CEAPS performances.

Visual reconstructions

The experiments show that situations with high noise, SNR < 20dB, are challenging since the samples for the widths have a too strong variance. Furthermore, these samples, used to compute the evidences, have a direct impact on the model selection process.

In practice, above a certain level of noise, the method tends to favor the PSD shapes with thicker tails, by considering that these thicker tails account for the noise. In fact, for high noise levels there is a smaller amount of information, thus more uncertainty in the estimation, which eventually triggers estimation errors for the PSD widths and even missclassifications. More specifically, the noise level is underestimated and either the PSD widths are overestimated, or a model with thicker tails is selected.

Using the samples employed to compute the evidences, we can also compute PM estimates for the texture parameters, the noise precision and the unknown image. Consequently, as an additional result, our algorithm provides a PM estimate of the original image, conditionally on the selected model. Fig. 5.4 shows examples of the reconstruction. We can observe situations (Figs. 5.4c, 5.4f and 5.4i) where CEAPS successfully restores the texture even if the observations are severely degraded. This illustrates the method's high capacity to handle the blur and the noise. This is due to the strength of the information given by the structure of the PSD and to the method's optimality from the classification and estimation risk point of view. Nevertheless, there are also situations, such as Fig. 5.4l, where the image is degraded to an extent that impairs a reconstruction, in most cases, this being due to a low information scenario.

Conclusion and Perspectives

In this chapter we have presented a method for texture model choice from indirect observations. The textured images are modeled by Scale Mixture of Gaussian Random Fields or Gaussian Random Fields with parametric Power Spectral Density. By applying a Bayesian formalism, we were able to determine the posterior model probabilities based on the evidences, this approach being optimal from the risk point of view. Moreover, the within-model simulation technique that is employed translates into a sweep of all possible models and the computation of the evidence for each model in the dictionary. This quantity can be determined only by numeric methods, since the integral in (5.2) is intractable. We have compared several methods for numerically computing the evidence based on samples from the a posteriori law and we have presented the performances of these methods, deciding on which is the most adapted method for our application.

As a secondary result, this approach provides chains of samples for the parameters, conditionally on each model M = k. These samples can be used to obtain estimates that are optimal from the mean square error point of view, by using the Posterior Mean estimator, the estimator that minimizes this error among all possible estimators.

The immediate perspective of this work is to compare our approach to other concurrent model choice methods.

Further developments include, but are not limited to, including the texture modeling related perspectives presented in Chapter 3, i.e., extending the texture model to dependent Fourier coefficients, and the class of compared PSD models. However, that will come at a cost, i.e., the coefficients dependence will be reflected in the computational complexity. Another idea is the use of multi-modal PSDs, in order to obtain more structured, quasi-periodic textures.

From a different perspective, related to the work presented in Chapter 4 and [START_REF]Bayesian Texture and Instrument Parameter Estimation From Blurred and Noisy Images Using MCMC[END_REF], the method can be adapted to deal with unknown non-parametric PSFs (blind) or parametric PSFs (semi-blind) and the estimation of their parameters. The most complex problem in this thesis is the one of textured image segmentation from indirect observations. Most existing approaches do not tackle the indirect observations issue and focus only on segmentation. Image segmentation is a computer vision problem consisting in partitioning an image into several groups of adjacent pixels that have a certain homogeneity property (gray level, color, texture, or other features) or that compose an object of interest.

State of the art

The literature in the field of image segmentation is extensive as this topic has been of great interest for decades. Classifying the existing approaches is not an easy task, however, from the beginning, a distinction can be made between the region-based and the contourbased methods. An important aspect related to the region-based approaches is that they always provide closed contours, as opposed to the contour-based methods, for which the contours are not necessarily closed.

The most straightforward image segmentation method is thresholding, however, this method is seldom applicable, since it is only adapted for piecewise constant images. In this context, a large variety of approaches are defined for segmentation.

Region growing methods are iterative algorithms that start from initial seeds and expand each region in its neighborhood based on a similarity criterion. [START_REF] Zhang | A diffusion approach to seeded image segmentation[END_REF] presents a seeded image segmentation based on a heat diffusion process, [GUSV + 09] describes an unsupervised region growing and multiresolution merging algorithm for color image segmentation, based on color edges and [START_REF] Alpert | Image Segmentation by Probabilistic Bottom-Up Aggregation and Cue Integration[END_REF] presents a bottom-up aggregation image segmentation.

Partial differential equations based techniques have also been employed for image segmentation, for instance in [START_REF] Chan | An Active Contour Model without Edges[END_REF] introducing an active contour without edges method for object segmentation, based on level sets, curve evolution and the Mumford Shah model.

A very popular approach to image segmentation is graph partitioning, where each pixel of the image is considered as the node of a graph and the goal is to connect the nodes having certain similarity properties. [START_REF] Felzenszwalb | Efficient Graph-Based Image Segmentation[END_REF] uses a graph based image model and measures the evidence for a boundary between two regions in order to achieve an image segmentation. [START_REF] Boykov | Graph Cuts and Efficient N-D Image Segmentation[END_REF] describes the basic framework for efficient object extraction from multi-dimensional image data using graph cuts. Also based on graph cuts, [START_REF] Barbu | Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities[END_REF] presents an image segmentation method based on a generalized Swendsen-Wang sampler. This method uses a Potts prior for the labels and, based on an adjacency graph, computes probabilities for each edge, performs graph clustering and graph flipping (instead of single vertex flipping as in the case of the Gibbs sampler).

The watershed approach is based on the topography of the gradient magnitude on the image. [START_REF] Malik | Contour and Texture Analysis for Image Segmentation[END_REF] presents a normalized cuts approach relying on a local measure of similarity of the textural characteristics in a neighborhood of the pixel, while [Gra06] uses a small number of predefined labels and computes for each unlabeled pixel a probability based on the speed of a random walk to reach a pre-labeled pixel. The final label assigned to a pixel is the one maximizing this probability. [START_REF] Sinop | A Seeded Image Segmentation Framework Unifying Graph Cuts And Random Walker Which Yields A New Algorithm[END_REF] unifies the graph cuts and the random walker methods in a common framework, based on L q norms minimization for seeded image segmentation.

One of the first approaches for textured image segmentation is presented in [START_REF] Tuceryan | Moment-based texture segmentation[END_REF], which is based on using the moments in small windows of the image as texture features. The topic of image segmentation for textured images has been previously addressed in [START_REF] Arivazhagan | Texture segmentation using wavelet transform[END_REF]. This method consists in computing features based on the Discrete Wavelet Transform coefficients of blocks of the image, evaluating the difference between these features on adjacent blocks and applying post-processing techniques in order to obtain a one-pixel thick contours. This method does not provide a label field, thus no information about which texture belongs to which region. Another method providing texture edges [START_REF] Wolf | Patch-based texture edges and segmentation[END_REF] uses active contours and the patch based approach for texture analysis. However, none of the aforementioned approaches is formulated in the context of indirect observations. Image segmentation for textured images is achieved in [START_REF] Lillo | Supervised Segmentation Based on Texture Signatures Extracted in the Frequency Domain[END_REF] based on features extracted from the Fourier transform of the learning textures. A significant method for image segmentation based on both gray level (intervening contour framework) and texture (textons) is presented in [START_REF] Malik | Contour and Texture Analysis for Image Segmentation[END_REF]. The segmented image is obtained using a normalized cuts approach. [MRY + 11] presents a method that models a homogeneous textured region of an image by a Gaussian distribution and the region boundaries by adaptive chain codes.

The segmentation is obtained using a clustering process. An image segmentation approach devoted to strongly resembling textures in given in [START_REF] Galun | Texture segmentation by multiscale aggregation of filter responses and shape elements[END_REF]. The goal is to accurately characterize the textures and this is achieved by combining a collection of statistics and filter responses. This local information is then used in an aggregation process to determine the segmentation. A three stage segmentation method is presented in [START_REF] Liu | Image and Texture Segmentation Using Local Spectral Histograms[END_REF] and relies on characterizing both textured and non-textured regions using local spectral histograms. Texel based image segmentation is achieved in [TA09] by identifying the modes in the pdf of region properties. Nevertheless, this method does not deal with the aspect of indirect observations either.

A very significant class of segmentation methods are those relying on a probabilistic formulation and using models for the image and the labels. [START_REF] Geman | Boundary Detection by Constrained Optimization[END_REF] presents an approach for image partitioning into homogeneous regions and for locating edges based on disparity measures. In [START_REF] Tu | Image segmentation by data driven Markov chain Monte Carlo[END_REF], an image segmentation method is developed based on MCMC and the K adventurers algorithm with the goal of preserving ambiguities in the segmentation process. This method integrates clustering and edge detection in the proposal probabilities. [START_REF] Deng | Unsupervised image segmentation using a simple MRF model with a new implementation scheme[END_REF] introduced a weighted MRF model that estimates the model parameters and thus performs unsupervised image segmentation. [START_REF] Mignotte | A Segmentation-Based Regularization Term for Image Deconvolution[END_REF] formulates a Bayesian approach to image deconvolution with a regularization term based on a segmentation map, estimated using a Potts prior for the labels and a Gaussian model for the image, conditionally on the labels.

One of the most commonly used model for the labels in the probabilistic approaches is the Potts model. [START_REF] Celeux | EM-based image segmentation using Potts models with external field[END_REF] explores the Potts model based image segmentation by introducing a site dependent external field in the Potts model, while [START_REF] Pereyra | Estimating the Granularity Coefficient of a Potts-Markov Random Field Within a Markov Chain Monte Carlo Algorithm[END_REF] proposes a method for jointly estimating the temperature parameter β of the Potts label prior using a likelihood free MH algorithm.

A very interesting work is the Bayesian method for image segmentation from indirect data (with known blur and Gaussian noise) in [START_REF] Ayasso | Joint NDT Image Restoration and Segmentation Using Gauss-Markov Potts Prior Models and Variational Bayesian Computation[END_REF]. This method is based on a Potts model for the labels, the pixels that belong to different regions being considered independent of each other. Within a region, the pixels are either independently Gaussian, or Markov-Gaussian. This type of image model makes the method suited mostly for the piecewise constant images. Nevertheless, the problem formulation in the Markov-Gaussian case is rather complicated and slightly unclear, resulting a difficult comprehension of the topic. This chapter addresses a similar problem to the one in [START_REF] Ayasso | Joint NDT Image Restoration and Segmentation Using Gauss-Markov Potts Prior Models and Variational Bayesian Computation[END_REF] and presents a probabilistic method for image deconvolution and segmentation, based on a Potts model for the labels. However, in our case, the pixels of each region are modeled by a GRF. This allows us to segment images containing different textures, with the same gray level, but different textural characteristics.

Problem Statement

The current chapter presents a segmentation method for textured images affected by blur and noise. To be more exact, y is the blurred and noisy observation of the original image x. The unobserved image x is composed of several regions, each of these regions consisting in a different texture. In our case, all the textures have the same mean gray level. This means that the information regarding the gray level in a neighborhood of a pixel is not pertinent in distinguishing the regions. In this context, the goal is to jointly segment and restore the image based on the textural characteristics.

To summarize, in this problem:

• the original image x -contains several regions, -each region consists in a patch of texture belonging to one of K classes, -each of the K textures is modeled by a GRF with parametric PSD. The parametric PSD has a known model driven by the unknown parameter set θ;

• the PSF is known;

• the hyperparameters (γ n -noise parameter, and γ k -scale parameters for the textures) are unknown.

Let us consider that each patch of texture is extracted from a full image x k following the GRF model from Chapter 3, driven by the parameters θ k and γ k . Then, the variable hierarchy of our segmentation problem from indirect data is graphically represented in Figure 6.1:

y γ n α n , β n x 1 θ 1 γ 1 θ m 1 , θ M 1 α 1 x , β 1 x x K θ K γ K θ m K , θ M K α K x , β K x z β Figure 6
.1: Hierarchical variable dependency for the deconvolution-segmentation of textured images.

Based on this variable dependency, the joint law for this problem can be expressed as:

f (y, z, x 1..K , γ n , γ 1..K ,θ 1..K ) = f (y|γ n , z, x 1..K ) • π(γ n ) • f (z|β) • K k=1 f (x k |θ k , γ k ) • K k=1 π(θ k ) • K k=1 π(γ k ) (6.1)
As can be noticed in the variable dependency and in the joint law, the image x is not directly probabilized and is not assigned a prior. Nevertheless, this variable represents a deterministic transformation of a series of probabilistic quantities. The process of obtaining the image x containing the K textures, starting from the full textures x k , k = 1..K and the labels z, can be mathematically formalized as follows:

x = S 1 x 1 + S 2 x 2 + ... + S K x K = K k=1 S k x k (6.2)
where S k = S k (z) are P × P size diagonal binary matrices obtained based on the labels z. These matrices extract from the single texture image x k the pixels with label k and replace the other pixels with 0. A more detailed description of this extraction and zero-padding process is given in Appendix D. The structure of the aforementioned matrices is:

S k = Diag {δ(z p , k), p = 1..P }
where the entries equal to 1 are on the positions p with z p = k.

This process can be illustrated by the schematic:

X 1 : S 1 X 1 : X 2 :
S 2 X 2 :

X 3 :

S 3 X 3 : ⇒ X :

Bayesian Formulation

By using the form in (6.2) for the original image in (2.3), we obtain the following likelihood:

f (y|z, x 1...K , γ n ) = (2π) -P • γ P n • exp -γ n y -H k S k x k 2 (6.3)
where H represents the circulant, known convolution matrix. We have chosen a parametric form for this filter, driven by the width parameter w f :

• h nm = exp - 1 2 w 2 f ν 2 n + ν 2 m (6.4)
This is an isotropic Gaussian filter, centered in the 0 frequency, and w f is the inverse width of the Fourier transform of the convolution filter, as already employed in Chapters 4 and 5.

A prior must be assigned to the hidden labels, which encodes all the information we possess about this hidden field. An appropriate choice is the Potts model, which captures the tendency of neighbor labels of having the same value. The strength of this tendency is parametrized by the inverse temperature parameter β:

f (z|β) = C z (β) • exp β r∼s δ(z r , z s ) (6.5)
Further details about this model are given in Appendix C.

As for the rest of the unknowns, for the textures x k we have employed the GRF model, previously described in Section 3.2.1 and used in Chapter 4. The γ n and γ k have been assigned Gamma priors, while for the shape parameters of the textures, θ k , we have chosen uniform priors.

Based on the likelihood law and the aforementioned priors, the joint law describing this problem (6.1) becomes:

f (y, z, x 1...K ,γ n , γ 1...K , θ 1...K ) = C y exp -γ n y -H k S k x k 2 • C γn γ αn+P -1 n exp (-γ n β n ) • C z (β) • exp β r∼s δ(z r , z s ) • C x k R k (θ k ) -1 exp -γ k x k 2 R k (θ k ) • k C γ k • γ α k +P -1 k • exp (-γ k β k ) • k U θ k m ,θ k M (θ k ) (6.6)
The a posteriori law, which is proportional to the joint law in (6.6), summarizes all the information about the unknowns contained by the data and the a priori models.

Estimators

Based on the aforementioned a posteriori law, parameter estimates can be obtained. The estimators we have chosen to use differ among the parameters:

• for the labels z we have chosen a marginalized MAP (mMAP) estimator;

• for γ n , γ k and θ k we have used the PM estimator;

• for the textures x k we have used the PM estimator, conditionally on the estimated labels ẑ.

An estimate of the original image x can be obtained based on the estimates for the labels ẑ and for the textures xk as follows:

x = Ŝ1 x1 + Ŝ2 x2 + ... + ŜK xK (6.7)
where the estimated extraction matrices Ŝk = Ŝk ( ẑ) are obtained based on ẑ.

Computing the Estimators -A Posteriori Conditionals

The complicated posterior law cannot be exploited directly to straightforwardly determine the estimates. Consequently, the information must be extracted using numerical methods. The Gibbs sampler is employed in this work. This algorithm sequentially samples the a posteriori conditional laws, which in this case write:

• the noise parameter γ n has a standard Gamma conditional, meaning that its sampling does not pose any difficulties;

f (γ n | * ) ∝ γ αn+P -1 n • exp -γ n y -H k S k x k 2 + β n (6.8)
• the PSD scale parameters γ k , k = 1..K, also have Gamma forms, which can be straightforwardly sampled;

f (γ k | * ) ∝ γ α k +P -1 k • exp -γ k x k 2 R k (θ k ) + β k (6.9)
• the conditional law of the PSD shape parameters θ k , k = 1..K, is highly nonlinear and has a non-standard form. Nevertheless, this law can be sampled using an RWMH step;

f (θ k | * ) ∝ p λ p (θ k ) • exp -γ k p λ p (θ k )| • x kp | 2 • U θ k m ,θ k M (θ k ) (6.10)
• the labels have a non-standard complicated law. The sampling of this law will be detailed in the next section.

f (z| * ) ∝ exp -γ n y -H k S k (z)x k 2 + β r∼s δ(z r , z s ) (6.11)
• the law for each complete texture x k , k = 1..K, has a Gaussian, non-separable form, of dimension P . The sampling of this law is explained in a subsequent section;

f (x k | * ) ∝ exp - 1 2 γ n y -H k S k x k 2 + γ k x k 2 R k (θ k ) (6.12)
Based on the values of the samples drawn from the posterior law via the Gibbs algorithm, the estimators are computed as follows:

-The mMAP labels estimation consists in selecting for each label the value that has been most often selected during the sampling process. In practice, this comes down to independently computing an a posteriori histogram of the sampled values for each label and selecting the value that corresponds to the maximum in the histogram. This value gives the estimate ẑ.

-The PM estimates of the rest of the variables consist in calculating the expectancy of the posterior law. Since this comes down to calculating an intractable integral, practically, the estimates are obtained numerically, by averaging the posterior samples.

Sampling Aspects

This section provides detailed explanations regarding the cumbersome task of sampling the labels, i.e., the f (z| * ) law in (6.11), and the full textures, i.e., the f (x k | * ) law in (6.12). These two sampling processes represent the major algorithmic challenges of our approach. Since the treated problem has not been tackled in the literature before, these aspects had never been encountered and addressed, either. For this reason, the sampling processes, especially that of the textures, have been conceptualized in several different forms that have allowed us to gradually mature our perspective and arrive at the version presented in this section.

Sampling the Labels

The conditional law for the labels (6.11) has a complicated and non separable form and thus its sampling is not an easy task. A first solution is to apply Gibbs sampling and thus sample the labels sequentially, each conditionally on the rest of the labels. Nevertheless, this solution implies a prohibitive computational time, thus the goal is to parallelize this process.

The prior Potts law has been defined using a 4-neighborhood. This prior field can be sampled using a Gibbs algorithm in only two steps. In each of these steps, we are sampling labels that are conditionally independent of each other. In this context, the label lattice can be viewed as a checkerboard [START_REF] Winkler | Image Analysis, Random Fields and Markov Chain Monte Carlo Methods: A Mathematical Introduction[END_REF], where the "white" labels are independent of each other, conditionally on the "black" labels. The same reasoning applies to the "black" labels, which are mutually independent, conditionally on the "white" ones.

However, these considerations are only related to the prior. The a posteriori conditional law for the labels also contains the likelihood term. This term can be expressed in a separable form with respect to the labels, thus the label sampling can be performed in a parallel, computationally efficient manner.

Remark: If this were not true, the sampling should have been performed sequentially, one label at a time, conditionally on the rest of the labels and variables. This would have translated into a prohibitively costly sampling process. Sampling a label, say z p comes down to computing all the probabilities:

Pr(z p = 1| * ) ∝ exp -γ n y -Hx p 1 2 + β r∼p δ(z r , 1) . . . Pr(z p = k| * ) ∝ exp -γ n y -Hx p k 2 + β r∼p δ(z r , k) . . . Pr(z p = K| * ) ∝ exp -γ n y -Hx p K 2 + β r∼p δ(z r , K) (6.13)
up to the multiplicative factor. This factor can be determined knowing that the sum of the probabilities is 1.

The x p k in (6.13) represents an image that has all the pixels identical to x, except for pixel p. The pixel at position p in x p k is obtained by extracting the pixel p from x k . In this manner, the data adequacy term measures which class is the most appropriate for describing pixel p in the sense that it provides the maximum level of coherence with the rest of the image.

The second term is the contribution of the prior. This term is meant to ensure a coherence of the current label with its neighbors. This is weighted by the inverse temperature parameter, β, which tunes the strength of the connection with the neighbor labels.

To compute the probabilities in (6.13), we must evaluate the value of the two terms at pixel p. The second term is easily computed by counting the neighbors of pixel p having label k and multiplying with β.

Let us now focus on the first term, that we will denote V p,k = y -Hx p k 2 . In order to write this term in a more convenient form, let us introduce the following auxiliary quantities:

-a vector 1 p in which only one element is different from 0. This element is equal to 1 and has the index p: 1 p = 0 0 . . . 0 1 0 . . . 0 t P -a scalar quantity δ kp ∈ C that records the difference between the pth pixels of the full texture image x k and of the image x:

δ kp = -x p + x kp
Then, using these notations, V p,k can be written as:

V p,k = y -H x -δ kp 1 p 2 = (y -Hx) ȳ -δ kp H1 p 2 = ȳ † ȳ + δ k 2 p 1 † p H † H 1 p H1p 2 -2δ kp 1 † p H † ȳ (6.14)
Let us analyze each of these terms:

1. The first term ȳ † ȳ does not depend on p or k. This means that its computation is not required in the process of label sampling, thus this term can be included in the multiplicative factor;

2. The H1 p 2 factor of the second term does not depend on p. This is due to the circulant form of the H matrix. For this reason, this norm only needs to be computed once for the entire problem. In fact, this norm amounts to the sum P p=1 |

• h p | 2 . Consequently, we only need to compute the δ kp , which can be done in parallel for all values of p and k and stored in a N × N × K block.

3. Finally, the third term can be written as:

δ kp 1 † p H † ȳ = δ kp 1 † p H † ȳ FFT select (6.15)
The product H † ȳ can be computed efficiently by FFT, once for all the iterations. The product with 1 † p consists in selecting the pixel p of H † ȳ.

Sampling the Image

The law in (6.12) is Gaussian, however, its sampling will pose problems due to the high dimension of the variable. Let us rewrite this expression as:

f (x k | * ) ∝ exp - 1 2 x k -m k 2 Σ k (6.16) where Σ k = γ n S † k H † HS k + γ k R -1 k (θ k ) -1 = Q -1 k m k = Σ k • γ n S † k H † ȳk ȳk = y -H l =k S l x l (6.17)
This ȳk represents in fact the observations from which we are subtracting the contribution of all the regions, except region k. In fact, ȳk collects the part of the observations that are relevant about texture k.

Remark: To improve the readability of the equations, in the following we will use the simplified notation R k = R k (θ k ).

The P -dimensional Gaussian law in (6.16) has mean m k and covariance matrix Σ k . The sampling of this high dimension Gaussian requires computing the matrix Σ k , which is a very costly or even impossible operation due to the fact that it implies determining the inverse of a P × P matrix (for instance, for a 256 × 256 pixel image, Σ k is of size 65536 × 65536). This can be performed efficiently only in the case where the covariance matrix is sparse or has a special structure (circulant).

In our case, the R k , H and, by extension, the H † H matrices are circulant. However, the presence of the S k matrix breaks the circularity of the first term, making it impossible to perform the sampling straightforwardly and efficiently by FFT. Moreover, since this sampling operation must be performed at each iteration, it is obvious that a direct computation would render the algorithm prohibitively expensive.

Nevertheless, the literature accounts for a method of sampling high dimension Gaussians that have the covariance matrix with a certain structure [START_REF] Orieux | Sampling High-Dimensional Gaussian Distributions for General Linear Inverse Problems[END_REF]. This method is called Sampling by Perturbation-Optimization and consists in building a "perturbed" criterion, based on the target law, and optimizing this criterion in order to obtain a sample of the target. The Sampling by Perturbation-Optimization can be applied in cases where the precision matrix and the mean can be written as a sum of the form:

Π = C c=1 M t c R -1 c M c µ = Π -1 C c=1 M t c R -1 c µ c (6.18)
Remark: This sampling method is only sure to provide samples of the target law if the optimization is perfect. The work in [START_REF] Gilavert | Rééchantillonnage gaussien en grande dimension pour les problèmes inverses[END_REF] provides a means of ensuring that the obtained samples are indeed from the target even in the case of an imperfect optimization.

By identifying Π = Σ -1 k , the most obvious choice is C = 2 and:

     M † 1 = S † k H † R 1 = 1/γ n I P µ 1 = ȳk      M t 2 = I P R 2 = 1/γ k R k µ 2 = O P
The perturbation phase of this algorithm consists in drawing the following samples:

ξ 1 ∼ N (µ 1 , R 1 ) ξ 2 ∼ N (µ 2 , R 2 ) (6.19)
The cost of sampling these variables is not prohibitive: ξ 1 corresponds to a decorrelated white noise of mean µ 1 and ξ 2 corresponds to a texture having the same model as the prior for x k . ξ 2 can be easily obtained by FFT.

In order to obtain a sample of the textured image x k , the following criterion must be optimized:

J(x k |ξ 1 , ξ 2 ) =γ n (ξ 1 -HS k x k ) † (ξ 1 -HS k x k ) + γ k (ξ 2 -x k ) † R -1 k (ξ 2 -x k ) =x † k Q k x k -2x † k q k + cte. (6.20) where q k = γ n S † k H † ξ 1 + γ k R -1 k ξ 2 .

Optimization algorithm

The advantage of this approach is that it does not need to store or compute the inverse of the matrix Q k , which is of size P × P . We have focused on optimal step gradientbased algorithms. The techniques we are employing are all guaranteed to converge. The advantage of these approaches is that the only necessary ingredients for performing the optimization are the methods to compute the product Q k x k and q k .

• The product can be written as:

Q k x k = γ n S † k H † H S k x k select FFT select (6.21)
and thus efficiently computed through a series of FFT and pixel selection. The selection process is achieved in the spatial domain, while the costly matrix multiplications are performed in the Fourier domain.

• The q k term writes:

q k = γ n S † k H † ξ 1 FFT select +γ k R -1 k ξ 2 FFT (6.22)
and can also be efficiently computed by FFT.

Theoretically, there is no constraint on the optimization technique to be used at this point. We have chosen the:

-gradient descent, -conjugate gradient.

The high dimension of the problem seems to favor the conjugate gradient optimization, since it is faster than the classical gradient descent. However, we have experienced convergence difficulties, i.e., the optimization algorithm did not converge in P iterations, making the overall algorithm very slow. In practice, we have noticed that the step length at each iteration was extremely low. For this reason, irrespectively of the initialization, the evolution in the parameter space was very slow. Consequently, the differences between the current values and the updates were almost insignificant. This represents the cause of the extremely reduced convergence speed. The solution for this problem was the use of the preconditioner:

M = γ n H † H + γ k R -1 k (6.23)
This form for the preconditioner has not been randomly chosen. It is in fact an approximation of the inverse Hessian of the problem that we are trying to solve. This approximation has a Circulant-block-Circulant form and was used instead of the actual inverse Hessian, which has no special form, for computational efficiency reasons. This preconditioner was used for both of the aforementioned optimization algorithms and we have compared the performances. Thusly, we have compared:

-a preconditioned gradient descent, -a preconditioned conjugate gradient. In this context, the two methods have yielded similar results, requiring only a few iterations to converge.

Finally, we have chosen to focus on the preconditioned gradient descent method. This algorithm requires at each iteration the computation of the descent direction. This is represented by the product Mg (i) k , where g

(i) k = Q k x (i)
k -q k is the gradient at iteration i. Since the matrix M is circulant, this product can be easily computed in the Fourier domain.

The second ingredient that is necessary is the step length at each iteration α i :

α i = g (i) k † M † g (i) k g (i) k † M † Q k Mg (i) k (6.24)
Figure 6.3 illustrates a set of full textures samples. A first remark that can be made is that these samples are not homogeneous. This is due to the data term for the texture x k , which contains information only on the positions having label k. For the rest of the pixels, the only contribution comes from the regularization term.

Approximating Q k

The most direct possibility to increase the computational efficiency in this case is to approximate Q k by a Circulant-block-Circulant matrix. This allows us the perform the sampling by FFT and yields a considerable gain in terms of speed.

The exact solution of (6.20) is:

xk = γ n S † k H † HS k + γ k R -1 k -1 γ n S † k H † ξ 1 + γ k R -1 k ξ 2 (6.25)

First circulant approximation

This approximation consists in eliminating the S k matrix from the expression of Q k , in order to keep only Circulant-block-Circulant matrices:

Qk = γ n H † H + γ k R -1 k (6.26)
The solution of the approximate system, corresponding to (6.25), is:

ẋk = γ n H † H + γ k R -1 k -1 γ n S † k H † ξ 1 + γ k R -1 k ξ 2 (6.27)
In this case, due to the Circulant-block-Circulant structure of the covariance matrix, the Fourier coefficients of x k are independent and thus can be sampled in parallel. Each of these coefficients, say

•

x kp has the form:

• x kp = γ n | • h p | 2 + γ k λ p (θ k ) -1 γ n • h † p • ξ 1p + γ k λ p (θ k ) • ξ 2p
(6.28)

Second circulant approximation

A second and more rude approximation for the matrix Q k is to completely eliminate the H:

Qk = γ n I + γ k R -1 k (6.29)
The solution of the new approximate system is:

ẍk = γ n I + γ k R -1 k -1 γ n S † k H † ξ 1 + γ k R -1 k ξ 2 (6.30)
Similarly to the case of the first approximation, the covariance matrix has a Circulantblock-Circulant structure, i.e., it is diagonalizable by FFT. Consequently, the Fourier coefficients of x k are independent and can be sampled in parallel. Each of these coefficients, say

•

x kp , has the form:

• x kp = [γ n + γ k λ p (θ k )] -1 γ n • h † p • ξ 1p + γ k λ p (θ k ) • ξ 2p (6.31)
Remark: This approximation does not increase the speed performance as compared to the first approximation, being in fact a special case of the latter.

Cost of the approximations 6.5. Results

87

The proposed approximations allow for a significant speed increase as compared to the "exact" method based on sampling by optimization. Nevertheless, this comes at a cost, since we are in fact modifying the criterion that must be optimized. Consequently, it is interesting to analyze the difference between the expected results. This can be done by comparing the expectation and the variance of the solutions for the three systems:

• Exact system E [ xk ] = γ n S † k H † HS k + γ k R -1 k -1 γ n S † k H † ȳk var [ xk ] = γ n S † k H † HS k + γ k R -1 k -1
(6.32)

• First approximation

E [ ẋk ] = γ n H † H + γ k R -1 k -1 γ n S † k H † ȳk var [ ẋk ] = γ n H † H + γ k R -1 k -1 γ n S † k H † HS k + γ k R -1 k • γ n H † H + γ k R -1 k -1 • Second approximation E [ ẍk ] = γ n I + γ k R -1 k -1 γ n S † k H † ȳk var [ ẍk ] = γ n I + γ k R -1 k -1 γ n S † k H † HS k + γ k R -1 k • γ n I + γ k R -1 k -1
By analyzing the statistics of these solutions, we can notice that they differ both in mean and in variance. From a theoretical point of view, this means that we are not actually sampling the same law. The following section will provide a comparison of the method performances in the three cases in order to evaluate the impact of the approximations on the estimation results.

Results

This chapter has presented so far our method for joint image segmentation and deconvolution with hyperparameters and texture parameters estimation. This is a very complex problem, which has confronted us with a series of algorithmic challenges, especially related to the sampling of the labels and of the textures. These challenges were overcome due to our adapted formulation of the problem and the resulting method provides not only an estimated label field, but also an estimate of the original image.

The problem of segmentation has a considerable degree of difficulty, especially in this case, where the regions consist in textures with the same mean and the data is affected by blur and noise.

The practical implementation of our deconvolution-segmentation method is described by Algorithm 3. Its implementation has lead us to a series of practical considerations:

• we propose to initialize the labels by a white noise with values in the range 1..K. Our tests have shown that the algorithm converges faster when using this initialization as compared to an initialization by a constant label field;

• the use of the approximations for the Q k matrix in the texture sampling phase yields a speed acceleration by a factor of 5. Whilst the exact sampling method takes about 5 minutes for a 200 × 200 pixel image, the method based on approximated sampling by FFT takes less than 1 minute;

• the preconditioned gradient descent and the preconditioned conjugate gradient methods have similar performances in this case. The non preconditioned versions of these algorithms are very slow to converge for this problem;

• using a convenient initialization for the textures in the optimization algorithm can significantly reduce the convergence time. We have tested several initializations ranging from a constant image to a white noise. Obviously, the best performances are obtained for an approximation of the sample:

x init k = γ n H † H + γ k R -1 k -1 γ n S † k H † ξ 1 + γ k R -1 k ξ 2 (6.33)
This section is structured as follows:

1. A performance evaluation of the exact method is provided

• for different image topologies,

• in various combinations of blur and noise.

2. The influence of the inverse temperature parameter β on the estimation is analyzed. This influence is illustrated on a specific image topology that is rather sensitive to the choice of this parameter.

3. The estimation results of the exact method are compared to those of the two proposed approximations. This analysis is performed in order to evaluate the cost in terms of estimation quality that must be paid for the speed gain.

Evaluation of the exact method

Let us start our analysis with the evaluation of the exact method performances and a detailed interpretation of the results. Using the exact method means in fact that the samples of the complete textures x k are obtained through the sampling by Perturbation-Optimization method with a preconditioned gradient descent optimizer. Further on, we will also present segmentation results for different image topologies. This will allow us to evaluate the method versatility, adaptability and identify the potential weak spots. Figure 6.4 illustrates the method's performance on a typical image topology, in an observation scenario with w f = 4.5 and γ n = 60, corresponding to an SNR = 18dB.

This first example consists in an image composed of 4 regions and containing K = 3 different classes of texture. In this case, the deconvolution-segmentation algorithm con-verges to a label configuration very similar to the true labels, with only 1.2% of missclassified pixels, despite the considerable degradation of the image. Moreover, the texture parameters relative estimation error is small, less than 5%. The full textures x k are also accurately estimated, having the same characteristics as the original textures and also preserving the correct phase. The blur and the noise are eliminated and the resulting textured image strongly resembles the original image.

Label analysis

Our method is region-based, meaning that it provides closed contours, unlike a great part of the existing works in textured image segmentation. This can be seen in Figure 6.5 where we show the contour representing the set of boundary pixels between different regions, i.e., the set of pixels for which at least one of its 4 neighbors has a different label. One of the main advantages of using a probabilistic approach is the fact that it not only provides estimates for the unknowns, but also information regarding the uncertainty associated to these estimates. Figures 6.6 and 6.7 illustrate our analysis on the label estimates and their probability, in a 2D form and as 1D cross sections, respectively. Figure 6.6a shows the probabilities of the selected labels ẑ = z MAP . As specified in the previous sections, the estimate is represented by the label having the maximum marginal posterior probability. Nevertheless, for certain pixels this maximum probability can have a small value, while for others a high value. Obviously, the estimate is more reliable in the second case. We can observe this fluctuation in Figure 6.6a, this probability field being obtained as the maximum of the three probability fields shown in Figure 6.6c.

The uncertainty associated to the chosen label is high at certain locations in the image and it is safe to assume that at these locations there is a higher chance of selecting a wrong label.

Remark: The colored lines in Figures 6.6a and 6.6c indicate the position of the 1D cross sections presented in Figure 6.7.

The analysis of the selected label posterior probability can be done on any image, even when the true label configuration is unknown. In order to verify if indeed we are more prone to error in the area with small posterior probability, we have compared our estimation result to the true labels in a case where they are known. We can immediately notice in Figure 6.6b that all of the miss-labeled pixels are in fact positioned in these areas of weaker posterior probability, shown in Figure 6.6a. This reinforces our statement concerning the utility of the probabilistic approach, due to its ability to identify and quantify the problem difficulty.

The first plot in Figure 6.7 shows the superposition of the true label field and the estimated labels. We can notice that the estimated labels are very close to the true labels, the only errors being recorded at the region edges, where the contours may be off by a few pixels. Nevertheless, there are no oscillations and the regions are correctly identified.

The second plot shows the probability of the aforementioned selected labels. For the majority of the selected labels, the probability is quite elevated, approaching 1. The regions where this probability drops superpose perfectly on the contours. This means that the less reliable estimations occur around the region boundaries.

The last plot illustrates the probabilities of each label, individually. The probabilities of the selected labels are simply the maxima of these k probabilities.

Restored image analysis

Figure 6.8 is dedicated to the analysis of the image restoration results and the posterior statistics of the estimate. This figure superposes a cross section of the original image and the estimated image. The first aspect that should be noticed is that the scale is correctly estimated. Also, overall, the textural characteristics are very accurately restored, i.e., the 94 Chapter 6. Deconvolution Segmentation for Textured Images θ k parameters estimates are accurate.

The dotted lines represent the -3σ and +3σ levels around the PM image estimate, where σ is the posterior standard deviation, computed empirically from the samples. This defines a confidence region around our estimate and we can notice that the true image is always confined inside this interval.

Other image topologies

In the case of the second image topology, illustrated in Figure 6.9, although the number of textures is more reduced K = 2, the task is more difficult due to the shape of the regions. The presence of a relatively thin, continuous structure makes the tuning of the β parameter more difficult. For this reason, the label estimation presents flaws, especially under the form of discontinuities in the original structure. The third image, represented in Figure 6.10, permits the observation of another interesting aspect. Despite the good overall label estimation, there is a region that disappears completely. The particularity of this region is that although it has a significant size in terms of pixels, it has an elongated form and it is very thin. The reason behind this estimation error is again the β parameter and the representation capabilities of the Potts model itself, which seems to be more appropriate for representing compact, blobby regions rather than elongated structures. This statement is supported by the fact that the structure represented by the point of the "i" is correctly estimated, although it is smaller.

Blur and noise influence

Further results show how the estimation is affected by various observation conditions. Figure 6.11 illustrates the method performances in a weaker convolution case w f = 2 and higher level of noise γ n = 40.

The method performs very well in this case, the estimated label field being very close to the true labels. The superior quality of the estimates as compared to the results in Figure 6.9 can be explained by the fact that the observations are less degraded (an SNR of 21dB as opposed to 19dB) despite the higher noise level.

Influence of the β parameter

Figure 6.12 illustrates the influence of the β parameter in the case where we have encountered the aforementioned difficulty regarding the elongated structure. This test is performed in a case with no convolution. It is obvious that the choice of β significantly affects the estimation results. A small value permits to identify all the regions, even the thin ones, however, at the cost of producing parasitic regions. For a too large value, there are no parasitic regions, nevertheless, the thin structures are lost. In this work, the value of β is set by manual tuning. In cases where β must be tuned automatically, it can be included in the estimation process, however, the results will probably be less accurate than in the manual tuning case.

Influence of the approximations

Figure 6.13 shows the estimation results in a case with no blur. Thus, this is a particular situation of denoising-segmentation. The "exact method", consisting in sampling the x k textures using the sampling by Perturbabtion-Optimization algorithm is compared to the two approximations that we have proposed in Section 6.4.2.2 in order to enhance the speed performances. In this case, the two approximations for sampling the textures by FFT are equivalent, since there is no blur. This is in fact confirmed by the estimation results, which are very similar for the two approximations. Moreover, when compared to the exact sampling method, the results provided by the approximations based methods are practically indistinguishable. This means that the speed gain provided by the use of the approximations does not come at the cost of poor quality estimation for the labels and the image.

Conclusion and Perspectives

This chapter has presented our method for joint image deconvolution and segmentation, dedicated to textured images. This is a very difficult problem due to the large amount of unknowns and the complicated dependencies of the variables. The formulation of the problem itself has not been a trivial task and has demanded a long reflection time in order to find the best manner to accurately express the hierarchical dependencies. In this context, the most adapted choice was to model K full images, corresponding to each class, rather than model the pixels of the x image. This has allowed us to obtain an expression of the posterior in which the dependency with respect to each unknown has a relatively convenient form.

Our Bayesian method relies on sampling the posterior conditional laws. Based on these samples, estimators for the labels and the image are built. Nevertheless, the sampling process for the full textures and for the labels has also proved to be rather challenging and has demanded a certain amount of ingenuity to overcome the impasses. The proposed methodological and algorithmic original aspects have contributed to developing a method that is both theoretically sound and efficient for solving this difficult problem never tackled before.

The previous section has presented the results of a series of tests performed on various convolution and noise conditions, for different image topologies and for several values of the Potts field temperature parameter. These results have shown that the method is able to accurately segment the image, provide a good estimation for the textures parameters and thus restore the original image (affected by noise and blur).

Our tests have shown however that the method is sensitive to the tuning of the temperature parameter β. This leads us to the future developments of this work:

• the first perspective of our current work is to automatically estimate β;

• a second future contribution is the use of our SMGRF model for the constituent textures x k . Nevertheless, this will add an extra layer of sampling for the auxiliary variables. Except for the extra computational time, this extension should be straightforward and should not pose any methodological difficulty;

• the third future development of this work aims at performing a myopic deconvolution, i.e., considering that w f , the width of the convolution filter, is unknown and estimating it along with the rest of the parameters.

As it can be seen from this brief listing of the perspectives, the work on this topic is far from being over. Nevertheless, even in its current form, the method presented in CHAPTER 7

Conclusion and Perspectives

This work has presented a wide range of inverse problems of central importance in image processing, translated in the case of textured images. In this context, a series of advances have been made on an algorithmic, modeling and methodological level.

The algorithmic contribution consists in developing an efficient version of the Metropolis-Hastings algorithm. The proposition law of this new sampler is based on a component resembling the Newton direction, in which the inverse Hessian is replaced by the Fisher matrix. As a result, the problems related to taking the inverse of the Hessian matrix are avoided and the directional component of the proposal is always sure to have the direction of gradient descent. Moreover, in the case of our particular a posteriori conditional law, this proposal can be expressed only based on the gradient of the law.

As for the perspectives of this work, this sampler can be integrated in a wide spectrum of applications in order to accelerate the sampling process. ********** The model-related contribution consists in developing a model for texture analysis and synthesis based on a Scale Mixture of Gaussian Random Fields with parametric Power Spectral Density. The parametric model for the Power Spectral Density and the values of its parameters encode the textural characteristics and allow us to classify the texture or to generate new texture samples having the same features as the original one. This non-Gaussian model is built using a set of auxiliary variables such that the law of the image, conditionally on the auxiliary variables, is Gaussian and, marginally with respect to these variables, it is no longer Gaussian. In the current work, for algorithm efficiency purposes, these auxiliary variables are considered independent.

The perspectives regarding this direction of research are to use correlated priors for the auxiliary variables. This means that the auxiliary variables will no longer be independent. This would probably trigger an increase of the representation capabilities of our texture model. Nevertheless, the side effect is an increased computational load for sampling a texture. Another possible development is the use of more complex shapes for the Power Spectral Density and investigating how the phase of the Fourier Transform coefficients can be used to obtain more complex models. ********** The methodological contribution is threefold, since we have addressed separately three inverse problems in the context of indirect observations of textured images.

Chapter 7. Conclusion and Perspectives

Firstly, we have devised a method for myopic deconvolution and parameter estimation for textured images corrupted by a blur and noise. These textured images are modeled by a Gaussian Random Field with parametric Power Spectral Density. The Point Spread Function also has a parametric exponential form. The aforementioned parameters are unknown and are estimated by the method. This work can be immediately extended by using more complicated form for the convolution filter. Another extension can be the use of the Scale Mixture of Gaussians model for the textured image.

***

The second methodological contribution is an evidence-based classification method for indirectly observed textured images. As in the previous case, the textured images are affected by a blur and by noise and the goal is to determine the model and the parameters of the texture's Power Spectral Density. In this case, the convolution filter is known, thus we are not in a myopic deconvolution case, however, the noise level is unknown and must be estimated. The textures are modeled either by Scale Mixture of Gaussian Random Fields or by Gaussian Random Fields and the method selects the most probable model to have generated the realization. The method relies on sampling and the algorithmic implementation embeds the efficient Fisher matrix based sampler. This work can be continued by extending the dictionary of shapes for the Power Spectral Density. Another future development is the estimation of the convolution filter, along with the rest of the unknowns.

***

The final methodological advance is our deconvolution-segmentation method for textured images. As in the classification case, the convolution filter is known and the noise level is estimated. This method achieves very satisfactory results in the case of a very complicated problem. Texture segmentation in itself is a very challenging task, whilst the literature holds no accounts of joint deconvolution-segmentation methods for textured images. Our method is based on a Potts model for the labels, with a temperature parameter that is manually tuned, and on a Gaussian Random Field model for the textures, with a given set of models for the Power Spectral Density, but unknown parameters.

Besides the significant theoretical challenges related to this problem, the implementation of this method has confronted us with the difficult task of sampling the labels and the pixels in an efficient manner, to avoid having a prohibitive computational load. This constraint has gradually guided us towards the final modeling of this problem that is presented in this manuscript. Nevertheless, this is the result of a series alternatives that we have explored. Among all these formulations, the current version is the only one that allowed us to accurately represent the variable dependencies and to obtain good speed performances. This topic can be further developed by integrating an automatic estimation of the temperature parameter, in order to permit the method to adapt to any image topology without user intervention. Another extension can be to also estimate the convolution filter. *** Finally, an ambitious project would be to combine all these problems in order to obtain a method that, starting from blurred and noisy observations of an image made up of several regions, each with its own texture, would be able to estimate the convolution filter, segment the image and determine the model of the Power Spectral Density and the values of its parameters for each of the constituent textures. ********** All the methods we have devised are based on optimal estimators, such as the Posterior Mean for the parameter estimation and the evidence-based classifier for the model selection. Consequently, the methods themselves are optimal from the Mean Square Error and the Mean Classification Risk point of view.

The goal of this work was to provide answers to questions that have not yet been posed or answered in the literature in the context of indirect observations of textured images.

Truncation Matrices

The extraction matrices S k from Chapter 6 perform at the same time two tasks:

• they extract the pixels having label k,

• they perform zero-padding, i.e., they replace the pixels having labels l = k with zero.

This can be formalized by passing through the intermediary matrices T k , for which:

S k = T t k T k (D.1)
The T k matrices are binary truncation matrices, based on the labels z. They are used to extract the pixels with label k from the full image x k . Let us consider that:

. . .

0 0 0 0 • • • 1      P          card I k
The special properties of these matrices are that they cover the entire pixel lattice and that they select each pixel only once:

T t 1 T 1 + T t 2 T 2 + ... + T t K T K = I P (D.2) T k T t k = I card I k (D.3)
The second step of the process is the zero-padding, i.e., arranging the selected pixels on their original positions in the image and replacing of the pixels corresponding to labels l = k by zeros. This can be achieved using the T k matrix as follows:

x k = T t k x I k = T t k T k x k (D.4)

  Influence of the approximations . . . . . . . . . . . . . . . . . . 6.6 Conclusion and Perspectives . . . . . . . . . . . . . . . . . . . . . . . . 7 Conclusion and Perspectives A Fisher information for indirectly observed GRF textures B Optimal Bayesian Estimation B.1 Posterior Mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Evidence based Classifier . . . . . . . . . . . . . . . . . . . . . . . . . .

Me2.

  le choix du modèle de la Densité Spectrale de Puissance d'une image texturée, à partir d'observations indirectes, Me3. la déconvolution-segmentation d'images texturées, avec estimation des hyperparamètres et paramètres de la texture, tandis que la dernière contribution est algorithmique : xii Contents A1. Développement d'une version efficace de l'algorithme d'échantillonnage de Metropolis-Hastings, basé sur la matrice d'information de Fisher, et son intégration dans nos algorithmes.

  Contents xvétiquettes et EAP pour le reste des paramètres.

  ********** Du point de vue méthodologique, notre contribution est triple, car nous traitons séparément trois problèmes inverses dans le contexte d'observations indirectes d'images tex-xvi Contents turées. D'abord, nous avons développé une méthode pour la déconvolution myope et l'estimation des paramètres pour des images texturées convoluées et bruitées. Ces images texturées sont modélisées par des champs aléatoires gaussiens avec des Densités Spectrales de Puissance paramétriques. La fonction d'étalement du point a une forme paramétrique aussi, appartenant à la famille exponentielle. Les paramètres de ces deux champs sont inconnus et sont estimés par notre méthode.

  projet très ambitieux serait de combiner tous ces problèmes pour obtenir une méthode qui, à partir d'observations convoluées et bruitées d'images texturées, arrive à estimer le filtre convolutif, segmenter l'image, sélectionner le modèle et estimer les paramètres des Densités Spectrales de Puissance de chaque texture présente dans l'image. ********** Toutes les méthodes qui ont été développées sont basées sur des estimateurs optimaux, comme la Moyenne A Posteriori, pour l'estimation des paramètres, et le classifieur basé sur le calcul de l'évidence, pour le choix du modèle. Par conséquent, les méthodes ellesmêmes sont optimales du point de vue de l'erreur quadratique moyenne et du risque de classification moyen. L'objectif principal de ce travail a été d'offrir des réponses à des questions qui n'avaient pas été encore posées ou résolues dans le contexte d'observations indirectes d'images texturées.

1. 1 .

 1 Figure 1.1: Texture realizations

Figure 2 . 2 :

 22 Figure2.2: The stabilization of the PM for several chains in the case of each of the four studied methods. Each chain corresponds to a different initialization. As expected, all the algorithms converge to the same value, but the FMH reaches equilibrium the fastest. Although superior in terms of computation speed per iteration, the IMH and RWMH require a longer interval to converge.

Figures 2 .

 2 Figures 2.2 and 2.3 illustrate the samples evolution in the case of the MH samplers

  Introduction and State of the art . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Texture Modeling by Random Fields . . . . . . . . . . . . . . . . . . . . 25 3.2.1 Gaussian Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.2 Non-Gaussian Model . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3 Conclusion and Perspectives . . . . . . . . . . . . . . . . . . . . . . . . 35

Figure 3

 3 Figure 3.1: Example of an isopole pair, presented in [Jul80].

Figure 3

 3 Figure 3.4: Texture realizations of a zero-mean GRF, with the PSD consisting in a single, continuous parametric shape.

Figure 3

 3 Figure 3.5: Texture realizations of a zero-mean GRF, with the PSD composed of different uniformly sampled continuous parametric shapes.

Figure 3

 3 Figure 3.6: Texture realizations for the GRF (left) and the GSMRF (right) models with the same PSD.

  Figure 3.7: Texture realizations for the GSMRF model with scale parameters having a Potts prior and a Gamma prior, respectively.

Figure 4 . 1 :

 41 Figure 4.1: Hierarchical variable dependencies for the textured image unsupervised myopic deconvolution.

  TF and PSD partial overlap

Figure 4 . 2 :

 42 Figure 4.2: Illustration of different relative positioning and widths for the TF and PSD, resulting in different situations information-wise.

Figure 4 . 3 :

 43 Figure 4.3: Fisher information for the noise parameter γ n in the case of different signal levels, filter widths, PSD widths and relative TF and PSD positioning.

Figure 4 . 4 :

 44 Figure 4.4: Fisher information for the PSD central frequency ν 0 x as a function of different parameters. The legends are the same as those from Figure 4.3.

Figure 4

 4 Figure 4.5: Fisher information for the PSD width u x as a function of the SNR, the filter width and the relative positioning of the TF and the PSD. The legends are the same as those from Figure 4.3.

Figure 4

 4 Figure 4.6: Illustrations of the algorithm performances for the cases presented in Figure 4.2. On the first row, the case with narrow, centered PSD and wide TF, on the second row, a PSD far from the null frequency and very narrow TF and on the third row, partial overlap.

Figure 4

 4 Figure 4.7: Example of observation system in the Fourier domain, with the corresponding observations and the deconvolution result, for various SNR values in a partial overlap case.

Contents 5. 1

 1 Model Choice -State of the art . . . . . . . . . . . . . . . . . . . . . . . 53 5.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.1 Texture coefficients law and PSD models . . . . . . . . . . . . . . 54 5.2.2 Probabilistic model choice . . . . . . . . . . . . . . . . . . . . . . 56 5.2.3 Joint law and priors for the model, image and parameters . . . . . . 57 5.3 Evidence calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.3.1 Evidence approximations based on posterior samples . . . . . . . . 58 5.3.2 Posterior sampling . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.3.3 Gibbs within-model posterior sampling . . . . . . . . . . . . . . . 61 5.3.4 Implementation issues . . . . . . . . . . . . . . . . . . . . . . . . 63 5.4 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.4.1 RWMH vs FMH . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.4.2 HMA and LMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4.3 CEAPS performances . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4.4 Visual reconstructions . . . . . . . . . . . . . . . . . . . . . . . . 68 5.5 Conclusion and Perspectives . . . . . . . . . . . . . . . . . . . . . . . . 70

Figure 5 . 1 :

 51 Figure 5.1: Hierarchical variable dependency for the texture model choice from indirect observations.

Figure 5

 5 Figure5.2: NLL evolution -posterior sampling with the two versions of the MH algorithm. On the abscissa it is represented the time in seconds. The two chains have the same initialization. The burn-in stage is considerably longer for the RWMH as compared to FMH.

Figure 5 . 4 :

 54 Figure 5.4: Reconstruction results: 4 cases to be read from left to right. On the first column the original, unobserved texture, x|θ * , M = k * , in the center the distorted observations, y, and on the right column the results of the deconvolution, for the selected model, x| θ, M = k.

Figure 6 . 2 :

 62 Figure 6.2: Image forming process based on the x k textures and the labels z.

Figure 6

 6 Figure 6.3: Samples of the full textures for a problem with K = 3.

  Figure 6.4: Segmentation results and reconstructed textured images. The temperature parameter is fixed to β = 1.

  Probability of selected labels (b) Miss-labeled pixels (c) Probability of each pixel of having label 1, 2 or 3, respectively.

Figure 6 . 6 :

 66 Figure 6.6: Link between the probability of the selected label and the estimation error. The probabilities of each possible label value.

Figure 6

 6 Figure6.7: Cross-section of the true and estimated label fields, the probability of the selected labels and the probabilities of every possible label value k.

Figure 6

 6 Figure 6.8: Cross-section of the true and estimated image. The dotted lines indicate the ±3σ interval around the PM estimate.

  Figure 6.9: Segmentation results and reconstructed textured images. The temperature parameter is fixed to β = 0.7.

  Figure 6.10: Segmentation results and reconstructed textured images. The temperature parameter is fixed to β = 0.9.

Figure 6 .

 6 Figure 6.11: Segmentation results in an observation case w f = 2 and γ n = 40 corresponding to a blurred SNR = 21.5dB.

  Figure 6.12: Segmentation results and reconstructed textured images for the two approximations in an observation case w f = 0 and γ n = 50 corresponding to an SNR = 24dB.

Figure 6

 6 Figure 6.13: Segmentation results and reconstructed textured images for the two approximations in an observation case w f = 0 and γ n = 10 corresponding to an SNR = 8dB.

I

  k = {i|z(i) = k} → T k : x I k card I k ×1Then the T k matrices are of size cardI k × P and have the structure:

  Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Addressed Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Efficient Metropolis-Hastings Samplers . . . . . . . . . . . . . . . . . . 13 2.4.1 Independent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

	Chapter 1. Introduction
	with the label field and the role of these matrices used in Chapter 6.	CHAPTER 2
	Inverse Problems, Bayesian Framework
	and Stochastic Sampling
	Contents
	2.1 Inverse Problems in Image Processing . . . . . . . . . . . . . . . . . . .	7
	2.1.1 Indirect Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	8
	2.1.2 Full Inversion Problem . . . . . . . . . . . . . . . . . . . . . . . . 10
	2.2 Bayesian Approach -Prior, Joint and Posterior Laws . . . . . . . . . . 11
	2.2.1 Uninformative Priors. Conjugacy . . . . . . . . . . . . . . . . . . 12
	2.2.2 2.4.2 Standard Random-Walk . . . . . . . . . . . . . . . . . . . . . . . 16
	2.4.3 Langevin adapted Random-Walk . . . . . . . . . . . . . . . . . . . 16
	2.4.4 Hessian adapted Random-Walk . . . . . . . . . . . . . . . . . . . 17
	2.4.5 Fisher adapted Random-Walk . . . . . . . . . . . . . . . . . . . . 18
	2.5
	• Appendix D -Truncation Matrices explains in detail the structure, the connection

  Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

	CHAPTER 4
	Unsupervised Myopic Deconvolution of
	a Textured Image
	Contents
	4.1

4.1.1 Deconvolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.2 Data, image and TF models . . . . . . . . . . . . . . . . . . . . . 41 4.2 Myopic deconvolution for textured images . . . . . . . . . . . . . . . . . 42 4.2.1 Information and qualitative estimation performance analysis . . . . 42 4.2.2 Bayesian setting: priors, posterior and conditional posteriors . . . . 46 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.4 Conclusion and Perspectives . . . . . . . . . . . . . . . . . . . . . . . . 51

  Table 4.1: Parameter estimation MSE. The first part of the table shows the performances at a fixed SNR value and for different relative positioning of the PSD and TF. The second part of the table shows the method's sensitivity to noise by quantifying the error in favorable situations information-wise, for different SNR levels.

	4.3. Results								49
	SNR	TF γ n	γ x	ν 0 x	ν 0 y	u x	u y	η
		Narrow 1.2 3.2	4	4.2 8.5	8	4.3
	20 dB	Wide 2.5 1.9	3	3.1 6.3 6.5	12
		Overlap 1.5 1.1 1.2 1.3 3.4 3.5	3
	15 dB Overlap 25 dB	3.2 10.8 12.5 11.8 15.9 17.3 14.5 1.9 1.6 1.1 0.9 3.1 3 2.7

Table 5

 5

	True model	Estimated model Lo GL Exp GG
	Lo	85 10	1	4
	GL	21 69	3	7
	La	2	4	87	7
	GG	4	8	16	72

.2: CEAPS model selection performance for GRF textures (correct classifications rate in %) for a partial overlap case.

Table 5 .

 5 4 lists the average classification success rate for the CEAPS and the GMLE, every method being tested on 20 texture realizations, with various parameter values for

				Algorithm Classification accuracy (%)
				CEAPS			89
				GMLE			86
	Table 5.4: Average model selection performance (classification success rate in %) com-
	parison between CEAPS and GMLE. The averaging is done over PSD model and PSD
	parameters.
	x 10 4						x 10 4
							-2.465				GaussGen
	-1			GaussGen Laplace		-2.47				Laplace Lorentz Student
	-1.5			Lorentz Student		-2.475
							-2.48
	-2
							-2.485
	9.95 -2.5	9.96	9.97	9.98	9.99	10	-2.49	9.96	9.97	9.98	9.99	10
					x 10 4					x 10 4
	(a) NLL chains for candidate models

  State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.3 Bayesian Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.3.1 Estimators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.3.2 Computing the Estimators -A Posteriori Conditionals . . . . . . . 79 6.4 Sampling Aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.4.1 Sampling the Labels . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.4.2 Sampling the Image . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.5.1 Evaluation of the exact method . . . . . . . . . . . . . . . . . . . . 89 6.5.2 Influence of the β parameter . . . . . . . . . . . . . . . . . . . . . 96 6.5.3 Influence of the approximations . . . . . . . . . . . . . . . . . . . 97 6.6 Conclusion and Perspectives . . . . . . . . . . . . . . . . . . . . . . . . 98

	CHAPTER 6
	Deconvolution Segmentation for
	Textured Images
	Contents
	6.1

Last, but not least, the algorithmic contribution is represented by the development of a new efficient version of the Metropolis Hastings algorithm, with a directional component of the proposal function based on the "Newton direction" and the Fisher information matrix. This particular directional component allows for an efficient exploration of the parameter space and, consequently, increases the convergence speed of the algorithm.To summarize, this work presents a series of methods to solve three image processing problems in the context of blurry and noisy textured images. Moreover, we present two connected contributions, one regarding the texture models and one meant to enhance the performances of the samplers employed for all of the three methods.

this chapter addresses a problem that had not been tackled so far, while achieving very satisfactory results.

Acknowledgements

Chapter 5. Model Choice for the Law and the PSD of a Textured Image estimator cannot be used in the case of incomplete data and thus it is not applicable for indirect observations or the non-Gaussian texture model (due to the presence of the auxiliary variables).

Among the most frequently employed classifiers is the well-known MLE-based An Information Criterion (AIC) [START_REF] Akaike | A new look at the statistical model identification[END_REF], with penalization term 2D k where D k is the dimension of model k. Equally popular is the Bayesian approach, with the Deviance Information Criterion [START_REF] Spiegelhalter | Bayesian measures of model complexity and fit[END_REF] and the Bayesian Information Criterion (BIC) [Sch78], also based on the maximized value of the likelihood and a penalization term of the form D k ln N , where N is number of observations. In the same class of Bayesian approaches are included the Bayes factors and the method used in this work, based on the evidence. This quantity is mostly known in the literature as the marginal likelihood, since it is in fact obtained by marginalization. However, for brevity, we will systematically be referring to it as the evidence (in favor of a certain model).

An important aspect is model complexity penalization. While the AIC tends to select too complex models, the BIC is prone to underfitting. Hence, from a methodological perspective, being able to quantify the performances of the model choice algorithm is crucial. For this reason, using a method that is optimal in some sense offers the certainty that, at least from that point of view, it cannot be outperformed by any concurrent approach. Moreover, all these criteria are based on approximations. The evidence based methods represent an alternative Bayesian approach [START_REF] Beal | Variational Algorithms for Approximate Bayesian Inference[END_REF] that does not rely on an approximation and can be formulated under the form of Bayes factors or the posterior model probabilities computation method presented in this work. This method is based on the formulation of an optimal decision function from the Bayes risk viewpoint. The form of the Bayes risk and details concerning the computation of this classifier are given in Appendix B. The optimality is achieved due to the fact that the cost is averaged over all possible data realizations, true models and parameter values. Moreover, for a binary cost function C(k, k * ) = 1 -δ(k, k * ), the method implicitly selects the model with the MAP probability.

Problem Statement

Our problem consists in texture model selection from an indirectly observed realization. More precisely, textures modeled by GRFs or SMGRFs, with PSDs taking various parametric forms, are observed via an imperfect system that introduces a blur and noise. Starting from these observations, the goal is to select the model that best describes the texture.

Texture coefficients law and PSD models

The texture model, say M = k * , is chosen among K possible models. The same notations used in the previous chapters are employed here for the unobserved image x, the observations y and the noise precision γ n . The texture parameters corresponding to model M = k are denoted by γ k , s k and θ k .

Chapter 6. Deconvolution Segmentation for Textured Images Algorithm 3: Deconvolution-Segmentation Algorithm for Textured Images input : Data y, number of textures K, temperature parameter for the Potts field β output: Samples for the labels z (t) , the textures x

1..K , β);

% count label occurrences for every pixel p LabelOcc(p, z p ) = LabelOcc(p, z p ) + 1; end % determine the parameter estimates; γn = PostMean(γ

Fisher information for indirectly observed GRF textures

The data may carry different amounts of information, this affecting directly the estimation performances. The amount of available information depends on the SNR = γ n /γ x or the particular parameter values. The mean available information on Ψ is quantified by the Fisher information matrix. We focus on each parameter ψ (any one of Ψ's elements), through its diagonal elements: the expectation of the log-likelihood's second derivative:

f (y|Ψ) can be obtained from (2.4) and the texture models:

where

After derivation, the expression contains first and second order derivatives of r p (Ψ) with respect to ψ. When taking the expectation, knowing that E y|ψ | • y p | 2 = r p (Ψ), the second order derivatives cancel out. Thus:

where r p is the derivative of r p w.r.t. the current parameter ψ. Then:

) is a decreasing function, hence, the smaller the γ n (the higher the noise level), the easier it is to estimate its level.

Similarly,

is also a decreasing function, i.e., the smaller the γ x (higher signal level), the easier its estimation.

Appendix A. Fisher information for indirectly observed GRF textures

For the texture parameters θ ∈ θ:

Another interesting case is the noiseless scenario (γ n = ∞):

The Fisher information regarding θ depends only on the PSD and its parameters (as in [START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF]). This means that the amount of information for the texture parameters does not depend on the form of the TF.

The same considerations hold for the TF parameter η:

thus, for low SNR, I (η) is small and for no noise, I (η) only depends on the TF:

Optimal Bayesian Estimation

A Bayes estimator is the optimal estimator with respect to a certain cost function. This estimator minimizes the expected value with respect to the parameter and the data of the cost function, also known as Bayesian risk.

The estimators employed in this work, the PM and the evidence based classifier, are optimal from the Mean Squared Error (MSE) and the Mean Classification Error (MCE) standpoint, respectively.

B.1 Posterior Mean

Let us consider the y be the observations and θ be the unknown parameter, with prior distribution π(θ).

Based on the observations, the estimator δ(y) = θ(y) can be built. A cost function C(θ, θ) will quantify the significance of being wrong. An optimal Bayesian estimator minimizes the Bayesian risk function, i.e., the posterior expected cost:

By averaging the cost function over all the possible observations and all the possible values of the unknown, we obtain the Bayes risk:

Then, the decision function that minimizes this risk is: The last equality holds true because the cost function is non-negative by definition, thus ρ(π(θ), θ(y)|y) ≥ 0 and minimizing the integral of a non-negative function is equivalent to minimizing the function at each point.

Let us now consider the cost function as being the squared error:

We can then determine the estimator that minimizes the MSE by computing:

which results in:

Consequently, the estimator that minimizes the MSE is the PM.

B.2 Evidence based Classifier

Let us now consider that the data y are realizations of the law f (y|M * ). The goal is to determine the value of M which is the most likely, from the Bayes risk point of view, to have generated the data. In this case, M is a discrete variable, with prior probability π(M ).

Similarly, a decision function ∆(y) = M (y) will select a certain model, based on the current observations. The cost function C(M * , M (y)) quantifies the cost of choosing the wrong model. The Bayes risk writes in this case:

The optimal decision function, i.e., the decision that minimizes the Bayes risk is:

Potts Model

The Potts model originates from statistical mechanics and was used to describe the particle interaction in a crystalline lattice. The model has later been introduced to image processing and more specifically to image segmentation, since it is appropriate for modeling the prior behavior of the labels.

In a Potts model with K states, each node of a graph is assigned an integer value between 1 and K, the assembly of all these assignments being called a configuration. The set of all possible configurations is {1...K} P , where P is the total number of pixels.

The energy of a configuration z is denoted as E(z), representing the number of edges in the graph with endpoints labeled differently. It has the form:

A probability distribution function, known as the Gibbs distribution is assigned to the set of states:

where C z (β) is the normalizing constant that makes g a probability distribution. Often, this quantity is also referred to as the partition function. The parameter β is a measure of the inverse temperature. The value of β influences the structure of the resulting Potts field as follows:

• for low β the Potts field will have small regions and an overall noise-like character, and is an intractable quantity, except for the special case K = 2, the Ising model [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF][START_REF] Giovannelli | Ising field parameter estimation from incomplete and noisy data[END_REF]. This intractability is due to the number of terms in the sum that increases exponentially with P .