Méthode proposée

Cependant, un réglage fin des hyper-paramètres de la fonction de coût (e.g., paramètre de régularisation) peut être fastidieux. Récemment, des approches dites supervisées basées sur l'apprentissage profond ont été proposées pour la restauration d'images. Elles reposent sur des modèles non-linéaires complexes dont les paramètres sont appris à partir de bases de données d'images. Ces méthodes sont flexibles et peuvent arriver à d'excellents résultats pourvu que la base de données soit de qualité et de dimension suffisante. De plus, malgré leur complexité, ces approches peuvent être implémentées efficacement sur des cartes GPU, grâce à des environnements de programmation adaptés (e.g., PyTorch). En plus du point mentionné ci-dessus lié à l'importance de la base de données, un autre inconvénient est le manque d'interprétabilité de ces méthodes. Enfin, aussi bien les techniques d'optimisation ou d'apprentissage, fournissent généralement une estimation ponctuelle de la solution. En particulier, elles permettent rarement d'extraire des statistiques d'ordre supérieur, et donc ne peuvent pas quantifier l'incertitude sur le résultat obtenu. Dans cette thèse, nous considérons une troisième famille d'approche, appelée l'approche bayésienne, pour traiter le problème de restauration d'image. Par rapport aux méthodes mentionnées précédemment, le cadre bayésien considère les paramètres recherchés (image/noyau dans notre cas) comme des grandeurs aléatoires plutôt que déterministes. L'objectif devient l'estimation de leur distribution a posteriori, sachant les données observées (i.e., l'image dégradée). L'application du cadre bayésien présente deux avantages principaux. Tout d'abord, la statistique du bruit, et des informations sur a priori sur les grandeurs inconnues peuvent être prises en compte explicitement, via la formule de Bayes. Deuxièmement, étant donné que toute la distribution a posteriori est estimée, il est facile d'obtenir n'importe quel moment, tels que la moyenne et la covariance, ou un intervalle de confiance pour une prise de décision ultérieure.

Néanmoins, une tâche complexe réside dans la mise en oeuvre de l'approche bayésienne, à savoir le calcul de la distribution a posteriori. En effet, dans la plupart des iv Contents formulations bayésiennes pour la restauration d'image, la constante de normalisation de cette distribution ne peut pas être calculée analytiquement même si le modèle a priori et le modèle d'observation sont simples (e.g., modèles gaussiens). De nombreuses méthodes algorithmiques ont été envisagées pour fournir une bonne approximation de la distribution a posteriori recherchée. Dans cette thèse, nous apportons des contributions à cette problématique, dans le cadre de deux problèmes issus de la restauration d'images : (i) l'identification de flou non stationnaire, (ii) la déconvolution aveugle. Pour le premier problème, nous proposons une approche de simulation stochastique, de type filtrage particulaire (Chapitre 3). Pour le second problème, nous proposons une méthode d'approximation par technique bayésienne variationnelle (Chapitre 4). Puis nous en proposons une version supervisée via une technique de déroulement de réseaux de neurones (Chapitre 5). Nos propositions méthodologiques sont étayées par un grand nombre de résultats numériques, impliquant des modèles de flou et bruit variés. Les contributions détaillées de chaque chapitre sont décrites ci-après.

Chapitre 3: Approche de filtrage particulaire pour l'identification de flou variant spatialement Les images produites par des instruments optiques souffrent souvent de flou causé par la diffraction de la lumière ou le mouvement de l'objet. La présence du flou provoque Contents v inconnus. Cela nous permet de traiter l'incertitude de manière systématique et l'inclusion de connaissances a priori sur les inconnues. En conséquence, nous sommes en mesure non seulement de fournir des mesures d'incertitude sur les paramètres à estimer, mais aussi de propager cette incertitude à des tâches où ces paramètres estimés sont utilisés (ici, la tâche de restauration). Dans notre approche, la densité de probabilité a priori fonction (pdf) des paramètres prend une structure markovienne, ce qui nous permet d'exploiter des approches d'inférence existantes pour l'inférence séquentielle. Les récursions bayésiennes obtenues ne donnent pas lieu à des formes explicites et nous recourons au filtrage particulaire ("particle filtering", PF) pour une résolution approchée du problème d'inférence. Nous commençons par la méthode classique bootstrap PF (BPF). Le BPF parvient à obtenir des résultats satisfaisants dans certains de nos exemples mais il présente certaines lacunes, par exemple dans des scénarios avec un faible bruit d'observation où se pose un problème de dégénérescence des particules. Nous proposons donc une nouvelle méthode appelée GIAnPF ("Generalized Interacting Annealed PF"). GIAnPF gère la dégénérescence des particules en considérant une séquence de distributions intermédiaires, qui sont des versions "recuites" de la distribution de filtrage initiale. Nous démontrons la bonne performance de l'algorithme proposé ainsi que sa robustesse à travers plusieurs expériences numériques, incluant trois formes différentes de flous paramétriques variables spatialement. Nous illustrons en Figure 1 les résultats d'identification obtenus par différentes méthodes à partir d'acquisitions simulées sur une image de microscopie décomposée en 64 patchs et avec un écart-type de bruit égal à 10 -2 . Sur cet exemple, les deux méthodes BPF et GIAnPF basée sur la méthodologie PF, fournissent de meilleurs résultats visuels, que les autres approches. De plus, en terme quantitatif, notre méthode GIAnPF atteint une plus faible erreur, au sens de la métrique RMSE. Nous avons également illustré la validité des estimations obtenues, et l'interêt de la quantification d'incertitude, pour résoudre des problèmes de restauration d'image. Le travail présenté dans ce chapitre a fait l'objet de deux publications [Huang et al. 2019, Huang et al. 2021a].

Chapitre 4: Approche bayésienne variationnelle pour la déconvolution aveugle Nous considérons le problème de déconvolution aveugle, où une image floutée et bruitée doit être restaurée sans connaissance du modèle de flou. Nous proposons dans ce chapitre un algorithme bayésien varationnel visant à fournir une approximation optimale (au sens de la divergence de Kullback-Leibler) de la distribution a posteriori des paramètres à estimer (ici, l'image, le flou, et des hyper-paramètres du modèle). Nous supposons la PSF stationnaire, et satisfaisant une contrainte d'égalité linéaire (e.g., somme à un, symmétrie axiale). Nous adoptons de plus un modèle a priori de type SAR sur la PSF, dont la moyenne est définie à l'avance, et une régularisation de type variation totale sur l'image. Nous cherchons ensuite à minimiser la divergence de Kullback-Leibler (KL) entre la distribution a posteriori obtenue et un modèle séparable (de type "champs moyens") de celle-ci, par un algorithme de minimisation alternée. Même dans le cas d'un bruit gaussien, la mise à jour n'a pas une forme explicite et nous proposons donc d'adopter une stratégie d'augmentation de données inspirée des méthodes semi-quadratiques en optimisation. Nous obtenons alors l'algorithme VBA (variational Bayesian algorithm), dont nous donnons une première forme adaptée à un bruit gaussien, puis une extension de celui-ci à un modèle de bruit plus général. L'algorithme obtenu est rapide et performant. Néanmoins, des tests préliminaires ont montrés une grande sensibilité des résultats au réglage d'hyper-paramètres tels que le poids de la régularisation sur le noyau, et la variance du bruit. Cela nous conduit à la contribution suivante, présentée dans le prochain chapitre.

Chapitre 5: Algorithme VBA supervisé par déroulement profond ("deep unrolling") Nous considérons ici le problème de déconvolution aveugle dans le cas d'un bruit gaussien. Nous proposons l'algorithme unfoldedVBA, qui applique le paradigme du déroulement profond, à l'approche VBA décrite dans le chapitre précédent. Les itérations VBA sont intégrées à une architecture de réseau neuronal profond, chaque itération devenant une couche du réseau. Cela permet (i) d'apprendre les hyperparamètres (en particulier, le niveau de bruit) de VBA de manière automatique et supervisée, (ii) d'améliorer la qualité des résultats en entraînant le modèle avec une fonction de coût liée à la métrique de qualité d'image, (iii) de mettre en oeuvre la méthode en tirant pleinement parti des ressources GPU disponibles, et donc de considérablement réduire le temps de restauration d'image pendant la phase de test. Contrairement aux méthodes standard d'apprentissage profond pour la déconvolution aveugle, notre méthode est l'interprétable, grâce à la technique de déroulement. Soulignons que les méthodes bayésiennes variationnelles ont déjà été employées dans un contexte d'apprentissage profond. Notamment, ce sont des méthodes standards pour l'entrainement des auto-encodeurs variationnels, et des réseaux de neurones bayésiens. Cependant, à notre connaissance, notre travail est le premier à étudier le Contents vii déroulement profond d'une technique bayésienne variationnelle. En complément des itérations VBA déroulées, nous intégrons à notre architecture des couches supplémentaires, permettant d'apprendre les hyper-paramètres et d'améliorer la qualité de l'image en sortie via un post-traitement. L'architecture complète est présentée en Figure 2, pour le cas d'images en niveaux de gris et couleur. Nous montrons la supériorité de notre approche par rapport à l'état de l'art sur plusieurs jeux de données, impliquant des flous variés (e.g., gaussien, mouvement). Nous montrons dans la Figure 3 un exemple de restauration pour une image couleur dégradée par un flou de focus. Nous voyons que l'approche proposée permet d'améliorer la qualité visuelle et quantitative de l'image restaurée, par rapport à l'algorithme initial VBA. Les contributions des chapitres 4 et 5 ont fait l'objet de la publication [Huang et al. 2021b]. x Contents spatialement étant donné une image originale et sa version dégradée. Nous construisons un modèle d'espace d'état probabiliste tenant compte de la régularité entre les noyaux de flou voisins et proposons un nouvel algorithme basé sur le filtrage de particules bootstrap pour échantillonner des particules pondérées, et ainsi construire la distribution postérieure recherchée. Des expériences numériques sur divers exemples de flou spatialement variants illustrent les avantages et la stabilité de notre approche. Ensuite, nous nous concentrons sur le problème de déconvolution aveugle d'image, consistant à estimer conjointement l'image et le noyau de convolution à partir d'une image bruitée et floue. Nous adoptons une approche bayésienne variationnelle, pour construire une approximation appropriée de la distribution postérieure. Nous introduisons des étapes de majoration pour permettre des mises à jour explicites des variables même pour des a prioris non conjugués et un bruit non gaussien. Cela donne l'algorithme bayésien variationnel. Nous proposons ensuite de dérouler cet algorithme sur des couches de réseaux de neurones, en suivant le paradigme de déroulement profond récemment introduit dans la littérature. L'algorithme obtenue bénéficie d'un nombre réduit de paramètres à régler, d'une exécution rapide sur carte GPU et fournit des meilleurs résultats de restauration. La supériorité du la méthode est illustrée sur trois jeux de données impliquant des images naturelles en couleur ou niveaux de gris et diverses formes de flous.

Entrée

Motivation

The problem of restoring an image from a degraded, blurred and noisy version thereof arises in many application frameworks, such as microscopy, medical imaging and astronomy. Non-blind image deblurring reads as a linear inverse problem requiring to invert the effect of a blur operator, in the presence of noisy measurements. Simple (pseudo)inversion of the blur degradation operator might seem as an appealing solution approach, but in most cases it involves computing the inverse of a very high-dimensional, ill-conditioned or even non-square matrix, and whenever it exists, the result of the direct inversion is rarely satisfactory. In the case of blind restoration, the problem is even more complicated since the blur operator is also unknown, and we observe that the degraded image from which the image and the blur are to be estimated. Many approaches have been proposed for image restoration. Traditional methods rely on optimization approaches, minimizing an appropriate cost function by a suitable iterative optimization algorithm. These approaches are reliable and demonstrate satisfactory results. However, finetuning the hyperparameters of the cost function (e.g., regularization parameter) or of the optimization solver (e.g., stepsize) can be tedious. Supervised approaches based on deep learning have also been proposed for image restoration. They rely on complex nonlinear models whose parameters are learned from image databases. These methods are flexible and can achieve excellent results assuming a training database of sufficient quality and size. Despite their inherent model complexity, these approaches can be implemented efficiently on GPU cards, thanks to suitable programming environments (for example, PyTorch). In addition to requiring large database, another disadvantage is the lack of interpretation of these methods. Both optimization and learning techniques generally offer a point estimate of the solution. In particular, they rarely allow to extract higher order statistics, and therefore cannot quantify the uncertainty on the result obtained.

In this thesis, we consider a third family of approaches, called the Bayesian approach, to deal with the image restoration problem. Compared to the methods mentioned previously, the Bayesian framework considers the sought parameters (image/blur in our case) as random quantities rather than deterministic ones. The goal becomes the estimation of their distribution a posteriori, given the observed data (i.e., the degraded image). The application of the Bayesian framework has two main advantages. First of all, the noise statistic and a priori information on the unknown quantities can be taken into account explicitly, via the Bayes' formula. Second, since the entire a posteriori distribution is estimated, it is easy to obtain any moment, such as the mean and covariance, or a confidence interval for further decision making.

Nevertheless, a complex task lies in the implementation of the Bayesian approach, namely the calculation of the a posteriori distribution. Indeed, in most Bayesian formulations for image restoration, the normalization constant of this distribution cannot be analytically computed even if the model a priori and the observation model are simple (e.g., Gaussian models). Many algorithmic methods have been considered to provide a good approximation of the desired a posteriori distribution. In this thesis, we make contributions to this question, in the context of two problems arising from image restoration, namely (i) non-stationary blur identification, and (ii) blind deconvolution. For the first problem, we propose a sequential stochastic simulation approach based on the particle filtering framework (Chapter 3). For the second problem, we propose an approximation method through a variational Bayesian technique (Chapter 4). Then we build a supervised version of this method using a neural network unfolding methodology (Chapter 5). Our methodological contributions are supported by a large number of numerical results, involving various models of blur and noise. The detailed contributions of each chapter are described below.

Main contributions

Chapter 3: Particle filtering method for space-varying blur identification Images produced by optical instruments often suffer from blur caused by light diffraction or object motion. The presence of the blur causes the spread of an infinitesimal point source in the acquired image defining the "Point Spread Function" (PSF). When the PSF is unknown, a blind deconvolution strategy can be used to restore the image. Another approach is to firstly identify the blur and then remove it from the degraded image by a non-blind restoration strategy. This method requires good precision for the estimation of the PSF. This can be effectively achieved by a prior step of acquiring standardized and calibrated objects, such as fluorescent spherical microbeads in microscopy or resolution diagrams for calibration of digital cameras. The identification of the PSF is typically formulated as a (constrained) least squares minimization problem. The identified PSF can be further used to fit a parametric model (often nonlinear) to determine the characteristics of an optical system. In most realistic scenarios, a stationary PSF model is not suitable due to depth of field, object or camera movement, anisotropic optical aberrations of the lens, or atmospheric turbulence. Such sources of image degradation give rise to so-called spatially variant blur. Identifying the PSF in such context requires estimating a map of several PSFs, describing

Main contributions

the blur kernel at each location in the optical plane. In this chapter, we address the problem of the Bayesian estimation of the parameters of a spatially variant PSF from calibrated image acquisitions. We consider a parametric sequential model for the map of PSFs which assumes small variations between PSFs acting on neighboring regions of the image. We formulate a state-space model where each time step corresponds to a patch location in the image. We adopt a probabilistic approach aiming to produce the a posteriori distribution of the unknown parameters. This allows us to deal with uncertainty in a systematic way as well as the introduction of a priori knowledge about the unknowns. Consequently, we are able not only to provide measures of uncertainty on the parameters to be estimated, but also to propagate this uncertainty to tasks where these estimated parameters are used (here, the restoration task). In our approach, the probability density a priori function of the parameters takes a Markovian structure, which allows us to exploit existing inference approaches for sequential inference. The Bayesian recursions obtained do not give rise to explicit forms and we resort to particle filtering (PF) for an approximate resolution of the inference problem. We start with the classic method bootstrap PF (BPF). BPF manages to reach satisfactory results in some of our examples but it has some shortcomings, for example it fails in scenarios with low observation noise where a problem of particle degeneracy arises. We therefore propose a new method called Generalized Interacting Annealed PF (GIAnPF). GIAnPF handles particle degeneracy by considering a sequence of intermediate distributions, which are "annealed" versions of the initial filtering distribution. We demonstrate the good performance of the proposed algorithm as well as its robustness through several numerical experiments, including three different forms of spatially variable parametric fuzziness. We also illustrate the validity of the estimates obtained, and the interest of uncertainty quantification, to solve image restoration problems.

Chapter 4: Variational Bayesian algorithm for blind image deconvolution under Gaussian and non-Gaussian noise

We consider the problem of blind deconvolution, where a blurred and noisy image must be restored without assuming the knowledge of the blur kernel. We propose in this chapter a variational Bayesian algorithm aiming to provide an optimal approximation (in the sense of the Kullback-Leibler divergence) of the a posteriori distribution of the parameters to be estimated (here, the image, the blur, and hyperparameters of the model). We assume the PSF is stationary and satisfies a linear equality constraint (e.g., sum-to-one, axial symmetry). We also adopt an a priori model of the SAR type on the PSF, whose mean is defined in advance, and a total-variation (TV) regularization on the image. We then seek to minimize the Kullback-Leibler (KL) divergence between the a posteriori distribution obtained and a separable model (of the "mean field" type) of it, by an alternating minimization algorithm. Even in the simpler case of Gaussian noise, the updates do not have closed forms and we therefore propose a data augmentation strategy inspired by semi-quadratic methods used in optimization. This yields our algorithm VBA (variational Bayesian algorithm), for which we give a first form adapted to a Gaussian noise, then an extension of it to a more general noise model. The algorithm obtained is fast and efficient. Nevertheless, preliminary tests have shown its high sensitivity to the tuning of hyperparameters such as the weight of the regularization on the kernel, and the variance of the noise. This leads us to our next contribution, presented in the next chapter.

Chapter 5: Unrolled variational Bayesian algorithm for blind image deconvolution

We focus here on the blind image deconvolution problem in the case of Gaussian noise. We propose the unfoldedVBA algorithm, which applies the deep unfolding paradigm, to the VBA approach described in the previous chapter. VBA iterations are embedded into a deep neural network architecture, with each iteration becoming a layer of the network. This allows (i) to learn the hyperparameters (in particular, the noise level) of VBA in an automatic and supervised way, (ii) to improve the quality of the results by training the model with a cost function linked to the metric image quality, and (iii) to implement the method by taking advantage of the available GPU resources, and therefore to considerably reduce the image restoration time during the test phase. Unlike standard deep learning methods for blind deconvolution, our method is interpretable, thanks to the unrolling paradigm. Note that variational Bayesian methods have already been used in a deep learning context. In particular, they are standard methods for training variational auto-encoders, and Bayesian neural networks. However, to our knowledge, our work is the first to investigate the deep unrolling of a variational Bayesian technique. In addition to the VBA layers, we integrate additional layers into our architecture, allowing us to learn the hyperparameters and to improve again the quality of the output image via post-processing. We show the superiority of our approach compared to the state-of-the-art on several datasets, involving various blurs (e.g., Gaussian blurs, real-world motion blurs), and grayscale/color natural images.

Publications

Outline

• Chapter 1 presents the motivation and context of our research, and introduces our contributions.

• Chapter 2 gives an introduction to the field of image restoration, with a specific focus on Bayesian techniques.

• Chapter 3 presents our first contribution, that is the proposition of a sequential particle filtering algorithm for the probabilistic identification of the parameters of spatially-variant blurs.

• Chapter 4 introduces our second contribution, that is the construction of two variational Bayesian algorithms, to deal with image blind deconvolution problems under Gaussian and non-Gaussian noises.

• Chapter 5 states our third contribution, which is the introduction of the unrolling paradigm within the VBA method proposed in the previous chapter, and its experimental validation onto various datasets of blind image deblurring.

• Chapter 6 concludes the thesis and raises some perspectives. The so-called Bayesian approach in statistical analysis is based on Bayes' rule. The main advantage of Bayesian approach is that instead of deriving only point estimate of unknowns like other methods (e.g., optimization-based), it provides a full distribution and thus yields, for instance, uncertainty quantification of the estimation [Jeffreys 1939, Robert et al. 2009, Mattei 2020]. Bayesian methods have been applied to different problems related to image processing. The rest of this chapter is organized as follows. In Section 2.1, we introduce the Bayesian paradigm, and provide typical Bayesian formulations encountered in image restoration. In Section 2.2 and Section 2.3, we review two main classes of Bayesian inference methods (static and dynamic ones) respectively, which will be at the core of this thesis developments. In Section 2.4, we conclude this chapter.

Bayesian framework for imaging science

Bayesian paradigm

In the field of statistical inference, there are two mainstreams called frequentist [Bailer-Jones 2017] and Bayesian [Bernardo & Smith 1994, Gelman et al. 2004, Robert 2007] statistics. For some unknown parameters to be estimated, the frequentist approach assumes them to be fixed and uses only the information from the current experiment. The Bayesian or probabilistic approach allows to quantify the uncertainty and to incorporate prior information. We first review some basic concepts of Bayesian approaches.

Chapter 2. Background

Let Θ be the collection of unknown parameters to be estimated, and y the given data. Both of them are random variables, thus we can get access to the posterior distribution of the unknowns p(Θ|y) according to the Bayes' rule, given the likelihood p(y|Θ) and the prior distribution p(Θ) as follows:

p(Θ|y) = p(Θ)p(y|Θ) p(y) . (2.1)
Hereabove,

p(y) = p(Θ)p(y|Θ)dΘ (2.2)
is the marginal density of the variable Θ and can be seen as a normalization constant of p(Θ|y). This means that we have the proportionality relation:

p(Θ|y) ∝ p(Θ)p(y|Θ). (2.3)
Eq. ( 2.3) is especially useful since in most cases the constant (2.2) cannot be computed analytically.

MAP/MMSE estimator

Let us review some popular point-wise estimators. Since the exact value of the parameter Θ is unknown, we need to resort to estimation methods. The relevant maximum likelihood estimate (MLE) estimator is derived as

ΘMLE = argmax Θ p(y|Θ), (2.4)
by obtaining the estimate of the unknown parameters Θ only based on the observation model p(y|Θ). Moreover, if we have access to the prior distribution p(Θ), we can then estimate it by the posterior distribution p(Θ|y).

When maximizing this posterior distribution p(Θ|y), we obtain the maximum a posteriori (MAP) estimator, ΘMAP = argmax Θ p(Θ|y).

(2.5)

Note that if the prior distribution, p(Θ), is chosen as uniform, the MAP estimator coincides with the MLE estimator since

ΘMAP = argmax Θ p(Θ|y) = argmax Θ p(y|Θ)p(Θ).
(2.6)

Equivalently, ΘMAP = argmin Θ (ln p(y|Θ)ln p(Θ)) .

(2.7)

However, the closed-form solution is intractable in most applications, so we resort to iterative optimization methods to compute the estimator numerically.

On the other hand, we can also choose to minimize the mean square error between the estimation Θ and the true value Θ,

ΘMMSE = argmin Θ E Θ [( Θ -Θ) 2 ].
(2.8)

Thus the Bayesian estimator is the mean of the posterior distribution of Θ, ΘMMSE = E Θ|y (Θ) = Θp(Θ|y)dΘ.

(2.9)

This posterior mean estimator is often called minimum mean squared error (MMSE) estimator of the parameter Θ. Note that if the posterior distribution is Gaussian, these two Bayesian estimators coincide since the mean and the mode of a normal distribution are identical. However, in most cases, as soon as the posterior distribution is not symmetric, MAP estimator and MMSE estimator would give different results. In general, both MAP estimator and MMSE estimator are point estimates instead of Bayes estimators unless Θ is discrete variable, which also demand a credible interval for the estimation. In conclusion, this shortage of point-wise estimators forces to find alternative Bayesian approaches with the ability to provide the uncertainty of the unknown parameters.

Image restoration

Image restoration is a key problematic of the field of computer vision. Blur is a main cause of image degradation that commonly appears during the process of acquisition of the images due to many reasons such as camera shake [Fergus et al. 2006], fast moving of the object [Rozumnyi et al. 2021], and out of focus of the object [Chan & Nguyen 2011]. Undoubtedly, the blur deteriorates valuable image information [Shi et al. 2014], thus image deblurring becomes a hot topic in many areas, for example, astronomical imaging [Murtagh et al. 2007], microscopy [START_REF] Holmes | [END_REF]] and medical imaging [START_REF] Michailovich | [END_REF], where the images are usually required to be in high quality. We display some examples of images degraded by different kinds of blurs in Fig. 2.1. On top of blur degradation, images also suffer from noise coming from the sensor [START_REF] Boyat | [END_REF], for example, the popular additive Gaussian noise [Russo 2003] is a side effect of sensor heat. The salt and pepper noise [Alqadi 2018], also known as impulse noise [Jayamanmadharao et al. 2010], is caused by sharp and sudden disturbances in the image signal. Poisson noise [START_REF] Bertero | [END_REF]] appears because of the effects of electromagnetic waves. Multiplicative speckle noise [Guadalupe Sanchez et al. 2012] is present in low level luminance images such as ultrasound images. We display some examples of images corrupted by different kinds of noise in Fig. 2.2.

In image deblurring (also called image deconvolution), given a degraded, blurred, and noisy image, the goal is to restore a clean image along with an estimation of corresponding blur kernel (if the latter is unknown). Mathematically, the observation model is defined as follows: where y ∈ R N is the degraded image that is observed; x ∈ R N is the clean image to be restored; n ∈ R N corresponds to some additive random noise (more sophisticated noise models could arise, see [START_REF] Marnissi | [END_REF], Dupe et al. 2008]); and H ∈ R N×N is a linear operator related to the blur kernel h ∈ R M . According to whether the blur kernel is known or unknown, one can divide the image restoration problem into two different categories: non-blind image deconvolution or blind image deconvolution. In the former case, the blur kernel is known and the image must be restored given the blurry noisy image and the blur kernel [START_REF] Marnissi | [END_REF], Sun et al. 2014, Schuler et al. 2013].

y = Hx + n, ( 2
As for the blind deconvolution problem, one needs to restore the image and estimate the blur kernel given the observed degraded image [START_REF] Kundur | [END_REF], Huang et al. 2021a, Likas & Galatsanos 2004]. Let us make a special focus on this problem that will be at the core of our thesis. There are mainly two approaches to tackle image blind deconvolution. The first one is called sequential identification process [Huang et al. 2021a, Carasso 2001] where the blur kernel is derived first through some kernel estimation methods applied on calibrated images, then non-blind image deconvolution is conducted to restore the image using the previously approximated blur kernel. This problem will be investigated in Chapter 3. A second method is the joint estimation approach [Babacan et al. 2009, Ren et al. 2020], which estimates the image and blur kernel simultaneously. Since this inverse problem is severely ill-posed, an infinite number of combinations of the pairs (H, x) might exist. This demands to adopt suitable prior assumptions on the sought image and on the sought blur kernel. We will tackle this problem in Chapters 4 and 5.

In Chapter 3 of this thesis, we will focus on the blur identification problem. This amounts to identifying an estimate ĥ of an original blur kernel h, the given y, and the original image x. In the case of stationary blur, the observation model is

y = X h + n, (2.11)
where X is the suitable block Toeplitz matrix related to the original image x, and n models some noise. In calibration, it is useful to introduce a parametric model for the blur kernel depending on unknown parameters Θ,

h = h( Θ).
(2.12)

The goal is thus to find estimate for Θ given the knowledge of the pair (x, y). In Chapter 3, we will consider the more challenging case of a spatially-variant blur, requiring the estimation of a set of some parameters Θ t , associated to image positions t = 1, • • • , T . A complete review of the state-of-the-art for this problem will be provided in Chapter 3.

In Chapter 4 and 5, we focus on the joint image blind deconvolution problem. Given the observation model (2.10), we seek for an estimate (x, ĥ) of the pair of variables (x, h) given y. Following a Bayesian approach, our target is to find an entire distribution describing the sought parameters (i.e., image, kernel, hyperparameters). Many approaches have been proposed to solve this task, a complete review will be provided in Chapter 4.

For both calibration and blind deconvolution problems, a key aspect is the definition of suitable priors for the unknowns. One can impose smoothness on the image by applying total-variation (TV) [START_REF] Marnissi | [END_REF], Babacan et al. 2009] or simultaneous autoregressive (SAR) prior [START_REF] Tzikas | [END_REF], Shi et al. 2017]. For specific images such as US images, one can rely on a sparsity prior [Zhao et al. 2016, Besson et al. 2017]. For the blur kernel, it is usually assumed to satisfy a normalization constraint (e.g., sum-to-one) [START_REF] Zhou | [END_REF] to avoid the scaling ambiguity, as well as a non-negativity constraint [Mateos et al. 2016]. Sometimes, for special shapes of blurs, symmetry might be also assumed, i.e., the values of the blur are symmetrical according to the anti-diagonal line of 2D blur kernel, as shown in Fig. 2.3. The estimation of spatially-variant blurs can also benefit from priors accounting for their smooth variations in the image space [Denis et al. 2015

Bayesian inference methods

In the following, we review two families of algorithms that perform approximate Bayesian inference, namely variational Bayesian algorithms (Section 2.2.1) and importance sampling (Section 2.2.2). Let us point out that another important family of methods are Markov chain Monte Carlo (MCMC) methods that approximate the target distribution by relying on Markov chain theory and Monte Carlo integration [Robert & Casella 2004[START_REF] Gamerman | [END_REF], Robert et al. 2018]. We refer the reader to [START_REF] Pustelnik | [END_REF] for a review of the usage of these techniques in the field of image restoration. Examples of applications in the context of image deblurring can be found in [Marnissi et al. 2020, Marnissi et al. 2018].

Variational Bayesian algorithm

The variational Bayesian algorithm (VBA) attempts to find a more tractable approximation p(Θ) of the unknown distribution p(Θ|y). This is a powerful tool especially when we want to marginalize the unknowns without explicitly computing them because the involved integral is intractable. The variational Bayesian paradigm mainly has two advantages, first

Bayesian inference methods

it provides us an analytical approximation of the posterior distribution of the unknowns.

Secondly it enables us to get access to the high-order moments of the unknowns that allows us to obtain the uncertainty of the estimation, which is essentially important for a Bayesian method. Many applications of variational Bayesian principle have been arising to solve image deblurring problems [START_REF] Marnissi | [END_REF], Likas & Galatsanos 2004, Babacan et al. 2009, Zheng et al. 2015c]. Among different types of approximation, mean field approximation is the most used one that makes use of a factorized structure of q(Θ) to approximate the posterior distribution,

p(Θ|y) ≈ q(Θ) = R ∏ r=1 q r (Θ r ), (2.13)
with R being the number of factorized unknowns in the set of Θ. Then, one searches for a distribution q(Θ) with form (2.13) that approximates the posterior distribution p(Θ|y) [START_REF] Blei | [END_REF]]. The quality of the approximation is evaluated using the Kullback-Leibler (KL) divergence, defined as

KL(q(Θ)||p(Θ|y)) = q(Θ) ln q(Θ) p(Θ|y) dΘ (2.14)
where the equality only holds when q(Θ) = p(Θ|y). It is a non-symmetric measure of the difference between the target posterior distribution and its approximation. The optimal distribution q(Θ) is the one that minimizes the KL divergence. One can show that the solution of each q r (Θ r ) minimizing the KL divergence reads, q r (Θ r ) ∝ exp < ln p(Θ, y) > ∏ i =r q i (Θ i ) .

(2.15)

The so-called VBA consists in alternating over the updates of the factors q r , using the expression (2.15), while keeping others unchanged. An application of this methodology to image blind deconvolution will be presented in Chapters 4 and 5.

Importance sampling and adaptive importance sampling

Importance sampling (IS) is a Monte Carlo methodology for the approximation of intractable integrals [Robert & Casella 2004, Liu 2004]. While IS is a useful approximation technique in many settings, here we focus again on the Bayesian inference problem. Following the already defined notation, the posterior distribution is

p(Θ|y) = p(y|Θ)p(Θ) p(y) ∝ π(Θ), (2.16)
where π(Θ) = p(y|Θ)p(Θ) and p(y) = π(Θ)dΘ is a constant. In many applications, we are interested in computing a moment w.r.t. the posterior distribution,

I = E Θ|y [g(Θ)] = g(Θ)p(Θ|y)dΘ (2.17) = 1 p(y) g(Θ)π(Θ)dΘ, (2.18)
where g is any integrable function w.r.t. p(Θ|y). However, Eq. (2.18) cannot be computed, since in most cases, this integral is intractable and p(y) is also not available. IS [Elvira & Martino 2021] is a family of methods that generate samples from the so-called proposal or importance distribution and assign each sample with an importance weight, accounting for the mismatch between target and proposal distributions. The main difference between this approach and direct Monte Carlo sampling is that, in IS, the particle approximation of the distribution is composed of weighted Dirac masses. The standard IS algorithm is shown in Table . 2.1.

Table 2.1: IS algorithm 1. Draw M i.i.d. samples from the importance distribution q(Θ),

Θ (m) ∼ q(Θ), m = 1, • • • , M (2.19)
2. Compute the associated importance weights

w (m) = π(Θ (m) ) q(Θ (m) ) , m = 1, • • • , M (2.20)
and compute the normalized weights over all the samples,

w (m) = w (m) ∑ M i=1 w (i) , m = 1, • • • , M (2.21)
3. Output: It provides M weighted particles {Θ (m) , w (m) } M m=1 for the approximation of the posterior distribution.

IS methods allow for the approximation of the expectation above via the selfnormalized importance sampling (SNIS) estimator,

Ĩ = M ∑ m=1 w (m) g(Θ (m) ).
(2.22)

IS-based methods are powerful approximation techniques, but their performance strongly relies on the choice of the proposal distribution. Therefore, many approaches are proposed in order to tackle this challenge. One alternative is the usage of several proposals, under the so-called multiple IS (MIS) methodology [START_REF][END_REF], Elvira et al. 2016a, Sbert et al. 2018, Elvira et al. 2019b]. Adaptive importance sampling (AIS) goes one step further and consists in iteratively adapting a set of proposals (see [START_REF] Bugallo | [END_REF] for a review). One of the main families of adaptive importance sampling methods is the population Monte Carlo (PMC) algorithms [Iba 2001, Cappé et al. 2004, Cappé et al. 2008, Elvira et al. 2017a, Elvira et al. 2017b]. In PMC, the location parameters of the proposals are adapted following resampling-based schemes. Few PMC methods adapt also the scale parameters, e.g., O-PMC [START_REF] Elvira | [END_REF]. Other families of AIS methods include the moment-matching AMIS algorithm [Cornuet et al. 2012, El-Laham et al. 2019b, El-Laham et al. 2019a] or the layered AIS algorithm (e.g., LAIS [Martino et al. 2017b, Martino et al. 2017a, Mousavi et al. 2021]).

As a result, AIS methods can approximate the expectation of any function g(•) as

Ĩt = t ∑ τ=1 M ∑ m=1 ρ (m) τ g(Θ (m) τ ).
(2.23)

where

ρ (n) t = w (n) t ∑ t τ=1 ∑ M m=1 w (m) τ , t = 1, • • • , T, n = 1, • • • , M.
(2.24)

Note that this estimator re-normalizes all importance weights over iterations. Other alternatives for combining the partial estimators of each iteration are provided in [Nguyen et al. 2015] and [START_REF] Cappé | [END_REF]].

Sequential Bayesian inference

In this section, we address the Bayesian inference problem in state-space models, a class of probabilistic models that allow for sequential processing of the observations. In the following, we describe this class of models and particle filtering (PF), a family of Monte Carlo methods for the approximation of posterior distributions in this setting.

State-space models

State-space models (SSMs) allow to describe dynamic systems that evolve over time/space. More precisely, a general probabilistic SSM is given by where p(θ 0:T ) is the prior distribution, p(y 1:T |θ 0:T ) is the likelihood distribution of the observation model and p(y 1:T ) is the normalization constant defined as p(y 1:T ) = p(y 1:T |θ 0:T )p(θ 0:T )dθ 0:T .

θ 0 ∼ p(θ 0 ), (
(2.29)

Chapter 2. Background

The previous problem is called the smoothing problem and requires batch processing of all observations. The main disadvantage of the full posterior approach in (2.28) is that (a) we have to reprocess all other observations when we obtain a new observation, and (b) the dimension of the inferential problem grows with t. Since computing the full joint posterior distribution at all time steps is computationally expensive and unnecessary in real-time applications, we consider building the sequence of the so-called filtering posterior distributions instead, defined as where g(•) : R n → R m is an arbitrary function.

p(θ t |y 1:t ), t = 1, • • • , T, ( 2

Particle filtering

The computation of these filtering distributions is called the optimal Bayesian filtering problem. Unfortunately, for most realistic state-space models, these distributions cannot be obtained analytically, thus we resort to their approximations. PFs, also called sequential Monte Carlo (SMC) methods, are Monte Carlo algorithms that are well suited for approximate Bayesian inference in SSMs. PFs are based on IS, reviewed in Chapter 3, and they approximate the distributions (and integrals) of interest through weighted particles. Table 2.2 describes the bootstrap PF (BPF) algorithm, the most standard and widely used PF.

Step 2(a) describes the simulation of the particles and Step 2(b) details the computation of the normalized weights. While not being essential, Step 2(c) implements a resampling step that attenuates the particle degeneracy problem (see for instance [START_REF] Djurić | [END_REF]).

0 } M m=1 from the prior p(θ 0 ). 2. Filtering step. For t = 1, ..., T :

(a) Simulate θ (m) t ∼ p(θ t |θ (m) t-1 ), m = 1, ..., M (2.32) (b) Compute the normalized weights by w (m) t ∝ p(y t |θ (m) t , X t ), m = 1, ..., M (2.33) (c) Resample M times from {θ (m) t } M m=1 with associated probabilities {w (m) t } M m=1 , i.e., for m = 1, ..., M, let θ (m) t = θ ( j) t with probability w ( j) t , j = 1, ..., M 3. Output: The filtering step provides M weighted particles {θ (m) t , w (m) t } M
m=1 at each t, for the approximation of the filtering distribution.

As a result of the BPF, for each time step t, we can approximate the filtering distribution p(θ t |y 1:t ) by the M generated particles (θ (m) t ) 1≤m≤M and associated normalized weights

(w (m) t ) 1≤m≤M computed in Step 2(a)-(b) as follows: p(θ t |y 1:t ) ≈ M ∑ m=1 w (m) t δ (θ t -θ (m) t ).
(2.34)

Then, the approximation of the expectation of any arbitrary function g(•) is

E[g(θ t )|y 1:t ] ≈ M ∑ m=1 w (m) t g(θ (m) t ).
(2.35)

For example, we can compute the mean and the second order moment of the posterior distribution by taking g(θ) as θ and θ 2 , respectively.

The BPF is just an instance of particle filtering, but many other filters have been proposed. For instance, auxiliary particle filters (APFs) were originally proposed [Pitt & Shephard 1999] and have gained popularity in the last years due to their flexibility to implicitly adapt the multiple importance sampling proposals [START_REF] Klaas | [END_REF], Elvira et al. 2018a, Elvira et al. 2019a[START_REF] Branchini | [END_REF]. Also, multiple PFs [Djurić et al. 2007], Gaussian PFs [Kotecha & Djurić 2003], and Rao-Blackwellized PFs [Särkkä et al. 2007] have shown to provide superior performance in several applications. In all cases, the performance of the filter depends on the number of particles, M. The adaptation of M has been recently proven to be essential to attain successful approximations in challenging problems [Elvira et al. 2016b, Elvira et al. 2021]. Finally, note that the filtering problem is challenging whenever the observations are very informative, which leads to a steep or peaky likelihood. In order to attenuate this issue, [Deutscher et al. 2000] proposed an annealed particle filter by performing one annealing run at each time step. Another related work is [START_REF] Gall | [END_REF], which combines a PF and the interacting annealing Chapter 2. Background algorithm. The application of this methodology to the problem of blur identification will be presented in Chapter 3.

Conclusion

In this chapter, we have provided an overview of the Bayesian paradigm and introduced an inference problem of interest in the field of image restoration. Then we have described two main classes of Bayesian inference strategies, namely static and dynamic ones, and the main methods for each of these classes. This thesis will mostly focus on VBA and BPF strategies. Our goal is to bring novel methodological contributions in these two methods, so as to derive efficient Bayesian-based strategies for solving two major problems of image restoration, namely blur identification and blind deconvolution.

Introduction

Images produced by optical instruments often suffer from blur caused by light diffraction or object motion. The presence of the blur provokes an infinitesimal pointsource to be spread in the acquired image defining the so-called Point Spread Function (PSF). When the PSF is unknown, one can resort to a blind deconvolution strategy [START_REF][END_REF], Carasso 2001, Kenig et al. 2010, Reeves & Mersereau 1992] to jointly retrieve the image and blur. Another approach consists in first identifying the blur and then removing it from the degraded images using a non-blind restoration strategy [START_REF] Pustelnik | [END_REF], Berisha & Nagy 2014]. This demands a good accuracy for Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification the PSF estimation, which can be efficiently reached by a preliminary acquisition step of normalized and calibrated objects, such as fluorescent spherical microbeads in microscopy [Kirshner et al. 2012b, Hortholary et al. 2021] or resolution charts in digital camera calibration [Mannan & Langer 2016, Bell et al. 2016]. The PSF identification problem is typically formulated as a least-squares one [Tezaur et al. 2015, Mourya 2016].

On top of serving for image restoration purposes, the identified PSF can also be fitted into a parametric non-linear model in order to determine characteristics of the optical system [Anthony & Granick 2009, Chouzenoux et al. 2019, Kirshner et al. 2012a, Xue & Blu 2014]. In most realistic scenarios, the stationary PSF model is not suitable due to extended depth of field [Trouvé et al. 2011, Zhang et al. 2016], moving objects or cameras [START_REF] Zhe | [END_REF], Zhang & Hirakawa 2016], anisotropic optical lens aberrations [START_REF] Jezierska | [END_REF], or atmospheric turbulence [START_REF] Li | [END_REF], Schmitz et al. 2020, Chimitt & Chan 2020]. Such image degradation sources give rise to a so-called spatially variant blur [Chakrabarti et al. 2010, Hirsch et al. 2010]. The PSF identification then requires the recovery of a PSF map, describing the blur kernel at each location of the spatial plane [Tezaur et al. 2015].

In this chapter, we address the problem of estimating in a sequential manner the parameters of spatially variant PSFs from calibrated image acquisitions. We consider a flexible piecewise constant parametric model for the space-varying PSF map that allows us to describe smooth variations among PSFs acting on neighbor regions. We formulate a statespace model where each time step corresponds to a different patch location. We adopt a probabilistic approach, aiming at producing posterior distributions of the unknowns. The probabilistic approach allows for dealing with the uncertainty in a systematic manner and the inclusion of prior knowledge about the unknowns. As a consequence, we are able not only to provide uncertainty measures on the unknown parameters, but also to propagate this uncertainty to useful tasks where those estimated parameters are used (e.g., in deblurring). In our approach, the prior probability density function (pdf) of the parameters has a Markovian structure, which allows us to inherit existing inference approaches for the sequential inference. Note that more complicated local dependencies could be modeled, e.g., via Markov random fields, at the expense of losing the aforementioned simplicity in the inference task. Our flexible formulation includes the consideration of potential non-standard observation and transition models. More precisely, we can operate virtually with any nonlinear and non-Gaussian model. Due to the intractability of the Bayesian recursions, we resort to PF for an approximate inference solution. We start by considering existing PF methods, such as the well-known BPF [Gordon et al. 1993]. While the BPF can obtain adequate results in some models, it presents some deficiencies, for instance in scenarios with low observation noise, where the peaky likelihood challenges the diversity among particles (see for instance [START_REF] Djurić | [END_REF]] about the particle degeneracy problem). Due to these limitations, we propose a new PF method that we call Generalized Interacting Annealed PF (GIAnPF). The new PF tackles the particle degeneracy by considering a sequence of intermediate distributions, that are annealed versions of the filtering distribution. The considered sequential Bayesian framework provides us with three key advantages, namely (i) low computational cost and limited memory load due to a sequential processing; (ii) a flexible choice of the state-space models that enables us to consider non-linear relations and non-Gaussian noise; and (iii) a measure of statistical uncertainty on the estimated parameters of kernels. The sequential approach can be exploited to build online implementations in applications involving very large scale images available in batches. Dealing with the image in a sequential manner also alleviates the need of large memory capacity, since the information of previously processed patches does not need to be stored. We demonstrate the good performance of the proposed approach as well as its robustness through several sets of numerical experiments for three representative parametric blur shapes. We also illustrate the validity of the resulting estimations, to solve image restoration problems.

The rest of Chapter 3 is organized as follows. Section 3.2 introduces the problem of spatially variant blur identification. We review the literature around this topic and discuss the construction of parametric blur models. Section 3.3 presents the state-space models for parametric PSF estimation and proposes our algorithm for Bayesian inference. Section 3.4 shows abundant experimental validation, and Section 3.5 closes the section with some final remarks.

Problem statement

Let us assume the clean image is x ∈ R N and its corresponding blurred and noisy image y ∈ R N that is corrupted by spatially variant blur and noise with given level. Each pair of image x, y ∈ R N are decomposed into a partition of T patches (x t ) T t=1 ∈ R P , PT = N, and (y t ) T t=1 ∈ R P . We assume that each patch of the observed image is given by:

y t = X t h t + n t . (3.1)
Hereabove, y t ∈ R P represents the blurry noisy patch, h t ∈ R L is the blur kernel to be estimated at patch t, and n t ∈ R P models an additive noise. Moreover, X t ∈ R P×L is a suitable matrix related to the blur of a patch t of original image x. The model in Eq. (3.1) identifies with the PSF-interpolation model from [Denis et al. 2015, Eq. (22)], where the interpolation strategy depends on the choice made for X t (see also [START_REF] Hadj | [END_REF], Hirsch et al. 2010]). In all our practical experiments, for simplicity and complexity reasons, we will adopt the piecewise constant PSF model from [Denis et al. 2015, Eq.20], that is the order 0 (i.e., nearest neighbor) PSF-interpolation model. Note that several other strategies are available for modeling the space-varying blur operator, for instance in [Escande & Weiss 2015, Denis et al. 2015, Sroubek et al. 2016, Hirsch et al. 2010, Nagy et al. 2004], with various modeling accuracy and computational costs. Our choice in Eq. (3.1) is particularly well suited for fast blur identification due to the separability over patches. This model was also considered in [START_REF] Bar | [END_REF], Hadj et al. 2014, Mourya 2016, Tezaur et al. 2015].

In many practical situations, a prior over the blur shape is available. It is thus convenient to introduce a parametric model for those blur kernels, with the advantage of reducing the number of unknowns to estimate in the blur identification task. Each kernel is then assumed to read:

h t = h(ρ t ), t = 1, . . . , T (3.2)
Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification where h is a known function, not depending on t, that describes the general shape of the kernels, specified by parameters (ρ t ) T t=1 ∈ R K , with K ≥ 1 typically much lower than L. In this section, we focus on the problem of sequential parametric blur estimation, which amounts to retrieve, for every t = 1, . . . , T , an estimate for ρ t , given the knowledge of the past and present observations (X i , y i ) t i=1 and the function h(•).The main underlying assumption in our approach is the smooth variation of the kernel shape parameters when progressing sequentially along the patch indexes t = 1, . . . , T (see discussion in Section 3.3.2 about patch ordering). Therefore, our approach is particularly well suited for space-varying blurs arising from optical aberrations (e.g., phase aberration [Schuler et al. 2012], atmospheric aberration [Chimitt & Chan 2020]), and smooth motion (e.g., camera motion blur [Šorel & Flusser 2008]). One can also refer to the smooth varying PSF maps used as illustrative examples in [Escande & Weiss 2015, Denis et al. 2015]. As described in Section 3.3, the spatial dependence is encoded in the transition density among consecutive patches, which acts as a prior pdf. We note that, as it is often the case in Bayesian inference, vague enough priors are usually adequate particularly when enough informative observations are available. In contrast, this section does not cover PSF maps with abrupt changes, such as those considered in [START_REF] Zhe | [END_REF], Trouvé et al. 2011]. As we will show below, we exploit the similarity of kernels that are spatially close during the inference/estimation task. More precisely, unlike in other approaches, here the data of one patch is implicitly used to better estimate the parameter of other patches, which explains the good performance of the proposed modeling and methodology.

Related literature and contributions

The problem of space-varying blur identification has been widely studied in the literature of image processing. It is important to distinguish two types of approaches. First, there exists a bunch of methods, for performing image deblurring in the presence of an unknown spacevarying blur degradation, thus corresponding to blind image restoration [Zhang et al. 2016, Hadj et al. 2014, Chakrabarti et al. 2010, Zhang & Hirakawa 2016] in the case of motion blur. We also refer the reader to the recent work [Alam et al. 2019] for a review on this topic. Let us also mention [Hirsch et al. 2010] for the case of multi-frame blind image restoration. Most of these methods are focused on the restoration task, and not on the quantitative estimation of the blur map itself. This is at the exception of [Zhang et al. 2016, Chakrabarti et al. 2010, Zhang & Hirakawa 2016], that jointly restore the image and a piecewise-constant map, obtained through a segmentation-based strategy, describing the non-stationary blur effects. Though image restoration is a problem of high interest, in certain applications such as microscopy or astronomical imaging, an accurate qualitative and quantitative knowledge of the blur effects is key for a better understanding and improving (e.g., through calibration) of the imaging device [START_REF] Jezierska | [END_REF], Schmitz et al. 2020, Pankajakshan et al. 2012, Lefort et al. 2020]. For performing blur identification, it is necessary to make structural assumptions on the blur map to be estimated. For e.g., [Trouvé et al. 2011] assumes a finite set of candidate defocus blurs, to be tested in each location of the image, or [START_REF] Rajagopalan | [END_REF] considers parametric blurs depending only on a single parameter. The blur estimation can also be facilitated by making use of a calibrated image (e.g., fluorescent bead in microscopy, distant star in astronomical imaging) in order to reduce the ill-posedness of the identification problem. When the PSF is stationary, and no parametric model of it is further assumed, its estimation from calibrated image acquisition can be easily solved by a penalized least-squared algorithm (see for instance [START_REF] Chouzenoux | [END_REF] for an efficient method in the case of large size images). The problem raised in this case is similar to the one arising in high dynamic range image fusion, for which waveletbased methods have shown their efficiency [Tico & Pulli 2009, Tico et al. 2010]. If, additionally, a parametric model of the PSF is available, non-linear least squares method, such as Levenberg-Macquart [Kirshner et al. 2012a], or more recent proximity-based algorithms [Chouzenoux et al. 2019], can be employed for estimating directly the sought parameters. However, when PSF is non stationary, the accurate identification of its shiftvariant evolution becomes much more challenging to resolve, even when the image is calibrated (i.e., known). This problem is typically addressed using optimization-based methods [Tezaur et al. 2015, Mourya 2016]. We also refer the reader to [Schmitz et al. 2020] and references therein, for a review of the problem and recent insights, in the particular context of astronomical imaging. In the aforementioned works, no parametric model of the PSF was assumed though it is highlighted in [Schmitz et al. 2020] as a promising research direction. Moreover, most available techniques address the problem in a batch manner, requiring the loading of the full image before the starting of the identification process, which can be at the price of a high memory cost in high resolution imaging or even incompatible with an on-the-fly image acquisition (e.g., satellite imaging). Furthermore, up to our knowledge, no Bayesian-based techniques have been proposed so far, so that available strategies provide only point-wise estimators for the PSF field. We conclude this stateof-the-art review by emphasizing that several works exploit the smooth variability of the PSFs within the field of view, with the aim to reduce the cost of the non-stationary blur operator [Denis et al. 2015, Escande & Weiss 2015] or to perform PSF field interpolation [Ngolé Mboula & Starck 2017, Chen et al. 2005]. In this chapter, we propose a Bayesian sampling method, for the identification of parametrized space-varying blurs in the context of calibrated images. Our contribution is threefold: (i) the particle filtering strategy employed for sampling leads to a sequential and fast estimation result; (ii) the proposed method is versatile and easily adapted to diverse parametric blur shapes; and (iii) the Bayesian framework allows us to provide a posterior estimation, including uncertainty quantification, for the PSF map parameters.

Relevant parametric PSF models

The general model of Eq. (3.2) enables to encompass various blur shapes. Hereafter, we present three interesting classes of parametric blur models, that will be used thereby in the experiments of this chapter. To simplify the notation, let us ignore the subscript for patch t in this subsection. Without loss of generality, we consider kernels with square support parametrized by a grid with length √ L (assumed to be an odd integer), centered and regularly spaced, denoted

G = {- √ L-1 2 , - √ L-1 2 + 1, . . . , 0, 1, . . . , √ L-1 2 }.
The blur shape is determined by the values of a given function h(c 1 , c 2 ) for (c 1 , c 2 ) ∈ G 2 . Otherwise stated, (c 1 , c 2 ) are the (signed) distances, in pixels, from the origin (i.e., the center position of the squared kernel) in the horizontal and vertical axis, respectively. The vector h(ρ) ∈ R L is then simply deduced by ordering the entries of matrix [h(c 1 , c

2 ; ρ)] (c 1 ,c 2 )∈G 2 ∈ R √ L× √ L ,
following the lexicographic order 1 .

Generalized Gaussian blur space-varying generalized Gaussian blur shapes have been employed for instance in [Liu et al. 2021b] in the context of out-of-focus image deblurring. Such blur is parametrized by ρ = (θ , s) with θ ∈ R and s = (s 1 , s 2 ) ∈ (0, +∞) 2 the orientation and width parameters, respectively. The parametric model then reads:

(∀(c 1 , c 2 ) ∈ G 2 ), h(c 1 , c 2 ; ρ) = λ exp - 1 2 ([c 1 c 2 ]Σ -1 (θ , s)[c 1 c 2 ] ) p α p , (3.3)
with p > 0 and α > 0 the power and the scale of generalized Gaussian model, respectively, and λ the normalization constant reads

λ = p πΓ( 1 p )2 1 p α 1 |Σ(θ , s)| 1 2 , (3.4)
where Γ(•) is the Gamma function. The support size is:

s max = √ L 2α(-2 ln(a)) 1/p -1/2 . (3.5) Moreover, Σ(θ , s) = R θ s 2 1 0 0 s 2 2 R θ (3.6)
with R θ ∈ R 2×2 the rotation matrix defined by the orientation angle θ ,

R θ = cos θ -sin θ sin θ cos θ . (3.7)
We can ensure that the kernel value in its support corner is less or equal than a h(0, 0; ρ), for some a ∈ (0, 1), as soon as s 1 ≤ s max , s 2 ≤ s max with s max given in (3.5) above. Note that when p = α = 1 in (3.3), we recover the common Gaussian blur model, used for instance in fluorescence microscopy in [Zhang et al. 2007]. An example is displayed in Fig. 3.1(left).

Defocus blur Defocus space-varying blurs have been considered for example in [START_REF] Trouvé-Peloux | [END_REF] in the context of depth estimation in digital camera images. Defocus blur kernel can be parametrized by ρ = (θ , s) with θ ∈ R and s ∈ (0, +∞) 2 . The parametric model is, with the normalization constant λ defined as λ = πs 1 s 2 L.

(∀(c 1 , c 2 ) ∈ G 2 ), h(c 1 , c 2 ; ρ) = λ if (c 1 , c 2 ) ∈ C(θ , s) 0 otherwise , (
(3.9)

The support size is:

s max = max G = √ L -1 2 . (3.10)
and the elliptical domain:

C(θ , s) = (c 1 , c 2 ) ∈ G 2 such that (c 1 cos θ + c 2 sin θ ) 2 s 2 1 + (c 1 sin θ -c 2 cos θ ) 2 s 2 2 ≤ 1 . (3.11)
The latter can easily be ensured to be included into the kernel support, under the condition that s 1 ≤ s max , s 2 ≤ s max with s max given above. An example of such blur is displayed in Fig. 3.1(middle).

Skew-normal blur

The skew-normal blur shape [Azzalini & Valle 1996] was first considered in [Hansen & Jensen 2008] to build synthetic image deblurring problems. Its ability for modeling shift-variant PSFs in real digital imaging sensors was then illustrated in [Simpkins & Stevenson 2014]. This blur model depends of five positive parameters ρ = (s 1 , s 2 , w, α 1 , α 2 ), as follows:

(∀(c 1 , c 2 ) ∈ G 2 ), h(c 1 , c 2 ; ρ) = λ exp - 1 2 γ 1 γ 2 S -1 (w) γ 1 γ 2 × α 1 γ 1 +α 2 γ 2 -∞ exp - z 2 2 dz,
(3.12)

where

γ 1 = c 1 s 1 , γ 2 = c 2 s 2 , S(w) = 1 w w 1 , (3.13)
and the normalization constant λ is given by,

λ = 2 (2π) 3 2 |S(w)| 1 2 . (3.14)
The support size is:

s max = √ L 2 -ln(a) . (3.15)
As shown in [Azzalini & Valle 1996, Eq. (2.4)], parameters (w, α 1 , α 2 ) must ensure that the matrix below is positive definite:

  1 w-δ 1 δ 2 (1-δ 2 1 ) 1 2 (1-δ 2 2 ) 1 2 w-δ 1 δ 2 (1-δ 2 1 ) 1 2 (1-δ 2 2 ) 1 2 1   . (3.16)
Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification which amounts to satisfying the following technical conditions:

δ 1 δ 2 -(1 -δ 2 1 )(1 -δ 2 2 ) < w < δ 1 δ 2 + (1 -δ 2 1 )(1 -δ 2 2 ), α 1 = δ 1 -δ 2 w (1 -w 2 )(1 -w 2 -δ 2 1 -δ 2 2 + 2δ 1 δ 2 w) , (3.17) α 2 = δ 2 -δ 1 w (1 -w 2 )(1 -w 2 -δ 2 1 -δ 2 2 + 2δ 1 δ 2 w)
, for some (δ 1 , δ 2 ) ∈ (-1, 1) 2 . Hereagain, the maximal proportion a ∈ (0, 1) reached on the corner of the kernel support can be imposed by setting s 1 and s 2 no greater than s max given in (3.15) above. The skew-normal blur reduces to the Gaussian blur when α 1 = α 2 = 0. An example is displayed in Fig. 3.1(right).

Space-varying blurs modeling and inference

In this section, we first present our novel modeling approach for space-varying blur maps based on state-space models. We describe the standard BPF for inference in such models. After pointing the limitations of the BPF in this context, we propose the GIAnPF. Finally, we discuss the properties of GIAnPF and the connections with other methods of the literature.

State-space modeling for blur identification

We start by considering a generative SSM, where the hidden state represents the hidden (hence unknown) parameters of the variant PSFs. This generative modeling allows for a systematic Bayesian estimation of the unknown parameters, which are considered hidden states in the SSM literature [START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF]]. It allows to explicitly model the spatial smoothness of the PSF field. For example, the variation among neighbor patches can be a small rotation and change of width (see hereafter for a discussion regarding the order of the patches). Let us consider a state-space model given by p(ρ 0 ), . In particular, one can obtain the so-called filtering posterior distribution p(ρ t |X 1:t , y 1:t ), i.e., the posterior on ρ t conditioning on all data up to t.

(

Patch ordering for a sequential processing

The proposed approach requires setting the order in which the data will be sequentially processed. The parameters of consecutively processed patches must keep certain similarity in such a way the spatial information can be exploited. The underlying assumption in our approach is that the blurs affecting consecutive patches in the sequence have similar shape parameters. Otherwise stating, the parameters of the blurs should not change abruptly when going from patch t to patch t -1. Smooth PSF maps can then be identified, as soon as consecutive patches are spatially close. This is actually not the case if the patches are numbered naively following the lexicographic order, as illustrated in Figure 3.2(left). More suitable ordering must be adopted. If no additional structural assumption is available on the PSF map, we would suggest the zig-zag order, reminiscent from the one used in DCTbased image compression [Rabbani & Jones 1991] (Fig. 3.2(middle)). Circular blur maps, as those modeling phase aberrations ( [Denis et al. 2015, Sec.3.4], [Schuler et al. 2012, Simpkins & Stevenson 2014], [Mourya 2016, Chap.4]), would be identified better by using a spiral ordering (Fig. 3.2(right)). Phase aberrations (and thus, the PSF width) tend to increase when getting further from the center. The spiral ordering is thus more suitable as it implicitly promotes changes of the PSF parameters depending mostly on the distance to the image center. The spiral ordering will be retained in most experiments of the section. Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification

The smooth variations between consecutive kernels can be easily encoded using a Gaussian distribution for p(ρ t |ρ t-1 ). In order to ensure that the sampled blurs do not have degenerate shapes, range constraints can furthermore be enforced on some parameters (e.g., strict positivity for the width parameter, bounded PSF support), simply by considering truncated Gaussian distributions (see more details in Section 3.4). Note that it would be also possible to learn an optimal pattern to scan/process the data. Intuitively, the modeling/inference results of our proposed method would benefit for finding an order where consecutive kernels would be as similar as possible. This could be done for instance by first running a fast (but inaccurate) optimization-based strategy to choose the ordering in which the SSM would be constructed.

Bayesian inference in SSMs through particle filtering

Linear-Gaussian state-space models are often used since they allow for closed-form filtering and smoothing distributions via the well-known Kalman filter (KF) [Kalman 1960]. However, in this section we focus on more complex models which broadens the flexibility and applicability of our approach. Unfortunately, in these models the targeted pdf are intractable and approximations are required. Extensions of the KF exist, e.g., the extended Kalman filter (EKF) or the unscented Kalman filter (UKF) although they can greatly deviate from the true solution when the model is heavily non-linear non-Gaussian. In order to tackle these limitations, we consider an approach based on Monte Carlo approximations, and in particular on PF [Gordon et al. 1993]. We start by considering the BPF in Table 3.1, the first and arguably the most relevant PF [Gordon et al. 1993]. The BPF departs by simulating M samples (or particles) from the prior distribution. Then, for each patch t, it simulates M particles {ρ

(m) t } M
m=1 from the transition model by conditioning on the previous particles (Step 2a). Then, normalized importance weights are computed as the evaluation of the likelihood at each ρ (m) t in such a way all weights sum up to one (Step 2b). A resampling step is performed (Step 2d) in order to avoid particle degeneracy (see more details in [START_REF] Tiancheng | [END_REF]). As a result, the BPF approximates the filtering posterior distribution p(ρ t |X 1:t , y 1:t ) for each patch t from the set of M weighted particles {ρ

(m) t , w (m) t } M m=1 .
Note that other existing PFs can be also used in this problem, e.g., the APFs [Pitt & Shephard 2001, Elvira et al. 2018a, Elvira et al. 2019a] or the Rao-Blackwellized PF [START_REF] Djurić | [END_REF], Doucet & Johansen 2009[START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF]]. The research on PF has been very active in the last decade, and new PFs have been proposed in order to overcome existing challenges. In the specific application we tackle here, we are facing the well-known particle degeneracy problem, appearing due to very informative observations (hence, peaky likelihood). In particular, typical suitable parametric PSF models are very sensitive w.r.t. some of their parameters and the variance of the observation noise can be relatively low. Particle degeneracy refers to the (undesirable) effect where few particles (or even only one) remain after the resampling step. This reduces the diversity of the particle approximation and endangers the estimation of the blur parameters in the following patches. In the next section, we present a new PF that addresses the shortcomings of off-the-shelf PFs in such challenging

GIAnPF

In this section we propose the novel GIAnPF algorithm. It implements an annealing mechanism (see for instance [Neal 2001]) that overcomes the aforementioned limitations. Instead of directly approximating the targeted distribution, a sequence of modified distributions is considered. In such way, the algorithm starts by considering a more convenient distribution (i.e., sufficiently spread) while the last distribution is the true posterior of interest (in our case, the filtering distribution). The GIAnPF algorithm is described in Table 3.2. It starts by simulating M particles from the prior distribution in Step 1. However, unlike in the BPF or the APF, GIAnPF considers a sequence of Q intermediate steps in the processing of the data at each patch t. In the prediction step, M particles of the previous annealing layer {ρ (m) t,Q+1 } M m=1 are generated from the transition model by conditioning on the particles at time t -1 (Step 2a). Considering a sequence 0

= β Q+1 < β Q < • • • < β 1 < β 0 = 1
and starting from q = Q with a relatively small β Q , we first compute the weights by evaluating the tempered likelihood at ρ (m) t,q+1 (Step 2b), then each particle is resampled from the whole set of M particles with probability equal to the associated normalized weights (Step 2c). Finally, the new set of particles is simulated according to a transition kernel (Step 2d). We perform one step of Metropolis-Hastings (M-H) algorithm to sample from the transition kernel T t,q+1 (•) with proposal N (•, κ 2 q+1 I), where κ 2 q+1 is the variance of the kernel [Robert & Casella 2004], as it is common in sequential Monte Carlo samplers [Del Moral et al. 2006]. Note that in the last step, with q = 0, the exponent is β 0 = 1, and therefore the M-H targets the true posterior pdf of the kernel parameters, which ensures the invariance of the particle approximation. Rationale and parameter selection. The rationale of GIAnPF is as follows. The first exponent, β Q , is chosen in such a way the targeted distribution is a sufficiently flattened version of the likelihood. Then, for each annealing layer q = Q, . . . , 1, a larger β q-1 than the Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification

ρ (m) t,Q+1 ∼ p(ρ t,Q+1 |ρ (m) t-1,0 ), m = 1, ..., M (3.23) For q = Q, Q -1, . . . , 0, b)
Compute the normalized tempered weights as:

ξ (m) ∝ p(y t |ρ (m) t,q+1 , X t ) β q -β q+1 , m = 1, . . . , M (3.24) where 0 = β Q+1 < β Q < • • • < β 1 < β 0 = 1 c) resample ρ (m) t,q+1 = ρ ( j) t,q+1 , with probability ξ ( j) , for each m = 1, ..., M d) sample as ρ (m) t,q ∼ T t,q+1 (ρ (m) t,q+1 ), m = 1, . . . , M (3.25)
3. Output: At each t, we provide a set of M unweighted particles {ρ

(m) t,0 } M m=1
, for the approximation of the filtering distribution.

previous one is adopted so that the likelihood used when evaluating the weights becomes increasingly closer to the original one. The final particles at this time step t, {ρ (m) t,0 } M m=1 , are then those used for the approximation of the filtering distribution. This annealing procedure is especially helpful when the likelihood varies a lot among the particles, so that only few particles would have been chosen if the exact likelihood was used to do the sampling. In Step 2b), we compute the tempered weights to deal with the particle degeneracy problem. In Step 2c), the particle ρ (m) t,q+1 is resampled from the set of simulated particles, {ρ

( j) t,q+1 } M j=1
, with associated probabilities {ξ ( j) } M j=1 . In Step 2d), we apply one step of the M-H to the re-sampled particles, to ensure that the particle approximation converges to the true targeted distribution (i.e., without the tempering exponent). Note that GIAnPF can be seen as a generalized version of the BPF, since it contains this algorithm as particular case when we set Q = 0 and the M-H iteration (Step 2d) is avoided, i.e., with κ q+1 = 0. Connections to the literature. The GIAnPF connects with several algorithms in the literature. As stated above, it can be seen as a generalization of the BPF [Gordon et al. 1993]. It also holds clear links with annealing schemes such as [Neal 2001, Deutscher et al. 2000, Filippone 2014, Ulker et al. 2011]. In the context of PF, [Liu et al. 2021b] conducts interacting MCMC sampling procedure with the particles obtained from the BPF as the starting point. In [Madapura & Li 2007], the authors incorporate the annealing strategy within a Rao-Blackwellised Particle Filter and update the sample size using Kullback-Leibler Di-vergence transformation. The algorithm proposed in [START_REF] Gall | [END_REF]] also bears some similarities with GIAnPF in its use of an annealing scheme in a PF. However, unlike in [START_REF] Gall | [END_REF], our proposed method keeps closer ties to the standard BPF by modifying the sampling/weighting and adding an M-H update, which allows the particle approximation to converge to the filtering distribution with M (for any value of Q), and hence constructing consistent estimators of the sought intractable integrals.

Experimental results

Construction of an experimental database

We first describe the experimental database. It is composed of a set of several pairs of clean/degraded images associated to various choices for the patch decomposition, shift variant blur maps, and noise levels.

Images

We consider three images of size N = 512 × 512 (displayed in Fig. 3.4), namely, Chart, Cells, and Hubble. These are representative images of calibrated objects (e.g., resolution charts, fluorescent beads and distant stars) typically used to estimate spatially variant blur parameters. In all our experiments, patches of size 64 × 64 are used and the blurry noisy version y is generated following Eqs. (3.1)-(3.2), using i.i.d. zero-mean Gaussian noise with standard deviation σ n > 0, so that p(y t |ρ t ) = N (X t h(ρ t ), σ 2 n I).

Shift-variant blur maps

Five realistic space-varying blur maps, associated to the three blur shapes presented in Section 3.2.2, are constructed. To this aim, we take inspiration from the examples provided in [Escande & Weiss 2015, Fig.8], [Berisha & Nagy 2014, Fig.7], [Berisha & Nagy 2014, Fig.7] and [Simpkins & Stevenson 2014, Fig.4]. In most of our models, we will consider the existence of an optical center for the device, and that the blur width/orientation depends on the relative position of the patch with respect to this center. Instead of considering deterministic maps, as in the aforementioned works, we include some randomized settings in the expressions so as to simulate more realistic scenarios of an optical device whose effect on the image may slightly vary with respect to external conditions (temperature, laser power, planarity of the disposal). For a pre-specified abscissa-ordinate system in the 2D image plan, we denote, for a given patch t, o t = (o 1,t , o 2,t ) the coordinates of the patch center, ô = ( ô1 , ô2 ) the coordinates of the chosen optical center of the device, and t = o t -ô the Euclidian distance between both. The associated parameters, gathered in vector ρ t , are defined as follows, for the five considered models. Figure 3.3 displays examples of three realizations of each of those, illustrating the variability among patches, and also among different realizations for the same model.

Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification

Map 1 We assume Gaussian blurs, with size L = 15 2 , with width and orientation parameters s t and θ t at each patch t. We opt for the spiral patch ordering and we set, for every t,

θ t = arctan o 1,t -ô1 o 2,t -ô2 + ε t , (3.26)
where

ε t ∼ N ( π 2 , σ 2 ε ), s 1,t = δ 1,t t , s 2,t = δ 2,t t . (3.27) where δ 1,t ∼ U([δ 1,min , δ 1,max ]), δ 2,t ∼ U([δ 2,min , δ 2,max ]) with (σ ε , δ 1,min , δ 1,max , δ 2,min , δ 2,max ) positive scalars.
Map 2 We consider Gaussian blurs with size L = 15 2 with width and orientation parameters s t and θ t at each patch t. We use the spiral patch ordering and, for every t, we set

θ t = arctan o 1,t -ô1 o 2,t -ô2 + ε t , (3.28)
where

ε t ∼ N ( π 2 , σ 2 ε ), s 1,t = 1 8 (δ 1,t |o 1,t -ô1 | + 1) s 2,t = 1 8 (δ 2,t |o 2,t -ô2 | + 1) (3.29) with δ 1,t ∼ U([δ 1,min , δ 1,max ]), δ 2,t ∼ U([δ 2,min , δ 2,max ]
) and (σ ε , δ 1,min , δ 1,max , δ 2,min , δ 2,max ) some positive scalars.

Map 3 We consider Gaussian blurs with size L = 15 2 with width and orientation parameters s t and θ t at each patch t. The patches are ordered in lexicographic order and, for every t, we set

θ t ∼ U([0, π 8 ]), (3.30)
and

s 1,t = 1 8 (δ 1,t -o 1,t ), s 2,t = 1 15 (δ 2,t -o 2,t ), (3.31) where δ 1,t ∼ U([δ 1,min , δ 1,max ]), δ 2,t ∼ U([δ 2,min , δ 2,max ]
), and (δ 1,min , δ 1,max , δ 2,min , δ 2,max ) are some positive scalars.

Map 4 We assume blurs with defocus shape, with size L = 25 2 parametrized by width and orientation parameters s t and θ t at each patch t. We order the patches in lexicographic order and, for every t, we set

θ t ∼ U([0, π 8 ]), (3.32)
and

s 1,t = 1 3 (δ 1,t -o 1,t ), s 2,t = 1 6 (δ 2,t -o 2,t ), (3.33) with δ 1,t ∼ U([δ 1,min , δ 1,max ]), δ 2,t ∼ U([δ 2,min , δ 2,max ]
), and (δ 1,min , δ 1,max , δ 2,min , δ 2,max ) some positive scalars.

Map 5 We specify here a generative model for the particular case of the skew-normal blur, with size L = 15 2 , with the aim to mimic the shift-variant blur map that was considered in the experiments of [Simpkins & Stevenson 2014]. We order the patches in spiral order and we set

ω t = sign o 1,t -ô1 o 2,t -ô2 1 -2 |o 1,t -ô1 | -|o 2,t -ô2 | . (3.34)
We then create a shift-variant map with skew normal shapes of parameters

                 s 1,t = 1 3 (|o 2,t -ô2 | + 1 2 )( t + 1 5 ) + δ 1,t , s 2,t = 1 3 (|o 1,t -ô1 | + 1 2 )( t + 1 5 ) + δ 2,t , w t = 7 5 ω t t + δ 3,t , α 1,t = δ 4,t (1 -t ), α 2,t = δ 5,t (1 -t ), (3.35) with δ 1,t ∼ U([0, δ 1,max ]), δ 2,t ∼ U([0, δ 2,max ]), δ 3,t ∼ U([0, δ 3,max ]), δ 4,t ∼ U([δ 4,min , δ 4,max ]), δ 5,t ∼ U([δ 5,min , δ 5,max ])
, where (δ 1,max , δ 2,max , δ 3,max , δ 4,min , δ 4,max , δ 5,min , δ 5,max ) are positive scalars. Table 3.3 shows the settings for the realistic blur maps considered in Section 3.4.1.2. Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification

Quantitative comparison methodology

We will run several experiments of blur map estimation, using BPF and GIAnPF methods, both on toy examples and on our realistic dataset. Table 3.4 lists the transition models assumed when running BPF and GIAnPF, where N [s min ,s max ] (s t-1 , σ 2 s ) denotes a truncated normal distribution, i.e., the pdf is proportional to a normal pdf with mean s t-1 and variance σ 2 s in the support [s min , s max ], and 0 otherwise. More precisely, Model A will be Parameters settings

Model A      θ t ∼ N (θ t-1 , σ 2 θ ) s 1,t ∼ N [s min ,smax] (s t-1 , σ 2 s ) s 2,t ∼ N [s min ,smax] (s t-1 , σ 2 s ) s min = 10 -2 , s max set as in (3.5). Model B      θ t ∼ N (θ t-1 , σ 2 θ ) s 1,t ∼ N [s min ,smax] (s t-1 , σ 2 s ) s 2,t ∼ N [s min ,smax] (s t-1 , σ 2 s ) s min = 4 • 10 -2 , s max as in (3.10). Model C                s 1,t ∼ N [s min ,smax] (s t-1 , σ 2 s ) s 2,t ∼ N [s min ,smax] (s t-1 , σ 2 s ) w ∼ N [w min ,wmax] (s t-1 , σ 2 w ) α 1,t ∼ N [α min ,αmax] (s t-1 , σ 2 α ) α 2,t ∼ N [α min ,αmax] (s t-1 , σ 2 α )
s min = 5 • 10 -2 , s max as in (3.15), {w min , w max } = {-0.9, 0.9}, {α min , α max } = {0.5, 2}, so that (3.17) holds. employed to retrieve Gaussian blur maps (e.g., those resulting from Maps 1, 2 and 3), Model B to retrieve defocus blur maps (e.g., in case of Map 4) and Model C to estimate skew normal blur maps (e.g., Map 5). The upper bounds s max defined are used, to preserve consistency of the estimated shapes with the considered blur support width L. An analysis of robustness to the setting for L will be presented in Section 3.4.4.4. In order to quantify the gain of exploiting the spatial structure of the problem, we will compare the filtering pdf approximated by BPF and GIAnPF algorithms, and the posterior distribution that only considers the data corresponding to each patch independently. In the latter case, we use IS [Elvira et al. 2019b] to approximate the intractable posterior. Note that both PF and IS have a comparable computational complexity per sample/particle, while PF allows to exploit the information of previously processed data (i.e., patches), due to an implicit sequential IS structure. We also perform comparisons with two optimization-based methods considering non-parametric models for the kernels. In the so-called NP formulation, the kernels (h t ) T t=1 are estimated by minimizing a least-squares function under simplex and smoothness constraints, in a fully parallel manner for each patch t. We also compared with a more sophisticated non-parametric formulation where we included, in addition to the NP cost function, the total-variation based spatial regularization among kernels of neighboring patches from [Mourya 2016, Chap.4], yielding NP + method. Note that both resulting constrained convex minimization problems are solved with the FISTA algorithm [START_REF] Beck | [END_REF].

In all experiments, the relative mean squared error (RMSE) of blur kernel (RMSE h ) averaged over patches is used to evaluate the numerical performance of the blur estimation, i.e., RMSE

h = 1 T ∑ T t=1 ||h t -h t || 2
||h t || 2 with h t and h t the original kernels and the estimators re-spectively at patch t. For methods providing a posterior estimation (i.e., BPF, GIAnPF and IS), we calculate h t = h(ρ t ) with ρ t the mean estimator of the unknown parameters, while for NP and NP + , h t is directly the solution to the optimization problem, in patch t. We also include the standard deviation as well as the 95% credible interval of the estimators of BPF, GIAnPF and IS methods, to evaluate the statistical accuracy of the results. All the presented results are averaged over 100 random runs, and all the hyperparameters are tuned so as to minimize the RMSE (averaged over 5 random trials). The best results will be marked in bold cases in the tables. The numerical experiments are conducted in a Matlab environment on a computer with an Xeon(R) W-2135 processor (3.7 GHz clock frequency) and 12 GB of RAM.

Validation of the proposed method

We first discuss the settings and properties of the proposed method in an illustrative toy synthetic example. Our tests are conducted on the image Chart decomposed into T = 64 patches, and corrupted by skew-normal blurs whose parameters (ρ t ) T t=1 are generated following the transition model C, with (σ s , σ w , σ α ) = (10 -1 , 2 • 10 -1 , 10 -1 ) and, except otherwise stated, a noise standard deviation σ n = 0.05. The blur parameter estimations are then conducted with BPF and GIAnPF, assuming transition Model C with known hyperparameters.

Setting particle and layer numbers The BPF method, as well as the proposed GIAnPF method, requires the setting of the number of particles. For GIAnPF, the number of layers Q plays also a role. This is necessary to reach a compromise between a precise target reconstruction and a reasonable computational time, for setting those parameters. We display in Fig. 3.5 the averaged RMSE on 100 runs with respect to the computational time, for different settings of the parameters of BPF and GIAnPF. The associated values for the parameters (β q ) 0≤q≤Q are chosen on a logarithmic grid between 1 and 10 -2 , and β Q+1 = 0. For a given number of particles M, the GIAnPF algorithm always generates a better approximation of the unknown blurs than BPF, although at the expense of a slightly increased. Note that, for the same time budget, GIAnPF outperforms BPF. In all the upcoming experiments, we will set M = 3000 for BPF, IS, and Q = 2 and M = 500 for GIAnPF, so that they have comparable time requirement.

Alleviating the particle degeneracy One of the advantages of the novel GIAnPF method is the promotion of the diversity among the particles. In Fig. 3.6, we display the normalized effective sample size (NESS), computed as NESS =

1 M ∑ M m=1 w (m) t
, where w (m) t are the normalized weights (i.e., 0 ≤ NESS ≤ 1). We average the effective sample size (ESS) (see [START_REF] Kong | A note on importance sampling using standardized weights[END_REF], Kong et al. 1994, Elvira et al. 2018b]) over 100 runs and all the patches, divided by the corresponding number of particles M, for various level of noise σ n . The results are also averaged over 10 independent generations of the data. It can be seen that GIAnPF algorithm has higher NESS, thus preserves much more diversity in the particles for all analyzed noise levels. Similarly, Fig. 3.7 displays, for the same experiment, the RMSE of the Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification blur estimate using the mean of the approximate posterior of each filter. We can see that a larger diversity (larger NESS) translates into a smaller error.

Comparative performance on the experimental database

In this section, we compare the performance of GIAnPF and other approaches in different scenarios arising from our experimental database.

Quantitative performance Table 3.6 displays the RMSE h obtained when estimating kernels of Maps 1 to 5, from observations of Cells image, degraded by two noise levels, namely σ n = 0.01 and σ n = 0.05. Illustrative examples of results are also provided in Figure 3.8. In all the Bayesian methods, we display the blur maps resulting from kernel parameters equal to the MMSE estimator, i.e., the mean of the posterior distribution given the available data. The methods IS and NP, which do not exploit the spatial smoothness among neighboring kernels, reach the worst performance. A visual inspection of Figure 3.8 confirms this result. Among the spatially regularized methods, i.e., NP+, BPF and GIAnPF, the former is still far from reaching the quality of estimation of two latter methods though it improves over NP. Such behavior is expected as BPF and GIAnPF incorporate explicit knowledge about the parametric shape of the blur, leading to less diversity and thus less error in the restored blur maps. This can be seen in Figure 3.8. An analysis of the sensitivity of the methods over an error in the assumed parametric shape will be discussed in Section 3.4.4.1. In all cases presented in Table 3.6, the proposed method GIAnPF performs similarly or better than the standard BPF. The benefits from the annealing procedure are particularly noticeable in the case of σ n = 0.01. This improved performance can be explained by the annealed approach of GIAnPF, which is particularly effective when the posterior pdf of the parameters has the probability mass concentrated in small regions of the space, as it happens with low values of σ n (due to a peaky likelihood).

Complexity comparison

We perform a scalability analysis of the different methods, by applying those to the estimation of Map 5 blur map, from resized versions of the large scale aerial image from [START_REF] Chouzenoux | [END_REF]. The same hyperparameter settings as in Table 3.6 were used. In Table 3.5, we report the averaged computational time over 10 random trials. For IS and NP, that are fully separable methods onto patches, we also report the time per patch, that could be reached using parallel implementation on T cores. Both of these methods benefit from fast complexity cost assuming parallel processing of the patches. It is worthy to notice that the cost per patch of NP still increases for the largest images, probably due to memory saturation issues, while IS has a constant time requirement per patch. However, despite their rapidity, let us recall that the qualitative results for both these methods were rather poor on our experiments. The complexity burden dramatically increases for the regularized optimization-based method NP + . In the latter, the data is processed in a batch way. Parallelization over patches for this method is not straightforward, up to our knowledge, due to the non-separable structure of the underlying optimization problem. In contrast, the complexity of the proposed method stays very reasonable even for the largest image example, though involving the probabilistic estimation of T × K = 5120 parameters. As in the previous experiment, BPF and GIAnPF again show a similar computational cost.

Stability to the noise level We display in Fig. 3.9 a comparison of the RMSE h reached by the methods, when estimating kernels from observations degraded by increasing noise level σ n . One can observe that NP and NP+ performance rapidly deteriorates when the noise increases. In contrast, IS, BPF and GIAnPF perform in a relatively stable way. Moreover, the proposed GIAnPF method again outperforms the other competitors. As expected, its superiority over the standard BPF is less visible as σ n becomes higher, since GIAnPF is particularly dedicated to problems with peaky likelihoods, which arise when the noise level is low.

Robustness analysis

In practical contexts, it is frequent that the assumed parametric PSF model is erroneous (e.g., over-simplified with wrong support size) and/or that the calibrated image is only partially known (e.g., in digital camera imaging, a noisy non-blurry version can be acquired using short exposure settings). Furthermore, the noise level in the blurred image is usually estimated and not known with perfect accuracy. It is thus of main importance to quantify the robustness of our method to such model mismatch. To this aim, we conducted four experiments, aiming at reproducing different realistic scenarios of imperfect knowledge of the acquisition model. Except if specified otherwise, in all the examples of this section, we used the image Cells. The degraded image is obtained by applying the space-varying Gaussian blur Map 1, and then adding a noise with level equals to σ n = 0.05.

Choice of the parametric blur model

We first evaluated the influence of a mismatch between the assumed blur shape, and its actual one. We ran the different methods assuming instead a generalized Gaussian blur shape with different powers p, where p = 1 corresponds to the ground truth. The retrieved average RMSE h of BPF, GIAnPF, IS when taking different values for p are displayed in Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification Fig. 3.10. We also displayed the results of NP and NP + approaches, that remain unchanged since no parametric model for the kernels is considered in those methods. As expected, the best performance are obtained for the correct p = 1 setting. More interestingly, we can observe that, even if a wrong parametric model is adopted, i.e., using p = 1, our proposed method GIAnPF still performs better than its competitors, and in particular the non parametric method NP + , for a wide range of values for p.

Noisy calibrated image

In the previous experiments, we assumed that the original image x is known in a perfect manner. However, in reality, it is commonly not the case. Here, we consider the situation when only a noisy version of x is provided in the estimation process, denoted x = x + w where w is an additive Gaussian noise follows N (0, σ 2 x I). Eleven different levels for the noise on x are considered, namely σ x = {0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10}. Note that this range of values covers from low noise level values (where traditionally the particle-based filters struggle to operate due to the mass being concentrated in small regions; see for instance [Del Moral et al. 2006]) to relatively large noise values (which makes the inference less accurate). We display in Fig. 3.11 the averaged RMSE over 10 noisy realizations for each σ x , for the different approaches. The improvement of our proposed method GIAnPF with respect to BPF is more significant at low noise level, since particle degeneracy mostly arises in such context. Moreover, Fig. 3.11 shows that GIAnPF still outperforms its competitors when the noise level on x increases.

Setting of noise level

In the previous analysis, the exact observation noise variance σ n was assumed to be known and applied explicitly in our blur identification method. In more realistic situations, such noise level would have to be learned, leading to an estimate σn . We ran the different blur identification methods for the five different values σn = {0.01, 0.025, 0.05, 0.07, 0.1} (recall that the true observation noise level here is σ n = 0.05). Note that the optimizationbased methods are not sensitive to such change, so that the results for NP and NP + are unchanged. We display in Fig. 3.12 the average RMSE, for the different values for σn . We can notice that GIAnPF gives the best performance even when the noise level is poorly estimated. As expected, it reaches minimal mean and variance values when the ground truth noise level is used in the estimation process.

Setting of kernel size

In blur estimation, especially when the blur is non-stationary, a challenge lies in the setting a priori, of the kernel support width L. In this experiment we propose to assess the different approaches in the case where an erroneous L is assumed. In particular, we estimate the blur from the image Cells, when it is degraded by Map 1, and a noise level σ n = 0.05. We run the different methods/models, by setting the size for the kernel width as L ∈ {5, 7, 9, 11, 13, 15, 19, 23, 27}, while the groundtruth is L = 15. All the true/estimated kernels have been extended to the largest tested size L = 27 by zero-padding, to allow valid computation of estimation errors. We display in Fig. 3.13 the obtained RMSE values. One can notice that all methods have stable estimation error, as long as the assumed L is greater or equal to L, which is indeed expected. In contrast, when L < L, the results deteriorate for all methods. We observe that the GIAnPF method outperforms its competitors in almost all tested values for L.

Image restoration

We conclude this experimental section by illustrating the ability of the obtained kernel estimates to lead to satisfying and interpretable image restoration results. We assume that we have access to a pair (x, y) of clean/degraded version of a given calibrated image. We focus on the estimation of blurs resulting from our synthetic blurred maps, and then the use of those results for the restoration of several non-calibrated images. Through this example, we also illustrate and discuss the advantage of providing probabilistic blur estimates. We consider the image Chart, and its degraded version corrupted by the space-varying Map 4 with T = 64 and noise with standard deviation σ n = 0.1. We ran GIAnPF method and NP + approaches, identified as the two best blur identification procedures in the previous section. As shown in Fig. 3.14(left), we then perform the restoration of a set of four natural images2 , from their degraded version in Fig. 3.14(right), following the same blur model, and noise level σ x = 0, 0.01, 0.05 and 0.1. The restoration is performed in a non-blind manner, by making use of the fast majorize-minimize memory-gradient (3MG) algorithm from [Chouzenoux et al. 2013], that was recently adapted to the problem of image retrieval under shift-varying blur [Cadoni et al. 2016[START_REF] Chalvidal | [END_REF]]. Note that the variational formulation adopted in 3MG is very similar to the one in [Denis et al. 2015]. The regularization parameters of 3MG are finetuned so as to maximize SNR x = 20log 10 ||x|| ||x-x|| , with x the original image, and x the restored one. We present the results in Table 3.7, in terms of SNR x between x and the restored images obtained with the blur estimates of either GIAnPF or NP + . In the case of GIAnPF, we ran the restoration with 100 samples from the estimated posterior distribution of the blur maps, for a single run of GIAnPF, which allows us to compute the mean and the 95% confidence intervals for both image quality metrics. In Fig. 3.15(left), we show the mean of the pixel-wise squared error between the true kernel and these 100 samples obtained from the posterior distribution. This metric, called Bayesian mean squared error (MSE), allows to assess the uncertainty of GIAnPF on each kernel estimate. One can identify patches with larger estimation errors, corresponding to flat zones in the Chart images, thus making the blur estimation less accurate. In contrast, the patches localized in zones with significant content (sharp edges, in particular) benefit from very good estimation quality for the associated blur (pixel-wise error lower than 10 -6 ). Standard deviation is shown in Fig. 3.15(right), illustrating that the estimated distribution is more spread for difficult patches, while it is rather peaky for the patches that are better estimated. One can see from Table 3.7 that GIAnPF outperforms NP + in almost all examples, with up to 0.4 improvement in signal-to-noise 40 Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification ratio (SNR) score. Examples of visual results obtained from a given GIAnPF sample, are also displayed in Fig. 3.14(right).

We then illustrate the usefulness of the probabilistic estimation with associated uncertainty quantification that GIAnPF provides, for the image restoration task. Specifically, we focus on the restoration of Boat image degraded by Map 4 and no noise (i.e., σ x = 0). We display in Fig. 3.16(left) the uncertainty quantification, defined as the standard deviation per pixel of the 100 restored images, obtained using the 100 samples from the estimated posterior of the blur maps, in a single run of GIAnPF. This uncertainty map can be compared with Fig. 3.16(right). The latter displays the mean square error map, obtained by averaging the results of 100 restorations performed with the MMSE estimator of the kernel parameters (mean of the approximated posterior), obtained by 100 different independent runs of GIAnPF. It is noticeable that our probabilistic approach with the inference performed by GIAnPF; is able to quantify large uncertainty (see Fig. 3.16(left)) in areas where the squared errors are also large (see Fig. 3.16(right)). Let us point out that the uncertainty map from Fig. 3.16(left), is obtained without the need of processing the data multiple times, and does not require the knowledge of the original boat image.

Conclusion

This chapter addressed the estimation of the PSF parameters for spatially-varying blurs from calibrated image acquisitions. We proposed an original statistical modeling of the problem, accounting for the spatial dependency among neighboring kernels, and we applied a sequential Bayesian inference technique in this context. In order to alleviate the particle degeneracy problem brought by the BPF in some cases, we also proposed a new sampling method called the GIAnPF. Our results in different scenarios illustrated the good performance of the approach, including a useful uncertainty quantification. The novel approach opens many possibilities beyond this work. For instance, different noise distributions could be immediately used. Moreover, other state-space models, not necessarily Markovian, could be considered. These perspectives will be discussed thoroughly in Chapter 6. Each column displays one realization from the generative model. One can notice that, for each given row (i.e., blur map), the three images displayed share the same global aspect (e.g., circular map in rows 1, 2 or 5) with some slight variability, for instance in the blur kernel widths and orientations.

Chapter 3. Probabilistic modeling and inference for sequential space-varying blur identification

Introduction

Image blind deconvolution problem arises in many fields of image processing such as astronomy [Murtagh et al. 2007], biology [START_REF] Holmes | [END_REF]] and medical imaging [START_REF] Michailovich | [END_REF]. Given a degraded, blurred and noisy image, the aim is to restore a clean image along with an estimate of the blur kernel. Blind deconvolution is a severely ill-posed problem as there exists an infinite number of pairs (image/blur) that lead to the same observed image. Blind deconvolution methods available in literature adopt either a sequential identification process [Carasso 2001], or a joint estimation approach [Levin et al. 2009]. In the former, the blur kernel is identified first, possibly through a calibration step [Bell et al. 2016, Huang et al. 2021a, Chouzenoux et al. 2019]. Then the unknown image is inferred using a non-blind image restoration method. This approach was the one considered in our Chapter 3. In the latter, the blur kernel and unknown image are simultaneously estimated. This joint approach is retained in this chapter, and Chapter 5.

Since the problem of joint estimation is highly ill-posed, it is mandatory to incorporate prior knowledge on the sought unknowns. The retained prior strongly influences the choice for the solver. Let us distinguish two main classes of joint blind deconvolution approaches. A first option consists of formulating the problem as the minimization of a cost function gathering a data fidelity term (e.g., least-squares discrepancy) and penalties/constraints acting on the image and kernel variables. In such a way, it is quite standard to impose normalization and sparsity enhancing constraints on the kernel coefficients to avoid scale ambiguity inherent to the blind deconvolution model [START_REF] Komodakis | [END_REF], Kotera et al. 2013, Krishnan et al. 2011]. One can also easily impose the smoothness of the image, by adopting total-variation based regularization [START_REF] Chan | [END_REF]]. Several other efficient choices have been proposed in the literature, along with suitable iterative optimization methods to solve the resulting problems [Levin 2006, Shan et al. 2008, Joshi et al. 2008, Abboud et al. 2019, Bolte et al. 2010, Liu et al. 2021a]. The main advantage of this family of methods is probably its flexibility. But this comes at the price of heavy parameter tuning. The second option is to resort to a Bayesian formulation to express the model and a priori knowledge on the variables. The estimates are then defined from the estimation of the moments (typically, the mean) of a posterior distribution given the observed data and prior. As this typically involves the evaluation of intractable integrals, sampling [Robert & Casella 2004, Luengo et al. 2020] or approximation [Evans & Swartz 1995] strategies are used. MCMC methods have been widely used for blind deconvolution involving 1D sparse signals [START_REF] Rosec | [END_REF], Ge et al. 2011, Kail et al. 2012], but it is up to our knowledge scarcely employed in large-scale problems [Bishop et al. 2008], probably for computational time reasons. Another family of approach consists in adopting the so-called variational Bayesian approximation paradigm [START_REF] Fox | [END_REF], Blei et al. 2017, Zheng et al. 2015a, Gharsalli et al. 2012]. Then, a simpler (usually separable) approximation to the posterior is built through the minimization of a suitable divergence. This approach leads to fast Bayesian-based algorithms, whose great performance has been assessed in the context of non-blind [START_REF] Marnissi | [END_REF], Zheng et al. 2015b] and blind [Babacan et al. 2009] image restoration. Bayesian-based techniques usually require less parameters than optimization-based ones. Moreover, they can provide higher-order moments estimates, such as covariance matrices, which are of high interest for assessing probabilistically the uncertainty of the results. However, dealing with complex noise models and priors in such methods may be tricky, and the algorithms may be quite computationally heavy, hereagain with hyperparameters to be tuned.

A recent trend is to insert optimization-based steps in Bayesian sampling/approximation methods for a more versatile modality and faster computations. See, for example, [START_REF] Marnissi | [END_REF], Pereyra et al. 2016, Marnissi et al. 2020] for applications of such ideas in the context of large-scale image processing. Following this strategy, in this chapter, we introduce an algorithm called VBA, based on the Bayesian variational strategy (see Section 2.2.1), enhanced by optimization-based ideas inherited from [START_REF] Marnissi | [END_REF], with the advantages to cope with a large set of priors on the kernel and the image, and to present a reduced computational cost. The rest of Chapter 4 is organized as follows. In Section 4.2, we introduce the image degradation model and introduce our Bayesian modeling, and provide the background for deriving our algorithm. Section 4.3 explicitly describes the iterative updates of the proposed VBA in the case of Gaussian noise. Section 4.4 explicitly describes the iterative updates of the proposed VBA in the case of non-Gaussian noise. Section 4.5 concludes this chapter. The next Chapter 5 will be dedicated to the extension of VBA using the deep unrolling paradigm, and to the numerical experiments.

Problem statement 4.2.1 Observation model

We focus on the restoration of an original image x ∈ R N , from a degraded version of it y ∈ R N , related to x according to the following model:

(∀i ∈ {1, ..., N}) y i = D([ H x] i ).
(4.1)

Hereabove, n ∈ R N models some additive random perturbation on the observation. Moreover, H ∈ R N×N is a linear operator modeling the effect of a blur kernel h ∈ R M . D(•) is a noise model that may depend on the data available. In the following paragraphs, we'll discuss the case of Gaussian noise and non-Gaussian noise respectively.

In this work, we focus on the generalized blind deconvolution problem where the matrix associated with a given kernel h = [h 1 , . . . , h M ] reads

H = M ∑ m=1 h m S m , (4.2)
with {S 1 , . . . , S M } is a set of M known sparse N × N real-valued matrices. This model allows to retrieve the standard image deblurring model, as a special case when H identifies with a 2D discrete convolution matrix with suitable padding. The considered problem amounts to retrieving an estimate ( x, h) of the pair of variables ( x, h) given y. Due to the ill-posedness of this inverse problem, assumptions are required on the sought image / kernel and on the noise statistics to reach satisfying results. Furthermore, we introduce a linear equality constraint on the blur kernel estimate h. A general expression of such a constraint is as follows:

h = Tz + t, (4.3)
where T = (T m,p ) 1≤m≤M,1≤p≤P ∈ R M×P is a matrix of rank P ∈ {1, . . . , M} and t = [t 1 , . . . ,t M ] ∈ R M is some vector of predefined constants. Vector z = [z 1 , . . . , z P ] ∈ R P becomes the new unknown of the problem, along with the image x. A typical linear equality constraint in such context is the sum-to-one constraint, i.e., ∑ M m=1 h m = 1. Other examples will be provided in the experimental Section 5.2. We can thus rewrite (4.2) as

H = P ∑ p=1 z p K p + K 0 = H(z), (4.4) with (∀p ∈ {1, . . . , P}) K p = M ∑ m=1 T m,p S m ∈ R N×N , (4.5)
and

K 0 = M ∑ m=1 t m S m ∈ R N×N . (4.6)

Hierarchical Bayesian Modeling

Let us now introduce the hierarchical Bayesian model on which our VBA method will be grounded.

Likelihood

First, we express the likelihood p(y|x, z) of the observed data, given the unknowns (x, z).

In general case, the likelihood can be expressed as follows:

p(y|x, z) ∝ exp - N ∑ i=1 φ i ([Hx] i ; y i ) = exp - N ∑ i=1 φ i ([H(z)x] i ; y i ) . (4.7)
According to different noise models D(•), φ i can have various forms (see [Marnissi et al. 2017, Tab I] for examples). For instance, it reduces to a least square function in the case of additive Gaussian noise.

Prior

As already mentioned, it is necessary to incorporate suitable prior knowledge on the sought quantities to limit the problem ill-posedness. We here consider a wide range of sparsity enhancing prior for the image x, by adopting the generic model,

p(x|γ) ∝ γ N 2κ exp -γ J ∑ j=1 ||D j x|| 2κ , (4.8)
with κ ∈ (0, 1] a scale parameter and (D j ) 1≤ j≤J ∈ (R S×N ) J both assumed to be known. For instance, an isotropic total variation prior is obtained by setting κ = 1/2, S = 2, J = N and for every j ∈ {1, . . . , N},

D j x = [[∇ h x] j , [∇ v x] j ] ∈ R 2
gathers the horizontal and vertical gradients of x at pixel j. Other relevant choices are discussed in [START_REF] Marnissi | [END_REF].

Hereabove, γ > 0 is a regularization hyperparameter that we incorporate in our hierarchical model. We assume a Gamma distribution on γ, p(γ) ∝ γ α-1 exp(-ηγ), (4.9)

where α ≥ 0 and η ≥ 0 are the (known) shape and inverse scale parameters of the Gamma distribution. Such choice for the hyperparameter is rather standard in the context of Bayesian image restoration.

Regarding the blur h, we adopt the so-called SAR model, successfully used for Bayesianbased blind deconvolution in [Babacan et al. 2009]. The model relies on the following Gaussian model, (4.10) where A ∈ R Q×M with Q ∈ N \ {0} denotes a matrix of rank M. m ∈ R M is the mean of the underlying Gaussian distribution, and ξ > 0 is such that ξ A A is its inverse covariance matrix. If h follows this distribution, the projection of h onto the affine space defined by (4.3) is also Gaussian as well as the vector z associated with each projected vector. More precisely, z follows a Gaussian distibution with mean µ = T -1 (m -t) and covariance ma-

p(h|ξ ) ∝ ξ M 2 exp - ξ 2 ||A(h -m)|| 2 ,
trix ξ -1 T -1 (A A) -1 (T -1 )
where T -1 is the left inverse of T, i.e., T -1 = (T T) -1 T . This yields the following prior for the variable of interest z:

p(z|ξ ) ∝ ξ P 2 exp - ξ 2 (z -µ) L(z -µ) , (4.11)
where L = T T(T A A) -1 T -1 T T. We will consider (L, µ) to be predefined by the user, so as to be adapted to the sought properties of the blur kernel to estimate. The hyperparameter ξ will be learned during a training phase, as we will explain in Section 5.1.

Hierarchical model

Let us assume that (x, γ) and z are mutually independent. According to Bayes formula, the posterior distribution of the unknowns Θ = (x, z, γ) given the observed data y is defined as

p(Θ|y) ∝ p(y|x, z)p(x|γ)p(z|ξ )p(γ), (4.12)
where the four factors on the right side have been defined above.

Variational Bayesian Inference

The Bayesian inference paradigm seeks for solving the blind restoration problem through the exploration of the posterior p(Θ|y). Typically, one would be interested in the posterior mean, its covariance, or its modes (i.e., maxima). Let us make (4.12) explicit:

p(Θ|y) ∝ exp -γ J ∑ j=1 ||D j x|| 2κ - N ∑ i=1 φ i ([H(z)x] i ; y i ) × γ N 2κ +α-1 exp(-ηγ)ξ P 2 exp - ξ 2 (z -µ) L(z -µ) . (4.13)
Unfortunately, neither p(Θ|y), nor its moments (e.g., mean, covariance), nor its mode positions have a closed form. In particular p(y), which acts as a normalization constant, cannot be calculated analytically. We thus resort to the variational Bayesian framework, introduced in Section 2.2.1, to approximate this distribution by a more tractable one, denoted by q(Θ), for which the estimators are easier to compute. The approximation is computed with the aim to minimize the KL divergence between the target posterior and its approximation, which amounts to determining q opt (Θ) = argmin q KL(q(Θ)||p(Θ|y)),

= argmin q q(Θ)ln q(Θ) p(Θ|y) dΘ, (4.14)

where the equality holds only when q(Θ) = p(Θ|y). In order to make the solution of the above minimization problem tractable, a typical strategy is to make use of a so-called mean field approximation of the posterior, combined with an alternating minimization procedure, which yields to the so-called VBA. More specifically, we consider the following factorization: q(Θ) = q X (x)q Z (z)q Γ (γ).

(4.15)

Due to the intricate form of the chosen prior on the image and of the likelihood (in the case of non-Gaussian noise), we introduce extra approximation steps, relying on a majorationminimization (MM) strategy, reminiscent from [START_REF] Marnissi | [END_REF]]. In addition, we propose a strategy to reduce the time complexity of VBA, so as to deal with medium to large size images. In Section 4.3 and 4.4, we describe our developments in the case of Gaussian and non-Gaussian noises, respectively. As we will emphasize, even for the simplest case of Gaussian noise, the method requires the setting of two cumbersome hyperparameters, namely the regularization weight ξ and the noise level denoted β . In Chapter 5.1, we will show how to unroll the VBA method as a neural network structure, so as to learn the parameters (ξ , β ) in a supervised fashion. All the numerical experiments are reported in Chapter 5.1.

Case of Gaussian noise

Let us first focus on the resolution of the inverse problem (4.1) when D(•) models additive Gaussian noise. First, we reformulate the observation model as

y = H x + n. (4.16)
We then express the corresponding likelihood p(y|x, z) of the observed data, given the unknowns (x, z). Since the noise is assumed to be Gaussian distributed, the likelihood can be expressed as follows:

p(y|x, z) = β N 2 exp - β 2 ||y -H(z)x|| 2 , (4.17)
where we recall that β denotes the inverse of the noise variance. In the sequel, we will assume that n is a realization of an additive Gaussian noise with zero mean and standard deviation σ . In the remainder of the section, it will be convenient to set β = σ -2 . We now describe our proposed implementation of the VBA when applied to the approximation to the posterior in (4.13). We first present an MM-based procedure to handle the complicated form of the prior term on variable x. Then, we give the explicit expressions of the updates performed in the alternating minimization method.

MM-based approximation

Let us focus on the prior term in (4.8). This distribution is difficult to deal with as soon as κ is different from 1 (in which case a Gaussian distribution is retrieved). We thus propose

Case of Gaussian noise 57

to construct a surrogate for the prior on x. We use the tangent inequality for concave functions, which yields the following majorant function for the κ -function with κ ∈ (0, 1]:

(∀u > 0)(∀v ≥ 0) v κ ≤ (1 -κ)u κ + κu κ-1 v. (4.18)
Let us introduce the vector of auxiliary positive variables λ = (λ j ) 1≤ j≤J . From the previous inequality, we then deduce the following majorant function for the negative logarithm of the prior distribution:

(∀x ∈ R N ) γ J ∑ j=1 ||D j x|| 2κ ≤ J ∑ j=1 F j (D j x, λ j ; γ), (4.19)
where, for every j ∈ {1, . . . , J},

F j (D j x, λ j ; γ) = γ κ||D j x|| 2 + (1 -κ)λ j λ 1-κ j . (4.20)
This majorant function can be understood as a Gaussian lower bound on the prior distribution on x, which will appear more tractable in the VBA implementation. We will also show that the update of the auxiliary variables remains rather simple, thus not impacting the complexity of the whole procedure. In a nutshell, using (4.13), and (4.19), we obtain the following inequality:

p(Θ|y) ≥ F(Θ|y; λ) (4.21)
where the lower bound on the posterior distribution is

F(Θ|y; λ) = Cγ N 2κ exp - β 2 ||y -H(z)x|| 2 -F(x, λ; γ) p(γ)p(z|ξ ). (4.22)
Hereabove we have introduced the shorter notation

F(x, λ; γ) = J ∑ j=1 F j (D j x, λ j ; γ), (4.23)
and C is a multiplicative constant independent from Θ. Inequality (4.21) leads to the following majorization of the KL divergence involved in (4.14):

KL(q(Θ)||p(Θ|y)) ≤ KL(q(Θ)||F(Θ|y; λ)). (4.24)
By minimizing the upper bound in (4.24) with respect to λ, we can keep it as tight as possible, so as to guarantee the good performance of the VBA. To summarize, we propose to solve Problem (4.14) through the following four iterative steps:

1. Minimizing KL(q(Θ)||F(Θ|y; λ)) w.r.t. q X (x).

2. Minimizing the upper bound KL(q(Θ)||F(Θ|y; λ)) in (4.24) w.r.t. q Z (z).

3. Update the auxiliary variables (λ j ) 1≤ j≤J to minimize KL(q(Θ)||F(Θ|y; λ)).

4. Minimizing KL(q(Θ)||F(Θ|y; λ)) w.r.t. q Γ (γ).

Subsequently, at a given iteration k of the proposed algorithm, the corresponding estimated variables will be indexed by k.

VBA updates

Let us now describe the four steps of the proposed VBA, starting from a given iteration k associated with the current approximated distributions q k X (x), q k Z (z), and q k Γ (γ), and the auxiliary parameter estimate λ k . We also denote by (x k , z k , γ k ) the estimates of the means of q k X , q k Z , and q k Γ , and (C k x , C k z ) the covariance estimates for q k X and q k Z .

4.3.2.1 Update of q X (x)
By definition,

q k+1 X (x) = argmin q X KL(q X (x)q k Γ (γ)q k Z (z)||F(Θ|y; λ k )). (4.25)
The standard solution provided by (2.15) remains valid, by replacing the joint distribution by a lower bound chosen proportional to F(Θ|y; λ k ):

q k+1 x (x) ∝ exp < ln F(x, z, γ | y; λ k ) > q k Γ (γ),q k Z (z) ∝ exp ln F(x, z, γ | y; λ k )q k Γ (γ)q k Z (z)dγdz . (4.26)
By decomposing the different terms and using (4.4),

q k+1 x (x) ∝ exp - 1 2 x β E q k Z (z) (H) E q k Z (z) (H) + P ∑ p=1 P ∑ q=1 e p cov q k Z (z) (z)e q K p K q + 2E q k Γ (γ) (γ)D Λ k D x + β x E q k Z (z) (H) y (4.27)
where

E q k Z (z) (H) = P ∑ p=1 e p E q k Z (z) (z)K p + K 0 , (4.28)
D = [D 1 , . . . , D J ] , (4.29)
Λ k is the block diagonal matrix whose diagonal elements are (κ(λ k j ) κ-1 I S ) 1≤ j≤J , and (e 1 , . . . , e P ) is the canonical basis of R P . We thus obtain a Gaussian distribution:

q k+1 X (x) = N (x; xk+1 , Čx k+1 ), (4.30)
parametrized by

( Čx k+1 ) -1 = β (H k ) H k + P ∑ p=1 P ∑ q=1 e p C k z e q K p K q + 2γ k D T Λ k D, (4.31) xk+1 = β Čx k+1 (H k ) y, (4.32)
with

H k = H(z k ).
In image restoration applications, dimension N can be rather large (typically greater than 4.3. Case of Gaussian noise 59 10 6 variables), so that the storage of the full covariance matrix Čx k+1 is neither desirable nor usually possible. We thus propose to resort to a diagonal approximation to this matrix when required, so that the update finally reads:

q k+1 X (x) = N (x; x k+1 , C k+1 x ), (4.33)
with

C k+1 x = Diag δ k+1
x (4.34)

x k+1 = CG ( Čx k+1 ) -1 , β (H k ) y , (4.35)
where δ k+1

x ∈ R N is the vector of the inverses of the diagonal elements of ( Čx k+1 ) -1 , and CG(A,b) denotes the application of a linear conjugate gradient solver to the linear system Ax=b.

4.3.2.2 Update of q Z (z)
According to the VBA principle,

q k+1 Z (z) = argmin q Z KL(q k+1 X (x)q k Γ (γ)q Z (z)||F(Θ|y; λ k )). (4.36)
Using (2.15) and the previously introduced bound F(Θ | y; λ k ), we have

q k+1 Z (z) ∝ exp ln F(x, z, γ | y; λ k )q k Γ (γ)q k+1 X (x)dγdx . (4.37)
Replacing the involved quantities by their expression yields

q k+1 Z (z) ∝ exp - 1 2 z β B k+1 + ξ L z + z β a k+1 + ξ Lµ , (4.38)
where a k+1 = (a k+1 p ) 1≤p≤P ∈ R P and B k+1 = (B k+1 p,q ) 1≤p,q≤P ∈ R P×P are such that, for every (p, q) ∈ {1, . . . , P} 2 ,

a k+1 p =E q k+1 X (x) (x) K p y -E q k+1 X (x K p K 0 x) = (x k+1 ) K p y -B k+1 p,0 , (4.39)
B k+1 p,q = E q k+1 X (x K p K q x) = trace K p C k+1 x K q + (x k+1 ) K p K q x k+1 (4.40) with B k+1 p,0 = E q k+1 X (x K p K 0 x) = trace K p C k+1 x K 0 + (x k+1 ) K p K 0 x k+1 . (4.41)
Thus, the update for the distribution q Z reads

q k+1 Z (z) = N (z; z k+1 , C k+1 z ), (4.42)
with

(C k+1 z ) -1 = β B k+1 + ξ L, (4.43)
z k+1 = C k+1 z β a k+1 + ξ Lµ . (4.44)

Update of λ

Let us now express the update of the auxiliary variable. We aim at finding

λ k+1 = argmin λ KL(q k+1 X (x)q k Γ (γ)q k+1 Z (z)||F(Θ|y; λ)). (4.45)
This amounts to finding, for every j ∈ {1, . . . , J},

λ k+1 j = argmin λ j ∈[0,+∞) q k+1 X (x)q k Γ (γ)q k+1 Z (z) × log q k+1 X (x)q k Γ (γ)q k+1 Z (z) F(Θ|y, λ) dΘ, = argmin λ j ∈[0,+∞) J ∑ j=1 q k+1 X (x)q k Γ (γ) × F j (D j x, λ j ; γ)dxdγ, = argmin λ j ∈[0,+∞) κE q k+1 x (x) ||D j x|| 2 + (1 -κ)λ j λ 1-κ j . (4.46)
The explicit solution to the above minimization problem yields the following update:

λ k+1 j = E q k+1 x (x) ||D j x|| 2 = ||D j x k+1 || 2 + trace D j D j C k+1 x .
(4.47)

4.3.2.4 Update of q Γ (γ)
Finally, the update related to the hyperparameter γ is expressed as

q k+1 Γ (γ) = argmin q Γ KL(q k+1 X (x)q Γ (γ)q k+1 Z (z)||F(Θ|y; λ k+1 )). (4.48)
Using (2.15), we have

q k+1 Γ (γ) ∝ exp ln F(x, z, γ | y; λ k+1 ) × q k+1 X (x)q k+1 Z (z) dxdz . (4.49)
The above integral has the following closed form expression:

q k+1 Γ (γ) ∝ γ N 2κ +α-1 exp(-ηγ) × exp -γ J ∑ j=1 κ E q k+1 x (x) ||D j x|| 2 + (1 -κ)λ k+1 j (λ k+1 j ) 1-κ .
(4.50)

Case of non-Gaussian noise 61

It thus follows from (4.47) that the update of q Γ is

q k+1 Γ (γ) = Γ(d, b k+1 ), (4.51)
that is the Gamma distribution with parameters

d = N 2κ + α, b k+1 = J ∑ j=1 (λ k+1 j ) κ + η. (4.52)
The mean of q k+1 Γ is finally given by

γ k+1 = d b k+1 . (4.53)
Note that parameter d is not iteration dependent and can thus be precomputed from the beginning of the VBA.

Overview of VBA with Gaussian noise

Algorithm 1 provides a summary of the resulting VBA for solving the blind deconvolution problem introduced in Section 4.2 under Gaussian noise. We also specify our initialization strategy. More practical details about the latter will be discussed in the experimental Section 5.2. As a result, the optimal posterior distributions for both variables x and z will be approximated as Gaussian distributions, while the one for hyperparameter γ is approximated by a Gamma distribution. In particular, after K iterations, it is direct to extract from VBA outputs an estimate for the posterior mean of the image and the kernel, through variable x K and Tz K + t. The associated covariance matrices are given by C K x and TC K z T . These matrices can be useful to perform uncertainty quantification of the results. The VBA also allows us to estimate easily the hyperparameter γ involved in the image prior.

Algorithm 1 VBA approach for image blind deconvolution under Gaussian noise Initialization. Set hyperparameters (ξ , β , α, η). Define initial values for (x 0 , C 0

x , z 0 , C 0 z ). Compute λ 0 and γ 0 using (4.47) and (4.53), respectively. Iterative steps. For k = 0, 1, . . . , K:

1: Update the mean x k+1 and the covariance matrix C k+1

x of q k+1 X (x) using (4.34)-(4.35). 2: Update the mean z k+1 and the covariance matrix C k+1 z of q k+1 Z (z) using (4.43)-(4.44). 3: Update λ k+1 j using (4.47), for every j ∈ {1, . . . , J}.

4: Update the mean γ k+1 of q k+1 Γ (γ) using (4.52)-(4.53).

Case of non-Gaussian noise

In this section, we extend our implementation of the VBA when the noise degrading the image is non-Gaussian. We can still adopt the MM-based approximation for the prior term on variable x. The novelty here is to also apply an MM-based strategy to approximate the likelihood, which becomes more difficult for non-Gaussian noise. We detail below this new approximation, and give the explicit expressions of the updates performed in the resulting VBA.

MM-based approximation

Likelihood

Let us present our majorizing strategy to cope with the non-Gaussian noise case. Similarly as in [START_REF] Marnissi | [END_REF], Ass.III.1], we assume that the noise model leads to a likelihood in (4.7) such that, for every i ∈ {1, • • • , N}, φ i is differentiable on R and there exists µ i (y i ) > 0 such that the function defined by

v → v 2 2 -φ i (v;y i ) µ i (y i ) in convex on R. This assumption is satisfied when, for every i ∈ {1, • • • , N}, φ i is β i (y i )-Lipschitz differentiable.
Examples of useful functions and noise models satisfying this assumption were displayed in [Marnissi et al. 2017, Tab.1]. It is worthy to note that functions (φ i ) 1≤i≤N , and thus the parameters (µ i (y i )) 1≤i≤N implicitly account for the noise level present in the image.

For every i ∈ {1, • • • , N} and v ∈ R, let us define the following function:

ς i (v; y i ) = sup t∈R - 1 2 (v -t) 2 + φ i (t; y i ) µ i (y i ) . (4.54)
Then, according to [Marnissi et al. 2017, Prop.III.1], for every i ∈ {1, • • • , N},

(v ∈ R) φ i (v; y i ) = inf w i ∈R T i (v, w i ; y i ) (4.55)
where, for every v ∈ R,

T i (v, w i ; y i ) = µ i (y i ) 1 2 (v -w i ) 2 + ς i (w i ; y i ) . (4.56)
Moreover, the unique minimizer of

w i → T i (v, w i ; y i ) reads ŵi (v) = v - 1 µ i (y i ) φ i (v; y i ). (4.57)
This leads to construct an upper bound for the likelihood function.

(∀x ∈ R N ) φ (x, y) = inf w∈R N T (x, z, w; y) = inf w∈R N N ∑ i=1 T i ([H(z)x] i , w i ; y i ) (4.58) with (∀i ∈ {1, .., N}) T i (v, w i ; y i ) = µ i (y i ) 1 2 (v -w i ) 2 + ς i (w i ; y i ) (4.59)
and ς i defined in (4.54). Using (4.13), and (4.58), we obtain the following inequality:

p(Θ|y) ≥ F(Θ|y; w, λ) (4.60)
where the lower bound on the posterior distribution is

F(Θ|y; w, λ) = Cγ N 2κ exp (-T (x, z, w; y) -F(x, λ; γ)) p(γ)p(z|ξ ), (4.61)
using again the notation (4.23). Inequality (4.60) leads to the majorization KL(q(Θ)||p(Θ|y)) ≤ KL(q(Θ)||F(Θ|y; w, λ)).

(4.62)

By minimizing the upper bound in (4.62) with respect to λ and w, we can keep it as tight as possible, so as to guarantee the good performance of the VBA. As we aforementioned, the other steps of VBA remain the same as in the Gaussian case. To summarize, we propose to solve Problem (4.14) through the following five iterative steps:

1. Minimizing the upper bound KL(q(Θ)||F(Θ|y; w, λ)) w.r.t. q X (x).

2. Minimizing the upper bound KL(q(Θ)||F(Θ|y; w, λ)) w.r.t. q Z (z).

3. Update the auxiliary variables (w i ) 1≤i≤N to minimize KL(q(Θ)||F(Θ|y; w, λ)).

4. Update the auxiliary variables (λ j ) 1≤ j≤J to minimize KL(q(Θ)||F(Θ|y; w, λ)).

5. Minimizing the upper bound KL(q(Θ)||F(Θ|y; w, λ)) w.r.t. q Γ (γ).

At a given iteration k of the proposed algorithm, the corresponding estimated variables will be indexed by k.

VBA updates

Let us now describe the steps of the proposed VBA, starting from a given iteration k associated with the current approximated distributions q k X (x), q k Z (z), and q k Γ (γ), and the auxiliary parameters λ k and w k . We also denote by (x k , z k , γ k ) the estimates of the means of q k X , q k Z , and q k Γ , and (C k x , C k z ) the covariance estimates for q k X and q k Z .

4.4.2.1 Update of q X (x)

By definition,

q k+1 X (x) = argmin q X KL(q X (x)q k Γ (γ)q k Z (z)||F(Θ|y; w k , λ k )). (4.63)
Once again, we can use the solution provided by (2.15), by replacing the joint distribution by a lower bound chosen proportional to F(Θ|y; w k , λ k ):

q k+1 x (x) ∝ exp < ln F(x, z, γ | y; w k , λ k ) > q k Γ (γ),q k Z (z) ∝ exp ln F(x, z, γ | y; w k , λ k )q k Γ (γ)q k Z (z)dγdz . (4.64)
By decomposing the different terms and using (4.4),

q k+1 x (x) ∝ exp - 1 2 x E q k Z (z) (H) Diag(µ(y))E q k Z (z) (H) + P ∑ p=1 P ∑ q=1 e p cov q k Z (z) (z)e q K p Diag(µ(y))K q + 2E q k Γ (γ) (γ)D Λ k D x + x E q k Z (z) (H) u (4.65)
where

E q k Z (z) (H) = P ∑ p=1 e p E q k Z (z) (z)K p + K 0 , (4.66) D = [D 1 , . . . , D J ] , (4.67)
Λ k is the block diagonal matrix whose diagonal elements are (κ(λ k j ) κ-1 I S ) 1≤ j≤J , and (e 1 , . . . , e P ) is the canonical basis of R P . µ(y) = [µ 1 (y 1 ), ..., µ N (y N )] , u is a N × 1 vector whose i-th component is given by u i = µ i (y i )w k i . We thus obtain a Gaussian distribution:

q k+1 X (x) = N (x; xk+1 , Čx k+1 ), (4.68)
parametrized by

( Čx k+1 ) -1 = (H k ) Diag(µ(y))H k + P ∑ p=1 P ∑ q=1 e p C k z e q K p Diag(µ(y))K q + 2γ k D T Λ k D, (4.69) xk+1 = Čx k+1 (H k ) u, (4.70)
with

H k = H(z k ).
Hereagain, we propose to resort to a diagonal approximation to this matrix when required, so that the update finally reads:

q k+1 X (x) = N (x; x k+1 , C k+1 x ), (4.71)
with

C k+1 x = Diag δ k+1
x (4.72)

x k+1 = CG ( Čx k+1 ) -1 , (H k ) u , (4.73)
where δ k+1

x ∈ R N is the vector of the inverses of the diagonal elements of ( Čx k+1 ) -1 .

4.4.2.2 Update of q Z (z)
According to the VBA principle,

q k+1 Z (z) = argmin q Z KL(q k+1 X (x)q k Γ (γ)q Z (z)||F(Θ|y; w k , λ k )). (4.74)
Using the VBA optimality condition (2.15) and the previously introduced bound F(Θ | y; w k , λ k ), we have

q k+1 Z (z) ∝ exp ln F(x, z, γ | y; w k , λ k )q k Γ (γ)q k+1 X (x)dγdx . (4.75)
Replacing the involved quantities by their expression yields

q k+1 Z (z) ∝ exp - 1 2 z B k+1 + ξ L z + z a k+1 + ξ Lµ , (4.76)
where a k+1 = (a k+1 p ) 1≤p≤P ∈ R P and B k+1 = (B k+1 p,q ) 1≤p,q≤P ∈ R P×P are such that, for every (p, q) ∈ {1, . . . , P} 2 ,

a k+1 p =E q k+1 X (x) (x) K p u -E q k+1 X (x K p Diag(µ(y))K 0 x) = (x k+1 ) K p u -B k+1 p,0 , (4.77)
B k+1 p,q = E q k+1 X (x K p Diag(µ(y))K q x) = trace Diag(µ(y))K p C k+1 x K q + (x k+1 ) K p Diag(µ(y))K q x k+1 (4.78) with B k+1 p,0 = E q k+1 X (x K p Diag(µ(y))K 0 x) = trace Diag(µ(y))K p C k+1 x K 0 + (x k+1 ) K p Diag(µ(y))K 0 x k+1 . (4.79)
Thus, the update for the distribution q Z reads

q k+1 Z (z) = N (z; z k+1 , C k+1 z ), (4.80)
with

(C k+1 z ) -1 = B k+1 + ξ L, (4.81)
z k+1 = C k+1 z a k+1 + ξ Lµ . (4.82)

Update of w

Let us now express the update of the auxiliary variable w. We aim at finding

w k+1 = argmin w KL(q k+1 X (x)q k Γ (γ)q k+1 Z (z)||F(Θ|y; w, λ k )). (4.83)
which leads to

w k+1 = argmin w q k+1 X (x)q k Γ (γ)q k+1 Z (z) × ln q k+1 X (x)q k Γ (γ)q k+1 Z (z) L(Θ|y; w, λ k ) dΘ (4.84) = argmin w q k+1 X (x)q k Γ (γ)q k+1 Z (z) -lnL(Θ|y; w, λ k ) dΘ (4.85) = argmin w q k+1 X (x)q k+1 Z (z) N ∑ i=1 T i ([Hx] i , w i ; y i )dxdz (4.86) = argmin w q k+1 X (x)q k+1 Z (z) N ∑ i=1 T i ([H(z)x] i , w i ; y i )dxdz (4.87) = argmin w q k+1 X (x)q k+1 Z (z) N ∑ i=1 T i P ∑ p=1 z p K p + K 0 x i , w i ; y i dxdz. (4.88) (4.89)
It follows that, for every i ∈ {1, . . . , N},

w k+1 i = argmin w i T i H k+1 x k+1 i , w i ; y i (4.90) = H k+1 x k+1 i - 1 µ i (y i ) φ i H k+1 x k+1 i ; y i . (4.91)

Update of λ

Let us now express the update of the auxiliary variable λ. We aim at finding

λ k+1 = argmin λ KL(q k+1 X (x)q k Γ (γ)q k+1 Z (z)||F(Θ|y; w k+1 , λ)). (4.92)
This amounts to finding, for every j ∈ {1, . . . , J},

λ k+1 j = argmin λ j ∈[0,+∞) q k+1 X (x)q k Γ (γ)q k+1 Z (z) × log q k+1 X (x)q k Γ (γ)q k+1 Z (z) F(Θ|y, w k+1 , λ) dΘ, = argmin λ j ∈[0,+∞) J ∑ j=1 q k+1 X (x)q k Γ (γ) × F j (D j x, λ j ; γ)dxdγ, = argmin λ j ∈[0,+∞) κE q k+1 x (x) ||D j x|| 2 + (1 -κ)λ j λ 1-κ j . (4.93)
The explicit solution to the above minimization problem yields the following update:

λ k+1 j = E q k+1 x (x) ||D j x|| 2 = ||D j x k+1 || 2 + trace D j D j C k+1 x .
(4.94)

4.4.2.5 Update of q Γ (γ)
Finally, the update related to the hyperparameter γ is expressed as

q k+1 Γ (γ) = argmin q Γ KL(q k+1 X (x)q Γ (γ)q k+1 Z (z)||F(Θ|y; w k+1 , λ k+1 )). (4.95)
Using (2.15), we have

q k+1 Γ (γ) ∝ exp ln F(x, z, γ | y; w k+1 , λ k+1 ) × q k+1 X (x)q k+1 Z (z) dxdz . (4.96)
The above integral has the following closed form expression:

q k+1 Γ (γ) ∝ γ N 2κ +α-1 exp(-ηγ) × exp -γ J ∑ j=1 κ E q k+1 x (x) ||D j x|| 2 + (1 -κ)λ k+1 j (λ k+1 j ) 1-κ .
(4.97)

It thus follows from (4.94) that the update of q Γ is

q k+1 Γ (γ) = Γ(d, b k+1 ), (4.98)
that is the Gamma distribution with parameters

d = N 2κ + α, b k+1 = J ∑ j=1 (λ k+1 j ) κ + η. (4.99)
The mean of q k+1 Γ is finally given by

γ k+1 = d b k+1 .
(4.100)

Overview of VBA with non-Gaussian noise

Algorithm 2 provides a summary of the resulting VBA for solving the blind deconvolution problem introduced in Section 4 under non-Gaussian noise. As a result, the optimal posterior distributions for both variables x and z will be approximated as Gaussian distributions, while the one for hyperparameter γ is approximated by a Gamma distribution. We also provide the update formula for the hyperparameter w.

Algorithm 2 VBA approach for image blind deconvolution under non-Gaussian noise

Initialization. Set hyperparameters (ξ , α, η). Define initial values for (x 0 , C 0 x , z 0 , C 0 z ). Compute λ 0 , γ 0 and w 0 using (4.94), (4.100) and (4.91), respectively. Iterative steps. For k = 0, 1, . . . , K:

1: Update the mean x k+1 and the covariance matrix C k+1

x of q k+1 X (x) using (4.72)-(4.73). 2: Update the mean z k+1 and the covariance matrix C k+1 z of q k+1 Z (z) using (4.81)-(4.82). 3: Update w k+1 i using (4.91), for every i ∈ {1, . . . , N}. 4: Update λ k+1 j using (4.94), for every j ∈ {1, . . . , J}.

5: Update the mean γ k+1 of q k+1 Γ (γ) using (4.99)-(4.100).

Conclusion

In this chapter, we showed how to address the image blind deblurring problem through the variational Bayesian paradigm. We introduced a general hierarchical Bayesian model involving suitable priors for the image and the kernel. We introduced majorant functions to allow closed-form updates for each variables, for both Gaussian and non-Gaussian noises. Nonetheless, our preliminary experiments showed that it appears difficult to find an efficient manner to estimate the hyperparameter ξ , as this value highly fluctuates from one image/kernel pair to the other so that a simple prior modeling does not appear obvious. Moreover, in both considered noise settings, our VBA requires the knowledge of the noise level, through β (in Gaussian case) or (µ i (y i )) 1≤i≤N (in non-Gaussian case). This is limiting, and one might prefer to have this quantity estimated in an automatic manner. We thus propose in the next chapter, to resort to a supervised learning strategy to learn both ξ and β along the iterates of VBA, in the spirit of recent works [Bertocchi et al. 2019] on the unrolling (also called unfolding) of iterative algorithms.

CHAPTER 5

Deep unrolled variational Bayesian algorithm for image blind deconvolution In Chapter 4, we address the problem of blind image deconvolution using the variational Bayesian strategy. This strategy, as well as the references cited in the introduction of Chapter 4, belong to the class of traditional restoration methods. They rely on explicit and explainable steps, grounded on solid optimization/Bayesian methodological tools. This somewhat contrasts with the recent trend in image restoration relying on the deep neural network machinery. The use of deep learning for image deblurring only goes back to the last decade [Kupyn et al. 2018, Dong et al. 2016, Albluwi et al. 2018]. These methods rely on the deep learning methodology. More precisely, a supervised learning strategy is adopted to learn (implicitly) some prior information on the image/kernel from a so-called training set. A highly non-linear and multi-layers architecture is built, and its parameters (i.e., neuron weights) are estimated by back-propagation to minimize a given loss function associated with the task at hand (e.g., image visual quality). Several recent works propose neural network architectures dedicated to the problem of image blind deconvolution. Let us mention DeblurGAN [Kupyn et al. 2018], based on conditional generative adversarial networks and a multi-component loss function, SRCNN [START_REF] Dong | [END_REF]] and its extended version, DBSRCNN [Albluwi et al. 2018], relying on a convolutional neural network (CNN) architecture. SelfDeblur [Ren et al. 2020] combines an optimization-based Chapter 5. Deep unrolled variational Bayesian algorithm for image blind deconvolution method with two generative networks for modeling deep priors on image and blur kernel, respectively. Other image deblurring problems have been studied using neural network techniques, for example [Wen et al. 2021] presents an unsupervised method based on multi-adversarial CycleGAN [Zhu et al. 2017] for high-resolution image generation, and [Mohan et al. 2021] proposes a dynamic scene deblurring method for unconstrained duallens (DL) cameras. Most of these methods are supervised. They can perform quite well, as long as the training set is large and representative enough. Moreover, they are well suited to GPU-based implementation. However, they have traditionally suffered from lack of interpretability and robustness [START_REF] Bietti | [END_REF]].

Both traditional (optimization/Bayesian) approaches and deep learning methods are able to provide good performance. For most traditional image restoration methods, the common important challenge is the (hyper)parameter tuning. The parameters can arise from the model itself (e.g., prior weight, noise level) or from the resolution method (e.g., algorithm stepsize). Several strategies are available, based on empirical search, or more advanced statistical tools. However, these methods are often time-consuming without guaranteed success. While for deep learning methods, they lack theoretical guarantee and sometimes they demand numerous trials to find an optimal network structure. An emerging set of methods, in the field of inverse problems in signal/image processing, performs algorithm unrolling [Monga et al. 2021]. An iterative method (e.g., an optimization algorithm) is unrolled as layers of a neural network. The reduced set of parameters of this network are learnt by supervised training. Promising results have been obtained in the context of image deconvolution in [START_REF] Li | [END_REF], Bertocchi et al. 2019, Pesquet et al. 2021, Galinier et al. 2020], as well as other image restoration applications [Zhang et al. 2020]. One key advantage is that cumbersome parameter tuning can be avoided, by simply embedding the sought parameters into the network and learning them at training phase. Theoretical results assessing the stability and robustness of unrolling techniques can be found in [Pesquet et al. 2021[START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF], Chouzenoux et al. 2021]. These methods are also closely related to plug-and-play techniques where a trained neural network is employed as the denoiser [Pesquet et al. 2021, Sun et al. 2021, Zhang et al. 2021]. Deep unrolling can also be combined with a bilevel optimization framework, still with the aim for efficient hyperparameter learning (see [Crockett & Fessler 2022] for a recent survey).

In this chapter, we propose a novel approach for blind image deconvolution, that aims at gathering the best of the aforementioned methods. To be specific, we apply the unrolling paradigm to the algorithm presented in Chapter 4. We focus here on the case of Gaussian noise only. We create a deep neural network architecture, where VBA iterations proposed in Section 4.3 are integrated as layers. This allows us to (i) learn the hyperparameters (in particular, the noise level) of VBA in an automatic supervised fashion, (ii) improve further the quality of the results by choosing a dedicated loss in the training phase, and (iii) implement the method by taking full advantage of possible GPU resources, thus considerably reducing the processing time during the test phase. In contrast to standard deep learning methods for blind deconvolution, all these benefits come along with a preservation of the interpretability of the method, thanks to the unrolling technique. Let us emphasize that variational Bayesian methods often appear in deep learning context. Indeed, they are backbones of variational autoencoders [Kingma & Welling 2014] and also constitute methods of choice for training Bayesian neural networks [START_REF] Jospin | [END_REF]. However, up to our knowledge, our work is the first to investigate the unrolling of a variational Bayesian technique. The rest of Chapter 5 is organized as follows. The unrolling of VBA is presented in Section 5.1. Numerical results, including comparisons with various methods (in particular, our algorithm from Section 4.3), are presented in Section 5.2. Section 5.3 concludes this chapter.

Deep unrolled architecture

Overview

We introduce a supervised learning strategy to estimate the hyperparameter ξ and the inverse of the noise variance β , that are required to run VBA from Section 4.3. We adopt the so-called unrolling (or unfolding) methodology [Monga et al. 2021]. The idea is to view each iteration of an iterative algorithm as one layer of a neural network structure. Each layer can be parametrized by some quantities that are learned from a training database so as to minimize a task-oriented loss function. The advantage of the unrolling approach is threefold: (i) each layer mimics one iteration of the algorithm and thus it is highly interpretable, (ii) the choice of the loss is directly related to the task at the end, which is beneficial for the quality of the results, (iii) once trained, the network can be applied easily and rapidly on a large set of test data without any further tuning. In particular, its implementation can make use of GPU-accelerated frameworks.

Several recent examples in the field of image processing have shown the benefits of unrolling [START_REF] Gilton | [END_REF], Li et al. 2021, Tolooshams et al. 2021, Nan et al. 2020] when compared to standard black-box deep learning techniques or more classical restoration methods based on Bayesian or optimization tools. Let us in particular mention the works [START_REF] Li | [END_REF], Li et al. 2019] for the application of unrolling in the context of blind image restoration.

Let us now specify the unrolling procedure in the context of VBA. Let K > 0 be the number of iterations of the VBA described in Algorithm 1, thus corresponding to K layers of a neural network architecture. Iteration k ∈ {0, . . . , K -1} of our unrolled VBA can be conceptually expressed as

(x k+1 , C k+1 x , z k+1 , C k+1 z , λ k+1 , γ k+1 ) = A(x k , C k x , z k , C k z , λ k , γ k , ξ k , β k ).
(5.1)

The initialization procedure for (x 0 , C 0 x , z 0 , C 0 z , λ 0 , γ 0 ) is detailed in Algorithm 1. For k ∈ {0, . . . , K -1}, the expressions of (x k+1 , C k+1

x , z k+1 , C k+1 z , γ k+1 , λ k+1 ) as a function of the input arguments of A(•) are given respectively by (4.72)-(4.73), (4.81)-(4.82), (4.94), and (4.99)-(4.100). Furthermore, (ξ k , β k ) 0≤k≤K-1 are now learned, instead of being constant and preset by the user. This leads to the unfoldedVBA architecture depicted in Fig. 5.1, which can be summarized into the composition of K layers L K-1 • • • • • L 0 . Each layer L k with k ∈ {0, . . . , K -1} is made of three main blocks, that are two neural networks, namely NN k σ and NN k ξ , and the core VBA block A(•). There remains to specify our strategy for building the two inner networks, with the aim to learn (ξ k , β k ) 0≤k≤K-1 . deconvolution

ξ k ξ k+1 k σk n A σk+1 n A β k β k+1 y • • • • • • x K z K x V BA z V BA Layer k + 1 Layer NN k+1 σ NN k+1 ξ NN k σ NN k ξ C k+1 x C k x x k x k+1 z k z k+1 C k z C k+1 z γ k γ k+1 λ k λ k+1

Learning hyperparameter ξ

For every k ∈ {0, ..., K -1}, neural network NN k ξ takes as input the current kernel estimate h k = Tz k + t and delivers ξ k as an output. The architecture of the neural network is shown in Fig. 5.2. Note that the Softplus function, defined as

(∀x ∈ R) Softplus(x) = ln(1 + exp (x)), (5.2)
is used as a last layer, in order to enforce the strict positivity of the output hyperparameter ξ k .

Learning noise parameter β

When the noise parameter β is unknown, it might be useful to include a procedure to learn it automatically, again in a supervised fashion. In this case, we propose to introduce simple nonlinear mappings such that, for every k ∈ {0, . . . , K -1},

σ k = NN k σ (y) = Softplus(ρ k ) σ (y) + Softplus(τ k ), (5.3)
and β k = (σ k ) -2 . Hereabove, y is the observed degraded image, from which we deduce the wavelet-based variance estimator (also used in [Bertocchi et al. 2019]),

σ (y) = median(|W H y|) 0.6745 , (5.4)
where |W H y| gathers the absolute value of the diagonal coefficients of the first level Haar wavelet decomposition of the degraded image y. Moreover, (ρ k , τ k ) 0≤k≤K-1 are two scalar parameters to be learned during the training phase.

Complete architecture

We now present our complete blind deconvolution architecture for grayscale images and color images in Fig. 5.3. First, let us notice that VBA and its unrolled variant is designed for grayscale images. We thus generalized the architecture from Fig. 5.3(top), to process color images. To this end, we first transform the input RGB image to its YUV representation, which takes human perception into consideration. The network NN k σ is first applied to the luminance part y Y of the image. After applying the unfoldedVBA network (see Fig. 5.1), we obtain z VBA and x VBA as outputs. The latter is a restored version of the luminance channel. The remaining (U,V) color channels are simply obtained by median filtering of (y U , y V ). Both architectures in Fig. 5.3 additionally involve post-processing layers. More precisely, we first include a linear layer so as to encode the linear transformation (4.3), and then deduce the estimated blur kernel h. Second, we also allow a post-processing layer L pp acting on the image, so as to reduce possible residual artifacts, finally yielding x. In the case of color images, the post-processing is applied on the RGB representation to avoid chromatic artifacts.

Training procedure

The training of both proposed architectures from Fig. 5.3 requires to define a loss function, measuring the discrepancy between the output ( x, h) and the ground truth (x, h), that we denote hereafter by ( x, h, x, h). In the blind deconvolution application, one can for instance consider a loss function related to the error reconstruction on the kernel, or to the image quality, or a combination of both. Two training procedures will be distinguished and discussed in our experimental section, namely: Greedy training The parameters of the unfolded VBA are learned in a greedy fashion so as to minimize the kernel reconstruction error at each layer. Then, the post-processing network is learned in a second step, so as to maximize an image quality metric such as the structural similarity index measure (SSIM) [Wang et al. 2004]. End-to-end training The parameters of the complete architecture are learned end-to-end so as to maximize the image quality metric. Whatever the chosen training procedure, it is necessary to make use of a back-propagation step, that is to differentiate the loss function with respect to all the parameters of the network. Most operations involved in Fig. 5.3 can be differentiated efficiently using standard auto-differentiation tools. However, we observed in our experiments that it is beneficial (and sometimes even necessary) for a stable training phase to avoid using such tools deconvolution for differentiating the VBA layer A(•) involved in Fig. 5.1. In practice, we used the explicit expressions for the partial derivatives of it. Note that we followed the approach from [Charlier et al. 2021] to obtain the expression of the derivatives for the CG solver.

Experimental results

Problem formulation and settings

Problem overview

We focus on the resolution of the blind image deconvolution problem, where x ∈ R N is an original image, either grayscale or colored one. We come back to the model presented in Section 4.2.1, where the linear operator H ∈ R N×N models the application of a blur kernel h ∈ R M to the image. The noise n is assumed to be an additive white Gaussian noise with zero mean and standard deviation σ . In the case of color images, we assume that the same kernel, and the same noise level, is applied to each of the three RGB channels.

Datasets

Let us introduce the three datasets used to train and test our network, and compare our approach with state-of-the-art techniques. In all the cases, the training set is made of 100 images from the COCO training set. The validation set contains 40 images from the BSD500 validation set. For Dataset 1 and Dataset 2, the test set consists of 30 images from the Flickr30 test set. For Dataset 3, we adopted the same testing strategy as in [Ren et al. 2020], using the 4 images from [Levin et al. 2009]. Each image is center-cropped using a window of size N = 256 × 256. Each original image x is associated to a degraded version of it, y, built from Model (4.1). Various blur kernels and noise levels are used, as detailed hereafter. Dataset 1: All involved images are converted in grayscale. Each image of the database is blurred with 10 randomly generated Gaussian blurs, and then corrupted by additive noise. Thus in total, we have 1000 (= 100 × 10) training images, 400 (= 40 × 10) validation images, and 300 (= 30 × 10) test images for Dataset 1. The Gaussian blurs are of size 9 × 9. Two of them are isotropic with standard deviation randomly generated following a uniform law within [0.2, 0.4]. Eight of them are anisotropic with orientation either π/4 or 3π/4 (with equal probability) and vertical/horizontal widths (i.e., standard deviations of the 2D Gaussian shape) uniformly drawn within [0.15, 0.4]. On this dataset, the noise standard deviation is set to σ = 0.01, and assumed to be known (so that blocks (NN k σ ) 1≤k≤K of our architecture are overlooked). Dataset 2: All the images are then colored ones. We degraded each of them with 15 different blurs, namely 10 Gaussian blurs (simulated using the same procedure as above), two uniform blurs with width 5 × 5 and 7 × 7 pixels, and 3 out-of-focus blurs. For the latter, the vertical and horizontal widths are set randomly within [0.2, 0.5], and the orientation is either π/4 or 3π/4 (with equal probability). Furthermore, for each blurred image, zero-mean Gaussian noise is added, with standard deviation σ randomly chosen, with uniform distribution over [0.005, 0.05]. The noise level is not assumed to be known and is estimated using the proposed NN σ architecture. In total, we have 1500 (= 100 × 15) training images, 600 (= 40 × 15) validation images, and 450 (= 30 × 15) test images, on this dataset. Examples of blurs involved in Dataset 2 are depicted in Fig. 5.4. Dataset 3: All the images are converted in grayscale. We degraded each of them with 8 different real-world motion blurs from1 [Levin et al. 2009]. As the kernels have different sizes ranging from 13 × 13 to 27 × 27, we pad them with 0 values, so to reach a common large size of 27 × 27. For each blurred image, zero-mean Gaussian noise is added, with standard deviation σ randomly chosen, following a uniform distribution within [5 × 10 -6 , 5 × 10 -5 ]. The noise level is not assumed to be known and is estimated using the proposed NN σ architecture. In total, we have 800 (= 100 × 8) training images, 320 (= 40 × 8) validation images, and 32 (= 4 × 8) test images, on this dataset. The motion blurs involved in Dataset 3 are illustrated in Fig. 5.5.

VBA settings

Different settings are adopted for the implementation of the VBA step on each dataset, that we specified hereafter. For all considered datasets, we seek for kernels whose entries satisfy a sum-to-one constraint. For the datasets involving synthetic blurs (namely, Dataset 1 and Dataset 2), we impose additionally an axial symmetry along the main diagonal axis. In all cases, the constraints can be easily translated into the affine constraint (4.3). For Dataset deconvolution -1, while we have P = M -1 for Dataset 3. Regarding the choice of the blur prior, hereagain, a distinction is made between the first two scenarios and the last one. The synthetic blurs from Dataset 1 and Dataset 2 are either very smooth (e.g., Gaussian blurs) or present large constant areas (e.g., defocus blurs). In order to incorporate such prior knowledge, for both these datasets, A ∈ R (2M+1)×M is set to the matrix that computes the horizontal and vertical differences between pixels, augmented with an additional first row corresponding to an averaging operation, i.e., [1, . . . , 1]/M. This choice allows us to promote smooth variations in the kernels, while satisfying the required full column rank assumption on A. A constant vector with entries equal to 1 M is set for the prior mean m. In contrast, the motion blurs of Dataset 3 depicted in Fig. 5.5 present complex shapes, with no specific property on their gradient. We thus simply set A to the identity matrix in the prior, while a rough estimation of the blur obtained from2 [Kotera et al. 2013] is used for m. Matrices (D j ) 1≤ j≤J and parameter κ, involved in (4.8), are set in such a way that the chosen prior on the image yields an isotropic total-variation regularization (see our comment in Section 4.2.2.2). We must now specify the initialization for VBA iterates/layers. Our initial guess x 0 for the image is the degraded one. The associated covariance matrix C 0

x is initialized using the identity matrix. For Dataset 1 and Dataset 2, the blur is initialized with a 5 × 5 uniform kernel, from which we deduce the corresponding z 0 , while we use the prior mean for initializing the method for Dataset 3. In all VBA experiments, the covariance matrix C 0 z is set to a multiple of identity matrix. As the blurs involved in the Dataset 3 are of large size, a diagonal approximation of C k z is used in (4.69) so as to cope with memory saturation issues. A similar strategy was adopted for instance in [START_REF] Marnissi | [END_REF] and was not observed to yield any degradation of the resulting quality. The hyperparameters (α, η) involved in the prior law on parameter γ are set to zero in practice which is equivalent to impose a non-informative Jeffrey improper prior. Finally, the conjugate gradient solver used for the update of the image is run over 10 iterations which appears sufficient to reach practical stability. The solver is initialized with the degraded image.

Training specifications

We present results obtained by adopting the two training strategies described in Section 5.1.5. For the greedy training, we make use of the mean squared error on the estimated kernel, as a loss function for the unfoldedVBA layers, defined as ( x, h, x, h) = hh 2 . The SSIM loss ( [Wang et al. 2004]), between the output image x and the ground truth x is used to train the post-processing layer L pp . For the end-to-end training, we use again SSIM between x and x. We use warm initialization for end-to-end training, that is we initialize with the weights learnt during the greedy training phase, associated with a weight decay procedure. ADAM optimizer, with mini-batch size equal to 10, is used for all the training procedures. Its parameters such as learning rate (lr), weight decay (wd) and epochs number are finetuned, so as to obtain stable performance on each validation set. The number of layers K (i.e., number of VBA iterations) is set during the greedy training, and kept the same for the end-to-end training. In practice, we increase K as long as a significant decrease in the averaged MSE over the training set was observed. We specify in Table 5.1 all the retained settings. The train/validation/test phase are conducted with a code implemented in Pytorch (version 1.7.0) under Python (version 3.6.10) environment, and run on an Nvidia DGX wokstation using one Tesla V100 SXM2 GPU (1290 MHz frequency, 32GB of RAM). Our code is made available at https://github.com/yunshihuang/unfoldedVBA.

Comparison to other methods

The proposed method is compared to several blind deconvolution approaches available in the literature: Optimization-based methods: We first evaluate the VBA described in Section 4.3, in the favorable situation where the noise level σ is assumed to be known, and parameter ξ is finetuned empirically (see more details hereafter). VBA is run until reaching practical convergence, i.e., when the relative squared distance between two consecutive image iterates gets lower than 10 -5 . We also test two optimization-based approaches for blind deconvolution. The first one is called deconv2D. It makes use of the proximal alternating algorithm from [Bolte et al. 2010], to minimize a least-squares data fidelity term combined with various priors, namely total variation and positivity constraint on the image, sum-to-one and deconvolution quadratic constraint on the kernel. This method is implemented in Matlab, and inherits some of the software accelerations discussed in [Abboud et al. 2019] for blind video deconvolution. The second competitor in this category is the blinddeconv approach3 from [Krishnan et al. 2011]. For the sake of fair comparisons, for both datasets, we finetune the hyperparameters of these three methods on 40% of the training set and apply an average of the found values on the test set. Moreover, following the use of these three methods, we perform a non-blind deconvolution step BM3D-DEB4 [Lebrun 2012], which uses their respective estimated blur kernel to restore the image.

Deep learning methods: We perform comparisons with three recent deep learning architectures for blind deconvolution. SelfDeblur5 [Ren et al. 2020] is an unsupervised approach able to jointly perform the image restoration and kernel estimation tasks. DB-SRCNN6 [Albluwi et al. 2018] and DeblurGAN7 [Kupyn et al. 2018] are two supervised deep learning techniques. In contrast with SelfDeblur, they both only provide the estimated image, but do not estimate the kernel. Both these methods have been retrained on each of our datasets, using the same settings as in their initial implementation. Moreover, we adapted DBSRCNN to color images using the same pipeline as for our method, that is applying DBSRCNN on the luminance channel while simply nonlinearly filtering the chrominance ones. Finally, the kernels provided as outputs of the SelfDeblur method are normalized as a post-processing step to satisfy the sum-to-one constraint, to make the comparisons more faithful.

Evaluation metrics

All the methods are evaluated in terms of their performance on the blur kernel estimation (when available) and on the image restoration. Different metrics are used. For the blur kernels, we evaluate (i) the MSE, (ii) the so-called H ∞ error defined as the ∞ norm of the difference between the 2D discrete Fourier coefficients (with suitable padding) of the estimated and of the true kernel, and (iii) the mean absolute error (MAE) defined as the

Experimental results

Dataset 1

In Tables 5.2 and 5.3, we report the results of kernel estimation and image restoration, computed on the test set, using the different methods. As could be expected, the greedy approach tends to give more weight to the kernel quality than the end-to-end training. Our two training approaches yield great performance, when compared to all the other tested approaches. One can notice that the VBA with finetuned value for ξ performs quite well, showing the validity of our Bayesian formulation. The proposed unrolled VBA technique allows us to avoid a manual tuning of this parameter, and further increases the resulting quality. This is a direct outcome of the supervised training procedure aiming at maximizing quality scores, and also to the introduction of a post-processing step on the images. DBSRCNN has a good performance in terms of image restoration in this dataset. However, it is not capable of estimating the blur kernel, which might be useful for various applications. We display two examples of results in Fig 5 .7, extracted from our test set. One can notice, by visual inspection of these results, the high quality of the restored images. No artifacts can be observed, which is confirmed by a low average value of the PieAPP index on the test set. Moreover, the kernels are generally estimated quite accurately, as shown by the low MSE score and the good retrieval of their general structure. In the few cases when the unfolded VBA algorithm fails to give a perfect recovery of the blur kernel as in Fig. 5.7(bottom), the estimation is still accurate enough to yield a good recovery of the image whatever greedy training or end-to-end training is used. One can also notice that our method tends to provide better contrasted images, compared to its closest competitor in the image restoration task that is DBSRCNN. We display in Fig. 5.6(left) the evolution of the SSIM loss during the end-to-end training of the proposed architecture, showing the increase of the loss, then its stabilization, for both training and validation set, thus confirming an appropriate setting of ADAM optimizer parameters. Finally, Table 5.4(left) displays the average test time for each methods, that is the computational time required to restore one example of the dataset, once the method is finetuned/trained. We displayed CPU time for a fair comparison between methods, for codes ran on a Dell workstation equipped with an Xeon(R) W-2135 processor (3.7 GHz clock frequency and 12 GB of RAM). GPU time is also indicated when available. DBSRCNN makes use of an optimized C implementation, for its test phase on CPU. De-blurGAN is also very fast, but the resulting quality was quite poor in our experiments. Our method reaches a reasonable computational time on CPU. It becomes quite competitive when making use of GPU implementation, as the unrolled architecture is well suited for that purpose. This allows to drop the test time per image to few seconds, making it advantageous, with the addition benefit of better quality results in average, and of an available kernel estimate.

Dataset 2

The results of kernel estimation and image restoration on Dataset 2 using the various methods are shown in Tables 5.5 and 5.6, respectively. This dataset is more challenging, as it includes color images, various blur shapes, and various noise levels. The latter are not assumed to be known anymore. Hereagain, we can observe that the greedy training yields the best performance in terms of kernel estimation for the three considered metrics. In contrast, end-to-end training tends to favor the restored image quality while still providing a good kernel quality compared to other methods. In this more complicated context, standard VBA does not perform very well, as setting ξ becomes tedious for such an heterogeneous dataset. Let us note that the noise level is assumed to be known for this particular method, putting it in a quite favorable situation, compared to the other competitors, including our proposed approach. DBSRCNN provides again a good image recovery, but our proposed approach still outperforms it for both SSIM and PieAPP metrics. We display two examples of restoration in Fig 5 .8, when the sought blur is uniform, and out-of-focus, respectively. Such blur shapes are challenging and the MSE on the estimated blur might appear not excellent. Nevertheless, our method remains the best among the compared ones. The visual quality of the image generated by the proposed method is also very satisfying. We display in Fig. 5.6(right) the evolution of the SSIM loss during the end-to-end training, witnessing the absence of any overfitting issue. Moreover, we present in Fig 5.9 the evolution of the MSE loss on the kernel estimate, along the K = 21 layers of the architecture trained in an end-to-end manner. The MSE was averaged on test set examples associated to either Gaussian or out-of-focus blurs, respectively. These plots show that, for our choice of K (finetuned on the validation set), the MSE values are close to minimal. Larger K implied an increase of memory and training time, while not necessarily improving the results quality. One can also notice more fluctuations in the case of out-of-focus blur, which turns out to be more challenging to restore. A similar curve was obtained for uniform blurs, not shown by lack of space. Finally,

Dataset 3

The results of kernel estimation and image restoration on Dataset 3 using the various methods are shown in Tables 5.7 and 5.8, respectively. The level of complexity of the blind deconvolution task again increases compared to the former experiments, as the blurs are of larger size and have less smooth kernel. The noise level also has to be estimated. Regarding the performance of kernel estimation, we can observe that the best performance are reached by the optimization-based method blinddeconv, followed closely by our method. SelfDeblur and deconv2D are slightly behind, in terms of quantitative metrics. Regarding the image restoration, the results are more contrasted. Hereagain, the proposed method with end-to-end training gives the best performance regarding the quality of the restored image, and at the same time keeps the quality of blur kernel estimation compared to greedy training. The standard VBA procedure reaches rather fair quality on this dataset. We must howewer point out that, as for Dataset 2, this method was run under the simplifying assumption that the noise level is known. Both deep learning methods DBSRCNN and DeblurGAN fail to recover a good image quality. We display two examples of restoration in Fig 5 .10 where the sought kernels correspond to two different motion blurs from the dataset. We can notice that most of the methods struggle in the estimation of the blur kernel. Although SelfDeblur yields fairly good results according to a visual inspection, this method seems to suffer from a shift ambiguity, typical in blind deconvolution, which explains its rather limited quantitative scores compared to the other methods. As for the image quality, our proposed method is the most satisfying visually, which is aligned with its best PieAPP score by far among all the methods. We display in Fig. 5.11(left) the evolution of the SSIM loss during the end-to-end training. The evolution of the MSE loss on the kernel estimate is displayed in Fig. 5.11(right), across the K = 10 layers of the architecture trained in an end-to-end manner. Slight fluctuations can be observed compared to the plots obtained for both synthetic blur datasets, which shows the high challenge of estimating such complex motion blur kernels. Table 5.4(right) displays the average test time for each methods. With the exception of DBSRCNN and DeblurGAN (but yielding rather low quality results), all the methods see their computational time increasing, compared to the other two datasets. This is directly related to the increase of the blur kernel size. This being said, the computational time for our proposed method remains competitive, especially when making use of GPU.

Conclusion

This chapter proposed a novel method for blind image deconvolution that combines the variational Bayesian algorithm from Chapter 4 with a neural network architecture. Our experiments illustrated the excellent performance of this method on three datasets, comprising grayscale and color images, and degraded with various synthetic and real-world kernel types. Compared to state-of-the-art variational and deep learning approaches, our method delivered a more accurate estimation of both the image and the blur kernels. It also included an automatic noise estimation step, so that it required little hyperparameter tuning, com-Chapter 5. Deep unrolled variational Bayesian algorithm for image blind deconvolution pared to its counterpart from Chapter 4. The proposed method is very competitive in terms of computational time during the test phase, while showing similar train time to its deep learning competitors. The main core of the proposed architecture is highly interpretable, as it implements unrolled iterates of a well sounded Bayesian-based blind deconvolution method. As a byproduct, it also outputs estimates for the covariance matrices of both sought quantities (image/kernel). This information could be of interest for uncertainty quantification and model selection tasks (see for instance [Huang et al. 2021a, Repetti et al. 2019]). More generally, our work demonstrates that unrolling VBA algorithms constitute a promising research direction for solving challenging problems arising in Data Science. Image deblurring is a challenging problem common to several areas from astronomy, biology, medical science to embedded camera systems. It becomes even more severe when the blur is unknown and/or spatially-variant. Most state-of-the-art deblurring methods, including traditional optimization-based methods and more recent deep learning strategies, rely on deterministic models and tools, providing point-wise image estimates. In contrast, in this thesis, we focus on image restoration methods relying on the Bayesian framework. The image degradation model, the priors and the resolution tools become probabilistic. The estimates are defined from the estimation of the moments (typically, the mean) of a posterior distribution given the observed data and prior. Higher order moments (typically, the covariance) can also be exhibited, allowing to quantify uncertainty about the estimates. Bayesian approaches are thus very appealing. However, they typically involve the evaluation of intractable integrals. This can be overcome through the use of sampling or approximation strategies. Among available Bayesian techniques, few can efficiently tackle large scale image deblurring problems involving complex (e.g., spatially-variant) blurs and non-Gaussian noise. Moreover, hyperparameter tuning can be cumbersome. This PhD thesis has brought novel contributions to these challenges by proposing advances along two classes of Bayesian methods, namely particle filtering and variational Bayesian algorithms.

6.1 Contributions 6.1.1 Chapter 3: A new BPF approach for spatially-variant PSFs identification

We addressed the problem of spatially-variant PSFs identification from calibrated images. We formulated a state-space model to describe the space-varying PSF map, accounting for the spatial dependency among PSFs of the neighboring regions. We then proposed the algorithm GIAnPF, relying on the BPF paradigm, to solve the blur identification task in a probabilistic manner. Our algorithm implements an annealing mechanism to overcome the limitations of traditional BPF. The validity of the proposed method has been validated on numerous examples involving various shapes of blur kernels and PSF maps. This work let to the publications [Huang et al. 2019, Huang et al. 2021a].

6.1.2 Chapter 4: Novel VBA methods for blind image deconvolution

In [START_REF] Marnissi | [END_REF]], an efficient VBA strategy was proposed to tackle the case of non-blind image deconvolution. In this chapter, we proposed an extension of this methodology to solve the problem of blind image deconvolution. We introduced a generic linear equality constraints (e.g., sum-to-one) and a SAR prior on the unknown kernel. We then showed how to build suitable majorizing approximations of the original complicated posterior distribution, so as to yield a VBA method with simple updates while preserving its versatility. We provided two algorithms, respectively for addressing the case of Gaussian and non-Gaussian noise.

Chapter 5: Unrolled VBA for blind image deconvolution

Preliminary tests on the algorithms proposed in Chapter 4 showed high sensitivity to some hyperparameters. In particular, the parameter related to the kernel prior distribution is cumbersome to set. Moreover, the algorithms both require the knowledge of the noise level at hand. Finally, although reasonably good restoration results have been achieved, the methods were still far from the quality of some recent deep-learning methods. To address this challenge, we built an unrolling strategy where the iterations of VBA are recast as layers of neural network (NN) structures. The unrolling paradigm allowed to propose a supervised training approach where the NN parameters are learned so as to maximize image and kernel quality metrics. On top of solving the hyperparameter tuning problem, unrolling also led to fast algorithms amenable to implementations on GPU tools. Two unfolded VBA architectures have been proposed, able to process grayscale or color images, as well as to estimate the noise level. Our experiments on three datasets showed the superiority of our approach compared to several approaches from the state-of-the-art.

Perspectives

Our work raises several short and long term research paths that we list in this section.

Chapter 3

In this chapter, a BPF method has been proposed for spatially-variant blur identification. Several perspectives are foreseen, that we list herebelow, ranked from shorter term to longer term and open problems.

• Although our method is rather versatile to the noise and evolutionary laws assumed in our state-space model, we only tested it in scenarios involving Gaussian (or truncated Gaussian) laws. A first short-term perspective is thus to perform a more extensive numerical analysis, to cope with more sophisticated laws.

• Real-data experiments would certainly enforce the impact of this work. We did preliminary attempts, but we faced the issue of a mismatch between the piecewise constant spatially-variant blur degradation model considered in our approach, and the blur actually occurring in real-world data. Adapting our method to more realistic spatially-variant models (see [Denis et al. 2015, Escande & Weiss 2015]) would be an interesting path.

• Among realistic scenarios, we have considered the extension and application of our approach to the problem of blur calibration in radio interferometry [START_REF] Birdi | [END_REF]]. Discussions on this topic have been initiated this year with Dr. A. Repetti from Heriot-Watt University, Edinburgh. The main challenges are (i) the very large scale of the data, (ii) the complex-valued acquisitions in the Fourier domain, and (iii) the construction of a suitable SSM mimicking the neighboring interactions between the blur parameters.

• The robustness experiment considered in Section 3.4.4.2 is actually reminiscent from the image fusion problem [Tico & Pulli 2009]. In this inverse problem, one observes a pair of images, one being degraded with high noise level but no blur, the other being blurry and noisy with low level noise. The goal is to provide an estimation of the original image, jointly with the blur identification. Such problem arises typically in satellite imaging and in digital camera sensing, when pairs of short/long exposure images of the same scene can be acquired. Mapping such problem within our framework would require to include in our observation model the noise degradation of the image. Then, our idea would be to solve the resulting problem in an alternating fashion, using our GIAnPF for the blur identification step, while an optimization-based (e.g., total least squares) strategy would be employed to estimate, in a point-wise fashion, the original image.

• The main weakness of the model considered in Chapter 3 is that it is unidimensional. An important direction for future work is thus its extension to a multidimensional SSM to better encompass the variations among the PSFs. The increase of the dimensionality in the evolution equation of SSM is a challenging topic, as it requires to define suitable Markovian structures to preserve somewhat the causality relationships and thus to keep a sequential modeling/inference. A promising idea that we did not pursue by lack of time is to rely on 2D hidden Markov models as the ones described in [START_REF] Champagnat | [END_REF]]. The inference could be performed using tools described in [Lindsten & Schön 2013].

Chapters 4 and 5

• A first short-term goal is to extend the experiments from Chapter 5 to an example involving non-Gaussian noise, as it was done in the work by [START_REF] Marnissi | [END_REF]].

The VBA method in that case has been described in Chapter 4, but our numerical experiments on that algorithm stayed incomplete and were thus not reported in the manuscript. The extension of the unrolled architecture to the non-Gaussian case would also be an interesting avenue for short-term research.

• Another natural extension of these chapters is to tackle a real-world application. To that end, we are starting a collaboration with Dr. Nora Ouzir (CentraleSupélec), for extending and applying our method to the blind deconvolution of ultrasound imaging. This requires several changes in the noise and prior models, leading to a new structure for the unfolded VBA approach. We obtained promising preliminary results, and are confident that we will be able to finalize a publication on that topic by the end of the year.

• Our contributions in Chapter 4 would gain impact if they were assessed by theoretical convergence analysis results. Specifically, we would be interested in analyzing the convergence properties of VBA iterations, using recent results from the field of fixed point analysis [START_REF] Combettes | Fixed Point Strategies in Data Science[END_REF]. To conduct such study, we would go back to simpler settings, such as non-blind image deblurring under Gaussian noise. Despite the high interest of this topic, the convergence properties of VBA schemes have been scarcely studied in the literature, and we believe that the analysis tools might be now mature enough to start this ambitious topic.

• Deep unrolled architectures are promoted as being easier to study theoretically. In particular, stability to noise in the NN inputs/parameters can be studied [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF], Bertocchi et al. 2019]. Then, it would be very interesting to perform such analysis to the proposed architectures from Chapter 5.

Figure 1 :

 1 Figure 1: Noyaux originaux et estimés avec les méthodes d'optimisation NP, NP + , et les méthodes d'échantillonnage IS, BPF, GIAnPF. La métrique RMSE (Relative Mean Square Error) vaut, respectivement, 0.2117, 0.1756, 0.1258, 0.1246 et 0.0680. Ce résultat est extrait de la publication [Huang et al. 2021a].

Figure 2 :

 2 Figure 2: Architecture proposée pour la déconvolution aveugle, dans le cas d'images en niveaux de gris (haut) et couleur (bas).

  vector of n rows R m×n set of real matrices of m rows and n columns (a, b) real interval excluding a and b s scalars will be denoted by lowercase letters v vectors will be denoted by lowercase bold letters M matrices will be denoted by uppercase bold letters M transportation of the matrix M M -1 inverse of the matrix M trace(M) trace of the matrix M I m square identity matrix in R m×m v i i-th entry of the vector v M i, j element in i-th row and j-th column of the matrix M |s| absolute value of the scalar s |M| determinant of the matrix M ||v|| 2 norm of the vector v ||M|| 2 norm of the matrix M diag(M) vector composed of the diagonal entries of the matrix M Diag(v) diagonal matrix whose entries are given by the vector v δ f (v) δ v i functional derivatives of the function f with respect to v i ∂ f (v) ∂ v i partial derivatives of the function f with respect to v i ln natural logarithm exp natural exponential KL( f ||g) Kullback Leibler divergence between functions f and g x t the state of x at time step t x (i) t the i-th sample of the state of x at time step t p(x) prior distribution of x p(x|y) conditional distribution of x given yE x ( f (x))expectation of the function f (x) with respect to p(x)E x|y ( f (x)) expectation of the function f (x) with respect to p(x|y) U([a, b]) uniform distribution over the interval [a, b] N (µ, Σ)Gaussian distribution with mean vector µ and covariance matrix Σ N [a,b] (µ, Σ) truncated Gaussian distribution with mean vector µ and covariance matrix Σ in the support [a, b] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

  Figure 2.1: Examples of images degraded by different blurs.

  Figure 2.2: Examples of images degraded by different noises.

  Figure 2.3: Three examples of blurs that follow the diagonal symmetry.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.1 Related literature and contributions . . . . . . . . . . . . . . . . . 22 3.2.2 Relevant parametric PSF models . . . . . . . . . . . . . . . . . . . 23 3.3 Space-varying blurs modeling and inference . . . . . . . . . . . . . . . 25 3.3.1 State-space modeling for blur identification . . . . . . . . . . . . . 26 3.3.2 Patch ordering for a sequential processing . . . . . . . . . . . . . . 26 3.3.3 Bayesian inference in SSMs through particle filtering . . . . . . . . 27 3.3.4 GIAnPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.4.1 Construction of an experimental database . . . . . . . . . . . . . . 30 3.4.2 Validation of the proposed method . . . . . . . . . . . . . . . . . . 34 3.4.3 Comparative performance on the experimental database . . . . . . 35 3.4.4 Robustness analysis . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4.5 Image restoration . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

  3.8) 1 The lexicographic ordered vector associated to a matrix a b c d is the column vector [a b c d] .

Figure 3

 3 Figure 3.1: (left) Gaussian blur with θ = 0, s = [0.2; 0.1] √ L, α = 1, p = 1 and L = 15 2 ; (middle) defocus blur with θ = 0, s = [0.2; 0.1] √ L and L = 25 2 ; (right) skew-normal blur with s 1 = 0.15 √ L, s 2 = 0.15 √ L, w = 0, α 1 = 1.5, α 2 = 1.5 and L = 15 2 .

Figure 3

 3 Figure 3.2: Patch ordering (in the case of 4 × 4 structure): (left) lexicographic ; (middle) zig-zag ; (right) spiral.

  0.01, δ 1,min = 0.595, δ 1,max = 0.605,δ 2,min = 0.295, δ 2,max = 0.305 Map 2 σ ε = 0.01, δ 1,min = 2.95, δ 1,max = 3.05, δ 2,min = 2.95, δ 2,max = 3.05 Map 3 δ 1,min = 2.45, δ 1,max = 2.55, δ 2,min = 1.95, δ 2,max = 2.05 Map 4 δ 1,min = 1.45, δ 1,max = 1.55, δ 2,min = 1.45, δ 2,max = 1.55 Map 5δ 1,max = 0.01, δ 2,max = 0.01, δ 3,max = 0.1, δ 4,min = 1.5, δ 4,max = 2, δ 5,min = 2, δ 5,max = 2.5 Table3.3: Numerical settings for the generation of the experimental shift-variant blur maps.

Figure 3 . 3 :

 33 Figure 3.3: Three realizations (from left to right) of the space-varying blur Maps 1 to 5 (from top to bottom).

Figure 3 Figure 3 Figure 3 . 6 :Figure 3

 33363 Figure 3.4: Test images Chart, Cells and Hubble.

Figure 3 Figure 3

 33 Figure 3.9: Mean/variance of the RMSE h reached by the different methods. Hubble image blurred by space-varying kernels generated following Map 5 and three different values for σ n .

Figure 3 .Figure 3

 33 Figure 3.14: (left) Original images ; (middle) blurry noisy versions with σ n = 0.01 ; (right) example of restored image using the mean blur estimate from one run of GIAnPF.

Figure 3

 3 Figure 3.16: Uncertainty quantification vs error in the estimation of the kernels by GIAnPF on the image boat ; (left) log 10 of the pixel-wise standard deviation of the restored image, when using 100 samples from the posterior of the blur parameters, approximated by a single run of GIAnPF. (right) log 10 of the mean pixel-wise quadratic error between the true and the restored image when using the posterior mean estimator of the parameters approximated by GIAnPF, averaged over 100 independent runs.

Contents 5 . 1

 51 Deep unrolled architecture . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.1.2 Learning hyperparameter ξ . . . . . . . . . . . . . . . . . . . . . 71 5.1.3 Learning noise parameter β . . . . . . . . . . . . . . . . . . . . . 72 5.1.4 Complete architecture . . . . . . . . . . . . . . . . . . . . . . . . 72 5.1.5 Training procedure . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2.1 Problem formulation and settings . . . . . . . . . . . . . . . . . . 73 5.2.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

Figure 5 . 1 :

 51 Figure 5.1: Architecture of unfoldedVBA network.

Figure 5 . 2 :

 52 Figure 5.2: Neural network architecture NN k ξ for estimating ξ k , for k ∈ {0, ..., K -1}.

  Figure 5.3: Proposed blind image restoration pipeline for grayscale (top) and color (bottom) images.

Figure 5

 5 Figure 5.4: Examples of blur kernels used to construct Dataset 2.

Figure 5

 5 Figure 5.5: Motion blur kernels from Dataset 3.

  epoch = 10 K = 21, epoch = 10 K = 10, epoch = 10 lr = 5 × 10 -3 lr = 5 × 10 -3 (for L 0 ) lr = 10 -4 (for L 0 ) lr = 10 -3 (for other layers) lr = 10 -3 (for other layers) Post-processing L pp et al. 2015] epoch = 30, lr = 10 -3 epoch = 200, lr = 10 -3 epoch = 200, lr = 10 -= 5 × 10 -5 lr = 5 × 10 -5 lr = 5 × 10 -5 wd = 10 -4 wd = 10 -4 wd = 10 -4

Figure 5

 5 Figure 5.6: Evolution of SSIM loss along epochs of end-to-end training phase, averaged either on training or on validation sets of Dataset 1 (left) and Dataset 2 (right).

Figure 5 Figure 5

 55 Figure 5.9: MSE loss along the layers for proposed method using end-to-end training. Average over Dataset

  Bayesian inference methods . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Variational Bayesian algorithm . . . . . . . . . . . . . . . . . . . . 2.2.2 Importance sampling and adaptive importance sampling . . . . . . 2.3 Sequential Bayesian inference . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Probabilistic modeling and inference for sequential space-varying blur identification 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Related literature and contributions . . . . . . . . . . . . . . . . . 3.2.2 Relevant parametric PSF models . . . . . . . . . . . . . . . . . . . 3.3 Space-varying blurs modeling and inference . . . . . . . . . . . . . . . . . 3.3.1 State-space modeling for blur identification . . . . . . . . . . . . . 3.3.2 Patch ordering for a sequential processing . . . . . . . . . . . . . . 3.3.3 Bayesian inference in SSMs through particle filtering . . . . . . . . 3.3.4 GIAnPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Le problème de la restauration d'une image à partir d'une version dégradée, floutée et bruitée, de celle-ci se pose dans de nombreux cadres applicatifs, tels que la microscopie, l'imagerie médicale et l'astronomie. La restauration non aveugle s'apparente à l'inversion d'un système linéaire, lié à l'opérateur de flou, en présence de mesures bruitées. Une simple inversion matricielle semble fournir une approche de résolution triviale, mais dans la plupart des cas, elle implique de calculer l'inverse d'une matrice de très grande dimension, mal conditionnée voir non carrée, et quand il existe, le résultat de l'inversion directe est rarement satisfaisant. Dans le cas de la restauration aveugle, le problème est encore plus compliqué puisque l'opérateur de flou est lui aussi inconnu, et on observe que l'image dégradée à partir de laquelle l'image et le flou sont à estimer. De nombreuses approches ont été proposées pour la restauration d'image. Par exemple, les méthodes traditionnelles basées sur des approches d'optimisation, minimisant une fonction de coût appropriée par un algorithme d'optimisation itératif adapté. Ces approches sont fiables et montrent des résultats satisfaisant.

		Contents
	Résumé étendu	iii
	Abstract/Résumé	ix
	Notation	xi
	List of acronyms	xiii
	1 General introduction
	1.1 Motivation .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Background 2.1 Bayesian framework for imaging science . . . . . . . . . . . . . . . . . . . 2.1.1 Bayesian paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 MAP/MMSE estimator . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Image restoration . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 2.3.1 State-space models . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Particle filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 3.4.1 Construction of an experimental database . . . . . . . . . . . . . . 3.4.2 Validation of the proposed method . . . . . . . . . . . . . . . . . .

  Bayesian framework for imaging science . . . . . . . . . . . . . . . . . 7 2.1.1 Bayesian paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 MAP/MMSE estimator . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.3 Image restoration . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Bayesian inference methods . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Variational Bayesian algorithm . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Importance sampling and adaptive importance sampling . . . . . . 13 2.3 Sequential Bayesian inference . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 State-space models . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.2 Particle filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

	CHAPTER 2
	Background
	Contents
	2.1

  ∈ R n is the state at time step t and y t ∈ R m is the observation, t = 1, • • • , T ; p(θ 0 ) is the prior distribution of the hidden state θ 0 at the initial time step t = 0; p(θ t |θ t-1 ) is the state model, and p(y t |θ t ) is the observation model given the current state. The three equations above are in general parametric, and the model parameters are either known or estimated through point-wise estimators[Chouzenoux & Elvira 2020] or via probabilistic techniques[Andrieu et al. 2010].The sequence of states {θ 0 , θ 1 , • • • , θ T } is unknown, but can be estimated through the observations {y 1 , y 2 , • • • , y T }. Using the Bayes' rule, we can compute the full joint posterior distribution of all the states given all the observations as follows:

				2.25)
	θ t ∼ p(θ t |θ t-1 ),		(2.26)
	y t ∼ p(y t |θ t ),		(2.27)
	where θ t p(θ 0:T |y 1:T ) =	p(y 1:T |θ 0:T )p(θ 0:T ) p(y 1:T )	,	(2.28)

Table 2

 2

	.2: BPF algorithm

  = 1, ..., T , where p(ρ 0 ) is the prior distribution, p(ρ t |ρ t-1 ) is the transition model that generates the patch parameters t given the previous patch t -1, and p(y t |ρ t , X t ) is the observation model of patch t (that can be seen as the likelihood function when y t is observed). The goal is to approximate probabilistically and sequentially the unknown parameters ρ t for each patch t given subsets of data {X t , y t } T t=1

		3.18)
	p(ρ t |ρ t-1 ),	(3.19)
	p(y t |ρ t , X t ),	(3.20)
	for t

Table 3 .

 3 1: BPF algorithm for space-varying blur identification. 1. Initialization. Draw M i.i.d. samples, {ρ

	(m) 0 } M m=1 from the prior p(ρ 0 ).

scenario.

Table 3 .

 3 2: Generalized Interacting Annealed PF for space-varying blur identification.

	1. Initialization. Draw M i.i.d. samples, {ρ	(m) 0,0 } M m=1 from the prior p(ρ 0 ). Set Q ∈ N.
	2. Bayesian recursion. For t = 1, ..., T :
	a) Propagate the particles as

Table 3

 3

.4: Considered transitions models as in Eq. (3.19).

  Table3.5: Averaged computational time in seconds, for the different methods, in the case of Map 5 (i.e., L = 225 and K = 5) applied to various images of size N, decomposed into T patches. Time per patch is also reported for IS and NP, as they process patches independently.

	Image size	N = 256 2	N = 512 2	N = 1024 2	N = 2048 2
	Patches	T = 16	T = 64	T = 256	T = 1024
	NP	36.16 (2.26) 135.5 (2.11) 4987 (19.48) 82273 (80.34)
	NP +	35.80	127.3	4473	71913
	IS	5.21 (0.32) 21.85 (0.34) 91.29 (0.36)	374.7 (0.37)
	BPF	5.09	20.15	85.53	345.4
	GIAnPF	4.67	19.03	78.81	329.9

Table 3

 3

		0.45
		0.4
		0.35
		0.3
		0.25
	RMSE h	0.2
		0.15	NP
			NP	+
			IS
			BPF
			GIAnPF
		0	0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09	0.1
			x

Figure 3.11: Mean/variance of the RMSE h reached by the different methods, for various noise level σ x deteriorating the input calibrated image. h reached by the different methods, when using an erroneous noise level estimate σn (the ground truth is σ n = 0.05).

Figure 3.13: Mean/variance of the RMSE for various kernel support width L (ground truth size: L = 15). .7: SNR x values of restored images using estimated blur derived from NP + and GIAnPF applied on the calibrated image Chart.

Table 5 .

 5 1: Settings for the training phases in our experiments

  Chapter 5. Deep unrolled variational Bayesian algorithm for image blind deconvolution

	Method	MSE	H ∞ error	MAE
	VBA	0.0017 (0.0022)	0.1674 (0.1003)	0.0472 (0.0317)
	deconv2D	0.0025 (0.0035)	0.1483 (0.1037)	0.0489 (0.0395)
	blinddeconv	0.0013 (0.0011)	0.1553 (0.0660)	0.0417 (0.0203)
	SelfDeblur	0.0143 (0.0181)	0.3253 (0.1105)	0.1350 (0.0663)
	Proposed (greedy)	0.0008 (0.0012)	0.1165 (0.0677)	0.0281 (0.0168)
	Proposed (end-to-end)	0.0009 (0.0013)	0.1188 (0.0672)	0.0289 (0.0170)

The fastest method is DBSRCNN, though we must emphasize that, in contrast with all the other methods based on Matlab/Python softwares,

Table 5 .

 5 2: Quantitative assessment of the restored kernels. Mean (standard deviation) values computed over the test sets of Dataset 1.

	Method	SSIM	PSNR	PieAPP
	Blurred	0.6542 (0.1072)	22.2254 (2.3779)	4.1794 (0.9005)
	VBA	0.7603 (0.0814)	23.7332 (2.5672)	1.5109 (0.6184)
	deconv2D	0.7668 (0.0912)	24.5970 (2.8656)	1.9289 (0.4959)
	blinddeconv	0.7528 (0.0963)	23.9347 (2.4299)	1.9170 (0.6630)
	SelfDeblur	0.6948 (0.1006)	22.2704 (2.1255)	3.3178 (0.7291)
	DBSRCNN	0.7780 (0.0895)	24.9561 (2.9800)	1.5959 (0.6463)
	DeblurGAN	0.6613 (0.0731)	22.4388 (2.4074)	1.8937 (0.7630)
	Proposed (greedy)	0.7945 (0.0890)	24.7093 (2.9351)	1.4047 (0.6437)
	Proposed (end-to-end)	0.7989 (0.0886)	24.6638 (3.0711)	1.1976 (0.5433)

Table 5

 5

.3: Quantitative assessment of the restored images. Mean (standard deviation) values computed over the test sets of Dataset 1.

Table 5 .

 5 4: Average test time per image, using CPU (resp. GPU).

  Table 5.4(right) presents the average test time of the different methods. Again, our method appears competitive in terms of running time. Chapter 5. Deep unrolled variational Bayesian algorithm for image blind deconvolution

	Method	MSE	H ∞ error	MAE
	VBA	0.0148 (0.0139)	0.4492 (0.1638)	0.1339 (0.0627)
	deconv2D	0.0099 (0.0160)	0.2796 (0.1692)	0.0869 (0.0576)
	blinddeconv	0.0245 (0.0264)	0.3113 (0.1409)	0.1596 (0.1106)
	SelfDeblur	0.4336 (0.3188)	0.8720 (0.2774)	0.6975 (0.3047)
	Proposed (greedy)	0.0037 (0.0079)	0.1888 (0.1061)	0.0570 (0.0414)
	Proposed (end-to-end)	0.0039 (0.0079)	0.1960 (0.1056)	0.0588 (0.0411)

Table 5 .

 5 5: Quantitative assessment of the restored kernels. Mean (standard deviation) values computed over the test sets of Dataset 2.

	Method	SSIM	PSNR	PieAPP
	Blurred	0.5427 (0.1150)	21.7994 (2.1679)	4.2378 (0.8539)
	VBA	0.4024 (0.1571)	16.0371 (4.1798)	2.4218 (0.5545)
	deconv2D	0.6880 (0.1065)	23.1940 (2.8986)	2.2245 (0.6721)
	blinddeconv	0.6961 (0.1034)	23.2663 (2.7229)	2.3259 (0.8080)
	SelfDeblur	0.5107 (0.1305)	19.9943 (2.1467)	5.9269 (1.4066)
	DBSRCNN	0.6948 (0.1688)	23.6041 (4.2073)	1.9474 (0.7171)
	DeblurGAN	0.3370 (0.0740)	17.2781 (1.2909)	3.6581 (1.0040)
	Proposed (greedy)	0.7454 (0.1015)	23.2169 (2.4442)	1.7250 (0.5324)
	Proposed (end-to-end)	0.7518 (0.1025)	23.5631 (2.5959)	1.7681 (0.5502)

Table 5 .

 5 6: Quantitative assessment of the restored images. Mean (standard deviation) values computed over the test sets of Dataset 2.

Table 5 .

 5 7: Quantitative assessment of the restored kernels. Mean (standard deviation) values computed over the test sets of Dataset 3.

	Method	SSIM	PSNR	PieAPP
	Blurred	0.5830 (0.1500)	21.1404 (2.8626)	3.4647 (0.6608)
	VBA	0.7344 (0.1242)	24.3958 (3.0088)	1.7028 (0.4914)
	deconv2D	0.5440 (0.1422)	20.4229 (2.6241)	2.2811 (0.6407)
	blinddeconv	0.4847 (0.1669)	16.0690 (2.7140)	0.6957 (0.5241)
	SelfDeblur	0.7381 (0.1897)	24.3973 (4.6249)	0.6893 (0.3982)
	DBSRCNN	0.5327 (0.1111)	19.7017 (2.2976)	2.3438 (0.5672)
	DeblurGAN	0.5536 (0.1727)	21.0929 (3.2151)	2.3692 (0.7155)
	Proposed (greedy)	0.9197 (0.0353)	29.5730 (1.9817)	0.1778 (0.2408)
	Proposed (end-to-end)	0.9295 (0.0342)	31.0749 (2.0658)	0.2227 (0.2309)

Table 5 .

 5 8: Quantitative assessment of the restored images. Mean (standard deviation) values computed over the test sets of Dataset 3.

la propagation d'une source ponctuelle infinitésimale dans l'image acquise définissant la "Point Spread Fonction" (PSF). Lorsque la PSF est inconnue, on peut recourir à une stratégie de déconvolution aveugle pour restorer l'image. Une autre approche consiste d'abord à identifier le flou puis à l'éliminer de l'image dégradée par une stratégie de restauration non aveugle. Cette méthode demande une bonne précision pour l'estimation de la PSF, qui peut être atteinte efficacement par une étape préalable d'acquisition d'objets normalisés et calibrés, tels que des microbilles sphériques fluorescentes en microscopie ou des diagrammes de résolution pour l'étalonnage de l'appareil photo numérique. L'identification de la PSF est typiquement formulé sous la forme des moindres carrés. En plus de servir pour la phase de restauration d'image, la PSF identifiée peut également être utilisée pour ajuster un modèle paramétrique (souvent non linéaire) afin déterminer les caractéristiques d'un système optique. Dans la plupart des scénarios réalistes, un modèle de PSF stationnaire n'est pas adapté en raison de la profondeur de champ, du mouvement des objets ou de la caméra, d'abérrations optiques anisotropes de la lentille, ou de la turbulence atmosphérique. De telles sources de dégradation d'image donnent lieu à ce que l'on appelle flou variable spatialement. L'identification de la PSF dans ce contexte nécessite d'estimer une carte de plusieurs PSFs, décrivant le noyau de flou à chaque emplacement du plan optique. Dans ce chapitre, nous abordons ce problème de l'estimation bayésienne des paramètres d'une PSF variable spatialement à partir d'acquisitions d'images calibrées. Nous considérons un modèle séquentiel paramétrique pour la carte des PSFs qui suppose des faibles variations entre les PSFs agissant sur des régions voisines de l'image. Nous formulons un modèle d'espace d'états où chaque pas de temps correspond à un emplacement de patch dans l'image. Nous adoptons une approche probabiliste visant à produire la distribution a posteriori des paramètres

http://sipi.usc.edu/database/

https://webee.technion.ac.il/people/anat.levin/

http://zoi.utia.cas.cz/deconv_sparsegrad

Matlab code: https://dilipkay.wordpress.com/blind-deconvolution/

Matlab code:https://webpages.tuni.fi/foi/GCF-BM3D/index.html#ref_ software

Python/Pytorch code: https://github.com/csdwren/SelfDeblur

Python/Pytorch code: https://github.com/Fatma-ALbluwi/DBSRCNN

Python/Pytorch (training) and Matlab C-mex (testing) codes: https://github.com/ KupynOrest/DeblurGAN

norm of the difference between h and h. For evaluating the image quality, we compute (i) the SSIM, (ii) the peak-signal-to-noise ratio (PSNR), and (iii) the PieAPP value[Prashnani et al

. 2018], between the estimated image x and the ground truth x.

test examples involving either Gaussian (left), or out-of-focus (right) blur shapes.