Anna Jezierska

Audrey Repetti

Mai Quyen Pham

Feriel Abboud

Yosra Marnisi

Emilie Chouzenoux

Aurélie Pirayre

Mounir Kaaniche

Florian Dupuy

Olena Tankyevych

Camille Couprie

Ravi Kiran

Anisia Florescu

Giovanni Petrazzuoli

Without

Patrice Herault

Aida Arida

Pascale Asmar

Mashhour Mohammad El Achkar

El Chakouch

Ahmad Diab

Karam Fayad

Rim Kaddah

Ahlem Khlass

Karim Youssef

Pierre Youssef

Kifah Sarraf

M El Gheche

C Chaux

J.-C Pesquet

B Pesquet-Popescu

J Farah

Convex optimization aims at searching for the minimum of a convex function over a convex set. While the theory of convex optimization has been largely explored for about a century, several related developments have stimulated a new interest in the topic. The first one is the emergence of efficient optimization algorithms, such as proximal methods, which allow one to easily solve large-size nonsmooth convex problems in a parallel manner. The second development is the discovery of the fact that convex optimization problems are more ubiquitous in practice than was thought previously.

In this thesis, we address two different problems within the framework of convex optimization. The first one is an application to computer stereo vision, where the goal is to recover the depth information of a scene from a pair of images taken from the left and right positions. The second one is the proposition of new mathematical tools to deal with convex optimization problems involving information measures, where the objective is to minimize the divergence between two statistical objects such as random variables or probability distributions.

We propose a convex approach to address the problem of dense disparity estimation under varying illumination conditions. A convex energy function is derived for jointly estimating the disparity and the illumination variation. The resulting problem is tackled in a set theoretic framework and solved using proximal tools. It is worth emphasizing the ability of this method to process multicomponent images under illumination variation. The conducted experiments indicate that this approach can effectively deal with the local illumination changes and yields better results compared with existing methods.

We then extend the previous approach to the problem of multi-view disparity estimation. Rather than estimating a single depth map, we estimate a sequence of disparity maps, one for each input image. We address this problem by adopting a discrete reformulation that can be efficiently solved through a convex relaxation. This approach offers the advantage of handling both convex and nonconvex similarity measures within the same framework. We have shown that the additional complexity required by the application of ix x Abstract our method to the multi-view case is small with respect to the stereo case.

Finally, we have proposed a novel approach to handle a broad class of statistical distances, called ϕ-divergences, within the framework of proximal algorithms. In particular, we have developed the expression of the proximity operators of several ϕ-divergences, such as Kulback-Leibler, Jeffrey-Kulback, Hellinger, Chi-Square, I α , and Renyi divergences. This allows proximal algorithms to deal with problems involving such divergences, thus overcoming the limitations of current state-of-the-art approaches for similar problems. The proposed approach is validated in two different contexts. The first is an application to image restoration that illustrates how to employ divergences as a regularization term, while the second is an application to image registration that employs divergences as a data fidelity term.

Resumé

Cette thèse s'inscrit dans le contexte de l'optimisation convexe. Elle apporte à ce domaine deux contributions principales. La première porte sur les méthodes d'optimisation convexe non lisse appliquées à la vision par ordinateur. Quant à la seconde, elle fournit de nouveaux résultats théoriques concernant la manipulation de mesures de divergences, telles que celles utilisées en théorie de l'information et dans divers problèmes d'optimisation.

Le principe de la stéréovision consiste à exploiter deux images d'une même scène prises sous deux points de vue, afin de retrouver les pixels homologues et de se ramener ainsi à un problème d'estimation d'un champ de disparité. Dans ce travail, le problème de l'estimation de la disparité est considéré en présence de variations d'illumination. Ceci se traduit par l'ajout, dans la fonction objective globale à minimiser, d'un facteur multiplicatif variant spatialement, estimé conjointement avec la disparité. Nous avons mis l'accent sur l'avantage de considérer plusieurs critères convexes et non-nécessairement différentiables, et d'exploiter des images multicomposantes (par exemple, des images couleurs) pour améliorer les performances de notre méthode. Le problème d'estimation posé est résolu en utilisant un algorithme parallèle proximal basé sur des développements récents en analyse convexe.

Dans une seconde partie, nous avons étendu notre approche au cas multi-vues qui est un sujet de recherche relativement nouveau. Cette extension s'avère particulièrement utile dans le cadre d'applications où les zones d'occultation sont très larges et posent de nombreuses difficultés. Pour résoudre le problème d'optimisation associé, nous avons utilisé des algorithmes proximaux en suivant des approches multi-étiquettes relaxés de manière convexe. Les algorithmes employés présentent l'avantage de pouvoir gérer simultanément un grand nombre d'images et de contraintes, ainsi que des critères convexes et non convexes. Des résultats sur des images synthétiques ont permis de valider l'efficacité de ces méthodes, pour différentes mesures d'erreur.

xi xii

Resumé

La dernière partie de cette thèse porte sur les problèmes d'optimisation convexe impliquant des mesures d'information (ϕ-divergences), qui sont largement utilisés dans le codage source et le codage canal. Ces mesures peuvent être également employées avec succès dans des problèmes inverses rencontrés dans le traitement du signal et de l'image. Les problèmes d'optimisation associés sont souvent difficiles à résoudre en raison de leur grande taille.

Dans ce travail, nous avons établi les expressions des opérateurs proximaux de ces divergences. En s'appuyant sur ces résultats, nous avons développé une approche proximale reposant sur l'usage de méthodes primales-duales. Ceci nous a permis de répondre à une large gamme de problèmes d'optimisation convexe dont la fonction objective comprend un terme qui s'exprime sous la forme de l'une de ces divergences.

Vue d'ensemble du travail effectu é

Optimisation convexe Le chapitre 2 a pour but d'introduire le lecteur aux problèmes inverses dans le cadre convexe. Notre exposé se limite à la description des outils mathématiques nécessaires à notre étude. Pour ce faire, nous avons fourni, dans un premier temps, un bref aperçu des problèmes inverses en nous basant sur des approches variationnelles, pour ensuite présenter les algorithmes d'optimisation convexe. Une attention particulière a été accordée à la projection sur des ensembles convexes, à la fin de ce chapitre.

Estimation de la disparité La stéréoscopie est l'un des problèmes fondamentaux de la vision par ordinateur. Le procédé de la stéréoscopie est calqué sur la perception humaine du relief grâce aux deux images planes que perçoit chaque oeil. Ce processus constitue un enjeu majeur en vision par ordinateur. Une étape cruciale dans la résolution de ce problème est la mise en correspondance stéréoscopique.

Le problème de la mise en correspondance est un procédé de reconstruction de la profondeur à partir de deux images de la même scène acquises avec des angles différents. En établissant que deux points correspondants, dans deux images acquises, représentent la projection d'un même point de la scène, on peut calculer par triangulation la position tridimensionnelle exacte de ce point dans l'espace.

L'un des problèmes les plus délicats en vision par ordinateur, notamment dans le cadre de la mise en correspondance stéréoscopique, est Resumé xiii celui de changements d'illumination qui peuvent affecter deux vues d'une même scène. Un autre problème critique provient de la présence d'objets dont la surface est uniforme ou de la présence de textures répétitives. En effet, les pixels dans ces régions ne représentent pas de primitives discriminantes, c'est à dire que plusieurs correspondants potentiels ayant la même valeur de luminance existent dans l'autre image.

Les régions homogènes sont donc difficiles à apparier. Par ailleurs, une autre difficulté rencontrée par les méthodes de mise en correspondance concerne les occultations. Des points de la scène peuvent être visibles dans une image de la paire stéréoscopique mais pas dans l'autre.

Pour pallier ces difficultés, les méthodes de mise en correspondance sont amenées à exploiter toutes les informations disponibles afin de faciliter la recherche et la détermination des correspondants. Les techniques les plus simples et les plus courantes se basent sur une recherche locale de similarité entre de petites régions, typiquement par optimisation d'une fonction de corrélation. On peut aussi segmenter les images et mettre en correspondance les régions obtenues. L'avantage de ces méthodes est que la mise en oeuvre est simple et efficace pour un temps de calcul relativement faible. Cependant, les techniques locales sont sensibles aux occultations, aux régions faiblement texturées et aux motifs répétitifs. Les approches globales tentent de pallier ces inconvénients en minimisant une fonction de coût globale qui fait intervenir l'ensemble des pixels des images, mais au prix d'une complexité de calcul plus élevée. Plusieurs méthodes d'optimisation ont été proposées dans la littérature. Nous pouvons distinguer les méthodes à variables discrètes, qui effectuent un appariement de couples de pixels, comme la programmation dynamique ou les coupures de graphes, des méthodes à variables continues, où la disparité estimée prend des valeurs réelles "quelconques", comme dans les approches variationnelles.

Dans notre travail, nous nous plaçons dans ce dernier cadre en employant des méthodes parallèles reposant sur des opérateurs proximaux. Dans le Chapitre 3, le problème de mise en correspondance est alors formulé comme un problème d'optimisation convexe sous contraintes. Une fonction objective convexe est minimisée sur l'intersection d'ensembles convexes. Ces ensembles sont associés à différentes contraintes modélisant des informations a priori sur les propriétés des champs à estimer. L'originalité de cette approche se traduit par plusieurs aspects:

• L'estimation de la carte de disparité est effectuée ici en présence de variations d'illumination dans la scène observée. Ceci se traduit par l'ajout, dans la fonction objective globale à minimiser, d'un xiv Resumé facteur multiplicatif variant spatialement, qui est estimé conjointement avec la disparité.

• Des images multi-composantes sont prises en considération. Cette technique est élaborée en prenant en compte les trois canaux couleurs (RGB, YUV, . . . ). La formulation globale n'étant pas convexe, on propose une méthodologie d'approximation convexe basée sur un développement de Taylor au premier ordre.

• Le critère convexe (distance) considéré peut être choisi de manière très générale et peut être non différentiable (une norme 1 par exemple).

• L'algorithme considéré offre une grande flexibilité quant au choix du critère à minimiser et à l'incorporation de contraintes multiples. Par ailleurs, il faut souligner que la méthode peut être implantée sur une architecture parallèle (multi-coeurs, multi-processeurs, GPU, ...).

Dans les résultats de simulation fournis, nous démontrons l'efficacité de notre approche.

Estimation de la disparité en multi-vues L'imagerie multi-vues a suscité un vif intérêt dans différentes applications telles que le cinéma numérique, la médecine, le trafic aérien, les technologies militaires, les jeux vidéos, etc. Le développement de la télévision 3D et des supports auto-stéréoscopiques a suscité de nombreuses études portant sur le traitement de séquences d'images. Il s'agit en particulier, d'extraire de l'information 3D à partir des vues disponibles, sous forme de cartes de disparités.

Dans le Chapitre 4, nous étendons l'approche présentée dans le Chapitre 3 au cas multi-vues. On dispose de caméras calibrées, dont les paramètres sont connus et les images acquises sont rectifiées. Nous proposons d'estimer les cartes de disparité de chaque image présente dans la séquence en exploitant toutes les données dans les vues adjacentes. La stratégie adoptée permet d'exploiter les zones occultées qui se trouvent dans certaines vues, ainsi que d'améliorer la précision des estimations (par rapport au cas stéréo).

Pour résoudre les problèmes d'optimisation associés, nous avons utilisé des algorithmes proximaux en suivant deux approches de relaxation convexe différentes:

• Premièrement, nous avons considéré une approche multi-étiquettes fondée sur des développements récents de l'analyse convexe. La technique employée présente l'avantage de pouvoir manipuler un grand nombre d'images avec une complexité faible par rapport Resumé xv au cas stéréo. Cet algorithme offre une grande flexibilité quant à l'incorporation de critères convexes mais aussi non convexes ( 1 -tronquée, 1 2 -tronquée, . . . )

• Deuxièmement, nous appliquons une linéarisation basée sur un développement de Taylor amélioré.

Des résultats sur des images ont permis, pour différentes formes de critère de mesure d'erreur, de valider l'efficacité de la première méthode par rapport à la seconde et d'autres approches existantes.

Opérateurs proximaux des ϕ-divergences

La théorie de l'information est une branche émergente des mathématiques appliquées et de l'informatique, issue de la volonté de quantifier les échanges d'information. Dans les années 40, Claude Shannon a développé la théorie de l'information qui décrit les aspects les plus fondamentaux des systèmes de transmission. Cette théorie s'intéresse à la construction et à l'étude de modèles mathématiques, en se basant essentiellement sur la théorie des probabilités. Elle s'est développée au fur et à mesure et est devenue aujourd'hui incontournable dans la conception des systèmes de communication. Elle a été étendue à d'autres champs d'applications, ou le codage, comme pour le codage source/canal.

Une fonction fondamentale de la mesure de l'information est la divergence de Kullback -Leibler, qui mesure la distance entre deux objets statistiques telles que des variables aléatoires par le biais de leurs distributions de probabilité. Une telle divergence est un cas particulier d'une classe plus large de distances appelées ϕ-divergences, incluant celles de Jeffreys-Kullback, Hellinger, Chi-Square, I α , et de Renyi. Nous présentons dans ce qui suit un bref résumé sur l'état de l'art concernant ces divergences:

• La divergence de Kullback-Leibler(KL) a un rôle de premier plan dans le domaine de codage source/canal pour le calcul de la capacité de transmission et de la fonction de débit-distorsion. Ces problèmes sont souvent résolus par des algorithmes basés sur la minimisation alternée proposés par Blahut et Arimoto. La symétrisation classique de Kullback-Leibler est la divergence de Jeffreys-Kullback, qui est connue dans plusieurs applications surtout dans les méthodes basées sur les algorithmes des k-moyennes.

• La divergence d'Hellinger a été initialement introduite par Beran et elle est connue sous différentes dénominations (telles que la xvi Resumé divergence de Jeffreys-Matusita). Dans le domaine de la théorie de l'information, la divergence d'Hellinger est couramment utilisée dans l'analyse des données, des statistiques et de l'apprentissage automatique.

• La divergence Chi-carré a été introduite par Pearson qui l'a utilisée pour évaluer quantitativement si un phénomène observé tend à confirmer une hypothèse donnée. Cette divergence a été appliquée avec succès dans différents contextes, tels que la théorie de l'information et le traitement du signal, comme une mesure de dissimilarité entre deux distributions de probabilités.

Les méthodes proximales constituent un cadre approprié pour les ϕdivergences, puisque ce sont des méthodes permettant de traiter des fonctions convexes non lisses. La difficulté provient du fait que le calcul de l'opérateur proximal des ces divergences n'est pas trivial. En effet, une ϕ-divergence est une somme de fonctions à deux variables non séparables. Par conséquent, le défi consiste à calculer des expressions explicites de ces opérateurs. Dans le Chapitre 5, nous avons établi les expressions des opérateurs proximaux de ces divergences. Ce calcul nous a permis de développer de nouvelles méthodes proximales qui permettent d'aborder des formes plus générales de problèmes d'optimisation. Les expressions des opérateurs proximaux présentées, directement liée à des fonctions convexes de deux variables réelles, ont enrichi la liste des fonctions pour lesquelles un tel opérateur proximal peut être facilement calculé. Optimization is a rich and thriving mathematical discipline. Historically, the first methodology for solving nonlinear optimization problems is due to Fermat (1601-1665). Because of this work, Lagrange (1736-1813) stated clearly that he considered Fermat to be the inventor of calculus (as opposed to Newton (1643-1727) and Leibnitz (1646-1716) who were later locked in a bitter struggle for this honor). At the beginning of the nineteenth century, the least-squares method was firmly established as a scientific approach by Gauss, Legendre, and Laplace within the space of a decade. Since then, optimization theory has grown ever more sophisticated. The powerful and elegant language of convex analysis unifies much of this theory, to which a strong contribution was brought, among others, by Moreau (1923Moreau ( -2014)). A large number of important analytic problems have illuminating optimization formulations and hence can be approached through convex analysis tools. In this thesis, we present two main contributions within the framework of convex optimization. The first one concerns an application of recent nonsmooth optimization methods to computer stereo vision, while the second one provides novel theoretical tools to deal with divergence measures in optimization problems.

Disparity estimation One of the main challenges in computer vision consists of reproducing human abilities by means of mathematical models involving geometry, statistics, and learning theory. A popular example is computer stereo vision, inspired from a biological process called stereopsis. Stereopsis is the impression of depth that is perceived when a scene is observed with two eyes. Indeed, binocular vision creates two slightly different images of the scene in the two eyes, due to the different positions of the eyes on the head. These differences, referred to as binocular disparity, provide information that the brain can use to calculate depth in the visual scene, providing a major means of depth perception.

Introduction

In stereo vision, binocular disparity refers to the difference in location of an object seen from two different viewpoints. Such a difference results from the horizontal separation of cameras (see Figure 1.1). In the same way as the brain computes the depth of the visual scene, the disparity serves to extract 3D information from 2D digital images, since the disparity is inversely proportional to the depth. Disparity map can be recovered by estimating the relative position of features in the stereo pair. This process is called stereo matching. Looking for correspondences between stereo images is a difficult task, because of the presence of hidden areas (i.e. occlusions) and because of the fact that the light is reflected differently depending on the viewing angle. A huge amount of solutions have been proposed to address stereo matching. In this thesis, we present original approaches based on convex optimization. Divergence measures Information theory is a branch of applied mathematics and computer science involving the quantization of information. It was developed by Shannon to find fundamental limits on communication and signal processing operations. Since its inception, it has broadened to find applications in many other areas, including statistical inference, natural language processing, cryptography, neurobiology, quantum computing, and other forms of data analysis. Important sub-fields of information theory are source/channel coding, algorithmic information theory, complexity theory, and measures of information (see Figure 1.2). A key measure of information is the Kullback-Leibler divergence, which quantifies the distance between two statistical objects such as random variables or probability distributions. Such a divergence is a special case of a broader class of statistical distances called ϕ-divergences. Thanks to their interesting properties, ϕ-divergences have been widely used in optimization problems as similarity criteria between the sought solution and a given observation. However, existing approaches involving these kinds of measures are often restricted to problems where the minimization is performed w.r.t. one of the arguments of the divergence, or by using alternating minimization approaches which require specific assumptions to be valid. In this thesis, we present an original convex optimization approach that circumvents the above limitations.

(p) = 1 1-α ( i p α i -1) Tα (p||q) = 1 1-α (1 - p α q α -1 )
Bregman divergences (1967) BF (p||q) = F (p) -F (q)-< p -q, F (q) > B r egm an-C siszár divergence (1991)

F α ( x ) = x -l o g x -1 α = 0 x l o g x -x + 1 α = 1 1 α ( 1 -α ) ( -x α + α x -α + 1 ) 0 < α < 1 C siszár ' f -divergence D f (p||q) = pf ( q p ) (Ali& Silvey 1966, Csiszár 1967) Amari α-divergence (1985) f α ( x ) = x l o g x α = 1 -l o g x α = -1 4 1 -α 2 ( 1 -x 1 + α 2 ) -1 < α < 1
Quant um ent ropy S(ρ) = -kTr(ρ log ρ) (Von Neumann 1927)

K olmogorov K (p||q) = |q -p| ( K o l m o g o r ov -Sm i r n o f f m a x | p -q | ) Hellinger H (p||q) = ( √ p - √ q) 2 = 2(1 - √ f g Cherno divergence (1952) Cα (p||q) = -ln p α q 1-α χ 2 t est χ 2 (p||q) = ( q-p

Introduction

Challenges

Estimating depth from stereo images is a central challenge in computer vision and pattern analysis. Recovering 3D information from a pair of stereo cameras has been a popular topic because the additional 3D information provided by this technology contains significantly more information than 2D clues produced by traditional cameras. However, this 2D to 3D evolution has always been facing many challenges, which can be grouped into two main categories: accuracy and efficiency. Accuracy becomes an important concern in applications such as precise 3D surface modeling, especially when dealing with object surfaces with reflectance behavior, rich geometric structure, significant amount of occlusion and poor texture. Efficiency is one of the main issues when the stereo system is employed in real-time applications such as robot navigation, video surveillance, and interactive user interfaces. Unfortunately, these challenges often conflict with each other: in order to improve the quality of stereo matching, people usually cast the problem as a global optimization problem, which results in a high computational cost and a poor efficiency. On the other hand, most efficient stereo matching algorithms are based on only local information, which leads to poor accuracy in some difficult situations. Another challenge in depth estimation from stereo images is the presence of hidden surfaces (occlusions), moving objects or surfaces that reflect the light differently depending on the viewing angle, therefore complicating the matching task.

The rapid development of computer technology in the last decades has changed the methodologies underlying the 3D reconstruction. In this context, a possible approach consists of first extracting the depth from multiple images, and then reconstructing 3D scene using the depth information. Therefore, multi-view depth estimation has become a crucial step in computer vision. However, this kind of methods have to face the challenge of efficiently handling a massive amount of multi-view data (w.r.t. the stereo case).

ϕ-divergences are often used as discrete measures in signal processing problems. They serve as dissimilarity functions in many information theoretic models, data recovery tasks, and machine learning. Their popularity stems from their intimate connections with information theory concepts such as the entropy and the mutual information. Being dissimilarity functions, they may serve as criteria to be minimized in optimization problems. However, ϕ-divergences are nonsmooth functions, which rules out smooth optimization techniques. Proximal methods represent an ideal framework for ϕ-divergences, as these methods are able to deal with nonsmooth convex functions. Here, the difficulty stems from the fact that computing the proximity operator of ϕ-divergences is not trivial since they are not separable functions of scalar variables. Consequently, the challenge is to derive closed form expressions of these operators.

Contributions

Applications to computer stereo vision:

• We propose a convex optimization framework to address the problem of dense disparity estimation under varying illumination conditions. A convex energy function that takes into account the illumination variation model is derived by resorting to a relaxation based on a first-order Taylor approximation around an initial estimate. This energy is then minimized while taking into consideration various convex constraints arising from prior knowledge and observed data (Chapter 3).

The proposed method allows us to incorporate various convex distances and it relies on the extension of the Parallel ProXimal Algorithm (PPXA).

• We show, for the same optimization problem, the ability to consider multicomponent images with illumination variation, and also the flexibility in minimizing various (possibly nonsmooth) convex similarity measures and combining them with various convex constraints.

• In the context of multi-view images, we address the problem of disparity estimation from rectified images. We estimate a sequence of disparity maps by exploiting, for each of these maps, information from all the available views.

• We investigate the possibility of handling nonconvex similarity measures by adopting a discrete formulation of the disparity estimation problem which is exactly solved by a tight convex relaxation (Chapter 4). We show that the additional complexity incurred by our method in the multi-view context is small with respect to the stereo case.

Theoretical tools for dealing with divergence measures in optimization:

• ϕ-divergences are often restricted to problems where the minimization is performed w.r.t. one of the arguments of the divergence, or they are based on an alternating minimization process which requires specific assumptions to converge. In this work, we develop a novel proximal method that allows us to address more general forms of optimization problems, by computing the expressions of the proximity operators of two-variable convex functions. The expressions we have derived for these divergences enrich the list of functions for which such proximity operators can be easily calculated (Chapter 5).

• We have provided two different applications for divergence proximity operators in signal processing:

We extended the Non Local Total Variation (NLTV) framework by using information divergences to build new sparsity regularization measures for signal recovery.

We have involved the divergence in the expression of the data fidelity term for image registration, more precisely for disparity estimation under illumination variation.

Organization and Summary

This thesis is organized as follows. Chapter 2 introduces the reader to the topic of convex optimization and the associated algorithmic solutions.

Then, Chapter 3 addresses the dense disparity estimation problem from multicomponent images, formulated as an inverse problem by minimizing a convex objective function under multiple convex constraints. These constraints arise from prior knowledge and rely on various properties of the field to be estimated. We address the problem of disparity estimation under varying illumination conditions, and develop a spatially varying multiplicative model that accounts for photometric changes between both images in the stereo pair. Results on synthetic and real stereo pairs demonstrate the efficiency of the proposed method in recovering illumination changes and disparity map jointly, making disparity estimation very robust w.r.t. such changes.

For large scale or complex scenes, we address, in Chapter 4, the problem of multi-view disparity estimation from rectified images. A sequence of disparity maps is estimated by using two different approaches, which differ in the type of convex relaxation employed. For each of these maps, information from all the available views are exploited without adding much computational cost with respect to the two-view case.

In Chapter 5, we focus on convex optimization problems involving information measures (ϕ-divergence), which have been extensively investigated in source and channel coding. These measures can also be successfully used in inverse problems encountered in signal and image processing. The related Introduction optimization problems are often challenging due to their large size. We derive closed-form expressions of the proximity operators of these divergences. Building upon these results, we develop an efficient proximal approach. This allows us to address, in Chapter 6, two applications of convex optimization problems involving one of these divergences.

Finally, we draw some conclusions and perspectives in Chapter 7.

-Chapter 2 -

Convex optimization background

"Everything is possible convex." P. L. Combettes and J.-C. Pesquet Optimization constitutes a wide and rich mathematical field. Many applications in computer design or analysis, in particular for real-time reactive or automatic systems, rely on a wealth of these thriving tools. In recent years, the theory underlying the computational techniques in convex and nonconvex optimization have been studied by many authors. However, the elegant and dynamical language of convex analysis enlights much more the potential of optimization theory. Convex analysis, its extensions and its applications have known many developments in the recent years. This chapter mainly focuses on the convex formulation of inverse problems. In this part, we aim at describing the tools relevant to our studies. Firstly, we begin with a brief overview of inverse problems addressed by variational approaches. Then, we provide a brief presentation of convex optimization algorithms. Finally, we discuss different methods for computing the projection onto a convex set. 10 Chapter 2. Convex optimization background § 2.1 inverse problems In many optimization problems, one is interested in measuring a physical signal which is not directly available. In this context, we often classify the problems along two axes:

• Direct problem: having a complete description of a physical system, we can estimate the outcome of some measurements. One can say it is a problem of inferring the consequences from given causes.

• Inverse problem: finding the unknown causes of known consequences (it is associated with the reversal of the cause-effect relationship).

Using the actual result of some measurements, we are able to infer the values of the parameters that characterize the system Bertero and Boccacci [1998].

The inherent inverse problem consists in recovering the target signal by making explicit any available a priori information on the model parameters.

There exist many applications to inverse problems which can be found in many fields such as satellite image restoration A popular inverse problem consists in recovering a signal (e.g. an image) which is as close as possible to a reference one. Because of sensor imperfections and acquisition errors, the measured data are often noisy and degraded by a linear operator. Let us consider the basic approach in restoration where the direct problem can be formulated as follows:

z = T x + w, (2.1)
where z = (z j ) 1≤j≤K ∈ R K is the degraded observation and x = (x i ) 1≤i≤N ∈ R N is the image to be recovered. The linear operator T = (T j,i ) 1≤j≤K,1≤i≤N ∈ R K×N describes the physical laws linking x to the measurements z (e.g. motion blur, defocussing blur, ...) and the components of w = (w j ) 1≤j≤K ∈ R K are realizations of a mutually independent zero-mean (for example Gaussian) noise process. Solving the inverse problem associated with (2.1) aims at recovering x from z. One of the most well-known approaches in this context is the least squares method [START_REF] Gauss | Theory of the Combination of Observations Least Subject to Errors[END_REF]; [START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF]. It consists in minimizing the quadratic distance between the observation and the image to be recovered: argmin

x∈R N T x -z 2 .
(2.2) The solution to Problem (2.2) may be highly unstable and sensitive to small changes in the data. In such a case, we can certainly see why an inverse problem would be "difficult". A more precise way to look at this might be defined by how "solvable" a given problem is, which leads to the notion of well-posed and ill-posed problems proposed by [START_REF] Hadamard | Sur les problèmes aux dérivés partielles et leur signification physique[END_REF].

Well-posed problems Hadamard proposed three properties that a problem must possess in order to be classified as a well-posed one: existence, uniqueness and stability. For our purposes, considering only finite dimensional linear operators, we can think of data and solutions as the input and output vectors of some linear transformations. The first two properties seem rather obvious: a solution x exists for each observation z and this solution is unique. The third property is a question of stability, requiring that small changes to the input do not produce arbitrarily large changes to the output. A problem that lacks any one of these properties is by definition ill-posed. Thus, to solve the problem in a reliable manner, it is necessary to incorporate additional information Tikhonov [1963]. One of the first approaches to get a stable and satisfactory solution consists in adding a quadratic penalization term in order to achieve a small energy:

argmin x∈R N T x -z 2 + λ x 2 , (2.3)
where λ > 0 is the regularization parameter which balances the solution between the data fidelity term T x -z 2 and the penalization x 2 . A more general formulation of the regularized problem will be presented afterwards.

Data fidelity and prior information

Since many inverse problems are inherently ill-posed, a Bayesian approach may help in finding a stable and satisfactory solution. Let x be a realization of a random variable X with prior probability density function P X and z a realization of a random variable Z with conditional likelihood P Z|X=x . One aims at finding the signal x having the maximum a posterior probability density function P X|Z=z (x), given the observation z, i.e.

x * = arg max x P X|Z=z (x).

(2.4)

Using Bayes rule, the posterior probability density function can be decomposed into likelihood and prior probability measures :

x * = arg max

x P Z|X=x (z) P X (x).

(2.5)

Taking the minus-logarithm of the above expression, eventually leads to the following minimization problem:

minimize x -log P Z|X=x (z) g(x)
+ (-log P X (x)

Ψ(x) ), (2.6)
where the function g is the data fidelity term and the function Ψ is the prior term. Note that the fidelity term typically presents a separable form, as it is common to assume the statistical conditional independence of the components of vector Z:

g(x) = N i=1 g (i) (x, z i ).
(2.7)

The terms g (i) (x, z i ) may be explicitly defined by a particular likelihood distribution P Z i |X=x (z i ) (as in denoising problems), or they may be implicitly derived from local matching scores (as in disparity estimation problems). From a practical standpoint, it is important to balance the importance of the data fidelity term g against the regularization Ψ, as it has an impact over the signal to be recovered. A possible solution is to introduce a regularization constant λ ∈]0, +∞[ into the minimization problem, leading to the following unconstrained formulation: minimize x g(x) + λΨ(x).

(2.8)

The alternative solution is to enforce an upper bound on the function g or Ψ, yielding the following alternative constrained formulations:

minimize x g(x) s. t. Ψ(x) ≤ τ (resp. minimize x Ψ(x) s. t. g(x) ≤ η), (2.9)
where τ and η are some positive constants which are chosen depending on some prior information. Indeed, under technical assumptions, both formulations are equivalent for some specific values of λ and τ (resp. η).

However, the constrained formulation may be considered more practical, as the solution is less sensitive to the choice of regularization parameter λ, and even the choice of τ (resp. η) may be easier than that of λ Convex data fidelity term In general, the function g may be either convex or nonconvex. Unfortunately, finding a global minimum to Problems (2.8) and (2.9), when g is nonconvex, turns out to be a very hard problem. In some cases, good results may be obtained by local minimization, starting from a good initialization, using for example PDEs (Partial Differential Equations) or other variational approaches. However, such methods cannot guarantee any form of quality of the result, as they may get stuck in local minima.

On the other hand, when g is convex, there exist efficient numerical methods to solve the above problems with guaranteed convergence (under weak conditions). Consequently, a variety of convex relaxation techniques have been developed to make the function g convex: first order Taylor approximation [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], convex-hull approximation Jiang et al. [2007], exact relaxation for binary problems [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF], exact relaxation for multi-label problems [START_REF] Pock | Global solutions of variational models with convex regularization[END_REF], and semidefinite programming [START_REF] Lieven | Semidefinite programming[END_REF].

Convex regularization:

The aim of regularization is to convey some prior knowledge about the signal to be recovered, in order to make the original problem well-posed. A priori knowledge should not reflect any specific measurements of the physical signal which is currently available and analyzed. In this regard, the more carefully the regularity is modeled (in terms of energy to be minimized), the better the quality of the estimated signal. Since the signals of interest are often piecewise smooth, popular regularization models tend to enforce a sparsity condition by penalizing an analysis representation of the target signal. In this context, Total Variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]; [START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF]; [START_REF] Peyré | Group sparsity with overlapping partition functions[END_REF]; [START_REF] Bayram | Directional total variation[END_REF] has emerged as a simple, yet successful, regularization, consisting in penalizing the gradient coefficients. However, TV fails to preserve textures, details, and fine structures, because they are hardly distinguishable from noise. To improve this behaviour, the TV model has been extended by using a non-locality principle [START_REF] Bibliography | Nonlocal operators with applications to image processing[END_REF]; [START_REF] Werlberger | Motion estimation with non-local total variation regularization[END_REF]; [START_REF] Peyré | A review of adaptive image representations[END_REF]; [START_REF] Chierchia | Epigraphical splitting for solving constrained convex optimization problems with proximal tools[END_REF]; [START_REF] Couprie | Dual constrained TV-based regularization on graphs[END_REF]. Another approach to overcome these limitations is to replace the gradient operator with a framebased transformation which may yield a more suitable sparse representation of the signals of interest [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The connections between these two different approaches have been studied in [START_REF] Cai | Image restoration: Total variation, wavelet frames and beyond[END_REF].

• Introducing a Total Variation regularization has an effect on smoothing homogeneous regions in the image while preserving edges. Initially introduced in image processing by [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], the TV regularity measure has been proven to be very useful in image recovery and denoising problems [START_REF] Bibliography | Image restoration subject to a total variation constraint[END_REF]; Aujol [2009] and variational stereo methods [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF]. Let ∇ (1) and ∇ (2) denote discrete horizontal and vertical gradients (for example, obtained by cyclic convolutions [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]; [START_REF] Pustelnik | Parallel ProXimal algorithm for image restoration using hybrid regularization[END_REF]). Then, a discrete version of the Total Variation is the following one:

(∀x ∈ R N ) TV(x) = s∈D | ∇ (1) x(s)| 2 + | ∇ (2) x(s)| 2 (2.10)
where D ⊂ Z 2 is the image domain and s ∈ R 2 is the spatial position in the image x, here column-wise reshaped as an N -dimensional vector.

The associated constraint set is

C = {x ∈ R N | TV(x) ≤ τ } (2.11)
where τ can be computed based on some prior information on the target solution (or approximated from an initial estimate).

C = {x ∈ R N | NLTV(x) ≤ τ } (2.
→]-∞, +∞] is dom f = {x ∈ R N |f (x) < +∞}. Γ 0 (R N ) is the class of lower semi- continuous convex functions from R N to ] -∞, +∞] such that dom f = ∅. Let f ∈ Γ 0 (R N ). The conjugate of f is the function f * ∈ Γ 0 (R N ) defined by f * : R N →] -∞, +∞] : u → sup x∈R N x u -f (x) (2.14)
and the subdifferential of f is the set-valued operator

∂f : R N → 2 R N : x → {u ∈ R N |(∀y ∈ R N )(y -x) u + f (x) ≤ f (y)}. (2.15) When f is Gâteaux-differentiable at y ∈ R N , ∂f (y) = {∇f (y)} where ∇f (y) is the gradient of f at y. A differentiable convex function has β-Lipschitz continuous gradient ∇f if (∀(x, y) ∈ R N × R N ) ∇f (x) -∇f (y) ≤ β x -y , (2.16)
where β ∈]0, +∞[. Let C be a nonempty subset of R N . The indicator function of C is

ι C : R N →] -∞, +∞] : x → 0, if x ∈ C +∞, if x / ∈ C.
(2.17)

Let S be a convex subset of R N . The relative interior of S, i.e. the set of points x ∈ S such that the cone generated by -x + S is a vector subspace of R N , is denoted by ri S.

The epigraph of f ∈ Γ 0 (R N ) is the nonempty closed convex subset of R N × R defined as epi f = {(x, ζ) ∈ R N × R | f (x) ≤ ζ}, (2.18)
and the lower level set of f at height ζ ∈ R is the nonempty closed convex subset of R N defined as

lev ≤ζ f = {x ∈ R N | f (x) ≤ ζ}. (2.19) Let f ∈ Γ 0 (R N ) and g ∈ Γ 0 (R N ).
We denote by f g the infimal convolution of f and g

f g : R N → [-∞, +∞] : x → inf y∈R N f (y) + g(x -y).
(2.20)

Problem formulation and algorithms

By generalizing (2.6), a large range of inverse problems can be formulated under the following form:

minimize x∈R N h(x) + R r=1 g r (T r x) s. t.        H 1 x ∈ C 1 , . . . H S x ∈ C S , (2.21)
where, for every r ∈ {1, . . . , R}, T r ∈ R Kr×N and g r ∈ Γ 0 (R Kr ), h : R N → (-∞, +∞] is a convex differentiable function with a µ-Lipschitzian gradient for some µ ∈ (0, +∞) and, for every s ∈ {1, . . . , S}, H s ∈ R Ms×N and C s ⊂ R Ms is a nonempty closed convex subset.

Gradient descent

(∀y ∈ R N ) prox ϕ (y) = argmin u∈R N 1 2 u -y 2 + ϕ(u).
(2.22)

The proximity operator can be interpreted as a sort of subgradient step for the function ϕ, as p = prox ϕ (y) is uniquely defined through the inclusion y -p ∈ ∂ϕ(p) [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. This explains why proximity operators, yielding first order methods, are useful in optimization. As a matter of fact, the proximity operator is the cornerstone of so-called proximal methods.

Proximity operators enjoy many interesting properties (2.24)

Forward-backward splitting aims at finding a point x ∈ R N satisfying the condition.

1 0 ∈ ∂g(x) + ∇h(x) ⇔ x = prox g x -∇h(x) . (2.25)
Indeed, under the assumption that lim x →+∞ g(x) + h(x) = +∞, it can be shown [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] that the sequence (x [n] ) n∈N defined as

x [n+1] = x [n] + λ n prox γng x [n] -γ n ∇h(x [n] ) -x [n]
(2.26) converges to a solution to Problem (2.24) for γ n ∈ ]0, 2/µ[ and λ n ∈ ]0, 1].

1 0 ∈ ∂g(x) + ∇h(x) ⇔ x -∇h(x) -x ∈ ∂g(x) ⇔ x = prox g x -∇h(x)
Douglas-Rachford splitting Let h = 0, R = 2, T 1 = T 2 = Id and S = 0 in (2.21). Then the problem reads minimize

x∈R N g 1 (x) + g 2 (x).
(2.27) Douglas-Rachford splitting aims at finding some points x, ū ∈ R N that satisfy the condition.

2 0 ∈ ∂g 1 (x)+∂g 2 (x) ⇔ ū -x ∈ ∂g 2 (x) x -ū ∈ ∂g 1 (x) ⇔ x = prox g 2 (ū) 0 = prox g 1 (2x -ū) -x.
(2.28) As a matter of fact, under the assumption (ri dom g 1 ) ∩ (ri dom g 2 ) = ∅, it can be shown Combettes and [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] that the sequence (x [n] ) n∈N defined as

x [n] = prox γg 2 (u [n] ),

u [n+1] = u [n] + λ n prox γg 1 (2 x [n] -u [n] ) -x [n] ,
(2.29) converges to a solution to Problem (2.27) for λ n ∈ ]0, 2[ and γ > 0.

PPXA+ Let h = 0 in Problem (2.21). Then, we have to minimize

x∈R N R r=1 g r (T r x) + S s=1
ι Cs (H s x).

(2.30)

Let L = [T 1 . . . T R H 1 . . . H S ] ∈ R Q×N with Q = R r=1 K r + S s=1 M s , ∀v = [v 1 . . . v R+S ] ∈ R Q g(v) = R r=1 g r (v r ) + S s=1
ι Cs (v R+s ), (2.31)

and D = v ∈ R Q v = Lx, x ∈ R N .
Then, Problem (2.30) can be cast over the Douglas-Rachford framework by setting g 1 = ι D and g 2 = g,

minimize v∈R Q ι D (v) + g(v). (2.32)
Moreover, it can be shown Combettes and [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] that, 3(2.29)

⇔            v [n] = prox γg (u [n] ) v [n] = P D (v [n] ) u [n+1] = u [n] + λ n 2 v [n] -u [n] -v [n] u [n+1] = u [n] + λ n ( v [n] -u [n] ), (2.33)
and the sequence ( u [n] ) n∈N converges to a solution of (2.32). Furthermore, the constraint D enforces ( u [n] ) n∈N = (L x [n] ) n∈N , for some sequences (x [n] ) n∈N . When L L is invertible, it can be proved Pesquet and Pustelnik [2012] that (x [n] ) n∈N converges to a solution of (2.30) and that P D (v [n] ) = L x [n] with x [n] = (L L) -1 L v [n] . As a result, if there exists an x ∈ R N such that Lx ∈ ri dom g, the sequence (x [n] ) n∈N generated by

           v [n] = prox γg (u [n] ) x [n] = (L L) -1 L v [n] u [n+1] = u [n] + λ n L(2 x [n] -x [n] ) -v [n] x [n+1] = x [n] + λ n ( x [n] -x [n] ),
(2. 34) converges to a solution to Problem (2.30). The detailed iterations are illustrated in Algorithm 1. Such an algorithm can be considered as a generalization of numerous methods [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems[END_REF]; [START_REF] Eckstein | Parallel alternating direction multiplier decomposition of convex programs[END_REF]

R+s ∈ R Ms x [0] ∈ R N
for n = 0, 1, . . .

(∀r ∈ {1, . . . , R}) v [n] r = prox γgr (u [n] r ) (∀s ∈ {1, . . . , S}) v [n] R+s = PC s (u [n] R+s ) x [n] = R r=1 T r Tr + S s=1 H s Hs -1 R r=1 T r v [n] r + S s=1 H s v [n] R+s (∀r ∈ {1, . . . , R}) u [n+1] r = u [n] r + λn Tr(2 x [n] -x [n] ) -v [n] r (∀s ∈ {1, . . . , S}) u [n+1] R+s = u [n] R+s + λn Hs(2 x [n] -x [n] ) -v [n] R+s x [n+1] = x [n] + λn ( x [n] -x [n] )

Primal-dual algorithms

The scope of primal methods is somewhat limited, as they require to determine the inverse of some linear operators (see Algorithm 1), which turns out to be burdensome or even impossible to compute in many practical cases. Primal-dual algorithms overcome the inversion issue by resorting to the Fenchel-Rockafellar duality framework. In Problem (2.21), let us add a

function f ∈ Γ 0 (R N ), let L = [T 1 . . . T R H 1 . . . H S ] ∈ R Q×N and define g ∈ Γ 0 (R Q ) as in (2.31).
Instead of seeking a solution to the primal problem minimize

x∈R N f (x) + g(Lx) + h(x), (2.35)
or a solution to the dual problem minimize

v∈R Q g * (v) + (f * h * )(-L v), (2.36)
one aims at finding a saddle-point solution to the primal-dual problem

min x∈R N max v∈R Q f (x) -g * (v) F 1 (x,v) + v Lx + h(x) F 2 (x,v)
.

(2.37)

The classical Karush-Kuhn-Tucker theory asserts that a solution to the primal-dual problem is a pair (x, v) satisfying the inclusion,

4 0 ∈ ∂f (x) + L v + ∇h(x) 0 ∈ ∂g * (v) -Lx ⇔ x = prox f x -L v -∇h(x) v = prox g * (v + Lx), (2.38
) where x and v are solutions, respectively, to the primal and dual problems.

The primal-dual inclusion in (2.38) falls in the forward-backward framework. Hence, it can be solved by resorting to the classical forward-backward splitting method or to the forward-backward-forward splitting method [START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF]. The detailed iterations associated with a primal-dual algorithm based on Tseng's splitting are illustrated in Algorithm 2. § 2.

Projections

As mentioned before, computing the projection P C onto a nonempty closed convex subset C of a real Hilbert space R N requires to solve a constrained quadratic minimization problem:

(∀y ∈ R N ) P C (y) = argmin u∈C u -y . (2.39)
The distance to C of every point y ∈ R N is then given by d C (y) = y-P C (y) . However, it turns out that a closed form expression of the solution to (2.39) is available in a limited number of instances. Such well-known examples are the projections onto a closed half-space and a hyperslab: Chapter 2. Convex optimization background Algorithm 2 M+LFBF Combettes and Pesquet [2012] for Problem (2.21).

β = µ + R r=1 Tr 2 + S s=1 Hs 2 1/2 (γn) n∈N ⊂]0, β -1 [ (∀r ∈ {1, . . . , R}) u [0] r ∈ R Kr (∀s ∈ {1, . . . , S}) v [0] s ∈ R Ms x [0] ∈ R N for n = 0, 1, . . . p [n] = ∇h(x [n] ) + R r=1 T r u [n] r + S s=1 H s v [n] s x [n] = prox γnf x [n] -γnp [n]
(∀r ∈ {1, . . . , R}) u

[n] r = prox γng * r u [n]
r + γnTrx [n] (∀r ∈ {1, . . . , R}) u

[n+1] r = u [n] r + γnTr( x [n] -x [n] ) (∀s ∈ {1, . . . , S}) v [n] s = prox γnι * Cs v [n] s + γnHsx [n] (∀s ∈ {1, . . . , S}) v [n+1] s = v [n] s + γnHs( x [n] -x [n] ) p [n] = ∇h( x [n] ) + R r=1 T r u [n] r + S s=1 H s v [n] s x [n+1] = x [n] -γn p [n] -p [n] Example 2.3.1 Let u ∈ R N , let η ∈ R, and let C = {x ∈ R N | x | u ≤ η}.
(2.40)

Then

• If u = 0 and η ≥ 0, in which case C = R N and P C = Id .

• If u = 0 and η < 0, in which case C = ∅.

• If u = 0, in which case C = ∅ and

(∀x ∈ R N ) P C (x) =    x, if x | u ≤ η, x + η -x | u u 2 u, otherwise.
(2.41)

Example 2.3.2 Let u ∈ R N , let (η 1 , η 2 ) ∈ R 2 ,

and let

C = {x ∈ R N | η 1 ≤ x | u ≤ η 2 }. (2.42)
Then exactly one of the following holds:

• u = 0 and η 1 ≤ 0 ≤ η 2 , in which case C = R N and P C = Id .

• u = 0 and η 1 > 0 (resp. η 2 < 0), in which case C = ∅.

• u = 0 and η 1 > η 2 , in which case C = ∅.

• u = 0 and η 1 ≤ η 2 , in which case C = ∅ and

(∀x ∈ R N ) P C (x) =            x + η 1 -x | u u 2 u, if x | u < η 1 , x, if η 1 ≤ x | u ≤ η 2 , x + η 2 -x | u u 2 u, otherwise. (2.43)
The projection onto the 2 -norm ball has also a closed-form expression:

Example 2.3.3 Let τ ∈ R, and set C = {x ∈ R N | x 2 ≤ τ }. Then (∀x ∈ R N ), P C (x) =    x, if x 2 ≤ τ, x √ τ x otherwise.
(2.44)

Subgradient projection When an expression of the direct projection is not available, a possible solution is to approximate the convex set C by a half-space, which leads to the concept of subgradient projection. Let us assume that C = w ∈ R N ϕ(w) ≤ η is nonempty, where η ∈ R and ϕ ∈ Γ 0 (R N ). Then, for any y ∈ R N and t ∈ ∂ϕ(y), the convex set C is a subset of the half-space

C y η = u ∈ R N u -y | t ≤ η -ϕ(y) .
(2.45)

If t = 0, the subgradient projection of y ∈ R N onto C reads P C y η (y) =    y, if ϕ(y) ≤ η, y + η -ϕ(y) t 2 t, otherwise.
(2.46)

Subgradient projection was used in several works. It is involved in Polyak's algorithm by alternating the exact projection and the subgradient projection [START_REF] Polyak | Minimization of unsmooth functionals[END_REF] (see also [START_REF] Shor | Minimization Methods for Non-differentiable Functions[END_REF] for alternative projected subgradient approaches). In [START_REF] Combettes | A block-iterative surrogate constraint splitting method for quadratic signal recovery[END_REF], a parallel block iterative algorithm was proposed in order to solve Problem (2.21). However, the main limitation of this method is that the objective function to be minimized must be quadratic and strictly convex. For recent works about subgradient projection methods, the reader may refer to [START_REF] Slavakis | The adaptive projected subgradient method over the fixed point set of strongly attracting nonexpansive mappings[END_REF]; [START_REF] Bouboulis | Adaptive learning in complex reproducing kernel Hilbert spaces employing Wirtinger's subgradients[END_REF].

Chapter 2. Convex optimization background

Epigraphical projection An alternative approach to circumvent the computation of a direct projection consists in resorting to a set of epigraphical projections [START_REF] Chierchia | Epigraphical splitting for solving constrained convex optimization problems with proximal tools[END_REF]. This approach applies to constraints that can be expressed as the lower level set of a sum of convex functions evaluated over different blocks of the signal:

(∀y = [y 1 , . . . , y L ] ∈ R M ) ϕ(y) = L =1 ϕ (y ) ≤ η, (2.47)
where, for every ∈ {1, . . . , L},

y ∈ R M and ϕ ∈ Γ 0 (R M ) is such that ri(dom ϕ ) = ∅.
The epigraphical-splitting technique consists in introducing an auxiliary vector ζ = ζ 1≤ ≤L ∈ R L into the minimization problem, so that Condition (2.47) can be equivalently rewritten as

       L =1 ζ ≤ η, (∀ ∈ {1, . . . , L}) ϕ (y ) ≤ ζ .
(2.48)

(2.49)

Dealing with additional constraints in the original minimization is not a problem for proximal splitting algorithms, as far as the projections onto the associated constraint sets can be computed. In the present case, the projection onto (2.48) is given by Example 2.3.1 with u = (1, . . . , 1) ∈ R L , whereas the projection onto (2.49) is given by

(∀(y, ζ) ∈ R M × R L ) (p, θ) = P epi ϕ 1 ו••×epi ϕ L (y, ζ), (2.50) where θ = (θ ) 1≤ ≤L , p = [p 1 , . . . , p L ] ∈ R M is blockwise decomposed as y, and (∀ ∈ {1, . . . , L}) (p , θ ) = P epi ϕ (y , ζ ). (2.51)
Therefore, the problem reduces to the lower-dimensional problem of the determination of the projection onto the convex subset epi ϕ of R M × R for each ∈ {1, . . . , L}. The key point of this approach is that epigraphical projections are closely related to proximity operators. Indeed, the epigraphical projector in (2.51) is given by Chierchia et al. [2013]:

p = prox 1 2 (max{ϕ -ζ ,0}) 2 (y ), θ = max{ϕ (p ), ζ }.
(2.52)

In [START_REF] Chierchia | Epigraphical splitting for solving constrained convex optimization problems with proximal tools[END_REF], the closed-form expression of the epigraphical projection has been given for several functions ϕ , such as the absolute power raised to a power β ∈ [1, +∞[, the distance to a convex set and the p -norm with p ∈ {1, 2, +∞}. Note also that, when the epigraphical projection has no explicit form, a polyhedral approximation can still be applied.

Numerical methods Beside the aforementioned approaches, specific numerical methods have been developed to directly solve the constrained quadratic minimization problem in (2.39) when

C = {x ∈ R N | ϕ(x) ≤ η}, (2.53)
with η > 0 and ϕ ∈ Γ 0 (R N ). Indeed, it can be shown that [Bauschke and Combettes, 2011, Prop. 28.30]

P C (x) = x, if ϕ(x) ≤ η, proxλ ϕ (x), otherwise, (2.54)
where λ is a solution in ]0, +∞[ to the equation:

ϕ(prox λϕ (x)) = η. (2.55)
In other words, the projection onto the lower level set of ϕ can be computed by using the proximity operator of ϕ, provided that the correct λ is found. In general, the solutions to (2.55) admit no closed-form expressions. However, one may find λ by standard numerical methods, as the left-hand side of (2.

Conclusion

In this chapter, we have presented a quick overview of modern convex optimization. Firstly, we have introduced inverse problems and discussed the advantages of formulating them as convex optimization problems. Then, we have emphasized the need to look for prior information which aims at incorporating some additional knowledge in the problem formulation. To maintain some smoothness, while preserving some discontinuities of the signal, we have used a regularization term (such as TV and NLTV) which can be used either as a convex constraint set or a penalization function.

Moreover, we have presented some algorithms that can efficiently deal with a large panel of convex problems. A special attention has been given to proximal algorithms, as they are suitable for large-size problems and are amenable to parallel implementations. Finally, to compute the projection onto some convex sets, we have detailed the existing approaches for the exact computation of projections and we have made an overview of the numerical methods in the case where the projection operator admits no closed-form expression.

In the next chapters, we will investigate some examples of constrained convex optimization problems in disparity estimation and multi-view reconstruction.

-Chapter 3 -A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation

"Two Eyes = Three" 3D Dimensions
This chapter deals with the problem of matching a pair of multi-component images by jointly estimating the disparity and the illumination variation.

Since the global formulation is nonconvex, the problem is addressed by solving a sequence of convex relaxations. Each convex relaxation is non trivial and involves many constraints aiming at imposing some regularity on the solution. Experiments demonstrate that the method is efficient and provides better results compared with other approaches. § 3.1 Introduction

The purpose of computer vision is to duplicate the abilities of human vision in order to perceive the observable world through digital images. Among the large panel of human abilities, depth perception is certainly one of the most appealing features to be electronically reproduced, as it may serve as a basis for more complex tasks, such as object recognition, robot navigation, reconstruction of three-dimensional shapes, and realistic scene visualization.

The depth of an object can be perceived thanks to binocular vision. As a simple experiment, hold your finger vertically in front of your eyes and close each eye alternately: you will notice that the finger jumps left and right relative to the background of the scene. Indeed, because of the horizontal separation, the eyes view an object from a slightly different angle, causing a horizontal displacement from the left eye to the right eye (Figure 3.1). The same mechanism underlies computer stereo vision, where the eyes are Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation replaced by two cameras viewing the scene from slightly different angles (Figure 3.2). The disparity is the difference between the positions of the same object in the two images. If an object is far away, the disparity will be small; otherwise, the disparity will be large if the object is close.

Stereo matching is the process of finding the pixels in two (or more) images that correspond to the same 3D point of the observed scene, in order to compute the disparity between such pixels (i.e. the difference in coordinates) and, hence, to infer the depth of the associated 3D point.

Several preliminary steps are required before that stereo matching can take place: image acquisition, calibration, and rectification. In the present work, we focus on stereo matching, but in order to provide a comprehensive introduction to stereo vision, we briefly explain these preliminary steps in the following.

Acquisition and calibration

The pin-hole model is one of the simplest tools to describe image acquisition. In this model, the image is formed by a projection onto a 2D plane:

P : R 3 -→ R 2 (3.1) (X, Y, Z) -→ (f X Z , f Y Z ) (3.2)
where f is the focal distance and (X, Y, Z) are the 3D coordinates. Calibration follows the image acquisition and consists of estimating the internal and the external geometry of the acquisition system. In particular, internal calibration aims at finding the focal point, the optical center of the camera and the size/angle of pixels, while external calibration aims at finding the rotations and translations of the camera with respect to an external reference.

The calibration problem is usually considered as a solved problem, but current research still focuses on this subject, because the calibration is a crucial step in 3D reconstruction.

Rectification

Let us consider a 3D point P that is visible from both cameras. Given its projection P L on the left image, its corresponding projection P R on the right image plane has to be located on the epipolar line. The epipolar line is the intersection between the epipolar plane (defined by the optical centers of the stereo cameras C L and C R and the 3D point P ) and the image plane 3.3. Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation Stereo matching uses triangulation based on epipolar geometry to find the correspondences between left and right images. Indeed, the epipolar line corresponding to a pixel in one image is used to constrain the search for the corresponding pixel in the other image. Rectification is a transformation process that puts the original images onto a common image plane, where the epipolar lines are parallel to the horizontal image axes. So doing, the corresponding search simplifies to a one-dimensional one. Note that the epipolar rectification is trivial when the calibration parameters are known; otherwise, it can be easily achieved by exploiting some corresponding points between the images.

Stereoscopic vision

In the following, we consider that the mapping is done from the left image (called "reference") to the right one. Each pixel in the left image is associated with a pixel in the right one which corresponds to the same point in the scene: they are homologous pixels. The difference between the pixel coordinates of the homologue points is called disparity.

The matching problem is then equivalent to searching for the disparity field u, which attributes a disparity for each pixel in the left image (x L , y L ) as follows:

u : R 2 → R 2 (x L , y L ) → (x L -x R , y L -y R ), (3.3)
where (x R , y R ) are the coordinates in the right image of the point corresponding to (x L , y L ). If the stereo images are rectified, the vertical components of the corresponding pixels are equal (y L = y R ). In this case, the disparity is a scalar and reads

u : R 2 → R (x L , y L ) → x L -x R . (3.4)
The 3D position (X, Y, Z) of a point P can be reconstructed from the perspective projection of P on the image planes of the cameras (Figure 3.4), once the relative position and orientation of the two cameras is known. Let us consider the optical setting in the figure, which is also called standard model.

• C L and C R are two pinhole cameras with parallel optical axes. Let f be the focal length of both cameras.

• The baseline is perpendicular to the optical axes. Let B be the distance between the two lens centers.

• the optical axes lie in the XZ plane, which is parallel to the image plane of both cameras, X axis equals the baseline, and the origin of (X, Y, Z) world reference system is the lens center of the left camera.

x L f = X Z and x R f = X -B Z . (3.5)
One can express the disparity u(x L , y L ) = x L -x R of a 3D point as a function of the depth Z between the corresponding 2D points as follows:

u(x L , y L ) = x L -x R = f B Z (3.6)
Using the relationship (3.6), we can draw two important conclusions:

• disparity is inversely proportional to the depth of a point Z, i.e. far points have low disparity (for example the horizon has disparity of zero) and close points have a high disparity;

Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation

• the disparity is proportional to the baseline B: the smaller the baseline, the lower the disparity. Stereo matching is a challenging task, because of the intrinsic difficulty in finding reliable correspondences between images in the presence of several sources of uncertainty, such as noise, illumination variations, and occlusions (which are pixels only visible from one view of the stereo images). A large number of algorithms have been proposed to perform stereo matching. Since the publication of the Middlebury database, which makes it possible to quantitatively compare the different algorithms, many solutions have been proposed to improve different aspects of stereo matching algorithms.

Related work

In recent years, much progress has been made in dense stereo matching, thanks to the development of In a majority of dense disparity estimation techniques, the scene is assumed to be Lambertian. However, in the presence of illumination variation (that often occurs in practice), this assumption is violated. Some algorithms have been proposed to compensate illumination changes by normalizing pixel values [START_REF] Buchsbaum | A spatial processor model for object colour perception[END_REF]; [START_REF] Finlayson | Comprehensive colour image normalization[END_REF]; van de Weijer et al. [2007]. These algorithm can be applied to compensate for illumination changes as a preprocessing step before the stereo matching. The grey-world assumption algorithm [START_REF] Buchsbaum | A spatial processor model for object colour perception[END_REF] normalizes intensities based on the averages. The comprehensive normalization algorithm [START_REF] Finlayson | Comprehensive colour image normalization[END_REF] attempts to remove the dependency on scene geometry as well as the effect of global illumination iteratively. The grey-edge algorithm van de Weijer et al. [2007] employs the average edge difference for the normalization.

(i) powerful
Other algorithms have been proposed to decrease the illumination variation effect by using invariant similarity measures [START_REF] Hirschmuller | Evaluation of cost functions for stereo matching[END_REF]; Heo et al. [2011] in the matching cost computation. The normalized cross correlation (NCC) can be used as a similarity measure, which minimizes the effect of illumination conditions by compensating the differences of the gain and the bias between stereo images Hirschmuller and Scharstein [2007]. The adaptive normalized cross correlation (ANCC) Heo et al. [2011] employs the NCC with adaptively support weights after log-chrominance transformation and a suitable normalization that revert the illumination changes to a linear model. In addition, methods based on nonparametric local transforms followed by NCC and rank-transform Zabih and Woodfill [1994]; Jung and Kim [2012] have been used.

Recently, the light field has been established as additional information to describe the visual appearance of a scene. One of these approaches uses the light detection to employ the structure tensor of an epipolar plane image to obtain a robust local disparity estimation Wanner and Goldluecke [2013]. In Heo et al. [2013], the authors proposed a stereo color histogram equalization to produce color-consistent stereo images from which they can get more accurate disparity maps. from multicomponent images under illumination variation On the other hand, several techniques have been proposed to model the illumination changes. In this context, Cox et al. [START_REF] Cox | Dynamic histogram warping of image pairs for constant image brightness[END_REF] proposed a dynamic histogram warping that consists of directly matching histogram values and performing a global optimization via dynamic programming. However, this approach is restricted to the case of a spatially invariant relationship between the intensities of the two images. In [START_REF] Davis | BRDF invariant stereo using light transport constancy[END_REF], a new method for stereo reconstruction called Light Transport Constancy was introduced to formulate a rank constraint and enable correspondence of non Lambertian surfaces.

The work in Gennert [1988] shows that the intensities of two corresponding points are related by a spatially varying multiplicative term. Following this model, a stereo algorithm based on a nonconvex cost function was developed. The minimization of the cost function from the associated Euler-Lagrange equations becomes a difficult task. In the same context, based on this multiplicative model, dense disparity map estimation is formulated in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] as a constrained optimization problem in which a strictly convex quadratic objective function is minimized under various convex constraints. The resulting optimization problem is solved via a parallel block iterative algorithm involving subgradient projections.

Chapter outline

In this chapter, we present a convex optimization approach based on a parallel proximal algorithms The remainder of the chapter is organized as follows. The notation, background, and considered model are presented in Section 3.2. Then, in Section 3.3, we describe the adopted parallel proximal algorithm that allows us to solve the derived convex minimization problem. Some simulation results are shown in Sections 3.4 and 3.5 on the MiddleBury data set and other real images. Finally, conclusions are drawn in Section 3.6. § 3.2 Stereo matching model

We consider two multicomponent images of the same scene acquired by a stereoscopic camera. The left (resp. right) view is given by the function

I 1 : R 2 → R K s → I 1 (s) = (I (1) 1 (s), I (2) 1 (s), ..., I (K) 1 (s)) (3.7) (resp. I 2 : R 2 → R K s → I 2 (s) = (I (1) 2 (s), I (2)
2 (s), ..., I

(K) 2 (s)) ).
For every k ∈ {1, . . . , K}, I

(k) 1 (s) (resp. I (k)
2 (s)) represents the k-th component of image I 1 (resp. I 2 ) at position s ∈ R 2 . Particular cases of this framework include color images by taking K = 3. Each component then corresponds to one of the color channels of a specific color system (RGB, YUV, ...). Throughout this chapter, it will be assumed that I 1 and I 2 are differentiable functions.

Corresponding points in the stereo images tend to have similar values. The pixel at position (i 1 , i 2 ) ∈ R 2 in the left image I 1 corresponds to a pixel at position (i 1 , i 2 ) ∈ R 2 in the right image I 2 , the disparity between these pixels being equal to (i 1 -i 1 , i 2 -i 2 ). As already mentioned, when stereo images are rectified [START_REF] Fusiello | A compact algorithm for rectification of stereo pairs[END_REF], the vertical component of the disparity vector vanishes. The disparity thus reduces to u(i 1 , i 2 ) = i 1 -i 1 . Finding for each pixel (i 1 , i 2 ) in the left image the corresponding pixel (i 1 -u(i 1 , i 2 ), i 2 ) in the right image constitutes the goal of stereo matching. The problem is then equivalent to finding a disparity field u : R 2 → [0, +∞) that minimizes a similarity measure. However, under varying illumination conditions, the corresponding points in a stereo pair do not have the same component values. The illumination variation model we will employ is grounded on the work by [START_REF] Gennert | Brightness-based stereo matching[END_REF] and takes the following form:

(∀(i 1 , i 2 ) ∈ R 2 ) I 2 (i 1 -u(i 1 , i 2 ), i 2 ) v(i 1 , i 2 )I 1 (i 1 , i 2 ) (3.8)
where v : R 2 → [0, +∞) represents the illumination field. The problem now is not only to estimate the disparity u but also to estimate the illumination field v. The next section aims at better formulating the considered problem.

Problem formulation

In what follows, we propose a variational approach to jointly estimate the disparity u and the illumination variation v. Based on Model (3.8), we can Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation formulate the stereo matching problem as the minimization of a similarity measure J:

J(u, v) = K k=1 s=(i 1 ,i 2 )∈A\O φ (k) v(s)I (k) 1 (s) -I (k) 2 (i 1 -u(s), i 2 ) (3.9)
where, for every k ∈ {1, . . . , K}, φ (k) is assumed to belong to Γ 0 (R), the class of proper lower semi-continuous convex functions from R to ] -∞, +∞]. In addition, A ⊂ Z 2 is the considered finite image domain, u = (u(s)) s∈A , and v = (v(s)) s∈A . The occlusion areas denoted by O correspond to pixels only visible from one view of the stereo pair, which should not be taken into account in the computation of the similarity measure. Consequently, they have been discarded in the expression of the cost function.

In this chapter, we will be mainly concerned with convex optimization approaches. Unfortunately, J is nonconvex with respect to the variable u.

To tackle this difficulty, we perform the first-order Taylor expansion of the disparity compensated right image around an initial value of the disparity ū. The latter can be derived from a rough estimation, e.g. by a block-based correlation method Tsai et al. [2003]. In practice, we can iteratively update our initial value in order to make our final solution weakly dependent on it. When the magnitude of the difference of the fields u and ū = (ū(s)) s∈A is small enough, we obtain the following expression: for every k ∈ {1, . . . , K} and s = (i 1 , i 2 ) ∈ A,

I (k) 2 (i 1 -u(s), i 2 ) I (k) 2 (i 1 -ū(s), i 2 ) -(u(s) -ū(s))∇ (1) I (k) 2 (i 1 -ū(s), i 2 ) (3.10)
where

∇ (1) I (k) 2
denotes the horizontal gradient of the k-th component of the right image. As a consequence of (3.9) and (3.10), we can approximate the cost function by

J(u, v) = K k=1 s∈A\O φ (k) T (k) 1 (s)u(s) + T (k) 2 (s)v(s) -r (k) (s) (3.11)
where, for every k ∈ {1, . . . , K} and s

= (i 1 , i 2 ) ∈ A,      T (k) 1 (s) = ∇ (1) I (k) 2 (i 1 -ū(s), i 2 ) T (k) 2 (s) = I (k) 1 (s) r (k) (s) = I (k) 2 (i 1 -ū(s), i 2 ) + ū(s)T (k) 1 (s).
(3.12) As previously mentioned, our objective is to jointly estimate u and v. Thus, by defining w = (u, v),

(∀s ∈ A) w(s) = u(s) v(s) (3.13)
and, for every k ∈ {1, . . . , K} and s

∈ R 2 , T (k) (s) = [T (k) 1 (s), T (k)
2 (s)], the above expression can be reexpressed more concisely as

J(w) = K k=1 J (k) (w) (3.14)
where, for every k ∈ {1, . . . , K},

J (k) (w) = s∈A\O φ (k) T (k) (s)w(s) -r (k) (s) . (3.15)
Despite the convexity of the function J, optimizing this criterion is an ill-posed problem, since we have two variables to estimate at each point s and the components of T (k) (s) may locally vanish for some k and s. Therefore, we need to incorporate additional prior information on the desired disparity and illumination variation fields.

Introducing prior information

Our objective here is to introduce prior information on our target solution. This can be done either by adding some regularization terms to J or by incorporating some convex constraints to the problem. In this chapter, we will follow the second approach. Let (S i ) 1≤i≤m denote the m nonempty closed convex sets modelling the constraints one wants to impose. The resulting optimization problem can then be formulated as:

Find w ∈ S = m i=1 S i such that J(w) = inf J(S). (3.16)
Constraint sets defined on the Hilbert space H = R |A| × R |A| can be described as lower level sets:

(∀i ∈ {1, ..., m}) S i = {w ∈ H | f i (w) ≤ δ i } , (3.17)
where, for every i ∈ {1, . . . , m}, f i : H → R is a convex function and δ i ∈ R. One of the potential advantages of this approach, with respect to the regularization formulation, is that it appears often easier, from a physical viewpoint, to set the constraint bounds than to determine the regularization parameters.

We will now review some of the constraints that can be applied to our problem. Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation

Constraints on the disparity

Range values As the minimum and maximum values u min and u max of the disparity field are often known, we propose to introduce a constraint related to these bounds. The corresponding constraint can be expressed as

S 1,1 = {(u(s)) s∈A ∈ R |A| | (∀s ∈ A) u min ≤ u(s) ≤ u max }. (3.18)
It can be noticed that the disparity is always nonnegative, so that u min ≥ 0.

First-order smoothness constraint A piecewise constant behaviour of the disparity field is often expected. Therefore, it is meaningful to introduce a constraint enforcing smoothness in homogeneous areas while preserving edges. One can think of using the

(∀u ∈ R |A| ) TV(u) = s∈A | ∇ (1) u(s)| 2 + | ∇ (2) u(s)| 2 . (3.19)
The associated constraint set is

S 1,2 = u ∈ R |A| TV(u) ≤ τ 2 (3.20)
with τ 2 > 0.

Frame analysis constraint Alternatively, we can adopt a frame analysis approach to construct a smoothness constraint. Such transform can be described by an analysis frame operator F : R |A| → R Q with Q ≥ |A|, which associates to u ∈ R |A| its frame coefficients (F u) q 1≤q≤Q . In particular, it is well-known [START_REF] Meyer | Ondelettes et opérateurs: Ondelettes. Actualités mathématiques[END_REF]

S 1,2 = u ∈ R |A| Q q=1 η q |(F u) q | ≤ τ 2 (3.21)
where (η q ) 1≤q≤Q ∈ [0, +∞[ Q and τ 2 > 0.

Second-order constraint Some recent works [START_REF] Lefkimmiatis | Hessian-based norm regularization for image restoration with biomedical applications[END_REF] have shown that considering the second-order derivatives of the target image can be helpful in data recovery problems. Let a discrete version of the Hessian of u ∈ R |A| at s ∈ A be

∇ 2 u(s) = ∇ (1,1) u(s) ∇ (1,2) u(s) ∇ (1,2) u(s) ∇ (2,2) u(s) . (3.22)
Consequently, a second-order variant of the discrete total variation can be defined as

(∀u ∈ R |A| ) TV 2 (u) = s∈A | ∇ (1,1) u(s)| 2 + | ∇ (2,2) u(s)| 2 + 2| ∇ (1,2) u(s)| 2 (3.23)
and the associated constraint set is

S 1,3 = u ∈ R |A| TV 2 (u) ≤ τ 3 (3.24)
with τ 3 > 0.

Constraints on the illumination field

Range values Similarly to the disparity field, estimations of the minimum and maximum values v min ≥ 0 and v max of the illumination field are often available. The corresponding constraint can be expressed as

S 2,1 = { v(s) s∈A ∈ R |A| | (∀s ∈ A) v min ≤ v(s) ≤ v max }. (3.25)
First-order smoothness constraint Since the illumination field is usually smoothly varying, we can add a quadratic constraint on the discrete gradient of the illumination field. This amounts to considering the following constraint set:

S 2,2 = v ∈ R |A| ∇v 2 2 ≤ κ 2 (3.26)
where κ 2 > 0 and

(∀v ∈ R |A| ) ∇v 2 = s∈A | ∇ (1) v(s)| 2 + | ∇ (2) v(s)| 2 1/2 . ( 3

.27) from multicomponent images under illumination variation

Second-order constraint We can also introduce a constraint on the discrete Hessian of the illumination field, thus yielding

S 2,3 = v ∈ R |A| ∇ 2 v 2 2 ≤ κ 3 (3.28)
where κ 3 > 0 and, for every v ∈ R |A| ,

∇ 2 v 2 = s∈A | ∇ (1,1) v(s)| 2 + | ∇ (2,2) v(s)| 2 + 2| ∇ (1,2) v(s)| 2 1/2 .
(3.29)

Resulting constraints applied to vector w

Based on the constraints presented above, four different constraint sets are considered in this work:

(i) The first one concerns the range of values taken by the disparity and the illumination fields. It is defined as

S 1 = S 1,1 × S 2,1 . (3.30)
(ii) The second set is based on the first-order regularity of the illumination and disparity fields:

S 2 = S 1,2 × S 2,2 (3.31)
(iii) The third set involves the frame analysis formulation instead of the total variation for the disparity, and reads:

S 2 = S 1,2 × S 2,2 . (3.32)
(iv) The last constraint set models the second-order smoothness constraints:

S 3 = S 1,3 × S 2,3 .
(3.33) § 3.3 Proximal approaches for convex optimization

Optimization background

The constraints proposed in Section 3.2.2.3 are separable, since each constraint set S i with i ∈ {1, . . . , m} can be expressed as S 1,i ×S 2,i . 1 In addition, 1 Constraint S 2 will be substituted for S2 in some of our experiments.

for every i ∈ {1, . . . , m}, one can express S 1,i (resp. S 2,i ) as

L -1 1,i (C 1,i ) (resp. L -1 2,i (C 2,i )) where C 1,i (resp. C 2,i ) is a nonempty closed convex subset of R N 1,i (resp. R N 2,i ) and L 1,i (resp. L 2,i ) is a matrix in R N 1,i ×|A| (resp. R N 2,i ×|A| ).
Hence, we have chosen L 1,1 and L 2,1 equal to the |A| × |A| identity matrix. By setting L 1,2 and L 2,2 equal to the concatenation of the horizontal and vertical gradient operators (N 1,2 = N 2,2 = 2|A|), C 1,2 and C 2,2 reduce to 2,1 balls Van Den Berg and Friedlander [2008]; [START_REF] Eldar | Robust recovery of signals from a structured union of subspaces[END_REF]. A similar choice can be made for constraint S 3 by considering linear operators related to the second-order derivatives (N 1,3 = N 2,3 = 3|A|). In the considered cases, the projection onto each convex set C 1,i or C 2,i takes a closed-form expression or it can be computed in a finite number of operations Van Den Berg and Friedlander [2008]. The considered optimization problem then takes the following generic form:

minimize L 1,i u ∈ C 1,i , L 2,i v ∈ C 2,i , i ∈ {1, . . . , m} J(w).
(3.34) mireille Numerical solutions to this problem can be provided by parallel proximal splitting methods. These methods basically consist of iterating computations of proximity operators, as we have shown in Chapter 2.

PPXA+ algorithm

As shown by Algorithm 3, we proceed by splitting function J into the K terms corresponding to each image component and by exploiting the separability of the closed convex constraint sets (C 1,i × C 2,i ) 1≤i≤m . The algorithm is initialized by setting the positive weights (ω 1 , ω 2 , . . . , ω m+K ) ∈]0, +∞[ m+K as well as the variables (z 1,i,0 ) 1≤i≤m+K and (z 2,i,0 ) 1≤i≤m+K associated with the convex constraints (C 1,1 ×C 2,1 , . . . , C 1,m ×C 2,m ) and the split data fidelity criterion J → w → (J (1) (w), . . . , J (K) (w)).

The main loop first consists of computing the proximity operators of the data fidelity terms, as well as the direct projections onto the different convex sets modeling prior information. By using basic properties of proximity operators [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], it can be shown that, for every k ∈ {1, . . . , K}, ω k > 0, and (z (3.35) where z 1 = (z 1 (s)) s∈A , z 2 = (z 2 (s)) s∈A , and Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation Algorithm 3 PPXA+

1 , z 2 ) ∈ R |A| × R |A| , prox J (k) ω k (z 1 , z 2 ) = (z 1 (s)+µ (k) (s)T (k) 1 (s), z 2 (s)+µ (k) (s)T (k) 2 (s)) s∈A
• N otation N1,m+1 = . . . = N1,m+K = |A| N2,m+1 = . . . = N2,m+K = |A| For k = 1, 2, . . . , K L 1,m+k = I , L 2,m+k = I • Initialization (ω1, ω2, . . . , ωm+K ) ∈]0, +∞[ m+K (z [0] 1,i ) 1≤i≤m+K ∈ R N 1,1 × R N 1,2 × • • • × R N 1,m+K (z [0] 2,i ) 1≤i≤m+K ∈ R N 2,1 × R N 2,2 × • • • × R N 2,m+K Q1 = m+K i=1 ωi(L1,i) L1,i -1 , Q2 = m+K i=1 ωi(L2,i) L2,i -1 u [0] = Q1 m+K i=1 ωi(L1,i) z [0] 1,i , v [0] = Q2 m+K i=1 ωi(L2,i) z [0] 2,i

• M ain loop

For n = 0, 1, . . .

For i = 1, 2, . . . , m p [n] 1,i = PC 1,i (z [n] 1,i ) and p [n] 2,i = PC 2,i (z [n] 2,i ) For k = 1, 2, . . . , K (p [n] 1,m+k , p [n] 2,m+k ) = prox J (k) ω m+k (z [n] 1,m+k , z [n] 2,m+k ) Averaging c [n] 1 = Q1 m+K i=1 ωi(L1,i) p [n] 1,i , c [n] 2 = Q2 m+K i=1 ωi(L2,i) p [n] 2,i Updates For i = 1, 2, . . . , m z [n+1] 1,i = z [n] 1,i + λ [n] (L1,i(2c [n] 1 -u [n] ) -p [n] 1,i ) z [n+1] 2,i = z [n] 2,i + λ [n] (L2,i(2c [n] 2 -v [n] ) -p [n] 2,i ) For k = 1, 2, . . . , K z [n+1] 1,m+k = z [n] 1,m+k + λ [n] (2c [n] 1 -u [n] -p [n] 1,m+k ) z [n+1] 2,m+k = z [n] 2,m+k + λ [n] (2c [n] 2 -v [n] -p [n] 2,m+k ) un+1 = u [n] + λ [n] (c [n] 1 -u [n] ) vn+1 = v [n] + λ [n] (c [n] 2 -v [n] ) µ (k) (s) =              1 γ (k) (s) prox γ (k) (s)φ (k) ω k (τ (k) (s)) -τ (k) (s) if γ (k) (s) = 0 and s ∈ O 0 otherwise, (3.36)
τ (k) (s) = T (k) 1 (s)z 1 (s) + T (k) 2 (s)z 2 (s) -r (k) (s) (3.37) γ (k) (s) = |T (k) 1 (s)| 2 + |T (k) 2 (s)| 2 . (3.38)
Subsequently, an averaging of the variables resulting from the computation of the projections and proximity operators is performed before the update equations. It must be emphasized that the matrix inversions required in the algorithm can be efficiently performed by making use of discrete Fourier diagonalization techniques [START_REF] Pustelnik | Relaxing tight frame condition in parallel proximal methods for signal restoration[END_REF].

To gain insights in the convergence properties of this algorithm, assume that, for every k ∈ {1, . . . , K}, φ (k) is a finite function. In our problem, the data fidelity term is convex with respect to the variable (u, v) after Taylor expansion and each individual constraint is also convex with respect to either u or v, as mentioned in Sections 3.2.1 and 3.2.2. Then, provided that the sequence of relaxation parameters (λ n ) n∈N involved in Algorithm 3 is such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2, where λ ∈]0, 2[, [Pesquet and Pustelnik, 2012, Prop. 5.2] allows us to guarantee that the sequence (u n , v n ) n∈N generated by the algorithm converges to a solution to problem (3.16). The value of this parameter is set to 1.5 for all iterations. As the relaxation parameter constitutes an internal parameter of the algorithm, its choice does not modify the optimality of the solution. However, its value plays a role in the convergence speed as illustrated in Table 3.1. According to this table, the value of λ n has been chosen so as to maximize the convergence speed. Finally, note that PPXA+ can be similarly applied in order to minimize J over S 1 ∩ S 2 or S 1 ∩ S 3 . Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation § 3.4 Results with local constraints

Error measures

First of all, we would like to define the quality measures we adopt for assessing the performance of our methods. Let z ∈ R |A| denote an estimate of z ∈ R |A| .

The measures below allow us to provide quantitative evaluations of our results:

(i) Mean Absolute Error (MAE) MAE(z, z) = 1 |A \ O| s∈A\O |z(s) -ẑ(s)|. ( 3.39)
(ii) Average percentage of bad pixels

Err T (z, z) = 1 |A \ O| s∈A\O 1 {|z(s)-ẑ(s)|>T } (3.40)
where T is some given positive threshold value chosen here equal to 2 and, for every s ∈ A,

1 {|z(s)-ẑ(s)|>T } = 1 if |z(s) -ẑ(s)| > T 0 otherwise. (3.41)

Initialization

As previously mentioned, initial estimates ū and v are obtained based on a block matching technique. Usually, two criteria are considered in the minimization performed in the block matching procedure: the sum of squared differences and the sum of absolute differences. In this work, we consider another popular measure, namely the Normalized Cross Correlation (NCC) similarity measure [START_REF] Tsai | The evaluation of normalized cross correlations for defect detection[END_REF]. Note that more sophisticated measures could be envisaged such as the one proposed in [START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF]. NCC is defined as follows:

(∀(i 1 , i 2 ) ∈ A)(∀u ∈ N) ψ i 1 ,i 2 (u) = K k=1 ψ (k) i 1 ,i 2 (u) (3.42)
where

ψ (k) i 1 ,i 2 (u) = (j 1 ,j 2 )∈B I (k) 1 (i 1 + j 1 , i 2 + j 2 )I (k) 2 (x -u + j 1 , i 2 + j 2 ) (j 1 ,j 2 )∈B (I (k) 1 (i 1 + j 1 , i 2 + j 2 )) 2 (j 1 ,j 2 )∈B (I (k) 2 (x -u + j 1 , i 2 + j 2 )) 2 (3.43)
and B corresponds to the matching block (here of size 11 × 11) centered at pixel (i 1 , i 2 ) in the reference image. In order to get the corresponding pixel (i 1 -ū(i 1 , i 2 ), i 2 ) in the right image, the matching process determines ū(i 1 , i 2 ) which maximizes the similarity measure

ψ i 1 ,i 2 : ū(i 1 , i 2 ) ∈ argmax u ∈ N∩[u min ,umax] ψ i 1 ,i 2 (u). (3.44) Once ū(i 1 , i 2 ) is found, v(i 1 , i 2
) is evaluated by the following weighted least squares estimation:

v(i 1 , i 2 ) = K k=1 θ (k) (j1,j2)∈B I (k) 1 (i 1 + j 1 , i 2 + j 2 )I (k) 2 (i 1 -ū(i 1 , i 2 ) + j 1 , i 2 + j 2 ) K k=1 θ (k) (j1,j2)∈B I (k) 1 (i 1 + j 1 , i 2 + j 2 ) 2 , (3.45)
where (θ (k) ) 1≤k≤K are nonnegative real weights. In our tests, the illumination variation is the same for all the components and, for YUV images, the best results have been obtained by choosing a weight equal to 1 for the Y component, and 0 for the U and V components. However, with the previous technique, some remaining artefacts can be observed in the disparity maps, e.g. contouring effects and oversmoothing in some areas. In order to reduce these artifacts and limit the influence of occlusion areas, a bidirectional matching process is performed. More precisely, we first apply the NCC method to find the left disparity map denoted by ū1 (i 1 , i 2 ), by proceeding as described above. Similarly, the same method is applied to find the right disparity map denoted by ū2 (i 1 , i 2 ) by taking now the right view as the reference image. Finally, we derive the initial disparity ū(i 1 , i 2 ) as follows:

ū(i 1 , i 2 ) = ū2 (i 1 -ū1 (i 1 , i 2 ), i 2 ). (3.46)
The upper and lower bounds in the constrained formulation are estimated using prior knowledge when available or from the initial estimates ū and v, otherwise.

Algorithm implementation

In this section, we summarize the algorithm steps and parameters that are required in our joint disparity and illumination variation estimation method.

Inputs:

-Occluded areas O estimated from the initial disparity map based on the ordering Yuille and Poggio [1984] and border constraints.

Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation

-(T (k) 1 , T (k) 2 , r (k) ) for k ∈ {1, 2, 3}.
-Minimum and maximum estimates of the fields to be estimated.

-Upper bounds for the constraints.

Parameter choice:

-For S 1 , ω 1 =100, ω 2 =100.

-For S 2 , S 3 or S 4 : ω 3 =200, ω 4 =200.

-For J (1) ,J (2) and J (3) : ω 5 =ω 6 =ω 7 =10.

-λ n ≡ 1.5.

The weights were chosen so that all the terms in the averaging step have a contribution of the same order.

Stopping criterion:

-The algorithm is stopped when the following condition is satisfied for 10 successive iterations: u n+1 -u n < u n where = 10 -5 .

Outputs:

-The final values of the disparity vector u n and the illumination variation vector v n .

Gray level images without illumination variation

Although the purpose of this work is to focus on situations where illumination variation arises, in this first experiment we aim at illustrating the performance of the proposed approach on grayscale images (K = 1) without any variation of the illumination field (for every s ∈ A, v(s) = 1). This allows us to compare our method with some recent state-of-the-art methods referenced in the Middlebury website. 2 In Table 3.2, we display the results obtained on 3 different stereo pairs (Teddy, Venus, and Cones) for which we compare the results obtained by:

• the correlation-based approach used to initialize our approach (NCC),

• a graph-cut based approach (Gc+Occ) [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF],

• a total variation based algorithm (ConvexTV) [START_REF] Bibliography | A convex formulation of continuous multilabel problems[END_REF],

• a method based on the knowledge of Ground Control Points (Global-GCP) Wang and Yang [2011],

• a histogram-based method (HistoAggr) Min et al.

[2011],

• a curvelet-based approach (CurveletSupWgt) Mukherjee et al.

[2010],

• the approach based on set theoretic estimation (subg. proj.) in Miled et al.

[2009a],

• our proposed method based on PPXA+ minimizing an 1 cost subject to constraints S 1,1 , S 1,2 , and S 1,2 .

Note that the frame employed to define constraint S 1,2 simply consists of an overcomplete Haar wavelet representation carried out over 1 resolution level. Other wavelet families have been tested without observing significant improvements. The proposed approach appears to be competitive with respect to the other ones, although most of them are not able to deal with illumination variation. Moreover, as already mentioned, the proposed method can be implemented on a GPU architecture, thus reducing drastically the computation time. Indeed, the proximity operators and the projections onto the different convex sets can be computed in parallel. For example, with a GPU implementation, when considering the range values S 1,1 and the total variation constraints S 1,2 , the following executing times per iteration were measured: for Teddy: 21.7ms, for Venus: 5.72ms, and for Cones: 21.98ms [START_REF] Gaetano | Parallel implementations of a disparity estimation algorithm based on a proximal splitting method[END_REF].

Let us now turn our attention to more challenging scenarios including illumination variation both for grayscale and color images. Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation

Gray level images with illumination variation

In this section, we evaluate our method on grayscale images (K = 1) in the presence of real and simulated illumination changes. The results are obtained on three stereo pairs depicted in Figure 3.6, named "Books", "Dolls", and "Parking meter". In a first experiment, we estimate the disparity map and the illumination variation (real case) on the Books pair downloaded from MiddleBury website. We take here φ (1) = | • | for the cost function and consider the two convex constraints S 1 and S 2 for which the upper bounds were evaluated on the ground truth maps. In Figure 3.7, we compare our results with the subgradient projection approach proposed in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], where a strictly convex quadratic function is used. We also compare our approach with two naive procedures for the disparity field estimation: the correlation based approach ignoring the illumination variation and a histogram equalization technique followed by a variational approach for disparity estimation, in the spirit of the approach we developed.

Concerning the choice of the criterion and of the constraint sets, we present in Figure 3.8 results obtained by using the following settings:

• for a)-b): S 1 , S 2 , and φ (1) = | • |; • for c)-d) same constraints as a)-b) and φ (1) = (•) 4 ;
• for e)-f): S 1 , S 2 , S 3 , and φ (1) = | • |;

• for g)-h) same constraints as e)-f) and φ

(1) = | • | 3 2 .
By comparing the obtained results, it turns out that the 1 cost is a good choice, and that the use of second-order constraints leads to a marginal improvement.

Results with local constraints 51

In a third experiment (see Figure 3.10), we show the results provided by the proposed method (with φ (1) = |•|) for a real image stereo pair from JISCT database,3 taken under simulated illumination variation. To introduce a significant illumination variation, the left image is kept unchanged and the right image is modified by multiplying it with the Gaussian profile g depicted in Figure 3.9 and defined as:

(∀(i 1 , i 2 ) ∈ A) g(i 1 , i 2 ) = α 1 exp - (i 1 -i 1,0 ) 2 + (i 2 -i 2,0 ) 2 2ξ 2 + α 2 (3.47)
Figure 3.9: Synthetic Gaussian profile.

where (i 1,0 , i 2,0 ) is the image center, ξ = 512 is the illumination deviation, α 1 and α 2 are given constants equal to 1.8 and -0.6 respectively. In this case, the constraint bounds in S 1 and S 2 are not known. Hence, they are computed in the first place using the initial disparity, and then divided by a factor 2, thus assuming that the expected result should be smoother than the initial estimate. [2010], the nonlocal total variation regularization is introduced using the second order Taylor expansion of the data term applied in motion estimation. The nonlocal total variation regularization incorporates this information directly into the energy definition yielding:

NLTV(u) = s∈A n∈Ns⊂Ws ω s,n |u(s) -u(n)| 2 , (3.48)
where W s is a set of positions n ∈ A \ {s} located into a Q × Q window centered at s, where Q ∈ N is odd, and N s is the actual neighborhood support of pixel s, which may contain a smaller subset of positions in W s . The associated constraint is expressed as

S 1,4 = u ∈ R |A| NLTV(u) ≤ τ 4 (3.49)
with τ 4 > 0.

Weight estimation for NLTV operator

For every s ∈ A, we design the support N s by selecting the pixels n ∈ W s that are most similar to s according to their weights ω s,n . For every n ∈ W s , the weight ω s,n depends on the similarity between patches built around the pixels s and n of the disparity map, as will be shown afterwards. In general, from multicomponent images under illumination variation such similarities can be estimated by following two different approaches. The first solution consists of computing the weights on the left image, assuming that object boundaries provide cues to localize discontinuities in the disparity map, as will be shown next. The alternative solution consists of a two-step approach, where an estimate of the disparity map is first computed with the TV regularizer and subsequently used to compute the weights. We follow the latter approach, as the former one suffers from the fact that similar patches in the left image may not necessarily correspond to regions with similar disparities.

In both cases, the weights are computed as follows

ω s,n = ω s exp -δ B s F s x -B n F n x 2 , (3.50)
where δ ∈ ]0, +∞[, ω s ∈ (0, +∞), B s (resp. B n ) selects a Q × Q patch centered at position s (resp. n), and F s (resp. F n ) is a linear processing of the image depending on the position s (resp. n). The constant ω s is set so as to normalize the weights (i.e. n∈Ns ω s,n = 1). In the present work, we consider the foveated self-similarity measure recently introduced in Foi and Boracchi [2012], due to its better performance in denoising. This approach can be derived from (3.50) by setting F s (resp. F n ) to a set of low-pass Gaussian filters whose variances increase as the spatial distance from the patch center s (resp. n) grows.

For every s ∈ A, the neighbourhood N s is built according to the procedure described in Gilboa and Osher [2007]. In practice, we limit the size of the neighbourhood |N s | to 14, and set Q = 11 and Q=3.

Results

To emphasize the potential of our approach, we realized some experiments on the already introduced Middlebury benchmark database. For the NLTV constraint, we have used an epigraphical projection in order to handle the projection onto 12 ball norm Chierchia et al. [2013]. Concerning the choice of the constraint sets, for NLTV (resp. TV) we use S 1 , S 2,2 , and S 4 (resp. S 1 and S 2 ). Note that the matrix inversion required by primal algorithms (such as PPXA+) can be performed by resorting to variable splitting as detailed in [START_REF] Peyré | Group sparsity with overlapping partition functions[END_REF]. In this respect, primal-dual methods often constitute a better choice, as they require nothing but the linear operators and their associated adjoint operators. For this reason, in the following we employ the 3.5. Experiments with non local constraints 55 M+LFBF algorithm(described in the previous chapter), which offers a good performance, robustness to numerical errors, and guaranteed convergence. The images used in our experiments are displayed in Figure 3.12. mireille mireille -Case 1: Gray images without illumination variation (Table 3.3). According to the percentage of bad pixels and the Mean Absolute Error, we can see that the NLTV approach performs better than the total variation one, with a gain of 0.83% in terms of bad pixels. -Case 2: Color images under illumination variation (Table 3.4). Also in this case, NLTV performs better than TV with a gain in the percentage of bad pixels of 0.51% and 2.01% for Cloth and Wood pairs respectively.

Chapter 3. A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation Finally, in Figure 3.13, we illustrate the recovered disparity maps using our method for color stereo pairs Cloth and Wood. These visual results clearly show the efficiency of our method in recovering the depth estimation in addition to illumination variations, for multicomponent images. § 3.

Conclusion

In this chapter, we have investigated the application of a parallel proximal algorithm to the dense disparity estimation problem, for multicomponent (e.g. color) stereo matching, under illumination variation conditions. The proposed approach is flexible, as it allows us to consider various convex objective functions and constraints. It is also able to exploit the potentials offered by multicore/GPU parallel system architectures. However, one of the current limitations of the proposed method is that it is limited to convex cost functions, thus requiring a linearization of the original disparity model and making the estimation potentially sensitive to high matching error values. In the next chapter, we extend this approach to the multi-view context, but we resort to the convex relaxation proposed in Cremers et al. [2011], which allows us to consider nonconvex cost functions with no approximation in the computed solution (except for a quantization error). -Chapter 4 -

Multi-view disparity estimation

"When you feel life at crossroads, you need higher perspective view." Toba Beta

In the previous chapter, we discussed the problem of depth estimation from two views, under illumination variation. However, using two views, we have a limited field of view and we can only see a portion of a scene. In particular, a global description of objects often cannot be reconstructed from only two viewpoints due to occlusions. For this purpose, multi-view imaging has been used for many vision tasks to make the objects more visible.

This chapter addresses the problem of multi-view disparity estimation from rectified images. A sequence of disparity maps is estimated by using a multilabel optimization approach. For each of these maps, information from all the available views are exploited. The joint estimation problem is solved in a convex manner without adding much complexity with respect to the stereo matching case. The good performance of the proposed method is illustrated through experiments. § 4.1 Introduction

We consider the classical problem of estimating a dense 3D structure of a scene from a collection of calibrated views. Being one of the fundamental problems in computer vision, this issue has lately gained much attention in this field and remains an active research area. There are different types of approaches to tackle this problem, depending on the exploited image information. All these methods aim at reversing the image formation process. However, an image is essentially formed by projecting the 3D scene onto a 2D plane (which causes the loss of depth) and, in a mathematical sense, the inverse projection mapping does not exist, since all 3D points along a visual 59 60 Chapter 4. Multi-view disparity estimation ray are projected onto the same image 2D point. This makes the problem of multi-view 3D reconstruction ill-posed. Thus, additional constraints are needed in order to make the reconstruction problem well-posed. Multi-view is known to produce detailed reconstructions with a quality close to that of laser scanned models [START_REF] Seitz | A comparison and evaluation of multi-view stereo reconstruction algorithms[END_REF].

Related Work

The problem of 3D modeling from multiple images or multi-view stereopsis has been widely studied in the field of computer vision, computer graphics, and computational photography. A large panel of variational methods for multi-view 3D reconstruction inherits the active contour framework originally proposed for image segmentation [START_REF] Kass | Snakes: Active contour models[END_REF]. In this context, the problem consists of modeling a continuous two-dimensional surface as the minimizer of an appropriate energy functional, by using geodesic active contour models [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Caselles | Minimal surfaces: A three dimensional segmentation approach[END_REF], level sets Faugeras and Keriven [1998], triangle meshes [START_REF] Hernandez | Silhouette and stereo fusion for 3d object modeling[END_REF]; [START_REF] Duan | Shape reconstruction from 3D and 2D data using PDE-based deformable surfaces[END_REF] or graph cuts Vogiatzis et al. [2007]. A generalization of this approach has been developed in [START_REF] Pons | Multi-view stereo reconstruction and scene flow estimation with a global image-based matching score[END_REF], which allows to replace the classical photo-consistency estimation with a global matching score on the entire image domain.

On the other side, the Patch-based Multi-view Stereo algorithm proposed by Furukawa and Ponce [2010] appears as the state-of-the-art method among those based on feature extraction and matching. However, methods based on feature extraction and matching may show significant errors caused by severe occlusions and highly reflective surfaces. Probabilistic models have been proposed as an alternative approach to handle complicated scene geometry. These models do not make hard decisions about surface geometry or appearance: they explicitly represent uncertainties by assigning probabilities to multiple hypotheses. Early works along this line of reasoning [START_REF] Broadhurst | A probabilistic framework for space carving[END_REF]; [START_REF] Bhotika | A probabilistic theory of occupancy and emptiness[END_REF] can be regarded as extensions of Space Carving [START_REF] Kutulakos | A theory of shape by space carving[END_REF] and, more recently, algorithms based on generative models for the reverse image formation process have been introduced Gargallo et al. [2007]; [START_REF] Liu | A complete statistical inverse ray tracing approach to multi-view stereo[END_REF]. Using a Bayesian approach, these algorithms infer the maximum a posteriori estimation from the joint probability of all the images.

A different approach amounts at recovering a set of consistent viewdependent depth maps from an image sequence, instead of reconstructing a complete 3D model. This is mainly motivated by applications such as view synthesis, depth-based segmentation, and video enhancement. Early approaches in depth estimation Okutomi and [START_REF] Okutomi | A multiple-baseline stereo[END_REF]; [START_REF] Collins | A space-sweep approach to true multi-image matching[END_REF] used local window-based methods and employed a local winner-takes-all strategy. Our work is somewhat related to that of [START_REF] Kang | Extracting view-dependent depth maps from a collection of images[END_REF]; Gargallo and Sturm [2005], which also aims at inferring consistent depth maps from multiple images. The authors in [START_REF] Gargallo | Bayesian 3D modeling from images using multiple depth maps[END_REF] formulated 3D modeling from images as a Bayesian MAP problem, and solved it using the expectation-maximization algorithm. They used the estimated depth map to determine the visibility prior. A multiple-depth-map prior was finally used to smooth and merge the depths while preserving discontinuities. [START_REF] Larsen | Temporally consistent reconstruction from multiple video streams using enhanced belief propagation[END_REF] presented an approach for 3D reconstruction from multiple synchronized video streams. In order to improve the final reconstruction quality, they used the optical flow to find corresponding pixels in the subsequent frames acquired by the same camera, and enforced the temporal consistency in reconstructing successive frames. With the observation that the depth error in conventional stereo methods grows quadratically with depth, Gallup et al. [2008] proposed a multi-baseline and multi-resolution stereo method to achieve constant depth accuracy by varying the baseline and resolution proportionally to depth.

In summary, although many approaches have been proposed to model 3D objects or to estimate depths using multiple input images, the problem of how to appropriately extract information and recover consistent depths from an image sequence remains challenging. In the present work, a global multilabel approach is proposed for multi-view disparity estimation, which is grounded on the use of a convex optimization strategy. A global solution to the problem is provided by minimizing a regularized similarity measure. The regularization enforces the smoothness of the sought disparity maps while preserving coherent contours in the set of images. The proposed formulation handles possibly nonconvex similarity measures. We also point out that the additional complexity required by the application of our method to the multi-view context is small with respect to the stereo case.

The chapter is organized as follows. In Section 4.2, we provide a general formulation of the multi-view disparity estimation problem. In Section 4.3, we describe a convex relaxation method based on the multilabel approach and we propose an algorithm to solve the resulting convex optimization problem. In Section 4.4, we describe an alternative convex relaxation method based on an improved first-order-Taylor approximation around an initial estimate. In Section 4.5, we present some numerical comparisons of the proposed multi-view disparity estimation methods and we show the effectiveness of the first one with respect to the standard stereoscopic approach (i.e. based on two views) and other multi-view estimation approach. Some conclusions are given in Section 4.6. § 4.2 Multi-view disparity estimation Although two images are enough to estimate the depth of a scene from point correspondences, having additional views presents two significant benefits. Firstly, the estimation process becomes more accurate: since there is more information about the same 3D point, this redundancy can be exploited in order to improve the depth estimation quality. Secondly, occluded pixels (i.e. pixels that have no match in some of the images) can be handled using the other images where a correspondence exists. The uncertainty related to occluded pixels can thus be reduced in multi-view systems.

This work investigates the disparity estimation problem for an arbitrary number of views acquired by rectified cameras [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]; [START_REF] Nozick | Multiple view image rectification[END_REF], with the camera centers aligned.

Problem formulation

We suppose that N rectified views (I n ) 1≤n≤N are taken along a straight line with equal distance between any two consecutive camera positions (Figure 4.1). n,m ) s∈A are proportional to the distance between the two cameras, the focal length f , and to the inverse depths of the objects (Z (s) ) in the scene:

(∀s ∈ A), u (s) n,m = f (m -n)B Z (s) . (4.1)
Using the same notation as in the previous chapter, s = (i 1 , i 2 ) ∈ A denotes a pixel coordinate and A ⊂ R 2 is the image support. Under classical epipolar geometry assumptions (see Chapter 3), the intensities of corresponding pixels in the rectified images are related by

I (i 1 ,i 2 ) n = I (i 1 -u (s) n,m ,i 2 ) m . (4.2)
The overall number of disparity fields is equal to N (N -1)/2, but these are obviously related since they correspond to the same 3D scene. In particular, for every n ∈ {1, . . . , N }, if k n = n is some given index in {1, . . . , N } (typically, k n = n + 1 is the right image index, when n < N and k n = N -1 is the left image index, when n = N ), we have

(∀(n, m) ∈ {1, . . . , N } 2 , n = m) u n,m = α n,m u n,kn (4.3)
where α n,m is a real-valued constant depending on the relative positions of the cameras. According to Equation (4.1), when the cameras are equidistant, α n,m = m-n kn-n . Hence, by using (4.3), the number of disparity fields to be estimated reduces to N . By taking advantage of the N views, the u n,kn field can be estimated with a variational approach that minimizes a criterion of the form

f n (u n,kn ) = N m=1 m =n s∈Dn,m ψ(I (s) n , I (i 1 -αn,mu (s) n,kn ,i 2 ) m ) (4.4)
where D n,m denotes some finite subset of A of unoccluded pixels (i.e. pixels in the n-th view that can be matched with pixels in the m-th view) and ψ : R 2 → R is a similarity measure. As explained in the following sections, the proposed approach is applicable to any kind of function of interest ( 2 norm, truncated 1 cost function, etc...). § 4.3 Approach based on multilabeling

Convex representation

In order to estimate u n,kn , we propose to adopt a multilabel approach [START_REF] Bibliography | A convex formulation of continuous multilabel problems[END_REF]; [START_REF] Hiltunen | Comparison of two proximal splitting algorithms for solving multilabel disparity estimation problems[END_REF]. More specifically, the disparity u n,kn is quantized over Q + 1 quantization levels r 0 , r 1 , . . . , r Q , such that r 0 < r 1 < • • • < r Q . We employ a bijection that maps the disparity map 64 Chapter 4. Multi-view disparity estimation u n,kn ∈ {r 0 , . . . , r Q } P , with P = card(A), to a multivariate binary field θ n = (θ n,1 , . . . , θ n,Q ) ∈ B with components

(∀ q ∈ {1, . . . , Q})(∀ s ∈ A) θ (s) n,q = 1 if u (s)
n,kn ≥ r q 0 otherwise where

B = {θ n ∈ ({0, 1} P ) Q | (∀ s ∈ A) 1 ≥ θ (s) n,1 ≥ • • • ≥ θ (s) n,Q ≥ 0}.
More precisely, there exists a bijection between the set of disparity fields u n,kn ∈ {r 0 , . . . , r Q } P and the set of multivariate label images θ n,q as defined above. We have then

(∀ s ∈ A) u (s) n,kn = r 0 + Q q=1 (r q -r q-1 )θ (s) n,q (4.5)
The estimation of the quantized disparity u n,kn is thus equivalent to finding the associated multivariate binary field θ n ∈ B. It can be further observed that the global error measure can be re-expressed as a linear function of θ n as follows:

f n (θ n ) = N m=1 m =n s∈Dn,m Q q=0 ψ(I (s) n , I (i 1 -αn,mrq,i 2 ) m )(θ (s) n,q -θ (s) n,q+1 ) (4.6)
where θ n,0 = 1 and θ n,Q+1 = 0. Consequently, the minimization problem can be expressed as:

minimize θn∈B f n (θ n ) + ρ(θ n ) (4.7)
where ρ is some regularization function used to promote the spatial regularity of the binary images. A typical choice for ρ is the isotropic discrete total variation which is denoted by tv. Then, our problem becomes: (r q -r q-1 ) tv(θ n,q ) (4.9)

minimize θn∈B f n (θ n ) + µ Q q=1 (r q -r q-1 ) tv(θ n,q ), µ > 0. ( 4

Approach based on multilabeling 65

where

R = {θ n ∈ ([0, 1] P ) Q |(∀ s ∈ A) 1 ≥ θ (s) n,1 ≥ • • • ≥ θ (s) n,Q ≥ 0}.
Although the convex function f n (θ n ) + ρ(θ n ) is minimized over a convex set R, the uniqueness of a solution is not guaranteed [START_REF] Ekeland | Convex analysis and variational problems[END_REF]. In order to solve Problem (4.9), we need to introduce some auxiliary functions.

Firstly, let g n : R → R be a linear function defined as

(∀ θ n ∈ R) g n (θ n ) = f n (θ n ) - N m=1 m =n s∈Dn,m ψ I (s) n , I (i 1 -αn,mr 0 ,i 2 ) m (4.10) = Q q=1 s∈A θ n,q ς n,q = θ n ς n , (4.11)
where • • is the standard scalar product in the space (R P ) Q and ς n = (ς n,1 , . . . , ς n,Q ) ∈ (R P ) Q . For every q ∈ {1, . . . , Q}, the components of vector ς n,q are given by

(∀ s ∈ A) ς (s) n,q = N m=1 m =n 1 n,m (s) ψ(I (s) n , I (i 1 -αn,mrq,i 2 ) m ) -ψ(I (s) n , I (i 1 -αn,mr q-1 ,i 2 ) m )
where, for every m ∈ {1, . . . , N } with n = m, 1 n,m (s) = 1 if s ∈ D n,m and 0 otherwise. In addition, the constraint set R can be decomposed as the intersection of two more tractable sets:

θ n ∈ R ⇔ (θ n ∈ E 1 and Lθ n ∈ E 2 ) (4.12)
where

E 1 = ([0, 1] P ) Q , E 2 = ([0, +∞[ P ) Q-1 , and L : (R P ) Q → (R P ) Q-1
is a linear operator, calculating the successive differences between the Q components of θ n , defined as

(∀ θ n ∈ (R P ) Q ) Lθ n = (ζ n,1 , . . . , ζ n,Q-1 ), (4.13)
and for all s ∈ A,   

ζ n,1 (s) . . .

ζ n,Q-1 (s)    = L    θ n,1 (s) . . . θ n,Q (s)    (4.14)
where L ∈ R (Q-1)×Q , is such that

L =          1 -1 0 . . . . . . 0 0
         . (4.15)
Finally, the total variation function can be expressed as

ρ(θ n ) = µ h(Dθ n ) (4.16)
where

D : (R P ) Q → (R 2P
) Q is a block diagonal operator whose q-th diagonal term with q ∈ {1, . . . , Q} corresponds to a discrete gradient operator applied to θ n,q computed in the horizontal and vertical directions. D is a linear operator given by

D =              D 1 0 . . . 0 D 2 0 . . . 0 D 1 0 D 2 . . . . . . . . . 0 0 0 0 D 1 0 0 . . . 0 D 2              (4.17)
where D 1 and D 2 are the spatial gradient operators operating in the horizontal/vertical directions with periodic extension at boundaries. These operators correspond to 2D filters with frequency responses 1 -exp(-2πν 1 ) and 1 -exp(-2πν 2 ), where (ν 1 , ν 2 ) are the 2D horizontal/vertical frequency variables ( 2 = -1). The above function h is a sum of 1,2 norms:

1 ∀ δ n = (δ n,1 , . . . , δ n,Q ) ∈ (R 2P ) Q h(δ n ) = Q q=1 (r q -r q-1 ) δ n,q 1,2 .
With these definitions, the optimization problem reads minimize

θn∈(R P ) Q g n (θ n ) + µ h(Dθ n ) s. t. θ n ∈ E 1 , Lθ n ∈ E 2 , (4.18)
and the corresponding constrained formulation reads minimize

θn∈(R P ) Q g n (θ n ) s. t.      θ n ∈ E 1 , Lθ n ∈ E 2 , h(Dθ n ) ≤ η, (4.19)
where η ≥ 0.

1 See Section 6.1 for the definition of the 1,2-norm.

Proximal solution

Proximity operators

Since g n is a linear function, it can be further deduced from the standard properties of the proximity operator [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] that the proximity operator of γg n + ι E 1 with γ ∈ 0, +∞] is given by:

prox γgn+ι E 1 = P E 1 (• -γς n ) (4.20)
For the proximity operator of h, needed to solve Problem (4.18), we have:

∀ δ n = (δ n,1 , . . . , δ n,Q ) ∈ (R 2P ) Q prox µh δ n = ξ n (4.21)
where, ξ n = (ξ n,1 , . . . , ξ n,Q ) ∈ (R 2P ) Q , for every q ∈ {1, . . . , Q} and s ∈ A, ξ (s) n,q = max 0, 1 -µ(r q -r q-1 ) Remark 4.3.1 The projection onto the polyhedron R may be computed in a finite number of iterations by solving the dual problem, which consists of a quadratic criterion subject to a positivity constraint. However, since at each iteration we need to compute a large number of projections (for each pixel in the image), this approach turns out to be slower than splitting the constraint R into

δ (s) n,q 1,2 δ (s) n,q . ( 4
E 1 + E 2 .

M+LFBF algorithm

We propose here to use a primal-dual proximal algorithm whose main advantage is the absence of any matrix inversion. More precisely, we use the M+LFBF algorithm proposed in [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallelsum type monotone operators[END_REF]. The algorithm is initialized by setting the constants β and γ, and the variables v The main loop consists of computing the proximity operator given in (4.20), as well as the direct projections onto the different convex sets modeling prior information.

Algorithm 4 M+LFBF Combettes and Pesquet [2012]

β = D 2 + L 2 1/2 ∈]0, (1 + β) -1 [, γ ∈ [ , (1 -)/β] v [0] 1 ∈ (R 2P ) Q , v [0] 2 ∈ (R P ) Q-1 θ [0] n ∈ (R P ) Q for t = 0, 1, . . . p [t] = D v [t] 1 + L v [t] 2 θ [t] n = PE 1 θ [t] n -γςn -γp [t] v [t] 1 = v [t] 1 + γDθ [t] n -γP lev ≤η h v [t] 1 +γDθ [t] n γ v [t] 2 = v [t] 2 + γLθ [t] n -γPE 2 v [t] 2 +γLθ [t] n γ v [t+1] 1 = v [t] 1 + γD( θ [t] n -θ [t] n ) v [t+1] 2 = v [t] 2 + γL( θ [t] n -θ [t] n ) p [t] = D v [t] 1 + L v [t] 2 θ [t+1] n = θ [t] n -γ p [t] -p [t] Remark 4.3.2 mireille 1-If γ ∈ [ , (1-) β ] where β = D 2 + L 2 1/2 and ∈]0, (1 + β) -1 [, then the sequence (θ [t]
n ) t∈N generated by the Algorithm 4 converges to a solution θ n to Problem 4.19 [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallelsum type monotone operators[END_REF]).

2-In practice, to accelerate the convergence, the step-size γ has to be chosen as large as possible (for example, = 0.001/(β + 1) and γ = (1 -)/β).

3-Due to the forms of operators L and D, we have

L = 2 (4. 23
)
(that is the maximum singular value of matrix L). In addition, according 4.4. Improved Taylor approximation 69 to (4.17), we have

D = D * D 1 2 = D * 1 D 1 + D * 2 D 2 1 2 = max ν 1 ,ν 2 (|1 -exp(2πν 1 )| 2 + |1 -exp(2πν 2 )| 2 ) 1 2 = 2 √ 2. (4.24)
This shows that β = 2 √ 3. § 4.4 Improved Taylor approximation

An alternative way to deal with multi-view disparity estimation consists of relaxing the data fidelity term by resorting to a linearization around an initial estimate. In Chapter 3, we presented an approach based on the first-order Taylor approximation. In the present section, we improve this approach by using a second-order Taylor approximation [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF], which allows us to tighten the relaxation of the data fidelity term in the absence of illumination variation.

Recall that (I n ) 1≤n≤N denote the N rectified views, and the intensities of corresponding pixels between two views are related by

I (i 1 ,i 2 ) n = I (i 1 -u (s) n,m ,i 2 ) m . (4.25)
When the magnitude of the difference between the fields u n,m and the initial estimate ūn,m is small enough, we obtain the following expression: ∀(n, m) ∈ {1, . . . , N } 2 , n = m and ∀s = (i 1 , i 2 ) ∈ D n,m

I (i 1 -u (s) n,m ,i 2 ) m I (i 1 -ū (s) n,m ,i 2 ) m -(u (s) n,m -ū(s) n,m )∇ (1) I (i 1 -ū (s) n,m ,i 2 ) m + 1 2 (u (s) n,m -ū(s) n,m ) 2 ∇ (1,1) I (i 1 -ū (s) n,m ,i 2 ) m + o(|u (s) n,m -ū(s) n,m | 2 ) O(|u (s) n,m -ū (s) n,m | 2 ) , (4.26)
where ∇ (1) (resp. ∇ (1,1) ) denotes the first-order (resp. the second order) derivative along the first space dimension (see Chapter 3, Equation 3.22). By deriving the two members of Equation (4.25), we obtain:

∇ (1) I (s) n = ∇ (1) (I (i 1 -u (s) n,m ,i 2 ) m ) = ∇ (1) I (i 1 -u (s) n,m ,i 2 ) m (1 -∇ (1) u (s) n,m ). (4.27) Since u (s)
n,m is locally smooth, it appears reasonable to assume that ∇ (1) u (s) n,m 0. The results in Section 4.5 confirm the validity of this assumption. By performing a second term Taylor expansion around the initial estimate ū, we Chapter 4. Multi-view disparity estimation obtain:

∇ (1) I (i 1 -u (s) n,m ,i 2 ) m ∇ (1) I (i 1 -ū (s) n,m ,i 2 ) m -(u (s) n,m -ū(s) n,m )∇ (1,1) I (i 1 -ū (s) n,m ,i 2 ) m + o(|u (s) n,m -ū(s) n,m | 2 ). (4.28)
Then, combining Equation (4.26) and (4.28), we obtain

I (i 1 -u (s) n,m ,i 2 ) m I (i 1 -ū (s) n,m ,i 2 ) m - 1 2 (u (s) n,m -ū(s) n,m )∇ (1) (I (i 1 -ū (s) n,m ,i 2 ) m + I (s) n ) + o(|u (s) n,m -ū(s) n,m | 2 ) (4.29)
Following the approach in Chapter 3, the criterion is reformulated as: To minimize the criterion J n , we need to incorporate additional prior information on the disparity field (see Chapter 3, Section 3.2.2.1), which yields a minimization problem expressed as: minimize

J n (u n,kn ) = N m=1 m =n s∈Dn,m ψ I (s) n -I (i 1 -αn,mu (s) n,kn ,i 2 ) m = N m=1 m =n s∈Dn,m ψ T (s) n,m u (s) n,kn -r (s) n,m (4
u n,kn ∈R P J n (u n,kn ) s. t. u n,kn ∈ S 1,1 , u n,kn ∈ S 1,2 , (4.32)
where

S 1,1 = {u n,kn ∈ R P | (∀s ∈ R P ) u min ≤ u (s) n,kn ≤ u max }, S 1,2 = {u n,kn ∈ R P | TV(u n,kn ) ≤ τ 2 }.
To solve Problem 4.32, we will use Algorithm 5, which is a modified version of PPXA+ employed in Chapter 3. 2 Algorithm 5 PPXA+

• Initialization For m = 1, . . . , N -1 Lm+2 = I (ω1, ω2, . . . , ωN+1) ∈]0, +∞[ N +1 (z [0] i ) 1≤i≤N +1 ∈ R P × R 2P × R P • • • × R P Q1 = N +1 i=1 ωi(Li) Li -1 u [0] n,kn = Q1 N +1 i=1 ωi(Li) zi,0
For l = 0, 1, . . .

p [l] 1 = PC 1,1 (z [l] 1 ), p [l] 2 = PC 1,2 (z [l] 2 ) For m = 1, . . . , N -1 p [l] m+2 = prox J (m) n ω m+2 (z [l] m+2 ) c [l] = Q1 N +1 i=1 ωi(Li) p [l] i z [l+1] 1 = z [l] 1 + λ [l] (L1(2c [l] -u [l] n,kn ) -p [l] 1 ) z [l+1] 2 = z [l] 2 + λ [l] (L2(2c [l] -u [l] n,kn ) -p [l] 2 ) For m = 1, . . . , N -1 z [l+1] m+2 = z [l] m+2 + λ [l] (2c [l] -u [l] n,kn -p [l] m+2 ) u [l+1] n,kn = u [l] n,kn + λ [l] (c [l] -u [l] n,kn ) § 4.

Numerical results

In this section, we present some numerical examples to illustrate the performance of the proposed method for multi-view disparity estimation. In this regard, we will show the improvement resulting from the use of more than two images for estimating the disparity map and the flexibility of our algorithm to handle both convex and nonconvex similarity measures ψ. The quality of the results is evaluated by computing the mean absolute error (MAE) and the percentage of unoccluded badly estimated pixels between the computed disparity u n,kn and the ground truth u n,kn . A pixel is considered as unoccluded in I n if it appears at least in one image of the considered image sequence (I m ) 1≤m≤N,m =n . Thus, the percentage of false pixels reads

Err T (u n,kn , u n,kn ) = 1 N m=1 m =n D n,m s∈ N m=1 m =n Dn,m 1 {|u (s) n,kn -û (s) n,kn |>T } , (4.33)
where T is a predefined threshold, as in Chapter 3. Since ground truth disparity maps for image sequences are hard to find, a set of experiments are conducted to evaluate the proposed framework.

Disparity map estimation from multiple images

Firstly, we limit ourselves to the dataset available at the Middlebury stereo vision website. Each dataset consists of 9 or 7 rectified views taken from equidistant cameras along a line. Ground truth disparity maps, created by using the structured lighting technique of [START_REF] Scharstein | High-accuracy stereo depth maps using structured light[END_REF], are only provided for two viewpoints in each dataset. In the case when we have 9 views, we estimate the disparity between view 2 and view 4, which is then evaluated using the provided ground truth divided by two. Finally, by adding more images, we demonstrate the efficiency of our method in recovering depth information from multiple images. Our analysis is conducted by estimating the disparity map u 2,4 on three image sequences (Sawtooth, Teddy, and Cones -Figure 4.2) and by comparing different similarity measures in the optimization criterion.

Multilabel approach

The main advantage of this method is the possibility to handle both convex and nonconvex error similarity measures ψ, namely the 4.1: The (MAE, Err 1 ) results on Sawtooth image, with Q=14, r0=2, rQ=9 and (∀ q ∈ {2, . . . , Q}) rq -rq-1 = 0.5, and η = 6.9 × 10 3 . Moreover, we illustrate in Figure 4.3 the estimated disparity maps. The red regions indicate the pixels associated with a wrong disparity value. 4.4, the improved first-order expansion reduces the error introduced by the Taylor approximation around the same initial estimate. Note that the improved approximation can only be applied to the disparity estimation problem for gray level images without illumination variation. This is the reason why it was not employed in Chapter 3, where the considered images were affected by illumination variations.

1 1, 1 2 1 2 , (I2, I4) (0.

Performance comparison of different methods

In this section, we compare the approaches based on multi-labeling and Taylor approximation with the method in Woodford et al. [2009] based on graph-cuts.

To do so, we employ some standard datasets taken from the Middlebury Database ("Sawtooth", "Teddy" and "Cones", Figure 4.2) and a synthetic dataset called "Sponza", which consists of five images with the corresponding ground-truth disparity maps (Figure 4.4) generated by using the photorealistic rendering Embree software. 4 Each disparity map is computed by exploiting depth information and image characteristics (Chapter 3, Equation 3.6). Table 4.5 collects the percentage of error pixels obtained with the improved-first-order Taylor approximation, the multilabel approach, and the method in Woodford et al. [2009]. One can see that the multilabel approach performs better than the others, as it employs a tighter convex relaxation of the data fidelity term. This can be also observed in Figure 4.6, where the disparity map of "Sponza" dataset estimated with the aforementioned methods is displayed. Nonetheless, the Taylor-approximation approach is much faster (see Chapter 3, Section 3.4), i.e. it trades off the precision with the speed.

Implementation parameters To perform the above experiments, we used the following setups for the considered algorithms:

• Multilabel. For Sawtooth, Teddy, and Cones, we fix the same parameters as in Tables 4.1, 4.2, and 4.3. For Sponza, sequence, we have: ψ = 1, ( =40), Q=80, r 0 = 11, r Q = 91, (∀ q ∈ {2, . . . , Q}), r q -r q-1 = 1, and η = 54 × 10 3 .

• Improved Taylor approximation For Sponza we have ψ = 1 , τ 2 = 11 × 10 4 , and [u min , u max ]= [START_REF] Hajlaoui | Satellite image restoration in the context of a spatially varying point spread function[END_REF]91]. (For Sawtooth, Teddy, and Cones, see Table 4.4).

• Graph cuts [START_REF] Woodford | Global stereo reconstruction under second-order smoothness priors[END_REF] The authors incorporated the second-order smoothness priors for a global optimization of the stereo reconstruction. They employed triple cliques to estimate depth by using graph cuts. We have used the code publicly available 5 for gray level images. The default setting has been unchanged (except for Sponza).

Disparity map sequence estimation

In this section, we show the reconstruction result of a 3D scene from a sequence of disparity maps. Figure 4.7 displays such a disparity sequence, where five images were used to estimate each map. The improvement is obvious with respect to the stereo case, using the 1 similarity measure. Finally, an extended disparity map is obtained from the disparity map sequence (see Figure 4.8) as well as a 3D reconstruction of the "Sponza" sequence (see Figure 4.9) performed from the estimated maps. Five views, of the 3D scene, are presented from different camera position: Left, Down, Center, Up and Right. § 4.

Conclusion

In this chapter, we have proposed a multilabel optimization approach for performing multi-view disparity estimation. Rather than estimating a single depth map, we associate a depth map with each input image (or a subset of those images). We have proposed two different approaches to estimate the depth map of each chosen reference frame, which differ in the type of convex relaxation employed. The first approach employs a relaxation based (MAE, Err 2 ) obtained using five images for each disparity map estimation, Q=80, r 0 =11, r Q =91 and (∀ q ∈ {2, . . . , Q}), r q -r q-1 = 1.

For comparison, we indicate the values of (MAE, Err 2 ) stereo obtained in the stereo case when u n,kn is estimated from (I n ,I kn ). ψ is the 1 -norm.

on multi-labeling, with the possibility of using nonconvex similarity measures. So doing, there is no more dependence of the method performance on the initial estimate, at the cost of an increased complexity. The second approach employs a relaxation based on an improved first-order Taylor approximation around an initial estimate. The numerical issues related to the computation of our disparity sequence have been addressed by using a proximal algorithm. The proposed method leads to improved results compared with other methods. It is also worth emphasizing that this method can be efficiently implemented on multicore architectures. -Chapter 5 -

Proximity operators of discrete information divergences

"Knowledge is of two kinds. We know a subject ourselves, or we know where we can find information upon it. "

Samuel Johnson

The notion of information measure plays a central role in the fields of probability, statistics, and information theory, as it allows one to assess how close two distributions are from each other. Among the large panel of available measures, a special attention has been paid to ϕ-divergences, e.g. Kullback-Leibler, Jeffreys-Kullback, Hellinger, Chi-Square, I α , and Renyi divergences. Indeed, these divergences have been studied extensively in many different contexts, including convex optimization. In the latter context, the solution of optimization problems may be challenging.

One of the main limitations of existing convex optimization methods is that they restrict the use of ϕ-divergences to problems in which the minimization is performed w.r.t. one of the divergence arguments. A possible solution to overcome this issue is to resort to alternating minimization approaches, which however require specific assumptions to be valid [START_REF] Tseng | Alternating projection-proximal methods for convex programming and variational inequalities[END_REF]; [START_REF] Bauschke | Joint minimization with alternating Bregman proximity operators[END_REF]. In the present chapter, we derive closed-form expressions of the proximity operators of the aforementioned ϕ-divergences. These novel results enrich the list of functions for which the proximity operator can be easily computed. Building upon these results, we are able to address a wider range of convex optimization problems by means of proximal methods.

84

Chapter 5. Proximity operators of discrete information divergences § 5.1 introduction Claude [START_REF] Shannon | A mathematical theory of communication[END_REF], in his seminal paper "A Mathematical Theory of Communication", introduced a powerful mathematical framework to quantify the notion of information, laying the foundations for the field of information theory and originating a major revolution in communications and related fields. The power of the framework introduced by Shannon is reflected in his two main results known as the Source Coding Theorem and the Channel Coding Theorem.

In the same spirit as [START_REF] Shannon | A mathematical theory of communication[END_REF], it was natural to investigate functionals which allow one to quantify how much information is shared between two probability distributions. As pointed out by [START_REF] Liese | Convex Statistical Distances[END_REF], the origins of these ideas can be traced back to the early twentieth century literature in the works by [START_REF] Pearson | On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can reasonable be supposed to have arisen from random sampling[END_REF] and [START_REF] Hellinger | Neue begrndung der theorie quadratischer formen von unendlichvielen vernderlichen[END_REF], although research in this area became much more prolific after the publication of Shannon's paper [START_REF] Shannon | A mathematical theory of communication[END_REF]. Motivated by this work, [START_REF] Kullback | On information and sufficiency[END_REF] introduced the information measure now known as the Kullback-Leibler divergence within the context of hypothesis testing.

Yet a higher level of generalization in the area of probabilistic divergences was achieved by the work by [START_REF] Csiszár | Eine informations theoretische ungleichung und ihre anwendung auf den beweis der ergodizität von markoffschen ketten[END_REF] (and independently also [START_REF] Ali | A general class of coefficients of divergence of one distribution from another[END_REF]; [START_REF] Cover | Elements of information theory[END_REF]), who introduced the notion of ϕ-divergences for probability distributions, a framework which encompasses a vast number of information measures currently used in the literature, including the Kullback-Leibler (KL), Jeffreys-Kullback (JK), Hellinger (Hel), Chi square, and I α divergences. The works by [START_REF] Arimoto | Information-theoretical considerations on estimation problems[END_REF]; [START_REF] Lecam | Convergence of estimates under dimensionality restrictions[END_REF]; [START_REF] Burbea | On the convexity of some divergence measures based on entropy functions[END_REF]; [START_REF] Basseville | Distance measures for signal processing and pattern recognition[END_REF]; [START_REF] Vajda | Theory of Statistical Inference and Information[END_REF]; [START_REF] Barron | Distribution estimates consistent in total variation and two types of information divergence[END_REF]; [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF]; [START_REF] Amari | Alpha-divergence is unique, belonging to both f-divergence and bregman divergence classes[END_REF]; Österreicher [2013]; [START_REF] Nielsen | On the chi square and higher-order chi distances for approximating f-divergences[END_REF]; [START_REF] Basseville | Divergence measures for statistical data processing[END_REF] have enriched this formalism and its applications in an extensive manner.

The KL divergence is known to play a prominent role in the computation of channel capacity and rate-distortion functions. One can address these problems with the celebrated alternating minimization algorithm proposed by [START_REF] Blahut | Computation of channel capacity and rate-distortion functions[END_REF]; [START_REF] Arimoto | An algorithm for computing the capacity of arbitrary discrete memoryless channels[END_REF]. However, other approaches based on geometric programming [START_REF] Chiang | Geometric programming duals of channel capacity and rate distortion[END_REF] may provide more efficient numerical solutions. The generalized KL divergence is also known as the I-divergence. It plays an important role in inverse problems for recovering the signal of interest in the presence of Poisson noise. In such a case, the I-divergence is usually employed as a data fidelity term. For instance, an alternating projection technique was proposed in Byrne [1993], where both the data fidelity term and the regularization term are based on the KL divergence. The problem was formulated in a similar manner in [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF], whereas in [START_REF] Fessler | Hybrid Poisson/polynomial objective functions for tomographic image reconstruction from transmission scans[END_REF] [2012] are grounded on proximal splitting methods. These proximal tools were recently shown to offer both efficient and flexible solutions to a wide class of possibly nonsmooth convex minimization problems (see [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]; [START_REF] Parikh | Proximal algorithms[END_REF] and references therein). Note that, in all of the aforementioned works, one of the two variables of the KL divergence is fixed.

As the KL divergence is a Bregman distance, optimization problems involving this function can also be addressed by using the alternating minimization approach proposed in Bauschke et al. [2006]. However, the required optimization steps may be difficult to implement and the convergence of the algorithm is only guaranteed under restrictive conditions. Moreover, a Kullback-Leibler proximal algorithm generalizing the EM algorithm was investigated in [START_REF] Chrétien | On EM algorithms and their proximal generalizations[END_REF]. The KL divergence is then used as a metric for the maximization of a log-likelihood function, rather than being one of the terms of the objective function. Recently, the authors in Vemuri et al. [2011]; [START_REF] Meizhu | Total bregman divergence and its applications to shape retrieval[END_REF] defined a new similarity measure called total Kullback Leibler divergence, which has the potential benefit of being invariant to axis rotations. On the other side, the classical symmetrization of Kullback-Leibler divergence (also known as Jeffreys-Kullback divergence Jeffreys [1946]) has been recently used in the k-means algorithm (as a replacement of the squared difference) [START_REF] Nielsen | Sided and symmetrized bregman centroids[END_REF]; [START_REF] Nielsen | On the symmetrical Kullback-Leibler Jeffreys centroids[END_REF], yielding analytical expressions of the centroids in terms of the Lambert W function.

text The Hellinger divergence was named after Ernst [START_REF] Hellinger | Neue begrndung der theorie quadratischer formen von unendlichvielen vernderlichen[END_REF] and later rediscovered in [START_REF] Beran | Minimum hellinger distance estimates for parametric models[END_REF]; [START_REF] Devijver | Pattern Recognition: A Statistical Approach[END_REF]; [START_REF] Gibbs | On choosing and bounding probability metrics[END_REF]; [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF]; [START_REF] Rauber | Probabilistic distance measures of the dirichlet and beta distributions[END_REF] under different names (such as Jeffreys-Matusita distance). In the field of information theory, the Hellinger divergence is commonly used in the analysis of nonparametric density estimation [START_REF] Lecam | Convergence of estimates under dimensionality restrictions[END_REF]; van de Geer [1993], statistics and Data Analytics Chang-Hwan [2012], and machine learning Cieslak et al. [2012].

text The Chi-square divergence can be traced back to [START_REF] Pearson | On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can reasonable be supposed to have arisen from random sampling[END_REF], where it was used to quantitatively assess whether an observed phenomenon tends to confirm or deny a given hypothesis. This work heavily contributed to the development of modern statistics: in 1984, the journal Science referred to [START_REF] Pearson | On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can reasonable be supposed to have arisen from random sampling[END_REF] as "one of the 20 most important scientific breakthroughs". Moreover, Chi-square divergence was also successfully applied in different contexts, such as information theory and signal processing, as a dissimilarity measure between two probability distributions [START_REF] Cover | Elements of information theory[END_REF]; [START_REF] Park | Estimation of symmetric chisquare divergence for point processes[END_REF].

The I α -divergence was originally proposed in [START_REF] Chernoff | A measure of asymptotic efficiency for tests of a hypothesis based on a sum of observations[END_REF] to statistically evaluate the efficiency of an hypothesis test. Subsequently, it was recognized as an instance of more general divergence classes [START_REF] Ali | A general class of coefficients of divergence of one distribution from another[END_REF], such as the class of ϕ-divergences [START_REF] Csiszár | Information measures: A critical survey[END_REF] and the class of Bregman divergences Amari [2009]. Indeed, I α -divergence has been extended by many researchers as in [START_REF] Liese | Convex Statistical Distances[END_REF]; [START_REF] Zhang | Divergence function, duality, and convex analysis[END_REF]; [START_REF] Minka | Divergence measures and message passing[END_REF]; Amari [2009], it has been considered in the context of Non-negative Matrix Factorization where the hyper-parameter α is associated with characteristics of a learning machine Cichocki et al. [2008].

text Rényi divergence was introduced in [START_REF] Rényi | On measures of entropy and information[END_REF] as a measure of information related to Rényi entropy. According to the definition given in [START_REF] Harremoës | Interpretations of rnyi entropies and divergences[END_REF], Rényi divergence measures "how much a probabilistic mixture of two codes can be compressed". It has been studied and applied in many areas [START_REF] Vajda | Theory of Statistical Inference and Information[END_REF]; Liese andVajda [1987, 2006] including image registration problems [START_REF] Hero | Applications of entropic spanning graphs[END_REF].

To the best of our knowledge, existing approaches for optimizing convex criteria involving these divergences are often restricted to problems where the minimization is performed w.r.t. one of the arguments of the divergence, or they are based on an alternating minimization process which requires specific assumptions to be valid. In the context of proximal methods, this lacuna can be explained by the fact that there exist few results concerning the proximity operator of non-separable convex functions (as opposed to separable ones, for which many results are available Chaux et al. In this work, we develop a novel proximal method that allows us to address more general forms of optimization problems. The expressions we derive for these divergences enrich the list of functions for which such proximity operators can be easily computed. In addition to its flexibility, the proposed approach leads to parallel proximal algorithms that can be efficiently implemented on multicore architectures.

The remaining of the chapter is organized as follows. In Section 5.2, we present the general form of the addressed optimization problem and introduce the notation used in this work. In Section 5.3, we study the proximity operators of ϕ-divergences and some of their properties. The aforementioned divergence proximity operators are described in Section 5.4. In Section 5.5, we explain how the forms of these proximity operators can be employed to deduce the expression of the projection onto the epigraph of some convex functions. Finally, Section 5.6 concludes the chapter. § 5.2 Problem formulation

Optimization problem

The objective of this chapter is to address convex minimization problems involving a discrete information divergence D, which can be formulated as follows:

Problem 5.2.1 Let D be a function in Γ 0 (R P × R P ). Let A and B be matrices in R P ×N , and let u and v be vectors in R P . For every s ∈ {1, . . . , S}, let R s be a function in Γ 0 (R Ks ) and T s ∈ R Ks×N . We want to minimize

x∈R N D(Ax + u, Bx + v) + S s=1 R s (T s x).
(

Hereabove Γ 0 (H) denotes the class of convex functions f defined on a real Hilbert space H and taking their values in ] -∞, +∞] which are lowersemicontinuous and proper (i.e. their domain dom f on which they take finite values is nonempty). Allowing functions to take the +∞ value makes Problem 5.2.1 a versatile one. For example, the involved functions may take finite values only for nonnegative-valued arguments. This also allows us to include in this formulation convex constrained problems by letting R s with s ∈ {1, . . . , S} be equal to the indicator function ι C of some nonempty closed convex set C. Recall that the indicator of a nonempty closed convex subset C of H is defined as

(∀x ∈ H) ι C (x) = 0 if x ∈ C +∞ otherwise. (5.2)
In applications to inverse problems, R s may also be some regularization function serving to enforce the smoothness of the sought solution or to model some additional prior information, e.g. its sparsity after some appropriate linear transform T s .

A special case of interest in information theory is obtained when u = v = 0, N = 2P , and

A = [I P 0], B = [0 I P ] (5.3)
where I P is the P × P identity matrix. By decomposing the vector x as x = [p q ] (5.4)

where (p, q) ∈ (R P ) 2 , and by setting (∀s ∈ {1, . . . , S})

T s = [U s V s ] (5.5)
where U s and V s are matrices in R Ks×P , Problem 5.2.1 becomes:

Problem 5.2.2 Let D be a function in Γ 0 (R P × R P ). For every s ∈ {1, . . . , S}, let R s be a function in Γ 0 (R Ks ), and let U s and V s be matrices in R Ks×P . We want to minimize

p∈R P , q∈R P D(p, q) + S s=1 R s (U s p + V s q). (5.6)
When the variables p and q correspond to probability mass distributions, a standard example of convex constraint subset of R 2P is the product of two unit simplexes of R P : C = (p (1) , . . . , p (P ) , q (1) , . . . , q (P ) ) ∈ [0, 1] 2P P i=1 p (i) = 1 and P i=1 q (i) = 1 .

(5.7)

Considered class of divergences

We will focus on additive information measures of the form ∀p = (p (i) ) 1≤i≤P ∈ R P ∀q = (q (i) ) 1≤i≤P ∈ R P D(p, q) = P i=1 Φ(p (i) , q (i) )

(5.8) where Φ is defined as follows D(p, q) = 0 ⇔ p = q.

∀(υ, ξ) ∈ R 2 Φ(υ, ξ) =                ξ ϕ υ ξ if υ ∈ [0,
(5.12)

Examples of ϕ-divergences will be provided in Sections 5.4.1, 5.4.2, 5.4.3, 5.4.4 and 5.4.6. For a thorough investigation of the rich properties of ϕ-divergences, the reader is refered to [START_REF] Csiszár | Eine informations theoretische ungleichung und ihre anwendung auf den beweis der ergodizität von markoffschen ketten[END_REF]; [START_REF] Ali | A general class of coefficients of divergence of one distribution from another[END_REF]; [START_REF] Basseville | Distance measures for signal processing and pattern recognition[END_REF]. Other divergences (e.g. Rényi divergence) are expressed as

∀(p, q) ∈ R P × R P D g (p, q) = g D(p, q) (5.13)
where g is some increasing function. Then, provided that g ϕ(1) = 0, for every (p, q) ∈ ([0, +∞[ P ) 2 such that [p q ] ∈ C where C is given by (5.7), we have D g (p, q) ≥ 0. Note that, from an optimization standpoint, minimizing D or D g (possibly subject to constraints) makes no difference, so that we will only address problems involving D in the rest of this chapter.

Proximity operators

Proposition 5.2. 3 Moreau [1965]; [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] Let f ∈ Γ 0 (H). Then, (i) for every x ∈ H, prox f x ∈ dom f .

(ii) For every (x, x) ∈ H 2 , x = prox f x ⇔ x -x ∈ ∂f (x), where ∂f (x) denotes the subdifferential of f at x. If f is Gâteaux differentiable at x, ∂f (x) = {∇f (x)} where ∇f (x) denotes the gradient of f at x.

(iii) For every (x, z) ∈ H 2 , prox f (•+z) x = prox f (x + z) -z. (5.14)
(iv) For every (x, z) ∈ H 2 and for every α ∈ R,

prox f +• z+α x = prox f (x -z). (5.15) (v) Let f * be the conjugate function of f (f * ∈ Γ 0 (H) defined by f * : u → sup x∈H ( x | u -f (x))). For every x ∈ H, prox f * x = x -prox f x. (5.16)
(vi) Let G be a real Hilbert space and let T : G → H be a bounded linear operator, the adjoint of which is denoted by T * . If T • T * = κ Id where κ ∈ ]0, +∞[, then

prox f •T = Id + 1 κ T * • (prox κf -Id ) • T.
(5.17)

Proximal splitting algorithms

As soon as we know how to calculate the proximity operators of the functions involved in Problem

∈ ]0, +∞[ S+1 , t0,0 ∈ R P , t1,0 ∈ R P , t2,0 ∈ R K 1 , . . . , tS+1,0 ∈ R K S x0 = (ω0A A + ω0B B + S s=1 ωsT s Ts) -1 (ω0A t0,0 + ω0B t1,0 + S s=1 ωsT s ts+1,0)
for n = 0, 1, . . . In this algorithm, ω 0 , . . . , ω S are weighting factors and (λ n ) n≥0 relaxation factors. Typical values for these parameters will be indicated later. For every n ≥ 0, the variables e 0,n ∈ (R P ) 2 , e 1,n ∈ R K 1 , . . . , e S,n ∈ R K S model possible errors in the computation of the proximity operators. For instance, these errors arise when the proximity operator is not available in a closed form and one needs to compute it through inner iterations. Under some technical conditions, the convergence of PPXA+ is guaranteed. (ii) There exists x ∈ R N such that

(r0,n, r1,n) = prox ω -1 0 D(•+u,•+v) (t0,
Ax + u ∈ ]0, +∞[ P , B x + v ∈ ]0, +∞[ P (∀s ∈ {1, . . . , S}) T s x ∈ ri(dom R s ). (5.25)
(iii) There exists λ ∈]0, 2[ such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2.

(iv) (∀s ∈ {0, . . . , S}) n∈N e s,n < +∞.

If the set of solutions to Problem (5.2.1) is nonempty, then any sequence (x n ) n∈N generated by Algorithm (6) converges to an element of this set.

It can be noticed that, at each iteration n, PPXA+ requires to solve a linear system in order to compute the intermediate variable y n . The computational cost of this operation may be high when N is large. Proximal primal-dual approaches Chen and Teboulle [1994]; [START_REF] Esser | A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science[END_REF]; Chambolle and Pock [2011]; Briceño-Arias and Combettes [2011]; [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallelsum type monotone operators[END_REF]; [START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF]; [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] allow us to circumvent this difficulty. An example of such an approach is the M+LFBF (Monotone+Lipschitz Forward Backward Forward) method [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallelsum type monotone operators[END_REF] which takes the following form: parallel. Second, they are tolerant with respect to errors in the computations of the proximity operators. 4

§ 5.3 Main result

As shown by Proposition 5.2.4, we need to compute the proximity operator of a scaled version of a function Φ ∈ Γ 0 (R 2 ) as defined in (5.9). In the following, Θ denotes a primitive on ]0, +∞[ of the function ζ → ζϕ (ζ -1 ). The following functions will subsequently play an important role:

ϑ -: ]0, +∞[ → R : ζ → ϕ (ζ -1 )
(5.27)

ϑ + : ]0, +∞[ → R : ζ → ϕ(ζ -1 ) -ζ -1 ϕ (ζ -1
).

(5.28)

A first technical result is as follows:

Lemma 5.3.1 Let γ ∈ ]0, +∞[, let (υ, ξ) ∈ R 2
, and define Proof. The derivative of ψ is

χ -= inf ζ ∈ ]0, +∞[ ϑ -(ζ) < γ -1 υ (5.29) χ + = sup ζ ∈ ]0, +∞[ ϑ + (ζ) < γ -1 ξ (5.
(∀ζ ∈ ]0, +∞[) ψ (ζ) = ϕ(ζ -1 ) -(ζ + ζ -1 )ϕ (ζ -1 ) + γ -1 υζ -γ -1 ξ = ζ γ -1 υ -ϑ -(ζ) + ϑ + (ζ) -γ -1 ξ. (5.32)
The function ϑ -is decreasing as the convexity of ϕ yields

(∀ζ ∈ ]0, +∞[) ϑ -(ζ) = -ζ -2 ϕ (ζ -1 ) ≤ 0.
(5.33)

Proof. If (χ -, χ + ) ∈ ]0, +∞[ 2 , we have already shown that Conditions (ii) and (iii) are satisfied if and only χ -< χ + . If χ -∈ ]0, +∞[ and χ + = +∞ (resp. χ -= 0 and χ + ∈ ]0, +∞[), we still have

ψ (χ -) = ϑ + (χ -) -γ -1 ξ < 0.
(5.42)

(resp. ψ (χ + ) = χ + γ -1 υ -ϑ -(χ + ) > 0) (5.43)
which shows that Condition (ii) (resp. Condition (iii)) is always satisfied.

By using the same expressions of χ -and χ + as in the previous lemmas, we obtain the following characterization of the proximity operator of any scaled version of Φ:

ψ ( ζ) = 0 ⇔ ζ γ -1 υ -ϑ -( ζ) = γ -1 ξ -ϑ + ( ζ) (5
∈ ]0, +∞[ is (∀(υ, ξ) ∈ R 2 ) prox γΦ (υ, ξ) = υ + γ ln ζ, ξ + γ( ζ -1 -1) if exp(υ/γ) > 1 -γ -1 ξ (0, 0)
otherwise (5.59) where, in the first case, ζ is the unique minimizer on ] exp(-υ/γ), +∞[ of the function:

ψ : ]0, +∞[ → R : ζ → ζ 2 2 -1 ln ζ + 1 2 γ -1 υ - 1 2 )ζ 2 + (1 -γ -1 ξ)ζ. (5.60)
Proof. For every (υ, ξ) ∈ R 2 , (υ, ξ) = prox γΦ (υ, ξ) is such that (υ, ξ) ∈ dom Φ Moreau [1962]. Let us first note that

υ ∈ ]0, +∞[ ⇔ (υ, ξ) ∈ ]0, +∞[ 2 .
(5.61)

We are now able to apply Proposition 5.3.3, where ψ is given by (5.60) and

(∀ζ ∈ ]0, +∞[) Θ(ζ) = ζ 2 2 1 2 -ln ζ -1 (5.62) ϑ -(ζ) = -ln ζ (5.63) ϑ + (ζ) = 1 -ζ -1 . (5.64)
In addition, χ -= exp(-υ/γ) (5.65) We have shown that υ > 0 ⇔ (5.69). So, υ = 0 when (5.69) is not satisfied. Then, the expression of ξ simply reduces to the asymmetric soft-thresholding rule Combettes and Pesquet [2007a]:

χ + = (1 -γ -1 ξ) -1 if ξ < γ +∞ otherwise. ( 5
ξ = ξ -γ if ξ > γ 0 otherwise.
(5.72)

However, exp(υ/γ) ≤ 1 -γ -1 ξ ⇒ ξ < γ, so that ξ is necessarily equal to 0.

Remark 5.4.2 (i) It can be noticed that

ψ ( ζ) = ζ ln ζ + γ -1 υ ζ - 1 ζ + 1 -γ -1 ξ = 0 ⇔ ζ -1 exp ζ -1 ( ζ -1 + γ -1 ξ -1) = exp(γ -1 υ). (5.73)
In the case where ξ = γ, the above equation reduces to

2 ζ -2 exp 2 ζ -2 = 2 exp(2γ -1 υ) ⇔ ζ = 2 W (2e 2υ/γ ) 1/2 (5.74)
where W is the Lambert W function [START_REF] Corless | On the Lambert W function[END_REF]. When ξ = γ, although a closed form expression of (5.73) is not available, efficient numerical methods to compute ζ can be developed as shown in the following.

(ii) For a given (υ, ξ) ∈ R 2 , minimizing ψ over ] exp(-υ/γ), +∞[ amounts to minimizing ψ 1 + ψ 2 over R, where which shows that ψ 1 and ψ 2 are strictly convex functions.

(∀ζ ∈ R) ψ 1 (ζ) = ζ 2 2 ln ζ + 1 2 γ -1 υ -1 2 ζ 2 if ζ ≥ exp(-υ/γ) +∞ otherwise (5.75) and (∀ζ ∈ R) ψ 2 (ζ) = -ln ζ + (1 -γ -1 ξ)ζ if ζ > 0 +∞ otherwise. ( 5
As ψ 1 and ψ 2 belong to Γ 0 (R), a possibility for performing this minimization is to employ the Douglas-Rachford algorithm. This requires calculating the proximity operators of τ ψ 1 and τ ψ 2 where τ > 0.

For every ζ ∈ R, prox τ ψ 1 ζ ∈ [exp(-υ/γ), +∞ [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]. We have then

ζ > exp(-υ/γ) ζ = prox τ ψ 1 ζ (5.79) ⇔ ζ > exp(-υ/γ) ζ -ζ + τ ζ ln ζ + γ -1 υζ = 0. (5.80) Since ζ → ζ + τ (ζ ln ζ + γ -1 υζ) is strictly increasing (hence, bijective) from [exp(-υ/γ), +∞[ to [exp(-υ/γ), +∞[, this is also equivalent to ζ > exp(-υ/γ) ζ -ζ + τ (ζ ln ζ + γ -1 υζ) = 0.
(5.81)

In addition, the latter equality can be rewritten as

ζ τ ζ exp ζ τ ζ = ζ τ exp(τ -1 + γ -1 υ) ⇔ ζ = ζ τ W τ -1 ζe τ -1 +γ -1 υ . (5.82)
In summary, (5.79) is equivalent to ζ > exp(-υ/γ) and (5.82). We deduce that

prox τ ψ 1 ζ =      ζ τ W τ -1 ζe τ -1 +γ -1 υ if ζ > exp(-υ/γ)
exp(-υ/γ) otherwise.

(5.83)

The proximity operator of τ ψ 2 follows from standard results on the calculation of proximity operators [START_REF] Pustelnik | Parallel ProXimal algorithm for image restoration using hybrid regularization[END_REF]:

prox τ ψ 2 ζ = 1 2 ζ + τ (γ -1 ξ -1) + |ζ + τ (γ -1 ξ -1)| 2 + 4τ . (5.84)
The Douglas-Rachford algorithm to compute the solution is then provided by:

where, in the first case, ζ is the unique minimizer on ]W (e 1-γ -1 υ ), +∞[ of the function:

ψ : ]0, +∞[ → R : ζ → ζ 2 2 + ζ -1 ln ζ + ζ 3 3 + 1 2 γ -1 υ - 3 2 ζ 2 -γ -1 ξζ.
(5.90)

Proof. We apply Proposition 5.3.3 where ψ is given by (5.90) and

(∀ζ ∈ ]0, +∞[) Θ(ζ) = ζ 2 3 4 - ζ 3 - 1 2 ln ζ (5.91) ϑ -(ζ) = ϑ + (ζ -1 ) = -ln ζ -ζ + 1.
(5.92)

The above equalities have been derived from (5.55) and (5.56). It can be deduced from (5.29), (5.30) and (5.92) that χ -+ ln χ -= 1 -γ -1 υ (5.93) W (e 1-γ -1 υ )W (e 1-γ -1 ξ ) < 1.

χ -1 + + ln(χ -1 + ) = 1 -γ -1 ξ (5.94) that is χ -= W (e 1-γ -1 υ ) (5.95) χ + = W (e 1-γ -1 ξ ) -1 . ( 5
(5.97) Under this assumption, the expression of the proximity operator follows from Proposition 5.3.3 and (5.92).

We have shown that prox γΦ (υ, ξ) ∈ ]0, +∞[ 2 ⇔ (5.97). Since prox γΦ (υ, ξ) ∈ dom Φ, we necessarily get prox γΦ (υ, ξ) = (0, 0), when (5.97) is not satisfied.

Remark 5.4.4 To minimize ψ, we need to find the zero on [χ -, χ + ] of the function:

ψ (ζ) = (ζ + 1) log ζ + ζ 2 + 1 - 1 ζ + ζ 2 + γ -1 υ - 3 2 ζ -γ -1 ξ. (5.98)
Thus can be performed by Algorithm 9, the convergence of which is proved in Section 5.7.

Algorithm 9 Newton algorithm

Hellinger divergence

Let us now consider the function of Γ 0 (R 2 ) given by

Φ : (υ, ξ) → ( √ υ - √ ξ) 2 if (υ, ξ) ∈ [0, +∞[ 2 +∞
otherwise.

(5.99)

This symmetric function satisfies (5.9) and (5.55) with

(∀ζ ∈ ]0, +∞[) ϕ(ζ) = ζ -ζ.
(5.100) Proposition 5.4.5 The proximity operator of γΦ with γ ∈ ]0, +∞[ is

(∀(υ, ξ) ∈ R 2 ) prox γΦ (υ, ξ) =            υ + γ(ρ -1), ξ + γ 1 ρ -1 if (υ < γ and (1 -γ -1 υ)(1 -γ -1 ξ) < 1) or υ ≥ γ (0, 0) otherwise (5.101)
where, in the first case, ρ is the unique solution on ] max(1 -γ -1 υ, 0), +∞[ of the equation:

ρ 4 + (γ -1 υ -1)ρ 3 + (1 -γ -1 ξ)ρ -1 = 0.
(5.102)

Proof. For every (υ, ξ) ∈ R 2 , (υ, ξ) = prox γΦ (υ, ξ) is such that (υ, ξ) ∈ [0, +∞[ 2 Moreau [1962]. By using the notation of Proposition 5.3.3 and by using Remark 5.3.4, we have:

(∀ζ ∈ ]0, +∞[) Θ(ζ) = ζ 2 2 - 2 5 ζ 5/2 + 1 (5.103) ϑ -(ζ) = ϑ + (ζ -1 ) = 1 -ζ (5.104) and χ -= (1 -γ -1 υ) 2 if υ < γ 0 otherwise (5.105) χ + = (1 -γ -1 ξ) -2 if ξ < γ +∞ otherwise. ( 5
ψ : ]0, +∞[ → R : ζ → 2 5 ζ 5/2 -2ζ 1/2 + γ -1 υ -1 2 ζ 2 + (1 -γ -1 ξ)ζ. (5.109)
This means that ζ is the unique solution on ]χ -, +∞[ of the equation: or, equivalently

ψ ( ζ) = ζ 3/2 -ζ -1/2 + (γ -1 υ -1) ζ + 1 -γ -1 ξ = 0. ( 5
υ < γ and (1 -γ -1 υ)(1 -γ -1 ξ) < 1 or υ ≥ γ. (5.112)
In turn, when (5.112) is not satisfied, we necessarily have υ = 0 or ξ = 0. In the first case, the expression of ξ is simply given by the asymmetric soft-thresholding rule in (5.72). Similarly, in the second case, we have

υ = υ -γ if υ > γ 0 otherwise.
(5.113) However, when υ > γ or ξ > γ, (5.111) is always satisfied, so that υ = ξ = 0. Altogether, the above results yield the expression of the proximity operator in (5.101).

Chi square divergence

Let us now consider the function of Γ 0 (R 2 ) given by

Φ : (υ, ξ) →          (υ -ξ) 2 ξ if υ ∈ [0, +∞[ and ξ ∈ ]0, +∞[ 0 if υ = ξ = 0 +∞ otherwise.
(5.114)

This function satisfies (5.9) with

(∀ζ ∈ ]0, +∞[) ϕ(ζ) = (ζ -1) 2 .
(5.115) onto the epigraph of a convex function is a useful tool for solving constrained convex optimization with the help of proximal algorithms. Epigraphical projection may be of main interest when dealing with constraint sets for which a closed form expression of the projection does not exist. In this part, we briefly explain how the forms of the proximity operator of divergences which have been derived in the previous sections can be employed to deduce the expression of the projection onto the epigraph of some convex functions. (5.150)

Let epi ϕ * = (υ * , ξ * ) ∈ R 2 ϕ * (υ * ) ≤ ξ * be the epigraph of ϕ * .
Then, the projection onto epi ϕ * is given by The above proposition is applicable to all the examples of ϕ-divergences wich have been presented in Section 5.3. In Table 5.1, we give the expressions of the corresponding functions ϕ * . In this chapter, we have shown how to deal with convex optimization problems involving discrete information divergences. Our analysis has emphasized that the difficulty lies in computing the proximity operator of two-variable real functions. To this end, we have carried out a systematic study of the proximity properties of ϕ-divergences, which has eventually led us to derive new closed-form expressions of proximity operators. These novel results have shed light on a new methodology that may replace the approaches usually followed to deal with problems involving these divergences (i.e. minimizing by fixing one of the two arguments or by alternating between the two variables).

∀(υ * , ξ * ) ∈ R 2 P epi ϕ * (υ * , ξ * ) = (υ * , -ξ * ) -prox Φ (υ * , -ξ * ). (5.151) Proof. The conjugate function of Φ is (∀(υ, ξ) ∈ R 2 ) Φ * (υ * , ξ * ) = sup (υ,ξ)∈R 2 υυ * + ξξ * -Φ(υ, ξ) . ( 5

Divergence

ϕ(ζ) ϕ * (ζ * ) ζ > 0 ζ * ∈ R Kullback-Leibler ζ ln ζ -ζ + 1 e ζ * -1 Jeffreys-Kullback (ζ -1) ln ζ W (e 1-ζ * ) + W (e 1-ζ * ) -1 + ζ * -2 Chi square (ζ -1) 2    ζ * (ζ * + 4) 4 if ζ * ≥ -2 -1 otherwise Hellinger 1 + ζ -2 √ ζ    ζ * 1 -ζ * if ζ * > -1 +∞ otherwise Renyi, α ∈]1, +∞[ ζ α    (α -1) ζ * α α α-1 if ζ * ≥ 0 0 otherwise Iα, α ∈]0, 1[ 1 -α + αζ -ζ α    (1 -α) 1 - ζ * α α α-1 -1 if ζ * ≤ α +∞ otherwise
In addition, we have provided a closed-form expression of the projection onto the epigraph of convex functions related to the considered ϕ-divergences.

In the next chapter, we will consider applications of these divergences, which may be useful in some problems, in the context of image restoration and image registration. § 5.7 Convergence proof of algorithm 9

We aim at finding the zero of the function ψ given below (along with its derivative ψ ):

∀ ζ ∈ ]W (e 1-γ -1 υ ), +∞[ ψ (ζ) = (ζ + 1) log ζ + ζ 2 + 1 - 1 ζ + ζ 2 + γ -1 υ - 3 2 ζ -γ -1 ξ, ψ (ζ) = log ζ + 1 ζ + 1 ζ 2 + 2ζ + γ -1 υ.
To do so, we employ the Newton algorithm, whose global convergence is guaranteed for any initial value by the following condition Thorlund-Petersen [2004]:

(∀a ∈]0, +∞[)(∀b ∈]a, +∞[) ψ (a) + ψ (b) > ψ (b) -ψ (a) b -a , which is equivalent to (b -a)ψ (a) + (b -a)ψ (b) -ψ (b) + ψ (a) > 0.
(5.157) Condition (5.157) can be rewritten as follows

(b -a)(log a + 1 a + 1 a 2 + 2a + γ -1 υ) +(b -a)(log b + 1 b + 1 b 2 + 2b + γ -1 υ) + (a + 1) log a + a 2 + 1 - 1 a + a 2 + γ -1 υ - 3 2 a -γ -1 ξ -(b + 1) log b + b 2 + 1 - 1 b + b 2 + γ -1 υ - 3 2 b -γ -1 ξ > 0,
which, after some simplification, boils down to To prove these two inequalities, we will study the derivative of g with respect to its arguments. The conditions in (5. The first and second partial derivatives of g w.r.t. x read

(b + 1) log a -(a + 1) log b + b a + b a 2 - 2 a - a b - a b 2 + 2 b -a 2 + b 2 +γ -1 υ(b -a) + 1 2 (b -a) > 0. ( 5
(∀y ∈]0, +∞[)(∀x ∈]0, y[) ∂g ∂x (x, y) = -log y - 1 y - 2y x 3 + 2 x 2 ∂ 2 g ∂x 2 (x, y) = 6y x 4 - 4 x 3 = 6y -4x a 4 > 0.
Since ∂ 2 g/∂x 2 is strictly positive, ∂g/∂x is strictly increasing w.r.t. x and lim

x→y ∂g ∂x (x, y) = -log y - 1 y = log 1 y - 1 y < 0.
Therefore, Condition (5.162) holds, and g is decreasing with respect to x.

The first and second partial derivatives of g w.r.t. y read

(∀x ∈]0, +∞[)(∀y ∈]x, +∞[) ∂g ∂y (x, y) = -x y - 1 y + x y 2 + 1 x 2 + 2y ∂ 2 g ∂y 2 (x, y) = x y 2 + 1 y 2 - 2x y 3 + 2. For every y ∈ [1, +∞[, (∀x ∈]0, y[) ∂ 2 g ∂y 2 (x, y) = x y 2 + 1 y 2 - 2x y 3 <1 +2 > 0,
and ∂g/∂y is strictly increasing w.r.t. y (since ∂ 2 g/∂y 2 is strictly positive) and

(∀x ∈ [1, +∞[) lim y→x ∂g ∂y (x, y) = -1 + 1 x 2 + 2x > 0, (∀x ∈]0, 1]) ∂g ∂y (x, 1) = 1 x 2 + 1 > 0.
For every y ∈ ]0, 1[, we have

(∀x ∈]0, y[) ∂g ∂y (x, y) = -x y + x y 2 - 1 y + 1 x 2 + 2y, = x -xy y 2 - 1 y + 1 x 2 + 2y, since x < y < 1, this implies that xy < x and 1 y < 1 x < 1 x 2 and (∀x ∈]0, y[) ∂g ∂y (x, y) = x -xy y 2 >0 - 1 y + 1 x 2 >0 +2y > 0.
As Condition (5.163) holds, g is increasing with respect to y.

-Chapter 6 -

2D Applications of divergence proximity operators

"Knowledge without application is like a book that is never read."

Christopher Crawford

Measures based on distance or divergence play a pivotal role in evaluating the dissimilarity between two objects (such as numbers, vectors, matrices, functions, images, . . . ). In this chapter, we investigate the use of divergences as a data fidelity term or a regularization term in two different image processing applications: restoration and registration.

As we have presented in the previous chapter, the objective function will have the following form:

x → D(Ax, Bx) (6.1)
where D is a function in Γ 0 (R P × R P ), and A and B are matrices in R P ×N . The most common assumption on the signal of interest x ∈ R N is that, by making an appropriate choice of matrices A and B, Ax and Bx are close in the metric induced by D.

Divergences have been widely used in the literature within many different contexts. Some measures have been used for regression analysis, some have been applied to compare two probability density functions, and some are used to measure the dissimilarity between two random vectors of the same distribution. In this chapter, we use information measures in two applications. Firstly, we extend the NLTV framework by using some information divergences to build new sparsity measures for signal recovery. In denoising and deconvolution examples, our approach is compared with other existing approaches. Secondly, an image registration application constitutes another example illustrating the use of the proposed divergence proximity operators. In Equation (3.48), we have defined NLTV as follows:

NLTV(x) = s∈A n∈Ns⊂Ws ω s,n |x (s) -x (n) | 2 , (6.2)
where W s is the set of positions n ∈ A \ {s} located into a Q × Q window centered at s, where Q ∈ N is odd. The idea behind NLTV is illustrated in Fig. 6.1. For each pixel in the image, we select some neighbours within a limited region and we compute the weights ω s,n according to a similarity measure between Q × Q patches built around pixels s and n. Large weights are desirable in flat areas in order to effectively suppress the noise, while small weights are desirable around object edges in order to preserve details and fine structures.

A key observation is that NLTV can be interpreted as a similarity measure. Indeed, the vector of weighted differences can be rewritten as:

(∀s ∈ A) [ω s,n (x (s) -x (n) )] n∈Ns = A s x -B s x, (6.3)
where

A s : x → [ω s,n x (s) ] n∈Ns , (6.4)
B s : x → [ω s,n x (n) ] n∈Ns . (6.5)
Consequently, the 1,p -NLTV regularization term can be expressed as follows b (s,m) , (6.9)

1,p -NLTV(x) = s∈A A s x -B s x p . ( 6
where ϕ : R → [0, +∞] belongs to Γ 0 (R) and is twice differentiable on ]0, +∞[. The aim of this part is to investigate the influence of the choice of the similarity measure on the performance of optimization approaches in the context of image recovery.

Numerical experiments

In order to assess the relevance of the proposed regularization, we present numerical experiments relative to an image restoration problem. The considered degradation model is z = Hx + w, where x ∈ R N denotes the image to be recovered, z ∈ R M are the observed data, H ∈ R M ×N corresponds to a matrix associated with a blur operator, and w ∈ R M is a realization of an additive zero-mean white Gaussian noise. The restoration process is based on a variational approach aiming at finding an image x ∈ R N as close as possible to x ∈ R N from the observation z. Indeed, as inverse problems are usually ill-posed, one needs some prior information about the estimate, such as the dynamic range of the pixel values and the smoothness of the estimate of x. Therefore, the restoration can be achieved by solving the following convex optimization problem:

minimize x∈R N 1 2λ Hx -z 2 + D(Ax, Bx) + ι C (x) (6.10)
where λ ∈]0, +∞[ is the regularization constant, and the operators A and B are chosen as explained before. The last term ι C constrains x to belong to the convex set C = [0, 255] N . Alternatively, we seek to describe the problem within a set theoretic framework. In several works [START_REF] Combettes | The use of noise properties in set theoretic estimation[END_REF], it was observed that an upper bound on the data fidelity term allows us to efficiently restrict the solution to vectors x such that:

Hx ∈ C = u ∈ R M | u -z 2 ≤ δM σ 2 (6.11)
where δ is a positive constant (usually close to 1) and σ 2 is the noise variance. This leads to the following variant of Problem (6.10):

minimize x∈R N D(Ax, Bx) + ι C (Hx) + ι C (x). ( 6
.12)

Note that, under technical assumptions (Chapter 2, Section 2.1.1), Problem (6.12) is equivalent to Problem (6.10). However the constrained formulation given above is usually considered to be more practical as the solution is less sensitive to the choice of δ than λ Chierchia et al. [2013].

The above problem can be solved using proximal optimization algorithms. Such methods require to compute the proximity operator of the divergence D, the projection onto the 2 ball, and the projection onto the hypercube [0, 255] N . These two projections are quite standard. A possible algorithm for solving Problem (6.12) is thus the M+LFBF primal-dual algorithm [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallelsum type monotone operators[END_REF]. The associated iterations are recalled in Algorithm 1, where at each iteration k, γ [k] is a step-size and e [k] ∈ (R P ) 2 The observed image is generated by degrading the original one with a convolution operator H, which is equal to identity for the denoising problem and corresponds to a truncated Gaussian point spread function with standard deviation 1.6 and kernel size 3 × 3, and 7 × 7 for the deconvolution problem. The noise variance is equal to 400 and 64 for the image denoising and restoration problems, respectively. The linear operators A and B associated with NLTV are computed from the TV image result obtained using the code in [START_REF] Foi | Foveated self-similarity in nonlocal image filtering[END_REF]. For the regularization parameters, we set Q = 11 and Q = 5.

The balance between the smoothness of the estimate and the data fidelity is controlled by the parameter δ tuned so as to maximize the Signal-to-Noise Ratio (SNR). The quality of the results, presented in Tables 6.1,6.2, and 6.3, is evaluated in terms of the SNR and the Structural Similarity index (SSIM) [START_REF] Malpica | Range image quality assessment by structural similarity[END_REF]. One can observe that the results obtained with D-NLTV improve upon the standard 1,2 -NLTV, for the different divergences with a gain up to 0.44 dB. Fig. 6.3 and 6.4 show that our method using Div-NLTV leads to significant visual improvements compared with 1,2 -TV and 1,2 -NLTV. Note that the considered divergences lead to similar results. We believe that this behavior can be explained by the fact that the ϕ-divergences present little difference when the two arguments are similar, i.e. the ratio is close to one. Indeed, the Taylor expansion around one is the same and equal to the quadratic difference as shown next. Recall that the considered ϕ-divergences are additive information measures of the form given by Equation (5.8).

Let us now perform the second-order Taylor expansion of the function Φ around an initial estimate (ῡ, ξ) ∈ R 2 as follows:

Φ(υ, ξ) Φ(ῡ, ξ) + (υ -ῡ) ∂Φ ∂υ (ῡ, ξ) + (ξ -ξ) ∂Φ ∂ξ (ῡ, ξ) + 1 2 (υ -ῡ) 2 ∂ 2 Φ ∂υ 2 (ῡ, ξ) + 1 2 (ξ -ξ) 2 ∂ 2 Φ ∂ξ 2 (ῡ, ξ) + (υ -ῡ)(ξ -ξ) ∂ 2 Φ ∂υ∂ξ (ῡ, ξ), (6.13)
where

∂Φ ∂υ (υ, ξ) = ϕ υ ξ ∂Φ ∂ξ (υ, ξ) = ϕ υ ξ - υ ξ ϕ υ ξ ∂ 2 Φ ∂υ 2 (υ, ξ) = 1 ξ ϕ υ ξ ∂ 2 Φ ∂ξ 2 (υ, ξ) = υ 2 ξ 3 ϕ υ ξ ∂ 2 Φ ∂υ∂ξ (υ, ξ) = - υ ξ 2 ϕ υ ξ . (6.14)
Since ϕ is such that

ϕ(1) = ϕ (1) = 0, (6.15)
we can deduce that when ξ → ῡ, the equations in (6.14) boil down to

∂Φ ∂υ (ῡ, ξ) = 0, ∂Φ ∂ξ (ῡ, ξ) = 0, ∂ 2 Φ ∂υ 2 (ῡ, ξ) = 1 ῡ ϕ (1), ∂ 2 Φ ∂ξ 2 (ῡ, ξ) = 1 ῡ ϕ (1), ∂ 2 Φ ∂υ∂ξ (ῡ, ξ) = -1 ῡ ϕ (1). (6.16) Thus (6.13) reads Φ(υ, ξ) 1 2ῡ ϕ (1) (υ -ῡ) 2 + (ξ -ξ) 2 -2(υ -ῡ)(ξ -ξ) 1 2ῡ ϕ (1)(υ -ξ) 2 , (6.17)
where ϕ (1) > 0 since ϕ is strictly convex.

So, as we can notice, divergences perform similarly to the squared 2norm when the two arguments are close. The latter computation is validated in the same context of image restoration using the NLTV regularization. Tables 6.4, 6.5, and 6.6 show that 2 -NLTV leads to close performance to the ϕ-divergences ("Best of Div." designates the best result obtained in restoration using ϕ-divergences, i.e. Tables 6.1, 6.2, and 6.3). where ∂F j ∂z j (I j , zj ) is a Jacobian matrix. With the linearization expressed in (6.19), the image registration problem can be reformulated under the form of Problem 5.1 (Chapter 5), where

A = ∂F 1 ∂z 1 (I 1 , z1 ) 0 , (6.20)
B = 0 ∂F 2 ∂z 2 (I 2 , z2 ) , (6.21)
and

u = F 1 (I 1 , z1 ) -z1 ∂F 1 ∂z 1 (I 1 , z1 ), (6.22) v = F 2 (I 2 , z2 ) -z2 ∂F 2 ∂z 2 (I 2 , z2 ). ( 6.23)
The additional regularization terms (R s •T s ) 1≤s≤S allow the incorporation of prior knowledge about the sought parameter vector x.

Disparity estimation under illumination variation

We conducted experiments for disparity estimation under illumination variation, using a stereoscopic pair of grayscale images I 1 and I 2 with size P 1 × P 2 (P = P 1 P 2 ). The related image mapping operators are given by F 1 (I 1 , z 1 ) = vec (z

(i 1 ,i 2 ) 1 I (i 1 ,i 2 ) 1 ) 1 ≤ i 1 ≤ P 1 1 ≤ i 2 ≤ P 2 (6.24) F 2 (I 2 , z 2 ) = vec (I (i 1 -z (i 1 ,i 2 ) 2 ,i 2 ) 2 ) 1 ≤ i 1 ≤ P 1 1 ≤ i 2 ≤ P 2 (6.25)
where

I (i 1 ,i 2 ) 1 (resp. I (i 1 ,i 2 ) 2
) is the intensity value in the left (resp. right) view at pixel (i 1 , i 2 ), z (i 1 ,i 2 ) 1 the illumination variation coefficient at (i 1 , i 2 ), and z (i 1 ,i 2 ) 2 the associated disparity.

The parameter vectors are here z 1 = vec (z

(i 1 ,i 2 ) 1
) 1≤i 1 ≤P 1 ,1≤i 2 ≤P 2 and z 2 = vec (z

(i 1 ,i 2 ) 2 ) 1≤i 1 ≤P 1 ,1≤i 2 ≤P 2 .
It can be observed that ∂F 1 ∂z 1 (I 1 , z1 ) = Diag (I (i 1 ,i 2 ) 1

) 1 ≤ i 1 ≤ P 1 1 ≤ i 2 ≤ P 2 (6.26) ∂F 2 ∂z 2 (I 2 , z2 ) = -Diag (∇ (1) I (i 1 -z

(i 1 ,i 2 ) 2 ,i 2 ) 2 ) 1 ≤ i 1 ≤ P 1 1 ≤ i 2 ≤ P 2 (6.27)
where ∇ (1) is the gradient operator w.r.t. the first spatial coordinate. It is useful to incorporate prior information about the solution so as to 6.2. Data fidelity based on ϕ-divergences 129 convert the original matching problem to a well-posed one. Similarly to what we have presented in Chapter 3, a first constraint set is introduced that takes into account the range of possible values for the disparity and the illumination variation:

C 1 = [z 1 , z 2 ] ∈ R 2P | (z 1 ) min ≤ z 1 ≤ (z 1 ) max , (z 2 ) min ≤ z 2 ≤ (z 2 ) max .
(6.28)

A second constraint is employed to enforce the smoothness of the estimated fields:

C 2 = {[z 1 , z 2 ] ∈ R 2P | ∇z 1 2 2 ≤ κ 1 , TV(z 2 ) ≤ κ 2 } (6.29)
where • 2 is the 2 norm, ∇ is the spatial gradient operator and TV is the discrete total variation semi-norm. Consequently, we have now S = 2 in Problem 5.1 (previous chapter), where R 1 and R 2 reduce to indicator functions of convex sets, T 1 = I 2P , and T 2 is a block diagonal matrix with diagonal terms set to gradient operators.

Experiments

We now illustrate the practical performance of our method on the "Cloth", "Baby 1", "Baby 2", and "Dolls" stereo pairs (Figure 6.5) downloaded from the MiddleBury website. 3The results are provided in Table 6.7 and Figure 6.6. The quality of the results was evaluated using three different metrics based on a ground truth: the MAE (Mean Absolute Error) evaluated between the computed maps (z j ) j∈{1,2} and the ground truth ( z j ) j∈{1,2} , the percentage of bad matching pixels in the disparity map Err T = 1 P

P 1 i 1 =1 P 2 i 2 =1
1(|z

(i 1 ,i 2 ) 2 -z (i 1 ,i 2 ) 2
)| > T ), Finally, in Figure 6.6, we illustrate the recovered disparity maps using JK and KL divergences for the gray stereo pair "Cloth". These visual results show the efficiency of the JK measure in recovering the depth estimation in the presence of the illumination variations. In this chapter, we have applied ϕ-divergences as a measures of dissimilarity in two different contexts:

• We have extended the NLTV framework by using information divergences to build new sparsity regularization measures for signal recovery.

• We have used the divergence as a data fidelity term for image registration, more precisely for disparity estimation under illumination variation.

However, as we have shown before, the performance achieved by different types of ϕ-divergence is similar. This behaviour may be explained by the fact that the values of these two-variable functions have the same second-order Taylor expansion when the two arguments are close to each others. The practical results obtained by implementing the proposed ϕ-divergence measures in an image processing context do not yield the expected performance gains with respect to standard measures such as the squared 2 -norm. However, it would be interesting to investigate the effectiveness of divergences in problems involving probability distribution. Indeed, it is our belief that these measures could constitute a successful element for other applications and that they are undoubtedly better adapted to some other problems. -Chapter 7 -

Conclusion

"You may not end up where you thought you'd be, but you'll end up right where you're meant to be."

In the first part of this thesis, we have investigated the problem of disparity estimation within the framework of convex optimization and extended it to the case of multiple stereo views. In the second part, we have provided novel tools to deal with convex divergence measures in optimization problems and we have applied them to image recovery problems. § 7.1 Contributions Disparity estimation under illumination variation We have proposed a convex optimization approach to deal with disparity estimation under illumination variation. Our solution consists of modeling the spatiallyvarying illumination variation by a multiplicative term. Since the global formulation is nonconvex, the problem is addressed by solving a sequence of convex relaxations.

• We have employed a relaxation based on a first-order Taylor approximation around an initial estimate to derive a convex energy function. This energy is then minimized subject to various convex constraints aiming at imposing some regularity on the solution.

The originality in our approach lies in the ability to deal with nonstrictly convex energies, which allows us to employ data fidelity terms other than quadratic ones.

• The problem of jointly estimating the disparity and the illumination variation is tackled in a set theoretic framework and solved using proximal tools. The original contributions of our work consist of the ability to consider multicomponent images in the presence of illumination variation, where each component corresponds to one of the channels of a specific color system (RGB, YUV, ...).

The proposed approach is also able to exploit the potentials offered by multicore/GPU parallel system architectures. The experiments indicate that this approach can effectively deal with the local illumination variation and yields better results compared with other approaches.

Disparity map sequence We have extended the previous approach to multi-view problems. Rather than estimating a single depth map, we have associated a depth map with each input image (or a subset of them). We have proposed two different approaches to estimate the depth map for each chosen reference frame, which differ in the type of convex relaxation employed. The first approach employs a relaxation based on multi-labeling, with the possibility of using nonconvex similarity measures. So doing, there is no more dependence with respect to the initial estimate, at the cost of an increased complexity. The second approach employs a relaxation based on an improved first-order Taylor approximation around an initial estimate. The experimental results illustrate, through a complex specular scene example, the effectiveness of our proposals.

ϕ-divergence proximity operator We have shown how to deal with convex optimization problems involving ϕ-divergences, when the minimization is jointly performed with respect to both divergence arguments. More specifically, we have derived new expressions for the proximity operator of several ϕ-divergences (e.g. Kullback-Leibler, Jeffreys-Kullback, Hellinger, Chi-Square, I α , and Renyi divergences), enriching the list of functions for which such operators can be easily computed. This has allowed us to address a generic class of convex optimization problems by means of efficient proximal algorithms, overcoming the limitations of current state-of-the-art approaches for such problems.

Applications for ϕ-divergence proximity operator We have investigated the use of divergences in two different image processing applications: restoration and registration. Firstly, we have extended the NLTV framework by using some information divergences to build new regularization terms for signal recovery. This regularity measure allows us to better exploit image features by making use of the relevant information along an edge or a regular texture pattern. Secondly, we have presented an image registration application to illustrate the use of divergences as data fidelity terms, where our objective is to determine spatial transformations that maximize a similarity metric between two images resulting from two different acquisitions.

Epigraphical projection Epigraphical projections have a great interest when dealing with constraint sets for which a closed form expression of the projection does not exist. Within this context, we have provided theoretical guidelines to show how the forms of the proximity operator of divergences can be employed to compute the expression of the projection onto the epigraph of a number of convex functions. § 7.2 Perspectives Extension to View synthesis application: The problem of estimating the disparity between two stereo images (see Chapter 3) can be extended to view synthesis. Because only a few cameras are available, additional views have to be synthesized from them. The main idea is to simultaneously estimate a new view and the corresponding disparity map. To do so, we need to synthesize a virtual depth map from the current view by applying a warping technique. We do not need any other side information to synthesize the virtual depth maps except camera parameters. In particular, it could be a good idea to investigate new ways of estimating the illumination variation for the new view synthesis by studying how the light could change from one image to another in multi-view datasets. This could be applied later in 3D tracking, stereo video coding, 3D scene interpretation, and 3D television.

Exploiting the dependence among the disparity sequence maps:

The multi-view methods presented in Chapter 4 can be improved in several ways. One of the most appealing ideas is to combine the disparity sequence estimation in a single optimization problem in order to exploit the huge correlation between the disparities of different views.

In this case, an additional piece of prior information may be included into the minimization problem: the 3D point projected in different views should have the same disparity value (since cameras are rectified and equidistant). The main challenges in this context are the occlusion areas and the non-convexity of the considered constraint, which needs to be relaxed using one of the techniques presented in Chapter 6.

Combining the discrete and continuous methods: There is also room for improving the potential of the multi-label approach for the disparity estimation sequence. The images of a scene under varying illuminations and from different viewpoints are highly interrelated, which make it possible to predict the object depth and the illumination variation. Illumination due to imperfectly calibrated cameras, different perspective projection directions, and different reflection effects varies from an object to the other. We could consider the local illumination variation that arises due to differences in camera positioning. Hence, it would be interesting to extend the multi-label approach in order to jointly estimate the disparity map and the illumination variation field. This would require combining the multilabel-based relaxation of Chapter 4 with the continuous approach of Chapter 3.

Disparity and Motion from a Multi-view video sequence: Stereo video sequences show the temporal evolution of a scene from two different viewpoints. In order to allow for an accurate stereo sequence analysis, it is essential to consider the spatio-temporal relationships that exist between the different images of the sequences. The stereo-motion consistency constraint can be exploited by jointly estimating the motion field of the left and right sequences, along with the disparity field of the current stereo pair. This approach can be efficiently extended to a multi-view sequence by using the techniques presented in Chapter 4.

ϕ-divergence in blind deconvolution ϕ-divergences can be applied in the context of blind deconvolution. Here, the degradation model is the same as in Equation (2.1):

z = T x + w,
but one aims at estimating both the target image x ∈ R N and the linear operator T ∈ R K×N . The common approach is to resort to an alternating minimization scheme, which however has no guaranteed convergence. An alternative approach consists of assuming that there exists an operator G such that GT = I, which allows one to reformulate the degradation model as

Gz = x + ω,
where ω = Gw denotes the transformed (and possibly amplified) noise. This leads to a convex optimization problem which aims at minimizing the divergence between Gz and x, while enforcing some kind of regularization on G and x. The presented approach constitutes a good example of minimization problems jointly performed on two variables.

ϕ-divergence for channel estimation It is worth mentioning that several applications in telecommunication systems could benefit from the proposed optimization approaches, e.g. for jointly estimating the transmitted data and the channel transfer function in the presence of noise.

In this context, ϕ-divergences play a central role in the formulation of the recovery problem. In addition, some constraints can be added to the problem in order to model the temporal evolution of the channel (due to Doppler effect). The resulting optimization problem is nonconvex and the existing methods generally follow an alternating minimization scheme (see e.g. [START_REF] Hu | A divergence minimization approach to joint multiuser decoding for coded CDMA[END_REF]). An alternative approach would consist of adopting the algorithms developed in Chapter 5 to address this nonconvex scenario, or by resorting to the convex reformulation provided above for blind deconvolution.

Epigraphical projection in video coding As we have seen in Section 5.5, the proximity operator of the considered ϕ-divergences allows us to compute the orthogonal projection onto the epigraph of functions listed in Table 5.1. Recently, epigraphical projections have been employed to address convex constraints for which the associated projection operator is not available in closed form [START_REF] Chierchia | Epigraphical splitting for solving constrained convex optimization problems with proximal tools[END_REF]. Our results thus substantially enrich the list of functions that can be handled by such an epigraphical approach. For example, the exponential function can be used in the context of video coding, where one needs to solve a rate allocation problem that aims at minimizing the rate budget while keeping the global distortion (expressed as a sum of exponentials composed by a linear operator) below a given bound. 4.7 (MAE, Err 2 ) obtained using five images for each disparity map estimation, Q=80, r 0 =11, r Q =91 and (∀ q ∈ {2, . . . , Q}), r qr q-1 = 1. For comparison, we indicate the values of (MAE, Err 2 ) stereo obtained in the stereo case when u n,kn is estimated from (I n ,I kn ). ψ is the 1 -norm.

List of Figures

Figure

  Figure 1.1

  Euclidean geomet r y I nfor mat ion geomet r ies Euclidean dist ance d2(p, q) = i (pi -qi ) 2 (Pyt hagoras' t heorem circa 500 BC) Minkowski dist ance (L k -norm) dk (p, q) = k i |pi -qi | k (H. Minkowski 1864-1909) Space-t im e geom et r y Manhat t an dist ance d1(p, q) = i |pi -qi | (city block-t axi cab) Mahalanobis met ric (1936) dΣ = (p -q) T Σ -1 (p -q) Quadrat ic dist ance dQ = (p -q) T Q(p -q) Riemannian met ric t ensor gi j d x i d s d x j d s ds (B. Riemann 1826-1866,) Physics ent ropy J K -1 -k p log p (Bolt zmann-Gibbs 1878) Informat ion ent ropy H (p) = -p log p (C. Shannon 1948) Fisher informat ion (local ent ropy) I (θ) = E[ ∂ ∂ θ ln p(X |θ) 2 ] (R. A. Fisher 1890-1962) K ullback-Leibler divergence K L(p||q) = p log p q = Ep [log P Q ] (relat ive ent ropy, 1951) -1) ln p α q 1-α (addit ive ent ropy) T sallis ent ropy (1998) (Non-addit ive ent ropy) Tα

G

  dist ance (1967) d(p, q) = -ln ( √ p -√ q) 2 N on-addit ive ent r opy cr o ss-en t r o p y co n d i t i o n a l en t r o p y m u t u a l i n f o r m a t i o n ( ch a i n r u l es) A ddit ive ent r opy Permissible Bregman divergences (Nock & Nielsen, 2007) N on-Euclidean geomet r ies St at ist ical geomet r y Je rey divergence (Jensen-Shannon) H ( p ) = K L ( p | | u ) ener alized P yt hagor as' t heor em (Generalized projection) I -project ion Quant um geomet r y ex p o n en t i a l f a m i l i es Log Det divergence D (P ||Q) = < P , Q -1 > -log det P Q -1 -dimP Von Neumann divergence D (P ||Q) = Tr(P (log P -log Q) -P + Q) It akura-Sait o divergence IS(p|q) = i ( pi qi -log pi qi -1) (Burg ent ropy) A lgor it hmic geomet r y? K olm ogor ov com plexit y K u l l b a ck -L ei b l er Distance between two algorithms ?

Figure 1 . 2 :

 12 Figure 1.2: Taxonomy of principal distances. c Frank Nielsen, 2007. Version 0.2 Last updated, May 2008

( a )Figure 2 . 1 :

 a21 Figure 2.1: Example of a degraded image resulting from the direct problem: z = T x + w.

Forward-

  Backward splitting Let R = 1, T 1 = Id and S = 0 in (2.21). Then the problem becomes minimize x∈R N h(x) + g(x).

  55) is decreasing in λ Parikh and Boyd [2013]. In sparse reconstruction applications, some interesting algorithms have been developed such as spectral projected gradient techniques Van Den Berg and Friedlander [2008], 5 1 -regularized least squares methods Koh et al. [2007], 6 sparse reconstruction by separable approximation Wright et al. [2009], 7 and many others Duchi et al. [2008], Schmidt et al. [2009], Weiss et al. [2009], Fadili and Peyré [2011]. Note that there also exist numerical methods that do not follow the aforementioned approach, such as Quattoni et al. [2009] and Kyrillidis et al. [2013]. § 2.

Figure 3 . 1 :

 31 Figure 3.1: If you focus on a tree, closing first your right eye and then your left eye, the tree will appear to move relatively to your finger.

Figure 3 . 2 :

 32 Figure 3.2: Cameras (at two different positions) capture two separate images.The disparity is the difference between the positions of the same object in the two images. If an object is far away, the disparity will be small; otherwise, the disparity will be large if the object is close.

Figure 3 . 3 :

 33 Figure 3.3: The epipolar lines are defined by the intersection of the left and right image planes with the epipolar plane (C L , C R , P ). The points e L and e R are the epipoles and are defined by the intersection of the image planes with the line connecting the two cameras: the baseline (C L , C R ).

Figure 3 . 4 :

 34 Figure 3.4: Two cameras (acquiring images of the same scene) have two different 2D representations of a common 3D point. With proper processing, the position and depth of the 3D point can be extracted from the images.

Figure 3 .

 3 Figure 3.5 depicts the camera configuration. By applying Thales theorem to the triangles of hypotenuses C L P and C R P , we get:

Figure 3 . 5 :

 35 Figure 3.5: The connection between disparity and depth.

  [START_REF] Figueiredo | An EM algorithm for wavelet-based image restoration[END_REF];[START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF];[START_REF] Raguet | Generalized forward-backward splitting[END_REF]; Briceño-Arias andCombettes [2011];[START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallelsum type monotone operators[END_REF];[START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. Similarly to the subgradient based methods, we use the multiplicative model by Gennert, which allows us to define a global similarity measure to be minimized over an intersection of convex constraint sets. The original contributions of our work consist of (i) the ability to consider multicomponent images in the presence of illumination variation (ii) the flexibility in minimizing various (possibly nonsmooth) convex similarity measures and combining them with various convex constraints.

Figure 3 . 6 :

 36 Figure 3.6: Stereo images: "Books" (left), "Dolls" (center), and "Parking meter" (right).

Figure 3 . 7 :Figure 3 . 8 :

 3738 Figure 3.7: Results for "Books" stereo pair: a)-b) ground truths, c)-d) subgradient projection method (c) disparity and d) illumination fields), e)f ) proposed approach with (u min , u max ) = (20, 75),(v min , v max ) = (0.1, 1.1), (τ 2 , κ 2 ) = (74000, 400) (e) disparity and f ) illumination fields), g) estimation of the disparity by NCC, h) estimation of the disparity after histogram equalization.

Figure 3 . 10 :

 310 Figure 3.10: Results for "Parking meter" stereo pair. a)-b) Subgradient projection method and c)-d) our approach with (u min , u max ) = (1, 8), (v min , v max ) = (0.1, 1.1), (τ 2 , κ 2 ) = (20000, 0.4). a)-c) is the measured illumination field, b)-d) the measured disparity field.

Figure 3 . 11 :

 311 Figure 3.11: Results for Dolls stereo pair: e)-f ) results from gray scale images and g)-h) from color (YUV) images. ((u min , u max ) = (20, 75), (v min , v max ) = (0.1, 1.1), (τ 2 , κ 2 ) = (92000, 240)).

  (a) Venus (b) Bull (c) Cloth (d) Wood

Figure 3 .

 3 Figure 3.12: Left images of the stereo pairs.

  Left image (b) True disparity (c) MAE=0.82, Err2=1.73 (d) Right image (e) True illumination (f ) MAE=0.067 (g) Left image (h) True disparity (i) MAE=2.53, Err2=8.10 (j) Right image (k) True illumination (l) MAE=0.015

Figure 3 . 13 :

 313 Figure 3.13: Visual results for with color images "Cloth" and "Wood". "Cloth" ( (a) to (f ) ) with (u min , u max ) = (15, 55),(v min , v max ) = (0.1, 0.6) and (τ 4 , κ 2 ) = (18800, 135). "Wood" ( (g) to (l) ) with (u min , u max ) = (21, 71),(v min , v max ) = (0.1, 0.6) and (τ 4 , κ 2 ) = (21000, 42).

Figure 4 .

 4 Figure 4.1: N images (view 1, view 2, . . . ) captured by N aligned and equidistant cameras (C 1 , C 2 , . . . ) at N different positions, resulting in a set of stereo pairs with equal baselines B.

  Guidelines for addressing the above optimization problems are provided in Chambolle and Pock [2011]; Briceño-Arias and Combettes [2011]; Combettes and Pesquet [2012]; Setzer et al. [2010]; Parikh and Boyd [2013] by employing algorithms involving computation of proximity operators.

2

  associated with the lower-level set constraint on h and E 2 . 68 Chapter 4. Multi-view disparity estimation

  .30) where, by using (4.3), (4.25), and (4.29), for every s = (i 1 , i 2 ) ∈ D n,m , = I (i 1 -αn,m ūn,kn (s),i 2 )

Figure 4 . 2 :

 42 Figure 4.2: Tested image sequences

Figure 4 . 3 :

 43 Figure 4.3: The disparity map u 2,4 and the indices (MAE, Err 1 ) obtained by using the proposed method with the truncated 1 -norm. From left to right: stereo, 3 views, 5 views. From top to button: Sawtooth, Teddy, and Cones.

Figure 4 . 4 :

 44 Figure 4.4: "Sponza" sequence. Each image is of size 480 × 960 (PPM grayscale format). Calibration parameters for the equidistant cameras are known.

Figure 4 . 5 :

 45 Figure 4.5: Image 3 of "Sponza" sequence.

Figure 4 . 6 :

 46 Figure 4.6: Visual comparison of "Sponza" disparity map u 3,4 estimated by the considered methods using 5 images.

u1, 2

 2 : (1.07, 4.43) u2,3: (0.94, 2.95) u3,4 : (0.94, 3.18) (1.56, 5.58)stereo (1.43, 5.17)stereo (1.55, 5.35)stereo u4,5: (1.01, 3.68) u5,4: (1.13, 5.58) (2.16, 6.69)stereo (2.33, 7.59)stereo

Figure 4 .

 4 Figure 4.7:(MAE, Err 2 ) obtained using five images for each disparity map estimation, Q=80, r 0 =11, r Q =91 and (∀ q ∈ {2, . . . , Q}), r q -r q-1 = 1. For comparison, we indicate the values of (MAE, Err 2 ) stereo obtained in the stereo case when u n,kn is estimated from (I n ,I kn ). ψ is the 1 -norm.

Figure 4 . 8 :

 48 Figure 4.8: The extended disparity map constructed using the sequence of disparity maps: (u 1,2 , u 2,3 , u 3,4 , u 4,5 , u 5,4 )

  ; Piro et al. [2008]; Dupé et al. [2009]; Zanella et al. [2009]; Piro et al. [2010]; Pustelnik et al. [2011]; Teuber et al. [2012]; Carlavan and Blanc-Féraud [2012]; Anthoine et al. [2012] more general forms of the regularization functions are considered. In particular, the developments in Dupé et al. [2009]; Pustelnik et al. [2011]; Teuber et al. [2012]; Anthoine et al.

  [2007]; Combettes and Pesquet [2011]). Some examples of such non-separable functions are: a quadratic function composed by a linear operator Combettes and Pesquet [2008], a separable function composed with an orthonormal or semi-orthogonal linear operator Combettes and Pesquet [2008], the 2 -norm Combettes and Pesquet [2008], the quadratic-over-linear function Benamou and Brenier [2000]; Benamou et al. [2002], the indicator function of a closed convex Combettes and Pesquet [2008] or the epigraph of a lower semi continuous convex function Chierchia et al. [2013], and a convex function composed with the distance to a closed convex set Combettes and Pesquet [2008].

  +∞[ and ξ ∈ ]0, +∞[ υ lim ζ→+∞ ϕ(ζ) ζ if υ ∈ ]0, +∞[ and ξ = 0 0 if υ = ξ = 0 +∞ otherwise (5.9) and ϕ : R → [0, +∞] belongs to Γ 0 (R) and is twice differentiable on ]0, +∞[. 1 So, Φ is the perspective function Bauschke and Combettes [2011] of ϕ on [0, +∞[ × ]0, +∞[ and D is a function in Γ 0 (R 2P ). If ϕ is strictly convex such that ϕ(1) = ϕ (1) = 0 (5.10) D belongs to the celebrated class of ϕ-divergences Csiszár [1963, 1967]. Then, ∀(p, q) ∈ [0, +∞[ P × [0, +∞[ P D(p, q) ≥ 0 (5.11)

  n, t1,n) + e0,n for s = 1, . . . , S rs+1,n = prox ω -1 s Rs ts+1,n + es,n yn = (ω0A A + ω0B B + S s=1 ωsT s Ts) -1 (ω0A r0,n + ω0B r1,n + S s=1 ωsT s rs+1,n) λn ∈ ]0, 2[ t0,n+1 = t0,n + λn A(2yn -xn) -r0,n t1,n+1 = t1,n + λn B(2yn -xn) -r1,n for s = 1, . . . , S ts+1,n+1 = ts+1,n + λn Ts(2yn -xn) -rs+1,n xn+1 = xn + λn(yn -xn).

Proposition 5 . 2 . 6 [

 526 Pesquet and Pustelnik, 2012, Corollary 5.3] Suppose that the following assumptions hold.(i) A A + B B + S s=1 T s T s is invertible.

  30)(with the usual convention inf ∅ = +∞ and sup ∅ = -∞). If χ -= +∞, the functionψ : ]0, +∞[ → R : ζ → ζϕ(ζ -1 ) -Θ(ζ) + γ -1 υ 2 ζ 2 -γ -1 ξζ (5.31)is strictly convex on ]χ -, +∞[. In addition, if(i) χ -= +∞ and χ + = -∞ (ii) lim ζ→χ - ζ>χ - ψ (ζ) < 0 (iii) lim ζ→χ + ψ (ζ) > 0then ψ admits a unique minimizer ζ on ]χ -, +∞[, and ζ < χ + .

Proposition 5 . 3 . 3

 533 Let γ ∈ ]0, +∞[ and (υ, ξ) ∈ R 2 . prox γΦ (υ, ξ) ∈ ]0, +∞[ 2 if and only if Conditions (i)-(iii) in Lemma 5.3.1 are satisfied. When these conditions hold,prox γΦ (υ, ξ) = υ -γ ϑ -( ζ), ξ -γ ϑ + ( ζ)(5.44)where ζ < χ + is the unique minimizer of ψ on ]χ -, +∞[.Proof. For every (υ, ξ) ∈ R 2 , such that Conditions (i)-(iii) in Lemma 5.3.1 hold, let υ = υ -γ ϑ -( ζ) (5.45) ξ = ξ -γ ϑ + ( ζ) (5.46)where the existence of ζ ∈]χ -, χ + [ is guaranteed by Lemma 5.3.1. As consequences of (5.34) and (5.36), υ and ξ are positive. In addition, since

. 47 )ξ

 47 we derive from (5.45) and (5.46) that ζ = ξ/υ > 0. This allows us to re-express (5.45) and (5.46) as υ --ξ + γ ∂Φ ∂ξ (υ, ξ) = 0. (5.51) Proposition 5.4.1 The proximity operator of γΦ with γ

  .66) According to (5.61) and Proposition 5.3.3, υ ∈ ]0, +∞[ if and only if Conditions (i)-(iii) in Lemma 5.3.1 hold. Since χ -∈ ]0, +∞[ and lim ζ→+∞ ψ (ζ) = +∞, Lemma 5.3.2 shows that these conditions are satisfied if and only ifξ < γ and exp(-υ/γ) < (1 -γ -1 ξ) exp(υ/γ) > 1 -γ -1 ξ.(5.69)Under this assumption, Proposition 5.3.3 leads to the following expressions of the proximity operator:υ = υ + γ ln ζ (5.70) ξ = ξ + γ 1 ζ -1 (5.71)where ζ is the unique minimizer on ] exp(-υ/γ), +∞[ of the function ψ.

  .76)We have then(∀ζ > exp(-υ/γ)) ψ 1 (ζ) = ln ζ + γ -1 υ + 1 > 0 (5

Fix

  ζ0 ∈ [χ-, χ+]for n = 0, 1, . . .ζn+1 = P [χ -,χ + ] ζn -ψ ( ζn)/ψ ( ζn) .

  .106) According to Proposition 5.3.3, (υ, ξ) ∈ ]0, +∞[ 2 if and only if Conditions (i)-(iii) in Lemma 5.3.1 hold. Under these conditions, Proposition 5.3.3 leads to υ = υ + γ( ζ 1/2 -1) (5.107) ξ = ξ + γ ζ -1/2 -1 (5.108) 5.4. Examples 105 where ζ is the unique minimizer on ]χ -, +∞[ of the function:

  .110) By setting ρ = ζ 1/2 , (5.102) is obtained. Since lim ζ→0 ζ>0 ψ (ζ) = -∞ and lim ζ→+∞ ψ (ζ) = +∞, Lemma 5.3.2 shows that Conditions (i)-(iii) are satisfied if and only ifυ < γ, ξ < γ, and (1 -γ -1 υ) 2 < (1 -γ -1 ξ) -2or υ < γ and ξ ≥ γ or υ ≥ γ and ξ < γ or υ ≥ γ and ξ ≥ γ (5.111)

Proposition 5 . 5 . 1

 551 Let ϕ : R → [0, +∞] be a function in Γ 0 (R) which is twice differentiable on ]0, +∞[. Let Φ be the function defined by (5.9) and ϕ * ∈ Γ 0 (R) the Fenchel-conjugate function of the restriction of ϕ on [0, +∞[, defined as(∀ζ * ∈ R) ϕ * (ζ) = sup ζ∈[0,+∞[ ζζ * -ϕ(ζ) .

  .152) From the definition of Φ, we deduce that(∀(υ, ξ) ∈ R 2 ) (5.153) Φ * (υ * , ξ * ) = sup sup (υ,ξ)∈[0,+∞[×]0,+∞[ υυ * + ξξ *ξ)∈[0,+∞[×]0,+∞[ υυ * + ξξ * -ξϕ υ ξ , 0 = sup{ι epi ϕ * (υ * , -ξ * ), 0} (5.155) = ι epi ϕ * (υ * , -ξ * ), (5.156) where the equality in (5.155) stems from [Bauschke and Combettes, 2011, Example 13.8]. Then, (5.151) follows from the conjugation property of the proximity operator (see Proposition 5.2.3 (v)).

+ y 2 .

 2 .158) We now show that Condition (5.158) holds true because b > a and it is a sum of two terms (b + 1) log a -(a + 1) log b + Proof. Equation (5.159) can be rewritten as (∀a ∈]0, +∞[)(∀b ∈]a, +∞[) g(a, b) -g(b, a) > 0 where g(x, y) = -(x + 1) log y -Therefore, we shall demonstrate that, for every b > a > 0, g is decreasing w.r.t. the first argument and increasing w.r.t. to the second argument, i.e. g(a, b) > g(b, b) and g(b, b) > g(b, a), (5.161) which implies that g(a, b) > g(b, a).

  161) are indeed equivalent to (∀y ∈]0, +∞[)(∀x ∈]0, y[) ∂g ∂x (x, y) < 0, (5.162) (∀x ∈]0, +∞[)(∀y ∈]x, +∞[) ∂g ∂y (x, y) > 0. (5.163)

§ 6 . 1

 61 Regularization based on ϕ-divergences In this section, we consider divergences for building non-local smoothness measures. The Non-Local Total Variation Gilboa and Osher [2009]; Buades et al. [2006] has been used as a popular and effective image prior model in regularization-based imaging problems Werlberger et al. [2010]; Peyré [2011]; Chierchia et al. [2013]; Couprie et al. [2013]; Teuber et al. [2013]. As we have mentioned in Chapter 2, NLTV is a powerful regularization operator. It makes full use of spatial information distributed over the different image regions.

. 6 )

 6 Let, for all s ∈ A,a (s) = A s x ∈ [0, +∞[ |Ns| and b (s) = B s x ∈ [0, +∞[ |Ns| . Thus, Equation (6.6) is equivalent to 1,p -NLTV(x) = s∈A a (s) -b (s) p . (6.7) (a) Fragment of the noisy image (Chapter 2) Image part highlighted by red color is illustrated in (b), (c) and (d). (b) Reference pixel and surrounding 5×5 bloc. (c) Candidates (d) Chosen neighboors

Figure 6 . 1 :

 61 Figure 6.1: Illustration of the idea behind nonlocal total variation.

Figure 6 . 2 :

 62 Figure 6.2: Example of angiographic (a), satellite (b,c,d,e), and medical (f,g,h) images.

Figure 6 . 3 :

 63 Figure 6.3: Angiographic + Satellite images: comparison of restoration results (SNR-SSIM) with σ 2 =64 and a blur operator of size 7×7. From the Left to the Right: Original image, degraded image, 12 -TV result, 12 -NLTV result and JK-NLTV result.

Figure 6 . 4 :

 64 Figure 6.4: Medical images: comparison of restoration results (SNR-SSIM) with σ 2 =64 and a blur operator of size 7×7. From the Left to the Right: Original image, degraded image, 12 -TV result, 12 -NLTV result and JK-NLTV result.

Figure 6 . 5 :

 65 Figure 6.5: Left images of the stereo pairs used in our simulations.

Figure 6 . 6 :

 66 Figure 6.6: Results for "Cloth" stereo pair: a)-b) stereo images, c)-d) ground truths, e)-f ) KL divergence and g)-h) JK divergence. Parameters: ((z 1 ) min , (z 1 ) max ) = (20, 75), ((z 2 ) min , (z 2 ) max ) = (0.1, 1.1), (κ 1 , κ 2 ) = (74000, 400).

  1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Taxonomy of principal distances. c Frank Nielsen, 2007. Version 0.2 Last updated, May 2008 . . . . . . . . . . . . . . . . 2.1 Example of a degraded image resulting from the direct problem: z = T x + w. . . . . . . . . . . . . . . . . . . . . . . . . 3.1 If you focus on a tree, closing first your right eye and then your left eye, the tree will appear to move relatively to your finger. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Cameras (at two different positions) capture two separate images. The disparity is the difference between the positions of the same object in the two images. If an object is far away, the disparity will be small; otherwise, the disparity will be large if the object is close. . . . . . . . . . . . . . . . . . . . . 3.3 The epipolar lines are defined by the intersection of the left and right image planes with the epipolar plane (C L , C R , P ). The points e L and e R are the epipoles and are defined by the intersection of the image planes with the line connecting the two cameras: the baseline (C L , C R ). . . . . . . . . . . . . . . 3.4 Two cameras (acquiring images of the same scene) have two different 2D representations of a common 3D point. With proper processing, the position and depth of the 3D point can be extracted from the images. . . . . . . . . . . . . . . . . . 3.5 The connection between disparity and depth. . . . . . . . . . 3.6 Stereo images: "Books" (left), "Dolls" (center), and "Parking meter" (right). . . . . . . . . . . . . . . . . . . . . . . . . . .

  139

  . . . . . . . . . . . . . . . . . . . . 4.8 The extended disparity map constructed using the sequence of disparity maps: (u 1,2 , u 2,3 , u 3,4 , u 4,5 , u 5,4 ) . . . . . . . . . . 4.9 3D reconstruction from the true disparity sequence (Left column) and from the proposed method (Right column). The camera position from the top to the buttom: Left, Down, Center, Up, and Right. . 81 6.1 Illustration of the idea behind nonlocal total variation. . . . . 6.2 Example of angiographic (a), satellite (b,c,d,e), and medical (f,g,h) images. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Angiographic + Satellite images: comparison of restoration results (SNR-SSIM) with σ 2 =64 and a blur operator of size 7×7. From the Left to the Right: Original image, degraded image, 12 -TV result, 12 -NLTV result and JK-NLTV result. 6.4 Medical images: comparison of restoration results (SNR-SSIM) with σ 2 =64 and a blur operator of size 7×7. From the Left to the Right: Original image, degraded image, 12 -TV result,

  Projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Related work . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Chapter outline . . . . . . . . . . . . . . . . . . . . . . 3.2 Stereo matching model . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Problem formulation . . . . . . . . . . . . . . . . . . . 3.2.2 Introducing prior information . . . . . . . . . . . . . . 3.2.2.1 Constraints on the disparity . . . . . . . . . 3.2.2.2 Constraints on the illumination field . . . . . 3.2.2.3 Resulting constraints applied to vector w . . 3.3 Proximal approaches for convex optimization . . . . . . . . .

	xx
	-Chapter 1 -
	Introduction
	"What is written without effort is in general read without pleasure."
	Samuel Johnson
	1 Introduction
	Introduction
	2 Convex optimization background
	2.1 inverse problems . 3 A parallel proximal splitting method for disparity estima-
	tion from multicomponent images under illumination varia-
	tion
	3.1 xix

. . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Data fidelity and prior information . . . . . . . . . . . 2.2 Algorithms in convex optimization . . . . . . . . . . . . . . . 2.2.0.1 Primal algorithms . . . . . . . . . . . . . . . 2.2.0.2 Primal-dual algorithms . . . . . . . . . . . . 2.3

  In the following, we denote by • the standard Euclidean norm, and by Id the identity matrix. The domain of a function f : R N

	13)
	For computing the projection onto the convex set C , specific numerical
	methods Chierchia et al. [2013] have been recently developed.
	§ 2.2 Algorithms in convex optimization

Notation:

  The operators (prox gr ) 1≤r≤R , (T r ) 1≤r≤R , (P Cs ) 1≤s≤S and (H s ) 1≤s≤S are used

	individually to solve Problem (2.21). This yields easily implementable algo-
	rithms, which circumvent the problem of computing the proximity operator
	of a sum of functions or a function composed with a linear operator, for
	which closed-form expressions are rarely available.
	Proximal algorithms can be classified in primal methods Daubechies
	et al. [2004]; Chaux et al. [2007]; Combettes and Pesquet [2007b]; Figueiredo
	et al. [2007]; Beck and Teboulle [2009]; Fornasier and Schönlieb [2009]; Steidl
	and Teuber [2010]; Combettes and Pesquet [2011]; Pesquet and Pustelnik
	[2012] and primal-dual methods Chen and Teboulle [1994]; Esser et al. [2010];
	Chambolle and Pock [2011]; Briceño-Arias and Combettes [2011]; Combettes
	and Pesquet [2012]; Vũ [2013]; Condat [2013]. The main difference is that
	primal algorithms generally require the computation of the inverse of the
	operator R r=1 T r T r + S s=1 H s H s (see Algorithm 1), while primal-dual
	ones only require the calculation of (T r ) 1≤r≤R , (H s ) 1≤s≤S and their adjoints
	(see Algorithm 2). Consequently, primal-dual methods are often easier to
	implement than primal ones, but their convergence may be slower Pustelnik
	et al. [2012]; Couprie et al. [2013].
	2.2.0.1 Primal algorithms
	Primal methods aim at directly solving Problem (2.21), without resorting
	to any duality arguments. In this context, forward-backward splitting and
	Douglas-Rachford splitting are popular, although they are limited to solving
	particular instances of Problem (2.21). Nonetheless, they are the building
	blocks of more sophisticated algorithms (such as PPXA+ and M+LFBF)
	that allow one to solve Problem (2.21) in its general form.
			Chaux et al. [2007];
	Combettes and Pesquet [2011]. In particular, they generalize the notion of
	projection onto a nonempty closed convex subset C of R N , in the sense that
	(∀y ∈ R N )	u∈C prox ι C y = P C (y) = argmin	u -y .	(2.23)
	Hence, proximal methods provide a unifying framework that allows one to
	address nonsmooth penalizations (R ≥ 1) and hard constraints (S ≥ 1).
	The key feature of proximal methods is that they proceed by splitting.

  In general, dense stereo techniques can be also classified as local[START_REF] Hirschmuller | Evaluation of cost functions for stereo matching[END_REF];[START_REF] Tsai | The evaluation of normalized cross correlations for defect detection[END_REF];[START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF] or global[START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF];[START_REF] Ishikawa | Exact optimization for markov random fields with convex priors[END_REF]; Boykov et al. [2011]; Klaus et al. [2006]; Yang et al. [2006]; Besse et al. [2012]; Deriche et al. [1996]; Miled et al. [2006, 2009a], depending on whether they rely on local window-based computations or on the minimization of a global energy function. In local methods, the disparity computation at a given point depends only on intensity values within a finite window. Clearly, these techniques assume that all pixels within the window have the same disparity; therefore, they are sensitive near object boundaries. Attempts to alleviate this problem include the use of adaptive windows Kanade and Okutomi [1994] and shiftable windows Kang et al. [2001].

global combinatorial optimization methods, such as graph-cuts Kolmogorov and Zabih [2001]; Ishikawa [2003]; Boykov et al. [2011]; Komodakis et al. [2011] and belief propagation Klaus et al. [2006]; Yang et al. [2006]; Besse et al. [2012]; (ii) variational approaches which have proven to be also very effective for globally solving the matching problem Deriche et al. [1996]; Miled et al. [2006, 2009a]; Pock et al. [2010]; Cremers et al. [2011]; Wanner and Goldluecke [2013].

Table 3 .

 3 1: Influence of λ n on execution time (in seconds) for PPXA+ algorithm. The stopping condition described in Section 3.4.3 is used. The machine is a quad core Intel Xeon Processor X5450, 3Ghz CPU and 8Gb RAM.

	λ n	0.5	0.7	0.9	1	1.1	1.3	1.5	1.7	1.9
	Teddy 18.16 12.97 10.09 9.09 8.26 6.98	6.29 6.77	7.35
	Cones 19.69 14.08 10.96 9.89 9.01 10.00 9.05 10.55 11.07

Table 3 .

 3 2: Comparative results based on MAE for Teddy, Venus, and Cones stereo pairs. The red superscript numbers represent method ranking.

	Method	Teddy Venus Cones
	NCC	1.169 8 0.411 8 1.134 8
	Gc+Occ Kolmogorov and Zabih [2001]	1.166 7 0.303 7 0.585 6
	ConvexTV Pock et al. [2008]	0.366 1 0.178 1 0.374 1
	GlobalGCP Wang and Yang [2011]	0.489 2 0.272 4 0.443 3
	HistoAggr Min et al. [2011]	0.788 4 0.278 5 0.434 2
	CurveletSupWgt Mukherjee et al. [2010] 0.831 5 0.299 6 0.831 7
	Subg. proj. Miled et al. [2009a]	0.890 6 0.207 2 0.551 5
	PPXA+	0.666 3 0.211 3 0.487 4

  3.5. Experiments with non local constraints 53 § 3.5 Experiments with non local constraintsIn this part, we exploit a recent extension of the total variation functional which is the Non Local TV (NLTV). Here, the nonlocality principal means that any point in the depth map can interact directly with any other point in the whole map. The main idea of nonlocal extension is based on the definition of nonlocal derivative operators, where the objective is to take into account additional prior information derived from the object itself. Such pixel grouping has already been considered in the field of image restoration, stereo matching, and motion estimation: when performing image restoration, the support weights are computed based on the fact that natural images exhibit a large amount of redundancy. This leads to the well-known neighborhood filters[START_REF] Buades | Nonlocal image and movie denoising[END_REF] which belong to the best performing methods. A variational interpretation of these ideas leads to the nonlocal regularizers Gilboa and Osher [2009]; Peyré [2011]; Chierchia et al. [2013]; Couprie et al.

[2013]

. In the context of stereo matching, Yoon and Kweon

[START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF]

computed adaptive support weights to locally adapt the local matching window to the image structure, which leads to a significant improvement of the method. In Xiao et al.

[2006]

, the authors proposed to incorporate the bilateral Gaussian filter into a variational motion estimation framework in order to improve the performance in occluded regions. In Werlberger et al.

Table 3 . 3 :

 33 Disparity estimation results obtained for the gray level images Venus and Bull, using an 1 cost.

		TV	NLTV	Parameters
	Venus MAE=0.21 MAE=0.18 τ 2 =14000,τ 4 =15600
		Err 2 =1.52	Err 2 =0.69 (u min , u max ) =(3,20)
	Bull	MAE=0.20 MAE= 0.18	τ 2 =2533,τ 4 =7458
		Err 2 =1.47	Err 2 = 1.16 (u min , u max ) =(3,20)

Table 3 . 4 :

 34 Disparity estimation results obtained for the color stereo pairs Cloth and Wood, under illumination variation using an 1 cost.

	TV	NLTV	Parameters
	Cloth MAE=0.84 MAE=0.82 τ 2 =45000, τ 4 =18800
	Err 2 =2.44	Err 2 =1.73 (u min , u max ) =(15,55)
			(v min , v max )=(0.1,0.6)
	Wood MAE=2.78 MAE=2.53 τ 2 =58000, τ 4 =21000
	Err 2 =10.47 Err 2 =8.10 (u min , u max ) =(21,71)
			(v min , v max )=(0.1,0.6)

  .8) Note that f n and tv are convex functions, but the problem is nonconvex due to the nonconvexity of the set B. We thus use a convex relaxation based on the approach of Cremers et al.[2011]. The problem is then equivalent to minimize

θn∈R f n (θ n ) + µ Q q=1

Table 4

 4

		1	1,	1 2	1 2 ,
	(I2, I4)	(0.64, 5.70)	=25, (0.56, 4.29) (0.60, 5.37)	=19, (0.64, 5.90)
	(I0, I2, I4)	(0.52, 5.33)	=23, (0.48, 4.08) (0.50, 5.02)	=19, (0.61, 5.67)
	(I0, I1, I2, I3, I4) (0.46, 4.59)	=23, (0.48, 3.82) (0.44, 4.33)	=19, (0.61, 5.64)

.2: The (MAE, Err 1 ) results on Teddy image, with Q=20, r0=6, rQ=26 and (∀ q ∈ {2, . . . , Q}) rq -rq-1 = 1, and η = 23 × 10 3 .

Table 4 . 3

 43

		1	1,	1 2	1 2 ,
	(I2, I4)	(0.49, 7.34)	=27, (0.41, 5.24) (0.49, 6.96)	=25, (0.59, 7.9)
	(I0, I2, I4)	(0.49, 4.71)	=27, (0.45, 3.72) (0.42, 4.44)	=25, (0.49, 5.23)
	(I0, I1, I2, I3, I4) (0.44, 4.04)	=27, (0.43, 3.68) (0.43, 3.81)	=25, (0.48, 5.19)

: The (MAE, Err 1 ) results on Cones image, with Q=26, r0=2, rQ=28 and (∀ q ∈ {2, . . . , Q}) rq -rq-1 = 1, and η = 27 × 10 3 .

Table 4 .

 4 4: (MAE, Err 1 ) results on u2,4 for three different stereo pairs using ψ = | • | and ψ = | • | 2 . Comparison is performed between the first-Order Taylor linearization and the improved one.

	Stereo pair	τ 2	[u min ,u max ] ψ = | • | p	1st-order	Improved 1st-order
	Sawtooth	14×10 3	[3, 18]	p = 1 p = 2	(0.29, 2.81) (0.34, 3.12)	(0.26, 2.79) (0.30, 3.10)
	Teddy	46×10 3	[0, 26]	p = 1 p = 2	(0.61, 7.81) (0.62, 7.93)	(0.60, 7.64) (0.62, 7.74)
	Cones	68×10 3	[0, 27]	p = 1 p = 2	(0.63, 7.84) (0.69, 7.80)	(0.66, 7.41) (0.62, 7.61)

Table 4 . 5 :

 45 Percentage of bad pixels obtained with the considered methods. , I 4 ) (I 0 , I 2 , I 4 ) (I 0 , I 1 , I 2 , I 3 , I 4 ) , I 4 ) (I 2 , I 3 , I 4 ) (I 1 , I 2 , I 3 , I 4 , I 5 )

	u 2,4 (I 2 Sawtooth Multilabel 1.23 Improved 1st-order 2.79	1.16 2.13	1.08 1.86
		Graph cuts	1.75	1.66	1.70
	Teddy	Multilabel Improved 1st-order	4.29 7.64	4.08 5.43	3.82 5.00
		Graph cuts	4.57	4.46	3.83
	Cones	Multilabel Improved 1st-order	5.24 7.41	3.72 5.4	3.68 4.65
		Graph cuts	7.02	6.02	4.65
	u 3,4 (I 3 Sponza Multilabel 5.4 Improved 1st-order 7.71	4.2 5.93	3.10 5.49
		Graph cuts	6.50	5.47	4.84

  5.2.1, various proximal methods can be employed to solve it numerically. Two examples of such methods are given subsequently. The first algorithm is PPXA+ Pesquet and Pustelnik [2012] which constitutes an extension of PPXA (Parallel ProXimal Agorithm) proposed in Combettes and Pesquet [2008]. As can be seen in Afonso et al. [2011]; Setzer et al. [2010], PPXA+ is also strongly related to augmented Lagrangian methods [Pesquet and Pustelnik, 2012, Section 6]. In the form provided below, this algorithm requires A A + B B + S s=1 T s T s to be invertible. 3

	Algorithm 6 PPXA+
	Initialization
	(ω0, . . . , ωS)

  .96) According to Proposition 5.3.3, prox γΦ (υ, ξ) ∈ ]0, +∞[ 2 if and only if Conditions (i)-(iii) in Lemma 5.3.1 hold. Lemma 5.3.2 shows that these conditions are satisfied if and only if

Table 5 .

 5

1: Conjugate function ϕ * of the restriction of ϕ to [0, +∞[. § 5.6 Conclusion

Table 6 .

 6 2: (SNR, SSIM) restoration results with σ 2 =64 and a blur operator of size 3×3.

Table 6 .

 6 7: Disparity estimation results (MAE, Err 2 ) obtained for gray level images. Hel (0.9057, 3.46) (0.9333, 5.57) (1.7631, 11.71) (1.2107, 13.87) Chi (0.8335, 3.43) (0.9328, 5.56) (1.7636, 11.71) (1.2127, 13.81) I0.2 (0.9058, 3.39) (0.9335, 5.55) (1.7741, 11.64) (1.2167, 13.83) § 6.3 Conclusion

		Cloth	Baby 1	Baby 2	Dolls
	2	(0.8385, 3.62) (0.9706, 5.80) (1.8336, 12.18) (1.2591, 14.41)
	KL	(0.8287, 3.36) (0.9301, 5.40) (1.6539, 11.56) (1.2107, 13.86)
	JK	(0.8348, 3.44) (0.9298, 5.32) (1.6503, 11.48) (1.2157, 13.36)

• La I α -divergence a été initialement proposée pour évaluer statistiquement l'efficacité d'un test d'hypothèse. Elle a été étendue par de nombreux chercheurs et est surtout considérée dans le contexte de la factorisation des matrices non-négatives.• La divergence de Rényi a été introduite par Alfréd Rényi comme une mesure de l'information liée à l'entropie Rényi. Elle a été étudiée et appliquée dans de nombreux domaines, y compris les problèmes de recalage d'image.Etant des fonctions de dissimilarité, les ϕ-divergences peuvent servir comme critères de minimisation dans les problèmes d'optimisation. Cependant, ce sont des fonctions non lisses, ce qui exclut les techniques d'optimisation lisses. De plus, les approches existantes impliquant ce type de mesures se limitent souvent à des problèmes où la minimisation est effectuée par rapport à l'un des arguments de la divergence, ou en utilisant des approches alternées qui nécessitent des hypothèses spécifiques pour être valides.

x -u ∈ ∂g1(x) ⇔ (2x -u) -x ∈ ∂g1(x) ⇔ x = prox g 1 (2x -u)

In (2.29) PD(2 v [n] -u [n] ) = 2 v [n] -u [n] , with v [n] = PD(v [n] ) and u [n] = PD(u [n] ) = PD u [n-1] + λn (2 v [n-1] -u [n-1] -v [n-1] ) = u [n-1] + λn( v [n-1] -u [n-1] ).

0 ∈ ∂g * (v) -Lx ⇔ (v + Lx) -v ∈ ∂g * (v) ⇔ v = prox g * (v + Lx).

http://www.cs.ubc.ca/labs/scl/SPGL11/

http://www.stanford.edu/ ~boyd/l1_ls/.

thttp://www.lx.it.pt/ ~mtf/SpaRSA/

http://vision.middlebury.edu/stereo/

http://vasc.ri.cmu.edu/idb/html/jisct/

In this section, the optimization is performed w.r.t. one variable instead of two. Recall that S1,1 (resp. S1,2) can be expressed as L1(C1,1) (resp. L2(C1,2)), where L1 is the identity matrix and L2 the gradient operator matrix in R 2P .

Note that, in Tables

4.1, 4.2, and 4.3, the five multi-view images (3rd line) contain the three multi-view images used in the previous experiment (2nd line), which in turn contain the left/right stereo images (1st line).

http://software.intel.com/en-us/articles/embree-photo-realistic-ray-tracing-kernels-0

http://www.robots.ox.ac.uk/∼ojw/software.htm

The existence of lim ζ→+∞ ϕ(ζ)/ζ in [0, +∞] is guaranteed [Hiriart-Urruty and Lemaréchal, 1996, Section

2.3].

Indeed, none of the considered ϕ-divergences can be expressed as a function of the difference between the two arguments.

A slightly more complex variant of PPXA+ exists where this assumption is not required Attouch and Soueycatt [2009]; Pesquet and Pustelnik [2012], Prop. 5.4.

Note also that we have employed PPXA+ (resp. M+LFBF) in Chapter 3 (resp. Chapter 4), where the optimization was performed with respect to one variable.

www.pcronline.com

www.cellimagelibrary.org/

http://vision.middlebury.edu/stereo/

-NLTV result and JK-NLTV result. . . . . . . . . . . . . . 6.5 Left images of the stereo pairs used in our simulations. . . . . 6.6 Results for "Cloth" stereo pair: a)-b) stereo images, c)-d) groundtruths, e)-f) KL divergence and g)-h) JK divergence. Parameters: ((z 1 ) min , (z 1 ) max ) = (20, 75), ((z 2 ) min , (z 2 ) max ) = (0.1, 1.1), (κ 1 , κ 2 ) = (74000, 400). . . . . . . . . . . . . . . . . . . . . . .

Acknowledgments

I would like to acknowledge to my supervisors Jean-Christophe Pesquet and

Chapter 5. Proximity operators of discrete information divergences

Numerous additional properties of proximity operators are mentioned in [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]; [START_REF] Parikh | Proximal algorithms[END_REF].

In this chapter, we will be mainly concerned with the determination of the proximity of D as introduced in Section 5.2.2 (in this case, H = R P ×R P ). The next result shows that the problem reduces to the determination of the proximity operator of a real function of two variables. Proposition 5.2. 4 Let D be defined by (5.8) where Φ ∈ Γ 0 (R 2 ) and let γ ∈ ]0, +∞[. Let u = (u (i) ) 1≤i≤P ∈ R P and v = (v (i) ) 1≤i≤P ∈ R P . Then, for every p = (p (i) ) 1≤i≤P ∈ R P and for every q = (q (i) ) 1≤i≤P ∈ R P , prox γD(•+u,•+v) (p, q) = (p (i) -u (i) ) 1≤i≤P , (q (i) -v (i) ) 1≤i≤P (5.18) where (∀i ∈ {1, . . . , P }) (p (i) , q (i) ) = prox γΦ (p (i) + u (i) , q (i) + v (i) ). (5.19) Proof. The result is a straightforward consequence of [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] and Proposition 5.2.3(iii) (by setting f = D and z = (u, v)).

Note that, although an extensive list of proximity operators of one-variable real functions can be found in [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], few results are available for real functions of two variables [START_REF] Benamou | A computational fluid mechanics solution to the monge-kantorovich mass transfer problem[END_REF]; [START_REF] Benamou | The monge-katonrovich mass transfer and its computational fluid mechanics formulation[END_REF]; [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]; [START_REF] Chierchia | Epigraphical splitting for solving constrained convex optimization problems with proximal tools[END_REF].

An example of such a result is provided below.

Proposition 5.2.5 Let ϕ ∈ Γ 0 (R) be an even differentiable function on

Then, for every

otherwise.

(5.20)

Proof. Let (ν, ξ) ∈ R 2 . From Proposition 5.2.3(i), we know that prox Φ (ν, ξ) ∈ [0, +∞[ 2 . By using Proposition 5.2.3(ii), we have the following equivalences:

Problem formulation

91

where Φ : (ν, ξ) → ϕ(ν -ξ). By using now Proposition 5.2.3(vi), we get

(5.22)

Similarly, we have

(5. 23) where the last equivalence results from the assumption that ϕ is even. Symmetrically,

(5.24)

The above two propositions provide a simple characterization of the proximity operators of some distances defined for nonnegative-valued vectors. However, the assumptions made in Proposition 5.2.5 are not satisfied by the class of functions Φ considered in Section 5.2.2. 2 In the next section, we will propose two algorithms for solving a general class of convex problems involving these functions Φ.

In this algorithm, (γ n ) n≥0 is a sequence of step-sizes, and e 0,n ∈ (R P ) 2 , e 1,n ∈ R K 1 , . . . , e S,n ∈ R K S still correspond to possible errors in the computation of the proximity operators. The convergence of the algorithm is secured by the following result. Proposition 5.2.7 [Combettes and Pesquet, 2012, Theorem 4.2] Suppose that the following assumptions hold.

If the set of solutions to Problem (5.2.1) is nonempty, then any sequence (x n ) n∈N generated by Algorithm (10) converges to an element of this set.

It is worth highlighting that these two algorithms share two interesting features: at first, many operations (e.g. the loops on s) can be implemented in 96 Chapter 5. Proximity operators of discrete information divergences This allows us to deduce that, if

(5.34)

Similarly, the function ϑ + is increasing as the convexity of ϕ yields

So, Conditions (ii) and (iii) are equivalent to

(5.40)

In view of (5.34) and (5.36), these inequalities are satisfied if and only if χ -< χ + . The inequality is also obviously satisfied if χ -= 0 or χ + = +∞.

In addition, we have The required assumptions in the previous lemma can often be simplified as stated below. Conversely, for every (υ, ξ) ∈ R 2 , let (υ, ξ) = prox γΦ (υ, ξ). If (υ, ξ) ∈ ]0, +∞[ 2 , (υ, ξ) satisfies (5.48) and (5.49). By setting ζ = ξ/υ > 0, after simple calculations, we find

(5.54)

According to (5.34) and (5.36), (5.52) and (5.53) imply that χ -= +∞, χ + = -∞, and ζ ∈]χ -, χ + [. In addition, according to Lemma 5.3.1, ψ is strictly increasing on ]χ -, +∞[ (since ψ is strictly convex on this interval). Hence, ψ has a limit at χ -(which may be equal to -∞ when χ -= -∞), and Condition (ii) is satisfied. Similarly, ψ has a limit at χ + (possibly equal to +∞ when χ + = +∞), and Condition (iii) is satisfied.

Remark 5.3.4 In (5.9), a special case arises when

where ϕ is a twice differentiable convex function on ]0, +∞[. Then Φ takes a symmetric form, so giving birth to L-divergences. It can then be deduced from (5.28) that

(5.56) § 5.4 Examples

Kullback-Leibler divergence

Let us now apply the results in the previous section to the function

(5.57) This is a function in Γ 0 (R 2 ) satisfying (5.9) with

(5.58) Chapter 5. Proximity operators of discrete information divergences Algorithm 8 Douglas-Rachford algorithm

(iii) More generally, we can derive the proximity operator of

where prox γΦ is provided by Proposition 5.4.1.

Jeffreys-Kullback divergence

Let us now consider the symmetrized form of (5.57) given by

(5.87)

This function belongs to Γ 0 (R 2 ) and satisfies (5.9) and (5.55) with

where, in the first case, ρ is the unique solution on ]0, 1 + γ -1 υ/2[ of the equation:

(5.117)

Proof. By proceeding similarly to the proof of Proposition 5.4.5. we have:

(5.120) and

+∞ otherwise

(5.121)

(

where ζ is the unique minimizer on ]χ -, +∞[ of the function:

Thus, ζ is the unique solution on ]χ -, +∞[ of the equation:

or υ > -2γ and ξ ≥ γ.

(5.128) When (5.128) does not hold, we necessarily have υ = 0. The end of the proof is similar to that of Proposition 5.4.5.

Renyi divergence

Let α ∈]1, +∞[ and let us now consider the function of Γ 0 (R 2 ) given by

which corresponds to the case when

(5.130) Note that the above function Φ allows us to generate the Rényi divergence up to a log transform and a multiplicative constant.

Proposition 5.4.7 The proximity operator of γΦ with γ ∈ ]0, +∞[ is

where, in the first case, ζ is the unique solution on ](αγυ -1 )

(5.132)

Chapter 5. Proximity operators of discrete information divergences

Proof. We proceed similarly to the previous examples by noticing that

and

(

Note that (5.132) becomes a polynomial equation when α is a rational number. In particular, when α = 2, it reduces to the cubic equation:

I α divergence

Let α ∈]0, 1[ and let us now consider the function of Γ 0 (R 2 ) given by

where, in the first case, ζ is the unique solution on max 1 -υ γα , 0

of the equation:

(5.143)

Proof. We have then (∀ζ ∈ ]0, +∞[)

and

(5.148)

+∞

otherwise.

( Algorithm 10

Results

In this section, we present the performance of the proposed D-NLTV regularization method in restoration experiments. In particular, the following choices for D are considered: squared 2 -norm, 1,2 -norm, Kullback-Leibler divergence, Jeffreys-Kullback, Hellinger divergence, Chi square divergence, and I α divergence (α=0.2). Note also that Rényi divergence is not suitable for the considered application, since its use is limited to probability distributions.

In our experiments, we have considered three types of images: angiographic, 1 medical, 2 and satellite images (Figure 6.2). § 6.2 Data fidelity based on ϕ-divergences

Image registration

The objective of image registration is to determine spatial transformations that maximize a similarity metric between two images resulting from two different acquisitions. Let the two original images be represented by data vectors I 1 ∈ R P and I 2 ∈ R P . For every j ∈ {1, 2}, let F j be an image mapping operator such that F j : R P × R N j → R P (I j , z j ) → F j (I j , z j ) (6.18) where F j (I j , z j ) is the mapped image, and z j ∈ R N j is a vector of mapping parameters He et al. [2003]; [START_REF] Hamza | Image registration and segmentation by maximizing the Jensen-Rnyi divergence[END_REF]; [START_REF] Zitovà | Image registration methods: a survey[END_REF]; [START_REF] Malviya | Wavelet based image registration using mutual information[END_REF]. When a Kullback metric is adopted, the image registration criterion to be minimized w.r.t. x = z 1 , z 2 reads D F 1 (I 1 , z 1 ), F 2 (I 2 , z 2 ) .

To determine the optimal parameter vector x ∈ R N with N = N 1 + N 2 , one can proceed by supposing that F 1 (resp. F 2 ) is differentiable and by performing a first-order Taylor expansion around an initial estimate z1 (resp. z2 ) as follows: (∀j ∈ {1, 2})

F j (I j , z j ) F j (I j , zj ) + ∂F j ∂z j (I j , zj )(z j -zj ) (6.19)