<figure><img src="image_1.png" /></figure>

HAL
open science

# Deep learning for inverse problems and application to omni-directional imaging

Rita Fermanian

► To cite this version:

Rita Fermanian. Deep learning for inverse problems and application to omni-directional imaging. Signal and Image processing. Université de Rennes, 2023. English. NNT: 2023URENS112. tel-04608348

HAL Id: tel-04608348

https://theses.hal.science/tel-04608348

Submitted on 11 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
## THÈSE DE DOCTORAT DE

L'UNIVERSITÉ DE RENNES

ÉCOLE DOCTORALE N° 601

Mathématiques, Télécommunications, Informatique, Signal, Systèmes, Électronique

Spécialité : Signal, Image, Vision

Par
Rita FERMANIAN

## Deep learning for inverse problems and application to omnidirectional imaging

Thèse présentée et soutenue à Rennes, le 6 Decembre 2023

Unité de recherche : INRIA Rennes - Bretagne Atlantique

## Rapporteurs avant soutenance :

Andrés Almansa Directeur de Recherche, CNRS-Université Paris Cité

Said Moussaoui Professeur, École centrale de Nantes

## Composition du Jury :

Président :

Examinateurs :

Andrés Almansa Directeur de Recherche, CNRS - Université Paris Cité

Charles Kervrann Directeur de Recherche, INRIA Rennes - Bretagne Atlantique

Charles Yaacoub Professeur, Université catholique de Lille

Nicolas Papadakis Directeur de Recherche, CNRS - Institut de Mathématiques de Bordeaux

Said Moussaoui Professeur, École centrale de Nantes

Dir. de thèse : Christine Guillemot

Directrice de Recherche, INRIA Rennes - Bretagne Atlantique

Enc. de thèse : Mikael Le Pendu

Chercheur, Interdigital Rennes
<figure><img src="image_2.png" /></figure>
# TABLE OF CONTENTS

Résumé 1

1 Introduction 7

1.1 Context 7

1.2 Motivation 8

1.3 Contributions 10

1.4 Outline 11

I Inverse problems in perspective imagery 13

2 An overview of inverse problems in perspective imagery 15

2.1 Inverse Problems 15

2.1.1 Introduction and problem statement 15

2.1.2 Inverse problems in imaging: examples and definitions 17

2.1.3 An ill-posed problem and need for regularization 21

2.2 Reconstruction approaches for solving inverse problems 22

2.2.1 Variational approach for solving inverse problems 22

2.2.2 Learning approach for solving inverse problems: Neural network regression 25

2.2.3 Bayesian approach for solving inverse problems 26

2.3 Optimization algorithms 29

2.3.1 Derivative-based optimization algorithms: first-order methods 29

2.3.2 Proximal algorithms 31

2.4 Learning to Regularize inverse problems 34

2.4.1 Regularizing by using the implicit prior captured by a neural network 35

2.4.2 Regularizing by leveraging the power of denoisers 35

2.4.3 Regularizing through generative models 40

2.4.4 End-to-end learning of the regularization function via deep unrolling 41
3 Regularization of the Deep Image Prior with a deep denoiser 45

<table><tr><td>3.1 Introduction</td><td>45</td></tr><tr><td>3.2 Related work: regularized Deep Image Prior</td><td>46</td></tr><tr><td>3.3 Regularizing the Deep Image Prior with a learned denoiser</td><td>47</td></tr><tr><td>3.4 Experimental results</td><td>49</td></tr><tr><td>3.5 Conclusion and discussion</td><td>52</td></tr></table>

4 PnP-ReG: learned regularizing gradient for Plug-and-play gradient descent 55

<table><tr><td>4.1 Introduction</td><td>55</td></tr><tr><td>4.2 Notations and problem statement</td><td>56</td></tr><tr><td>4.3 Training of the gradient of a regularizer</td><td>58</td></tr><tr><td>4.3.1 Mathematical derivations</td><td>58</td></tr><tr><td>4.3.2 Training framework for the regularizing gradient network G</td><td>59</td></tr><tr><td>4.3.3 Training details</td><td>61</td></tr><tr><td>4.4 Experimental results</td><td>62</td></tr><tr><td>4.4.1 Plug-and-play gradient descent</td><td>62</td></tr><tr><td>4.4.2 Unrolled gradient descent with G</td><td>67</td></tr><tr><td>4.4.3 Analysis of the joint training</td><td>72</td></tr><tr><td>4.5 Discussion</td><td>75</td></tr><tr><td>4.6 Conclusion</td><td>76</td></tr></table>

II Inverse problems in omnidirectional imagery 81

5 A Comprehensive Review of omnidirectional images: processing methodologies and inverse problems 83

<table><tr><td>5.1 Introduction to omnidirectional imaging</td><td>83</td></tr><tr><td>5.2 Acquisition of omnidirectional images</td><td>85</td></tr><tr><td>5.3 Projecting spherical data to the Euclidean space</td><td>85</td></tr><tr><td>5.3.1 Equirectangular projection (ERP)</td><td>85</td></tr><tr><td>5.3.2 Cubemap projection (CMP)</td><td>86</td></tr><tr><td>5.4 Convolution on omnidirectional images</td><td>87</td></tr><tr><td>5.5 Inverse problems in omnidirectional images</td><td>89</td></tr><tr><td>5.5.1 Super-resolution</td><td>89</td></tr></table>
<table><tbody><tr><td>5.5.2</td><td>Image denoising</td><td>90</td></tr><tr><td><b>6</b></td><td><b>SphereDRUNet: A spherical denoiser for omnidirectional images</b></td><td><b>91</b></td></tr><tr><td>6.1</td><td>Introduction</td><td>91</td></tr><tr><td>6.2</td><td>Related work</td><td>92</td></tr><tr><td>6.2.1</td><td>HEALPix sampling</td><td>92</td></tr><tr><td>6.2.2</td><td>On-the-Sphere Learning for Omnidirectional Images (OSLO)</td><td>93</td></tr><tr><td>6.3</td><td>Spherical image denoising</td><td>94</td></tr><tr><td>6.3.1</td><td>Training details</td><td>96</td></tr><tr><td>6.3.2</td><td>Comparison framework</td><td>97</td></tr><tr><td>6.4</td><td>Experimental results</td><td>98</td></tr><tr><td>6.4.1</td><td>Spherical convolution: OSLO-based vs. graph-based (DeepSphere)</td><td>99</td></tr><tr><td>6.4.2</td><td>Ablation study: training a denoiser on the equirectangular projection</td><td>101</td></tr><tr><td>6.5</td><td>Conclusion</td><td>103</td></tr><tr><td><b>7</b></td><td><b>Conclusion</b></td><td><b>105</b></td></tr><tr><td><b>A</b></td><td><b>Parameter settings of Chapter 4</b></td><td><b>109</b></td></tr><tr><td>A.1</td><td>Plug-and-play ADMM: formulation and parameter settings</td><td>109</td></tr><tr><td>A.2</td><td>Parameter settings for the other methods</td><td>111</td></tr><tr><td><b>Bibliography</b></td><td></td><td><b>113</b></td></tr></tbody></table>
LIST OF FIGURES

<table><tr><td>2.1</td><td>Degradation of an image by different types of noise.</td><td>17</td></tr><tr><td>2.2</td><td>PDF of Gaussian noise.</td><td>18</td></tr><tr><td>2.3</td><td>Inverse problem of denoising.</td><td>19</td></tr><tr><td>2.4</td><td>Inverse problem of deblurring.</td><td>20</td></tr><tr><td>2.5</td><td>Inverse problem of super-resolution.</td><td>20</td></tr><tr><td>2.6</td><td>Inverse problem of pixel-wise inpainting.</td><td>21</td></tr><tr><td>2.7</td><td>Unrolled gradient descent framework.</td><td>42</td></tr><tr><td>3.1</td><td>Architecture of the Deep Image Prior.</td><td>46</td></tr><tr><td>3.2</td><td>Set of images used for experiments in Chapter 3.</td><td>49</td></tr><tr><td>3.3</td><td>Visual comparison of super-resolution results for the DIP regularization methods</td><td>53</td></tr><tr><td>3.4</td><td>Visual comparison of super-resolution results for the DIP regularization methods</td><td>53</td></tr><tr><td>4.1</td><td>Framework for joint training of D and G.</td><td>60</td></tr><tr><td>4.2</td><td>Architecture of G.</td><td>61</td></tr><tr><td>4.3</td><td>Visual comparison of super-resolution results for ReG.</td><td>66</td></tr><tr><td>4.4</td><td>Visual comparison of deblurring results for ReG.</td><td>70</td></tr><tr><td>4.5</td><td>Visual comparison of pixel-wise inpainting results for ReG.</td><td>71</td></tr><tr><td>4.6</td><td>Visual comparison of super-resolution results for ReG in unrolled-GD.</td><td>72</td></tr><tr><td>4.7</td><td>Visual comparison of deblurring results for ReG in unrolled-GD.</td><td>73</td></tr><tr><td>4.8</td><td>Comparison of ReG's performance under different training strategies of ReG.</td><td>77</td></tr><tr><td>4.9</td><td>Comparison of the performances of the original and the updated denoisers in PnP-ADMM.</td><td>78</td></tr><tr><td>4.10</td><td>Comparison of the convergence of PnP-ADMM, RED, GS-PnP and PnP-ReG for super-resolution.</td><td>78</td></tr><tr><td>4.11</td><td>Comparison of the convergence of PnP-ADMM, RED, GS-PnP and PnP-ReG for deblurring.</td><td>79</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>5.1</td>
      <td>Image processing pipeline for 2D images</td>
      <td>84</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Image processing pipeline for omnidirectional imagery</td>
      <td>84</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Equirectangular projection</td>
      <td>86</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Cubemap projection</td>
      <td>87</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Visualization of a fish-eye image and its different planar projections</td>
      <td>88</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Visualization of HEALPix sampling over the sphere.</td>
      <td>93</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Spherical convolution in the pixel domain proposed by OSLO.</td>
      <td>94</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Architecture of the DRUNet network.</td>
      <td>95</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Architecture of the SphereDRUNet network.</td>
      <td>96</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Visual comparisons of denoising the sphere with SphereDRUNet and its projection with a DRUNet that was re-trained on ERP images.</td>
      <td>99</td>
    </tr>
    <tr>
      <td>6.6</td>
      <td>Visual comparisons of denoising the sphere with SphereDRUNet and its projection with a DRUNet that was re-trained on ERP images.</td>
      <td>100</td>
    </tr>
    <tr>
      <td>6.7</td>
      <td>Visual comparisons of denoising the sphere with SphereDRUNet and its projection with a DRUNet that was re-trained on ERP images.</td>
      <td>101</td>
    </tr>
    <tr>
      <td>6.8</td>
      <td>Illustration of the distortion caused by cubemap projection</td>
      <td>102</td>
    </tr>
    <tr>
      <td>6.9</td>
      <td>Comparison between OSLO-based SphereDRUNet and graph-based Sphere-DRUNet.</td>
      <td>103</td>
    </tr>
  </tbody>
</table>
LIST OF TABLES

<table><tbody><tr><td>3.1</td><td>Parameters of the proposed DIP-denoiser-ADMM</td><td>50</td></tr><tr><td>3.2</td><td>Results of denoising for the DIP regularization methods</td><td>51</td></tr><tr><td>3.3</td><td>Results of super-resolution for the DIP regularization methods</td><td>52</td></tr><tr><td>4.1</td><td>Parameters used for our PnP-ReG method.</td><td>65</td></tr><tr><td>4.2</td><td>Quantitative results of super-resolution (x2) for the PnP approaches in terms of PSNR.</td><td>67</td></tr><tr><td>4.3</td><td>Quantitative results of super-resolution (x3) for the PnP approaches in terms of PSNR.</td><td>68</td></tr><tr><td>4.4</td><td>Quantitative results of deblurring for the PnP approaches in terms of PSNR.</td><td>69</td></tr><tr><td>4.5</td><td>Quantitative results of pixel-wise inpainting for the PnP approaches in terms of PSNR.</td><td>72</td></tr><tr><td>4.6</td><td>Parameters used for unrolled gradient descent optimization.</td><td>74</td></tr><tr><td>4.7</td><td>Quantitative results of super-resolution for unrolled gradient descent in terms of PSNR.</td><td>75</td></tr><tr><td>4.8</td><td>Quantitative results of deblurring for unrolled gradient descent in terms of PSNR.</td><td>76</td></tr><tr><td>6.1</td><td>Quantitative results of denoising omnidirectional data (Spherical format, ERP and CMP) in terms of WS-PSNR.</td><td>101</td></tr><tr><td>6.2</td><td>Quantitative results of denoising omnidirectional data (Spherical format, ERP and CMP) in terms of S-PSNR.</td><td>102</td></tr><tr><td>6.3</td><td>Comparison of denoising equirectangular images in terms of WS-PSNR.</td><td>103</td></tr><tr><td>A.1</td><td>Parameters used for the PnP-ADMM in Chapter 4.</td><td>110</td></tr><tr><td>A.2</td><td>Parameters used for One-Net in Chapter 4.</td><td>111</td></tr><tr><td>A.3</td><td>Parameters used for RED in Chapter 4.</td><td>112</td></tr></tbody></table>
# RÉSUMÉ

Les images omnidirectionnelles, ou images sphériques, sont des données visuelles qui couvrent un champ de vision à 360°, capturant le champ lumineux convergeant vers un seul point depuis toutes les directions. Les images omnidirectionnelles ont fait l'objet d'une attention particulière ces dernières années, notamment dans les domaines de réalité augmentée (RA) et la robotique. Ces caméras souffrent souvent d'une distorsion en raison de leur large champ de vision et de leurs caractéristiques optiques uniques. Par conséquent, la qualité des images omnidirectionnelles capturées peut être degradée, ce qui nécessite le développement de solutions aux problèmes inverses spécifiques pour ces données.

D'une manière générale, les problèmes inverses font référence à la tâche consistant à trouver les entrées qui ont produit un ensemble d'observations. Plus particulièrement dans le domaine de l'imagerie, les limitations des systèmes d'acquisition, de transmission, ou encore de stockage peuvent conduire à des observations dégradées. Parfois, l'acquisition peut être coûteuse ou même nocive (par exemple, exposer le corps humain à des rayonnements dans le cas d'un scanner), il est donc souhaité de mesurer le moins de données possible. Ainsi, les problèmes inverses en imagerie concernent la récupération d'une image à partir des mesures dégradées que nous observons.

Les problèmes inverses rencontrés sont souvent mal-posés, ce qui signifie qu'ils n'ont pas de solution unique. Afin de résoudre le défi posé par la nature mal posée, des connaissances a priori supplémentaires sont utilisées pour restreindre l'ensemble des solutions admissibles.

Cette thèse est dédiée à contribuer aux méthodes de régularisation basées sur l'apprentissage profond pour les problèmes inverses en imagerie perspective et omnidirectionnelle.

## Motivation

Au cours des dernières décennies, des recherches approfondies ont été consacrées à la résolution de problèmes inverses, aboutissant à diverses approches différentes pour leur résolution. Par exemple, l'approche bayésienne consiste à considérer une perspective probabiliste et à définir l'a priori et le terme d'attache aux données par des distributions
de probabilité. La distribution a posteriori est exprimée à l'aide de la loi de Bayes :

$$p(x|y) = \frac{p(y|x)p(x)}{p(y)}, \qquad (1)$$

où $p(y|x)$ et $p(x)$ représentent respectivement la *vraisemblance* des observations et la distribution *a priori*.

La distribution a posteriori représente toutes les solutions possibles au problème inverse en question, compte tenu des données mesurées. Elle peut être évaluée à l'aide de différentes estimations, comme le MAP ou le MMSE.

Pour une bonne reconstruction, le choix du prior (la terme de régularisation) joue un rôle important. Bien que les méthodes de régularisation créées manuellement aient prouvé leur efficacité, il n'est pas facile de représenter entièrement la nature complexe des images naturelles à l'aide de telles distributions créées manuellement.

L'apprentissage profond a considérablement révolutionné les sciences de l'imagerie, et le domaine des problèmes inverses n'est pas resté insensible à cette influence. Notre attention se porte sur les progrès en matière de régularisation apprise des problèmes inverses, un domaine qui a récemment suscité beaucoup d'intérêt. Certaines de ces approches comprennent :

* **Régularisation en utilisant le prior implicitement capturé par l'architecture du réseau** : Le Deep Image Prior (DIP) [1] est une approche unique, où les auteurs montrent que l'architecture du réseau d'un modèle génératif peut elle-même être utilisée pour régulariser les problèmes inverses. Le réseau générateur n'est pas appris sur un ensemble de données. Il fonctionne en prenant un vecteur d'entrée aléatoire et en ajustant itérativement ses poids initialisés de manière aléatoire pour minimiser l'erreur quadratique moyenne entre sa sortie et l'observation dégradée. Cela signifie que le réseau capture implicitement un prior au cours du processus d'optimisation, régularisant efficacement les tâches de restauration d'image.

* **Régularisation en utilisant la puissance des débruiteurs** : L'idée d'exploiter les débruiteurs pour régulariser différents problèmes inverses a été proposée pour la première fois en 2013 par Venkatakrishnan et al. dans un framework appelé *Plug-and-play*, où l'opérateur proximal d'un régulariseur est remplacé par un débruiteur. La méthode a été initialement proposée avec un débruiteur classique, mais l'approche a rapidement évolué et est devenue un sous-domaine de recherche important. Les approches Plug-and-play avec des débruiteurs appris sont des méthodes
de pointe pour résoudre les problèmes inverses [2], [3]. Une autre catégorie de méthodes exprime le gradient de la fonction de régularisation comme étant le résidu du débruitage de l'image estimée [4]. Les débruiteurs ont également été utilisés pour remplacer la fonction de score dans les approches d'échantillonnage de pos- térieur [5]-[8]. Ainsi, les approches basées sur le débriteur constituent sans aucun doute une partie importante de la résolution des problèmes inverses d'une manière générique. Un inconvénient commun à ces approches est que la régularisation est implicitement ou explicitement définie par un débriteur. Cela implique un hyper- paramètre supplémentaire à affiner pour chaque application spécifique, alors qu'en théorie, un régulariseur devrait représenter l'a priori indépendamment du problème inverse spécifique en question.

* **Algorithmes de déroulement profond** : Les méthodes de déroulement profond consistent à itérer un algorithme d'optimisation pour un nombre prédéterminé d'ité- rations en remplaçant la fonction de régularisation par un réseau de neurones. Le réseau est entraîné de bout en bout pour un problème inverse spécifique. Ces mé- thodes donnent d'excellents résultats. Cependant, leur principal inconvénient est qu'elles ne sont pas génériques.

D'autre part, un nombre limité de travaux ont abordé le thème des problèmes inverses dans les images omnidirectionnelles, y compris la tâche de débruitage. Cela limite notre capacité à étendre facilement les méthodes de pointe à la sphère.

De plus, les images omnidirectionnelles sont souvent projetées sur des représentations planaires bidimensionnelles (par example en utilisant la projection équirectangulaire [9]) afin de profiter des outils et des méthodes développées pour les images en 2D. Cependant, ces projections entraînent d'importantes distorsions et une résolution spatiale non uni- forme. Par conséquent, l'utilisation d'algorithmes conçus pour les images 2D sans tenir compte de la géométrie sphérique conduit à des résultats limités. Une autre approche consiste à travailler directement sur la sphère. Cependant, définir des outils de traitement sur la sphère n'est pas une tâche simple.

L'objectif de cette thèse est de contribuer aux méthodes de régularisation basées sur l'apprentissage profond. Nous nous concentrons principalement sur les approches Plug-and-play, motivées par leur caractère générique et leurs performances exceptionnelles. Nous étudions également l'effet de la combinaison de différentes approches de régularisa- tion. En ce qui concerne les images omnidirectionnelles, nous nous concentrons sur la tâche spécifique du débruitage, étant donné que c'est l'élément clé des méthodes Plug-and-play.
# Contributions

Dans cette section, nous discutons des principales contributions de cette thèse et ré-
pertorions les publications associées.

Régularisation du Deep Image Prior avec un débriteur profond

Nous proposons une méthode d'optimisation combinant un débriteur appris avec le mo-
dèle génératif non entraîné, le Deep Image Prior (DIP), dans le cadre de la *Alternating
Direction Method of Multipliers* (ADMM). Nous comparons également différents régula-
risateurs de l'optimisation de DIP, pour des problèmes inverses en imagerie. L'objectif
est d'étudier l'effet de la combinaison du DIP non entraîné et d'un régulariseur générique
appris de manière supervisée à partir d'une grande collection d'images. Le débriteur
est utilisé comme un opérateur proximal dans le cadre de l'ADMM et peut être appris
indépendamment du problème inverse considéré. La régularisation proposée basée sur le
débriteur présente des avantages par rapport aux a priori fabriqués à la main et reste
générique, car elle peut être appliquée à différents problèmes inverses.

PnP-ReG : Gradient de régularisation appris pour la descente de
gradient Plug-and-play

Le cadre Plug-and-play (PnP) permet d'intégrer des priors avancés de débruitage d'images
dans des algorithmes d'optimisation afin de résoudre efficacement diverses tâches de
restauration d'images. Cependant, l'algorithme PnP original ne s'applique qu'aux algo-
rithmes proximaux. De plus, un hyperparamètre supplémentaire doit être ajusté pour
chaque application lorsque la régularisation est définie par un débriteur. Nous nous
appuyons sur l'hypothèse qu'un débriteur représente l'opérateur proximal d'un régula-
risateur différentiable, et nous établissons une relation entre le débriteur et le gradient
du régularisateur correspondant. Nous utilisons cette relation pour entraîner un réseau
modélisant directement le gradient du régularisateur, tout en apprenant conjointement le
débriteur correspondant. Nous utilisons ce réseau dans un algorithme de PnP avec la
descente de gradient et obtenons de meilleurs résultats par rapport à d'autres approches
génériques. Nous montrons également que le réseau de régularisation peut être utilisé
comme un réseau pré-entraîné pour la descente de gradient déroulée. Enfin, nous montrons
que le débriteur résultant permet une meilleure convergence de l'ADMM Plug-and-play
en pratique.

# SphereDRUNet : un débruiteur sphérique pour les images omnidirectionnelles

Nous abordons le problème du débruitage d'images omnidirectionnelles et nous visons à étudier l'avantage de débruiter directement l'image sphérique plutôt que sa projection. Nous introduisons un nouveau réseau appelé SphereDRUNet pour débruiter des images sphériques en utilisant des outils d'apprentissage profond sur un échantillonnage sphérique. Nous montrons que le débruitage direct de la sphère à l'aide de notre réseau donne de meilleures performances par rapport au débruitage des images équirectangulaires avec un modèle appris de manière similaire. Nous comparons également notre SphereDRUNet à un réseau de convolution basé sur des graphes et confirmons que l'approche de convolution que nous utilisons est plus efficace.

## Publications

* R. Fermanian, T. Maugey and C. Guillemot. SphereDRUNet : A Spherical Denoiser for Omnidirectional Images. IEEE 22nd International symposium on mixed and augmented reality adjunct (ISMAR-Adjunct), Sydney, Australia, 2023.

* R. Fermanian, M. Le Pendu and C. Guillemot. PnP-ReG : Learned Regularizing Gradient for Plug-and-Play Gradient Descent. SIAM Journal on Imaging Sciences 16.2 (2023) : 585-613.

* R. Fermanian, M. Le Pendu and C. Guillemot. Regularizing the Deep Image Prior with a Learned Denoiser for Linear Inverse Problems. IEEE 23rd International Workshop on Multimedia Signal Processing (MMSP), Tampere, Finland, 2021.

## Contenu

La suite de cette thèse est divisé en 6 chapitres.

Dans le Chapitre 2, nous élaborons le contexte des problèmes inverses pour l'imagerie. Nous commençons par introduire des problèmes inverses et illustrons plusieurs exemples. Nous discutons ensuite de différentes méthodologies de reconstruction et approches de régularisation.
Dans le chapitre 3, nous proposons une méthode combinant un débriteur appris et le Deep Image Prior dans le cadre de l'ADMM. Nous discutons de différentes approches existantes qui régularisent le DIP et montrons l'avantage de notre proposition de régulation appris.

Dans le chapitre 4, nous proposons une nouvelle méthode pour entraîner un réseau modélisant le gradient d'un régulariseur, conjointement avec un débriteur. En partant de l'hypothèse qu'un débriteur représente l'opérateur proximal d'un régulariseur différentiable (définissant l'image a priori), nous dérivons une fonction de perte qui relie le débriteur et le gradient du régulariseur correspondant. Nous utilisons ce réseau dans un algorithme de Plug-and-play avec de la descente de gradient et obtenons de meilleures performances par rapport à d'autres méthodes génériques. Nous utilisons ensuite notre réseau comme stratégie de pré-entraînement dans un algorithme de descente de gradient déroulé.

Dans le Chapitre 5, nous présentons les principaux concepts du traitement d'images omnidirectionnelles. Nous discutons des défis liés à l'acquisition, à la représentation et au traitement des images omnidirectionnelles, et présentons enfin la littérature sur les problèmes inverses.

Dans le chapitre 6, nous proposons SphereDRUNET, un nouveau réseau de débruitage d'images omnidirectionnel, en transférant le DRUNet [2] sur la sphère. Nous montrons que le débruitage d'image donne de meilleures performances lorsqu'il est effectué directement sur la sphère plutôt que via une projection.

Le chapitre 7 conclut cette thèse et propose des perspectives d'études ultérieures.
# INTRODUCTION

In this chapter, we present the main problem of interest of this thesis: inverse problems in imaging, with particular attention to omnidirectional imaging. We elaborate on the context and discuss the motivation behind our work. Then, we list our contributions and publications, and outline the remaining of the manuscript.

## 1.1 Context

Omnidirectional images, also known as 360° images or spherical images, are high-resolution visual data that cover a 360° field of view capturing the light field converging to a single point from all directions. From Augmented Reality (AR) to robotics, omnidirectional images have gained a particular attention in the recent years. These cameras often suffer from inherent distortion and increased noise due to their wide field of view and unique optical characteristics. Consequently, the quality of the captured omnidirectional images can be compromised, limiting the accuracy of subsequent computer vision tasks in AR/MR and necessitating the development of solutions for inverse problems in omnidirectional images.

In a general sense, inverse problems refer to the task of finding the input or underlying factors that have produced a particular set of data or observations. More particularly in the area of imaging, limitations in acquisition systems, transmission or even storage can lead to degraded measurements. Sometimes, acquisition can be costly or even harmful (i.e. exposing the human body to radiation in the case of CT scan), so it is desired to measure as few data points as possible. Hence, inverse problems in imaging concern recovering an image from its degraded measurements that we observe.

Unfortunately, solving inverse problems is not a straightforward task. In fact, the presence of noise and the incomplete measurements lead to challenges in the reconstruction. Moreover, these problems are *ill-posed*, meaning that they do not have a unique solution that is consistent with the observation. In order to address the ill-posed nature,
additional prior knowledge (i.e. regularization) is used to restrict the set of admissible solutions.

This thesis is dedicated to contribute to deep learning-based regularization methods for inverse problems in perspective and omnidirectional imagery.

## 1.2 Motivation

Over the past few decades, extensive research has been dedicated to solve inverse problems, resulting in various different approaches for their resolution. For instance, the Bayesian approach consists in considering a probabilistic perspective and defining the prior and the data term with probability distributions. The posterior distribution is the probability distribution of the unknown parameters x we wish to recover, conditional on the observed measurement y. It is expressed using Bayes' theorem as:

$$p(x|y) = \frac{p(y|x)p(x)}{p(y)}, \qquad (1.1)$$

where $p(y|x)$ and $p(x)$ represent respectively the likelihood of the observations and the prior. The posterior distribution represents all potential solutions to the inverse problem in hand, given the measurement. It can be evaluated using different point estimates, such as the MAP or the MMSE.

For a successful reconstruction, the choice of the prior plays a crucial role. Early methods used hand-crafted priors incorporating desired properties of the solution such as Total Variation (TV) [10] or wavelet decomposition [11]. Although these priors have demonstrated significant efficiency, it is not trivial to fully represent the complex nature of natural images using hand-crafted priors.

The breakthrough of deep learning has brought a remarkable shift in different areas of imaging sciences, including the solution of inverse problems. Our focus lies in recent advances on deep regularization of inverse problems, an area that has recently gained a lot of attention. Some of these approaches include:

* Regularizing using the implicit prior captured by the network architecture: The Deep Image Prior (DIP) [1] is a unique approach, where the authors show that the network architecture of a generative model itself can be used to regularize inverse problems. The generator network is not trained on a dataset. Instead, it operates by taking a random input vector and iteratively adjusting its randomly
initialized weights to minimize the mean square error between its output and the
degraded observation. This means that the network inherently captures an im-
plicit prior during the optimization process, effectively regularizing image restora-
tion tasks. In subsequent studies, researchers have aimed to improve the Deep
Image Prior’s performance in various ways, particularly by regularizing it using
hand-crafted regularizers.

*   **Regularizing using the power of denoising engines:** The idea of leveraging denoisers to regularize different inverse problems was first proposed in 2013 by Venkatakrishnan et al. in a framework called *Plug-and-play*, where the proximal operator of a regularizer is replaced by a denoiser. The method was initially proposed with a classical denoiser, but the approach quickly evolved and became an attractive sub-field of research. Deep denoiser-based Plug-and-play approaches are state-of-the-art methods for solving inverse problems [2], [3]. Another category of methods expresses the gradient of the regularization function to be the denoising residual of the estimated image [4]. Denoisers have also been used to replace the *score function* in posterior sampling approaches [5]–[8]. Thus, denoiser-based approaches are unarguably an important part of solving inverse problems in a generic way. A drawback shared by these approaches is that the regularization is implicitly or explicitly defined by a denoiser. This involves an additional hyper-parameter to be fine-tuned for each specific application, whereas in theory, a regularizer should represent the image prior independently of the specific inverse problem at hand.

*   **Deep unrolling:** It consists in iterating an optimization algorithm for a pre-determined number of iterations by replacing the regularization function by a neural network. The network is trained end-to-end for a specific inverse problem. Deep unrolling methods perform impressively well. However, their main drawback is that they are not generic.

On the other hand, it is worth noting that a limited number of works have addressed the topic of inverse problems in omnidirectional images, including the task of denoising. This limits our capability to easily extend state-of-the-art methods to the sphere.

Moreover, omnidirectional images are often mapped to two-dimensional planar rep-
resentations (e.g. using equirectangular projection [9]) in order to take advantage of
computing tools and methods that have been developed for 2D images. However, these
mappings cause significant distortions and non-uniform spatial resolution. Hence, using
algorithms designed for 2D images without considering the spherical geometry leads to
limited results. Another approach is to work directly on the sphere. However, defining
processing tools on the sphere is not a straightforward task, as it comes with the challenges
of spherical geometry.

The goal of this thesis is to contribute to deep learning-based regularization meth-
ods. We mainly focus on Plug-and-play approaches, motivated by their genericity and
outstanding performance. We also study the effect of combining different regularization
approaches. For omnidirectional images, we focus on the specific task of denoising, given
that it is the key element of Plug-and-play methods.

1.3 Contributions

In this section, we discuss the main contributions of this thesis and list the associated
publications.

Regularization of the Deep Image Prior with a deep denoiser

We propose an optimization method coupling a learned denoiser with the untrained gener-
ative model deep image prior (DIP) in the framework of the Alternating Direction Method
of Multipliers (ADMM). We also compare different regularizers of DIP optimization, for
inverse problems in imaging. The goal is to study the effect of combining the untrained
DIP and a generic regularizer learned in a supervised manner from a large collection of
images. When placed in the ADMM framework, the denoiser is used as a proximal oper-
ator and can be learned independently of the considered inverse problem. The proposed
denoiser-based regularization compares favourably with handcrafted priors and remains
generic, since it can be applied to different inverse problems.

PnP-ReG: learned regularizing gradient for Plug-and-play gradi-
ent descent

The Plug-and-play (PnP) framework makes it possible to integrate advanced image de-
noising priors into optimization algorithms, to efficiently solve a variety of image restora-
tion tasks. However, the original PnP algorithm only applies to proximal algorithms.
When working on the previous contribution, we observed that the parameter tuning of
the ADMM algorithm is not a straightforward task. Also, an additional hyper-parameter
must be tuned per application when the regularization is defined by a denoiser. We rely
on the assumption that a denoiser represents the proximal operator of a differentiable
regularizer, and derive a relation between the denoiser and the corresponding regular-
izer's gradient. We use this relation to train a network directly modeling the gradient of
the regularizer, while jointly training the corresponding denoiser. We use this network in
a PnP gradient descent algorithm and obtain better results comparing to other generic
approaches. We also show that the regularizer can be used as a pre-trained network for
unrolled gradient descent. Lastly, we show that the resulting denoiser allows for a better
convergence of the Plug-and-play ADMM.

SphereDRUNet: a spherical denoiser for omnidirectional images

We address the problem of omnidirectional image denoising and we aim to study the
advantage of denoising the spherical image directly rather than its mapping. We introduce
a novel network called SphereDRUNet to denoise spherical images using deep learning
tools on a spherical sampling. We show that denoising directly the sphere using our
network gives better performance, compared to denoising the projected equirectangular
images with a similarly learned model. We also compare our SphereDRUNet to a graph-
based convolution network and confirm that the convolution approach that we use is more
efficient.

List of publications

* R. Fermanian, T. Maugey and C. Guillemot. SphereDRUNet: A Spherical Denoiser for Omnidirectional Images. IEEE 22nd International symposium on mixed and augmented reality adjunct (ISMAR-Adjunct), Sydney, Australia, 2023.

* R. Fermanian, M. Le Pendu and C. Guillemot. PnP-ReG: Learned Regularizing Gradient for Plug-and-Play Gradient Descent. SIAM Journal on Imaging Sciences 16.2 (2023): 585-613.

* R. Fermanian, M. Le Pendu and C. Guillemot. Regularizing the Deep Image Prior with a Learned Denoiser for Linear Inverse Problems. IEEE 23rd International Workshop on Multimedia Signal Processing (MMSP), Tampere, Finland, 2021.

1.4 Outline

The remaining of this thesis is divided into 6 chapters.
In Chapter 2, we elaborate on the background of inverse problems for perspective imagery. We start by introducing inverse problems and illustrate several examples. We then discuss different reconstruction methodologies and regularization approaches.

In Chapter 3, we propose a method coupling a learned denoiser with the Deep Image Prior in the ADMM framework. We discuss different existing approaches that regularize the DIP and show the advantage of our proposed deep regularization.

In Chapter 4, we propose a novel method to train a network modeling the gradient of a regularizer, jointly with a deep denoiser. Based on the assumption that a denoiser represents the proximal operator of an underlying differentiable regularizer (defining the image prior), we derive a loss function that links the denoiser and the gradient of the corresponding regularizer. We use this network in a Plug-and-play gradient descent algorithm and obtain better performance comparing to other generic methods. We then use our network as a pre-training strategy in an unrolled gradient descent algorithm.

In Chapter 5, we present the main concepts of omnidirectional image processing. We discuss the acquisition, representation and processing challenges of omnidirectional images, and finally present the literature of inverse problems.

In Chapter 6, we propose SphereDRUNET, a novel omnidirectional image denoising network, by transferring the DRUNet [2] to the sphere. We show that image denoising gives better performance when it is performed directly on the sphere rather than via a mapping, even though it raises many challenges.

Chapter 7 concludes this thesis and proposes perspectives for further studies.
# PART I

## Inverse problems in perspective imagery
<figure><img src="image_3.png" /></figure>
# AN OVERVIEW OF INVERSE PROBLEMS IN PERSPECTIVE IMAGERY

Picture a doctor seeking to understand the inner workings of the human body from a limited set of diagnostic measurements or an astronomer striving to discover the characteristics of celestial objects, relying only on the faint traces of light that reach its telescope: these are scenarios illustrating instances that represent *inverse problems*. In most real-life applications, we do not have direct measurements of the signal of interest, but only a corrupted version of the original signal. Several phenomena may induce these distortions in practice, such as the acquisition system that is used or the communication channel through which the signal is transmitted. In such cases, we aim to reconstruct the original unknown signal from the distorted observations that we have.

Inverse problems have been broadly studied in the last few decades and several different sort of approaches have been proposed to solve them. This chapter aims, in a first place, to introduce inverse problems in signal and image processing, and secondly, to present an overview of the research and the methodologies that have been proposed to solve them in the literature.

## 2.1 Inverse Problems

### 2.1.1 Introduction and problem statement

Inverse problems refer to the task of recovering an image $x \in \mathbb{R}^n$ from its degraded measurements $y \in \mathbb{R}^m$ obtained by a certain degradation model:

$$y = \mathcal{H}(Ax), \qquad (2.1)$$
where $A \in \mathbb{R}^{m \times n}$ represents a linear degradation operator depending on the inverse problem and $\mathcal{H}$ is a non-linear model representing noise degradation.

Various different types of non-linear degradation may destroy the signal of interest. For example, the signal may suffer from multiplicative noise [12]–[14], such as speckle noise [15], [16]. The latter can be found in a wide range of systems, including synthetic aperture radar (SAR) images, ultrasound imaging, and many more. Moreover, the statistical nature of electromagnetic waves such as x-rays, visible lights and gamma rays causes shot noise. These ray sources emit a number of photons per unit time and have a random fluctuation of photons. The image has thus a spatial and temporal randomness. The resulting noise is signal-dependent and non-additive and it can be modeled by a Poisson distribution [17]–[20]. Furthermore, acquisition devices usually encounter internal fluctuations. This causes an additive noise that can be represented by a Gaussian distribution [21], [22]. Another possibility is the combination of Poisson and Gaussian noise (called Poisson-Gaussian noise) [23]–[26]. Figure 2.1 shows an image degraded with the aforementioned noise types.

In this work, we consider the case where the noise is additive (independent from the signal itself) and has a Gaussian distribution. Thus, the degraded signal can be obtained by:

$$y = Ax + \eta, \qquad (2.2)$$

where $\eta \in \mathbb{R}^m$ represents Additive White Gaussian Noise (AWGN). The probability density function (PDF) of a Gaussian noise follows a normal distribution and is given by:

$$\varphi(\eta) = \frac{1}{\sqrt{2\pi}\sigma^2} \exp\left(-\frac{(\eta - \mu)^2}{2\sigma^2}\right) \qquad (2.3)$$

where $\mu$ is the mean of the Gaussian distribution and $\sigma$ is its standard deviation (i.e. $\sigma^2$ is its variance). Figure 2.2 illustrates the PDF of a Gaussian noise.

The rationale behind commonly restricting the noise to Additive White Gaussian Noise (AWGN) is as follows: one might argue that the Poisson distribution is more suited to model imaging noise given that imaging sensors inherently quantify photons. Although this argument is correct, it is worth noting that in situations of high photon counts, the Poisson distribution becomes a Gaussian one [28]. As a result, when dealing with high photon counts, the measurements can be converted into a form akin to Gaussian contamination by using a variance-stabilizing transformation, such as the *Anscombe transform* [29]. Finally, since the formulations become simpler when using the Gaussian distribution,
<figure><img src="image_4.png" /><figcaption>Figure 2.1 – (a) Grayscale version of the pepper image from the Set14 dataset [27] (pixel intensities are ranged from 0 to 255). The original image is degraded with (b) zero mean Speckle noise of standard deviation 0.1, (c) Poisson noise with average intensity of photons equal to 1, (d) additive zero mean Gaussian noise of standard deviation 25 and (e) Poisson-Gaussian noise (poisson noise with average intensity of photons of 1, and additive zero mean Gaussian noise of standard deviation of 25).</figcaption></figure>

researchers have considered the AWGN when designing solutions to inverse problems.

## 2.1.2 Inverse problems in imaging: examples and definitions

In the image restoration tasks that we consider, we aim to recover the original image $x$, which is degraded by a certain operation that is represented by the degradation operator $A$. The latter depends on the inverse problem in hand. In this section, we define some inverse problems in imaging that we will solve in this work, such as denoising, super-resolution, deblurring and pixel-wise inpainting. Solution methodologies will be discussed in the upcoming sections.
<figure><img src="image_5.png" /><figcaption>Figure 2.2 – Probability Density Function of Gaussian noise of mean $\mu$ and standard deviation $\sigma$.</figcaption></figure>

### 2.1.2.1 Denoising

Image denoising is unarguably the most extensively studied among the inverse problems in imaging. The task of denoising refers to the process of removing unwanted noise from an observed image to enhance its quality and improve visual clarity. The presence of noise in images can be attributed to various factors such as sensor limitations, transmission errors, or environmental conditions during image acquisition. The different types of noise were discussed in the previous section.

In the case where the observed image is degraded with Additive White Gaussian noise of variance $\sigma^2$, Eq. 2.2 reduces to $y = x + \eta$, where $y$ is the noisy observation, $x$ is the original image and $\eta$ is an AWGN. Therefore, the degradation matrix $A$ is simply the identity matrix. Figure 2.3 illustrates the inverse problem of image denoising.

### 2.1.2.2 Deblurring

Image deblurring refers to the task of recovering a sharp, undistorted version of a blurred image. Blur in images can occur due to various factors, such as camera movement, defocusing, or atmospheric conditions. The objective of image deblurring is to estimate
<figure><img src="image_6.png" /><figcaption>Figure 2.3 – Inverse problem of denoising an image corrupted by a White Gaussian noise of variance σ² = 30.</figcaption></figure>

the original, sharp image that might have led to the blurry observation. Figure 2.4 illustrates the inverse problem of deblurring.

The degradation matrix *A* represents a blurring operator (i.e. convolution), typically followed by the addition of Gaussian noise. The blurring filter used in image deblurring is commonly known as the Point Spread Function (PSF). The PSF characterizes how light from each point in the original scene is spread or blurred across adjacent pixels in the resulting image due to factors like camera shake, defocus, or motion. The PSF is usually represented by a rectangular or a square matrix, with its size increased by the extent of blurring. For instance, in the case of motion blur caused by camera movement, the PSF might be a horizontal line kernel indicating the motion direction. For defocus blur, the PSF may have a circular shape. The PSF of the blurring process must be known or estimated accurately in order to perform image deblurring. Blind deconvolution techniques attempt to estimate the PSF simultaneously with the sharp image, making the reconstruction more complex and challenging.

### 2.1.2.3 Super-resolution

Super-resolution is the process of enhancing the resolution and quality of a given low-resolution image or video to obtain a higher-resolution version with more details and clarity. The goal is to reconstruct missing high-frequency details, edges, and textures that are lost during the downscaling process of a given image.

Generally, the low-resolution image is obtained after a blurring and a downsampling operation of the high-resolution image. In practice, noise may also be present in the acquired low-resolution. Hence, the degradation operator includes blurring, downsampling,
<figure><img src="image_7.png" /><figcaption>Figure 2.4 – Inverse problem of deblurring. The blurred image is generated using an isotropic Gaussian kernel followed by adding Gaussian noise of standard deviation $\sigma = \sqrt{2}$.</figcaption></figure>

and possibly noise. Figure 2.5 shows an example of super-resolution.

<figure><img src="image_8.png" /><figcaption>Figure 2.5 – Inverse problem of super-resolution. The low-resolution image is obtained by performing Gaussian blur over the original image and downscaling it with a factor 2.</figcaption></figure>

### 2.1.2.4 Pixel-wise inpainting

Pixel-wise inpainting, also known as image completion, is a task in image processing that aims to fill in missing pixels in an image. The goal is to generate a complete and visually plausible image from the partially observed or degraded version.

The missing pixels are often represented by masked pixels, where the pixel values are unknown or marked as invalid. The task is to estimate the missing pixel values and complete the image in a way that blends seamlessly with the surrounding content. The
degraded image is typically obtained by multiplying the original image with a binary mask. This binary mask is used to indicate the pixels that need to be inpainted or completed. The mask is a matrix of the same size as the original image, with values of 0 or 1. A value of 0 in the mask denotes a missing or corrupted pixel, while a value of 1 represents a valid pixel that should be kept unchanged. Figure 2.6 shows an example of the inverse problem of completion.

<figure><img src="image_9.png" /><figcaption>Figure 2.6 – Inverse problem of pixel-wise inpainting. The degraded image has a known pixel rate of $p = 10\%$.</figcaption></figure>

### 2.1.3 An ill-posed problem and need for regularization

The task of solving inverse problems defined by Eq. 2.2 is equivalent to finding an optimal solution $\hat{x}$ that best represents the original signal of interest $x$ using the measurement $y$. In practice, even in relatively simpler applications, obtaining a restored image of good quality is not as easy as it seems. This is due to the ill-posed nature of inverse problems, that we will elaborate in this section.

A problem is considered to be mathematically well-posed according to Hadamard [30] if the solution $\hat{x}$ satisfies all three conditions below:

* A solution exists

* The solution is unique

* The solution is stable (the solution depends continuously on the measurement)

The first condition requires the system to be consistent. The second condition restricts the solution to be unique. It ensures that there is only one possible solution that can be
determined from the available data. In order to meet this condition, the system should not be under-determined (e.g. the degradation matrix should not have $m < n$). Moreover, the stability of the solution with respect to the input data ensures that minor changes (i.e. noise) in the observed data do not drastically alter the solution. Stability is essential to prevent the problem from becoming ill-posed, where small errors in the data could lead to significant inaccuracies or even to the non-existence of the solution.

Inverse problems are ill-posed because they do not satisfy all of the aforementioned conditions. To address this ill-posed nature, regularization techniques are employed. Regularization involves incorporating additional prior knowledge or constraints into the problem to guide the solution process. By imposing constraints, such as smoothness or sparsity on the unknowns, regularization helps to narrow down the solution space, making it more stable and robust. The topic of regularization will be further detailed in the following sections.

## 2.2 Reconstruction approaches for solving inverse problems

To this end, we have discussed inverse problems without actually diving into the actual task of their reconstruction. So, how can we solve inverse problems? From classical approaches such as simple interpolation and filtering to deep learning methods, the answer to this question has evolved over time with the advancements in signal and image processing. In this section, we elaborate the most common reconstruction approaches that have been adopted in the literature.

### 2.2.1 Variational approach for solving inverse problems

A natural and naive approach to solve inverse problems is to apply a function enforcing similarity between the observation and the searched solution (e.g. the least squares method), as follows:

$$ \hat{x} = \underset{x \in \mathbb{R}^n}{\operatorname{argmin}} \mathcal{L}(Ax, y), \qquad (2.4) $$

However, since inverse problems are ill-posed, we use a regularizing function as discussed before. The variational approach for solving inverse problem amounts to solving an opti-
mization problem defined by:

$$ \hat{x} \in \operatorname*{argmin}_{x \in \mathbb{R}^n} \mathcal{L}(Ax, y) + \lambda \phi(x), \qquad (2.5) $$

where $\mathcal{L}(Ax, y)$ represents the data-fidelity term, $\phi(x)$ is the regularization function and $\lambda > 0$ balances the trade-off between the two terms. This parameter should be tuned in order to compromise between the two terms.

Note that the most common choice for $\mathcal{L}(Ax, y)$ is the $\mathcal{L}_2$ distance form, thus we generally get:

$$ \hat{x} \in \operatorname*{argmin}_{x \in \mathbb{R}^n} \|y - Ax\|_2^2 + \lambda \phi(x). \qquad (2.6) $$

The problem in Eq. 2.5 usually does not have a closed-form solution except in some simple cases, for instance when both the data-fidelity term and the regularization are quadratic. Consequently, the variational approach involves solving an optimization problem by using iterative schemes. Section 2.3 provides an overview of some iterative algorithms that are commonly used to solve variational problems.

The choice of the regularization function plays a crucial role in the solution of inverse problems. The prior should balance between incorporating desired properties of the solution and avoiding an excessive computational burden. In the following, we provide a non exhaustive discussion about some of the traditional priors (often referred to as "hand-crafted" priors) that have been designed throughout years.

* **Tikhonov Regularization:** Tikhonov regularization [31] is probably one of the most commonly used classical regularization methods to solve inverse problems. The regularization expression is given as $\phi(x) = \|\Gamma x\|_2^2$, where $\Gamma$ is the Tikhonov operator promoting some desirable properties. For example, if the signal is believed to be smooth, $\Gamma$ can assume the role of a high-pass operator like a Laplacian filter (i.e. $\phi(x) = \|Lx\|_2^2$). A particular case of Tikhonov regularization is the energy regularization ($L_2$ regularization) when $\Gamma = I$. Thus $\phi$ reduces to $\phi(x) = \|x\|_2^2$ encouraging solutions with smaller norms.

* **Weighted Smoothness:** Imposing a uniform constraint on the smoothness of the solution may result in overly-smoothed regions in the image. To overcome this, a weighted smoothing regularization can be considered, such as in [32], [33]. By assigning varying weights, the regularization process can effectively adapt to the inherent structure of the data: a high weight is assigned to regions where the
assumption of smoothness is believed to be valid, while a lower weight is allocated
to textured regions or edges. Forming a diagonal weight matrix W containing the
assigned weights in the diagonal, and considering the Laplacian regularization, the
regularization function becomes $\phi(x) = \|Lx\|_W^2$.

*   **Total-Variation:** Total-Variation [10] is a regularization function quantifying the intensity variation or discontinuity across neighboring pixels, hence promoting gradient sparsity. In the anisotropic form, total variation computes the $L_1$ norm of the image gradient: $TV_{aniso}(u) = \sum_{i=1}^{n} |(D_h u)_i| + |(D_v u)_i|$, while isotropic TV corresponds to the $L_2$ norm of the gradient vector: $TV_{iso}(u) = \sum_{i=1}^{n} \sqrt{(D_h u)_i^2 + (D_v u)_i^2}$. The variables $D_h$ and $D_v$ denote the first order finite difference discrete operators along the horizontal and vertical axes. Extensions of this line of research can be found in [11], [34].

*   **Sparsity methods:** A sparse signal refers to a type of signal in which most of the values are either zero or very close to zero, while only a few values are non-zero and carry a significant information. For example, wavelets provide a compact representation of data by using a minimal number of coefficients, thereby introducing the concept of sparsity, which conveniently links to the regularization $\phi(x) = \|Wx\|_1$ encouraging a reduction in the presence of non-zero Wavelet coefficients [35]. The concept of sparsity gained further interest as subsequent research explored redundant and learned representations. Ideal images are assumed to be linear combinations of atoms from a pre-defined dictionary $D$, denoted as $x = D\alpha$. Sparsity can hence be enforced through $\phi(x) = \|\alpha\|_0$. Different algorithms have been developed along these lines considering global and patch-based dictionaries [36]–[40].

*   **Self-similarity:** The concept of self-similarity in image processing is rooted in the idea that small patches within an image often exhibit strong similarities to other patches in the overall image content. Hence, treating these similar patches collectively is believed to contribute to enhanced restoration. For example, $\phi(x) = \sum_i \sum_{j \in \Lambda(i)} d\{R_i x, R_j x\}$ promotes proximity between $R_i x$ and $R_j x$, where $R_i x$ denotes the patch extracted from the image at location $i$, $\Lambda(i)$ refers to the set of indices in the neighborhood of $i$ and $d$ denotes a dissimilarity score.

*   **Low-rank assumption:** This approach is very similar to self-similarity regularization, but rather than imposing proximity between a set of closely related patches, a low-rank structure is forced on a matrix having these related patches as columns. This infers that these patches are spanned by a few main directions. For instance,
let $Z_{\Lambda(k)}$ be a matrix denoting a group of similar patches extracted from the image, where $\Lambda(.)$ is the collection of related patches and $k$ indicates the specific patch group. Let $\|\cdot\|_*$ denote the nuclear norm, which is used as a measure of the rank of the matrix. The low-rank assumption can be imposed by $\phi(x) = \sum_k \|Z_{\Lambda(k)}\|_*$ [41], [42].

It's important to note that the landscape of regularization techniques is vast and continually evolving. Hand-crafted regularization approaches have historically demonstrated significant efficiency in addressing specific challenges. However, in pursuit of enhanced performance, current research has shifted its focus towards learning prior knowledge. We detail the idea of leveraging deep learning tools to regularize inverse problems in Section 2.4.3.

### 2.2.2 Learning approach for solving inverse problems: Neural network regression

With the advancement of deep learning and the success of convolutional neural networks in different computer vision tasks, researchers started to develop learning-based solutions for inverse problems. The simplest approach is to train problem-specific neural networks that approximate the mapping between the space of measurements (i.e. degraded images) and the corresponding solution space (i.e. ground-truth images), without having knowledge of the degradation operator $A$ (an agnostic learner). In fact, using a large amount of training data, the network is able to learn the patterns and relationships between the distorted inputs and the corresponding clean outputs. As a result, the network can effectively reconstruct the high-quality outputs from given corrupted inputs. This approach yielded a significant performance improvement in the field of image restoration.

A wide range of network architectures have been designed and trained for different inverse problems, such as DnCNN [43], NLRN [44], FOCNet [45] and DRUNet [2] for image denoising, SRCNN [46], SRGAN [47], ESRGAN [48] and RCAN [49] for Super-Resolution and DCNN [50] for blind deconvolution.

However, the success of these networks highly depend on the complexity of the forward operator $A$ and the available amount of training data. For instance, in the medical imaging field, obtaining large datasets can be particularly challenging due to privacy concerns and ethical considerations. Another drawback of this end-to-end learning procedure is that
the network has to be re-trained again if a single parameter of the degradation process changes. Therefore, these limitations highlight the necessity for more adaptable and generalized methods in handling diverse degradation scenarios.

### 2.2.3 Bayesian approach for solving inverse problems

A common approach to solve inverse problems consists in considering a probabilistic perspective. Rather than treating $x$ and $y$ as constant values, it is assumed that the noise in the measured data causes some randomness, which further induces uncertainty in the restoration. It can be inferred that $x$ and $y$ are random variables having probability distributions. In this setting, the reconstructed image from the measurements can be calculated through its posterior density $p_{x|y}(x|y)$. For notation simplicity, we will denote $p_u(u)$ as $p(u)$ in the remainder of the manuscript. This approach is the core of this thesis.

Bayes' theorem expresses the posterior density $p(x|y)$ using the likelihood and the prior density:

$$p(x|y) = \frac{p(y|x)p(x)}{p(y)}, \qquad (2.7)$$

where $p(y|x)$ is the probability density of the measurement $y$ given $x$, representing the likelihood of the observations and $p(x)$ represents the prior distribution. $p(y)$ is the distribution of $y$ and can be considered as a normalization constant for $p(x|y)$, as it is independent of $x$.

In the assumption of an Additive White Gaussian Noise $\eta$ with standard deviation $\sigma$, the noise is assumed to be independent and identically distributed across pixels, which means that each pixel $y_i$ follows a Gaussian distribution with mean $Ax_i$ and standard deviation $\sigma$. Since the noise in each pixel is independent, the joint probability density function of all noisy pixels $y$ given the true image $x$ is the product of the PDFs of each pixel. Therefore, the likelihood can be expressed as:

$$ \begin{aligned} p(y|x) &= \prod_{i=0}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(y_i-Ax_i)^2}{2\sigma^2}} \\ &= \frac{1}{(2\pi\sigma^2)^{n/2}} e^{-\frac{\sum_{i=0}^{n} (y_i-Ax_i)^2}{2\sigma^2}} \\ &= \frac{1}{(2\pi\sigma^2)^{n/2}} e^{-\frac{1}{2\sigma^2} \|y-Ax\|_2^2}. \end{aligned} \qquad (2.8) $$

It therefore follows that $p(y|x) \propto e^{-\frac{1}{2\sigma^2} \|y-Ax\|_2^2}$.
Once the likelihood is derived and the prior is chosen, the posterior density can be evaluated using various point estimates, such as Maximum A Posteriori (MAP) and the Minimum Mean Square Error (MMSE).

The bayesian estimator can be obtained by calculating the minimum expectation of a certain cost function $C(\hat{x} - x)$, where the latter assesses the accuracy of the estimated $\hat{x}$ in comparison to the unknown signal $x$ with posterior distribution $p(x|y)$:

$$ \begin{aligned} \hat{x}_C &\in \underset{\hat{x}}{\operatorname{argmin}} \mathbb{E}_{x|y}[C(\hat{x} - x)], \\ \hat{x}_C &\in \underset{\hat{x}}{\operatorname{argmin}} \int C(\hat{x} - x)p(x|y) dx. \end{aligned} \qquad (2.9) $$

## Minimum Mean Squared Error estimator (MMSE)

For the cost function $C(\hat{x} - x) = \|\hat{x} - x\|_2^2$, solving Eq. 2.9 amounts to finding the Minimum Mean Squared Error (MMSE) estimate. Having the following minimization problem:

$$ \hat{x}_{\text{MMSE}} = \underset{\hat{x}}{\operatorname{argmin}} \int \|\hat{x} - x\|_2^2 p(x|y)dx, \qquad (2.10) $$

Let us note $Q(\hat{x}) = \int \|\hat{x} - x\|_2^2 p(x|y)dx$. We can thus state that finding the minimizer in Eq. 2.10 amounts to solving $\frac{\partial Q(\hat{x})}{\partial \hat{x}} = 0$.

$$ \begin{aligned} \frac{\partial Q(\hat{x})}{\partial \hat{x}} &= \int \frac{\partial \|\hat{x} - x\|_2^2}{\partial \hat{x}} p(x|y)dx, \\ &= \int 2(\hat{x} - x)p(x|y)dx, \\ &= 2\hat{x} - 2 \int xp(x|y)dx = 0. \end{aligned} \qquad (2.11) $$

Therefore, the MMSE estimator of the posterior mean is given by:

$$ \hat{x}_{\text{MMSE}} = \int xp(x|y)dx. \qquad (2.12) $$

While the formulation seems simple and clear, computing the MMSE estimate can be challenging, as it consists in estimating high dimensional integrals. Sampling algorithms are often used in order to perform numerical integration.
Maximum A Posteriori estimator (MAP)

Maximum A Posteriori estimator is derived from choosing the *Hit-or-miss* cost function:

$$ \mathcal{C}(\hat{x} - x) = \begin{cases} 0 & \text{if } \|\hat{x} - x\|_2 < \delta \\ 1 & \text{if } \|\hat{x} - x\|_2 \ge \delta, \end{cases} \qquad (2.13) $$

$\delta > 0$ being a predefined threshold.

In this case, the estimation in Eq. 2.9 reduces to:

$$ \hat{x}_{MAP} \in \underset{\hat{x}}{\operatorname{argmin}} \left(1 - \int_{\|\hat{x} - x\|_2 < \delta} p(x|y)dx.\right) \qquad (2.14) $$

As $\delta \to 0$ we can observe that this minimization is obtained when the posterior distribution is maximized. Therefore, the estimation is computed by maximizing the probability density of the posterior $p(x|y)$:

$$ \begin{align} & \underset{x}{\operatorname{argmax}} p(x|y) \nonumber \\ & \underset{x}{\operatorname{argmax}} \frac{p(y|x)p(x)}{p(y)} \tag{2.15} \\ & \underset{x}{\operatorname{argmin}} -\log p(y|x) - \log p(x) \nonumber \end{align} $$

As established in Eq. 2.8, $p(y|x) \propto e^{-\frac{1}{2\sigma^2}\|y-Ax\|_2^2}$ in the case of AWGN. Eq. 2.15 is therefore equivalent to:

$$ \underset{x}{\operatorname{argmin}} -\log\left(e^{-\frac{1}{2\sigma^2}\|y-Ax\|_2^2}\right) + \phi(x) \qquad (2.16) $$

We finally get the MAP estimation:

$$ \hat{x}_{MAP} \in \underset{x}{\operatorname{argmin}} \frac{1}{2} \|y - Ax\|_2^2 + \sigma^2 \phi(x). \qquad (2.17) $$

By combining the likelihood (i.e. data fidelity term) and the prior (i.e. regularization term), the MAP estimate effectively balances the trade-off between fitting the data and conforming to prior knowledge, leading to a reliable restoration despite the degraded observation. One key point is that the term $-\log p(x)$ should assume a closed-form formulation that facilitates feasible numerical optimization. Consequently, most attempts at describing $p(x)$ have opted for the Gibbs distribution representation [51], $p(x) = c.e^{-\phi(x)}$,
thereby redirecting our attention towards the energy function $\phi(x)$ rather than the prior distribution $p(x)$. It's worth noting that within the research community, there's a prevalent simplification where the term 'prior' is sometimes used interchangeably to refer to the regularization function (or even a denoiser which is also related to the prior distribution).

In most cases, Eq. 2.17 does not have a closed-form solution, therefore the estimation is computed using iterative optimization methods, as in the variational case. In Section 2.3, we discuss some optimization methods that can be used to solve common variational problems.

## 2.3 Optimization algorithms

This section aims to provide a brief overview of some of the optimization methods used to solve common variational problems. While numerous categories of methods exist, our focus here is primarily on commonly employed algorithms in the recent literature of inverse problems, notably first-order optimization methods and proximal methods.

### 2.3.1 Derivative-based optimization algorithms: first-order methods

Let's consider minimizing functions that are smooth and differentiable. In this case, derivative-based algorithms, such as gradient descent, are effective tools for solving such optimization problems.

#### 2.3.1.1 Gradient Descent

Let us consider minimizing the following optimization problem:

$$ \hat{x} = \operatorname*{argmin}_{x \in \mathbb{R}^n} f(x). \qquad (2.18) $$

If $f$ is a convex differentiable function with a gradient that is $\beta$-Lipschitz, then we can solve Eq. 2.18 by finding $\hat{x}$ such that:

$$ \nabla f(\hat{x}) = 0. \qquad (2.19) $$

Thus, we can easily solve 2.18 in the case where $\nabla f$ is sufficiently straightforward to
derive a closed-form solution from Eq. 2.19 which is computationally tractable. However, finding such a solution can become challenging when $f$ has a more complex form and we can no longer solve 2.19. In this case, the solution of the minimization can be defined as the limit of a sequence constructed iteratively through some algorithm. One category of methods to find an iterative solution of the minimization problem defined in Eq. 2.18 is the gradient-based algorithms [52].

Gradient methods are based on the principle that for any $x$, the negative gradient $-\nabla f(x)$ points in the direction of the steepest decrease, making it a favorable choice for finding the minimizer of $f$. Let $\gamma_k$ be the step-size of the gradient descent at iteration $k$ and $K \in \mathbb{N}$ be the total number of iterations until convergence. Starting with an initial point $x_0 \in \mathbb{R}^n$, we can derive $\hat{x}$ by applying the Algorithm 1:

**Algorithm 1 Gradient Descent**

1: Initialize $x_0 \in \mathbb{R}^n$

2: **for** $k \leftarrow 0$ to $K-1$ **do**

3: $x_{k+1} \leftarrow x_k - \gamma_k \nabla f(x_k)$

4: **end for**

### 2.3.1.2 Projected Gradient Descent

Let us now consider a constrained optimization problem, where the objective function needs to be minimized while satisfying a certain constraint.

$$ \hat{x} = \underset{x \in \mathbb{R}^n}{\operatorname{argmin}} f(x), \qquad (2.20) \\ \text{s.t. } g(x)=0 $$

where $f$ is the original objective function to be minimized and $g$ is a convex function defining the constraint or the feasible set of the solution.

Projected gradient descent consists in projecting the current iterate onto the feasible set at each iteration, to ensure that the constraints are satisfied.

**Algorithm 2 Projected Gradient Descent**

1: Initialize $x_0 \in \mathbb{R}^n$

2: **for** $k \leftarrow 0$ to $K-1$ **do**

3: $x_{k+\frac{1}{2}} \leftarrow x_k - \gamma_k \nabla f(x_k)$

4: $x_{k+1} \leftarrow \text{Proj}_{g(x)}(x_{k+\frac{1}{2}})$

5: **end for**
Algorithm 2 describes the method of projected gradient descent. The projection operator is defined as

$$ \mathrm{Proj}_{g(x)}(x_k) = \underset{\hat{x} \in g(x)}{\operatorname{argmin}} \|\hat{x} - x_k\|_2, \qquad (2.21) $$

where the feasible set is defined by the function g. The projection itself is an optimization problem to be solved.

## 2.3.2 Proximal algorithms

### Proximal Operator

Before diving into proximal algorithms, let us first define a proximal operator. Consider $V$ to be a real Hilbert space and $\Gamma_0(V)$ to represent the set of lower semi-continuous convex functions from $V$ to $[-\infty, +\infty]$. Let $h \in \Gamma_0(V)$ be a closed proper convex function. Jean Jacques Moreau [53] proposed that for each $x \in V$, the minimization on the right-hand side of Eq. 2.22 has a unique solution and is denoted by $\mathrm{prox}_{\lambda h}(u)$.

$$ \mathrm{prox}_{\lambda h}(u) = \underset{x \in V}{\operatorname{argmin}} h(x) + \frac{1}{2\lambda} \|x - u\|_2^2. \qquad (2.22) $$

The term $\mathrm{prox}_{\lambda h}(u)$ is called the proximal operator of the scaled function $\lambda h$, or the proximal operator of $h$ with parameter $\lambda$. Let us note that, in practice, if $h$ is not convex, there may be an abuse of notation, as Eq. 2.22 will have multiple solutions.

If $h$ is differentiable and $\lambda$ is small, the proximal operator of $\lambda h$ can also be evaluated by doing a gradient step on $h$:

$$ \mathrm{prox}_{\lambda h}(u) \approx u - \lambda \nabla h(u). \qquad (2.23) $$

Proximal algorithms are advantageous tools to solve non-smooth minimization problems. The base operation of a proximal algorithm is to compute the proximal operator of a function, which reduces to solve a simple convex optimization sub-problem. The latter often has a closed-form solution or can be easily solved with simple methods. Examples of proximal algorithms are presented next.
### 2.3.2.1 Proximal Gradient Descent

Let's consider a minimization optimization as follows:

$$ \hat{x} = \underset{x \in \mathbb{R}^n}{\operatorname{argmin}} f(x) + g(x), \qquad (2.24) $$

where $f$ is a smooth and differentiable function whereas $g$ is a non-smooth, possibly non-differentiable function that enforces some constraint.

Proximal gradient descent (also called Forward-Backward Splitting algorithm) consists in computing the gradient of the smooth part of the objective function at the current point $x_k$ (i.e. $\nabla f(x_k)$) and updating the latter, then evaluating the proximal operator at the updated point $x_k - \gamma_k \nabla f(x_k)$. The proximal operator is used to handle the non-smooth regularization term $g$. These 2 steps are performed repeatedly at each iteration until reaching convergence. Algorithm 3 details the proximal gradient descent algorithm.

#### Algorithm 3 Proximal Gradient Descent

1: Initialize $x_0 \in \mathbb{R}^n$

2: **for** $k \leftarrow 0$ to $K-1$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp; $x_{k+\frac{1}{2}} \leftarrow x_k - \gamma_k \nabla f(x_k)$

4: &nbsp;&nbsp;&nbsp;&nbsp; $x_{k+1} \leftarrow \text{Prox}_{\gamma_k g}(x_{k+\frac{1}{2}})$

5: **end for**

### 2.3.2.2 Half-Quadratic Splitting

Half-Quadratic Splitting (HQS) [54] is a simple yet effective splitting method for solving optimization problems. The basic idea behind HQS is to split the original optimization problem into two sub-problems, each of which addressing one part of the original problem. The algorithm first reformulates the problem by introducing an auxiliary variable to the Eq. 2.24, as follows:

$$ \begin{array}{ll} \underset{x,z \in \mathbb{R}^n}{\operatorname{argmin}} & f(x) + g(z). \qquad (2.25) \\ \text{subject to} & x = z \end{array} $$

The HQS iteration consists in alternatively applying the proximal operators of $f$ and $g$, as follows:
$$x^{k+1} = \operatorname*{argmin}_{x} f(x) + \frac{\rho^k}{2} \|x - z^k\|_2^2 = \operatorname*{prox}_{\frac{1}{\rho^k} f}(z^k), \qquad (2.26)$$

$$z^{k+1} = \operatorname*{argmin}_{z} g(z) + \frac{\rho^k}{2} \|z - x^{k+1}\|_2^2 = \operatorname*{prox}_{\frac{1}{\rho^k} g}(x^{k+1}), \qquad (2.27)$$

where $\rho^k$ is a positive penalty parameter for weighting the penalty terms. In order to enforce the constraint $x = z$, it is necessary to progressively increase the penalty parameter $\rho^k$ towards infinity within the HQS iterations.

### 2.3.2.3 Alternating Direction Method of Multipliers

Given the optimization problem given in Eq. 2.24, the Alternating Direction Method of Multipliers (ADMM) algorithm aims to solve the problem by first decoupling the data term and the prior term by adding an auxiliary variable $z$ just like the HQS, yielding a constrained optimization problem which is equivalent to 2.25.

In order to deal with the constraint $x = z$, the ADMM introduces a Lagrange multiplier: the original problem is reformulated into an augmented Lagrangian form, which includes the Lagrange multiplier (dual variable) associated with the constraints. The augmented Lagrangian of the optimization problem with the equality constraint $x = z$ is given by:

$$\begin{aligned} \mathcal{L}(x, z, l) &= f(x) + g(z) + l^T(x - z) + \frac{\rho}{2} \|x - z\|_2^2, \\ &= f(x) + g(z) + \frac{\rho}{2} \|x - z + \frac{l}{\rho}\|_2^2 - \frac{1}{2\rho} \|l\|_2^2, \end{aligned} \qquad (2.28)$$

where $l$ is the dual variable and $\rho$ is a positive penalty parameter.

Each ADMM iteration then consists in an alternate minimization over $x$ and $z$, as follows:

$$x^{k+1} = \arg \min_{x} \mathcal{L}(x, z^k, l^k) = \operatorname*{prox}_{\frac{1}{\rho^k} f}(z^k - \frac{l^k}{\rho^k}) \qquad (2.29)$$

$$z^{k+1} = \arg \min_{z} \mathcal{L}(x^{k+1}, z, l^k) = \operatorname*{prox}_{\frac{1}{\rho^k} g}(x^{k+1} + \frac{l^k}{\rho^k}) \qquad (2.30)$$

$$l^{k+1} = l^k + \rho^k (x^{k+1} - z^{k+1}), \qquad (2.31)$$

where the dual variable is typically zero-initialized. The penalty parameter $\rho^k$ can be
adjusted during the iterations to balance between convergence speed and accuracy. In Appendix A.1, we will take a closer look at the ADMM algorithm, providing detailed insights and exploring parameter settings in the context of the Plug-and-play framework.

## Closed-form solution of the data-fidelity optimization sub-problem

Both ADMM and HQS split the problem into two sub-problems. Hence, they can be used to solve the MAP estimation for linear inverse problems in Eq. 2.17, by taking $f(x) = \frac{1}{2} \|y - Ax\|_2^2$ and $g(x) = \sigma^2\phi(x)$. Therefore, the data-fidelity term $f$ and the prior term $g$ are handled as independent sub-problems that consist in evaluating the corresponding proximal operator $\text{prox}_{\frac{1}{\rho^k}f}$ (Eq. 2.26 and Eq. 2.29) and the regularization term $\text{prox}_{\frac{1}{\rho^k}g}$ (Eq. 2.27 and Eq. 2.30).

In particular, for the data fidelity term, the proximal operator is expressed as:

$$ \text{prox}_{\frac{1}{\rho^k}f}(u) = \underset{x}{\operatorname{argmin}} \mathcal{Q}(x, u) \qquad (2.32) $$

$$ \text{with } \mathcal{Q}(x, u) = \frac{1}{2} \|y - Ax\|_2^2 + \frac{\rho^k}{2} \|x - u\|_2^2. \qquad (2.33) $$

One can thus derive a general closed-form solution by finding $x$ that cancels out $\nabla_x \mathcal{Q}(x, u)$. This gives:

$$ \text{prox}_{\frac{1}{\rho^k}f}(u) = (A^T A + \rho^k I)^{-1}(A^T y + \rho^k u). \qquad (2.34) $$

In general, calculating the inverse directly might be computationally expensive, especially for large matrices. However, for many specific problems such as super-resolution, deblurring, and denoising, the matrix $A$ can often be represented in a particular form (e.g. diagonal matrix or block-diagonal matrix) which simplifies the computation and make the closed-form solution in Eq. 2.34 feasible.

## 2.4 Learning to Regularize inverse problems

The choice of the regularization function plays a crucial role in the efficient reconstruction of inverse problems. The advent of deep learning has remarkably transformed various domains, including the field of inverse problems and its regularization. In the last years, researchers have elaborated different ways to use learned priors, be it through generative models, deep denoisers or other approaches. In this section, we elaborate some
of these methodologies.

### 2.4.1 Regularizing by using the implicit prior captured by a neural network

A unique approach called the Deep Image Prior (DIP) was introduced by Ulyanov et al. [1], where it was shown that the CNN structure itself can inherently regularize image restoration without having to be trained on a specific dataset. DIP represents an untrained generative model, that takes a random input vector (i.e. noise) and adjusts its randomly initialized weights to minimize the mean square error (MSE) between its output and the degraded observation. The regularization typically employed in conventional methods is replaced by the implicit prior captured by the architecture of the generator network. The DIP formulation can be expressed as follows:

$$ \Theta^* \in \underset{\Theta}{\operatorname{argmin}} E(f_{\Theta}(z); x_0), \quad \text{s.t. } x^* = f_{\Theta^*}(z) \qquad (2.35) $$

where $f_{\Theta}$ is the generator network with random input $z$ and $\Theta$ its parameters to be optimized. $x_0$ is the degraded observation. The goal is to find the optimal set of parameters $\Theta$ that minimize the task-dependent data term $E(f_{\Theta}(z); x_0)$.

To effectively control the regularization in DIP, early stopping is commonly employed during the optimization process, which prevents the network from overfitting to the observed degradation. In practice, this early stopping is hard to control.

Subsequent works have aimed to enhance the performance of the Deep Image Prior by different methods. We will dedicate our next chapter to discuss these advancements and to study the effect of combining the DIP regularization with a deep regularization.

### 2.4.2 Regularizing by leveraging the power of denoisers

Image denoising has garnered significant research attention in the literature. Clarivate Web-Of-Science (WoS) identifies nearly 30,000 papers targeting image denoising published in the last 25 years, with an important increase in the number of publications over time. This growth in interest can be attributed not only to the crucial role of image denoising in pre-processing, but also to the recent revelation that denoisers can be exploited for broader applications, such as solving general inverse problems or synthesizing high-quality images. In the following, we discuss the ability to deploy denoising engines for regularizing inverse
problems.

### 2.4.2.1 Plug-and-play (PnP) framework

The Plug-and-play framework revolutionized traditional inverse problem-solving by seamlessly incorporating the capabilities of advanced denoisers. First, let's recall the equation of the proximal operator of a function $\phi$ with parameter $\lambda$:

$$ \text{prox}_{\lambda h}(u) = \underset{x \in \mathbb{R}^n}{\operatorname{argmin}} \frac{1}{2} \|x - u\|_2^2 + \lambda h(x). \quad (2.36) $$

Observing Eq. 2.17, we can state that the proximal operator in Eq. 2.36 can be seen as a particular case of inverse problem where the degradation operator is the identity matrix, and the degradation only consists in the addition of White Gaussian Noise of standard deviation $\sqrt{\lambda}$.

Inspired by this perspective, and motivated by the significant advancements in the field of image denoising, Venkatakrishnan et al. [55] introduced the Plug-and-play (PnP) framework, which stimulated an important sub-field of research in image restoration. The authors suggest to solve different types of inverse problems by replacing the proximal operator of the regularization term by a powerful denoiser in an ADMM framework.

In the context of MAP estimation, the constrained optimization can be formulated as:

$$ \begin{array}{ll} \hat{\mathbf{x}} = \underset{\mathbf{x},z}{\operatorname{argmin}} & \frac{1}{2} \|\mathbf{A}\mathbf{x} - \mathbf{y}\|_2^2 + \sigma^2 \phi(\mathbf{z}), \quad (2.37) \\ \text{subject to} & \mathbf{x} = \mathbf{z}. \end{array} $$

Each ADMM iteration then consists in an alternate minimization over $\mathbf{x}$ and $\mathbf{z}$ with the introduction of an additional variable $\mathbf{l}$ called dual variable as well as a penalty parameter $\rho$, as follows:

$$ \mathbf{x}^{k+1} = \underset{\mathbf{x}}{\operatorname{argmin}} \ \|\mathbf{A}\mathbf{x} - \mathbf{y}\|_2^2 + \rho^k \left\| \mathbf{x} - \left(\mathbf{z}^k - \frac{\mathbf{l}^k}{\rho^k}\right) \right\|_2^2, \quad (2.38) $$

$$ \mathbf{z}^{k+1} = \underset{\mathbf{z}}{\operatorname{argmin}} \frac{1}{2} \left\| \mathbf{z} - \left( \mathbf{x}^{k+1} + \frac{\mathbf{l}^k}{\rho^k} \right) \right\|_2^2 + \frac{\sigma^2}{\rho^k} \phi(\mathbf{z}) = \underset{\mathbf{z}}{\operatorname{prox}}_{\frac{\sigma^2}{\rho^k} \phi} (\mathbf{x}^{k+1} + \mathbf{l}^k / \rho), \quad (2.39) $$

$$ \mathbf{l}^{k+1} = \mathbf{l}^k + \rho^k (\mathbf{x}^{k+1} - \mathbf{z}^{k+1}), \quad (2.40) $$
where the dual variable is typically zero-initialized.

The Plug-and-play ADMM thus replaces the z-update in 2.39 by the following denoising step:

$$z^{k+1} = D_s(x^{k+1} + l^k/\rho), \qquad (2.41)$$

where $D$ is a denoiser suitable for Gaussian noise of standard deviation $s$, with $s = \sqrt{\sigma^2/\rho} = \sigma/\sqrt{\rho}$.

While early works used traditional denoising methods such as BM3D [55], [56] or non-local means [57], recent research extended the PnP framework to powerful deep CNNs such as DnCNN [43] or DRUNet [2]. The latter leads to a State-Of-the-Art image restoration technique. Also, the Plug-and-play framework was applied using different proximal algorithms such as HQS [2], [3] and Proximal Gradient Descent [58].

An advantage of PnP methods is that, unlike agnostic learning approaches where the training takes place for each inverse problem separately, integrating deep denoisers in PnP frameworks decouples the learning phase from the inverse problem in hand, which makes it generic.

However, a limitation of proximal PnP algorithms lies in the fact that the regularizer is only defined implicitly via the denoiser, meaning that the objective function is not clearly defined. This limits the interpretability of the algorithm, since it becomes challenging to keep tractability of the optimization problem. Moreover, for most denoisers, there is no theoretical guarantee that there exists a potential regularization function whose proximal operator can be accurately represented by the denoiser [59]. This induces convergence issues in proximal PnP algorithms.

Although still an open question, designing convergence proofs for the Plug-and-play framework has become a major focus in the recent literature [3], [56], [58], [60]–[64]. The theoretical convergence analysis of these algorithms being outside the scope of this thesis, we refer interested readers to the review papers [65], [66]. In summary, most of the proposed solutions come with certain limitations, such as imposing strong assumptions on the denoiser (often not suitable for deep CNNs) or necessitating performance compromises.

### 2.4.2.2 Regularization by Denoising (RED)

In order to overcome some limitations arising from the absence of a clear objective function, Romano et al. [4] introduced an alternative approach for denosing-based priors called Regularization by Denoising (RED). Unlike some other methods where the prior
is implicitly defined, RED takes a different approach by explicitly expressing the regularization function $\phi$. The latter is defined to be proportional to the inner product between the estimated image and its denoising residual, which can be written as follows:

$$ \phi(x) = \frac{1}{2} x^T [x - D(x)]. \qquad (2.42) $$

where $D$ is an off-the-shelf denoiser. This explicit formulation results in an image-adaptive Laplacian-based regularization term, effectively promoting smoothness in the restored image. Under the assumption that the denoiser $D$ is locally homogeneous and its Jacobian is symmetric, the gradient of the regularization term defined in Eq. 2.42 is computed as:

$$ \nabla \phi(x) = x - D(x). \qquad (2.43) $$

Using this formulation, RED proposes to solve Eq. 2.17 by using different optimization algorithms such as gradient descent or ADMM. Unlike Plug-and-play methods, RED is not restricted by the choice of the algorithm, making it adaptable to different optimization strategies. Another key advantage of the RED algorithm is its explicit formulation of the regularization function $\phi(x)$ as given in Eq. 2.42. This explicit expression allows for better understanding of the overall Bayesian objective function.

However, Reehorst and Schniter [59] proved later that the gradient expression proposed in RED is not justified with denoisers that lack Jacobian symmetry, which excludes most practical denoisers such as block-matching and 3D filtering (BM3D) [67], Non-local means (NLM) [68] or state-of-the-art deep neural networks.

Another extension to the RED algorithm proposed by Cohen et al. [69] explores a projected RED via Fixed-Point Projection approach (RED-PRO). Rather than explicitly defining the regularization term, the fixed-point set of a demi-contractive denoiser $D\{x \in \mathbb{R}^n: x = D(x)\}$ is considered as a prior for the image restoration. However, in practice, verifying the assumption that the denoiser is demi-contractive is not straightforward [70].

### 2.4.2.3 The score function in the context of Plug-and-play MAP estimators

A different approach leading to the same formulation as RED uses Tweedie's identity to match the score function with an MMSE denoiser.
Starting from the MAP estimation for solving an inverse problem, we have:

$$ \hat{x}_{MAP} = \underset{x}{\operatorname{argmin}} \frac{1}{2} \|y - Ax\|_2^2 - \sigma^2 \log(p(x)). \quad (2.44) $$

Solving Eq. 2.44 by Steepest Descent (SD) amounts to solving the iterative formula:

$$ \hat{x}_{k+1} = \hat{x}_k - \gamma_k [A^T(Ax_k - y) - \sigma^2 \nabla_x \log(p(x))|_{\hat{x}_k}], \quad (2.45) $$

where $\nabla_x \log(p(x))$ is defined as the score function in statistical literature. Moreover, a well-known mathematical finding commonly credited to Tweedie [71], [72], Stein [73] or Miyasawa [74] outlines:

$$ \nabla_y \log(p_{\sigma_0}(y)) = \frac{D(y, \sigma_0) - y}{\sigma_0^2}, \quad (2.46) $$

where $y = x + v$ is a degraded version of $x$, corrupted with AWGN $v \sim \mathcal{N}(0, \sigma_0^2 I)$, and $D(y, \sigma_0)$ is the optimal MMSE denoiser for a Gaussian noise of standard deviation $\sigma_0$. The proof of this result is relatively straightforward and can be found in [7].

However, the relevant score function is that of the noisy observation density $p_{\sigma_0}(y)$, rather than the one corresponding to the prior density desired in Eq. 2.45. The noisy observation density $p_{\sigma_0}(y)$ is the result of a convolution between the prior density $p(x)$ and the noise distribution: $p_{\sigma_0}(y) = p(x) \otimes \mathcal{N}(0, \sigma_0^2 I)$. It is assumed, that when $\sigma_0$ is sufficiently small, $p(x)$ can be approximated by a slightly blurry PDF, and the update in Eq. 2.45 becomes:

$$ \hat{x}_{k+1} = \hat{x}_k - \gamma_k [A^T(Ax_k - y) + c(\hat{x}_k - D(\hat{x}_k, \sigma_0)]. \quad (2.47) $$

Here, $c = \frac{\sigma_0^2}{\sigma_0^2}$, and this is equivalent to the RED algorithm when solved with Steepest Descent.

In practice, we don't know the real prior distribution of natural images, hence we cannot obtain the MMSE denoiser. However, deep denoisers can often approximate the MMSE denoisers very well, given that they are typically learned using an MSE loss on a large set of natural images.
### 2.4.2.4 Posterior sampling: The score function in the context of MMSE estimators

Different from the aforementioned approaches, we elaborate here on the computation of MMSE estimators. We recall that the MMSE estimator of the posterior mean is given by:

$$ \hat{x}_{MMSE} = \int xp(x|y)dx. \qquad (2.48) $$

This computation can be challenging, as it consists in estimating high dimensional integrals. Therefore, sampling algorithms such as the Monte Carlo integration method can be used in order to perform numerical integration. Posterior sampling has been considered in [5]–[8].

When solving 2.48, sampling algorithms consist in taking different samples from the posterior distribution. As the number of samples increases, the average of those samples approaches the expected value of the random variable. In the context of integration, this means that as we consider more samples, the average of the function values over those samples converges to the true integral value.

Posterior sampling approaches are largely based on the annealed Langevin dynamics algorithm [6], a Monte Carlo method with the subsequent transition rule:

$$ \begin{aligned} x_{k+1} &= x_k + \gamma \nabla \log p(x_k|y) + \sqrt{2\gamma} z_k \\ &= x_k + \gamma \nabla \log p(y|x_k) + \gamma \nabla \log p(x_k) + \sqrt{2\gamma} z_k, \end{aligned} \qquad (2.49) $$

where $z_k \sim N(0, I)$ and $\gamma > 0$ is a step-size. The term $\nabla \log p(y|x_k)$ reflects the data fidelity and $\nabla \log p(x_k)$ is the score function which can be approximated by an MMSE denoiser as discussed in Section 2.4.2.3. Finally, Eq. 2.49 is an iterative algorithm that converges towards a solution that reflects the posterior distribution.

### 2.4.3 Regularizing through generative models

Another category of methods define the image prior by generative models such as a Generative Adversarial Network (GAN) [75] or a Variational Audoencoder (VAE) [76]. A generative model $G: \mathbb{R}^k \to \mathbb{R}^n$ learns to take a low-dimensional latent variable $z \in \mathbb{R}^k$ and map it to a high dimensional sample space $x \in \mathbb{R}^n$ such that $x = G(z)$. The idea is to solve the MAP estimator on the latent variable $z$ and restrict the output to be in the
range of the generator:

$$ \hat{z} = \underset{z}{\operatorname{argmin}} -\log p(y|G(z)) - \log p(z), \qquad (2.50) $$

The reconstructed image is obtained by $\hat{x} = G(\hat{z})$. This approach was first introduced by Bora et al. [77], motivated by the fact that well-trained generative models approximate the space of natural signals or images.

Subsequently, a lot of research has been done in this direction, resulting in using generative models for regularizing inverse problems either implicitly by restricting the solution in the range of a generative model, or explicitly by learning a CNN regularizer using generative models.

Shah et al. [78] later proposed to solve the approach proposed by [77] with a projected gradient descent algorithm and discussed convergence guarantees.

Furthermore, rather than restricting $\hat{x}$ to fall within the range of the generator and minimizing on the latent variable $z$, Gonzalez et al. [79] optimize the joint posterior $p(x, z|y)$ of the space variable $x$ and the latent variable $z$, given the observation $y$. They evaluate the joint $MAP_{x,z}$ estimator (JPMAP) using a VAE prior. The authors also emphasize the importance of training the VAE with a denoising criterion (denoising VAE) proposed by [80], which results in a more robust encoder that can generalize well for inputs that do not lie in the manifold of the generator network. However, a limitation of JPMAP lies in the fact that it is tailored for VAEs with a fixed Gaussian prior distribution over the latent space, which confines its applicability to toy examples.

In a subsequent work [81], Prost et al. generalize the JPMAP framework to regularize inverse problems with Hierarchical Variational Autoencoders (HVAE). In this context, the regularization strength can be controlled using the temperature of the prior in the latent space (the temperature refers to a hyper-parameter that controls the smoothness of the generated samples). This approach also overcomes the limitation of the previous JPMAP [79] to simple VAEs.

### 2.4.4 End-to-end learning of the regularization function via deep unrolling

Introduced by Gregor et al. [82], the idea of deep unrolling quickly grabbed the interest of researchers and made significant progress in subsequent works [83]–[88]. In the context of deep unrolling, often referred to as algorithm unrolling, the regularizing function is
learned end-to-end for a specific inverse problem. The idea is to iterate an optimization algorithm for a pre-determined number of iterations, and to replace the regularization function with a CNN to be learned specifically for a particular problem. The learned network replaces the gradient of the regularizer in the case of a gradient-based algorithm, or its proximal operator if the iterative algorithm is a proximal algorithm. These networks are often referred to as Artifact Removal operators trained to remove artifacts specific to the forward model.

For a more comprehensive understanding of deep unrolling, we illustrate the case of the algorithm of unrolled gradient descent. Solving Eq. 2.17 with the gradient descent algorithm amounts to the iterative process described by:

$$x_{k+1} = x_k - \gamma[A^T(Ax_k - y) + \nabla\phi(x_k)], \quad (2.51)$$

where $\gamma > 0$ denotes the step size of the optimization algorithm. In the context of the unrolled gradient descent, this process is carried out for a specified number of iterations (referred to as blocks), which we will denote by $B$. The latter remains constant throughout both the training and testing phases.

<figure><img src="image_10.png" /><figcaption>Figure 2.7 – Unrolled gradient descent framework. Here, $\nabla\phi$ is a trained neural network. The resulting output $\hat{x}$ is obtained after $B$ iterations of gradient descent with a fixed step size of $\gamma$.</figcaption></figure>

Figure 2.7 provides an illustrative representation. Here, the neural network parameterizes the gradient of the regularizer $\nabla\phi$, which is learned in conjunction with the forward model $A$ and training observations $y_i$. The learning process involves minimizing a loss function quantifying the discrepancy between the ground truth images $x_i$ and the output of the full network $\hat{x}(y_i)$. Therefore, this end-to-end learning approach tailors the learned component specifically to the underlying forward model.

As a result of the prior being specifically learned for the inverse problem in hand, deep unrolling methods yield significantly improved performance compared to Plug-and-play methods where the prior is generic. However, these end-to-end approaches lose the genericity of Plug-and-play algorithms. Furthermore, it is not clear what interpretation
can be given to the trained network which does not only learn prior knowledge on images,
but also task-specific features.
<figure><img src="image_11.png" /></figure>
# REGULARIZATION OF THE DEEP IMAGE
PRIOR WITH A DEEP DENOISER

**3.1 Introduction**

As deep learning evolved in the past decade, the utilization of end-to-end training methods brought about a significant breakthrough in solving inverse problems. Neural networks are typically used to learn to reconstruct degraded images by using a large amount of data, without any knowledge of the degradation operator.

As previously discussed in Section 2.4.1, an exceptional work by Ulyanov et al. [1] proposed to leverage the power of deep neural networks in a different way. The authors present the Deep Image Prior (DIP), where they show that the architecture of a generator network itself can be used as a regularization for inverse problems. DIP is an un-trained generative model that fits its parameters for the degraded image. Therefore, the prior is essentially hand-crafted, as no component of the model is learned beforehand from data.

Despite its simple formulation, the Deep Image Prior has demonstrated successful performance in solving various problems, including denoising, super-resolution, and inpainting. However, its results still lag behind those achieved by other state-of-the-art methods. This motivates to enhance the DIP by integrating explicit regularization to boost the implicit one.

Interestingly, prior research has explored the regularization of the DIP, yet these studies exclusively employed handcrafted regularization techniques, without exploiting the potential of learned regularizers.

In this chapter, we delve into the regularization of the estimate produced by the untrained generative DIP. We elaborate a method coupling a learned denoiser with the Deep Image Prior in the ADMM framework. The goal is to couple advantages of the DIP only optimized on the input image, with those of a generic prior learned from a large collection of natural images.
## 3.2 Related work: regularized Deep Image Prior

In this section, we elaborate details about the Deep Image Prior and discuss some regularization methods that have been applied to boost its performance in the literature. The deep generative model $x = f_{\Theta}(z)$ is learned by mapping a random code vector $z \in \mathbb{R}^N$ to an image x, i.e., by solving

$$ \Theta^* \in \underset{\Theta}{\operatorname{argmin}} \quad \|Af_{\Theta}(z) - y\|_2^2 \quad \text{s.t. } x^* = f_{\Theta^*}(z) \qquad (3.1) $$

where $\Theta$ represents the network parameters. The generator is randomly initialized with variables $\Theta$, which are optimized iteratively in a way that the output of the network is close to the target measurements. In most of the applications, a U-Net type architecture with skip-connections is used, having over 2 million parameters (Figure 3.1). Details of the network are adapted to the specific application.

<figure><img src="image_12.png" /><figcaption>Figure 3.1 – Architecture of the Deep Image Prior. An "hourglass" architecture is used. $n_u[i]$, $n_d[i]$, $n_s[i]$ represent the number of feature maps at the $i^{th}$ layer for the upsampling, downsampling and skip-connections respectively and $k_u[i]$, $k_d[i]$, $k_s[i]$ denote the kernel sizes. Skip connections and parameters are defined depending on the application. The figure is extracted from [1].</figcaption></figure>
## TV-regularized Deep Image Prior

Liu et al. [89] propose to use a total-variation regularization to enhance the performance of the Deep Image Prior. They consider the variational problem defined as:

$$ \Theta^* \in \underset{\Theta}{\operatorname{argmin}} \quad \|Af_{\Theta}(z) - y\|_2^2 + \lambda_{TV}(f_{\Theta}(z)) \quad \text{s.t. } x^* = f_{\Theta^*}(z), \qquad (3.2) $$

and define an $l_1$-based anisotropic TV given by $TV_{aniso}(u) = \sum_{i=1}^n |(D_h u)_i| + |(D_v u)_i|$, $D_h u$ and $D_v u$ being the first order finite difference discrete operators respectively along the horizontal and vertical axes. The objective function is minimized using a gradient descent optimization method.

Subsequently, Cascarano et al. [90] propose to solve this objective function with the ADMM. Later, the authors propose to add a weight on the TV regularization in [91]. Another difference between the DIP-TV proposed in [89] and the ADMM-DIP-TV [90], [91] works is that in DIP-TV, anisotropic TV is used, decoupling the contribution of both horizontal and vertical gradient components, whereas ADMM-DIP-TV uses isotropic TV, jointly considering the gradient components.

## Denoiser-based DIP regularization

Another approach to boost the performance of the DIP is proposed in [92], bringing-in the concept of Regularization by Denoising (RED) [4], which uses existing denoisers for regularizing inverse problems. By merging the DIP with the RED, the objective function of the resulting DeepRED becomes:

$$ \underset{\Theta}{\operatorname{argmin}} \quad \frac{1}{2} \|Af_{\Theta}(z) - y\|_2^2 + \frac{\lambda}{2} [f_{\Theta}(z)]^T ([f_{\Theta}(z)] - D(f_{\Theta}(z))), \qquad (3.3) $$

where $D$ is a denoiser applied on the output of the generative network, that they replace with the BM3D [67]. Here, the optimization is solved using ADMM.

## 3.3 Regularizing the Deep Image Prior with a learned denoiser

We propose to regularize the estimate produced by the DIP with a learned denoiser. With this combination, we aim to combine the benefits of a generic prior learned on a
large set of natural images and those of the Deep Image Prior optimized only on the input image, where the network architecture itself plays the role of the regularization.

Let us consider the problem of regularizing the DIP optimization:

$$ 
\begin{aligned}
\Theta^* \in \operatorname*{argmin}_{\Theta} & \frac{1}{2} \|Af_{\Theta}(z) - y\|_2^2 + \sigma^2 \phi(x) \\
\text{s.t.} \quad & x^* = f_{\Theta^*}(z),
\end{aligned} 
\quad (3.4) $$

where we want to replace $\phi$ with a learned regularizer. For simplicity, we will replace $f_{\Theta}(z)$ by $t$ wherever convenient. The augmented lagrangian form of the Eq. 3.4 can be written as:

$$ L(x, t, l) = \frac{1}{2} \|At - y\|_2^2 + \sigma^2 \phi(x) + \frac{\rho}{2} \left\| x - t + \frac{l}{\rho} \right\|_2^2 - \frac{1}{2\rho} \|l\|_2^2, \quad (3.5) $$

with $\rho$ being a positive penalty parameter of the constraint, and $l$ the dual variable. This minimization problem is solved using the iterative ADMM method. By alternatively optimizing $x$, $t$ and $u$ of the augmented Lagrangian $L(x, t, l)$, the ADMM iterate reads as:

$$ t^{k+1} = f_{\Theta^{k+1}}(z) \quad \text{with} $$

$$ \Theta^{k+1} = \operatorname*{argmin}_{\Theta} \frac{1}{2} \|Af_{\Theta}(z) - y\|_2^2 + \rho \left\| x^k - f_{\Theta}(z) + \frac{l^k}{\rho} \right\|_2^2, \quad (3.6) $$

$$ x^{k+1} = \operatorname*{argmin}_{x} \frac{1}{2} \left\| x - t^{k+1} + \frac{l^k}{\rho} \right\|_2^2 + \frac{\sigma^2}{\rho} \phi(x), \quad (3.7) $$

$$ l^{k+1} = l^k + \rho(x^{k+1} - t^{k+1}). \quad (3.8) $$

As discussed in Section 2.4.2.1, the $x$-update can be rewritten as $x^{k+1} = \operatorname*{prox}_{\frac{\sigma^2}{\rho}\phi}(t^{k+1} - u^k)$ and can be replaced by a denoiser. Hence, inspired by [55], we propose to replace the $x$-update by a learned denoiser [2]. We thus get:

$$ x^{k+1} = D_s(t^{k+1} - u^k), \quad (3.9) $$

with $D_s$ being a learned denoiser and $s = \frac{\sigma}{\sqrt{\rho}}$. We use the DRUNet [2], a state-of-the-art deep denoiser, which takes as input the noisy image concatenated in the channel dimension with a noise level map. It is worth noting that since ADMM separates the regularization
<figure><img src="image_13.png" /><figcaption>(a) Zebra</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) House</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Lighthouse</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Cameraman</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) Butterfly</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(f) Hill</figcaption></figure>

<caption>Figure 3.2 – The set of images used for the numerical experiments.</caption>

term $\phi$ from the degradation matrix $A$, the learned proximal operator can be used with any linear operator.

In practice, Eq. 3.6 can't be solved exactly due to the non-linearity of $f_{\Theta}(z)$. Furthermore, the DIP should not overfit the exact solution in order to provide an additional regularization effect in complement to the denoiser. Therefore, we solve Eq. 3.6 in exactly by performing a fixed number $n_{DIP}$ of gradient descent iterations within the ADMM iterations. For the initialization, we optimize the DIP using the original method [1] which performs gradient descent to optimize Eq. 3.1. A number $n_{init}$ of iterations is used for this initialization, followed by $n_{ADMM}$ iterations of the ADMM scheme.

We note that this approach can also be seen as a sort of a Plug-and-play ADMM, where the data-term is given by the generative Deep Image Prior.

**3.4 Experimental results**

We have reproduced the results of DIP [1], DIPTV [89], ADMM-DIPTV [90] and DeepRED [92] for the problems of denoising and super-resolution and we have compared them with our proposed method. In order to evaluate the input of the DIP in the PnP-ADMM algorithm, we have also compared the results with the PnP-ADMM with the denoiser of [2] for super-resolution, and with the denoiser for the inverse problem of denoising.

The test dataset (Figure 3.2) was formed by taking the 5 natural images from the
ADMM-DIPTV paper [90], and the zebra image used in [1], which contains interesting high frequency patterns.

In [90], the authors use isotropic TV and analyze the addition of ADMM by comparing with anisotropic TV in [89]. For the sake of a fair comparison between [89] and [90], we have reproduced the anisotropic version of both of the methods. For DeepRED, we have used the code provided by the authors which applies the BM3D denoiser to regularize the DIP optimization.

For fair comparisons, we have used the same network architecture for the DIP in all the compared methods.

<table><caption>Table 3.1 – Parameters used for the proposed method for (i) Denoising with noise standard deviations of ($\sigma$ = 20 and $\sigma$ = 30), (ii) Super-Resolution of factor 4 and 8. $n_{init}$: DIP iterations at initialisation, $n_{ADMM}$: ADMM iterations, $n_{DIP}$: number of iterations of Eq. 3.6 in each iteration of ADMM, $s$: noise level assumed by the learned denoiser in Eq. 3.9. The same step size is used for all gradient descent steps.</caption><thead><tr><th rowspan="2"></th><th colspan="2">(i) Denoising</th><th colspan="2">(ii) Super-resolution</th></tr><tr><th>σ = 20</th><th>σ = 30</th><th>Factor 4</th><th>Factor 8</th></tr></thead><tbody><tr><td>n<sub>init</sub></td><td>1750</td><td>1750</td><td>2500</td><td>4000</td></tr><tr><td>n<sub>ADMM</sub></td><td>100</td><td>100</td><td>100</td><td>100</td></tr><tr><td>n<sub>DIP</sub></td><td>250</td><td>250</td><td>200</td><td>150</td></tr><tr><td>s</td><td>20</td><td>30</td><td>41</td><td>45</td></tr><tr><td>Step size</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr></tbody></table>

Table 3.1 shows the parameters used for the proposed method for each of the cases of denoising and super-resolution. Tables 3.2 and 3.3 show the PSNR [dB] results for each of the compared methods for denoising and super-resolution respectively. The PSNR measures presented in this thesis are computed on the RGB channels in the following way:

$$ \text{PSNR} = 10 \cdot \log_{10} \left( \frac{\text{MAX}_{I}^{2}}{\text{MSE}} \right), \qquad (3.10) $$

where $MAX_I^2$ is the maximum possible pixel value in the image (i.e. 1 or 255 depending whether the image is normalized or not), and MSE is the Mean Squared Error calculated as:

$$ \text{MSE} = \frac{1}{W \cdot H} \sum_{x=0}^{W-1} \sum_{y=0}^{H-1} ((\hat{I}(x, y) - I^{gt}(x, y))^2 \qquad (3.11) $$
<table><caption>Table 3.2 – PSNR [dB] of denoised images obtained with TV, DIP[1], DIPTV[89], ADMM-DIPTV[90], the denoiser of [2] and our proposed method, for noise standard deviations of (i) 20 and (ii) 30.</caption><thead><tr><th></th><th></th><th>DIP</th><th>DIPTV</th><th>ADMM-DIPTV</th><th>denoiser</th><th>DIP-denoiser-ADMM</th></tr></thead><tbody><tr><td rowspan="6">(i)</td><td>House</td><td>27.93</td><td>27.81</td><td>27.96</td><td>31.55</td><td>31.35</td></tr><tr><td>Zebra</td><td>29.98</td><td>29.63</td><td>29.5</td><td>31.21</td><td>31.1</td></tr><tr><td>Lighthouse</td><td>29.03</td><td>29.55</td><td>29.3</td><td>32.18</td><td>31.96</td></tr><tr><td>Hill</td><td>31.29</td><td>32.53</td><td>32.07</td><td>34.02</td><td>33.93</td></tr><tr><td>Butterfly</td><td>30.07</td><td>30.52</td><td>30.62</td><td>33.23</td><td>32.96</td></tr><tr><td>Cameraman</td><td>30.01</td><td>31.4</td><td>31.04</td><td>33.7</td><td>33.67</td></tr><tr><td></td><td>Average</td><td>29.72</td><td>30.24</td><td>30.08</td><td>32.64</td><td>32.49</td></tr><tr><td rowspan="6">(ii)</td><td>House</td><td>25.75</td><td>26.18</td><td>25.72</td><td>29.28</td><td>29.15</td></tr><tr><td>Zebra</td><td>27.34</td><td>27.52</td><td>27.58</td><td>28.93</td><td>28.87</td></tr><tr><td>Lighthouse</td><td>26.7</td><td>27.48</td><td>27.2</td><td>29.96</td><td>29.75</td></tr><tr><td>Hill</td><td>29.61</td><td>30.81</td><td>30.28</td><td>32.62</td><td>32.51</td></tr><tr><td>Butterfly</td><td>28.16</td><td>28.48</td><td>28.52</td><td>30.85</td><td>30.83</td></tr><tr><td>Cameraman</td><td>27.79</td><td>28.88</td><td>28.41</td><td>31.15</td><td>31.14</td></tr><tr><td></td><td>Average</td><td>27.55</td><td>28.23</td><td>27.95</td><td>30.46</td><td>30.37</td></tr></tbody></table>

where *H* and *W* represent the number of rows and columns in the image respectively and $\hat{I}(x, y)$ and $I^{gt}(x, y)$ represent the pixel value of the reconstructed and original images at position $(x, y)$ respectively.

For denoising (Table 3.2), we can see that, our proposed method significantly improves the result of the DIP. Hence, as expected, learned regularizers outperform handcrafted regularizers when regularizing the DIP optimization. However, the learned denoiser itself performs better alone, and the DIP does not have an added value in this case. This was expected in the case of denoising, since the network is learned end-to-end for the denoising task.

However, for the task of super resolution (Table 3.3), we have a different outcome. In fact, we can see again that regularizing with a learned denoiser over the DIP optimization gives a better performance than using handcrafted priors. But also, as opposed to what we had for the denoising, combining the DIP with the denoiser improves the performance of the learned regularizer alone. This means that combining a prior learned on a large set of natural images with the untrained generative DIP is advantageous in this case.

The degradation filter we have used for the task of super resolution is a Gaussian filter of standard deviation $\sigma_f = 1$ for both x4 and x8 super-resolution factors. Visual comparisons are shown in Figures 3.3 and 3.4 for super-resolution of factor 4. More visual results
<table><thead><tr><th></th><th></th><th>Bicubic</th><th>DIP</th><th>DIPTV</th><th>ADMM-DIPTV</th><th>DeepRED</th><th>Denoiser-ADMM</th><th>DIP-denoiser-ADMM</th></tr></thead><tbody><tr><td rowspan="6">(i) factor 4</td><td>House</td><td>19.42</td><td>19.72</td><td>19.85</td><td>19.62</td><td>19.88</td><td>20.06</td><td>20.14</td></tr><tr><td>Zebra</td><td>23.00</td><td>24.13</td><td>24.31</td><td>23.84</td><td>24.69</td><td>24.71</td><td>24.9</td></tr><tr><td>Lighthouse</td><td>23.09</td><td>23.05</td><td>23.37</td><td>23.17</td><td>23.36</td><td>23.46</td><td>23.45</td></tr><tr><td>Hill</td><td>26.37</td><td>28.44</td><td>27.61</td><td>27.01</td><td>27.81</td><td>28.75</td><td>28.97</td></tr><tr><td>Butterfly</td><td>20.95</td><td>24.23</td><td>23.66</td><td>22.47</td><td>23.79</td><td>24.34</td><td>24.42</td></tr><tr><td>Cameraman</td><td>22.58</td><td>23.23</td><td>23.27</td><td>22.99</td><td>23.34</td><td>23.34</td><td>23.51</td></tr><tr><td></td><td>Average</td><td>22.58</td><td>23.80</td><td>23.67</td><td>23.18</td><td>23.81</td><td>24.11</td><td>24.23</td></tr><tr><td rowspan="6">(ii) factor 8</td><td>House</td><td>16.25</td><td>17.16</td><td>17.32</td><td>17.11</td><td>17.22</td><td>16.83</td><td>17.21</td></tr><tr><td>Zebra</td><td>16.84</td><td>18.01</td><td>18.15</td><td>17.75</td><td>18.05</td><td>17.37</td><td>17.84</td></tr><tr><td>Lighthouse</td><td>19.87</td><td>20.6</td><td>20.65</td><td>20.45</td><td>20.65</td><td>20.47</td><td>20.57</td></tr><tr><td>Hill</td><td>20.75</td><td>22.45</td><td>22.1</td><td>22.01</td><td>22.43</td><td>23.46</td><td>23.58</td></tr><tr><td>Butterfly</td><td>15.38</td><td>17.04</td><td>16.62</td><td>16.63</td><td>16.80</td><td>17.11</td><td>17.4</td></tr><tr><td>Cameraman</td><td>18.94</td><td>20.2</td><td>20.1</td><td>19.95</td><td>20.15</td><td>19.88</td><td>20.19</td></tr><tr><td></td><td>Average</td><td>18.00</td><td>19.24</td><td>19.15</td><td>18.98</td><td>19.22</td><td>19.18</td><td>19.46</td></tr></tbody></table>

Table 3.3 – PSNR [dB] of upsampled images obtained with bicubic interpolation, DIP[1], DIPTV[89], ADMM-DIPTV[90], DeepRED[92], PnP-ADMM with the denoiser of [2] and our proposed method, for (i) SR factor 4 and (ii) SR factor 8.

can be found on the web page: http://clim.inria.fr/DeepCIM/Regul-DIP/index.html.

## 3.5 Conclusion and discussion

In this chapter, we studied different approaches of regularizing the generative Deep Image Prior. We proposed to regularize the DIP with a learned denoiser, in order to combine the benefits of a generic prior learned on a large set of natural images and the regularization provided by the network architecture of the DIP. The proposed denoiser-based regularization compares favourably with handcrafted priors and can be applied to different inverse problems. Moreover, the use of DIP improves the performance of the PnP-ADMM in the case of super-resolution.

However, several aspects were observed while we were working on this contribution. The parameter tuning of the ADMM algorithm is not a straightforward task. Also, employing a denoiser adds the complexity of fine-tuning an extra hyper-parameter for each application, while theoretically, a regularizer should be able to determine the image prior regardless of the inverse problem in hand. Our experimental implementations led us to the reflection that solving the DIP optimization could have been significantly simplified
<figure><img src="image_19.png" /><figcaption>Figure 3.3 – Visual comparison of super-resolution results on the Butterfly image. Low resolution image is generated with a Gaussian kernel of standard deviation $\sigma_f = 1$ followed by a downsampling by a factor of 4.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 3.4 – Visual comparison of super-resolution results on the Cameraman image. Low resolution image is generated with a Gaussian kernel of standard deviation $\sigma_f = 1$ followed by a downsampling by a factor of 4.</figcaption></figure>
by employing a network that directly models the regularizer within a gradient descent
algorithm. In the upcoming chapter, we present a methodology to learn a network mod-
eling the gradient of a regularizer, by using a deep denoiser. Given that a denoiser can be
seen as an implicit prior, we aim to learn a network that models this prior. The resulting
network enables solving inverse problems using the prior embedded in a denoiser with a
simple gradient descent algorithm, avoiding the complexity of the parameter tuning in
proximal algorithms.
# PnP-REG: LEARNED REGULARIZING GRADIENT FOR PLUG-AND-PLAY GRADIENT DESCENT

## 4.1 Introduction

In this chapter, we present a novel approach to regularize inverse problems using deep neural networks. We propose a framework to train a network that models the gradient of a regularizer, which can be used in a Plug-and-play gradient descent algorithm to solve any inverse problem.

The Plug-and-play (PnP) framework, introduced in [55], makes it possible to integrate advanced image denoising priors into optimization algorithms, to efficiently solve a variety of image restoration tasks generally formulated as Maximum A Posteriori (MAP) estimation problems. However, as already discussed in Chapters 2 and 3, PnP methods remain limited to proximal algorithms, which makes the task of parameter tuning quite challenging. This introduces complexity and reduces the algorithm's comprehensibility and ease of interpretation. Subsequent work, such as RED [4] and PnP-SGD [64], used an explicit regularization to circumvent this issue. These methods use the denoising residual as the gradient of the regularization and are not always restricted to proximal algorithms. A recent method called GS-PnP [3] plugs a denoiser trained to perform an exact gradient step on a regularization function represented by a CNN, leading to convergence guarantees. Another line of work considers score matching methods that solve inverse problems in the MMSE framework. However, these methods approximate the true signal density by using a MMSE denoiser's residual which is proportional to the gradient of the log of the noisy signal density.

A common issue shared by all these methods is that the regularization term is (implicitly or explicitly) defined by a denoiser. This involves an additional noise level hyper-
parameter that must be tuned per-application. However, in theory, a regularizer (and thus its gradient) should fully determine the image prior, regardless of the task to be solved.

We propose a novel approach to train a network modeling the gradient of a regularizer. Our method makes use of a second network pre-trained for the denoising task. Our goal is to transfer the prior implicitly defined by a denoiser to a regularizing network. Hence, we aim to use our network in a simple gradient descent algorithm, to avoid the challenging parameter tuning of proximal PnP algorithms. Based on the assumption that the denoiser represents the proximal operator of an underlying differentiable regularizer (defining the image prior), we derive a loss function that links the denoiser and the regularizer's gradient. However, since there is no guarantee that this assumption is mathematically valid for a denoising neural network, we propose an approach where the pre-trained denoiser is modified jointly with the training of our regularizing network. This approach encourages the denoiser to be consistent with the definition of a proximal operator of a differentiable regularizer, and significantly improves our results in comparison to keeping the denoiser fixed.

We use our network to solve different inverse problems such as super-resolution, deblurring and pixel-wise inpainting in a simple gradient-based algorithm, and obtain better results when comparing to other generic methods. We also show that our training method can advantageously serve as a pre-training strategy, later facilitating a per-application tuning of the regularization network in the framework of unrolled gradient descent.

## 4.2 Notations and problem statement

In this section, we quickly revisit some previously established notations to ensure reader clarity and precision, and we introduce the problem that we treat in this chapter.

In this thesis, we consider linear inverse problems of the form:

$$y = Ax + \eta, \quad (4.1)$$

where $A \in \mathbb{R}^{m \times n}$ represents the degradation operator depending on the inverse problem and $\eta \sim \mathcal{N}(0, \sigma^2 I)$ typically represents Additive White Gaussian Noise (AWGN) of standard deviation $\sigma$. The reconstruction can be treated using Bayesian estimation that uses the posterior conditional probability $p(x|y)$. Maximum a posteriori probability (MAP) is
the most popular estimator in this scheme, where the estimation task is modeled as the optimization problem:

$$ \hat{x}_{MAP} = \underset{x}{\operatorname{argmin}} \frac{1}{2} \|y - Ax\|_2^2 + \sigma^2 \phi(x), \quad (4.2) $$

where $f(x) = \frac{1}{2} \|y - Ax\|_2^2$ is the data fidelity term whereas the $\phi(x)$ is the regularization term. Details of this notation can be found in Section 2.2.3. The problem in Eq. 4.2 does not have a closed-form solution in general. Therefore, it must be solved using different optimization algorithms. The Plug-and-play framework typically considers proximal splitting algorithms which decompose the problem into two sub-problems (one for each term in Eq. 4.2) and solve them alternately. In these algorithms, the regularization sub-problem consists in evaluating the proximal operator of the regularization term defined as:

$$ \begin{gathered} \operatorname{prox}_{\sigma^2\phi}(z) = \underset{x}{\operatorname{argmin}} \mathcal{F}_{\phi}(x, z, \sigma), \quad (4.3) \\ \text{with } \mathcal{F}_{\phi}(x, z, \sigma) = \frac{1}{2} \|x - z\|_2^2 + \sigma^2 \phi(x). \end{gathered} $$

This sub-problem can be replaced by a state-of-the-art Gaussian denoiser in a Plug-and-play proximal algorithm, since the proximal operator in Eq. 4.3 can be seen as a particular case of inverse problem where the degradation operator $A$ is the identity matrix, and the degradation only consists in the addition of White Gaussian Noise of standard deviation $\sigma$.

However, this approach does not directly generalize to the gradient descent algorithm where the update formula for the minimization in Eq. 4.2 is expressed as:

$$ \begin{aligned} \hat{x}_{k+1} &= \hat{x}_k - \mu \left[ \nabla f(\hat{x}_k) + \sigma^2 \nabla \phi(\hat{x}_k) \right], \\ &= \hat{x}_k - \mu \left[ A^T (A\hat{x}_k - y) + \sigma^2 \nabla \phi(\hat{x}_k) \right], \end{aligned} \quad (4.4) $$

$\mu$ being the step size. Here, instead of the proximal operator of the regularizer $\phi$, we need its gradient $\nabla\phi$, which cannot be replaced by a denoiser. In this chapter, we propose to train a network $\mathcal{G}$ that can serve as the gradient of the regularization term in the Plug-and-play gradient descent (Algorithm 4).
<caption>**Algorithm 4** Plug-and-play gradient descent</caption>

**Input:** $y, \sigma, \mu, N$
**Output:** $x$
1: $x_0 \leftarrow y$
2: **for** $k \leftarrow 0$ to $N-1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Apply $G_R \leftarrow \mathcal{G}(x_k)$
4: &nbsp;&nbsp;&nbsp;&nbsp;$$x_{k+1} \leftarrow x_k - \mu[A^T(A \cdot x_k - y) + \sigma^2 G_R]$$
5: **end for**

## 4.3 Training of the gradient of a regularizer

### 4.3.1 Mathematical derivations

We show in the following that it is mathematically possible to train a network that models the gradient of a regularizer by using a deep denoiser. Let us consider a denoiser $\mathcal{D}_\sigma$ defined as the proximal operator in Eq. 4.3:

$$ \mathcal{D}_{\sigma}(z) = \underset{x}{\operatorname{argmin}} \mathcal{F}_{\phi}(x, z, \sigma). \qquad (4.5) $$

Hence, for $\sigma$ and $z$ fixed, the denoised image $x = \mathcal{D}_\sigma(z)$ minimizes $\mathcal{F}_\phi(x, z, \sigma)$. Therefore, we have:

$$ \left. \frac{\partial \mathcal{F}_{\phi}}{\partial x} \right|_{x=\mathcal{D}_{\sigma}(z)} = 0, \qquad (4.6) $$

Furthermore, $\frac{\partial \mathcal{F}_\phi}{\partial x}$ can be computed as:

$$ \frac{\partial \mathcal{F}_{\phi}}{\partial x} = \frac{\partial \left[ \frac{1}{2} \|x-z\|_{2}^{2} + \sigma^{2} \phi(x) \right]}{\partial x}, \qquad (4.7) $$

$$ = x - z + \sigma^2 \cdot \frac{\partial \phi(x)}{\partial x}. \qquad (4.8) $$

Evaluating at the denoised image $x = \mathcal{D}_\sigma(z)$ thus gives:

$$ \left. \frac{\partial \mathcal{F}_{\phi}}{\partial x} \right|_{x=\mathcal{D}_{\sigma}(z)} = \mathcal{D}_{\sigma}(z) - z + \sigma^2 \cdot \left. \frac{\partial \phi(x)}{\partial x} \right|_{x=\mathcal{D}_{\sigma}(z)}. \qquad (4.9) $$
Using Eq. 4.6 and Eq. 4.9, we obtain:

$$ \mathcal{D}_{\sigma}(z) - z + \sigma^2 \frac{\partial \phi(x)}{\partial x} \Big|_{x=\mathcal{D}_{\sigma}(z)} = 0, \qquad (4.10) $$

$$ \sigma^2 \frac{\partial \phi(x)}{\partial x} \Big|_{x=\mathcal{D}_{\sigma}(z)} = z - \mathcal{D}_{\sigma}(z), \qquad (4.11) $$

$$ \sigma^2 \nabla \phi(\mathcal{D}_{\sigma}(z)) = z - \mathcal{D}_{\sigma}(z). \qquad (4.12) $$

Using Eq. 4.12, we can thus train a network $\mathcal{G}$ to model $\nabla\phi$ (i.e. gradient of the regularizer with respect to its input image) using the loss function:

$$ \mathcal{L}_{\mathcal{G}} = \left\| \sigma^2 [\mathcal{G}(\mathcal{D}_{\sigma}(z))] - (z - \mathcal{D}_{\sigma}(z)) \right\|_2^2. \qquad (4.13) $$

This requires the knowledge of the corresponding denoiser $\mathcal{D}_{\sigma}$. Note that Eq. 4.12 is valid for any value of $\sigma$ regardless of the degradation in $z$. Hence, $\sigma$ can be seen as a free parameter of our loss $\mathcal{L}_{\mathcal{G}}$. For small values of $\sigma$, the input $\mathcal{D}_{\sigma}(z)$ of the network $\mathcal{G}$ will be close to the degraded image $z$. Hence $\mathcal{G}$ will be trained to fit the artifacts in the degraded images (e.g. noise). On the other hand, for high values of $\sigma$, the input of $\mathcal{G}$ will be a strongly denoised image, with reduced artifacts but less details. Hence, $\mathcal{G}$ will be trained to recover the missing details. During the training, we vary the value of this parameter so that our network can recover details while also removing artifacts (see details in Section 4.3.2).

Also note that in practice, since $\phi$ is meant to be used in gradient-based algorithms, we only need the gradient $\nabla\phi$ rather than an explicit definition of $\phi$. Hence, we propose in what follows a framework for training $\mathcal{G}$ jointly with the denoiser $\mathcal{D}$, so that $\mathcal{G}$ can then be used in place of $\nabla\phi$ in Plug-and-play gradient descent. In the remainder of the paper, we refer to $\mathcal{G}$ as the regularizing gradient network (ReG).

### 4.3.2 Training framework for the regularizing gradient network

$\mathcal{G}$

The training framework is depicted in Figure 4.1. Let $\nu \to \mathcal{N}(0, \sigma_0)$ be a white Gaussian noise of standard deviation $\sigma_0$ that we use to corrupt the ground truth images $x_0$ of the training dataset to produce degraded images $z = x_0 + \nu$. Let $\sigma$ be a standard deviation value used as a parameter of our loss $\mathcal{L}_{\mathcal{G}}$, as defined in Section 4.3.1. In order to
<figure><img src="image_21.png" /><figcaption>Figure 4.1 – Framework for joint training of $\mathcal{D}$ and $\mathcal{G}$.</figcaption></figure>

handle different values of $\sigma$ in Eq. 4.13, $\mathcal{D}_\sigma$ is modeled as a non-blind deep denoiser that takes as input a noise level map (i.e. each pixel of the noise level map being equal to $\sigma$) concatenated with the noisy image $z$. This approach has been previously suggested in [2], [93], [94]. Note that if the denoiser $D_\sigma$ does not satisfy the formal definition of a MAP Gaussian denoiser for some differentiable prior, there may not exist a function $\nabla\phi$ that satisfies Eq. 4.12. We prevent this issue by starting from a pre-trained denoiser which we jointly update along with the training of the ReG network $\mathcal{G}$. In order to preserve the denoising performance of $D_\sigma$ during the training, we use an additional denoising loss $\mathcal{L}_{D_\sigma}$ defined as:

$$ \mathcal{L}_{D_\sigma} = \|\mathcal{D}_\sigma(z) - x_0\|_1. \qquad (4.14) $$

Note that Eq. 4.14 is a suitable loss for the denoiser only when $\sigma = \sigma_0$ since the non-blind denoiser must be parameterized with the true noise level $\sigma_0$ of the noisy input. The denoised output $\mathcal{D}_\sigma(z)$ is then inputted to the network modeling the gradient of the regularizer in order to train it using the loss $\mathcal{L}_G$ defined in Eq. 4.13.

Hence, our goal is to minimize the global loss defined as:

$$ \mathcal{L} = \delta \mathcal{L}_{D_\sigma} + \lambda \mathcal{L}_G, \qquad (4.15) $$

where $\lambda > 0$ and $\delta = \begin{cases} 1 & \text{if } \sigma_0 = \sigma \\ 0 & \text{otherwise.} \end{cases}$

For training the deep denoiser network, we should set the noise level $\sigma$ inputted to the network equal to the actual noise level $\sigma_0$ used for generating $\nu$. However, as explained in Section 4.3.1, for the loss $\mathcal{L}_G$, it is preferable to select $\sigma$ independently of $\sigma_0$. Hence, the input $\mathcal{D}_\sigma(z)$ of the ReG network $\mathcal{G}$ can cover a wide range of alterations, including
<figure><img src="image_22.png" /><figcaption>Figure 4.2 – Architecture of $\mathcal{G}$. The architecture is similar to the DRUNet architecture [2], with the input channel set to 3 instead of 4 (the ReG network does not need a noise level map as additional input as it does not depend on noise level).</figcaption></figure>

images with remaining noise (i.e. $\sigma < \sigma_0$) or with too strong denoising, and thus less details (i.e. $\sigma > \sigma_0$). We therefore choose to alternate during the training between either selecting independently $\sigma$ and $\sigma_0$, or setting $\sigma = \sigma_0$ in order to keep $\mathcal{D}$ faithful to the data. Furthermore, since the denoiser loss $\mathcal{L}_{\mathcal{D}_\sigma}$ is only valid when $\sigma = \sigma_0$, we omit this loss when $\sigma \neq \sigma_0$ by setting $\delta = 0$.

Note that we only need to alternate between setting $\sigma = \sigma_0$ and $\sigma \neq \sigma_0$ because we choose to train the denoiser jointly with the ReG network. A simpler training strategy would consist in first training the denoiser separately (only with the loss $\mathcal{L}_{\mathcal{D}_\sigma}$ and with $\sigma = \sigma_0$), and then training the regularizer (only with the loss $\mathcal{L}_{\mathcal{G}}$ and with $\sigma \neq \sigma_0$). However, our experimental results in Section 4.4.3 clearly show the advantage of our joint training scheme, which confirms that the separately trained denoiser does not satisfy the assumption of a Gaussian MAP denoiser for a differentiable prior.

### 4.3.3 Training details

For the training, we have used a state-of-the-art deep denoiser architecture in order to train our ReG network. We have chosen to work with the DRUNet proposed in [2] which is a combination of U-Net [95] and ResNet [96]. DRUNet additionally uses a bias-free network architecture, which has been shown to allow for good generalization of denoisers over various noise levels, even if they were not seen during training [97]. Given that it takes as input the noisy image concatenated in the channel dimension with a noise level map, it can suitably represent the non-blind denoiser $\mathcal{D}_\sigma$.
The architecture of $\mathcal{G}$ is shown in Figure 4.2. It is the same architecture as the bias-free DRUNet denoiser, with the only difference that it does not take a noise level map as additional input.

We have initialized $\mathcal{D}_\sigma$ using the pre-trained DRUNet denoiser (which we have re-produced based on the work in [2]). Then, we have trained $\mathcal{G}$ while jointly updating $\mathcal{D}_\sigma$, following the proposed framework in Section 4.3.2. The weight $\lambda$ of the loss $\mathcal{L}_\mathcal{G}$ in Eq. 4.15 has been set equal to 0.004. The parameters $\sigma$ and $\sigma_0$ have been selected following the alternating strategy described in Section 4.3.2: for half of the training iterations, we have used $\sigma = \sigma_0$ with a value chosen randomly with uniform distribution in [0, 50]; otherwise, $\sigma$ and $\sigma_0$ have been chosen independently with the same uniform distribution.

The remaining training details are similar to the ones presented in [2] for the DRUNet training: the same large dataset of 8694 images composed of images from the Waterloo Exploration Database [98], the Berkeley Segmentation Database [99], the DIV2K dataset [100] and the Flickr2K dataset [101] has been used. 16 patches of 128x128 have been randomly sampled from the training dataset for each iteration. We have used the ADAM optimizer [102] to minimize the loss $\mathcal{L}$ defined in Eq. 4.15. The learning rate has been initially set to 1e-4, and decreased by half every 100,000 iterations until reaching 5e-7, where the training stops.

## 4.4 Experimental results

### 4.4.1 Plug-and-play gradient descent

In this section, we evaluate the performance of the Plug-and-play gradient descent based on our ReG network. We refer to this approach as Plug-and-play Regularizing Gradient (PnP-ReG). We perform the evaluation on several inverse problems: super-resolution, deblurring and pixel-wise inpainting. Evaluations in this section are performed over the Set5 [103] and the CBSD68 [104] test datasets.

As the main goal of this approach is to solve inverse problems using simple gradient-based algorithms with a generic regularizer, we compare ourselves to algorithms that are designed to solve different inverse problems using a single regularization network in a Plug-and-play framework. Hence, we compare the performance of our PnP-ReG method to the PnP-ADMM with the DRUNet denoiser from [2] (see Appendix A.1 for the classical ADMM formulation and parameterization); the RED method [4] in gradient
descent with the same DRUNet denoiser; GS-PnP [3] where the network architecture used
in the denoising function is the DRUNet as well; and Chang's projection operator (One-
Net) [105] used in an ADMM framework. We also include a comparison with the Implicit
Prior of [7] only for pixel-wise inpainting, since the method is based on the assumption
that the kernel matrix is semi-orthogonal, which is not the case for the super-resolution
and deblurring applications in our work.

For fair comparisons, we have reproduced all the results under the same conditions, i.e. using the same initialization and the same degradation operator *A* for each application as described in the following subsection. We have tuned the parameters of each method on the Set5 dataset for each application to obtain the best results, and we have used the same parameters for the CBSD68 dataset.

4.4.1.1 Parameters setting and implementation details

In this section, we give further implementation and parameterization details for our
method as well as for the reference methods.

For the experimental results, we have used the ADAM optimizer [102] instead of the
simple gradient step (i.e. line 4 in Algorithm 4) to solve Eq. 4.2 in the Plug-and-play
framework. Similarly to [2], we have applied a periodical geometric self-ensemble data-
augmentation method. This involves one transformation (e.g. flipping, rotation) on the
input of the network and the counterpart inverse transformation on the output. Table 4.1
shows the parameters used during testing for PnP-ReG. The number of iterations *N* is
set to 1500. In theory, the parameter *σ* in Eq. 4.2 should be equal to the true standard
deviation *σ*<sub>*n*</sub> of the Gaussian noise added on the degraded image. However, when *σ*<sub>*n*</sub> = 0
(e.g. noiseless super-resolution, pixel-wise inpainting), choosing *σ* = 0 would completely
remove the regularization term. For these cases, we choose a small non-zero value of *σ*
depending on the application.

To reproduce the results of [105] we have used the model trained by the authors,
which takes input images of size 64x64. Hence we have applied the network on quarter-
overlapping sample patches in order to enhance the results by avoiding block artifacts.
Table A.2 in Appendix A.2 lists the tuned hyper-parameters for the reproduction of [105].

We have implemented RED [4] with gradient descent using the DRUNet denoiser.
For fair comparisons with our method, we have used the ADAM optimizer to perform
the gradient update step. The number of iterations is set to 300 for the tasks of Super-
resolution and Deblurring, and 800 in the case of pixel-wise inpainting. Table A.3 in
Appendix A.2 shows the tuned parameters for the implementation of RED.

For GS-PnP [3], we have used the parameters described in the paper, except for the standard deviation parameter of the denoiser in Super-resolution (noiseless or low-noise cases) and pixel-wise inpainting. In general, for a noise standard deviation $\sigma_n$ in the degradation, the authors of [3] suggest to parameterize the denoiser with a standard deviation of $2 * \sigma_n$. However, this is not suitable for the noiseless applications where $\sigma_n = 0$ (e.g. noiseless super-resolution, pixel-wise inpainting). Therefore, we have tuned this parameter and have set a denoiser standard deviation of 6/255 for Super-resolution in our paper. For pixel-wise inpainting, we have set this parameter to 50/255 for the first 20 iterations and to 20/255 for subsequent iterations. The remaining parameters were chosen as described in [3].

For the pixel-wise inpainting results of the Implicit Prior in [7], we have used the parameters suggested by the authors in the paper. We have computed the results by averaging 10 samples obtained with their stochastic method, as done in the original paper.

Finally, we have compared with the PnP-ADMM using DRUNet [2]. The PnP-ADMM formulation and its parameterization is elaborated in Appendix A.1. Furthermore, we have applied the periodical geometric self-ensemble data augmentation during the test as done in [2]. Table A.1 in Appendix A.1 shows the parameters that yielded the best results for the PnP-ADMM. We note that in the cases where we have added noise, fixing the standard deviation parameter of the denoiser gave better results than decreasing it throughout the iterations.

### 4.4.1.2 Super-resolution

First, we test our approach on the task of super-resolution. Low resolution images have been generated by applying a convolution kernel followed by a downsampling by a factor $t$. We have evaluated our method with bicubic and Gaussian convolutional kernels, with both 2x and 3x downsampling scales and Gaussian noise with 3 different noise levels $\sigma_n = \{0.0, \sqrt{2}, 2.55\}/255$. The Gaussian kernel has a standard deviation of $\sigma_b = 0.5 \cdot t$ (i.e. $\sigma_b = 1$ for x2 and $\sigma_b = 1.5$ for x3). In all the cases, the gradient descent has been initialized with a high resolution image obtained by the bicubic upsampling of the degraded image.

Tables 4.2 and 4.3 give a numerical comparison of our method with the aforementioned generic approaches for super-resolution of factor 2 and 3 respectively, for both bicubic and Gaussian filters. Numerical comparison gives higher values for PnP-ReG compared to the
<table><caption>Table 4.1 – Parameters used for our PnP-ReG method. $\mu$: gradient step size, $\sigma$: weight of the regularization, $\sigma_n$: standard deviation of the AWGN added on the degraded image.</caption><thead><tr><th></th><th></th><th>σ<sub>n</sub> * 255</th><th>μ</th><th>σ * 255</th></tr></thead><tbody><tr><td rowspan="3">Super-resolution x2</td><td rowspan="3">Bicubic</td><td>0</td><td>0.008</td><td>1.2</td></tr><tr><td>√2</td><td>0.002</td><td>√2</td></tr><tr><td>2.55</td><td>0.002</td><td>2.55</td></tr><tr><td rowspan="3">Gaussian</td><td rowspan="3">Bicubic</td><td>0</td><td>0.008</td><td>0.8</td></tr><tr><td>√2</td><td>0.002</td><td>√2</td></tr><tr><td>2.55</td><td>0.002</td><td>2.55</td></tr><tr><td rowspan="3">Super-resolution x3</td><td rowspan="3">Bicubic</td><td>0</td><td>0.002</td><td>0.9</td></tr><tr><td>√2</td><td>0.002</td><td>√2</td></tr><tr><td>2.55</td><td>0.002</td><td>2.55</td></tr><tr><td rowspan="6">Deblurring</td><td rowspan="3">Gaussian</td><td>0</td><td>0.004</td><td>0.4</td></tr><tr><td>√2</td><td>0.002</td><td>√2</td></tr><tr><td>2.55</td><td>0.002</td><td>2.55</td></tr><tr><td rowspan="3">Pixel-wise inpainting</td><td>0.1</td><td>0</td><td>0.025</td><td>3.6/255</td></tr><tr><td>0.2</td><td>0</td><td>0.01</td><td>1/255</td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

existing generic Plug-and-play approaches. Figure 4.3 shows a visual comparison of the results for a noiseless degradation with a bicubic kernel and a downsampling by a factor of 2. We observe sharper images with less aliasing artifacts produced by our approach.

### 4.4.1.3 Deblurring

We have degraded our images using two 25x25 isotropic Gaussian blur kernels of standard deviations of 1.6 and 2.0, as well as two anisotropic kernels that have been used in [106]. We have considered White Gaussian noise with 3 noise levels $\sigma_n = \{\sqrt{2}, 2.55, 7.65\}/255$. The blurred image is directly used as the initialization of the Plug-and-play gradient descent.

Table 4.4 gives the PSNR results [dB] obtained with our method when deblurring
<figure><img src="image_23.png" /><figcaption>Figure 4.3 – Visual comparison of super-resolution results obtained with the projection operator (One-Net) [105], GS-PnP [3], RED [4], PnP-ADMM with the denoiser of [2] and our PnP-ReG method. Low resolution images generated with a bicubic kernel followed by a downsampling by a factor of 2.</figcaption></figure>

images which have been degraded with isotropic and anisotropic kernels. We observe
higher PSNR values with respect to the other generic methods for most cases when the
added noise levels are σ_n = {√2, 2.55}/255. However, by increasing the noise level to
7.65/255, the best results have been obtained with GS-PnP. The visual comparison in
Figure 4.4 shows that our approach successfully recovers the details without increasing
the noise.
<table><caption>Table 4.2 – Super-resolution results (measured in PSNR [dB]) obtained with our PnP-ReG method; the input images have been corrupted using bicubic and Gaussian kernels followed by downsampling by a factor of 2 and adding Gaussian noise with 3 different noise levels $\sigma_n = \{0.0, \sqrt{2}, 2.55\}/255$. Comparison with the projection operator (One-Net) [105], GS-PnP [3], RED [4] and the PnP-ADMM with the denoiser of [2].</caption><thead><tr><th rowspan="2"></th><th rowspan="2">σ<sub>n</sub> * 255</th><th colspan="3">(i) Bicubic</th><th colspan="3">(ii) Gaussian</th></tr><tr><th>0.0</th><th>√2</th><th>2.55</th><th>0.0</th><th>√2</th><th>2.55</th></tr></thead><tbody><tr><td rowspan="5">Set5</td><td>OneNet</td><td>33.22</td><td>33.70</td><td>33.20</td><td>32.67</td><td>33.08</td><td>32.45</td></tr><tr><td>GS-PnP</td><td>34.58</td><td>34.49</td><td><b>34.31</b></td><td>33.98</td><td>33.88</td><td><b>33.59</b></td></tr><tr><td>RED</td><td>35.05</td><td>34.49</td><td>33.78</td><td>34.99</td><td>33.80</td><td>32.84</td></tr><tr><td>PnP-ADMM</td><td>35.20</td><td>34.42</td><td>33.80</td><td>35.14</td><td>33.69</td><td>32.74</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>35.34</b></td><td><b>34.90</b></td><td>34.29</td><td><b>35.30</b></td><td><b>34.33</b></td><td>33.41</td></tr><tr><td rowspan="5">CBSD68</td><td>OneNet</td><td>29.48</td><td>29.80</td><td>29.56</td><td>29.08</td><td>29.48</td><td>29.06</td></tr><tr><td>GS-PnP</td><td>29.95</td><td>29.93</td><td>29.81</td><td>29.57</td><td>29.51</td><td>29.39</td></tr><tr><td>RED</td><td>30.45</td><td>30.14</td><td>29.73</td><td>30.39</td><td>29.61</td><td>28.99</td></tr><tr><td>PnP-ADMM</td><td>30.48</td><td>30.17</td><td>29.92</td><td>30.42</td><td>29.78</td><td>29.23</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>30.55</b></td><td><b>30.34</b></td><td><b>30.06</b></td><td><b>30.50</b></td><td><b>30.01</b></td><td><b>29.44</b></td></tr></tbody></table>

### 4.4.1.4 Pixel-wise inpainting

The degradation in the pixel-wise inpainting task consists of multiplying the ground-truth image by a binary mask. We have tested our approach with both 20% and 10% of known pixels rates. For the initialization image $\hat{x}^0$, we have set the color of the unknown pixels to grey.

Table 4.5 shows a numerical evaluation of our method for the application of pixel-wise inpainting in terms of PSNR [dB]. We observe significant performance gains of the PnP-ReG method compared to the other methods, especially in the most challenging case where the known pixel rate is only 10%. Visual comparisons are given in Figure 4.5.

## 4.4.2 Unrolled gradient descent with $G$

Aside from the Plug-and-play gradient descent, our approach for training $G$ can also serve as a pre-training strategy for unrolled gradient descent. In unrolled optimization methods, the regularization network is trained for each inverse problem such that applying a fixed number of iterations of the algorithm (e.g. Eq. 4.4 for gradient descent) best approximates the ground truth image. We describe the unrolled training approach in
<table><caption>Table 4.3 – Super-resolution results (measured in PSNR [dB]) obtained with our PnP-ReG method; the input images have been corrupted using bicubic and Gaussian kernels followed by downsampling by a factor of 3 and adding Gaussian noise with 3 different noise levels $\sigma_n = \{0.0, \sqrt{2}, 2.55\}/255$. Comparison with the projection operator (One-Net) [105], GS-PnP [3], RED [4] and the PnP-ADMM with the denoiser of [2].</caption><thead><tr><th rowspan="2"></th><th rowspan="2">σ<sub>n</sub> * 255</th><th colspan="3">(i) Bicubic</th><th colspan="3">(ii) Gaussian</th></tr><tr><th>0.0</th><th>√2</th><th>2.55</th><th>0.0</th><th>√2</th><th>2.55</th></tr></thead><tbody><tr><td rowspan="5">Set5</td><td>OneNet</td><td>29.65</td><td>30.18</td><td>29.79</td><td>29.52</td><td>29.71</td><td>29.05</td></tr><tr><td>GS-PnP</td><td>31.48</td><td>31.43</td><td><b>31.24</b></td><td>30.91</td><td>30.84</td><td><b>30.59</b></td></tr><tr><td>RED</td><td>31.47</td><td>31.22</td><td>30.78</td><td>31.44</td><td>30.77</td><td>30.05</td></tr><tr><td>PnP-ADMM</td><td>31.49</td><td>31.05</td><td>30.39</td><td>31.45</td><td>30.22</td><td>29.17</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>31.75</b></td><td><b>31.46</b></td><td>31.13</td><td><b>31.60</b></td><td><b>31.09</b></td><td>30.39</td></tr><tr><td rowspan="5">CBSD68</td><td>OneNet</td><td>26.17</td><td>26.89</td><td>26.67</td><td>25.21</td><td>26.69</td><td>26.33</td></tr><tr><td>GS-PnP</td><td>27.11</td><td>27.08</td><td>26.99</td><td>26.80</td><td>26.73</td><td>26.64</td></tr><tr><td>RED</td><td>27.50</td><td>27.33</td><td>27.11</td><td>27.40</td><td>27.06</td><td>26.24</td></tr><tr><td>PnP-ADMM</td><td>27.51</td><td>27.33</td><td>26.97</td><td><b>27.47</b></td><td>26.83</td><td>26.20</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>27.55</b></td><td><b>27.41</b></td><td><b>27.23</b></td><td>27.46</td><td><b>27.17</b></td><td><b>26.77</b></td></tr></tbody></table>

Algorithm 5 for a gradient descent optimization. While this end-to-end training strategy loses the genericity of the Plug-and-play approach, it typically improves the performance.

However, to facilitate the training, it is generally required to initialize the network weights with a generically pre-trained version. When unrolling proximal algorithms such as ADMM, a pre-trained deep denoiser can be used since it can be interpreted as the proximal operator of a generic regularization function. On the other hand, the gradient descent requires instead the gradient of a regularizer. Hence, a denoiser cannot be directly used as a pre-trained network. By transferring the image prior implicitly represented by the denoiser $\mathcal{D}$ to our ReG network $\mathcal{G}$, our method thus provides a satisfying pre-trained network for unrolled gradient descent.

We have tested this approach for noiseless super-resolution of factor 2 and 3, as well as for deblurring using 2 isotropic Gaussian kernels (with standard deviations of 1.6 and 2.0) and additive Gaussian noise of standard deviation $\sqrt{2}/255$. We have compared our results with those obtained when using a network learned end-to-end in an unrolled environment without the pre-training (i.e. random weight initialization).

Both versions (pre-trained and not pre-trained) have been unrolled in the same conditions. Table 4.6 gives the training parameters for each of the different tasks. For all 4
<table><caption>Table 4.4 – Deblurring results (measured in PSNR [dB]) obtained with our PnP-ReG method. The input blurred images have been generated using two isotropic Gaussian kernels of respective standard deviations 1.6 and 2.0, and two anisotropic Gaussian kernels from [2] followed by adding Gaussian noise with 3 different noise levels $\sigma_n = \{\sqrt{2}, 2.55, 7.65\}/255$. We compare with the projection operator (One-Net) [105], GS-PnP [3], RED [4] and the PnP-ADMM with the denoiser of [2].</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th rowspan="2">σ<sub>n</sub> * 255</th><th colspan="3">Image</th><th colspan="3">Image</th></tr><tr><th>√2</th><th>2.55</th><th>7.65</th><th>√2</th><th>2.55</th><th>7.65</th></tr></thead><tbody><tr><td rowspan="5">Set5</td><td>OneNet</td><td>32.18</td><td>31.54</td><td>27.86</td><td>30.41</td><td>29.52</td><td>27.20</td></tr><tr><td>GS-PnP</td><td>32.84</td><td>32.31</td><td><b>30.93</b></td><td>31.04</td><td>29.98</td><td><b>29.76</b></td></tr><tr><td>RED</td><td>32.63</td><td>31.76</td><td>29.83</td><td>31.25</td><td>30.62</td><td>29.21</td></tr><tr><td>PnP-ADMM</td><td>32.83</td><td>32.06</td><td>30.43</td><td>31.55</td><td>30.88</td><td>29.30</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>33.40</b></td><td><b>32.51</b></td><td>30.63</td><td><b>31.96</b></td><td><b>31.19</b></td><td>29.46</td></tr><tr><td rowspan="5">CBSD68</td><td>OneNet</td><td>28.47</td><td>28.27</td><td>25.86</td><td>26.91</td><td>26.97</td><td>25.29</td></tr><tr><td>GS-PnP</td><td>28.98</td><td><b>28.64</b></td><td><b>27.64</b></td><td>27.5</td><td>27.33</td><td><b>26.57</b></td></tr><tr><td>RED</td><td>29.02</td><td>28.16</td><td>26.91</td><td>27.63</td><td>27.24</td><td>26.20</td></tr><tr><td>PnP-ADMM</td><td>29.21</td><td>28.54</td><td>27.20</td><td>27.91</td><td>27.40</td><td>26.25</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>29.38</b></td><td>28.62</td><td>27.07</td><td><b>27.99</b></td><td><b>27.41</b></td><td>26.21</td></tr><tr><td rowspan="5"></td><td rowspan="5">σ<sub>n</sub> * 255</td><td colspan="3">Image</td><td colspan="3">Image</td></tr><tr><th>√2</th><th>2.55</th><th>7.65</th><th>√2</th><th>2.55</th><th>7.65</th></tr><tr><td>OneNet</td><td>29.26</td><td>28.90</td><td>26.71</td><td>28.93</td><td>28.18</td><td>26.76</td></tr><tr><td>GS-PnP</td><td>30.18</td><td>29.30</td><td><b>29.13</b></td><td>29.70</td><td>29.07</td><td><b>28.55</b></td></tr><tr><td>RED</td><td>30.53</td><td>29.98</td><td>28.67</td><td>29.91</td><td>29.53</td><td>28.24</td></tr><tr><td rowspan="5">Set5</td><td>OneNet</td><td>29.26</td><td>28.90</td><td>26.71</td><td>28.93</td><td>28.18</td><td>26.76</td></tr><tr><td>GS-PnP</td><td>30.18</td><td>29.30</td><td><b>29.13</b></td><td>29.70</td><td>29.07</td><td><b>28.55</b></td></tr><tr><td>RED</td><td>30.53</td><td>29.98</td><td>28.67</td><td>29.91</td><td>29.53</td><td>28.24</td></tr><tr><td>PnP-ADMM</td><td>31.21</td><td>30.56</td><td>28.95</td><td>30.33</td><td>29.07</td><td>28.17</td></tr><tr><td>PnP-ReG (Ours)</td><td><b>31.26</b></td><td><b>30.57</b></td><td>28.84</td><td><b>30.67</b></td><td><b>29.95</b></td><td>28.31</td></tr><tr><td rowspan="5">CBSD68</td><td>OneNet</td><td>26.52</td><td>26.57</td><td>25.04</td><td>25.88</td><td>25.94</td><td>24.77</td></tr><tr><td>GS-PnP</td><td>27.02</td><td>27.04</td><td><b>26.32</b></td><td>26.27</td><td>26.41</td><td><b>25.72</b></td></tr><tr><td>RED</td><td>27.36</td><td>26.99</td><td>26.00</td><td>26.61</td><td>26.36</td><td>25.53</td></tr><tr><td>PnP-ADMM</td><td><b>27.83</b></td><td><b>27.32</b></td><td>26.16</td><td>26.9</td><td>26.46</td><td>25.46</td></tr><tr><td>PnP-ReG (Ours)</td><td>27.71</td><td>27.16</td><td>25.96</td><td><b>26.93</b></td><td><b>26.47</b></td><td>25.47</td></tr></tbody></table>

cases, we have trained the networks using the DIV2K dataset [100] of 800 images, over 600 epochs by randomly taking 48x48 patches from the dataset. In addition, we include a comparison with Total Deep Variation (TDV) method [88]. TDV also performs unrolled
<figure><img src="image_24.png" /><figcaption>Figure 4.4 – Visual comparison of deblurring results obtained with the projection operator (One-Net) [105], GS-PnP [3], RED [4], PnP-ADMM with the denoiser of [2] and our PnP-ReG method. The blurred images have been generated using an isotropic Gaussian kernel of standard deviation σ_b = 1.6 followed by adding Gaussian noise of standard deviation σ_n = √2/255.</figcaption></figure>

optimization, where the network represents the gradient of the regularization function. In [88], the network is not pre-trained. Instead, the training starts with a small number of unrolled iterations *N* = 2, and *N* is incremented every 700 epochs. Note that the authors originally trained the TDV network for *N* = 10 iterations of an unrolled proximal gradient descent algorithm. However, for fair comparisons with our approach, we re-trained it with *N* = 6 iterations.

Tables 4.7 and 4.8 show the PSNR results [dB] for both networks (pre-trained and not pre-trained) as well as for the TDV, for super-resolution and deblurring respectively, using Set5 [103], Set14 [27] and BSDS100[104]. As expected, using $\mathcal{G}$ for weight initialization improves the results of the unrolled gradient descent for the 4 tested cases (with up to 0.9 dB).

Some visual comparisons of super-resolution results with magnifying factor of 3, and of deblurring images corrupted by a Gaussian kernel of standard deviation 2.0 are respec-
<figure><img src="image_25.png" /><figcaption>Figure 4.5 – Visual comparison of pixel-wise inpainting results with known pixel rate of *p* = 20%, obtained with the projection operator (One-Net) [105], the Implicit Prior (IP) of [7], GS-PnP [3], RED [4], PnP-ADMM with the denoiser of [2] and our PnP-ReG method.</figcaption></figure>
<table><caption>Table 4.5 – Pixel-wise inpainting results (measured in PSNR [dB]) obtained with our PnP-ReG method; the corrupted images have been generated by keeping 20% and 10% of the known pixels. Comparison with the projection operator (One-Net) [105], the Implicit Prior (IP) of [7], GS-PnP [3], RED [4] and PnP-ADMM with the denoiser of [2].</caption><thead><tr><th rowspan="2"></th><th colspan="2">Set5</th><th colspan="2">CBSD68</th></tr><tr><th>(i) 10%</th><th>(ii) 20%</th><th>(i) 10%</th><th>(ii) 20%</th></tr></thead><tbody><tr><td>OneNet</td><td>17.20</td><td>24.06</td><td>14.72</td><td>20.46</td></tr><tr><td>IP (avg)</td><td>25.28</td><td>28.91</td><td>24.12</td><td>26.55</td></tr><tr><td>GS-PnP</td><td>25.08</td><td>28.70</td><td>23.58</td><td>25.91</td></tr><tr><td>RED</td><td>22.75</td><td>27.17</td><td>22.24</td><td>23.22</td></tr><tr><td>PnP-ADMM</td><td>26.20</td><td>30.20</td><td>24.06</td><td>26.75</td></tr><tr><td>PnP-ReG (Ours)</td><td><strong>26.94</strong></td><td><strong>30.36</strong></td><td><strong>24.43</strong></td><td><strong>27.09</strong></td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>Figure 4.6 – Visual comparison of super-resolution results between unrolled gradient descent with and without pre-training. Low resolution images have been generated with a bicubic kernel followed by a downsampling by a factor of 3.</figcaption></figure>

tively shown in Figures 4.6 and 4.7. The visual results confirm that better reconstruction of the details is obtained when our pre-training is used. While the TDV results display even sharper edges, the PSNR remains lower because of exaggerated sharpness in comparison to the ground truth.

### 4.4.3 Analysis of the joint training

In this section, we experimentally verify the advantages of updating the denoising network within the training of our ReG network $\mathcal{G}$, compared to letting the denoiser fixed. In the latter case, $\delta$ is always set to 0 (i.e. $\mathcal{L} = \lambda \mathcal{L}_{\mathcal{G}}$).

First, we compare the performance of our regularizing gradient network trained in
<caption>**Algorithm 5** Unrolled gradient descent with $\mathcal{G}$</caption>

1: initialize $\mathcal{G}$, $\sigma_n$, $\mu$, $A$, $B$
2: **for** each epoch **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;**for** each batch **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\eta \leftarrow \mathcal{N}(0, \sigma_n)$
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_{gt} \leftarrow$ ground-truth batch
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$y \leftarrow Ax_{gt} + \eta$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_0 \leftarrow y$
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $k \leftarrow 0$ to $N-1$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_{k+1} \leftarrow x_k - \mu[A^T(Ax_k - y) + \sigma^2\mathcal{G}(x_k)]$
10: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
11: &nbsp;&nbsp;&nbsp;&nbsp;loss $\leftarrow \|x_N - x_{gt}\|_2^2$
12: &nbsp;&nbsp;&nbsp;&nbsp;Update weights of $\mathcal{G}$ with back-propagation
13: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
14: **end for**

<figure><img src="image_27.png" /><figcaption>Figure 4.7 – Visual comparison of deblurring results between unrolled gradient descent with and without pre-training. The blurred images have been generated using an isotropic Gaussian kernel of standard deviation 2.0.</figcaption></figure>
Table 4.6 – Parameters used for unrolled gradient descent optimization for both the pre-trained and not pre-trained versions (i) Super-resolution of factor 2 and 3 with input images corrupted using a bicubic filter (ii) Deblurring images corrupted using an isotropic Gaussian kernel of standard deviation $\sigma_b$ of 1.6 and 2.0. $\sigma_n$: Standard deviation of the White Gaussian noise added on the corrupted image, $\sigma$: weight of the regularization term, $\mu$: gradient step size, $N$: number of unrolled iterations. The training is performed over the DIV2K dataset [100].

<table><thead><tr><th rowspan="2"></th><th colspan="2">(i) Super-resolution</th><th colspan="2">(ii) Deblurring</th></tr><tr><th>x2</th><th>x3</th><th>σ<sub>b</sub> = 1.6</th><th>σ<sub>b</sub> = 2.0</th></tr></thead><tbody><tr><td>σ<sub>n</sub> * 255</td><td>0</td><td>0</td><td>√2</td><td>√2</td></tr><tr><td>σ * 255</td><td>1.2</td><td>1.6</td><td>√2</td><td>√2</td></tr><tr><td>μ</td><td>0.008</td><td>0.1</td><td>0.004</td><td>0.004</td></tr><tr><td>N</td><td>6</td><td>6</td><td>6</td><td>6</td></tr></tbody></table>

both scenarios. An example of deblurring results with the Plug-and-play gradient descent is shown in Figure 4.8. It is clear that leaving the denoiser fixed to its pre-trained state degrades the performance of $G$: the reconstructed image in Figure 4.8 (a) remains more blurry than that in Figure 4.8 (b) and it also presents colored fringes artifacts. On the other hand, when the denoiser is updated when training $G$, the convergence of the Plug-and-play gradient descent is significantly improved as well as the visual result, as shown in Figure 4.8 (b,c).

A possible explanation of the worse results when fixing the denoiser in our training is that the pre-training of $D_\sigma$ does not guarantee that there exists a differentiable regularizer $\phi$ for which $\text{prox}_{\sigma^2\phi} = D_\sigma$ for every value of $\sigma$. In other words, the assumption that $D_\sigma$ is a MAP Gaussian denoiser for a differentiable prior may not be satisfied. However, by jointly updating the denoiser with the ReG network, the modified denoiser better represents such a MAP Gaussian denoiser for the corresponding regularizer.

In a second experiment, we compare the performances of the updated denoiser and the original DRUNet when used in the Plug-and-play ADMM algorithm. Figure 4.9 shows for each ADMM iteration the PSNR of the reconstructed images and the MSE of the difference between two consecutive iterations for both versions of the denoiser. We can observe that although the original DRUNet can obtain better PSNR performances when stopping the ADMM after a few iterations (see subfigure (a)), the algorithm does not converge, and may even strongly diverge after a sufficiently large number of iterations.
Table 4.7 – Super-resolution results (measured in PSNR [dB]) obtained with unrolled gradient descent (the input images have been corrupted using a bicubic kernel and a downsampling factor of 2 and 3). The evaluation has been performed using the Set5, Set14 and BSDS100 datasets. The restoration has been performed using unrolled gradient descent initialized with our pre-trained network $\mathcal{G}$ and without weight initialization, as well as TDV[88].

<table><thead><tr><th></th><th></th><th>pre-trained<br>(ours)</th><th>not<br>pre-trained</th><th>TDV</th></tr></thead><tbody><tr><td rowspan="3">(i) x2</td><td>Set5</td><td>35.61</td><td>35.42</td><td>34.57</td></tr><tr><td>Set14</td><td>31.18</td><td>30.93</td><td>30.31</td></tr><tr><td>BSDS100</td><td>30.77</td><td>30.68</td><td>30.23</td></tr><tr><td rowspan="3">(i) x3</td><td>Set5</td><td>32.10</td><td>31.88</td><td>31.47</td></tr><tr><td>Set14</td><td>28.00</td><td>27.79</td><td>27.40</td></tr><tr><td>BSDS100</td><td>27.68</td><td>27.63</td><td>27.34</td></tr></tbody></table>

On the other hand our modified denoiser allows for a better convergence of the Plug-and-play ADMM.

## 4.5 Discussion

In this section, we discuss the advantages and the limitations of our PnP-ReG method. First, we point out that the proposed approach only requires the tuning of the gradient step $\mu$ when the noise level is known, while the reference PnP methods heavily rely on the tuning of several hyper-parameters, sometimes even including the number of iterations. This makes our method easier to use since parameter tuning can take a lot of effort.

Furthermore, the comparison of the convergence of the different PnP algorithms in Figures 4.10 and 4.11 show that although our method requires more iterations than PnP-ADMM and GS-PnP to reach its highest PSNR result, it converges to a fixed point that provides the best quality. On the other hand, PnP-ADMM reaches its highest PSNR after a few iterations, but then diverges, which requires a careful tuning of the number of iterations for each application in practice. The convergence of our PnP-ReG method is more comparable to the RED-GD method (although significant PSNR gains are obtained with our method). This may be explained by the fact that both RED-GD and PnP-ReG are based on gradient descent while PnP-ADMM and GS-PnP use proximal algorithms.

Aside from the convergence analysis, one may notice that, even though we have de-
<table><caption>Table 4.8 – Deblurring results (measured in PSNR [dB]) obtained with unrolled gradient descent (the input images have been corrupted using isotropic Gaussian kernels with standard deviation 1.6 and 2.0 and Gaussian noise of standard deviation $\sqrt{2}/255$). The evaluation has been performed using the Set5, Set14 and BSDS100 datasets. The restoration has been performed using unrolled gradient descent initialized with our pre-trained network $\mathcal{G}$ and without weight initialization, as well as [88].</caption><thead><tr><th></th><th></th><th>pre-trained<br>(ours)</th><th>not<br>pre-trained</th><th>TDV</th></tr></thead><tbody><tr><td rowspan="3">(i) $\sigma_0 = 1.6$</td><td>Set5</td><td>33.51</td><td>32.57</td><td>32.32</td></tr><tr><td>Set14</td><td>30.14</td><td>29.33</td><td>29.17</td></tr><tr><td>BSDS100</td><td>29.80</td><td>29.12</td><td>29.01</td></tr><tr><td rowspan="3">(ii) $\sigma_0 = 2$</td><td>Set5</td><td>31.63</td><td>30.94</td><td>30.52</td></tr><tr><td>Set14</td><td>28.29</td><td>27.70</td><td>27.46</td></tr><tr><td>BSDS100</td><td>28.08</td><td>27.52</td><td>27.51</td></tr></tbody></table>

signed our loss from Eq. 4.12 so that $\mathcal{G}$ models $\nabla\phi$, we did not enforce the network to have a symmetric Jacobian matrix. As a result, our training framework does not guarantee that $\mathcal{G}$ can be interpreted mathematically as a conservative vector field, and thus, as the gradient of a scalar potential. However, using a non-conservative vector field for the gradient update step may not necessarily degrade the gradient descent performance. For example, advanced gradient-based algorithms typically alter the gradient (e.g. with momentum), in a way that loses the conservativeness property of the original gradient, while improving the algorithm’s performances (more robust to local minima, faster convergence...). Note also that it would be possible to enforce the Jacobian symmetry property in our method by using a network directly modeling $\phi$ and by explicitly computing its gradient, as done in [3], [107]. However, the explicit gradient computation has the same complexity as the network $\phi$, hence doubling the computation cost. A possible direction for future work would be to study whether such a constraint could improve the Plug-and-play gradient descent without sacrificing the computational complexity.

## 4.6 Conclusion

In this chapter, we have proposed a novel framework for regularizing linear inverse problems. Our approach makes it possible to solve Plug-and-play algorithms using gradient descent, where the gradient of the regularizer is required rather than its proximal
<figure><img src="image_28.png" /><figcaption>Figure 4.8 – Comparison of the performances of PnP-ReG, when the ReG network is trained (a) with a fixed denoiser and (b) while jointly updating the denoiser. (c) PSNR [dB] values over the gradient descent iterations. The results are shown for the deblurring problem (the blurred images have been generated using isotropic Gaussian kernels of standard deviation 1.6 followed by adding Gaussian noise of standard deviation $\sqrt{2}/255$).</figcaption></figure>

operator. We have proved that it is mathematically possible to train a network that models the gradient of a regularizer, jointly with a denoising neural network.

The results have demonstrated that the joint training of our network with a DRUNet denoiser has several advantages. First, our network can be used to regularize the gradients in a Plug-and-play gradient descent algorithm and can outperform other generic approaches in different inverse problems such as super-resolution, deblurring and pixel-
<figure><img src="image_29.png" /><figcaption>Figure 4.9 – Comparison of the performances of the original and the updated denoisers in PnP-ADMM for deblurring, using the Set5 dataset (the blurred images have been generated using an isotropic Gaussian kernel of standard deviation 1.6 followed by adding Gaussian noise of standard deviation $\sqrt{2}/255$): (a) Average PSNR [dB] over the ADMM iterations (b) MSE of the difference between two consecutive iterations (in log scale).</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 4.10 – Comparison of the convergence of PnP-ADMM with the denoiser of [2], RED [4], GS-PnP [3] and PnP-ReG for Super-resolution. The input low resolution images have been generated using bicubic downsampling with a factor 2. The results are averaged over the images of the Set5 dataset. For PnP-ADMM, we have used the setting with variable parameter $s^k$ until the 25th iteration for which the best results are obtained (see Table A.1), and let $s^k$ fixed afterwards.</figcaption></figure>
<figure><img src="image_31.png" /></figure>

Figure 4.11 – Comparison of the convergence of PnP-ADMM with the denoiser of [2], RED [4], GS-PnP [3] and PnP-ReG for Deblurring. The degraded images have been generated using an isotropic Gaussian kernel of standard deviation 2.0 followed by adding Gaussian noise of standard deviation $\sigma_n = 0.01$. The results are averaged over the images of the Set5 dataset. For PnP-ADMM, we have used the setting with fixed parameter $s^k = 30/255$ for which the best results are obtained (see Table A.1).

wise inpainting. Second, our network can also serve as a pre-training strategy for unrolled gradient descent and yield a significant improvement. Lastly, the joint training of the denoiser with the regularizing gradient network makes the former match better the definition of a proximal operator compared to the original pre-trained DRUNet.
<figure><img src="image_32.png" /></figure>
# Inverse problems in omnidirectional imagery
<figure><img src="image_33.png" /></figure>
# A COMPREHENSIVE REVIEW OF OMNIDIRECTIONAL IMAGES: PROCESSING METHODOLOGIES AND INVERSE PROBLEMS

The second part of this thesis is dedicated to advancing solutions for inverse problems in omnidirectional images. In this chapter, we provide a comprehensive review on omnidirectional images, exploring their acquisition, representation, challenges and the existing literature on inverse problems.

## 5.1 Introduction to omnidirectional imaging

Omnidirectional images, also known as spherical images or 360° images, are high-resolution visual contents that cover a 360° field of view capturing the light field converging to a single point. Thanks to their human-friendly nature, they have gained a particular attention in the domains of Virtual Reality (VR), Augmented Reality (AR) and robotics.

Nevertheless, the acquisition and the display of omnidirectional images present complex challenges. Also, traditional 2D processing tools can not be directly applied to spherical data. In fact, when dealing with omnidirectional images, the geometry is fundamentally different. A sphere has a curved surface, and points on its surface are not represented by simple *x* and *y* coordinates as in a 2D plane. Instead, spherical coordinates (i.e. polar angle $\Theta \in [0, \pi]$ and azimuth angle $\phi \in [-\pi, \pi]$) are used to indicate locations on the surface of the sphere. Hence, processing omnidirectional images is not a straightforward task.

Processing 2D images is relatively straightforward: we visualize and process the captured image after applying minor pre-processing steps (Figure 5.1). However, this sim-
<figure><img src="image_34.png" /><figcaption>Figure 5.1 – Image processing pipeline for 2D images. Figure extracted from [108].</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 5.2 – Image processing pipeline for omnidirectional imagery. Figure extracted from [108].</figcaption></figure>

plicity doesn't extend to omnidirectional imaging, where we don't directly process and visualize the captured data. Figure 5.2 represents a typical processing pipeline for omnidirectional imagery. First, it is not feasible, using a single camera, to capture the light ray coming from all directions effectively covering the entire sphere. Multiple captures should be merged in order to build an omnidirectional image that covers the entire sphere surface. Then, a mapping is typically used to process the spherical image. Finally, given the finite field of view in human vision, only a portion of the spherical image becomes visible. This necessitates projecting the relevant section onto a *viewport*. In other words, omnidirectional image processing is not a straightforward task and demands specialized processing tools and techniques. In the next section, we briefly discuss about the acquisition process of 360° images, then elaborate different representation models in Section 5.3.
## 5.2 Acquisition of omnidirectional images

Capturing spherical data can be done in different ways. It is important to note, that achieving a full 360° coverage using a single camera is not possible. Instead, several captures are typically combined using stitching methods, in order to build one spherical image. One approach involves the utilization of multiple perspective cameras, each covering a small angle of view. This approach is expensive, since it requires a lot of processing and a significant amount of perspective cameras to obtain one spherical image.

Moreover, some cameras can cover a wider angle of view, hence a fewer number of cameras are needed to cover the whole sphere. Two popular such devices are catadoptric and fish-eye based cameras. The former employs a curved mirror in its optical setup, which reflects incoming light rays from the scene onto the camera sensor. However, due to the substantial space required for a second symmetric mirror, it is constrained to capturing only one hemisphere. On the other hand, fish-eye cameras have gained significant popularity in recent years. These cameras use a special fish-eye lens that distorts the incoming light to cover a wide field of view. Fish-eye cameras equip nowadays most budget-friendly cameras available on the public market.

For more details on omnidirectional image acquisition, we refer interested readers to the Chapter 4 of Tzafestas’ book [109], as well as to [108], [110].

## 5.3 Projecting spherical data to the Euclidean space

Omnidirectional images rely on a spherical geometry, which is not an Euclidean space where 2D processing tools are well defined. In order to take advantage of computing tools and methods that have been developed for 2D images, spherical images are often mapped to two-dimensional planar representations such as equirectangular projection (ERP) [9], cubemap projection (CMP) [111] or rhombic dodecahedron projection [112]. While a vast number of such mappings have been designed<fnref n="1" />, in this section, we discuss two popular mapping methods that have been commonly adopted in the literature.

### 5.3.1 Equirectangular projection (ERP)

Equirectangular projection, also known as panorama projection, is the most popular mapping used to represent spherical data. It projects the omnidirectional image to

<footnote n="1">1. https://map-projections.net</footnote>
the Euclidean space by mapping each point on the sphere's surface to a corresponding
point on the grid, using a linear relationship between longitude (azimuthal angle) and
$x$-coordinate, and a linear relationship between latitude (polar angle) and $y$-coordinate.
However, this sampling induces significant distortions as we move towards the poles. It
can be seen in Figure 5.3, the red and blue portions are mapped to two equal regions on
the projection, whereas in reality, their sizes on the sphere are very different. Moreover,
this mapping causes strong deformation of objects (Figure 5.5(b)).

<figure><img src="image_36.png" /><figcaption>Figure 5.3 – Equirectangular projection. Figure extracted from [108].</figcaption></figure>

**5.3.2 Cubemap projection (CMP)**

Cubemap projection is a mapping technique that projects the spherical image onto the six faces of a cube (Figure 5.4). The projection on each of the faces imitates the capture of a perspective camera with a $\frac{\pi}{2} \times \frac{\pi}{2}$ field of view. An illustration of cubemap projection can be seen in Figure 5.5(c). The cube faces have natural image statistics, however, several discontinuities arise when using such mappings, which generates misleading artifacts (i.e. fake corners, cut objects...).

Projecting the spherical image into a 2D plane induces significant distortions, that can
lead to misrepresentations and inaccuracies when processing the data. These distortions
arise due to the inherent differences in the geometry between the spherical surface and
the flat plane. As a consequence, researchers have proposed direct processing of the
sphere itself (Figure 5.5 (d)), bypassing the need for intermediary mapping techniques.
Nevertheless, as already mentioned, processing on the sphere is not a straightforward task
and requires a deep understanding of the underlying non-Euclidean geometry.

Moreover, even elementary processing operations, such as translation or convolution,
are not easily defined. In the following section, we briefly review different approaches that
have been used to perform convolution on omnidirectional images.
<figure><img src="image_37.png" /><figcaption>Figure 5.4 – Cubemap projection. Figure extracted from [108].</figcaption></figure>

**5.4 Convolution on omnidirectional images**

The implementation of spherical convolutions comes with many challenges. First, considering registered images, a convolution is desired to be rotation equivariant for any rotation around the polar axis, meaning that the convolution operation should remain consistent under rotations of the input signal around the north pole. Consequently, the convolution filters of the CNN can be applied consistently all over the sphere. Second, the filters should be able to increase their convolutional range and have anisotropic responses to ensure expressiveness. Finally, computational efficiency in spherical convolutions is very important for real-time processing and memory limitation.

A popular approach to process the spherical image is to map it on a 2D planar grid and apply traditional 2D CNN tools. Primitive existing strategies consider a CNN trained on perspective images and either (i) apply it directly to an equirectangular projection [113], [114] or (ii) apply it to the tangent planes that they sample all over the sphere [115]–[118]. While the former is inaccurate because of the distortions induced by the perspective projection, the latter solution is exact but computationally inefficient. Another solution proposed by Su et al. [119] is based on knowledge distillation where the network takes an equirectangular image as input and for every image position, learns to mimic the output result of applying the corresponding tangent plane to a CNN learned on perspective images. However, they learn kernels of different sizes for every row of the equirectangular projection to ensure that the kernels cover the same surface as the spherical image. This comes with a significant model size and memory footprint, which limits not only the training, but also the usage of the learned model on devices with constrained computational resources. The authors later propose the KTN [120], [121] which solves the problem of memory bloat while achieving high accuracy in object detection, yet still suffers from a significant computational cost because the kernels are generated at run time.

A second category of methods involves directly implementing convolution on the sphere. One approach withing this category involves using the spectral domain [122]–[125]. Spherical harmonics form a complete orthonormal basis on the surface of the
<figure><img src="image_38.png" /><figcaption>Figure 5.5 – (a) Image acquired with double fish-eye camera and its corresponding (b) equirectangular projection, (c) cubemap projection and (d) spherical representation. Figure extracted from [108].</figcaption></figure>

sphere. Hence, in the spectral domain, spherical convolutions can be performed by the
simple product of the spherical harmonics coefficients of the spherical data and the ker-
nel function. While the convolution is consistent and expressive, these methods are very
complex and have significant memory overhead, thus are not suitable for real-world high-
resolution omnidirectional images.

Another approach is to model the sphere as a graph and to approximate the convolu-
tion using a Chebyshev polynomial representation [126], [127]. Although the convolution
is rotation equivariant and consistent over the sphere, the learned filters are not expres-
sive because a single coefficient by neighborhood can be learned, therefore they can only
be isotropic due to the missing directionality in a graph. To address this drawback of
graph-based methods, multi-graph learning was proposed by Khasanova et al. [128] where
multiple directed graphs are used to represent different directions in order to maintain the
orientation data and ensure filter expressiveness. Nevertheless, a large storage is required
to construct all the graphs.

Finally, to overcome the limitations of existing methods, Bidgoli et al.[129] have pro-
posed a novel framework for spherical processing (OSLO) where they reach the same filter
expressiveness, consistency and complexity of 2D CNNs. For this purpose, they sample
the sphere using HEALPix [130] and define the convolution in the pixel domain. We
elaborate OSLO and HEALPix in the upcoming chapter, as they play a pivotal role in
shaping our contribution.

**5.5 Inverse problems in omnidirectional images**

In contrast to the extensive exploration of inverse problems in 2D imaging, the field
of inverse problems in omnidirectional images has received comparatively less attention.
Most of the existing work has focused on super-resolution, since omnidirectional images
must require very high resolution in order to ensure an acceptable quality of user experi-
ence. In this section, we summarize the different approaches that have been proposed to
solve omnidirectional inverse problems. In particular, we focus on super-resolution and
denoising.

**5.5.1 Super-resolution**

Early methods for super-resolving omnidirectional images used a sequence of low-resolution images and assembled them to get high resolution outputs [131]–[133]. With the rise of deep learning, Ozcinar et al. [134] designed a generative adversarial network to super-resolve equirectangular projections. Along with the adversarial loss, they use a latitude related loss function to deal with the distortion caused by the projection. Recently, a significant number of end-to-end networks have been introduced to tackle omnidirectional image super-resolution. Deng et al. [135] proposed a latitude adaptive super-resolution network named LAU-Net+ for ERP images, where different upscaling factors can be used for different latitude bands. Furthermore, among the emerging methods, we note TCCL-Net [136], OSRT [137] and OPDN [138]. Finally, it is worth highlighting that prevailing methodologies in super-resolution mainly center around the equirectangular projection.
5.5.2 Image denoising

Although omnidirectional images have gained a lot of attention in the past years, only a few works have been dedicated to denoising them. Bigot et al. [139] adapt the Wiener filter and Tikhonov regularization to spherical images for the denoising application. Demonceaux et al. [140] denoise catadioptric inputs corrupted by White Gaussian noise by proposing a new neighborhood for Markov random fields (MRF) adapted for omnidirectional images. Furthermore, Alibouch et al. [141] adapt the Stein block thresholding method to omnidirectional images to remove additive white Gaussian noise. Finally, Phan et al. [142] solve poisson denoising by using a space-variant total-variation regularization on catadioptric images. They adapt the fidelity and the regularization terms with weighted functions based on the mean curvature of an image surface.
# SPHERE DRUNET: A SPHERICAL DENOISER FOR OMNIDIRECTIONAL IMAGES

## 6.1 Introduction

From wavelet thresholding to advanced deep-learning architectures, image denoising is a classical yet still one of the most widely explored topics in image processing. Despite the massive work that has been done to denoise traditional 2D images, very little effort has been dedicated to omnidirectional image denoising.

Most of the existing methods use conventional techniques and are primarily tailored to address the challenges presented by catadioptric inputs, thus they depend on the nature of the acquisition.

We dedicate this chapter to omnidirectional image denoising, as denoisers play an important role in today's era of regularizing inverse problems. Furthermore, a commonly adopted approach for processing spherical images involves applying two-dimensional processing techniques to their corresponding mappings, typically using the equirectangular projection. We propose an approach that is not sensitive to sampling distortion, by working directly on the sphere, and we aim to show that image denoising gives better performance when it is performed directly on the sphere rather than via a mapping, even though it raises many challenges.

Moreover, we introduce a novel spherical CNN denoiser (SphereDRUNet) by transforming the architecture of the DRUNet [2] using the elementary convolutional modules proposed by OSLO [129]. Our network takes directly spherical inputs rather than equirectangular projections. We evaluate the performance of SphereDRUNet in comparison with the initial DRUNet retrained and applied on equirectangular data. Since AWGN is a noise per pixel degradation, one can think that denoising the ERP could yield satisfac-
tory results. However, our results show that working on the sphere leads to a significant
performance gain for the task of image denoising. This also confirms that the OSLO
solution can be successfully extended to inverse problems such as denoising. We further
compare our SphereDRUNet to a graph-based convolution network and confirm that the
convolution approach that we use is more efficient.

6.2 Related work

As previously covered in Chapter 5, traditional CNN operations such as convolution and translation can not be straightforwardly extended to the spherical domain. Certain requirements should be respected while constructing specific spherical CNN architectures: the convolution filters should be rotation equivariant, highly expressive and computationally efficient. Although different approaches for spherical convolutions have been proposed such as spectral learning, 2D processing and graph-based representations, it hasn't been possible to assure these 3 key properties simultaneously to guarantee a proper spherical convolution.

To overcome the drawbacks of prevailing methods, Bidgoli et al. [129] have proposed a novel framework for spherical processing (OSLO) where all aforementioned properties are successfully satisfied. For this purpose, they sample the sphere using HEALPix [130] and define the convolution in the pixel domain. We choose to build our work using OSLO because of its efficient ability to conduct spherical convolution. Next, we elaborate on HEALPix and OSLO.

6.2.1 HEALPix sampling

Hierarchical Equal Area isoLatitude Pixelation (HEALPix) [130] is a sampling method
for spherical data that produces a configuration where the pixels are arranged in a
diamond-shaped pattern. HEALPix starts by forming the base resolution and divides
the surface of the sphere into 12 equal-area regions, each representing a specific area on
the sphere. Then, each of these regions is iteratively partitioned into 2x2 equal-area sub-
pixels, until reaching the target resolution. The HEALPix *k*<sup>*th*</sup> resolution is defined by a
parameter *N*<sub>side</sub>=2<sup>*k*</sup> and has a total of *N*<sub>pix</sub>=12*N*<sub>side</sub><sup>2</sup> = 12 x 2<sup>*k*</sup> x 2<sup>*k*</sup> pixels (Figure 6.1).
Each pixel in the final arrangement has eight adjacent neighbors (only 24 pixels have seven
neighbors, which is far from having a significant effect taking into account the resolution
for an omnidirectional image). This ensures the regularity of the neighborhood, which is a fundamental condition for an expressive and effective convolution. Furthermore, the relative distance and orientation between neighboring pixels remains consistent throughout the entire sphere. Thus, HEALPix offers a rigid structure that was further proved in [129]. Both the regularity and the rigidity of the sampling method are essential properties to ensure an expressive and effective convolution.

<figure><img src="image_39.png" /><figcaption>Figure 6.1 – Visualization of HEALPix sampling over the Sphere. (a) Base resolution. (b) First resolution. (c) Third resolution. (d) Fourth resolution</figcaption></figure>

## 6.2.2 On-the-Sphere Learning for Omnidirectional Images (OSLO)

On-the-sphere Learning for Omnidirectional Images (OSLO) proposes a new framework that defines all necessary operators to build a CNN on the sphere. The method samples the sphere using HEALPix and performs convolutions in the pixel domain. Pooling, stride, skip-connections and patching are defined as well. The advantage of OSLO is that the convolutions are rotation equivariant, expressive and efficient at the same time.

In fact, let $N_i(k)$ be the index assigned to the $k^{th}$ neighbor for node $i$ ($k = 1, ..., 8$ for a total of 8 adjacent neighbors of a vertex). Let $L_{in}$ and $L_{out}$ denote respectively the number of input and output features in the convolution. The convolution operation for an output feature $l$ ($1 \le l \le L_{out}$) is defined as:

$$x_i^l = \langle \theta_0, x_i \rangle + \sum_{k=1}^{8} \langle \theta_k, x_{N_i(k)} \rangle \cdot w_{N_i(k),i} \quad (6.1)$$

where $\theta_k$ represents the learnable filter weight at $k$ and $x_i$ is the input data/feature at point $i$. $w_{N_i(k),i}$ is set to 0 when the neighbor $N_i(k)$ is missing and to 1 otherwise, in order to deal with the 24 exceptions of the pixels that have a missing adjacent neighbor.

The described convolution is anisotropic and consistent all over the sphere. The same
weights are applied to compute the convolution output regardless of the position of the kernel on the sphere (Figure 6.2). Yet, it only supports 1-hop neighborhood. In order to increase the potential size of kernels and extend the local support to n-hop neighborhood, the authors propose an iterative computation of the 1-hop convolution with a proper aggregation method such as concatenation, addition or max aggregation. This makes the convolution highly expressive. Furthermore, the convolution being in the pixel domain, it simply consists in translating the filter over the sphere, hence the complexity linearly increases with the number of pixels.

<figure><img src="image_40.png" /><figcaption>Figure 6.2 – Spherical convolution in the pixel domain proposed by [129]. Each color represents a weight corresponding to a certain orientation. Filter weights remain the same, regardless of the position of the kernel on the sphere.</figcaption></figure>

Moreover, the authors use this framework for the task of compression and confirm
that this on-the-sphere solution outperforms similarly learned models applied to equirect-
angular images.

For these reasons, we choose OSLO to work on the sphere for omnidirectional image denosing. Further information about the other operations can be found in [129].

**6.3 Spherical image denoising**

Our main objective is to explore whether image denoising is more efficient when per-
formed directly on the sphere, compared to using a projection. We choose the DRUNet[2]
CNN denoiser which is a state-of-the-art Gaussian denoiser for perspective images, and
use OSLO to transfer its architecture to the sphere. Hence, a novel spherical denoiser
<figure><img src="image_41.png" /><figcaption>Figure 6.3 – Architecture of the DRUNet network [2]. The network takes as input the noisy image concatenated with a noise level map.</figcaption></figure>

SphereDRUNet is obtained that can be used to denoise spherical images corrupted with
any noise level.

The network architecture of DRUNet is shown in Figure 6.3. DRUNet is a bias-free model, which has been observed to be effective for generalization of denoisers over un-seen noise levels. The model takes as input a noise level map and integrates Residual blocks (ResNet [96]) into U-Net [95]. It consists of four scales having 64, 128, 256 and 512 number of channels in each of the layers. Each scale in the downscaling consists of 4 residual blocks followed by a 2x2 strided convolution (SConv) whereas a 2x2 trans-posed convolution (TConv) is followed by 4 residual blocks during the upscaling. Four skip connections are set between 2x2 strided convolutions (SConv) downscaling and 2x2 transposed convolutions (TConv) blocks at each of the four scales.

We transfer DRUNet to the sphere using the OSLO framework. The architecture of the
resulting SphereDRUNet is shown in Figure 6.4. We use 1-hop neighborhood convolution
for all the layers because all convolutions in DRUNet are of size 3x3. 2D stride operations
of size 2x2 amount to stride 4 on the sphere. For implementation purposes, we replace
the 2x2 transpose convolution by the sub-pixel convolution of size 4 (pixel-shuffle) since
both unpooling methods have equivalent performance.
<figure><img src="image_42.png" /><figcaption>Figure 6.4 – Architecture of the SphereDRUNet network with OSLO-based convolution.</figcaption></figure>

### 6.3.1 Training details

For the purpose of a fair comparison between denoising a spherical image and its projection, we adapt our training framework to that of the DRUNet. We use the SUN360 equirectangular image dataset [143] and consider 2105 images for training, 22 images for validation and 43 images for testing. All images in the initial dataset have a spatial size of 9104x4552. For each epoch of the training, we degrade the training ERP images by AWGN with a random noise level $\sigma$ chosen from [0,50], and concatenate a uniform map filled with $\sigma$ having the same size as the ERP image. We perform HEALPix sampling of resolution 10 (i.e. 12,582,912 pixels) of the original image, the noisy image and the concatenation. SphereDRUNet takes as input the concatenated data in the spherical form, and learns to denoise the image by minimizing the $L_1$ loss between the estimated clean image and the ground-truth image using the ADAM optimization. The learning rate is set to 5e-5 and decreased by half every 50,000 iterations. The gradient update is performed once every 8 patches, patches being of resolution 8 ($2^8$ x $2^8$ pixels). Finally, we train our network for a total of 25 epochs. Note that most of the modifications to the training framework compared to that of the DRUNet are essentially driven by the higher resolution of spherical images in contrast to perspective images.

Furthermore, since our goal is to study denoising directly on the sphere compared to denoising the equirectangular image, applying the DRUNet trained on perspective images to ERP projections might not seem fair. Therefore, we re-train the DRUNet on ERP
images in the same way as we trained the SphereDRUNet. In order to reduce the unfair bias caused by the presence of ground-truth images in ERP format, we resize the ERP images to 5056x2528 (12,781,568 pixels) so that both ERP and HEALPix have almost the same number of pixels.

## 6.3.2 Comparison framework

As discussed earlier, we proposed SphereDRUNet, a novel spherical denoising CNN and we re-trained the DRUNet with ERP images in order to apply it on spherical projections. We want to compare the performance of these two approaches in order to study whether spherical denoising is more efficient on the sphere or on a corresponding projection. We also include a comparison with a cubemap projection: we map the equirectangular image to the six faces of a cube and each face is denoised with the DRUNet. We then reconstruct the ERP image by combining the six planes.

To avoid a biased comparison, for testing both spherical and mapping-based approaches, we first corrupt the high-resolution equirectangular input with a Gaussian noise. Then, we downsample the corrupted image to get a HEALPix sampling of resolution 10 (12,582,912 pixels with $N_{side} = 2^{10}$) and resize the ERP image to 5056x2528 (12,781,568 pixels) in order to match the resolution of the HEALPix sampling. Thus, both test images are generated from the same corrupted input.

For quantitative evaluation, we calculate the Spherical PSNR (S-PSNR) [dB] [144] and the Weighted to Spherically uniform PSNR (WS-PSNR) [dB] [145]. The former maps the output and the ground-truth images (ERP or sphere) to a sphere by uniformly sampling 655,362 points, and computes the mean error between them to simulate the PSNR on the sphere. The sampling process for the S-PSNR calculation is different than the HEALPix sampling, hence the metric is not biased. S-PSNR is seen to be an estimation of the overall quality experienced by viewers across all potential views. On the other hand, in the calculation of the WS-PSNR, the mean squared error is weighted by the size of each pixel, as follows:

$$WS-PSNR = 10 \cdot \log_{10} \left( \frac{\text{MAX}_I^2}{\text{W-MSE}} \right), \qquad (6.2)$$

where MAXI is the maximum pixel intensity in the image, and W-MSE is given as:

$$\text{W-MSE} = \frac{\sum_{x=0}^{W-1} \sum_{y=0}^{H-1} ((\hat{I}(x, y) - I^{gt}(x, y))^2 \cdot w(y)}{\sum_{x=0}^{W-1} \sum_{y=0}^{H-1} w(y)} \qquad (6.3)$$
$W$ and $H$ denote the width and height of the image, $\hat{I}$ and $I^{gt}$ represent respectively the reconstructed and the ground-truth images and $w(y)$ represents the weight of the corresponding pixel in the weight matrix. We also show mollweide projections of the denoised images for some visual comparison.

## 6.4 Experimental results

We perform spherical denoising on images corrupted with an Additive White Gaussian noise (AWGN) of noise level $\sigma * 255 \in [10, 20, 30, 40, 50]$. Tables 6.1 and 6.2 show respectively the WS-PSNR [dB] and the S-PSNR [dB] of denoised equirectangular, cubemap and spherical images of the testing dataset. We can see that the on-the-sphere solution with our SphereDRUNet significantly outperforms denoising the equirectangular and the cubemap projections. This gain is even more important when the noise level is high. Figures 6.5, 6.6 and 6.7 show visual comparisons of the denoised ERP and spherical images for degradations of noise level $\sigma= 30/255$ and $\sigma= 50/255$. We observe that the DRUNet smoothes the equirectangular images and can not recover the high frequency details while the images denoised using SphereDRUNet on the sphere have less distortion and do not lose information (i.e. recover high frequency details such as texture). For instance, in Figure 6.5, we observe that the denoised equirectangular image has lost high frequency elements such as the details on the leaves or the stains on the train, whereas denoising the sphere successfully removes noise while also recovering details. This is probably because in order to recover an equirectangular image with spatial distortions, a network needs a larger amount of weights to reconstruct high frequency components compared to recovering the same high frequency element on the sphere. For instance, a line on the sphere remains a line regardless of its position, whereas it can take different forms on the ERP depending on where it lies, so a network needs more weights to recover the same high frequency components on the ERP compared to the sphere.

While the pixel distribution in the CMP is more uniform than that in the ERP, the mapping between the equirectangular projection and the cubemap projection induces a distortion and a significant shift in the pixels, especially because of the border discontinuities. We show this in Figure 6.8 where we map an ERP image into a CMP image (since the SUN360 dataset is in ERP format) and simply convert it back to ERP. A degraded image with an important shift is observed, with a WS-PSNR of 24.95 dB between the original ERP image and the re-converted projection. This further underlines the advantage